Refactor/histograms #203

Merged
altavir merged 6 commits from refactor/histograms into dev 2021-02-18 18:05:04 +03:00
23 changed files with 533 additions and 422 deletions

View File

@ -41,6 +41,7 @@
- Refactor histograms. They are marked as prototype
- `Complex` and related features moved to a separate module `kmath-complex`
- Refactor AlgebraElement
- Add `out` projection to `Buffer` generic
### Deprecated

View File

@ -76,6 +76,14 @@ benchmark {
iterationTimeUnit = "ms" // time unity for iterationTime, default is seconds
include("DotBenchmark")
}
configurations.register("expressions") {
warmups = 1 // number of warmup iterations
iterations = 3 // number of iterations
iterationTime = 500 // time in seconds per iteration
iterationTimeUnit = "ms" // time unity for iterationTime, default is seconds
include("ExpressionsInterpretersBenchmark")
}
}
kotlin.sourceSets.all {

View File

@ -10,7 +10,6 @@ import kscience.kmath.linear.Matrix
import kscience.kmath.linear.MatrixContext
import kscience.kmath.linear.inverseWithLup
import kscience.kmath.linear.real
import kscience.kmath.operations.invoke
import org.openjdk.jmh.annotations.Scope
import org.openjdk.jmh.annotations.State
import kotlin.random.Random
@ -34,14 +33,14 @@ internal class LinearAlgebraBenchmark {
@Benchmark
fun cmLUPInversion() {
CMMatrixContext {
with(CMMatrixContext) {
inverse(matrix)
}
}
@Benchmark
fun ejmlInverse() {
EjmlMatrixContext {
with(EjmlMatrixContext) {
inverse(matrix)
}
}

View File

@ -2,7 +2,6 @@ package kscience.kmath.benchmarks
import kscience.kmath.nd.*
import kscience.kmath.operations.RealField
import kscience.kmath.operations.invoke
import kscience.kmath.structures.Buffer
import org.openjdk.jmh.annotations.Benchmark
import org.openjdk.jmh.annotations.Scope
@ -12,7 +11,7 @@ import org.openjdk.jmh.annotations.State
internal class NDFieldBenchmark {
@Benchmark
fun autoFieldAdd() {
autoField {
with(autoField) {
var res: NDStructure<Double> = one
repeat(n) { res += one }
}
@ -20,7 +19,7 @@ internal class NDFieldBenchmark {
@Benchmark
fun specializedFieldAdd() {
specializedField {
with(specializedField) {
var res: NDStructure<Double> = one
repeat(n) { res += 1.0 }
}
@ -29,7 +28,7 @@ internal class NDFieldBenchmark {
@Benchmark
fun boxingFieldAdd() {
genericField {
with(genericField) {
var res: NDStructure<Double> = one
repeat(n) { res += 1.0 }
}

View File

@ -2,7 +2,6 @@ package kscience.kmath.benchmarks
import kscience.kmath.nd.*
import kscience.kmath.operations.RealField
import kscience.kmath.operations.invoke
import kscience.kmath.viktor.ViktorNDField
import org.jetbrains.bio.viktor.F64Array
import org.openjdk.jmh.annotations.Benchmark
@ -21,7 +20,7 @@ internal class ViktorBenchmark {
@Benchmark
fun automaticFieldAddition() {
autoField {
with(autoField) {
var res: NDStructure<Double> = one
repeat(n) { res += 1.0 }
}
@ -29,7 +28,7 @@ internal class ViktorBenchmark {
@Benchmark
fun realFieldAddition() {
realField {
with(realField) {
var res: NDStructure<Double> = one
repeat(n) { res += 1.0 }
}
@ -37,7 +36,7 @@ internal class ViktorBenchmark {
@Benchmark
fun viktorFieldAddition() {
viktorField {
with(viktorField) {
var res = one
repeat(n) { res += 1.0 }
}

View File

@ -2,7 +2,6 @@ package kscience.kmath.benchmarks
import kscience.kmath.nd.*
import kscience.kmath.operations.RealField
import kscience.kmath.operations.invoke
import kscience.kmath.viktor.ViktorNDField
import org.jetbrains.bio.viktor.F64Array
import org.openjdk.jmh.annotations.Benchmark
@ -22,7 +21,7 @@ internal class ViktorLogBenchmark {
@Benchmark
fun realFieldLog() {
realField {
with(realField) {
val fortyTwo = produce { 42.0 }
var res = one
repeat(n) { res = ln(fortyTwo) }
@ -31,7 +30,7 @@ internal class ViktorLogBenchmark {
@Benchmark
fun viktorFieldLog() {
viktorField {
with(viktorField) {
val fortyTwo = produce { 42.0 }
var res = one
repeat(n) { res = ln(fortyTwo) }

View File

@ -21,8 +21,17 @@ public class MstExpression<T, out A : Algebra<T>>(public val algebra: A, public
null
} ?: arguments.getValue(StringSymbol(value))
override fun unaryOperationFunction(operation: String): (arg: T) -> T = algebra.unaryOperationFunction(operation)
override fun binaryOperationFunction(operation: String): (left: T, right: T) -> T = algebra.binaryOperationFunction(operation)
override fun unaryOperation(operation: String, arg: T): T =
algebra.unaryOperation(operation, arg)
override fun binaryOperation(operation: String, left: T, right: T): T =
algebra.binaryOperation(operation, left, right)
override fun unaryOperationFunction(operation: String): (arg: T) -> T =
algebra.unaryOperationFunction(operation)
override fun binaryOperationFunction(operation: String): (left: T, right: T) -> T =
algebra.binaryOperationFunction(operation)
@Suppress("UNCHECKED_CAST")
override fun number(value: Number): T = if (algebra is NumericAlgebra<*>)

View File

@ -4,7 +4,7 @@ public abstract interface class kscience/kmath/domains/Domain {
}
public final class kscience/kmath/domains/HyperSquareDomain : kscience/kmath/domains/RealDomain {
public synthetic fun <init> ([D[DLkotlin/jvm/internal/DefaultConstructorMarker;)V
public fun <init> (Lkscience/kmath/structures/Buffer;Lkscience/kmath/structures/Buffer;)V
public fun contains (Lkscience/kmath/structures/Buffer;)Z
public fun getDimension ()I
public fun getLowerBound (I)Ljava/lang/Double;

View File

@ -16,6 +16,7 @@
package kscience.kmath.domains
import kscience.kmath.linear.Point
import kscience.kmath.structures.Buffer
import kscience.kmath.structures.RealBuffer
import kscience.kmath.structures.indices
@ -25,20 +26,20 @@ import kscience.kmath.structures.indices
*
* @author Alexander Nozik
*/
public class HyperSquareDomain(private val lower: RealBuffer, private val upper: RealBuffer) : RealDomain {
public class HyperSquareDomain(private val lower: Buffer<Double>, private val upper: Buffer<Double>) : RealDomain {
public override val dimension: Int get() = lower.size
public override operator fun contains(point: Point<Double>): Boolean = point.indices.all { i ->
point[i] in lower[i]..upper[i]
}
public override fun getLowerBound(num: Int, point: Point<Double>): Double? = lower[num]
public override fun getLowerBound(num: Int, point: Point<Double>): Double = lower[num]
public override fun getLowerBound(num: Int): Double? = lower[num]
public override fun getLowerBound(num: Int): Double = lower[num]
public override fun getUpperBound(num: Int, point: Point<Double>): Double? = upper[num]
public override fun getUpperBound(num: Int, point: Point<Double>): Double = upper[num]
public override fun getUpperBound(num: Int): Double? = upper[num]
public override fun getUpperBound(num: Int): Double = upper[num]
public override fun nearestInDomain(point: Point<Double>): Point<Double> {
val res = DoubleArray(point.size) { i ->

View File

@ -78,8 +78,7 @@ public interface NDStructure<T> {
strides: Strides,
bufferFactory: BufferFactory<T> = Buffer.Companion::boxing,
initializer: (IntArray) -> T,
): NDBuffer<T> =
NDBuffer(strides, bufferFactory(strides.linearSize) { i -> initializer(strides.index(i)) })
): NDBuffer<T> = NDBuffer(strides, bufferFactory(strides.linearSize) { i -> initializer(strides.index(i)) })
/**
* Inline create NDStructure with non-boxing buffer implementation if it is possible
@ -87,15 +86,13 @@ public interface NDStructure<T> {
public inline fun <reified T : Any> auto(
strides: Strides,
crossinline initializer: (IntArray) -> T,
): NDBuffer<T> =
NDBuffer(strides, Buffer.auto(strides.linearSize) { i -> initializer(strides.index(i)) })
): NDBuffer<T> = NDBuffer(strides, Buffer.auto(strides.linearSize) { i -> initializer(strides.index(i)) })
public inline fun <T : Any> auto(
type: KClass<T>,
strides: Strides,
crossinline initializer: (IntArray) -> T,
): NDBuffer<T> =
NDBuffer(strides, Buffer.auto(type, strides.linearSize) { i -> initializer(strides.index(i)) })
): NDBuffer<T> = NDBuffer(strides, Buffer.auto(type, strides.linearSize) { i -> initializer(strides.index(i)) })
public fun <T> build(
shape: IntArray,
@ -106,8 +103,7 @@ public interface NDStructure<T> {
public inline fun <reified T : Any> auto(
shape: IntArray,
crossinline initializer: (IntArray) -> T,
): NDBuffer<T> =
auto(DefaultStrides(shape), initializer)
): NDBuffer<T> = auto(DefaultStrides(shape), initializer)
@JvmName("autoVarArg")
public inline fun <reified T : Any> auto(
@ -120,8 +116,7 @@ public interface NDStructure<T> {
type: KClass<T>,
vararg shape: Int,
crossinline initializer: (IntArray) -> T,
): NDBuffer<T> =
auto(type, DefaultStrides(shape), initializer)
): NDBuffer<T> = auto(type, DefaultStrides(shape), initializer)
}
}

View File

@ -21,7 +21,7 @@ public typealias MutableBufferFactory<T> = (Int, (Int) -> T) -> MutableBuffer<T>
*
* @param T the type of elements contained in the buffer.
*/
public interface Buffer<T> {
public interface Buffer<out T> {
/**
* The size of this buffer.
*/

View File

@ -1,4 +1,10 @@
plugins { id("ru.mipt.npm.mpp") }
plugins {
id("ru.mipt.npm.mpp")
}
kscience {
useAtomic()
}
kotlin.sourceSets {
commonMain {

View File

@ -1,20 +1,49 @@
package kscience.kmath.histogram
/*
import kotlinx.atomicfu.atomic
import kotlinx.atomicfu.getAndUpdate
import kscience.kmath.operations.RealField
import kscience.kmath.operations.Space
/**
* Common representation for atomic counters
* TODO replace with atomics
*/
public expect class LongCounter() {
public fun decrement()
public fun increment()
public fun reset()
public fun sum(): Long
public fun add(l: Long)
public interface Counter<T : Any> {
public fun add(delta: T)
public val value: T
public companion object{
public fun real(): ObjectCounter<Double> = ObjectCounter(RealField)
}
}
public expect class DoubleCounter() {
public fun reset()
public fun sum(): Double
public fun add(d: Double)
public class IntCounter : Counter<Int> {
private val innerValue = atomic(0)
Zelenyy commented 2021-02-15 12:19:09 +03:00 (Migrated from github.com)
Review

Why not inline class?

Why not inline class?
altavir commented 2021-02-15 12:24:35 +03:00 (Migrated from github.com)
Review

It is an interface. Atomics should not be exposed to the outside world. Also, inlining does not give anything here.

It is an interface. Atomics should not be exposed to the outside world. Also, inlining does not give anything here.
override fun add(delta: Int) {
innerValue += delta
}
override val value: Int get() = innerValue.value
}
public class LongCounter : Counter<Long> {
private val innerValue = atomic(0L)
override fun add(delta: Long) {
innerValue += delta
}
override val value: Long get() = innerValue.value
}
public class ObjectCounter<T : Any>(public val space: Space<T>) : Counter<T> {
private val innerValue = atomic(space.zero)
override fun add(delta: T) {
innerValue.getAndUpdate { space.run { it + delta } }
}
override val value: T get() = innerValue.value
}

View File

@ -6,18 +6,16 @@ import kscience.kmath.structures.ArrayBuffer
import kscience.kmath.structures.RealBuffer
/**
* The bin in the histogram. The histogram is by definition always done in the real space
* The binned data element. Could be a histogram bin with a number of counts or an artificial construct
*/
public interface Bin<T : Any> : Domain<T> {
/**
* The value of this bin.
*/
public val value: Number
public val center: Point<T>
}
public interface Histogram<T : Any, out B : Bin<T>> : Iterable<B> {
public interface Histogram<T : Any, out B : Bin<T>> {
/**
* Find existing bin, corresponding to given coordinates
*/
@ -27,28 +25,31 @@ public interface Histogram<T : Any, out B : Bin<T>> : Iterable<B> {
* Dimension of the histogram
*/
public val dimension: Int
public val bins: Iterable<B>
}
public interface MutableHistogram<T : Any, out B : Bin<T>> : Histogram<T, B> {
public fun interface HistogramBuilder<T : Any> {
/**
* Increment appropriate bin
*/
public fun putWithWeight(point: Point<out T>, weight: Double)
public fun putValue(point: Point<out T>, value: Number)
public fun put(point: Point<out T>): Unit = putWithWeight(point, 1.0)
}
public fun <T : Any> MutableHistogram<T, *>.put(vararg point: T): Unit = put(ArrayBuffer(point))
public fun <T : Any, B : Bin<T>> HistogramBuilder<T>.put(point: Point<out T>): Unit = putValue(point, 1.0)
public fun MutableHistogram<Double, *>.put(vararg point: Number): Unit =
public fun <T : Any> HistogramBuilder<T>.put(vararg point: T): Unit = put(ArrayBuffer(point))
public fun HistogramBuilder<Double>.put(vararg point: Number): Unit =
put(RealBuffer(point.map { it.toDouble() }.toDoubleArray()))
public fun MutableHistogram<Double, *>.put(vararg point: Double): Unit = put(RealBuffer(point))
public fun <T : Any> MutableHistogram<T, *>.fill(sequence: Iterable<Point<T>>): Unit = sequence.forEach { put(it) }
public fun HistogramBuilder<Double>.put(vararg point: Double): Unit = put(RealBuffer(point))
public fun <T : Any> HistogramBuilder<T>.fill(sequence: Iterable<Point<T>>): Unit = sequence.forEach { put(it) }
/**
* Pass a sequence builder into histogram
*/
public fun <T : Any> MutableHistogram<T, *>.fill(block: suspend SequenceScope<Point<T>>.() -> Unit): Unit =
public fun <T : Any> HistogramBuilder<T>.fill(block: suspend SequenceScope<Point<T>>.() -> Unit): Unit =
fill(sequence(block).asIterable())

View File

@ -0,0 +1,76 @@
package kscience.kmath.histogram
import kscience.kmath.domains.Domain
import kscience.kmath.linear.Point
import kscience.kmath.misc.UnstableKMathAPI
import kscience.kmath.nd.NDSpace
import kscience.kmath.nd.NDStructure
import kscience.kmath.nd.Strides
import kscience.kmath.operations.Space
import kscience.kmath.operations.SpaceElement
import kscience.kmath.operations.invoke
/**
* A simple histogram bin based on domain
*/
public data class DomainBin<T : Comparable<T>>(
public val domain: Domain<T>,
public override val value: Number,
) : Bin<T>, Domain<T> by domain
@OptIn(UnstableKMathAPI::class)
public class IndexedHistogram<T : Comparable<T>, V : Any>(
override val context: IndexedHistogramSpace<T, V>,
public val values: NDStructure<V>,
) : Histogram<T, Bin<T>>, SpaceElement<IndexedHistogram<T, V>, IndexedHistogramSpace<T, V>> {
override fun get(point: Point<T>): Bin<T>? {
val index = context.getIndex(point) ?: return null
return context.produceBin(index, values[index])
}
override val dimension: Int get() = context.strides.shape.size
override val bins: Iterable<Bin<T>>
get() = context.strides.indices().map {
context.produceBin(it, values[it])
}.asIterable()
}
/**
* A space for producing histograms with values in a NDStructure
*/
public interface IndexedHistogramSpace<T : Comparable<T>, V : Any> : Space<IndexedHistogram<T, V>> {
//public val valueSpace: Space<V>
public val strides: Strides
public val histogramValueSpace: NDSpace<V, *> //= NDAlgebra.space(valueSpace, Buffer.Companion::boxing, *shape),
/**
* Resolve index of the bin including given [point]
*/
public fun getIndex(point: Point<T>): IntArray?
/**
* Get a bin domain represented by given index
*/
public fun getDomain(index: IntArray): Domain<T>?
public fun produceBin(index: IntArray, value: V): Bin<T>
public fun produce(builder: HistogramBuilder<T>.() -> Unit): IndexedHistogram<T, V>
override fun add(a: IndexedHistogram<T, V>, b: IndexedHistogram<T, V>): IndexedHistogram<T, V> {
require(a.context == this) { "Can't operate on a histogram produced by external space" }
require(b.context == this) { "Can't operate on a histogram produced by external space" }
return IndexedHistogram(this, histogramValueSpace.invoke { a.values + b.values })
}
override fun multiply(a: IndexedHistogram<T, V>, k: Number): IndexedHistogram<T, V> {
require(a.context == this) { "Can't operate on a histogram produced by external space" }
return IndexedHistogram(this, histogramValueSpace.invoke { a.values * k })
}
override val zero: IndexedHistogram<T, V> get() = produce { }
}

View File

@ -1,165 +0,0 @@
package kscience.kmath.histogram
import kscience.kmath.linear.Point
import kscience.kmath.nd.DefaultStrides
import kscience.kmath.nd.NDStructure
import kscience.kmath.operations.SpaceOperations
import kscience.kmath.operations.invoke
import kscience.kmath.structures.*
import kotlin.math.floor
public data class MultivariateBinDefinition<T : Comparable<T>>(
public val space: SpaceOperations<Point<T>>,
public val center: Point<T>,
public val sizes: Point<T>,
) {
public fun contains(vector: Point<out T>): Boolean {
require(vector.size == center.size) { "Dimension mismatch for input vector. Expected ${center.size}, but found ${vector.size}" }
val upper = space { center + sizes / 2.0 }
val lower = space { center - sizes / 2.0 }
return vector.asSequence().mapIndexed { i, value -> value in lower[i]..upper[i] }.all { it }
}
}
public class MultivariateBin<T : Comparable<T>>(
public val definition: MultivariateBinDefinition<T>,
public val count: Long,
public override val value: Double,
) : Bin<T> {
public override val dimension: Int
get() = definition.center.size
public override val center: Point<T>
get() = definition.center
public override operator fun contains(point: Point<T>): Boolean = definition.contains(point)
}
/**
* Uniform multivariate histogram with fixed borders. Based on NDStructure implementation with complexity of m for bin search, where m is the number of dimensions.
*/
public class RealHistogram(
private val lower: Buffer<Double>,
private val upper: Buffer<Double>,
private val binNums: IntArray = IntArray(lower.size) { 20 },
) : MutableHistogram<Double, MultivariateBin<Double>> {
private val strides = DefaultStrides(IntArray(binNums.size) { binNums[it] + 2 })
private val counts: NDStructure<LongCounter> = NDStructure.auto(strides) { LongCounter() }
private val values: NDStructure<DoubleCounter> = NDStructure.auto(strides) { DoubleCounter() }
public override val dimension: Int get() = lower.size
private val binSize = RealBuffer(dimension) { (upper[it] - lower[it]) / binNums[it] }
init {
// argument checks
require(lower.size == upper.size) { "Dimension mismatch in histogram lower and upper limits." }
require(lower.size == binNums.size) { "Dimension mismatch in bin count." }
require(!(0 until dimension).any { upper[it] - lower[it] < 0 }) { "Range for one of axis is not strictly positive" }
}
/**
* Get internal [NDStructure] bin index for given axis
*/
private fun getIndex(axis: Int, value: Double): Int = when {
value >= upper[axis] -> binNums[axis] + 1 // overflow
value < lower[axis] -> 0 // underflow
else -> floor((value - lower[axis]) / binSize[axis]).toInt() + 1
}
private fun getIndex(point: Buffer<out Double>): IntArray = IntArray(dimension) { getIndex(it, point[it]) }
private fun getCount(index: IntArray): Long = counts[index].sum()
public fun getCount(point: Buffer<out Double>): Long = getCount(getIndex(point))
private fun getValue(index: IntArray): Double = values[index].sum()
public fun getValue(point: Buffer<out Double>): Double = getValue(getIndex(point))
private fun getBinDefinition(index: IntArray): MultivariateBinDefinition<Double> {
val center = index.mapIndexed { axis, i ->
when (i) {
0 -> Double.NEGATIVE_INFINITY
strides.shape[axis] - 1 -> Double.POSITIVE_INFINITY
else -> lower[axis] + (i.toDouble() - 0.5) * binSize[axis]
}
}.asBuffer()
return MultivariateBinDefinition(RealBufferFieldOperations, center, binSize)
}
public fun getBinDefinition(point: Buffer<out Double>): MultivariateBinDefinition<Double> = getBinDefinition(getIndex(point))
public override operator fun get(point: Buffer<out Double>): MultivariateBin<Double>? {
val index = getIndex(point)
return MultivariateBin(getBinDefinition(index), getCount(index),getValue(index))
}
// fun put(point: Point<out Double>){
// val index = getIndex(point)
// values[index].increment()
// }
public override fun putWithWeight(point: Buffer<out Double>, weight: Double) {
val index = getIndex(point)
counts[index].increment()
values[index].add(weight)
}
public override operator fun iterator(): Iterator<MultivariateBin<Double>> =
strides.indices().map { index->
MultivariateBin(getBinDefinition(index), counts[index].sum(), values[index].sum())
}.iterator()
/**
* NDStructure containing number of events in bins without weights
*/
public fun counts(): NDStructure<Long> = NDStructure.auto(counts.shape) { counts[it].sum() }
/**
* NDStructure containing values of bins including weights
*/
public fun values(): NDStructure<Double> = NDStructure.auto(values.shape) { values[it].sum() }
public companion object {
/**
* Use it like
* ```
*FastHistogram.fromRanges(
* (-1.0..1.0),
* (-1.0..1.0)
*)
*```
*/
public fun fromRanges(vararg ranges: ClosedFloatingPointRange<Double>): RealHistogram = RealHistogram(
ranges.map(ClosedFloatingPointRange<Double>::start).asBuffer(),
ranges.map(ClosedFloatingPointRange<Double>::endInclusive).asBuffer()
)
/**
* Use it like
* ```
*FastHistogram.fromRanges(
* (-1.0..1.0) to 50,
* (-1.0..1.0) to 32
*)
*```
*/
public fun fromRanges(vararg ranges: Pair<ClosedFloatingPointRange<Double>, Int>): RealHistogram =
RealHistogram(
ListBuffer(
ranges
.map(Pair<ClosedFloatingPointRange<Double>, Int>::first)
.map(ClosedFloatingPointRange<Double>::start)
),
ListBuffer(
ranges
.map(Pair<ClosedFloatingPointRange<Double>, Int>::first)
.map(ClosedFloatingPointRange<Double>::endInclusive)
),
ranges.map(Pair<ClosedFloatingPointRange<Double>, Int>::second).toIntArray()
)
}
}

View File

@ -0,0 +1,121 @@
package kscience.kmath.histogram
import kscience.kmath.domains.Domain
import kscience.kmath.domains.HyperSquareDomain
import kscience.kmath.nd.*
import kscience.kmath.structures.*
import kotlin.math.floor
public class RealHistogramSpace(
private val lower: Buffer<Double>,
private val upper: Buffer<Double>,
private val binNums: IntArray = IntArray(lower.size) { 20 },
) : IndexedHistogramSpace<Double, Double> {
init {
// argument checks
require(lower.size == upper.size) { "Dimension mismatch in histogram lower and upper limits." }
require(lower.size == binNums.size) { "Dimension mismatch in bin count." }
require(!lower.indices.any { upper[it] - lower[it] < 0 }) { "Range for one of axis is not strictly positive" }
}
public val dimension: Int get() = lower.size
private val shape = IntArray(binNums.size) { binNums[it] + 2 }
override val histogramValueSpace: RealNDField = NDAlgebra.real(*shape)
override val strides: Strides get() = histogramValueSpace.strides
private val binSize = RealBuffer(dimension) { (upper[it] - lower[it]) / binNums[it] }
/**
* Get internal [NDStructure] bin index for given axis
*/
private fun getIndex(axis: Int, value: Double): Int = when {
value >= upper[axis] -> binNums[axis] + 1 // overflow
value < lower[axis] -> 0 // underflow
else -> floor((value - lower[axis]) / binSize[axis]).toInt()
}
override fun getIndex(point: Buffer<Double>): IntArray = IntArray(dimension) {
getIndex(it, point[it])
}
override fun getDomain(index: IntArray): Domain<Double> {
val lowerBoundary = index.mapIndexed { axis, i ->
when (i) {
0 -> Double.NEGATIVE_INFINITY
strides.shape[axis] - 1 -> upper[axis]
else -> lower[axis] + (i.toDouble()) * binSize[axis]
}
}.asBuffer()
val upperBoundary = index.mapIndexed { axis, i ->
when (i) {
0 -> lower[axis]
strides.shape[axis] - 1 -> Double.POSITIVE_INFINITY
else -> lower[axis] + (i.toDouble() + 1) * binSize[axis]
}
}.asBuffer()
return HyperSquareDomain(lowerBoundary, upperBoundary)
}
override fun produceBin(index: IntArray, value: Double): Bin<Double> {
val domain = getDomain(index)
return DomainBin(domain, value)
}
override fun produce(builder: HistogramBuilder<Double>.() -> Unit): IndexedHistogram<Double, Double> {
val ndCounter = NDStructure.auto(strides) { Counter.real() }
val hBuilder = HistogramBuilder<Double> { point, value ->
val index = getIndex(point)
ndCounter[index].add(1.0)
}
hBuilder.apply(builder)
val values: NDBuffer<Double> = ndCounter.mapToBuffer { it.value }
return IndexedHistogram(this, values)
}
public companion object {
/**
* Use it like
* ```
*FastHistogram.fromRanges(
* (-1.0..1.0),
* (-1.0..1.0)
*)
*```
*/
public fun fromRanges(vararg ranges: ClosedFloatingPointRange<Double>): RealHistogramSpace = RealHistogramSpace(
ranges.map(ClosedFloatingPointRange<Double>::start).asBuffer(),
ranges.map(ClosedFloatingPointRange<Double>::endInclusive).asBuffer()
)
/**
* Use it like
* ```
*FastHistogram.fromRanges(
* (-1.0..1.0) to 50,
* (-1.0..1.0) to 32
*)
*```
*/
public fun fromRanges(vararg ranges: Pair<ClosedFloatingPointRange<Double>, Int>): RealHistogramSpace =
RealHistogramSpace(
ListBuffer(
ranges
.map(Pair<ClosedFloatingPointRange<Double>, Int>::first)
.map(ClosedFloatingPointRange<Double>::start)
),
ListBuffer(
ranges
.map(Pair<ClosedFloatingPointRange<Double>, Int>::first)
.map(ClosedFloatingPointRange<Double>::endInclusive)
),
ranges.map(Pair<ClosedFloatingPointRange<Double>, Int>::second).toIntArray()
)
}
}

View File

@ -1,8 +1,8 @@
package scietifik.kmath.histogram
import kscience.kmath.histogram.RealHistogram
import kscience.kmath.histogram.fill
import kscience.kmath.histogram.RealHistogramSpace
import kscience.kmath.histogram.put
import kscience.kmath.operations.invoke
import kscience.kmath.real.RealVector
import kscience.kmath.real.invoke
import kotlin.random.Random
@ -11,12 +11,14 @@ import kotlin.test.*
internal class MultivariateHistogramTest {
@Test
fun testSinglePutHistogram() {
val histogram = RealHistogram.fromRanges(
val hSpace = RealHistogramSpace.fromRanges(
(-1.0..1.0),
(-1.0..1.0)
)
histogram.put(0.55, 0.55)
val bin = histogram.find { it.value.toInt() > 0 } ?: fail()
val histogram = hSpace.produce {
put(0.55, 0.55)
}
val bin = histogram.bins.find { it.value.toInt() > 0 } ?: fail()
assertTrue { bin.contains(RealVector(0.55, 0.55)) }
assertTrue { bin.contains(RealVector(0.6, 0.5)) }
assertFalse { bin.contains(RealVector(-0.55, 0.55)) }
@ -24,7 +26,7 @@ internal class MultivariateHistogramTest {
@Test
fun testSequentialPut() {
val histogram = RealHistogram.fromRanges(
val hSpace = RealHistogramSpace.fromRanges(
(-1.0..1.0),
(-1.0..1.0),
(-1.0..1.0)
@ -34,12 +36,45 @@ internal class MultivariateHistogramTest {
fun nextDouble() = random.nextDouble(-1.0, 1.0)
val n = 10000
histogram.fill {
val histogram = hSpace.produce {
repeat(n) {
yield(RealVector(nextDouble(), nextDouble(), nextDouble()))
put(nextDouble(), nextDouble(), nextDouble())
}
}
assertEquals(n, histogram.sumBy { it.value.toInt() })
assertEquals(n, histogram.bins.sumBy { it.value.toInt() })
}
@Test
fun testHistogramAlgebra() {
val hSpace = RealHistogramSpace.fromRanges(
(-1.0..1.0),
(-1.0..1.0),
(-1.0..1.0)
).invoke {
val random = Random(1234)
fun nextDouble() = random.nextDouble(-1.0, 1.0)
val n = 10000
val histogram1 = produce {
repeat(n) {
put(nextDouble(), nextDouble(), nextDouble())
}
}
val histogram2 = produce {
repeat(n) {
put(nextDouble(), nextDouble(), nextDouble())
}
}
val res = histogram1 - histogram2
assertTrue {
strides.indices().all { index ->
res.values[index] <= histogram1.values[index]
}
}
assertTrue {
res.bins.count() >= histogram1.bins.count()
}
assertEquals(0.0, res.bins.sumByDouble { it.value.toDouble() })
}
}
}

View File

@ -1,37 +0,0 @@
package kscience.kmath.histogram
public actual class LongCounter {
private var sum: Long = 0L
public actual fun decrement() {
sum--
}
public actual fun increment() {
sum++
}
public actual fun reset() {
sum = 0
}
public actual fun sum(): Long = sum
public actual fun add(l: Long) {
sum += l
}
}
public actual class DoubleCounter {
private var sum: Double = 0.0
public actual fun reset() {
sum = 0.0
}
public actual fun sum(): Double = sum
public actual fun add(d: Double) {
sum += d
}
}

View File

@ -1,7 +0,0 @@
package kscience.kmath.histogram
import java.util.concurrent.atomic.DoubleAdder
import java.util.concurrent.atomic.LongAdder
public actual typealias LongCounter = LongAdder
public actual typealias DoubleCounter = DoubleAdder

View File

@ -0,0 +1,155 @@
package kscience.kmath.histogram
import kscience.kmath.domains.UnivariateDomain
import kscience.kmath.misc.UnstableKMathAPI
import kscience.kmath.operations.Space
import kscience.kmath.structures.Buffer
import java.util.*
import kotlin.math.abs
import kotlin.math.sqrt
private fun <B : ClosedFloatingPointRange<Double>> TreeMap<Double, B>.getBin(value: Double): B? {
// check ceiling entry and return it if it is what needed
val ceil = ceilingEntry(value)?.value
if (ceil != null && value in ceil) return ceil
//check floor entry
val floor = floorEntry(value)?.value
if (floor != null && value in floor) return floor
//neither is valid, not found
return null
}
@UnstableKMathAPI
public class TreeHistogram(
override val context: TreeHistogramSpace,
private val binMap: TreeMap<Double, out UnivariateBin>,
) : UnivariateHistogram {
override fun get(value: Double): UnivariateBin? = binMap.getBin(value)
override val dimension: Int get() = 1
override val bins: Collection<UnivariateBin> get() = binMap.values
}
/**
* A space for univariate histograms with variable bin borders based on a tree map
*/
@UnstableKMathAPI
public class TreeHistogramSpace(
public val binFactory: (Double) -> UnivariateDomain,
) : Space<UnivariateHistogram> {
private class BinCounter(val domain: UnivariateDomain, val counter: Counter<Double> = Counter.real()) :
ClosedFloatingPointRange<Double> by domain.range
public fun produce(builder: UnivariateHistogramBuilder.() -> Unit): UnivariateHistogram {
val bins: TreeMap<Double, BinCounter> = TreeMap()
val hBuilder = object : UnivariateHistogramBuilder {
fun get(value: Double): BinCounter? = bins.getBin(value)
fun createBin(value: Double): BinCounter {
val binDefinition = binFactory(value)
val newBin = BinCounter(binDefinition)
synchronized(this) { bins[binDefinition.center] = newBin }
return newBin
}
/**
* Thread safe put operation
*/
override fun putValue(at: Double, value: Double) {
(get(at) ?: createBin(at)).apply {
counter.add(value)
}
}
override fun putValue(point: Buffer<Double>, value: Number) {
put(point[0], value.toDouble())
}
}
hBuilder.apply(builder)
val resBins = TreeMap<Double, UnivariateBin>()
bins.forEach { key, binCounter ->
val count = binCounter.counter.value
resBins[key] = UnivariateBin(binCounter.domain, count, sqrt(count))
}
return TreeHistogram(this, resBins)
}
override fun add(
a: UnivariateHistogram,
b: UnivariateHistogram,
): UnivariateHistogram {
require(a.context == this) { "Histogram $a does not belong to this context" }
require(b.context == this) { "Histogram $b does not belong to this context" }
val bins = TreeMap<Double, UnivariateBin>().apply {
(a.bins.map { it.domain } union b.bins.map { it.domain }).forEach { def ->
put(def.center,
UnivariateBin(
def,
value = (a[def.center]?.value ?: 0.0) + (b[def.center]?.value ?: 0.0),
standardDeviation = (a[def.center]?.standardDeviation
?: 0.0) + (b[def.center]?.standardDeviation ?: 0.0)
)
)
}
}
return TreeHistogram(this, bins)
}
override fun multiply(a: UnivariateHistogram, k: Number): UnivariateHistogram {
val bins = TreeMap<Double, UnivariateBin>().apply {
a.bins.forEach { bin ->
put(bin.domain.center,
UnivariateBin(
bin.domain,
value = bin.value * k.toDouble(),
standardDeviation = abs(bin.standardDeviation * k.toDouble())
)
)
}
}
return TreeHistogram(this, bins)
}
override val zero: UnivariateHistogram = produce { }
public companion object {
/**
* Build and fill a [UnivariateHistogram]. Returns a read-only histogram.
*/
public fun uniform(
binSize: Double,
start: Double = 0.0,
): TreeHistogramSpace = TreeHistogramSpace { value ->
val center = start + binSize * Math.floor((value - start) / binSize + 0.5)
UnivariateDomain((center - binSize / 2)..(center + binSize / 2))
}
/**
* Create a histogram with custom cell borders
*/
public fun custom(borders: DoubleArray): TreeHistogramSpace {
val sorted = borders.sortedArray()
return TreeHistogramSpace { value ->
when {
value < sorted.first() -> UnivariateDomain(
Double.NEGATIVE_INFINITY..sorted.first()
)
value > sorted.last() -> UnivariateDomain(
sorted.last()..Double.POSITIVE_INFINITY
)
else -> {
val index = sorted.indices.first { value > sorted[it] }
val left = sorted[index]
val right = sorted[index + 1]
UnivariateDomain(left..right)
}
}
}
}
}
}

View File

@ -1,76 +1,38 @@
package kscience.kmath.histogram
import kscience.kmath.linear.Point
import kscience.kmath.domains.UnivariateDomain
import kscience.kmath.misc.UnstableKMathAPI
import kscience.kmath.operations.Space
import kscience.kmath.operations.SpaceElement
import kscience.kmath.structures.Buffer
import kscience.kmath.structures.asBuffer
import kscience.kmath.structures.asSequence
import java.util.*
import kotlin.math.floor
//TODO move to common
public val UnivariateDomain.center: Double get() = (range.endInclusive - range.start) / 2
/**
* A univariate bin based an a range
* @param value The value of histogram including weighting
* @param standardDeviation Standard deviation of the bin value. Zero or negative if not applicable
*/
public class UnivariateBin(
public val position: Double,
public val size: Double,
) : Bin<Double> {
//internal mutation operations
internal val counter: LongCounter = LongCounter()
internal val weightCounter: DoubleCounter = DoubleCounter()
public val domain: UnivariateDomain,
override val value: Double,
public val standardDeviation: Double,
) : Bin<Double>, ClosedFloatingPointRange<Double> by domain.range {
/**
* The precise number of events ignoring weighting
*/
public val count: Long get() = counter.sum()
/**
* The value of histogram including weighting
*/
public override val value: Double get() = weightCounter.sum()
public override val center: Point<Double> get() = doubleArrayOf(position).asBuffer()
public override val dimension: Int get() = 1
public operator fun contains(value: Double): Boolean = value in (position - size / 2)..(position + size / 2)
public override fun contains(point: Buffer<Double>): Boolean = contains(point[0])
public override fun contains(point: Buffer<Double>): Boolean = point.size == 1 && contains(point[0])
}
/**
* Univariate histogram with log(n) bin search speed
*/
@OptIn(UnstableKMathAPI::class)
public abstract class UnivariateHistogram protected constructor(
protected val bins: TreeMap<Double, UnivariateBin> = TreeMap(),
) : Histogram<Double, UnivariateBin>, SpaceElement<UnivariateHistogram, UnivariateHistogramSpace> {
public operator fun get(value: Double): UnivariateBin? {
// check ceiling entry and return it if it is what needed
val ceil = bins.ceilingEntry(value)?.value
if (ceil != null && value in ceil) return ceil
//check floor entry
val floor = bins.floorEntry(value)?.value
if (floor != null && value in floor) return floor
//neither is valid, not found
return null
}
public override operator fun get(point: Buffer<out Double>): UnivariateBin? = get(point[0])
public override val dimension: Int get() = 1
public override operator fun iterator(): Iterator<UnivariateBin> = bins.values.iterator()
public interface UnivariateHistogram : Histogram<Double, UnivariateBin>,
SpaceElement<UnivariateHistogram, Space<UnivariateHistogram>> {
public operator fun get(value: Double): UnivariateBin?
public override operator fun get(point: Buffer<Double>): UnivariateBin? = get(point[0])
public companion object {
/**
* Build a histogram with a uniform binning with a start at [start] and a bin size of [binSize]
*/
public fun uniformBuilder(binSize: Double, start: Double = 0.0): UnivariateHistogramBuilder =
UnivariateHistogramSpace { value ->
val center = start + binSize * floor((value - start) / binSize + 0.5)
UnivariateBin(center, binSize)
}.builder()
/**
* Build and fill a [UnivariateHistogram]. Returns a read-only histogram.
*/
@ -78,35 +40,7 @@ public abstract class UnivariateHistogram protected constructor(
binSize: Double,
start: Double = 0.0,
builder: UnivariateHistogramBuilder.() -> Unit,
): UnivariateHistogram = uniformBuilder(binSize, start).apply(builder)
/**
* Create a histogram with custom cell borders
*/
public fun customBuilder(borders: DoubleArray): UnivariateHistogramBuilder {
val sorted = borders.sortedArray()
return UnivariateHistogramSpace { value ->
when {
value < sorted.first() -> UnivariateBin(
Double.NEGATIVE_INFINITY,
Double.MAX_VALUE
)
value > sorted.last() -> UnivariateBin(
Double.POSITIVE_INFINITY,
Double.MAX_VALUE
)
else -> {
val index = sorted.indices.first { value > sorted[it] }
val left = sorted[index]
val right = sorted[index + 1]
UnivariateBin((left + right) / 2, (right - left))
}
}
}.builder()
}
): UnivariateHistogram = TreeHistogramSpace.uniform(binSize, start).produce(builder)
/**
* Build and fill a histogram with custom borders. Returns a read-only histogram.
@ -114,48 +48,26 @@ public abstract class UnivariateHistogram protected constructor(
public fun custom(
borders: DoubleArray,
builder: UnivariateHistogramBuilder.() -> Unit,
): UnivariateHistogram = customBuilder(borders).apply(builder)
): UnivariateHistogram = TreeHistogramSpace.custom(borders).produce(builder)
}
}
public class UnivariateHistogramBuilder internal constructor(
override val context: UnivariateHistogramSpace,
) : UnivariateHistogram(), MutableHistogram<Double, UnivariateBin> {
private fun createBin(value: Double): UnivariateBin = context.binFactory(value).also {
synchronized(this) { bins[it.position] = it }
}
@UnstableKMathAPI
public interface UnivariateHistogramBuilder : HistogramBuilder<Double> {
/**
* Thread safe put operation
*/
public fun put(value: Double, weight: Double = 1.0) {
(get(value) ?: createBin(value)).apply {
counter.increment()
weightCounter.add(weight)
}
}
public fun putValue(at: Double, value: Double = 1.0)
override fun putWithWeight(point: Buffer<out Double>, weight: Double) {
put(point[0], weight)
}
/**
* Put several items into a single bin
*/
public fun putMany(value: Double, count: Int, weight: Double = count.toDouble()) {
(get(value) ?: createBin(value)).apply {
counter.add(count.toLong())
weightCounter.add(weight)
}
}
override fun putValue(point: Buffer<Double>, value: Number)
}
@UnstableKMathAPI
public fun UnivariateHistogramBuilder.fill(items: Iterable<Double>): Unit = items.forEach(::put)
public fun UnivariateHistogramBuilder.fill(items: Iterable<Double>): Unit = items.forEach(this::putValue)
@UnstableKMathAPI
public fun UnivariateHistogramBuilder.fill(array: DoubleArray): Unit = array.forEach(::put)
public fun UnivariateHistogramBuilder.fill(array: DoubleArray): Unit = array.forEach(this::putValue)
@UnstableKMathAPI
public fun UnivariateHistogramBuilder.fill(buffer: Buffer<Double>): Unit = buffer.asSequence().forEach(::put)
public fun UnivariateHistogramBuilder.fill(buffer: Buffer<Double>): Unit = buffer.asSequence().forEach(this::putValue)

View File

@ -1,25 +0,0 @@
package kscience.kmath.histogram
import kscience.kmath.operations.Space
public class UnivariateHistogramSpace(public val binFactory: (Double) -> UnivariateBin) : Space<UnivariateHistogram> {
public fun builder(): UnivariateHistogramBuilder = UnivariateHistogramBuilder(this)
public fun produce(builder: UnivariateHistogramBuilder.() -> Unit): UnivariateHistogram = builder().apply(builder)
override fun add(
a: UnivariateHistogram,
b: UnivariateHistogram,
): UnivariateHistogram {
require(a.context == this){"Histogram $a does not belong to this context"}
require(b.context == this){"Histogram $b does not belong to this context"}
TODO()
}
override fun multiply(a: UnivariateHistogram, k: Number): UnivariateHistogram {
TODO("Not yet implemented")
}
override val zero: UnivariateHistogram = produce { }
}