diff --git a/CHANGELOG.md b/CHANGELOG.md index a19b1f467..bc5447449 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -48,6 +48,8 @@ - Operations -> Ops - Default Buffer and ND algebras are now Ops and lack neutral elements (0, 1) as well as algebra-level shapes. - Tensor algebra takes read-only structures as input and inherits AlgebraND +- `UnivariateDistribution` renamed to `Distribution1D` +- Rework of histograms. ### Deprecated - Specialized `DoubleBufferAlgebra` diff --git a/build.gradle.kts b/build.gradle.kts index 445976853..1e703c5e1 100644 --- a/build.gradle.kts +++ b/build.gradle.kts @@ -11,7 +11,7 @@ allprojects { } group = "space.kscience" - version = "0.3.0-dev-20" + version = "0.3.0-dev-21" } subprojects { diff --git a/gradle.properties b/gradle.properties index a7cd2f876..80108737e 100644 --- a/gradle.properties +++ b/gradle.properties @@ -6,6 +6,8 @@ kotlin.code.style=official kotlin.jupyter.add.scanner=false kotlin.mpp.stability.nowarn=true kotlin.native.ignoreDisabledTargets=true +#kotlin.incremental.js.ir=true + org.gradle.configureondemand=true org.gradle.jvmargs=-XX:MaxMetaspaceSize=1G org.gradle.parallel=true diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/domains/Domain1D.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/domains/Domain1D.kt new file mode 100644 index 000000000..3d531349c --- /dev/null +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/domains/Domain1D.kt @@ -0,0 +1,59 @@ +/* + * Copyright 2018-2021 KMath contributors. + * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file. + */ + +package space.kscience.kmath.domains + +import space.kscience.kmath.linear.Point +import space.kscience.kmath.misc.UnstableKMathAPI + +@UnstableKMathAPI +public abstract class Domain1D>(public val range: ClosedRange) : Domain { + override val dimension: Int get() = 1 + + public operator fun contains(value: T): Boolean = range.contains(value) + + override operator fun contains(point: Point): Boolean { + require(point.size == 0) + return contains(point[0]) + } +} + +@UnstableKMathAPI +public class DoubleDomain1D( + @Suppress("CanBeParameter") public val doubleRange: ClosedFloatingPointRange, +) : Domain1D(doubleRange), DoubleDomain { + override fun getLowerBound(num: Int): Double { + require(num == 0) + return range.start + } + + override fun getUpperBound(num: Int): Double { + require(num == 0) + return range.endInclusive + } + + override fun volume(): Double = range.endInclusive - range.start + + override fun equals(other: Any?): Boolean { + if (this === other) return true + if (other == null || this::class != other::class) return false + + other as DoubleDomain1D + + if (doubleRange != other.doubleRange) return false + + return true + } + + override fun hashCode(): Int = doubleRange.hashCode() + + override fun toString(): String = doubleRange.toString() + + +} + +@UnstableKMathAPI +public val Domain1D.center: Double + get() = (range.endInclusive + range.start) / 2 diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/domains/HyperSquareDomain.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/domains/HyperSquareDomain.kt index bd5514623..485416a69 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/domains/HyperSquareDomain.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/domains/HyperSquareDomain.kt @@ -7,18 +7,28 @@ package space.kscience.kmath.domains import space.kscience.kmath.linear.Point import space.kscience.kmath.misc.UnstableKMathAPI import space.kscience.kmath.structures.Buffer +import space.kscience.kmath.structures.DoubleBuffer import space.kscience.kmath.structures.indices /** - * - * HyperSquareDomain class. - * - * @author Alexander Nozik + * A hyper-square (or hyper-cube) real-space domain. It is formed by a [Buffer] of [lower] boundaries + * and a [Buffer] of upper boundaries. Upper should be greater or equals than lower. */ @UnstableKMathAPI -public class HyperSquareDomain(private val lower: Buffer, private val upper: Buffer) : DoubleDomain { +public class HyperSquareDomain(public val lower: Buffer, public val upper: Buffer) : DoubleDomain { + init { + require(lower.size == upper.size) { + "Domain borders size mismatch. Lower borders size is ${lower.size}, but upper borders size is ${upper.size}." + } + require(lower.indices.all { lower[it] <= upper[it] }) { + "Domain borders order mismatch. Lower borders must be less or equals than upper borders." + } + } + override val dimension: Int get() = lower.size + public val center: DoubleBuffer get() = DoubleBuffer(dimension) { (lower[it] + upper[it]) / 2.0 } + override operator fun contains(point: Point): Boolean = point.indices.all { i -> point[i] in lower[i]..upper[i] } diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/domains/UnivariateDomain.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/domains/UnivariateDomain.kt deleted file mode 100644 index 9020ef8cb..000000000 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/domains/UnivariateDomain.kt +++ /dev/null @@ -1,33 +0,0 @@ -/* - * Copyright 2018-2021 KMath contributors. - * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file. - */ - -package space.kscience.kmath.domains - -import space.kscience.kmath.linear.Point -import space.kscience.kmath.misc.UnstableKMathAPI - -@UnstableKMathAPI -public class UnivariateDomain(public val range: ClosedFloatingPointRange) : DoubleDomain { - override val dimension: Int get() = 1 - - public operator fun contains(d: Double): Boolean = range.contains(d) - - override operator fun contains(point: Point): Boolean { - require(point.size == 0) - return contains(point[0]) - } - - override fun getLowerBound(num: Int): Double { - require(num == 0) - return range.start - } - - override fun getUpperBound(num: Int): Double { - require(num == 0) - return range.endInclusive - } - - override fun volume(): Double = range.endInclusive - range.start -} diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/misc/sorting.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/misc/sorting.kt index a144e49b4..dc5421136 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/misc/sorting.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/misc/sorting.kt @@ -3,43 +3,71 @@ * Use of this source code is governed by the Apache 2.0 license that can be found in the LICENSE file. */ + package space.kscience.kmath.misc -import kotlin.comparisons.* import space.kscience.kmath.structures.Buffer +import space.kscience.kmath.structures.VirtualBuffer /** - * Return a new list filled with buffer indices. Indice order is defined by sorting associated buffer value. - * This feature allows to sort buffer values without reordering its content. + * Return a new array filled with buffer indices. Indices order is defined by sorting associated buffer value. + * This feature allows sorting buffer values without reordering its content. * - * @return List of buffer indices, sorted by associated value. + * @return Buffer indices, sorted by associated value. */ -@PerformancePitfall @UnstableKMathAPI -public fun > Buffer.permSort() : IntArray = _permSortWith(compareBy { get(it) }) +public fun > Buffer.indicesSorted(): IntArray = permSortIndicesWith(compareBy { get(it) }) + +/** + * Create a zero-copy virtual buffer that contains the same elements but in ascending order + */ +@OptIn(UnstableKMathAPI::class) +public fun > Buffer.sorted(): Buffer { + val permutations = indicesSorted() + return VirtualBuffer(size) { this[permutations[it]] } +} -@PerformancePitfall @UnstableKMathAPI -public fun > Buffer.permSortDescending() : IntArray = _permSortWith(compareByDescending { get(it) }) +public fun > Buffer.indicesSortedDescending(): IntArray = + permSortIndicesWith(compareByDescending { get(it) }) + +/** + * Create a zero-copy virtual buffer that contains the same elements but in descending order + */ +@OptIn(UnstableKMathAPI::class) +public fun > Buffer.sortedDescending(): Buffer { + val permutations = indicesSortedDescending() + return VirtualBuffer(size) { this[permutations[it]] } +} -@PerformancePitfall @UnstableKMathAPI -public fun > Buffer.permSortBy(selector: (V) -> C) : IntArray = _permSortWith(compareBy { selector(get(it)) }) +public fun > Buffer.indicesSortedBy(selector: (V) -> C): IntArray = + permSortIndicesWith(compareBy { selector(get(it)) }) + +@OptIn(UnstableKMathAPI::class) +public fun > Buffer.sortedBy(selector: (V) -> C): Buffer { + val permutations = indicesSortedBy(selector) + return VirtualBuffer(size) { this[permutations[it]] } +} -@PerformancePitfall @UnstableKMathAPI -public fun > Buffer.permSortByDescending(selector: (V) -> C) : IntArray = _permSortWith(compareByDescending { selector(get(it)) }) +public fun > Buffer.indicesSortedByDescending(selector: (V) -> C): IntArray = + permSortIndicesWith(compareByDescending { selector(get(it)) }) + +@OptIn(UnstableKMathAPI::class) +public fun > Buffer.sortedByDescending(selector: (V) -> C): Buffer { + val permutations = indicesSortedByDescending(selector) + return VirtualBuffer(size) { this[permutations[it]] } +} -@PerformancePitfall @UnstableKMathAPI -public fun Buffer.permSortWith(comparator : Comparator) : IntArray = _permSortWith { i1, i2 -> comparator.compare(get(i1), get(i2)) } +public fun Buffer.indicesSortedWith(comparator: Comparator): IntArray = + permSortIndicesWith { i1, i2 -> comparator.compare(get(i1), get(i2)) } -@PerformancePitfall -@UnstableKMathAPI -private fun Buffer._permSortWith(comparator : Comparator) : IntArray { - if (size < 2) return IntArray(size) +private fun Buffer.permSortIndicesWith(comparator: Comparator): IntArray { + if (size < 2) return IntArray(size) { 0 } - /* TODO: optimisation : keep a constant big array of indices (Ex: from 0 to 4096), then create indice + /* TODO: optimisation : keep a constant big array of indices (Ex: from 0 to 4096), then create indices * arrays more efficiently by copying subpart of cached one. For bigger needs, we could copy entire * cached array, then fill remaining indices manually. Not done for now, because: * 1. doing it right would require some statistics about common used buffer sizes. @@ -53,3 +81,12 @@ private fun Buffer._permSortWith(comparator : Comparator) : IntArray */ return packedIndices.sortedWith(comparator).toIntArray() } + +/** + * Checks that the [Buffer] is sorted (ascending) and throws [IllegalArgumentException] if it is not. + */ +public fun > Buffer.requireSorted() { + for (i in 0..(size - 2)) { + require(get(i + 1) >= get(i)) { "The buffer is not sorted at index $i" } + } +} diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/AlgebraND.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/AlgebraND.kt index a071c1eb3..a9712e870 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/AlgebraND.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/AlgebraND.kt @@ -194,7 +194,7 @@ public interface RingOpsND> : RingOps>, Gro override fun multiply(left: StructureND, right: StructureND): StructureND = zip(left, right) { aValue, bValue -> multiply(aValue, bValue) } - //TODO move to extensions after KEEP-176 + //TODO move to extensions with context receivers /** * Multiplies an ND structure by an element of it. diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/BufferND.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/BufferND.kt index 539499794..2401f6319 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/BufferND.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/BufferND.kt @@ -32,18 +32,23 @@ public open class BufferND( /** * Transform structure to a new structure using provided [BufferFactory] and optimizing if argument is [BufferND] */ -public inline fun StructureND.mapToBuffer( - factory: BufferFactory = Buffer.Companion::auto, +public inline fun StructureND.mapToBuffer( + factory: BufferFactory, crossinline transform: (T) -> R, -): BufferND { - return if (this is BufferND) - BufferND(this.indices, factory.invoke(indices.linearSize) { transform(buffer[it]) }) - else { - val strides = DefaultStrides(shape) - BufferND(strides, factory.invoke(strides.linearSize) { transform(get(strides.index(it))) }) - } +): BufferND = if (this is BufferND) + BufferND(this.indices, factory.invoke(indices.linearSize) { transform(buffer[it]) }) +else { + val strides = DefaultStrides(shape) + BufferND(strides, factory.invoke(strides.linearSize) { transform(get(strides.index(it))) }) } +/** + * Transform structure to a new structure using inferred [BufferFactory] + */ +public inline fun StructureND.mapToBuffer( + crossinline transform: (T) -> R, +): BufferND = mapToBuffer(Buffer.Companion::auto, transform) + /** * Represents [MutableStructureND] over [MutableBuffer]. * diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/Buffer.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/Buffer.kt index 58c6d5ded..a1b0307c4 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/Buffer.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/Buffer.kt @@ -105,6 +105,16 @@ public interface Buffer { */ public val Buffer<*>.indices: IntRange get() = 0 until size +public fun Buffer.first(): T { + require(size > 0) { "Can't get the first element of empty buffer" } + return get(0) +} + +public fun Buffer.last(): T { + require(size > 0) { "Can't get the last element of empty buffer" } + return get(size - 1) +} + /** * Immutable wrapper for [MutableBuffer]. * diff --git a/kmath-core/src/commonTest/kotlin/space/kscience/kmath/misc/PermSortTest.kt b/kmath-core/src/commonTest/kotlin/space/kscience/kmath/misc/PermSortTest.kt index 0a2bb9138..4a724ac5f 100644 --- a/kmath-core/src/commonTest/kotlin/space/kscience/kmath/misc/PermSortTest.kt +++ b/kmath-core/src/commonTest/kotlin/space/kscience/kmath/misc/PermSortTest.kt @@ -6,14 +6,13 @@ package space.kscience.kmath.misc import space.kscience.kmath.misc.PermSortTest.Platform.* -import kotlin.random.Random -import kotlin.test.Test -import kotlin.test.assertEquals -import kotlin.test.assertTrue - import space.kscience.kmath.structures.IntBuffer import space.kscience.kmath.structures.asBuffer +import kotlin.random.Random +import kotlin.test.Test import kotlin.test.assertContentEquals +import kotlin.test.assertEquals +import kotlin.test.assertTrue class PermSortTest { @@ -29,9 +28,9 @@ class PermSortTest { @Test fun testOnEmptyBuffer() { val emptyBuffer = IntBuffer(0) {it} - var permutations = emptyBuffer.permSort() + var permutations = emptyBuffer.indicesSorted() assertTrue(permutations.isEmpty(), "permutation on an empty buffer should return an empty result") - permutations = emptyBuffer.permSortDescending() + permutations = emptyBuffer.indicesSortedDescending() assertTrue(permutations.isEmpty(), "permutation on an empty buffer should return an empty result") } @@ -47,25 +46,25 @@ class PermSortTest { @Test fun testPermSortBy() { - val permutations = platforms.permSortBy { it.name } + val permutations = platforms.indicesSortedBy { it.name } val expected = listOf(ANDROID, JS, JVM, NATIVE, WASM) assertContentEquals(expected, permutations.map { platforms[it] }, "Ascending PermSort by name") } @Test fun testPermSortByDescending() { - val permutations = platforms.permSortByDescending { it.name } + val permutations = platforms.indicesSortedByDescending { it.name } val expected = listOf(WASM, NATIVE, JVM, JS, ANDROID) assertContentEquals(expected, permutations.map { platforms[it] }, "Descending PermSort by name") } @Test fun testPermSortWith() { - var permutations = platforms.permSortWith { p1, p2 -> p1.name.length.compareTo(p2.name.length) } + var permutations = platforms.indicesSortedWith { p1, p2 -> p1.name.length.compareTo(p2.name.length) } val expected = listOf(JS, JVM, WASM, NATIVE, ANDROID) assertContentEquals(expected, permutations.map { platforms[it] }, "PermSort using custom ascending comparator") - permutations = platforms.permSortWith(compareByDescending { it.name.length }) + permutations = platforms.indicesSortedWith(compareByDescending { it.name.length }) assertContentEquals(expected.reversed(), permutations.map { platforms[it] }, "PermSort using custom descending comparator") } @@ -75,7 +74,7 @@ class PermSortTest { println("Test randomization seed: $seed") val buffer = Random(seed).buffer(bufferSize) - val indices = buffer.permSort() + val indices = buffer.indicesSorted() assertEquals(bufferSize, indices.size) // Ensure no doublon is present in indices @@ -87,7 +86,7 @@ class PermSortTest { assertTrue(current <= next, "Permutation indices not properly sorted") } - val descIndices = buffer.permSortDescending() + val descIndices = buffer.indicesSortedDescending() assertEquals(bufferSize, descIndices.size) // Ensure no doublon is present in indices assertEquals(descIndices.toSet().size, descIndices.size) diff --git a/kmath-for-real/build.gradle.kts b/kmath-for-real/build.gradle.kts index 4cccaef5c..18c2c50ad 100644 --- a/kmath-for-real/build.gradle.kts +++ b/kmath-for-real/build.gradle.kts @@ -21,19 +21,22 @@ readme { feature( id = "DoubleVector", - description = "Numpy-like operations for Buffers/Points", ref = "src/commonMain/kotlin/space/kscience/kmath/real/DoubleVector.kt" - ) + ){ + "Numpy-like operations for Buffers/Points" + } feature( id = "DoubleMatrix", - description = "Numpy-like operations for 2d real structures", ref = "src/commonMain/kotlin/space/kscience/kmath/real/DoubleMatrix.kt" - ) + ){ + "Numpy-like operations for 2d real structures" + } feature( id = "grids", - description = "Uniform grid generators", ref = "src/commonMain/kotlin/space/kscience/kmath/structures/grids.kt" - ) + ){ + "Uniform grid generators" + } } diff --git a/kmath-histograms/build.gradle.kts b/kmath-histograms/build.gradle.kts index 7e511faa0..51271bb0a 100644 --- a/kmath-histograms/build.gradle.kts +++ b/kmath-histograms/build.gradle.kts @@ -17,6 +17,8 @@ kotlin.sourceSets { commonTest { dependencies { implementation(project(":kmath-for-real")) + implementation(projects.kmath.kmathStat) + implementation("org.jetbrains.kotlinx:kotlinx-coroutines-test:1.6.0") } } } diff --git a/kmath-histograms/src/commonMain/kotlin/space/kscience/kmath/histogram/Counter.kt b/kmath-histograms/src/commonMain/kotlin/space/kscience/kmath/histogram/Counter.kt index 291284444..fe3278026 100644 --- a/kmath-histograms/src/commonMain/kotlin/space/kscience/kmath/histogram/Counter.kt +++ b/kmath-histograms/src/commonMain/kotlin/space/kscience/kmath/histogram/Counter.kt @@ -18,7 +18,8 @@ public interface Counter { public val value: T public companion object { - public fun double(): ObjectCounter = ObjectCounter(DoubleField) + public fun ofDouble(): ObjectCounter = ObjectCounter(DoubleField) + public fun of(group: Group): ObjectCounter = ObjectCounter(group) } } diff --git a/kmath-histograms/src/commonMain/kotlin/space/kscience/kmath/histogram/DoubleHistogramSpace.kt b/kmath-histograms/src/commonMain/kotlin/space/kscience/kmath/histogram/DoubleHistogramSpace.kt deleted file mode 100644 index 61d0b9f33..000000000 --- a/kmath-histograms/src/commonMain/kotlin/space/kscience/kmath/histogram/DoubleHistogramSpace.kt +++ /dev/null @@ -1,130 +0,0 @@ -/* - * Copyright 2018-2021 KMath contributors. - * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file. - */ - -package space.kscience.kmath.histogram - -import space.kscience.kmath.domains.Domain -import space.kscience.kmath.domains.HyperSquareDomain -import space.kscience.kmath.misc.UnstableKMathAPI -import space.kscience.kmath.nd.* -import space.kscience.kmath.operations.DoubleField -import space.kscience.kmath.structures.* -import kotlin.math.floor - -public class DoubleHistogramSpace( - private val lower: Buffer, - private val upper: Buffer, - private val binNums: IntArray = IntArray(lower.size) { 20 }, -) : IndexedHistogramSpace { - - 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 - - override val shape: IntArray = IntArray(binNums.size) { binNums[it] + 2 } - override val histogramValueSpace: DoubleFieldND = DoubleField.ndAlgebra(*shape) - - private val binSize = DoubleBuffer(dimension) { (upper[it] - lower[it]) / binNums[it] } - - /** - * Get internal [StructureND] 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): IntArray = IntArray(dimension) { - getIndex(it, point[it]) - } - - @OptIn(UnstableKMathAPI::class) - override fun getDomain(index: IntArray): Domain { - val lowerBoundary = index.mapIndexed { axis, i -> - when (i) { - 0 -> Double.NEGATIVE_INFINITY - shape[axis] - 1 -> upper[axis] - else -> lower[axis] + (i.toDouble()) * binSize[axis] - } - }.asBuffer() - - val upperBoundary = index.mapIndexed { axis, i -> - when (i) { - 0 -> lower[axis] - 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 { - val domain = getDomain(index) - return DomainBin(domain, value) - } - - override fun produce(builder: HistogramBuilder.() -> Unit): IndexedHistogram { - val ndCounter = StructureND.auto(shape) { Counter.double() } - val hBuilder = HistogramBuilder { point, value -> - val index = getIndex(point) - ndCounter[index].add(value.toDouble()) - } - hBuilder.apply(builder) - val values: BufferND = ndCounter.mapToBuffer { it.value } - return IndexedHistogram(this, values) - } - - override fun IndexedHistogram.unaryMinus(): IndexedHistogram = this * (-1) - - public companion object { - /** - * Use it like - * ``` - *FastHistogram.fromRanges( - * (-1.0..1.0), - * (-1.0..1.0) - *) - *``` - */ - public fun fromRanges(vararg ranges: ClosedFloatingPointRange): DoubleHistogramSpace = DoubleHistogramSpace( - ranges.map(ClosedFloatingPointRange::start).asBuffer(), - ranges.map(ClosedFloatingPointRange::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, Int>): DoubleHistogramSpace = - DoubleHistogramSpace( - ListBuffer( - ranges - .map(Pair, Int>::first) - .map(ClosedFloatingPointRange::start) - ), - - ListBuffer( - ranges - .map(Pair, Int>::first) - .map(ClosedFloatingPointRange::endInclusive) - ), - - ranges.map(Pair, Int>::second).toIntArray() - ) - } -} \ No newline at end of file diff --git a/kmath-histograms/src/commonMain/kotlin/space/kscience/kmath/histogram/Histogram.kt b/kmath-histograms/src/commonMain/kotlin/space/kscience/kmath/histogram/Histogram.kt index 4e803fc63..12edafd81 100644 --- a/kmath-histograms/src/commonMain/kotlin/space/kscience/kmath/histogram/Histogram.kt +++ b/kmath-histograms/src/commonMain/kotlin/space/kscience/kmath/histogram/Histogram.kt @@ -13,14 +13,23 @@ import space.kscience.kmath.structures.asBuffer /** * The binned data element. Could be a histogram bin with a number of counts or an artificial construct. */ -public interface Bin : Domain { +public interface Bin : Domain { /** * The value of this bin. */ - public val value: Number + public val binValue: V } -public interface Histogram> { +/** + * A simple histogram bin based on domain + */ +public data class DomainBin, D : Domain, out V>( + public val domain: D, + override val binValue: V, +) : Bin, Domain by domain + + +public interface Histogram> { /** * Find existing bin, corresponding to given coordinates */ @@ -32,29 +41,38 @@ public interface Histogram> { public val dimension: Int public val bins: Iterable + + public companion object { + //A discoverability root + } } -public fun interface HistogramBuilder { +public interface HistogramBuilder { /** - * Increment appropriate bin + * The default value increment for a bin */ - public fun putValue(point: Point, value: Number) + public val defaultValue: V + + /** + * Increment appropriate bin with given value + */ + public fun putValue(point: Point, value: V = defaultValue) } -public fun > HistogramBuilder.put(point: Point): Unit = putValue(point, 1.0) +public fun HistogramBuilder.put(point: Point): Unit = putValue(point) -public fun HistogramBuilder.put(vararg point: T): Unit = put(point.asBuffer()) +public fun HistogramBuilder.put(vararg point: T): Unit = put(point.asBuffer()) -public fun HistogramBuilder.put(vararg point: Number): Unit = +public fun HistogramBuilder.put(vararg point: Number): Unit = put(DoubleBuffer(point.map { it.toDouble() }.toDoubleArray())) -public fun HistogramBuilder.put(vararg point: Double): Unit = put(DoubleBuffer(point)) -public fun HistogramBuilder.fill(sequence: Iterable>): Unit = sequence.forEach { put(it) } +public fun HistogramBuilder.put(vararg point: Double): Unit = put(DoubleBuffer(point)) +public fun HistogramBuilder.fill(sequence: Iterable>): Unit = sequence.forEach { put(it) } /** * Pass a sequence builder into histogram */ -public fun HistogramBuilder.fill(block: suspend SequenceScope>.() -> Unit): Unit = +public fun HistogramBuilder.fill(block: suspend SequenceScope>.() -> Unit): Unit = fill(sequence(block).asIterable()) diff --git a/kmath-histograms/src/commonMain/kotlin/space/kscience/kmath/histogram/Histogram1D.kt b/kmath-histograms/src/commonMain/kotlin/space/kscience/kmath/histogram/Histogram1D.kt new file mode 100644 index 000000000..710e4478b --- /dev/null +++ b/kmath-histograms/src/commonMain/kotlin/space/kscience/kmath/histogram/Histogram1D.kt @@ -0,0 +1,64 @@ +/* + * Copyright 2018-2021 KMath contributors. + * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file. + */ + +package space.kscience.kmath.histogram + +import space.kscience.kmath.domains.Domain1D +import space.kscience.kmath.domains.center +import space.kscience.kmath.linear.Point +import space.kscience.kmath.misc.UnstableKMathAPI +import space.kscience.kmath.operations.asSequence +import space.kscience.kmath.structures.Buffer + + +/** + * A univariate bin based on a range + * + * @property binValue The value of histogram including weighting + */ +@UnstableKMathAPI +public data class Bin1D, out V>( + public val domain: Domain1D, + override val binValue: V, +) : Bin, ClosedRange by domain.range { + + override val dimension: Int get() = 1 + + override fun contains(point: Buffer): Boolean = point.size == 1 && contains(point[0]) +} + +@OptIn(UnstableKMathAPI::class) +public interface Histogram1D, V> : Histogram> { + override val dimension: Int get() = 1 + public operator fun get(value: T): Bin1D? + override operator fun get(point: Buffer): Bin1D? = get(point[0]) +} + +public interface Histogram1DBuilder : HistogramBuilder { + /** + * Thread safe put operation + */ + public fun putValue(at: T, value: V = defaultValue) + + override fun putValue(point: Point, value: V) { + require(point.size == 1) { "Only points with single value could be used in Histogram1D" } + putValue(point[0], value) + } +} + +@UnstableKMathAPI +public fun Histogram1DBuilder.fill(items: Iterable): Unit = + items.forEach(this::putValue) + +@UnstableKMathAPI +public fun Histogram1DBuilder.fill(array: DoubleArray): Unit = + array.forEach(this::putValue) + +@UnstableKMathAPI +public fun Histogram1DBuilder.fill(buffer: Buffer): Unit = + buffer.asSequence().forEach(this::putValue) + +@OptIn(UnstableKMathAPI::class) +public val Bin1D.center: Double get() = domain.center \ No newline at end of file diff --git a/kmath-histograms/src/commonMain/kotlin/space/kscience/kmath/histogram/HistogramND.kt b/kmath-histograms/src/commonMain/kotlin/space/kscience/kmath/histogram/HistogramND.kt new file mode 100644 index 000000000..68b24db5d --- /dev/null +++ b/kmath-histograms/src/commonMain/kotlin/space/kscience/kmath/histogram/HistogramND.kt @@ -0,0 +1,76 @@ +/* + * Copyright 2018-2021 KMath contributors. + * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file. + */ + +package space.kscience.kmath.histogram + +import space.kscience.kmath.domains.Domain +import space.kscience.kmath.linear.Point +import space.kscience.kmath.nd.DefaultStrides +import space.kscience.kmath.nd.FieldOpsND +import space.kscience.kmath.nd.Shape +import space.kscience.kmath.nd.StructureND +import space.kscience.kmath.operations.Group +import space.kscience.kmath.operations.ScaleOperations +import space.kscience.kmath.operations.invoke + +/** + * @param T the type of the argument space + * @param V the type of bin value + */ +public class HistogramND, D : Domain, V : Any>( + public val group: HistogramGroupND, + internal val values: StructureND, +) : Histogram> { + + override fun get(point: Point): DomainBin? { + val index = group.getIndexOrNull(point) ?: return null + return group.produceBin(index, values[index]) + } + + override val dimension: Int get() = group.shape.size + + override val bins: Iterable> + get() = DefaultStrides(group.shape).asSequence().map { + group.produceBin(it, values[it]) + }.asIterable() +} + +/** + * A space for producing histograms with values in a NDStructure + */ +public interface HistogramGroupND, D : Domain, V : Any> : + Group>, ScaleOperations> { + public val shape: Shape + public val valueAlgebraND: FieldOpsND //= NDAlgebra.space(valueSpace, Buffer.Companion::boxing, *shape), + + /** + * Resolve index of the bin including given [point]. Return null if point is outside histogram area + */ + public fun getIndexOrNull(point: Point): IntArray? + + /** + * Get a bin domain represented by given index + */ + public fun getDomain(index: IntArray): Domain? + + public fun produceBin(index: IntArray, value: V): DomainBin + + public fun produce(builder: HistogramBuilder.() -> Unit): HistogramND + + override fun add(left: HistogramND, right: HistogramND): HistogramND { + require(left.group == this && right.group == this) { + "A histogram belonging to a different group cannot be operated." + } + return HistogramND(this, valueAlgebraND { left.values + right.values }) + } + + override fun scale(a: HistogramND, value: Double): HistogramND { + require(a.group == this) { "A histogram belonging to a different group cannot be operated." } + return HistogramND(this, valueAlgebraND { a.values * value }) + } + + override val zero: HistogramND get() = produce { } +} + diff --git a/kmath-histograms/src/commonMain/kotlin/space/kscience/kmath/histogram/IndexedHistogramSpace.kt b/kmath-histograms/src/commonMain/kotlin/space/kscience/kmath/histogram/IndexedHistogramSpace.kt deleted file mode 100644 index 9275c1c5e..000000000 --- a/kmath-histograms/src/commonMain/kotlin/space/kscience/kmath/histogram/IndexedHistogramSpace.kt +++ /dev/null @@ -1,82 +0,0 @@ -/* - * Copyright 2018-2021 KMath contributors. - * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file. - */ - -package space.kscience.kmath.histogram - -import space.kscience.kmath.domains.Domain -import space.kscience.kmath.linear.Point -import space.kscience.kmath.misc.UnstableKMathAPI -import space.kscience.kmath.nd.DefaultStrides -import space.kscience.kmath.nd.FieldND -import space.kscience.kmath.nd.Shape -import space.kscience.kmath.nd.StructureND -import space.kscience.kmath.operations.Group -import space.kscience.kmath.operations.ScaleOperations -import space.kscience.kmath.operations.invoke - -/** - * A simple histogram bin based on domain - */ -public data class DomainBin>( - public val domain: Domain, - override val value: Number, -) : Bin, Domain by domain - -@OptIn(UnstableKMathAPI::class) -public class IndexedHistogram, V : Any>( - public val context: IndexedHistogramSpace, - public val values: StructureND, -) : Histogram> { - - override fun get(point: Point): Bin? { - val index = context.getIndex(point) ?: return null - return context.produceBin(index, values[index]) - } - - override val dimension: Int get() = context.shape.size - - override val bins: Iterable> - get() = DefaultStrides(context.shape).asSequence().map { - context.produceBin(it, values[it]) - }.asIterable() -} - -/** - * A space for producing histograms with values in a NDStructure - */ -public interface IndexedHistogramSpace, V : Any> - : Group>, ScaleOperations> { - //public val valueSpace: Space - public val shape: Shape - public val histogramValueSpace: FieldND //= NDAlgebra.space(valueSpace, Buffer.Companion::boxing, *shape), - - /** - * Resolve index of the bin including given [point] - */ - public fun getIndex(point: Point): IntArray? - - /** - * Get a bin domain represented by given index - */ - public fun getDomain(index: IntArray): Domain? - - public fun produceBin(index: IntArray, value: V): Bin - - public fun produce(builder: HistogramBuilder.() -> Unit): IndexedHistogram - - override fun add(left: IndexedHistogram, right: IndexedHistogram): IndexedHistogram { - require(left.context == this) { "Can't operate on a histogram produced by external space" } - require(right.context == this) { "Can't operate on a histogram produced by external space" } - return IndexedHistogram(this, histogramValueSpace { left.values + right.values }) - } - - override fun scale(a: IndexedHistogram, value: Double): IndexedHistogram { - require(a.context == this) { "Can't operate on a histogram produced by external space" } - return IndexedHistogram(this, histogramValueSpace { a.values * value }) - } - - override val zero: IndexedHistogram get() = produce { } -} - diff --git a/kmath-histograms/src/commonMain/kotlin/space/kscience/kmath/histogram/UniformHistogram1D.kt b/kmath-histograms/src/commonMain/kotlin/space/kscience/kmath/histogram/UniformHistogram1D.kt new file mode 100644 index 000000000..e13928394 --- /dev/null +++ b/kmath-histograms/src/commonMain/kotlin/space/kscience/kmath/histogram/UniformHistogram1D.kt @@ -0,0 +1,156 @@ +/* + * Copyright 2018-2021 KMath contributors. + * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file. + */ + +package space.kscience.kmath.histogram + +import space.kscience.kmath.domains.DoubleDomain1D +import space.kscience.kmath.misc.UnstableKMathAPI +import space.kscience.kmath.operations.Group +import space.kscience.kmath.operations.Ring +import space.kscience.kmath.operations.ScaleOperations +import space.kscience.kmath.operations.invoke +import space.kscience.kmath.structures.Buffer +import kotlin.math.floor + +@OptIn(UnstableKMathAPI::class) +public class UniformHistogram1D( + public val group: UniformHistogram1DGroup, + internal val values: Map, +) : Histogram1D { + + private val startPoint get() = group.startPoint + private val binSize get() = group.binSize + + private fun produceBin(index: Int, value: V): Bin1D { + val domain = DoubleDomain1D((startPoint + index * binSize)..(startPoint + (index + 1) * binSize)) + return Bin1D(domain, value) + } + + override val bins: Iterable> get() = values.map { produceBin(it.key, it.value) } + + override fun get(value: Double): Bin1D? { + val index: Int = group.getIndex(value) + val v = values[index] + return v?.let { produceBin(index, it) } + } +} + +/** + * An algebra for uniform histograms in 1D real space + */ +public class UniformHistogram1DGroup( + public val valueAlgebra: A, + public val binSize: Double, + public val startPoint: Double = 0.0, +) : Group>, ScaleOperations> where A : Ring, A : ScaleOperations { + + override val zero: UniformHistogram1D = UniformHistogram1D(this, emptyMap()) + + /** + * Get index of a bin + */ + @PublishedApi + internal fun getIndex(at: Double): Int = floor((at - startPoint) / binSize).toInt() + + override fun add( + left: Histogram1D, + right: Histogram1D, + ): UniformHistogram1D = valueAlgebra { + val leftUniform = produceFrom(left) + val rightUniform = produceFrom(right) + val keys = leftUniform.values.keys + rightUniform.values.keys + UniformHistogram1D( + this@UniformHistogram1DGroup, + keys.associateWith { + (leftUniform.values[it] ?: valueAlgebra.zero) + (rightUniform.values[it] ?: valueAlgebra.zero) + } + ) + } + + override fun Histogram1D.unaryMinus(): UniformHistogram1D = valueAlgebra { + UniformHistogram1D(this@UniformHistogram1DGroup, produceFrom(this@unaryMinus).values.mapValues { -it.value }) + } + + override fun scale( + a: Histogram1D, + value: Double, + ): UniformHistogram1D = UniformHistogram1D( + this@UniformHistogram1DGroup, + produceFrom(a).values.mapValues { valueAlgebra.scale(it.value, value) } + ) + + /** + * Fill histogram. + */ + public inline fun produce(block: Histogram1DBuilder.() -> Unit): UniformHistogram1D { + val map = HashMap() + val builder = object : Histogram1DBuilder { + override val defaultValue: V get() = valueAlgebra.zero + + override fun putValue(at: Double, value: V) { + val index = getIndex(at) + map[index] = with(valueAlgebra) { (map[index] ?: zero) + one } + } + } + builder.block() + return UniformHistogram1D(this, map) + } + + /** + * Re-bin given histogram to be compatible if exiting bin fully falls inside existing bin, this bin value + * is increased by one. If not, all bins including values from this bin are increased by fraction + * (conserving the norming). + */ + @OptIn(UnstableKMathAPI::class) + public fun produceFrom( + histogram: Histogram1D, + ): UniformHistogram1D = if ((histogram as? UniformHistogram1D)?.group == this) { + histogram + } else { + val map = HashMap() + histogram.bins.forEach { bin -> + val range = bin.domain.range + val indexOfLeft = getIndex(range.start) + val indexOfRight = getIndex(range.endInclusive) + val numBins = indexOfRight - indexOfLeft + 1 + for (i in indexOfLeft..indexOfRight) { + map[indexOfLeft] = with(valueAlgebra) { + (map[indexOfLeft] ?: zero) + bin.binValue / numBins + } + } + } + UniformHistogram1D(this, map) + } +} + +public fun Histogram.Companion.uniform1D( + valueAlgebra: A, + binSize: Double, + startPoint: Double = 0.0, +): UniformHistogram1DGroup where A : Ring, A : ScaleOperations = + UniformHistogram1DGroup(valueAlgebra, binSize, startPoint) + +@UnstableKMathAPI +public fun UniformHistogram1DGroup.produce( + buffer: Buffer, +): UniformHistogram1D = produce { fill(buffer) } + +/** + * Map of bin centers to bin values + */ +@OptIn(UnstableKMathAPI::class) +public val UniformHistogram1D.binValues: Map + get() = bins.associate { it.center to it.binValue } + + +//TODO add normalized values inside Field-based histogram spaces with context receivers +///** +// * Map of bin centers to normalized bin values (bin size as normalization) +// */ +//@OptIn(UnstableKMathAPI::class) +//public val UniformHistogram1D.binValuesNormalized: Map +// get() = group.valueAlgebra { +// bins.associate { it.center to it.binValue / group.binSize } +// } \ No newline at end of file diff --git a/kmath-histograms/src/commonMain/kotlin/space/kscience/kmath/histogram/UniformHistogramGroupND.kt b/kmath-histograms/src/commonMain/kotlin/space/kscience/kmath/histogram/UniformHistogramGroupND.kt new file mode 100644 index 000000000..90ec29ce3 --- /dev/null +++ b/kmath-histograms/src/commonMain/kotlin/space/kscience/kmath/histogram/UniformHistogramGroupND.kt @@ -0,0 +1,163 @@ +/* + * Copyright 2018-2021 KMath contributors. + * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file. + */ + +@file:OptIn(UnstableKMathAPI::class) + +package space.kscience.kmath.histogram + +import space.kscience.kmath.domains.HyperSquareDomain +import space.kscience.kmath.linear.Point +import space.kscience.kmath.misc.UnstableKMathAPI +import space.kscience.kmath.nd.* +import space.kscience.kmath.operations.DoubleField +import space.kscience.kmath.operations.Field +import space.kscience.kmath.operations.invoke +import space.kscience.kmath.structures.* +import kotlin.math.floor + +public typealias HyperSquareBin = DomainBin + +/** + * Multivariate histogram space for hyper-square real-field bins. + * @param bufferFactory is an optional parameter used to optimize buffer production. + */ +public class UniformHistogramGroupND>( + override val valueAlgebraND: FieldOpsND, + private val lower: Buffer, + private val upper: Buffer, + private val binNums: IntArray = IntArray(lower.size) { 20 }, + private val bufferFactory: BufferFactory = Buffer.Companion::boxing, +) : HistogramGroupND { + + 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 + + override val shape: IntArray = IntArray(binNums.size) { binNums[it] + 2 } + + private val binSize = DoubleBuffer(dimension) { (upper[it] - lower[it]) / binNums[it] } + + /** + * Get internal [StructureND] 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 getIndexOrNull(point: Buffer): IntArray = IntArray(dimension) { + getIndex(it, point[it]) + } + + override fun getDomain(index: IntArray): HyperSquareDomain { + val lowerBoundary = index.mapIndexed { axis, i -> + when (i) { + 0 -> Double.NEGATIVE_INFINITY + shape[axis] - 1 -> upper[axis] + else -> lower[axis] + (i.toDouble()) * binSize[axis] + } + }.asBuffer() + + val upperBoundary = index.mapIndexed { axis, i -> + when (i) { + 0 -> lower[axis] + 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: V): HyperSquareBin { + val domain = getDomain(index) + return DomainBin(domain, value) + } + + + override fun produce(builder: HistogramBuilder.() -> Unit): HistogramND { + val ndCounter = StructureND.buffered(shape) { Counter.of(valueAlgebraND.elementAlgebra) } + val hBuilder = object : HistogramBuilder { + override val defaultValue: V get() = valueAlgebraND.elementAlgebra.one + + override fun putValue(point: Point, value: V) = with(valueAlgebraND.elementAlgebra) { + val index = getIndexOrNull(point) + ndCounter[index].add(value) + } + } + hBuilder.apply(builder) + val values: BufferND = ndCounter.mapToBuffer(bufferFactory) { it.value } + return HistogramND(this, values) + } + + override fun HistogramND.unaryMinus(): HistogramND = + this * (-1) +} + +/** + * Use it like + * ``` + *FastHistogram.fromRanges( + * (-1.0..1.0), + * (-1.0..1.0) + *) + *``` + */ +public fun > Histogram.Companion.uniformNDFromRanges( + valueAlgebraND: FieldOpsND, + vararg ranges: ClosedFloatingPointRange, + bufferFactory: BufferFactory = Buffer.Companion::boxing, +): UniformHistogramGroupND = UniformHistogramGroupND( + valueAlgebraND, + ranges.map(ClosedFloatingPointRange::start).asBuffer(), + ranges.map(ClosedFloatingPointRange::endInclusive).asBuffer(), + bufferFactory = bufferFactory +) + +public fun Histogram.Companion.uniformDoubleNDFromRanges( + vararg ranges: ClosedFloatingPointRange, +): UniformHistogramGroupND = + uniformNDFromRanges(DoubleFieldOpsND, *ranges, bufferFactory = ::DoubleBuffer) + + +/** + * Use it like + * ``` + *FastHistogram.fromRanges( + * (-1.0..1.0) to 50, + * (-1.0..1.0) to 32 + *) + *``` + */ +public fun > Histogram.Companion.uniformNDFromRanges( + valueAlgebraND: FieldOpsND, + vararg ranges: Pair, Int>, + bufferFactory: BufferFactory = Buffer.Companion::boxing, +): UniformHistogramGroupND = UniformHistogramGroupND( + valueAlgebraND, + ListBuffer( + ranges + .map(Pair, Int>::first) + .map(ClosedFloatingPointRange::start) + ), + ListBuffer( + ranges + .map(Pair, Int>::first) + .map(ClosedFloatingPointRange::endInclusive) + ), + ranges.map(Pair, Int>::second).toIntArray(), + bufferFactory = bufferFactory +) + +public fun Histogram.Companion.uniformDoubleNDFromRanges( + vararg ranges: Pair, Int>, +): UniformHistogramGroupND = + uniformNDFromRanges(DoubleFieldOpsND, *ranges, bufferFactory = ::DoubleBuffer) \ No newline at end of file diff --git a/kmath-histograms/src/commonTest/kotlin/space/kscience/kmath/histogram/MultivariateHistogramTest.kt b/kmath-histograms/src/commonTest/kotlin/space/kscience/kmath/histogram/MultivariateHistogramTest.kt index 923cc98de..ca7c2f324 100644 --- a/kmath-histograms/src/commonTest/kotlin/space/kscience/kmath/histogram/MultivariateHistogramTest.kt +++ b/kmath-histograms/src/commonTest/kotlin/space/kscience/kmath/histogram/MultivariateHistogramTest.kt @@ -3,8 +3,11 @@ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file. */ +@file:OptIn(UnstableKMathAPI::class) + package space.kscience.kmath.histogram +import space.kscience.kmath.misc.UnstableKMathAPI import space.kscience.kmath.nd.DefaultStrides import space.kscience.kmath.operations.invoke import space.kscience.kmath.real.DoubleVector @@ -14,14 +17,14 @@ import kotlin.test.* internal class MultivariateHistogramTest { @Test fun testSinglePutHistogram() { - val hSpace = DoubleHistogramSpace.fromRanges( + val hSpace = Histogram.uniformDoubleNDFromRanges( (-1.0..1.0), (-1.0..1.0) ) val histogram = hSpace.produce { put(0.55, 0.55) } - val bin = histogram.bins.find { it.value.toInt() > 0 } ?: fail() + val bin = histogram.bins.find { it.binValue.toInt() > 0 } ?: fail() assertTrue { bin.contains(DoubleVector(0.55, 0.55)) } assertTrue { bin.contains(DoubleVector(0.6, 0.5)) } assertFalse { bin.contains(DoubleVector(-0.55, 0.55)) } @@ -29,7 +32,7 @@ internal class MultivariateHistogramTest { @Test fun testSequentialPut() { - val hSpace = DoubleHistogramSpace.fromRanges( + val hSpace = Histogram.uniformDoubleNDFromRanges( (-1.0..1.0), (-1.0..1.0), (-1.0..1.0) @@ -44,12 +47,12 @@ internal class MultivariateHistogramTest { put(nextDouble(), nextDouble(), nextDouble()) } } - assertEquals(n, histogram.bins.sumOf { it.value.toInt() }) + assertEquals(n, histogram.bins.sumOf { it.binValue.toInt() }) } @Test fun testHistogramAlgebra() { - DoubleHistogramSpace.fromRanges( + Histogram.uniformDoubleNDFromRanges( (-1.0..1.0), (-1.0..1.0), (-1.0..1.0) @@ -77,7 +80,7 @@ internal class MultivariateHistogramTest { assertTrue { res.bins.count() >= histogram1.bins.count() } - assertEquals(0.0, res.bins.sumOf { it.value.toDouble() }) + assertEquals(0.0, res.bins.sumOf { it.binValue.toDouble() }) } } } \ No newline at end of file diff --git a/kmath-histograms/src/commonTest/kotlin/space/kscience/kmath/histogram/UniformHistogram1DTest.kt b/kmath-histograms/src/commonTest/kotlin/space/kscience/kmath/histogram/UniformHistogram1DTest.kt new file mode 100644 index 000000000..09bf3939d --- /dev/null +++ b/kmath-histograms/src/commonTest/kotlin/space/kscience/kmath/histogram/UniformHistogram1DTest.kt @@ -0,0 +1,56 @@ +/* + * Copyright 2018-2021 KMath contributors. + * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file. + */ + +package space.kscience.kmath.histogram + +import kotlinx.coroutines.ExperimentalCoroutinesApi +import kotlinx.coroutines.test.runTest +import space.kscience.kmath.distributions.NormalDistribution +import space.kscience.kmath.misc.UnstableKMathAPI +import space.kscience.kmath.operations.DoubleField +import space.kscience.kmath.stat.RandomGenerator +import space.kscience.kmath.stat.nextBuffer +import kotlin.native.concurrent.ThreadLocal +import kotlin.test.Test +import kotlin.test.assertEquals + +@OptIn(ExperimentalCoroutinesApi::class, UnstableKMathAPI::class) +internal class UniformHistogram1DTest { + + @Test + fun normal() = runTest { + val distribution = NormalDistribution(0.0, 1.0) + with(Histogram.uniform1D(DoubleField, 0.1)) { + val h1 = produce(distribution.nextBuffer(generator, 10000)) + + val h2 = produce(distribution.nextBuffer(generator, 50000)) + + val h3 = h1 + h2 + + assertEquals(60000, h3.bins.sumOf { it.binValue }.toInt()) + } + } + + @Test + fun rebinDown() = runTest { + val h1 = Histogram.uniform1D(DoubleField, 0.01).produce(generator.nextDoubleBuffer(10000)) + val h2 = Histogram.uniform1D(DoubleField,0.03).produceFrom(h1) + + assertEquals(10000, h2.bins.sumOf { it.binValue }.toInt()) + } + + @Test + fun rebinUp() = runTest { + val h1 = Histogram.uniform1D(DoubleField, 0.03).produce(generator.nextDoubleBuffer(10000)) + val h2 = Histogram.uniform1D(DoubleField,0.01).produceFrom(h1) + + assertEquals(10000, h2.bins.sumOf { it.binValue }.toInt()) + } + + @ThreadLocal + companion object{ + private val generator = RandomGenerator.default(123) + } +} \ No newline at end of file diff --git a/kmath-histograms/src/jvmMain/kotlin/space/kscience/kmath/histogram/TreeHistogramGroup.kt b/kmath-histograms/src/jvmMain/kotlin/space/kscience/kmath/histogram/TreeHistogramGroup.kt new file mode 100644 index 000000000..6bec01f9b --- /dev/null +++ b/kmath-histograms/src/jvmMain/kotlin/space/kscience/kmath/histogram/TreeHistogramGroup.kt @@ -0,0 +1,179 @@ +/* + * Copyright 2018-2021 KMath contributors. + * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file. + */ + +@file:OptIn(UnstableKMathAPI::class) + +package space.kscience.kmath.histogram + +import space.kscience.kmath.domains.DoubleDomain1D +import space.kscience.kmath.domains.center +import space.kscience.kmath.misc.UnstableKMathAPI +import space.kscience.kmath.misc.sorted +import space.kscience.kmath.operations.Group +import space.kscience.kmath.operations.Ring +import space.kscience.kmath.operations.ScaleOperations +import space.kscience.kmath.structures.Buffer +import space.kscience.kmath.structures.first +import space.kscience.kmath.structures.indices +import space.kscience.kmath.structures.last +import java.util.* + +private fun > TreeMap.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 +} + +//public data class ValueAndError(val value: Double, val error: Double) +// +//public typealias WeightedBin1D = Bin1D + +/** + * A histogram based on a tree map of values + */ +public class TreeHistogram( + private val binMap: TreeMap>, +) : Histogram1D { + override fun get(value: Double): Bin1D? = binMap.getBin(value) + override val bins: Collection> get() = binMap.values +} + +/** + * A space for univariate histograms with variable bin borders based on a tree map + */ +public class TreeHistogramGroup( + public val valueAlgebra: A, + @PublishedApi internal val binFactory: (Double) -> DoubleDomain1D, +) : Group>, ScaleOperations> where A : Ring, A : ScaleOperations { + + internal inner class DomainCounter(val domain: DoubleDomain1D, val counter: Counter = Counter.of(valueAlgebra)) : + ClosedRange by domain.range + + @PublishedApi + internal inner class TreeHistogramBuilder : Histogram1DBuilder { + + override val defaultValue: V get() = valueAlgebra.one + + private val bins: TreeMap = TreeMap() + + private fun createBin(value: Double): DomainCounter { + val binDefinition: DoubleDomain1D = binFactory(value) + val newBin = DomainCounter(binDefinition) + synchronized(this) { + bins[binDefinition.center] = newBin + } + return newBin + } + + /** + * Thread safe put operation + */ + override fun putValue(at: Double, value: V) { + (bins.getBin(at) ?: createBin(at)).counter.add(value) + } + + fun build(): TreeHistogram { + val map = bins.mapValuesTo(TreeMap>()) { (_, binCounter) -> + Bin1D(binCounter.domain, binCounter.counter.value) + } + return TreeHistogram(map) + } + } + + public inline fun produce(block: Histogram1DBuilder.() -> Unit): TreeHistogram = + TreeHistogramBuilder().apply(block).build() + + override fun add( + left: TreeHistogram, + right: TreeHistogram, + ): TreeHistogram { + val bins = TreeMap>().apply { + (left.bins.map { it.domain } union right.bins.map { it.domain }).forEach { def -> + put( + def.center, + Bin1D( + def, + with(valueAlgebra) { + (left[def.center]?.binValue ?: zero) + (right[def.center]?.binValue ?: zero) + } + ) + ) + } + } + return TreeHistogram(bins) + } + + override fun scale(a: TreeHistogram, value: Double): TreeHistogram { + val bins = TreeMap>().apply { + a.bins.forEach { bin -> + put( + bin.domain.center, + Bin1D(bin.domain, valueAlgebra.scale(bin.binValue, value)) + ) + } + } + + return TreeHistogram(bins) + } + + override fun TreeHistogram.unaryMinus(): TreeHistogram = this * (-1) + + override val zero: TreeHistogram = produce { } +} + + +///** +// * Build and fill a histogram with custom borders. Returns a read-only histogram. +// */ +//public inline fun Histogram.custom( +// borders: DoubleArray, +// builder: Histogram1DBuilder.() -> Unit, +//): TreeHistogram = custom(borders).fill(builder) +// +// +///** +// * Build and fill a [DoubleHistogram1D]. Returns a read-only histogram. +// */ +//public fun uniform( +// binSize: Double, +// start: Double = 0.0, +//): TreeHistogramSpace = TreeHistogramSpace { value -> +// val center = start + binSize * floor((value - start) / binSize + 0.5) +// DoubleDomain1D((center - binSize / 2)..(center + binSize / 2)) +//} + +/** + * Create a histogram group with custom cell borders + */ +public fun Histogram.Companion.custom1D( + valueAlgebra: A, + borders: Buffer, +): TreeHistogramGroup where A : Ring, A : ScaleOperations { + val sorted = borders.sorted() + + return TreeHistogramGroup(valueAlgebra) { value -> + when { + value <= sorted.first() -> DoubleDomain1D( + Double.NEGATIVE_INFINITY..sorted.first() + ) + + value > sorted.last() -> DoubleDomain1D( + sorted.last()..Double.POSITIVE_INFINITY + ) + + else -> { + val index = sorted.indices.first { value <= sorted[it] } + val left = sorted[index - 1] + val right = sorted[index] + DoubleDomain1D(left..right) + } + } + } +} \ No newline at end of file diff --git a/kmath-histograms/src/jvmMain/kotlin/space/kscience/kmath/histogram/TreeHistogramSpace.kt b/kmath-histograms/src/jvmMain/kotlin/space/kscience/kmath/histogram/TreeHistogramSpace.kt deleted file mode 100644 index 0853615e6..000000000 --- a/kmath-histograms/src/jvmMain/kotlin/space/kscience/kmath/histogram/TreeHistogramSpace.kt +++ /dev/null @@ -1,171 +0,0 @@ -/* - * Copyright 2018-2021 KMath contributors. - * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file. - */ - -package space.kscience.kmath.histogram - -import space.kscience.kmath.domains.UnivariateDomain -import space.kscience.kmath.misc.UnstableKMathAPI -import space.kscience.kmath.operations.Group -import space.kscience.kmath.operations.ScaleOperations -import space.kscience.kmath.structures.Buffer -import java.util.* -import kotlin.math.abs -import kotlin.math.floor -import kotlin.math.sqrt - -private fun > TreeMap.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( - private val binMap: TreeMap, -) : UnivariateHistogram { - override fun get(value: Double): UnivariateBin? = binMap.getBin(value) - override val dimension: Int get() = 1 - override val bins: Collection get() = binMap.values -} - -@OptIn(UnstableKMathAPI::class) -@PublishedApi -internal class TreeHistogramBuilder(val binFactory: (Double) -> UnivariateDomain) : UnivariateHistogramBuilder { - - internal class BinCounter(val domain: UnivariateDomain, val counter: Counter = Counter.double()) : - ClosedFloatingPointRange by domain.range - - private val bins: TreeMap = TreeMap() - - 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, value: Number) { - require(point.size == 1) { "Only points with single value could be used in univariate histogram" } - putValue(point[0], value.toDouble()) - } - - fun build(): TreeHistogram { - val map = bins.mapValuesTo(TreeMap()) { (_, binCounter) -> - val count = binCounter.counter.value - UnivariateBin(binCounter.domain, count, sqrt(count)) - } - return TreeHistogram(map) - } -} - -/** - * A space for univariate histograms with variable bin borders based on a tree map - */ -@UnstableKMathAPI -public class TreeHistogramSpace( - @PublishedApi internal val binFactory: (Double) -> UnivariateDomain, -) : Group, ScaleOperations { - - public inline fun fill(block: UnivariateHistogramBuilder.() -> Unit): UnivariateHistogram = - TreeHistogramBuilder(binFactory).apply(block).build() - - override fun add( - left: UnivariateHistogram, - right: 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().apply { - (left.bins.map { it.domain } union right.bins.map { it.domain }).forEach { def -> - put( - def.center, - UnivariateBin( - def, - value = (left[def.center]?.value ?: 0.0) + (right[def.center]?.value ?: 0.0), - standardDeviation = (left[def.center]?.standardDeviation - ?: 0.0) + (right[def.center]?.standardDeviation ?: 0.0) - ) - ) - } - } - return TreeHistogram(bins) - } - - override fun scale(a: UnivariateHistogram, value: Double): UnivariateHistogram { - val bins = TreeMap().apply { - a.bins.forEach { bin -> - put( - bin.domain.center, - UnivariateBin( - bin.domain, - value = bin.value * value, - standardDeviation = abs(bin.standardDeviation * value) - ) - ) - } - } - - return TreeHistogram(bins) - } - - override fun UnivariateHistogram.unaryMinus(): UnivariateHistogram = this * (-1) - - override val zero: UnivariateHistogram by lazy { fill { } } - - 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 * 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) - } - } - } - } - } -} \ No newline at end of file diff --git a/kmath-histograms/src/jvmMain/kotlin/space/kscience/kmath/histogram/UnivariateHistogram.kt b/kmath-histograms/src/jvmMain/kotlin/space/kscience/kmath/histogram/UnivariateHistogram.kt deleted file mode 100644 index ac0576a8e..000000000 --- a/kmath-histograms/src/jvmMain/kotlin/space/kscience/kmath/histogram/UnivariateHistogram.kt +++ /dev/null @@ -1,79 +0,0 @@ -/* - * Copyright 2018-2021 KMath contributors. - * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file. - */ - -package space.kscience.kmath.histogram - -import space.kscience.kmath.domains.UnivariateDomain -import space.kscience.kmath.misc.UnstableKMathAPI -import space.kscience.kmath.operations.asSequence -import space.kscience.kmath.structures.Buffer - - -@UnstableKMathAPI -public val UnivariateDomain.center: Double - get() = (range.endInclusive + range.start) / 2 - -/** - * A univariate bin based on a range - * - * @property value The value of histogram including weighting - * @property standardDeviation Standard deviation of the bin value. Zero or negative if not applicable - */ -@UnstableKMathAPI -public class UnivariateBin( - public val domain: UnivariateDomain, - override val value: Double, - public val standardDeviation: Double, -) : Bin, ClosedFloatingPointRange by domain.range { - - override val dimension: Int get() = 1 - - override fun contains(point: Buffer): Boolean = point.size == 1 && contains(point[0]) -} - -@OptIn(UnstableKMathAPI::class) -public interface UnivariateHistogram : Histogram { - public operator fun get(value: Double): UnivariateBin? - override operator fun get(point: Buffer): UnivariateBin? = get(point[0]) - - public companion object { - /** - * Build and fill a [UnivariateHistogram]. Returns a read-only histogram. - */ - public inline fun uniform( - binSize: Double, - start: Double = 0.0, - builder: UnivariateHistogramBuilder.() -> Unit, - ): UnivariateHistogram = TreeHistogramSpace.uniform(binSize, start).fill(builder) - - /** - * Build and fill a histogram with custom borders. Returns a read-only histogram. - */ - public inline fun custom( - borders: DoubleArray, - builder: UnivariateHistogramBuilder.() -> Unit, - ): UnivariateHistogram = TreeHistogramSpace.custom(borders).fill(builder) - - } -} - -@UnstableKMathAPI -public interface UnivariateHistogramBuilder : HistogramBuilder { - /** - * Thread safe put operation - */ - public fun putValue(at: Double, value: Double = 1.0) - - override fun putValue(point: Buffer, value: Number) -} - -@UnstableKMathAPI -public fun UnivariateHistogramBuilder.fill(items: Iterable): Unit = items.forEach(this::putValue) - -@UnstableKMathAPI -public fun UnivariateHistogramBuilder.fill(array: DoubleArray): Unit = array.forEach(this::putValue) - -@UnstableKMathAPI -public fun UnivariateHistogramBuilder.fill(buffer: Buffer): Unit = buffer.asSequence().forEach(this::putValue) \ No newline at end of file diff --git a/kmath-histograms/src/jvmTest/kotlin/space/kscience/kmath/histogram/TreeHistogramTest.kt b/kmath-histograms/src/jvmTest/kotlin/space/kscience/kmath/histogram/TreeHistogramTest.kt index 28a1b03cb..c4eeb53cc 100644 --- a/kmath-histograms/src/jvmTest/kotlin/space/kscience/kmath/histogram/TreeHistogramTest.kt +++ b/kmath-histograms/src/jvmTest/kotlin/space/kscience/kmath/histogram/TreeHistogramTest.kt @@ -6,19 +6,24 @@ package space.kscience.kmath.histogram import org.junit.jupiter.api.Test +import space.kscience.kmath.operations.DoubleField +import space.kscience.kmath.real.step import kotlin.random.Random +import kotlin.test.assertEquals import kotlin.test.assertTrue class TreeHistogramTest { @Test fun normalFill() { - val histogram = UnivariateHistogram.uniform(0.1) { + val random = Random(123) + val histogram = Histogram.custom1D(DoubleField, 0.0..1.0 step 0.1).produce { repeat(100_000) { - putValue(Random.nextDouble()) + putValue(random.nextDouble()) } } - assertTrue { histogram.bins.count() > 10 } + assertTrue { histogram.bins.count() > 8} + assertEquals(100_000, histogram.bins.sumOf { it.binValue }.toInt()) } } \ No newline at end of file diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/distributions/Distribution.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/distributions/Distribution.kt index 3d3f95f8f..8dbcb3367 100644 --- a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/distributions/Distribution.kt +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/distributions/Distribution.kt @@ -27,7 +27,7 @@ public interface Distribution : Sampler { public companion object } -public interface UnivariateDistribution> : Distribution { +public interface Distribution1D> : Distribution { /** * Cumulative distribution for ordered parameter (CDF) */ @@ -37,7 +37,7 @@ public interface UnivariateDistribution> : Distribution { /** * Compute probability integral in an interval */ -public fun > UnivariateDistribution.integral(from: T, to: T): Double { +public fun > Distribution1D.integral(from: T, to: T): Double { require(to > from) return cumulative(to) - cumulative(from) } diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/distributions/NormalDistribution.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/distributions/NormalDistribution.kt index 66e041f05..cfc6eba04 100644 --- a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/distributions/NormalDistribution.kt +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/distributions/NormalDistribution.kt @@ -14,14 +14,9 @@ import space.kscience.kmath.stat.RandomGenerator import kotlin.math.* /** - * Implements [UnivariateDistribution] for the normal (gaussian) distribution. + * Implements [Distribution1D] for the normal (gaussian) distribution. */ -public class NormalDistribution(public val sampler: GaussianSampler) : UnivariateDistribution { - public constructor( - mean: Double, - standardDeviation: Double, - normalized: NormalizedGaussianSampler = ZigguratNormalizedGaussianSampler, - ) : this(GaussianSampler(mean, standardDeviation, normalized)) +public class NormalDistribution(public val sampler: GaussianSampler) : Distribution1D { override fun probability(arg: Double): Double { val x1 = (arg - sampler.mean) / sampler.standardDeviation @@ -43,3 +38,9 @@ public class NormalDistribution(public val sampler: GaussianSampler) : Univariat private val SQRT2 = sqrt(2.0) } } + +public fun NormalDistribution( + mean: Double, + standardDeviation: Double, + normalized: NormalizedGaussianSampler = ZigguratNormalizedGaussianSampler, +): NormalDistribution = NormalDistribution(GaussianSampler(mean, standardDeviation, normalized)) diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/Sampler.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/Sampler.kt index 0dd121f3b..c96ecdc8c 100644 --- a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/Sampler.kt +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/Sampler.kt @@ -67,3 +67,12 @@ public fun Sampler.sampleBuffer(generator: RandomGenerator, size: Int): @JvmName("sampleIntBuffer") public fun Sampler.sampleBuffer(generator: RandomGenerator, size: Int): Chain> = sampleBuffer(generator, size, ::IntBuffer) + + +/** + * Samples a [Buffer] of values from this [Sampler]. + */ +public suspend fun Sampler.nextBuffer(generator: RandomGenerator, size: Int): Buffer = + sampleBuffer(generator, size).first() + +//TODO add `context(RandomGenerator) Sampler.nextBuffer \ No newline at end of file diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/UniformDistribution.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/UniformDistribution.kt index 970a3aab5..20cc0e802 100644 --- a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/UniformDistribution.kt +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/UniformDistribution.kt @@ -8,9 +8,9 @@ package space.kscience.kmath.stat import space.kscience.kmath.chains.Chain import space.kscience.kmath.chains.SimpleChain import space.kscience.kmath.distributions.Distribution -import space.kscience.kmath.distributions.UnivariateDistribution +import space.kscience.kmath.distributions.Distribution1D -public class UniformDistribution(public val range: ClosedFloatingPointRange) : UnivariateDistribution { +public class UniformDistribution(public val range: ClosedFloatingPointRange) : Distribution1D { private val length: Double = range.endInclusive - range.start override fun probability(arg: Double): Double = if (arg in range) 1.0 / length else 0.0