Merge pull request #476 from mipt-npm/refactor/histogram

Complete refactor of histograms API
This commit is contained in:
Iaroslav Postovalov 2022-04-11 00:59:04 +07:00 committed by GitHub
commit ff58985d78
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
31 changed files with 947 additions and 582 deletions

View File

@ -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`

View File

@ -11,7 +11,7 @@ allprojects {
}
group = "space.kscience"
version = "0.3.0-dev-20"
version = "0.3.0-dev-21"
}
subprojects {

View File

@ -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

View File

@ -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<T : Comparable<T>>(public val range: ClosedRange<T>) : Domain<T> {
override val dimension: Int get() = 1
public operator fun contains(value: T): Boolean = range.contains(value)
override operator fun contains(point: Point<T>): Boolean {
require(point.size == 0)
return contains(point[0])
}
}
@UnstableKMathAPI
public class DoubleDomain1D(
@Suppress("CanBeParameter") public val doubleRange: ClosedFloatingPointRange<Double>,
) : Domain1D<Double>(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<Double>.center: Double
get() = (range.endInclusive + range.start) / 2

View File

@ -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<Double>, private val upper: Buffer<Double>) : DoubleDomain {
public class HyperSquareDomain(public val lower: Buffer<Double>, public val upper: Buffer<Double>) : 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<Double>): Boolean = point.indices.all { i ->
point[i] in lower[i]..upper[i]
}

View File

@ -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<Double>) : DoubleDomain {
override val dimension: Int get() = 1
public operator fun contains(d: Double): Boolean = range.contains(d)
override operator fun contains(point: Point<Double>): 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
}

View File

@ -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 <V: Comparable<V>> Buffer<V>.permSort() : IntArray = _permSortWith(compareBy<Int> { get(it) })
public fun <V : Comparable<V>> Buffer<V>.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 <V : Comparable<V>> Buffer<V>.sorted(): Buffer<V> {
val permutations = indicesSorted()
return VirtualBuffer(size) { this[permutations[it]] }
}
@PerformancePitfall
@UnstableKMathAPI
public fun <V: Comparable<V>> Buffer<V>.permSortDescending() : IntArray = _permSortWith(compareByDescending<Int> { get(it) })
public fun <V : Comparable<V>> Buffer<V>.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 <V : Comparable<V>> Buffer<V>.sortedDescending(): Buffer<V> {
val permutations = indicesSortedDescending()
return VirtualBuffer(size) { this[permutations[it]] }
}
@PerformancePitfall
@UnstableKMathAPI
public fun <V, C: Comparable<C>> Buffer<V>.permSortBy(selector: (V) -> C) : IntArray = _permSortWith(compareBy<Int> { selector(get(it)) })
public fun <V, C : Comparable<C>> Buffer<V>.indicesSortedBy(selector: (V) -> C): IntArray =
permSortIndicesWith(compareBy { selector(get(it)) })
@OptIn(UnstableKMathAPI::class)
public fun <V, C : Comparable<C>> Buffer<V>.sortedBy(selector: (V) -> C): Buffer<V> {
val permutations = indicesSortedBy(selector)
return VirtualBuffer(size) { this[permutations[it]] }
}
@PerformancePitfall
@UnstableKMathAPI
public fun <V, C: Comparable<C>> Buffer<V>.permSortByDescending(selector: (V) -> C) : IntArray = _permSortWith(compareByDescending<Int> { selector(get(it)) })
public fun <V, C : Comparable<C>> Buffer<V>.indicesSortedByDescending(selector: (V) -> C): IntArray =
permSortIndicesWith(compareByDescending { selector(get(it)) })
@OptIn(UnstableKMathAPI::class)
public fun <V, C : Comparable<C>> Buffer<V>.sortedByDescending(selector: (V) -> C): Buffer<V> {
val permutations = indicesSortedByDescending(selector)
return VirtualBuffer(size) { this[permutations[it]] }
}
@PerformancePitfall
@UnstableKMathAPI
public fun <V> Buffer<V>.permSortWith(comparator : Comparator<V>) : IntArray = _permSortWith { i1, i2 -> comparator.compare(get(i1), get(i2)) }
public fun <V> Buffer<V>.indicesSortedWith(comparator: Comparator<V>): IntArray =
permSortIndicesWith { i1, i2 -> comparator.compare(get(i1), get(i2)) }
@PerformancePitfall
@UnstableKMathAPI
private fun <V> Buffer<V>._permSortWith(comparator : Comparator<Int>) : IntArray {
if (size < 2) return IntArray(size)
private fun <V> Buffer<V>.permSortIndicesWith(comparator: Comparator<Int>): 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 <V> Buffer<V>._permSortWith(comparator : Comparator<Int>) : IntArray
*/
return packedIndices.sortedWith(comparator).toIntArray()
}
/**
* Checks that the [Buffer] is sorted (ascending) and throws [IllegalArgumentException] if it is not.
*/
public fun <T : Comparable<T>> Buffer<T>.requireSorted() {
for (i in 0..(size - 2)) {
require(get(i + 1) >= get(i)) { "The buffer is not sorted at index $i" }
}
}

View File

@ -194,7 +194,7 @@ public interface RingOpsND<T, out A : RingOps<T>> : RingOps<StructureND<T>>, Gro
override fun multiply(left: StructureND<T>, right: StructureND<T>): StructureND<T> =
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.

View File

@ -32,18 +32,23 @@ public open class BufferND<out T>(
/**
* Transform structure to a new structure using provided [BufferFactory] and optimizing if argument is [BufferND]
*/
public inline fun <T, reified R : Any> StructureND<T>.mapToBuffer(
factory: BufferFactory<R> = Buffer.Companion::auto,
public inline fun <T, R : Any> StructureND<T>.mapToBuffer(
factory: BufferFactory<R>,
crossinline transform: (T) -> R,
): BufferND<R> {
return if (this is BufferND<T>)
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<R> = if (this is BufferND<T>)
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 <T, reified R : Any> StructureND<T>.mapToBuffer(
crossinline transform: (T) -> R,
): BufferND<R> = mapToBuffer(Buffer.Companion::auto, transform)
/**
* Represents [MutableStructureND] over [MutableBuffer].
*

View File

@ -105,6 +105,16 @@ public interface Buffer<out T> {
*/
public val Buffer<*>.indices: IntRange get() = 0 until size
public fun <T> Buffer<T>.first(): T {
require(size > 0) { "Can't get the first element of empty buffer" }
return get(0)
}
public fun <T> Buffer<T>.last(): T {
require(size > 0) { "Can't get the last element of empty buffer" }
return get(size - 1)
}
/**
* Immutable wrapper for [MutableBuffer].
*

View File

@ -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)

View File

@ -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"
}
}

View File

@ -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")
}
}
}

View File

@ -18,7 +18,8 @@ public interface Counter<T : Any> {
public val value: T
public companion object {
public fun double(): ObjectCounter<Double> = ObjectCounter(DoubleField)
public fun ofDouble(): ObjectCounter<Double> = ObjectCounter(DoubleField)
public fun <T: Any> of(group: Group<T>): ObjectCounter<T> = ObjectCounter(group)
}
}

View File

@ -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<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
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<Double>): IntArray = IntArray(dimension) {
getIndex(it, point[it])
}
@OptIn(UnstableKMathAPI::class)
override fun getDomain(index: IntArray): Domain<Double> {
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<Double> {
val domain = getDomain(index)
return DomainBin(domain, value)
}
override fun produce(builder: HistogramBuilder<Double>.() -> Unit): IndexedHistogram<Double, Double> {
val ndCounter = StructureND.auto(shape) { Counter.double() }
val hBuilder = HistogramBuilder<Double> { point, value ->
val index = getIndex(point)
ndCounter[index].add(value.toDouble())
}
hBuilder.apply(builder)
val values: BufferND<Double> = ndCounter.mapToBuffer { it.value }
return IndexedHistogram(this, values)
}
override fun IndexedHistogram<Double, Double>.unaryMinus(): IndexedHistogram<Double, Double> = this * (-1)
public companion object {
/**
* Use it like
* ```
*FastHistogram.fromRanges(
* (-1.0..1.0),
* (-1.0..1.0)
*)
*```
*/
public fun fromRanges(vararg ranges: ClosedFloatingPointRange<Double>): DoubleHistogramSpace = DoubleHistogramSpace(
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>): DoubleHistogramSpace =
DoubleHistogramSpace(
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

@ -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<in T : Any> : Domain<T> {
public interface Bin<in T : Any, out V> : Domain<T> {
/**
* The value of this bin.
*/
public val value: Number
public val binValue: V
}
public interface Histogram<in T : Any, out B : Bin<T>> {
/**
* A simple histogram bin based on domain
*/
public data class DomainBin<in T : Comparable<T>, D : Domain<T>, out V>(
public val domain: D,
override val binValue: V,
) : Bin<T, V>, Domain<T> by domain
public interface Histogram<in T : Any, out V, out B : Bin<T, V>> {
/**
* Find existing bin, corresponding to given coordinates
*/
@ -32,29 +41,38 @@ public interface Histogram<in T : Any, out B : Bin<T>> {
public val dimension: Int
public val bins: Iterable<B>
public companion object {
//A discoverability root
}
}
public fun interface HistogramBuilder<in T : Any> {
public interface HistogramBuilder<in T : Any, V : Any> {
/**
* Increment appropriate bin
* The default value increment for a bin
*/
public fun putValue(point: Point<out T>, value: Number)
public val defaultValue: V
/**
* Increment appropriate bin with given value
*/
public fun putValue(point: Point<out T>, value: V = defaultValue)
}
public fun <T : Any, B : Bin<T>> HistogramBuilder<T>.put(point: Point<out T>): Unit = putValue(point, 1.0)
public fun <T : Any> HistogramBuilder<T, *>.put(point: Point<out T>): Unit = putValue(point)
public fun <T : Any> HistogramBuilder<T>.put(vararg point: T): Unit = put(point.asBuffer())
public fun <T : Any> HistogramBuilder<T, *>.put(vararg point: T): Unit = put(point.asBuffer())
public fun HistogramBuilder<Double>.put(vararg point: Number): Unit =
public fun HistogramBuilder<Double, *>.put(vararg point: Number): Unit =
put(DoubleBuffer(point.map { it.toDouble() }.toDoubleArray()))
public fun HistogramBuilder<Double>.put(vararg point: Double): Unit = put(DoubleBuffer(point))
public fun <T : Any> HistogramBuilder<T>.fill(sequence: Iterable<Point<T>>): Unit = sequence.forEach { put(it) }
public fun HistogramBuilder<Double, *>.put(vararg point: Double): Unit = put(DoubleBuffer(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> HistogramBuilder<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,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<T : Comparable<T>, out V>(
public val domain: Domain1D<T>,
override val binValue: V,
) : Bin<T, V>, ClosedRange<T> by domain.range {
override val dimension: Int get() = 1
override fun contains(point: Buffer<T>): Boolean = point.size == 1 && contains(point[0])
}
@OptIn(UnstableKMathAPI::class)
public interface Histogram1D<T : Comparable<T>, V> : Histogram<T, V, Bin1D<T, V>> {
override val dimension: Int get() = 1
public operator fun get(value: T): Bin1D<T, V>?
override operator fun get(point: Buffer<T>): Bin1D<T, V>? = get(point[0])
}
public interface Histogram1DBuilder<in T : Any, V : Any> : HistogramBuilder<T, V> {
/**
* Thread safe put operation
*/
public fun putValue(at: T, value: V = defaultValue)
override fun putValue(point: Point<out T>, value: V) {
require(point.size == 1) { "Only points with single value could be used in Histogram1D" }
putValue(point[0], value)
}
}
@UnstableKMathAPI
public fun Histogram1DBuilder<Double, *>.fill(items: Iterable<Double>): Unit =
items.forEach(this::putValue)
@UnstableKMathAPI
public fun Histogram1DBuilder<Double, *>.fill(array: DoubleArray): Unit =
array.forEach(this::putValue)
@UnstableKMathAPI
public fun <T : Any> Histogram1DBuilder<T, *>.fill(buffer: Buffer<T>): Unit =
buffer.asSequence().forEach(this::putValue)
@OptIn(UnstableKMathAPI::class)
public val Bin1D<Double, *>.center: Double get() = domain.center

View File

@ -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<T : Comparable<T>, D : Domain<T>, V : Any>(
public val group: HistogramGroupND<T, D, V>,
internal val values: StructureND<V>,
) : Histogram<T, V, DomainBin<T, D, V>> {
override fun get(point: Point<T>): DomainBin<T, D, V>? {
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<DomainBin<T, D, V>>
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<T : Comparable<T>, D : Domain<T>, V : Any> :
Group<HistogramND<T, D, V>>, ScaleOperations<HistogramND<T, D, V>> {
public val shape: Shape
public val valueAlgebraND: FieldOpsND<V, *> //= 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<T>): IntArray?
/**
* Get a bin domain represented by given index
*/
public fun getDomain(index: IntArray): Domain<T>?
public fun produceBin(index: IntArray, value: V): DomainBin<T, D, V>
public fun produce(builder: HistogramBuilder<T, V>.() -> Unit): HistogramND<T, D, V>
override fun add(left: HistogramND<T, D, V>, right: HistogramND<T, D, V>): HistogramND<T, D, V> {
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<T, D, V>, value: Double): HistogramND<T, D, V> {
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<T, D, V> get() = produce { }
}

View File

@ -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<in T : Comparable<T>>(
public val domain: Domain<T>,
override val value: Number,
) : Bin<T>, Domain<T> by domain
@OptIn(UnstableKMathAPI::class)
public class IndexedHistogram<T : Comparable<T>, V : Any>(
public val context: IndexedHistogramSpace<T, V>,
public val values: StructureND<V>,
) : Histogram<T, Bin<T>> {
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.shape.size
override val bins: Iterable<Bin<T>>
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<T : Comparable<T>, V : Any>
: Group<IndexedHistogram<T, V>>, ScaleOperations<IndexedHistogram<T, V>> {
//public val valueSpace: Space<V>
public val shape: Shape
public val histogramValueSpace: FieldND<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(left: IndexedHistogram<T, V>, right: IndexedHistogram<T, V>): IndexedHistogram<T, V> {
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<T, V>, value: Double): IndexedHistogram<T, V> {
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<T, V> get() = produce { }
}

View File

@ -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<V : Any>(
public val group: UniformHistogram1DGroup<V, *>,
internal val values: Map<Int, V>,
) : Histogram1D<Double, V> {
private val startPoint get() = group.startPoint
private val binSize get() = group.binSize
private fun produceBin(index: Int, value: V): Bin1D<Double, V> {
val domain = DoubleDomain1D((startPoint + index * binSize)..(startPoint + (index + 1) * binSize))
return Bin1D(domain, value)
}
override val bins: Iterable<Bin1D<Double, V>> get() = values.map { produceBin(it.key, it.value) }
override fun get(value: Double): Bin1D<Double, V>? {
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<V : Any, A>(
public val valueAlgebra: A,
public val binSize: Double,
public val startPoint: Double = 0.0,
) : Group<Histogram1D<Double, V>>, ScaleOperations<Histogram1D<Double, V>> where A : Ring<V>, A : ScaleOperations<V> {
override val zero: UniformHistogram1D<V> = 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<Double, V>,
right: Histogram1D<Double, V>,
): UniformHistogram1D<V> = 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<Double, V>.unaryMinus(): UniformHistogram1D<V> = valueAlgebra {
UniformHistogram1D(this@UniformHistogram1DGroup, produceFrom(this@unaryMinus).values.mapValues { -it.value })
}
override fun scale(
a: Histogram1D<Double, V>,
value: Double,
): UniformHistogram1D<V> = UniformHistogram1D(
this@UniformHistogram1DGroup,
produceFrom(a).values.mapValues { valueAlgebra.scale(it.value, value) }
)
/**
* Fill histogram.
*/
public inline fun produce(block: Histogram1DBuilder<Double, V>.() -> Unit): UniformHistogram1D<V> {
val map = HashMap<Int, V>()
val builder = object : Histogram1DBuilder<Double, V> {
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<Double, V>,
): UniformHistogram1D<V> = if ((histogram as? UniformHistogram1D)?.group == this) {
histogram
} else {
val map = HashMap<Int, V>()
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 <V : Any, A> Histogram.Companion.uniform1D(
valueAlgebra: A,
binSize: Double,
startPoint: Double = 0.0,
): UniformHistogram1DGroup<V, A> where A : Ring<V>, A : ScaleOperations<V> =
UniformHistogram1DGroup(valueAlgebra, binSize, startPoint)
@UnstableKMathAPI
public fun <V : Any> UniformHistogram1DGroup<V, *>.produce(
buffer: Buffer<Double>,
): UniformHistogram1D<V> = produce { fill(buffer) }
/**
* Map of bin centers to bin values
*/
@OptIn(UnstableKMathAPI::class)
public val <V : Any> UniformHistogram1D<V>.binValues: Map<Double, V>
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 <V : Any> UniformHistogram1D<V>.binValuesNormalized: Map<Double, V>
// get() = group.valueAlgebra {
// bins.associate { it.center to it.binValue / group.binSize }
// }

View File

@ -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<V> = DomainBin<Double, HyperSquareDomain, V>
/**
* Multivariate histogram space for hyper-square real-field bins.
* @param bufferFactory is an optional parameter used to optimize buffer production.
*/
public class UniformHistogramGroupND<V : Any, A : Field<V>>(
override val valueAlgebraND: FieldOpsND<V, A>,
private val lower: Buffer<Double>,
private val upper: Buffer<Double>,
private val binNums: IntArray = IntArray(lower.size) { 20 },
private val bufferFactory: BufferFactory<V> = Buffer.Companion::boxing,
) : HistogramGroupND<Double, HyperSquareDomain, V> {
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<Double>): 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<V> {
val domain = getDomain(index)
return DomainBin(domain, value)
}
override fun produce(builder: HistogramBuilder<Double, V>.() -> Unit): HistogramND<Double, HyperSquareDomain, V> {
val ndCounter = StructureND.buffered(shape) { Counter.of(valueAlgebraND.elementAlgebra) }
val hBuilder = object : HistogramBuilder<Double, V> {
override val defaultValue: V get() = valueAlgebraND.elementAlgebra.one
override fun putValue(point: Point<out Double>, value: V) = with(valueAlgebraND.elementAlgebra) {
val index = getIndexOrNull(point)
ndCounter[index].add(value)
}
}
hBuilder.apply(builder)
val values: BufferND<V> = ndCounter.mapToBuffer(bufferFactory) { it.value }
return HistogramND(this, values)
}
override fun HistogramND<Double, HyperSquareDomain, V>.unaryMinus(): HistogramND<Double, HyperSquareDomain, V> =
this * (-1)
}
/**
* Use it like
* ```
*FastHistogram.fromRanges(
* (-1.0..1.0),
* (-1.0..1.0)
*)
*```
*/
public fun <V : Any, A : Field<V>> Histogram.Companion.uniformNDFromRanges(
valueAlgebraND: FieldOpsND<V, A>,
vararg ranges: ClosedFloatingPointRange<Double>,
bufferFactory: BufferFactory<V> = Buffer.Companion::boxing,
): UniformHistogramGroupND<V, A> = UniformHistogramGroupND(
valueAlgebraND,
ranges.map(ClosedFloatingPointRange<Double>::start).asBuffer(),
ranges.map(ClosedFloatingPointRange<Double>::endInclusive).asBuffer(),
bufferFactory = bufferFactory
)
public fun Histogram.Companion.uniformDoubleNDFromRanges(
vararg ranges: ClosedFloatingPointRange<Double>,
): UniformHistogramGroupND<Double, DoubleField> =
uniformNDFromRanges(DoubleFieldOpsND, *ranges, bufferFactory = ::DoubleBuffer)
/**
* Use it like
* ```
*FastHistogram.fromRanges(
* (-1.0..1.0) to 50,
* (-1.0..1.0) to 32
*)
*```
*/
public fun <V : Any, A : Field<V>> Histogram.Companion.uniformNDFromRanges(
valueAlgebraND: FieldOpsND<V, A>,
vararg ranges: Pair<ClosedFloatingPointRange<Double>, Int>,
bufferFactory: BufferFactory<V> = Buffer.Companion::boxing,
): UniformHistogramGroupND<V, A> = UniformHistogramGroupND(
valueAlgebraND,
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(),
bufferFactory = bufferFactory
)
public fun Histogram.Companion.uniformDoubleNDFromRanges(
vararg ranges: Pair<ClosedFloatingPointRange<Double>, Int>,
): UniformHistogramGroupND<Double, DoubleField> =
uniformNDFromRanges(DoubleFieldOpsND, *ranges, bufferFactory = ::DoubleBuffer)

View File

@ -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() })
}
}
}

View File

@ -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)
}
}

View File

@ -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 <B : ClosedRange<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
}
//public data class ValueAndError(val value: Double, val error: Double)
//
//public typealias WeightedBin1D = Bin1D<Double, ValueAndError>
/**
* A histogram based on a tree map of values
*/
public class TreeHistogram<V : Any>(
private val binMap: TreeMap<Double, Bin1D<Double, V>>,
) : Histogram1D<Double, V> {
override fun get(value: Double): Bin1D<Double, V>? = binMap.getBin(value)
override val bins: Collection<Bin1D<Double, V>> get() = binMap.values
}
/**
* A space for univariate histograms with variable bin borders based on a tree map
*/
public class TreeHistogramGroup<V : Any, A>(
public val valueAlgebra: A,
@PublishedApi internal val binFactory: (Double) -> DoubleDomain1D,
) : Group<TreeHistogram<V>>, ScaleOperations<TreeHistogram<V>> where A : Ring<V>, A : ScaleOperations<V> {
internal inner class DomainCounter(val domain: DoubleDomain1D, val counter: Counter<V> = Counter.of(valueAlgebra)) :
ClosedRange<Double> by domain.range
@PublishedApi
internal inner class TreeHistogramBuilder : Histogram1DBuilder<Double, V> {
override val defaultValue: V get() = valueAlgebra.one
private val bins: TreeMap<Double, DomainCounter> = 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<V> {
val map = bins.mapValuesTo(TreeMap<Double, Bin1D<Double, V>>()) { (_, binCounter) ->
Bin1D(binCounter.domain, binCounter.counter.value)
}
return TreeHistogram(map)
}
}
public inline fun produce(block: Histogram1DBuilder<Double, V>.() -> Unit): TreeHistogram<V> =
TreeHistogramBuilder().apply(block).build()
override fun add(
left: TreeHistogram<V>,
right: TreeHistogram<V>,
): TreeHistogram<V> {
val bins = TreeMap<Double, Bin1D<Double, V>>().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<V>, value: Double): TreeHistogram<V> {
val bins = TreeMap<Double, Bin1D<Double, V>>().apply {
a.bins.forEach { bin ->
put(
bin.domain.center,
Bin1D(bin.domain, valueAlgebra.scale(bin.binValue, value))
)
}
}
return TreeHistogram(bins)
}
override fun TreeHistogram<V>.unaryMinus(): TreeHistogram<V> = this * (-1)
override val zero: TreeHistogram<V> = produce { }
}
///**
// * Build and fill a histogram with custom borders. Returns a read-only histogram.
// */
//public inline fun Histogram.custom(
// borders: DoubleArray,
// builder: Histogram1DBuilder<Double, Double>.() -> 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 <V : Any, A> Histogram.Companion.custom1D(
valueAlgebra: A,
borders: Buffer<Double>,
): TreeHistogramGroup<V, A> where A : Ring<V>, A : ScaleOperations<V> {
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)
}
}
}
}

View File

@ -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 <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(
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
}
@OptIn(UnstableKMathAPI::class)
@PublishedApi
internal class TreeHistogramBuilder(val binFactory: (Double) -> UnivariateDomain) : UnivariateHistogramBuilder {
internal class BinCounter(val domain: UnivariateDomain, val counter: Counter<Double> = Counter.double()) :
ClosedFloatingPointRange<Double> by domain.range
private val bins: TreeMap<Double, BinCounter> = 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<Double>, 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<Double, UnivariateBin>()) { (_, 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<UnivariateHistogram>, ScaleOperations<UnivariateHistogram> {
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<Double, UnivariateBin>().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<Double, UnivariateBin>().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)
}
}
}
}
}
}

View File

@ -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<Double>, ClosedFloatingPointRange<Double> by domain.range {
override val dimension: Int get() = 1
override fun contains(point: Buffer<Double>): Boolean = point.size == 1 && contains(point[0])
}
@OptIn(UnstableKMathAPI::class)
public interface UnivariateHistogram : Histogram<Double, UnivariateBin> {
public operator fun get(value: Double): UnivariateBin?
override operator fun get(point: Buffer<Double>): 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<Double> {
/**
* Thread safe put operation
*/
public fun putValue(at: Double, value: Double = 1.0)
override fun putValue(point: Buffer<Double>, value: Number)
}
@UnstableKMathAPI
public fun UnivariateHistogramBuilder.fill(items: Iterable<Double>): Unit = items.forEach(this::putValue)
@UnstableKMathAPI
public fun UnivariateHistogramBuilder.fill(array: DoubleArray): Unit = array.forEach(this::putValue)
@UnstableKMathAPI
public fun UnivariateHistogramBuilder.fill(buffer: Buffer<Double>): Unit = buffer.asSequence().forEach(this::putValue)

View File

@ -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())
}
}

View File

@ -27,7 +27,7 @@ public interface Distribution<T : Any> : Sampler<T> {
public companion object
}
public interface UnivariateDistribution<T : Comparable<T>> : Distribution<T> {
public interface Distribution1D<T : Comparable<T>> : Distribution<T> {
/**
* Cumulative distribution for ordered parameter (CDF)
*/
@ -37,7 +37,7 @@ public interface UnivariateDistribution<T : Comparable<T>> : Distribution<T> {
/**
* Compute probability integral in an interval
*/
public fun <T : Comparable<T>> UnivariateDistribution<T>.integral(from: T, to: T): Double {
public fun <T : Comparable<T>> Distribution1D<T>.integral(from: T, to: T): Double {
require(to > from)
return cumulative(to) - cumulative(from)
}

View File

@ -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<Double> {
public constructor(
mean: Double,
standardDeviation: Double,
normalized: NormalizedGaussianSampler = ZigguratNormalizedGaussianSampler,
) : this(GaussianSampler(mean, standardDeviation, normalized))
public class NormalDistribution(public val sampler: GaussianSampler) : Distribution1D<Double> {
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))

View File

@ -67,3 +67,12 @@ public fun Sampler<Double>.sampleBuffer(generator: RandomGenerator, size: Int):
@JvmName("sampleIntBuffer")
public fun Sampler<Int>.sampleBuffer(generator: RandomGenerator, size: Int): Chain<Buffer<Int>> =
sampleBuffer(generator, size, ::IntBuffer)
/**
* Samples a [Buffer] of values from this [Sampler].
*/
public suspend fun Sampler<Double>.nextBuffer(generator: RandomGenerator, size: Int): Buffer<Double> =
sampleBuffer(generator, size).first()
//TODO add `context(RandomGenerator) Sampler.nextBuffer

View File

@ -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<Double>) : UnivariateDistribution<Double> {
public class UniformDistribution(public val range: ClosedFloatingPointRange<Double>) : Distribution1D<Double> {
private val length: Double = range.endInclusive - range.start
override fun probability(arg: Double): Double = if (arg in range) 1.0 / length else 0.0