diff --git a/build.gradle.kts b/build.gradle.kts index 6ab33d31c..6d102a77a 100644 --- a/build.gradle.kts +++ b/build.gradle.kts @@ -1,8 +1,8 @@ plugins { - id("scientifik.publish") version "0.4.2" apply false + id("scientifik.publish") apply false } -val kmathVersion by extra("0.1.4-dev-4") +val kmathVersion by extra("0.1.4-dev-7") val bintrayRepo by extra("scientifik") val githubProject by extra("kmath") diff --git a/docs/codestyle.md b/docs/codestyle.md new file mode 100644 index 000000000..53789f7b2 --- /dev/null +++ b/docs/codestyle.md @@ -0,0 +1,22 @@ +# Local coding conventions + +Kmath and other `scientifik` projects use general [kotlin code conventions](https://kotlinlang.org/docs/reference/coding-conventions.html), but with a number of small changes and clarifications. + +## Utility class names +File name should coincide with a name of one of the classes contained in the file or start with small letter and describe its contents. + +The code convention [here](https://kotlinlang.org/docs/reference/coding-conventions.html#source-file-names) says that file names should start with capital letter even if file does not contain classes. Yet starting utility classes and aggregators with a small letter seems to be a good way to clearly visually separate those files. + +This convention could be changed in future in a non-breaking way. + +## Private variable names +Private variable names could start with underscore `_` in case the private mutable variable is shadowed by the public read-only value with the same meaning. + +Code convention do not permit underscores in names, but is is sometimes useful to "underscore" the fact that public and private versions define the same entity. It is allowed only for private variables. + +This convention could be changed in future in a non-breaking way. + +## Functions and properties one-liners +Use one-liners when they occupy single code window line both for functions and properties with getters like `val b: String get() = "fff"`. The same should be done with multiline expressions when they could be cleanly separated. + +There is not general consensus whenever use `fun a() = {}` or `fun a(){return}`. Yet from reader perspective one-lines seem to better show that the property or function is easily calculated. \ No newline at end of file diff --git a/examples/build.gradle.kts b/examples/build.gradle.kts index 8853b78a5..2fab47ac0 100644 --- a/examples/build.gradle.kts +++ b/examples/build.gradle.kts @@ -27,6 +27,7 @@ dependencies { implementation(project(":kmath-core")) implementation(project(":kmath-coroutines")) implementation(project(":kmath-commons")) + implementation(project(":kmath-prob")) implementation(project(":kmath-koma")) implementation(project(":kmath-viktor")) implementation(project(":kmath-dimensions")) @@ -57,6 +58,6 @@ benchmark { tasks.withType { kotlinOptions { - jvmTarget = Scientifik.JVM_VERSION + jvmTarget = Scientifik.JVM_TARGET.toString() } } \ No newline at end of file diff --git a/examples/src/main/kotlin/scientifik/kmath/commons/prob/DistributionBenchmark.kt b/examples/src/main/kotlin/scientifik/kmath/commons/prob/DistributionBenchmark.kt new file mode 100644 index 000000000..b060cddb6 --- /dev/null +++ b/examples/src/main/kotlin/scientifik/kmath/commons/prob/DistributionBenchmark.kt @@ -0,0 +1,71 @@ +package scientifik.kmath.commons.prob + +import kotlinx.coroutines.Dispatchers +import kotlinx.coroutines.async +import kotlinx.coroutines.runBlocking +import org.apache.commons.rng.sampling.distribution.ZigguratNormalizedGaussianSampler +import org.apache.commons.rng.simple.RandomSource +import scientifik.kmath.chains.BlockingRealChain +import scientifik.kmath.prob.* +import java.time.Duration +import java.time.Instant + + +private suspend fun runChain(): Duration { + val generator = RandomGenerator.fromSource(RandomSource.MT, 123L) + + val normal = Distribution.normal(NormalSamplerMethod.Ziggurat) + val chain = normal.sample(generator) as BlockingRealChain + + val startTime = Instant.now() + var sum = 0.0 + repeat(10000001) { counter -> + + sum += chain.nextDouble() + + if (counter % 100000 == 0) { + val duration = Duration.between(startTime, Instant.now()) + val meanValue = sum / counter + println("Chain sampler completed $counter elements in $duration: $meanValue") + } + } + return Duration.between(startTime, Instant.now()) +} + +private fun runDirect(): Duration { + val provider = RandomSource.create(RandomSource.MT, 123L) + val sampler = ZigguratNormalizedGaussianSampler(provider) + val startTime = Instant.now() + + var sum = 0.0 + repeat(10000001) { counter -> + + sum += sampler.sample() + + if (counter % 100000 == 0) { + val duration = Duration.between(startTime, Instant.now()) + val meanValue = sum / counter + println("Direct sampler completed $counter elements in $duration: $meanValue") + } + } + return Duration.between(startTime, Instant.now()) +} + +/** + * Comparing chain sampling performance with direct sampling performance + */ +fun main() { + runBlocking(Dispatchers.Default) { + val chainJob = async { + runChain() + } + + val directJob = async { + runDirect() + } + + println("Chain: ${chainJob.await()}") + println("Direct: ${directJob.await()}") + } + +} \ No newline at end of file diff --git a/examples/src/main/kotlin/scientifik/kmath/commons/prob/DistributionDemo.kt b/examples/src/main/kotlin/scientifik/kmath/commons/prob/DistributionDemo.kt index 3c5f53e13..e059415dc 100644 --- a/examples/src/main/kotlin/scientifik/kmath/commons/prob/DistributionDemo.kt +++ b/examples/src/main/kotlin/scientifik/kmath/commons/prob/DistributionDemo.kt @@ -5,10 +5,11 @@ import scientifik.kmath.chains.Chain import scientifik.kmath.chains.collectWithState import scientifik.kmath.prob.Distribution import scientifik.kmath.prob.RandomGenerator +import scientifik.kmath.prob.normal data class AveragingChainState(var num: Int = 0, var value: Double = 0.0) -fun Chain.mean(): Chain = collectWithState(AveragingChainState(),{it.copy()}){ chain-> +fun Chain.mean(): Chain = collectWithState(AveragingChainState(), { it.copy() }) { chain -> val next = chain.next() num++ value += next diff --git a/kmath-commons/src/main/kotlin/scientifik/kmath/commons/linear/CMMatrix.kt b/kmath-commons/src/main/kotlin/scientifik/kmath/commons/linear/CMMatrix.kt index 72e5fb95a..a17effccc 100644 --- a/kmath-commons/src/main/kotlin/scientifik/kmath/commons/linear/CMMatrix.kt +++ b/kmath-commons/src/main/kotlin/scientifik/kmath/commons/linear/CMMatrix.kt @@ -5,6 +5,7 @@ import org.apache.commons.math3.linear.RealMatrix import org.apache.commons.math3.linear.RealVector import scientifik.kmath.linear.* import scientifik.kmath.structures.Matrix +import scientifik.kmath.structures.NDStructure class CMMatrix(val origin: RealMatrix, features: Set? = null) : FeaturedMatrix { @@ -19,6 +20,16 @@ class CMMatrix(val origin: RealMatrix, features: Set? = null) : CMMatrix(origin, this.features + features) override fun get(i: Int, j: Int): Double = origin.getEntry(i, j) + + override fun equals(other: Any?): Boolean { + return NDStructure.equals(this, other as? NDStructure<*> ?: return false) + } + + override fun hashCode(): Int { + var result = origin.hashCode() + result = 31 * result + features.hashCode() + return result + } } fun Matrix.toCM(): CMMatrix = if (this is CMMatrix) { diff --git a/kmath-commons/src/main/kotlin/scientifik/kmath/commons/prob/CMRandomGeneratorWrapper.kt b/kmath-commons/src/main/kotlin/scientifik/kmath/commons/prob/CMRandomGeneratorWrapper.kt deleted file mode 100644 index 74e035ecb..000000000 --- a/kmath-commons/src/main/kotlin/scientifik/kmath/commons/prob/CMRandomGeneratorWrapper.kt +++ /dev/null @@ -1,32 +0,0 @@ -package scientifik.kmath.commons.prob - -import org.apache.commons.math3.random.JDKRandomGenerator -import scientifik.kmath.prob.RandomGenerator -import org.apache.commons.math3.random.RandomGenerator as CMRandom - -inline class CMRandomGeneratorWrapper(val generator: CMRandom) : RandomGenerator { - override fun nextDouble(): Double = generator.nextDouble() - - override fun nextInt(): Int = generator.nextInt() - - override fun nextLong(): Long = generator.nextLong() - - override fun nextBlock(size: Int): ByteArray = ByteArray(size).apply { generator.nextBytes(this) } - - override fun fork(): RandomGenerator { - TODO("not implemented") //To change body of created functions use File | Settings | File Templates. - } -} - -fun CMRandom.asKmathGenerator(): RandomGenerator = CMRandomGeneratorWrapper(this) - -fun RandomGenerator.asCMGenerator(): CMRandom = - (this as? CMRandomGeneratorWrapper)?.generator ?: TODO("Implement reverse CM wrapper") - -val RandomGenerator.Companion.default: RandomGenerator by lazy { JDKRandomGenerator().asKmathGenerator() } - -fun RandomGenerator.Companion.jdk(seed: Int? = null): RandomGenerator = if (seed == null) { - JDKRandomGenerator() -} else { - JDKRandomGenerator(seed) -}.asKmathGenerator() \ No newline at end of file diff --git a/kmath-commons/src/main/kotlin/scientifik/kmath/commons/prob/CommonsDistribution.kt b/kmath-commons/src/main/kotlin/scientifik/kmath/commons/prob/CommonsDistribution.kt deleted file mode 100644 index 94f8560a4..000000000 --- a/kmath-commons/src/main/kotlin/scientifik/kmath/commons/prob/CommonsDistribution.kt +++ /dev/null @@ -1,82 +0,0 @@ -package scientifik.kmath.commons.prob - -import org.apache.commons.math3.distribution.* -import scientifik.kmath.prob.Distribution -import scientifik.kmath.prob.RandomChain -import scientifik.kmath.prob.RandomGenerator -import scientifik.kmath.prob.UnivariateDistribution -import org.apache.commons.math3.random.RandomGenerator as CMRandom - -class CMRealDistributionWrapper(val builder: (CMRandom?) -> RealDistribution) : UnivariateDistribution { - - private val defaultDistribution by lazy { builder(null) } - - override fun probability(arg: Double): Double = defaultDistribution.probability(arg) - - override fun cumulative(arg: Double): Double = defaultDistribution.cumulativeProbability(arg) - - override fun sample(generator: RandomGenerator): RandomChain { - val distribution = builder(generator.asCMGenerator()) - return RandomChain(generator) { distribution.sample() } - } -} - -class CMIntDistributionWrapper(val builder: (CMRandom?) -> IntegerDistribution) : UnivariateDistribution { - - private val defaultDistribution by lazy { builder(null) } - - override fun probability(arg: Int): Double = defaultDistribution.probability(arg) - - override fun cumulative(arg: Int): Double = defaultDistribution.cumulativeProbability(arg) - - override fun sample(generator: RandomGenerator): RandomChain { - val distribution = builder(generator.asCMGenerator()) - return RandomChain(generator) { distribution.sample() } - } -} - - -fun Distribution.Companion.normal(mean: Double = 0.0, sigma: Double = 1.0): UnivariateDistribution = - CMRealDistributionWrapper { generator -> NormalDistribution(generator, mean, sigma) } - -fun Distribution.Companion.poisson(mean: Double): UnivariateDistribution = CMIntDistributionWrapper { generator -> - PoissonDistribution( - generator, - mean, - PoissonDistribution.DEFAULT_EPSILON, - PoissonDistribution.DEFAULT_MAX_ITERATIONS - ) -} - -fun Distribution.Companion.binomial(trials: Int, p: Double): UnivariateDistribution = - CMIntDistributionWrapper { generator -> - BinomialDistribution(generator, trials, p) - } - -fun Distribution.Companion.student(degreesOfFreedom: Double): UnivariateDistribution = - CMRealDistributionWrapper { generator -> - TDistribution(generator, degreesOfFreedom, TDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY) - } - -fun Distribution.Companion.chi2(degreesOfFreedom: Double): UnivariateDistribution = - CMRealDistributionWrapper { generator -> - ChiSquaredDistribution(generator, degreesOfFreedom) - } - -fun Distribution.Companion.fisher( - numeratorDegreesOfFreedom: Double, - denominatorDegreesOfFreedom: Double -): UnivariateDistribution = - CMRealDistributionWrapper { generator -> - FDistribution(generator, numeratorDegreesOfFreedom, denominatorDegreesOfFreedom) - } - -fun Distribution.Companion.exponential(mean: Double): UnivariateDistribution = - CMRealDistributionWrapper { generator -> - ExponentialDistribution(generator, mean) - } - -fun Distribution.Companion.uniform(a: Double, b: Double): UnivariateDistribution = - CMRealDistributionWrapper { generator -> - UniformRealDistribution(generator, a, b) - } \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/BufferMatrix.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/BufferMatrix.kt index 624cd44d4..73b18b810 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/BufferMatrix.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/BufferMatrix.kt @@ -27,7 +27,7 @@ class BufferMatrixContext>( @Suppress("OVERRIDE_BY_INLINE") object RealMatrixContext : GenericMatrixContext { - override val elementContext = RealField + override val elementContext get() = RealField override inline fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> Double): Matrix { val buffer = DoubleBuffer(rows * columns) { offset -> initializer(offset / columns, offset % columns) } @@ -101,8 +101,15 @@ infix fun BufferMatrix.dot(other: BufferMatrix): BufferMatrix.unsafeArray(): DoubleArray = if (this is DoubleBuffer) { + array + } else { + DoubleArray(size) { get(it) } + } + + val a = this.buffer.unsafeArray() + val b = other.buffer.unsafeArray() for (i in (0 until rowNum)) { for (j in (0 until other.colNum)) { diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgebra.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgebra.kt index 0456ffebb..abd8603f4 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgebra.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgebra.kt @@ -1,12 +1,8 @@ package scientifik.kmath.linear -import scientifik.kmath.operations.Field -import scientifik.kmath.operations.Norm -import scientifik.kmath.operations.RealField import scientifik.kmath.structures.Buffer import scientifik.kmath.structures.Matrix import scientifik.kmath.structures.VirtualBuffer -import scientifik.kmath.structures.asSequence typealias Point = Buffer @@ -19,8 +15,6 @@ interface LinearSolver { fun inverse(a: Matrix): Matrix } -typealias RealMatrix = Matrix - /** * Convert matrix to vector if it is possible */ diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/MatrixBuilder.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/MatrixBuilder.kt index 466dbea6e..516f65bb8 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/MatrixBuilder.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/MatrixBuilder.kt @@ -1,14 +1,46 @@ package scientifik.kmath.linear +import scientifik.kmath.structures.Buffer +import scientifik.kmath.structures.BufferFactory import scientifik.kmath.structures.Structure2D import scientifik.kmath.structures.asBuffer -class MatrixBuilder(val rows: Int, val columns: Int) { - operator fun invoke(vararg elements: T): FeaturedMatrix { +class MatrixBuilder(val rows: Int, val columns: Int) { + operator fun invoke(vararg elements: T): FeaturedMatrix { if (rows * columns != elements.size) error("The number of elements ${elements.size} is not equal $rows * $columns") val buffer = elements.asBuffer() return BufferMatrix(rows, columns, buffer) } + + //TODO add specific matrix builder functions like diagonal, etc } -fun Structure2D.Companion.build(rows: Int, columns: Int): MatrixBuilder = MatrixBuilder(rows, columns) \ No newline at end of file +fun Structure2D.Companion.build(rows: Int, columns: Int): MatrixBuilder = MatrixBuilder(rows, columns) + +fun Structure2D.Companion.row(vararg values: T): FeaturedMatrix { + val buffer = values.asBuffer() + return BufferMatrix(1, values.size, buffer) +} + +inline fun Structure2D.Companion.row( + size: Int, + factory: BufferFactory = Buffer.Companion::auto, + noinline builder: (Int) -> T +): FeaturedMatrix { + val buffer = factory(size, builder) + return BufferMatrix(1, size, buffer) +} + +fun Structure2D.Companion.column(vararg values: T): FeaturedMatrix { + val buffer = values.asBuffer() + return BufferMatrix(values.size, 1, buffer) +} + +inline fun Structure2D.Companion.column( + size: Int, + factory: BufferFactory = Buffer.Companion::auto, + noinline builder: (Int) -> T +): FeaturedMatrix { + val buffer = factory(size, builder) + return BufferMatrix(size, 1, buffer) +} diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/BigInt.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/BigInt.kt index 1661170d3..5e7d49b1c 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/BigInt.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/BigInt.kt @@ -2,6 +2,7 @@ package scientifik.kmath.operations import scientifik.kmath.operations.BigInt.Companion.BASE import scientifik.kmath.operations.BigInt.Companion.BASE_SIZE +import scientifik.kmath.structures.* import kotlin.math.log2 import kotlin.math.max import kotlin.math.min @@ -482,3 +483,18 @@ fun String.parseBigInteger(): BigInt? { } return res * sign } + +inline fun Buffer.Companion.bigInt(size: Int, initializer: (Int) -> BigInt): Buffer = + boxing(size, initializer) + +inline fun MutableBuffer.Companion.bigInt(size: Int, initializer: (Int) -> BigInt): MutableBuffer = + boxing(size, initializer) + +fun NDAlgebra.Companion.bigInt(vararg shape: Int): BoxingNDRing = + BoxingNDRing(shape, BigIntField, Buffer.Companion::bigInt) + +fun NDElement.Companion.bigInt( + vararg shape: Int, + initializer: BigIntField.(IntArray) -> BigInt +): BufferedNDRingElement = + NDAlgebra.bigInt(*shape).produce(initializer) \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/NumberAlgebra.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/NumberAlgebra.kt index 9639e4c28..2b6a92f14 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/NumberAlgebra.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/NumberAlgebra.kt @@ -102,7 +102,7 @@ object IntRing : Ring, Norm { override val zero: Int = 0 override inline fun add(a: Int, b: Int) = a + b override inline fun multiply(a: Int, b: Int) = a * b - override inline fun multiply(a: Int, k: Number) = (k * a) + override inline fun multiply(a: Int, k: Number) = k.toInt() * a override val one: Int = 1 override inline fun norm(arg: Int) = abs(arg) @@ -124,7 +124,7 @@ object ShortRing : Ring, Norm { override val zero: Short = 0 override inline fun add(a: Short, b: Short) = (a + b).toShort() override inline fun multiply(a: Short, b: Short) = (a * b).toShort() - override inline fun multiply(a: Short, k: Number) = (a * k) + override inline fun multiply(a: Short, k: Number) = (a * k.toShort()).toShort() override val one: Short = 1 override fun norm(arg: Short): Short = if (arg > 0) arg else (-arg).toShort() @@ -146,7 +146,7 @@ object ByteRing : Ring, Norm { override val zero: Byte = 0 override inline fun add(a: Byte, b: Byte) = (a + b).toByte() override inline fun multiply(a: Byte, b: Byte) = (a * b).toByte() - override inline fun multiply(a: Byte, k: Number) = (a * k) + override inline fun multiply(a: Byte, k: Number) = (a * k.toByte()).toByte() override val one: Byte = 1 override fun norm(arg: Byte): Byte = if (arg > 0) arg else (-arg).toByte() @@ -168,7 +168,7 @@ object LongRing : Ring, Norm { override val zero: Long = 0 override inline fun add(a: Long, b: Long) = (a + b) override inline fun multiply(a: Long, b: Long) = (a * b) - override inline fun multiply(a: Long, k: Number) = (a * k) + override inline fun multiply(a: Long, k: Number) = a * k.toLong() override val one: Long = 1 override fun norm(arg: Long): Long = abs(arg) @@ -180,4 +180,4 @@ object LongRing : Ring, Norm { override inline fun Long.minus(b: Long) = (this - b) override inline fun Long.times(b: Long) = (this * b) -} \ No newline at end of file +} diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/BufferedNDElement.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/BufferedNDElement.kt index 04049368a..be80f66bf 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/BufferedNDElement.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/BufferedNDElement.kt @@ -3,10 +3,10 @@ package scientifik.kmath.structures import scientifik.kmath.operations.* /** - * Base interface for an element with context, containing strides + * Base class for an element with context, containing strides */ -interface BufferedNDElement : NDBuffer, NDElement> { - override val context: BufferedNDAlgebra +abstract class BufferedNDElement : NDBuffer(), NDElement> { + abstract override val context: BufferedNDAlgebra override val strides get() = context.strides @@ -16,7 +16,7 @@ interface BufferedNDElement : NDBuffer, NDElement> { class BufferedNDSpaceElement>( override val context: BufferedNDSpace, override val buffer: Buffer -) : BufferedNDElement, SpaceElement, BufferedNDSpaceElement, BufferedNDSpace> { +) : BufferedNDElement(), SpaceElement, BufferedNDSpaceElement, BufferedNDSpace> { override fun unwrap(): NDBuffer = this @@ -29,7 +29,7 @@ class BufferedNDSpaceElement>( class BufferedNDRingElement>( override val context: BufferedNDRing, override val buffer: Buffer -) : BufferedNDElement, RingElement, BufferedNDRingElement, BufferedNDRing> { +) : BufferedNDElement(), RingElement, BufferedNDRingElement, BufferedNDRing> { override fun unwrap(): NDBuffer = this @@ -42,7 +42,7 @@ class BufferedNDRingElement>( class BufferedNDFieldElement>( override val context: BufferedNDField, override val buffer: Buffer -) : BufferedNDElement, FieldElement, BufferedNDFieldElement, BufferedNDField> { +) : BufferedNDElement(), FieldElement, BufferedNDFieldElement, BufferedNDField> { override fun unwrap(): NDBuffer = this diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt index f02fd8dd0..613a0d7ca 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt @@ -71,7 +71,7 @@ interface Buffer { fun Buffer.asSequence(): Sequence = Sequence(::iterator) -fun Buffer.asIterable(): Iterable = asSequence().asIterable() +fun Buffer.asIterable(): Iterable = Iterable(::iterator) val Buffer<*>.indices: IntRange get() = IntRange(0, size - 1) @@ -161,36 +161,7 @@ class ArrayBuffer(private val array: Array) : MutableBuffer { override fun copy(): MutableBuffer = ArrayBuffer(array.copyOf()) } -fun Array.asBuffer() = ArrayBuffer(this) - -inline class DoubleBuffer(val array: DoubleArray) : MutableBuffer { - override val size: Int get() = array.size - - override fun get(index: Int): Double = array[index] - - override fun set(index: Int, value: Double) { - array[index] = value - } - - override fun iterator() = array.iterator() - - override fun copy(): MutableBuffer = DoubleBuffer(array.copyOf()) -} - -@Suppress("FunctionName") -inline fun DoubleBuffer(size: Int, init: (Int) -> Double) = DoubleBuffer(DoubleArray(size) { init(it) }) - -/** - * Transform buffer of doubles into array for high performance operations - */ -val Buffer.array: DoubleArray - get() = if (this is DoubleBuffer) { - array - } else { - DoubleArray(size) { get(it) } - } - -fun DoubleArray.asBuffer() = DoubleBuffer(this) +fun Array.asBuffer(): ArrayBuffer = ArrayBuffer(this) inline class ShortBuffer(val array: ShortArray) : MutableBuffer { override val size: Int get() = array.size diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/DoubleBuffer.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/DoubleBuffer.kt new file mode 100644 index 000000000..c0b7f713b --- /dev/null +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/DoubleBuffer.kt @@ -0,0 +1,34 @@ +package scientifik.kmath.structures + +inline class DoubleBuffer(val array: DoubleArray) : MutableBuffer { + override val size: Int get() = array.size + + override fun get(index: Int): Double = array[index] + + override fun set(index: Int, value: Double) { + array[index] = value + } + + override fun iterator() = array.iterator() + + override fun copy(): MutableBuffer = + DoubleBuffer(array.copyOf()) +} + +@Suppress("FunctionName") +inline fun DoubleBuffer(size: Int, init: (Int) -> Double): DoubleBuffer = DoubleBuffer(DoubleArray(size) { init(it) }) + +@Suppress("FunctionName") +fun DoubleBuffer(vararg doubles: Double): DoubleBuffer = DoubleBuffer(doubles) + +/** + * Transform buffer of doubles into array for high performance operations + */ +val MutableBuffer.array: DoubleArray + get() = if (this is DoubleBuffer) { + array + } else { + DoubleArray(size) { get(it) } + } + +fun DoubleArray.asBuffer() = DoubleBuffer(this) \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDStructure.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDStructure.kt index 808f970c5..236a6f366 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDStructure.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDStructure.kt @@ -14,15 +14,25 @@ interface NDStructure { fun elements(): Sequence> + override fun equals(other: Any?): Boolean + + override fun hashCode(): Int + companion object { fun equals(st1: NDStructure<*>, st2: NDStructure<*>): Boolean { - return when { - st1 === st2 -> true - st1 is BufferNDStructure<*> && st2 is BufferNDStructure<*> && st1.strides == st2.strides -> st1.buffer.contentEquals( - st2.buffer - ) - else -> st1.elements().all { (index, value) -> value == st2[index] } + if (st1 === st2) return true + + // fast comparison of buffers if possible + if ( + st1 is NDBuffer && + st2 is NDBuffer && + st1.strides == st2.strides + ) { + return st1.buffer.contentEquals(st2.buffer) } + + //element by element comparison if it could not be avoided + return st1.elements().all { (index, value) -> value == st2[index] } } /** @@ -177,15 +187,25 @@ class DefaultStrides private constructor(override val shape: IntArray) : Strides } } -interface NDBuffer : NDStructure { - val buffer: Buffer - val strides: Strides +abstract class NDBuffer : NDStructure { + abstract val buffer: Buffer + abstract val strides: Strides override fun get(index: IntArray): T = buffer[strides.offset(index)] override val shape: IntArray get() = strides.shape - override fun elements() = strides.indices().map { it to this[it] } + override fun elements(): Sequence> = strides.indices().map { it to this[it] } + + override fun equals(other: Any?): Boolean { + return NDStructure.equals(this, other as? NDStructure<*> ?: return false) + } + + override fun hashCode(): Int { + var result = strides.hashCode() + result = 31 * result + buffer.hashCode() + return result + } } /** @@ -194,34 +214,12 @@ interface NDBuffer : NDStructure { class BufferNDStructure( override val strides: Strides, override val buffer: Buffer -) : NDBuffer { - +) : NDBuffer() { init { if (strides.linearSize != buffer.size) { error("Expected buffer side of ${strides.linearSize}, but found ${buffer.size}") } } - - override fun get(index: IntArray): T = buffer[strides.offset(index)] - - override val shape: IntArray get() = strides.shape - - override fun elements() = strides.indices().map { it to this[it] } - - override fun equals(other: Any?): Boolean { - return when { - this === other -> true - other is BufferNDStructure<*> && this.strides == other.strides -> this.buffer.contentEquals(other.buffer) - other is NDStructure<*> -> elements().all { (index, value) -> value == other[index] } - else -> false - } - } - - override fun hashCode(): Int { - var result = strides.hashCode() - result = 31 * result + buffer.hashCode() - return result - } } /** @@ -245,7 +243,7 @@ inline fun NDStructure.mapToBuffer( class MutableBufferNDStructure( override val strides: Strides, override val buffer: MutableBuffer -) : NDBuffer, MutableNDStructure { +) : NDBuffer(), MutableNDStructure { init { if (strides.linearSize != buffer.size) { diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Structure1D.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Structure1D.kt index df56017a3..9974538ec 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Structure1D.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Structure1D.kt @@ -17,7 +17,7 @@ interface Structure1D : NDStructure, Buffer { /** * A 1D wrapper for nd-structure */ -private inline class Structure1DWrapper(val structure: NDStructure) : Structure1D { +private inline class Structure1DWrapper(val structure: NDStructure) : Structure1D{ override val shape: IntArray get() = structure.shape override val size: Int get() = structure.shape[0] diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Structure2D.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Structure2D.kt index e736f84a0..bfaf8338d 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Structure2D.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Structure2D.kt @@ -14,7 +14,6 @@ interface Structure2D : NDStructure { return get(index[0], index[1]) } - val rows: Buffer> get() = VirtualBuffer(rowNum) { i -> VirtualBuffer(colNum) { j -> get(i, j) } @@ -58,22 +57,4 @@ fun NDStructure.as2D(): Structure2D = if (shape.size == 2) { error("Can't create 2d-structure from ${shape.size}d-structure") } -/** - * Represent this 2D structure as 1D if it has exactly one column. Throw error otherwise. - */ -fun Structure2D.as1D() = if (colNum == 1) { - object : Structure1D { - override fun get(index: Int): T = get(index, 0) - - override val shape: IntArray get() = intArrayOf(rowNum) - - override fun elements(): Sequence> = elements() - - override val size: Int get() = rowNum - } -} else { - error("Can't convert matrix with more than one column to vector") -} - - typealias Matrix = Structure2D \ No newline at end of file diff --git a/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/MatrixTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/MatrixTest.kt index 7d1209963..dcd510e32 100644 --- a/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/MatrixTest.kt +++ b/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/MatrixTest.kt @@ -17,7 +17,7 @@ class MatrixTest { @Test fun testBuilder() { - val matrix = Matrix.build(2, 3)( + val matrix = Matrix.build(2, 3)( 1.0, 0.0, 0.0, 0.0, 1.0, 2.0 ) diff --git a/kmath-core/src/jvmMain/kotlin/scientifik/kmath/operations/BigNumbers.kt b/kmath-core/src/jvmMain/kotlin/scientifik/kmath/operations/bigNumbers.kt similarity index 56% rename from kmath-core/src/jvmMain/kotlin/scientifik/kmath/operations/BigNumbers.kt rename to kmath-core/src/jvmMain/kotlin/scientifik/kmath/operations/bigNumbers.kt index f44d15042..76ca199c5 100644 --- a/kmath-core/src/jvmMain/kotlin/scientifik/kmath/operations/BigNumbers.kt +++ b/kmath-core/src/jvmMain/kotlin/scientifik/kmath/operations/bigNumbers.kt @@ -5,7 +5,7 @@ import java.math.BigDecimal import java.math.BigInteger import java.math.MathContext -object BigIntegerRing : Ring { +object JBigIntegerField : Field { override val zero: BigInteger = BigInteger.ZERO override val one: BigInteger = BigInteger.ONE @@ -14,9 +14,11 @@ object BigIntegerRing : Ring { override fun multiply(a: BigInteger, k: Number): BigInteger = a.multiply(k.toInt().toBigInteger()) override fun multiply(a: BigInteger, b: BigInteger): BigInteger = a.multiply(b) + + override fun divide(a: BigInteger, b: BigInteger): BigInteger = a.div(b) } -class BigDecimalField(val mathContext: MathContext = MathContext.DECIMAL64) : Field { +class JBigDecimalField(val mathContext: MathContext = MathContext.DECIMAL64) : Field { override val zero: BigDecimal = BigDecimal.ZERO override val one: BigDecimal = BigDecimal.ONE @@ -27,19 +29,4 @@ class BigDecimalField(val mathContext: MathContext = MathContext.DECIMAL64) : Fi override fun multiply(a: BigDecimal, b: BigDecimal): BigDecimal = a.multiply(b, mathContext) override fun divide(a: BigDecimal, b: BigDecimal): BigDecimal = a.divide(b, mathContext) -} - -inline fun Buffer.Companion.bigInt(size: Int, initializer: (Int) -> BigInteger): Buffer = - boxing(size, initializer) - -inline fun MutableBuffer.Companion.bigInt(size: Int, initializer: (Int) -> BigInteger): MutableBuffer = - boxing(size, initializer) - -fun NDAlgebra.Companion.bigInt(vararg shape: Int): BoxingNDRing = - BoxingNDRing(shape, BigIntegerRing, Buffer.Companion::bigInt) - -fun NDElement.Companion.bigInt( - vararg shape: Int, - initializer: BigIntegerRing.(IntArray) -> BigInteger -): BufferedNDRingElement = - NDAlgebra.bigInt(*shape).produce(initializer) \ No newline at end of file +} \ No newline at end of file diff --git a/kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/chains/BlockingIntChain.kt b/kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/chains/BlockingIntChain.kt new file mode 100644 index 000000000..6ec84d5c7 --- /dev/null +++ b/kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/chains/BlockingIntChain.kt @@ -0,0 +1,12 @@ +package scientifik.kmath.chains + +/** + * Performance optimized chain for integer values + */ +abstract class BlockingIntChain : Chain { + abstract fun nextInt(): Int + + override suspend fun next(): Int = nextInt() + + fun nextBlock(size: Int): IntArray = IntArray(size) { nextInt() } +} \ No newline at end of file diff --git a/kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/chains/BlockingRealChain.kt b/kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/chains/BlockingRealChain.kt new file mode 100644 index 000000000..6b69d2734 --- /dev/null +++ b/kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/chains/BlockingRealChain.kt @@ -0,0 +1,12 @@ +package scientifik.kmath.chains + +/** + * Performance optimized chain for real values + */ +abstract class BlockingRealChain : Chain { + abstract fun nextDouble(): Double + + override suspend fun next(): Double = nextDouble() + + fun nextBlock(size: Int): DoubleArray = DoubleArray(size) { nextDouble() } +} \ No newline at end of file diff --git a/kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/chains/Chain.kt b/kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/chains/Chain.kt index 1c2872d17..5635499e5 100644 --- a/kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/chains/Chain.kt +++ b/kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/chains/Chain.kt @@ -38,13 +38,16 @@ interface Chain: Flow { */ fun fork(): Chain - @InternalCoroutinesApi + @OptIn(InternalCoroutinesApi::class) override suspend fun collect(collector: FlowCollector) { - kotlinx.coroutines.flow.flow { while (true) emit(next()) }.collect(collector) + kotlinx.coroutines.flow.flow { + while (true){ + emit(next()) + } + }.collect(collector) } companion object - } @@ -68,6 +71,8 @@ class MarkovChain(private val seed: suspend () -> R, private val ge private var value: R? = null + fun value() = value + override suspend fun next(): R { mutex.withLock { val newValue = gen(value ?: seed()) @@ -97,6 +102,8 @@ class StatefulChain( private var value: R? = null + fun value() = value + override suspend fun next(): R { mutex.withLock { val newValue = state.gen(value ?: state.seed()) diff --git a/kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/streaming/BufferFlow.kt b/kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/streaming/BufferFlow.kt index cf4d4cc17..bef21a680 100644 --- a/kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/streaming/BufferFlow.kt +++ b/kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/streaming/BufferFlow.kt @@ -2,9 +2,11 @@ package scientifik.kmath.streaming import kotlinx.coroutines.FlowPreview import kotlinx.coroutines.flow.* +import scientifik.kmath.chains.BlockingRealChain import scientifik.kmath.structures.Buffer import scientifik.kmath.structures.BufferFactory import scientifik.kmath.structures.DoubleBuffer +import scientifik.kmath.structures.asBuffer /** * Create a [Flow] from buffer @@ -45,20 +47,28 @@ fun Flow.chunked(bufferSize: Int, bufferFactory: BufferFactory): Flow< */ fun Flow.chunked(bufferSize: Int): Flow = flow { require(bufferSize > 0) { "Resulting chunk size must be more than zero" } - val array = DoubleArray(bufferSize) - var counter = 0 - this@chunked.collect { element -> - array[counter] = element - counter++ - if (counter == bufferSize) { - val buffer = DoubleBuffer(array) - emit(buffer) - counter = 0 + if (this@chunked is BlockingRealChain) { + //performance optimization for blocking primitive chain + while (true) { + emit(nextBlock(bufferSize).asBuffer()) + } + } else { + val array = DoubleArray(bufferSize) + var counter = 0 + + this@chunked.collect { element -> + array[counter] = element + counter++ + if (counter == bufferSize) { + val buffer = DoubleBuffer(array) + emit(buffer) + counter = 0 + } + } + if (counter > 0) { + emit(DoubleBuffer(counter) { array[it] }) } - } - if (counter > 0) { - emit(DoubleBuffer(counter) { array[it] }) } } diff --git a/kmath-coroutines/src/jvmMain/kotlin/scientifik/kmath/structures/LazyNDStructure.kt b/kmath-coroutines/src/jvmMain/kotlin/scientifik/kmath/structures/LazyNDStructure.kt index 784b7cd10..e65940739 100644 --- a/kmath-coroutines/src/jvmMain/kotlin/scientifik/kmath/structures/LazyNDStructure.kt +++ b/kmath-coroutines/src/jvmMain/kotlin/scientifik/kmath/structures/LazyNDStructure.kt @@ -30,6 +30,20 @@ class LazyNDStructure( } return res.asSequence() } + + override fun equals(other: Any?): Boolean { + return NDStructure.equals(this, other as? NDStructure<*> ?: return false) + } + + override fun hashCode(): Int { + var result = scope.hashCode() + result = 31 * result + shape.contentHashCode() + result = 31 * result + function.hashCode() + result = 31 * result + cache.hashCode() + return result + } + + } fun NDStructure.deferred(index: IntArray) = diff --git a/kmath-for-real/src/commonMain/kotlin/scientifik/kmath/linear/RealVector.kt b/kmath-for-real/src/commonMain/kotlin/scientifik/kmath/linear/RealVector.kt deleted file mode 100644 index d3cf07e79..000000000 --- a/kmath-for-real/src/commonMain/kotlin/scientifik/kmath/linear/RealVector.kt +++ /dev/null @@ -1,48 +0,0 @@ -package scientifik.kmath.linear - -import scientifik.kmath.operations.Field -import scientifik.kmath.operations.Norm -import scientifik.kmath.operations.RealField -import scientifik.kmath.operations.SpaceElement -import scientifik.kmath.structures.Buffer -import scientifik.kmath.structures.DoubleBuffer -import scientifik.kmath.structures.asBuffer -import scientifik.kmath.structures.asSequence - -fun DoubleArray.asVector() = RealVector(this.asBuffer()) -fun List.asVector() = RealVector(this.asBuffer()) - - -object VectorL2Norm : Norm, Double> { - override fun norm(arg: Point): Double = - kotlin.math.sqrt(arg.asSequence().sumByDouble { it.toDouble() }) -} - -inline class RealVector(val point: Point) : - SpaceElement, RealVector, VectorSpace>, Point { - override val context: VectorSpace get() = space(point.size) - - override fun unwrap(): Point = point - - override fun Point.wrap(): RealVector = RealVector(this) - - override val size: Int get() = point.size - - override fun get(index: Int): Double = point[index] - - override fun iterator(): Iterator = point.iterator() - - companion object { - - private val spaceCache = HashMap>() - - inline operator fun invoke(dim:Int, initalizer: (Int)-> Double) = RealVector(DoubleBuffer(dim, initalizer)) - - operator fun invoke(vararg values: Double) = values.asVector() - - fun space(dim: Int) = - spaceCache.getOrPut(dim) { - BufferVectorSpace(dim, RealField) { size, init -> Buffer.real(size, init) } - } - } -} \ No newline at end of file diff --git a/kmath-for-real/src/commonMain/kotlin/scientifik/kmath/real/DoubleMatrixOperations.kt b/kmath-for-real/src/commonMain/kotlin/scientifik/kmath/real/DoubleMatrixOperations.kt deleted file mode 100644 index 7eeba3031..000000000 --- a/kmath-for-real/src/commonMain/kotlin/scientifik/kmath/real/DoubleMatrixOperations.kt +++ /dev/null @@ -1,146 +0,0 @@ -package scientifik.kmath.real - -import scientifik.kmath.linear.MatrixContext -import scientifik.kmath.linear.RealMatrixContext.elementContext -import scientifik.kmath.linear.VirtualMatrix -import scientifik.kmath.operations.sum -import scientifik.kmath.structures.Buffer -import scientifik.kmath.structures.Matrix -import scientifik.kmath.structures.asSequence -import kotlin.math.pow - -/* - * Functions for convenient "numpy-like" operations with Double matrices. - * - * Initial implementation of these functions is taken from: - * https://github.com/thomasnield/numky/blob/master/src/main/kotlin/org/nield/numky/linear/DoubleOperators.kt - * - */ - -/* - * Functions that help create a real (Double) matrix - */ - -fun realMatrix(rowNum: Int, colNum: Int, initializer: (i: Int, j: Int) -> Double) = - MatrixContext.real.produce(rowNum, colNum, initializer) - -fun Sequence.toMatrix() = toList().let { - MatrixContext.real.produce(it.size, it[0].size) { row, col -> it[row][col] } -} - -fun Matrix.repeatStackVertical(n: Int) = VirtualMatrix(rowNum*n, colNum) { - row, col -> get(if (row == 0) 0 else row % rowNum, col) -} - -/* - * Operations for matrix and real number - */ - -operator fun Matrix.times(double: Double) = MatrixContext.real.produce(rowNum, colNum) { - row, col -> this[row, col] * double -} - -operator fun Matrix.plus(double: Double) = MatrixContext.real.produce(rowNum, colNum) { - row, col -> this[row, col] + double -} - -operator fun Matrix.minus(double: Double) = MatrixContext.real.produce(rowNum, colNum) { - row, col -> this[row, col] - double -} - -operator fun Matrix.div(double: Double) = MatrixContext.real.produce(rowNum, colNum) { - row, col -> this[row, col] / double -} - -operator fun Double.times(matrix: Matrix) = MatrixContext.real.produce(matrix.rowNum, matrix.colNum) { - row, col -> this * matrix[row, col] -} - -operator fun Double.plus(matrix: Matrix) = MatrixContext.real.produce(matrix.rowNum, matrix.colNum) { - row, col -> this + matrix[row, col] -} - -operator fun Double.minus(matrix: Matrix) = MatrixContext.real.produce(matrix.rowNum, matrix.colNum) { - row, col -> this - matrix[row, col] -} - -// TODO: does this operation make sense? Should it be 'this/matrix[row, col]'? -//operator fun Double.div(matrix: Matrix) = MatrixContext.real.produce(matrix.rowNum, matrix.colNum) { -// row, col -> matrix[row, col] / this -//} - -/* - * Per-element (!) square and power operations - */ - -fun Matrix.square() = MatrixContext.real.produce(rowNum, colNum) { - row, col -> this[row, col].pow(2) -} - -fun Matrix.pow(n: Int) = MatrixContext.real.produce(rowNum, colNum) { - i, j -> this[i, j].pow(n) -} - -/* - * Operations on two matrices (per-element!) - */ - -operator fun Matrix.times(other: Matrix) = MatrixContext.real.produce(rowNum, colNum) { - row, col -> this[row, col] * other[row, col] -} - -operator fun Matrix.plus(other: Matrix) = MatrixContext.real.add(this, other) - -operator fun Matrix.minus(other: Matrix) = MatrixContext.real.produce(rowNum, colNum) { - row, col -> this[row,col] - other[row,col] -} - -/* - * Operations on columns - */ - -inline fun Matrix.appendColumn(crossinline mapper: (Buffer) -> Double) = - MatrixContext.real.produce(rowNum,colNum+1) { - row, col -> - if (col < colNum) - this[row, col] - else - mapper(rows[row]) - } - -fun Matrix.extractColumns(columnRange: IntRange) = MatrixContext.real.produce(rowNum, columnRange.count()) { - row, col -> this[row, columnRange.first + col] -} - -fun Matrix.extractColumn(columnIndex: Int) = extractColumns(columnIndex..columnIndex) - -fun Matrix.sumByColumn() = MatrixContext.real.produce(1, colNum) { _, j -> - val column = columns[j] - with(elementContext) { - sum(column.asSequence()) - } -} - -fun Matrix.minByColumn() = MatrixContext.real.produce(1, colNum) { - _, j -> columns[j].asSequence().min() ?: throw Exception("Cannot produce min on empty column") -} - -fun Matrix.maxByColumn() = MatrixContext.real.produce(1, colNum) { - _, j -> columns[j].asSequence().max() ?: throw Exception("Cannot produce min on empty column") -} - -fun Matrix.averageByColumn() = MatrixContext.real.produce(1, colNum) { - _, j -> columns[j].asSequence().average() -} - -/* - * Operations processing all elements - */ - -fun Matrix.sum() = elements().map { (_, value) -> value }.sum() - -fun Matrix.min() = elements().map { (_, value) -> value }.min() - -fun Matrix.max() = elements().map { (_, value) -> value }.max() - -fun Matrix.average() = elements().map { (_, value) -> value }.average() diff --git a/kmath-for-real/src/commonMain/kotlin/scientifik/kmath/real/RealVector.kt b/kmath-for-real/src/commonMain/kotlin/scientifik/kmath/real/RealVector.kt new file mode 100644 index 000000000..ff4c835ed --- /dev/null +++ b/kmath-for-real/src/commonMain/kotlin/scientifik/kmath/real/RealVector.kt @@ -0,0 +1,59 @@ +package scientifik.kmath.real + +import scientifik.kmath.linear.BufferVectorSpace +import scientifik.kmath.linear.Point +import scientifik.kmath.linear.VectorSpace +import scientifik.kmath.operations.Norm +import scientifik.kmath.operations.RealField +import scientifik.kmath.operations.SpaceElement +import scientifik.kmath.structures.Buffer +import scientifik.kmath.structures.DoubleBuffer +import scientifik.kmath.structures.asBuffer +import scientifik.kmath.structures.asIterable +import kotlin.math.sqrt + +fun DoubleArray.asVector() = RealVector(this.asBuffer()) +fun List.asVector() = RealVector(this.asBuffer()) + + +object VectorL2Norm : Norm, Double> { + override fun norm(arg: Point): Double = sqrt(arg.asIterable().sumByDouble { it.toDouble() }) +} + +inline class RealVector(private val point: Point) : + SpaceElement, RealVector, VectorSpace>, Point { + + override val context: VectorSpace + get() = space( + point.size + ) + + override fun unwrap(): Point = point + + override fun Point.wrap(): RealVector = + RealVector(this) + + override val size: Int get() = point.size + + override fun get(index: Int): Double = point[index] + + override fun iterator(): Iterator = point.iterator() + + companion object { + + private val spaceCache = HashMap>() + + inline operator fun invoke(dim: Int, initializer: (Int) -> Double) = + RealVector(DoubleBuffer(dim, initializer)) + + operator fun invoke(vararg values: Double): RealVector = values.asVector() + + fun space(dim: Int): BufferVectorSpace = + spaceCache.getOrPut(dim) { + BufferVectorSpace( + dim, + RealField + ) { size, init -> Buffer.real(size, init) } + } + } +} \ No newline at end of file diff --git a/kmath-for-real/src/commonMain/kotlin/scientifik/kmath/real/realBuffer.kt b/kmath-for-real/src/commonMain/kotlin/scientifik/kmath/real/realBuffer.kt new file mode 100644 index 000000000..d9ee4d90b --- /dev/null +++ b/kmath-for-real/src/commonMain/kotlin/scientifik/kmath/real/realBuffer.kt @@ -0,0 +1,8 @@ +package scientifik.kmath.real + +import scientifik.kmath.structures.DoubleBuffer + +/** + * Simplified [DoubleBuffer] to array comparison + */ +fun DoubleBuffer.contentEquals(vararg doubles: Double) = array.contentEquals(doubles) \ No newline at end of file diff --git a/kmath-for-real/src/commonMain/kotlin/scientifik/kmath/real/realMatrix.kt b/kmath-for-real/src/commonMain/kotlin/scientifik/kmath/real/realMatrix.kt new file mode 100644 index 000000000..813d89577 --- /dev/null +++ b/kmath-for-real/src/commonMain/kotlin/scientifik/kmath/real/realMatrix.kt @@ -0,0 +1,161 @@ +package scientifik.kmath.real + +import scientifik.kmath.linear.MatrixContext +import scientifik.kmath.linear.RealMatrixContext.elementContext +import scientifik.kmath.linear.VirtualMatrix +import scientifik.kmath.operations.sum +import scientifik.kmath.structures.Buffer +import scientifik.kmath.structures.DoubleBuffer +import scientifik.kmath.structures.Matrix +import scientifik.kmath.structures.asIterable +import kotlin.math.pow + +/* + * Functions for convenient "numpy-like" operations with Double matrices. + * + * Initial implementation of these functions is taken from: + * https://github.com/thomasnield/numky/blob/master/src/main/kotlin/org/nield/numky/linear/DoubleOperators.kt + * + */ + +/* + * Functions that help create a real (Double) matrix + */ + +typealias RealMatrix = Matrix + +fun realMatrix(rowNum: Int, colNum: Int, initializer: (i: Int, j: Int) -> Double): RealMatrix = + MatrixContext.real.produce(rowNum, colNum, initializer) + +fun Sequence.toMatrix(): RealMatrix = toList().let { + MatrixContext.real.produce(it.size, it[0].size) { row, col -> it[row][col] } +} + +fun Matrix.repeatStackVertical(n: Int): RealMatrix = + VirtualMatrix(rowNum * n, colNum) { row, col -> + get(if (row == 0) 0 else row % rowNum, col) + } + +/* + * Operations for matrix and real number + */ + +operator fun Matrix.times(double: Double): RealMatrix = + MatrixContext.real.produce(rowNum, colNum) { row, col -> + this[row, col] * double + } + +operator fun Matrix.plus(double: Double): RealMatrix = + MatrixContext.real.produce(rowNum, colNum) { row, col -> + this[row, col] + double + } + +operator fun Matrix.minus(double: Double): RealMatrix = + MatrixContext.real.produce(rowNum, colNum) { row, col -> + this[row, col] - double + } + +operator fun Matrix.div(double: Double): RealMatrix = + MatrixContext.real.produce(rowNum, colNum) { row, col -> + this[row, col] / double + } + +operator fun Double.times(matrix: Matrix): RealMatrix = + MatrixContext.real.produce(matrix.rowNum, matrix.colNum) { row, col -> + this * matrix[row, col] + } + +operator fun Double.plus(matrix: Matrix): RealMatrix = + MatrixContext.real.produce(matrix.rowNum, matrix.colNum) { row, col -> + this + matrix[row, col] + } + +operator fun Double.minus(matrix: Matrix): RealMatrix = + MatrixContext.real.produce(matrix.rowNum, matrix.colNum) { row, col -> + this - matrix[row, col] + } + +// TODO: does this operation make sense? Should it be 'this/matrix[row, col]'? +//operator fun Double.div(matrix: Matrix) = MatrixContext.real.produce(matrix.rowNum, matrix.colNum) { +// row, col -> matrix[row, col] / this +//} + +/* + * Per-element (!) square and power operations + */ + +fun Matrix.square(): RealMatrix = MatrixContext.real.produce(rowNum, colNum) { row, col -> + this[row, col].pow(2) +} + +fun Matrix.pow(n: Int): RealMatrix = MatrixContext.real.produce(rowNum, colNum) { i, j -> + this[i, j].pow(n) +} + +/* + * Operations on two matrices (per-element!) + */ + +operator fun Matrix.times(other: Matrix): RealMatrix = + MatrixContext.real.produce(rowNum, colNum) { row, col -> + this[row, col] * other[row, col] + } + +operator fun Matrix.plus(other: Matrix): RealMatrix = + MatrixContext.real.add(this, other) + +operator fun Matrix.minus(other: Matrix): RealMatrix = + MatrixContext.real.produce(rowNum, colNum) { row, col -> + this[row, col] - other[row, col] + } + +/* + * Operations on columns + */ + +inline fun Matrix.appendColumn(crossinline mapper: (Buffer) -> Double) = + MatrixContext.real.produce(rowNum, colNum + 1) { row, col -> + if (col < colNum) + this[row, col] + else + mapper(rows[row]) + } + +fun Matrix.extractColumns(columnRange: IntRange): RealMatrix = + MatrixContext.real.produce(rowNum, columnRange.count()) { row, col -> + this[row, columnRange.first + col] + } + +fun Matrix.extractColumn(columnIndex: Int): RealMatrix = + extractColumns(columnIndex..columnIndex) + +fun Matrix.sumByColumn(): DoubleBuffer = DoubleBuffer(colNum) { j -> + val column = columns[j] + with(elementContext) { + sum(column.asIterable()) + } +} + +fun Matrix.minByColumn(): DoubleBuffer = DoubleBuffer(colNum) { j -> + columns[j].asIterable().min() ?: throw Exception("Cannot produce min on empty column") +} + +fun Matrix.maxByColumn(): DoubleBuffer = DoubleBuffer(colNum) { j -> + columns[j].asIterable().max() ?: throw Exception("Cannot produce min on empty column") +} + +fun Matrix.averageByColumn(): DoubleBuffer = DoubleBuffer(colNum) { j -> + columns[j].asIterable().average() +} + +/* + * Operations processing all elements + */ + +fun Matrix.sum() = elements().map { (_, value) -> value }.sum() + +fun Matrix.min() = elements().map { (_, value) -> value }.min() + +fun Matrix.max() = elements().map { (_, value) -> value }.max() + +fun Matrix.average() = elements().map { (_, value) -> value }.average() diff --git a/kmath-for-real/src/commonTest/kotlin/scientific.kmath.real/RealMatrixTest.kt b/kmath-for-real/src/commonTest/kotlin/scientific.kmath.real/RealMatrixTest.kt index 31b8b5252..8918fb300 100644 --- a/kmath-for-real/src/commonTest/kotlin/scientific.kmath.real/RealMatrixTest.kt +++ b/kmath-for-real/src/commonTest/kotlin/scientific.kmath.real/RealMatrixTest.kt @@ -1,11 +1,12 @@ package scientific.kmath.real -import scientifik.kmath.real.* import scientifik.kmath.linear.VirtualMatrix import scientifik.kmath.linear.build +import scientifik.kmath.real.* import scientifik.kmath.structures.Matrix import kotlin.test.Test import kotlin.test.assertEquals +import kotlin.test.assertTrue class RealMatrixTest { @Test @@ -19,72 +20,72 @@ class RealMatrixTest { fun testSequenceToMatrix() { val m = Sequence { listOf( - DoubleArray(10) { 10.0 }, - DoubleArray(10) { 20.0 }, - DoubleArray(10) { 30.0 }).iterator() + DoubleArray(10) { 10.0 }, + DoubleArray(10) { 20.0 }, + DoubleArray(10) { 30.0 }).iterator() }.toMatrix() assertEquals(m.sum(), 20.0 * 30) } @Test fun testRepeatStackVertical() { - val matrix1 = Matrix.build(2, 3)( - 1.0, 0.0, 0.0, - 0.0, 1.0, 2.0 + val matrix1 = Matrix.build(2, 3)( + 1.0, 0.0, 0.0, + 0.0, 1.0, 2.0 ) - val matrix2 = Matrix.build(6, 3)( - 1.0, 0.0, 0.0, - 0.0, 1.0, 2.0, - 1.0, 0.0, 0.0, - 0.0, 1.0, 2.0, - 1.0, 0.0, 0.0, - 0.0, 1.0, 2.0 + val matrix2 = Matrix.build(6, 3)( + 1.0, 0.0, 0.0, + 0.0, 1.0, 2.0, + 1.0, 0.0, 0.0, + 0.0, 1.0, 2.0, + 1.0, 0.0, 0.0, + 0.0, 1.0, 2.0 ) assertEquals(VirtualMatrix.wrap(matrix2), matrix1.repeatStackVertical(3)) } @Test fun testMatrixAndDouble() { - val matrix1 = Matrix.build(2, 3)( - 1.0, 0.0, 3.0, - 4.0, 6.0, 2.0 + val matrix1 = Matrix.build(2, 3)( + 1.0, 0.0, 3.0, + 4.0, 6.0, 2.0 ) val matrix2 = (matrix1 * 2.5 + 1.0 - 2.0) / 2.0 - val expectedResult = Matrix.build(2, 3)( - 0.75, -0.5, 3.25, - 4.5, 7.0, 2.0 + val expectedResult = Matrix.build(2, 3)( + 0.75, -0.5, 3.25, + 4.5, 7.0, 2.0 ) assertEquals(matrix2, expectedResult) } @Test fun testDoubleAndMatrix() { - val matrix1 = Matrix.build(2, 3)( - 1.0, 0.0, 3.0, - 4.0, 6.0, 2.0 + val matrix1 = Matrix.build(2, 3)( + 1.0, 0.0, 3.0, + 4.0, 6.0, 2.0 ) val matrix2 = 20.0 - (10.0 + (5.0 * matrix1)) //val matrix2 = 10.0 + (5.0 * matrix1) - val expectedResult = Matrix.build(2, 3)( - 5.0, 10.0, -5.0, - -10.0, -20.0, 0.0 + val expectedResult = Matrix.build(2, 3)( + 5.0, 10.0, -5.0, + -10.0, -20.0, 0.0 ) assertEquals(matrix2, expectedResult) } @Test fun testSquareAndPower() { - val matrix1 = Matrix.build(2, 3)( - -1.0, 0.0, 3.0, - 4.0, -6.0, -2.0 + val matrix1 = Matrix.build(2, 3)( + -1.0, 0.0, 3.0, + 4.0, -6.0, -2.0 ) - val matrix2 = Matrix.build(2, 3)( - 1.0, 0.0, 9.0, - 16.0, 36.0, 4.0 + val matrix2 = Matrix.build(2, 3)( + 1.0, 0.0, 9.0, + 16.0, 36.0, 4.0 ) - val matrix3 = Matrix.build(2, 3)( - -1.0, 0.0, 27.0, - 64.0, -216.0, -8.0 + val matrix3 = Matrix.build(2, 3)( + -1.0, 0.0, 27.0, + 64.0, -216.0, -8.0 ) assertEquals(matrix1.square(), matrix2) assertEquals(matrix1.pow(3), matrix3) @@ -92,51 +93,61 @@ class RealMatrixTest { @Test fun testTwoMatrixOperations() { - val matrix1 = Matrix.build(2, 3)( - -1.0, 0.0, 3.0, - 4.0, -6.0, 7.0 + val matrix1 = Matrix.build(2, 3)( + -1.0, 0.0, 3.0, + 4.0, -6.0, 7.0 ) - val matrix2 = Matrix.build(2, 3)( - 1.0, 0.0, 3.0, - 4.0, 6.0, -2.0 + val matrix2 = Matrix.build(2, 3)( + 1.0, 0.0, 3.0, + 4.0, 6.0, -2.0 ) val result = matrix1 * matrix2 + matrix1 - matrix2 - val expectedResult = Matrix.build(2, 3)( - -3.0, 0.0, 9.0, - 16.0, -48.0, -5.0 + val expectedResult = Matrix.build(2, 3)( + -3.0, 0.0, 9.0, + 16.0, -48.0, -5.0 ) assertEquals(result, expectedResult) } @Test fun testColumnOperations() { - val matrix1 = Matrix.build(2, 4)( - -1.0, 0.0, 3.0, 15.0, - 4.0, -6.0, 7.0, -11.0 + val matrix1 = Matrix.build(2, 4)( + -1.0, 0.0, 3.0, 15.0, + 4.0, -6.0, 7.0, -11.0 ) - val matrix2 = Matrix.build(2, 5)( - -1.0, 0.0, 3.0, 15.0, -1.0, - 4.0, -6.0, 7.0, -11.0, 4.0 + val matrix2 = Matrix.build(2, 5)( + -1.0, 0.0, 3.0, 15.0, -1.0, + 4.0, -6.0, 7.0, -11.0, 4.0 ) - val col1 = Matrix.build(2, 1)(0.0, -6.0) - val cols1to2 = Matrix.build(2, 2)( - 0.0, 3.0, - -6.0, 7.0 + val col1 = Matrix.build(2, 1)(0.0, -6.0) + val cols1to2 = Matrix.build(2, 2)( + 0.0, 3.0, + -6.0, 7.0 ) + assertEquals(matrix1.appendColumn { it[0] }, matrix2) assertEquals(matrix1.extractColumn(1), col1) assertEquals(matrix1.extractColumns(1..2), cols1to2) - assertEquals(matrix1.sumByColumn(), Matrix.build(4, 1)(3.0, -6.0, 10.0, 4.0)) - assertEquals(matrix1.minByColumn(), Matrix.build(4, 1)(-1.0, -6.0, 3.0, -11.0)) - assertEquals(matrix1.maxByColumn(), Matrix.build(4, 1)(4.0, 0.0, 7.0, 15.0)) - assertEquals(matrix1.averageByColumn(), Matrix.build(4, 1)(1.5, -3.0, 5.0, 2.0)) + //equals should never be called on buffers + assertTrue { + matrix1.sumByColumn().contentEquals(3.0, -6.0, 10.0, 4.0) + } //assertEquals(matrix1.sumByColumn(), DoubleBuffer(3.0, -6.0, 10.0, 4.0)) + assertTrue { + matrix1.minByColumn().contentEquals(-1.0, -6.0, 3.0, -11.0) + } //assertEquals(matrix1.minByColumn(), DoubleBuffer(-1.0, -6.0, 3.0, -11.0)) + assertTrue { + matrix1.maxByColumn().contentEquals(4.0, 0.0, 7.0, 15.0) + } //assertEquals(matrix1.maxByColumn(), DoubleBuffer(4.0, 0.0, 7.0, 15.0)) + assertTrue { + matrix1.averageByColumn().contentEquals(1.5, -3.0, 5.0, 2.0) + } //assertEquals(matrix1.averageByColumn(), DoubleBuffer(1.5, -3.0, 5.0, 2.0)) } @Test fun testAllElementOperations() { - val matrix1 = Matrix.build(2, 4)( - -1.0, 0.0, 3.0, 15.0, - 4.0, -6.0, 7.0, -11.0 + val matrix1 = Matrix.build(2, 4)( + -1.0, 0.0, 3.0, 15.0, + 4.0, -6.0, 7.0, -11.0 ) assertEquals(matrix1.sum(), 11.0) assertEquals(matrix1.min(), -11.0) diff --git a/kmath-for-real/src/commonTest/kotlin/scientifik/kmath/linear/VectorTest.kt b/kmath-for-real/src/commonTest/kotlin/scientifik/kmath/linear/VectorTest.kt index 0e73ee4a5..28e62b066 100644 --- a/kmath-for-real/src/commonTest/kotlin/scientifik/kmath/linear/VectorTest.kt +++ b/kmath-for-real/src/commonTest/kotlin/scientifik/kmath/linear/VectorTest.kt @@ -1,5 +1,6 @@ package scientifik.kmath.linear +import scientifik.kmath.real.RealVector import kotlin.test.Test import kotlin.test.assertEquals diff --git a/kmath-histograms/src/commonMain/kotlin/scientifik/kmath/histogram/RealHistogram.kt b/kmath-histograms/src/commonMain/kotlin/scientifik/kmath/histogram/RealHistogram.kt index 85f078fda..4438f5d60 100644 --- a/kmath-histograms/src/commonMain/kotlin/scientifik/kmath/histogram/RealHistogram.kt +++ b/kmath-histograms/src/commonMain/kotlin/scientifik/kmath/histogram/RealHistogram.kt @@ -1,7 +1,7 @@ package scientifik.kmath.histogram import scientifik.kmath.linear.Point -import scientifik.kmath.linear.asVector +import scientifik.kmath.real.asVector import scientifik.kmath.operations.SpaceOperations import scientifik.kmath.structures.* import kotlin.math.floor diff --git a/kmath-histograms/src/commonTest/kotlin/scietifik/kmath/histogram/MultivariateHistogramTest.kt b/kmath-histograms/src/commonTest/kotlin/scietifik/kmath/histogram/MultivariateHistogramTest.kt index 4dc4dfc74..5edecb5a5 100644 --- a/kmath-histograms/src/commonTest/kotlin/scietifik/kmath/histogram/MultivariateHistogramTest.kt +++ b/kmath-histograms/src/commonTest/kotlin/scietifik/kmath/histogram/MultivariateHistogramTest.kt @@ -3,7 +3,7 @@ package scietifik.kmath.histogram import scientifik.kmath.histogram.RealHistogram import scientifik.kmath.histogram.fill import scientifik.kmath.histogram.put -import scientifik.kmath.linear.RealVector +import scientifik.kmath.real.RealVector import kotlin.random.Random import kotlin.test.Test import kotlin.test.assertEquals diff --git a/kmath-histograms/src/jvmMain/kotlin/scientifik/kmath/histogram/UnivariateHistogram.kt b/kmath-histograms/src/jvmMain/kotlin/scientifik/kmath/histogram/UnivariateHistogram.kt index 90b9aff5e..dcc5ac0eb 100644 --- a/kmath-histograms/src/jvmMain/kotlin/scientifik/kmath/histogram/UnivariateHistogram.kt +++ b/kmath-histograms/src/jvmMain/kotlin/scientifik/kmath/histogram/UnivariateHistogram.kt @@ -1,7 +1,7 @@ package scientifik.kmath.histogram -import scientifik.kmath.linear.RealVector -import scientifik.kmath.linear.asVector +import scientifik.kmath.real.RealVector +import scientifik.kmath.real.asVector import scientifik.kmath.structures.Buffer import java.util.* import kotlin.math.floor diff --git a/kmath-koma/src/commonMain/kotlin/scientifik.kmath.linear/KomaMatrix.kt b/kmath-koma/src/commonMain/kotlin/scientifik.kmath.linear/KomaMatrix.kt index 74681ac48..10deabd73 100644 --- a/kmath-koma/src/commonMain/kotlin/scientifik.kmath.linear/KomaMatrix.kt +++ b/kmath-koma/src/commonMain/kotlin/scientifik.kmath.linear/KomaMatrix.kt @@ -4,6 +4,7 @@ import koma.extensions.fill import koma.matrix.MatrixFactory import scientifik.kmath.operations.Space import scientifik.kmath.structures.Matrix +import scientifik.kmath.structures.NDStructure class KomaMatrixContext( private val factory: MatrixFactory>, @@ -85,6 +86,18 @@ class KomaMatrix(val origin: koma.matrix.Matrix, features: Set ?: return false) + } + + override fun hashCode(): Int { + var result = origin.hashCode() + result = 31 * result + features.hashCode() + return result + } + + } class KomaVector internal constructor(val origin: koma.matrix.Matrix) : Point { diff --git a/kmath-memory/src/commonMain/kotlin/scientifik/memory/Memory.kt b/kmath-memory/src/commonMain/kotlin/scientifik/memory/Memory.kt index f9f938dcc..7c6886202 100644 --- a/kmath-memory/src/commonMain/kotlin/scientifik/memory/Memory.kt +++ b/kmath-memory/src/commonMain/kotlin/scientifik/memory/Memory.kt @@ -69,3 +69,9 @@ inline fun Memory.write(block: MemoryWriter.() -> Unit) { * Allocate the most effective platform-specific memory */ expect fun Memory.Companion.allocate(length: Int): Memory + +/** + * Wrap a [Memory] around existing [ByteArray]. This operation is unsafe since the array is not copied + * and could be mutated independently from the resulting [Memory] + */ +expect fun Memory.Companion.wrap(array: ByteArray): Memory diff --git a/kmath-memory/src/jsMain/kotlin/scientifik/memory/DataViewMemory.kt b/kmath-memory/src/jsMain/kotlin/scientifik/memory/DataViewMemory.kt index 843464ab9..59945efb9 100644 --- a/kmath-memory/src/jsMain/kotlin/scientifik/memory/DataViewMemory.kt +++ b/kmath-memory/src/jsMain/kotlin/scientifik/memory/DataViewMemory.kt @@ -2,14 +2,7 @@ package scientifik.memory import org.khronos.webgl.ArrayBuffer import org.khronos.webgl.DataView - -/** - * Allocate the most effective platform-specific memory - */ -actual fun Memory.Companion.allocate(length: Int): Memory { - val buffer = ArrayBuffer(length) - return DataViewMemory(DataView(buffer, 0, length)) -} +import org.khronos.webgl.Int8Array class DataViewMemory(val view: DataView) : Memory { @@ -88,4 +81,17 @@ class DataViewMemory(val view: DataView) : Memory { override fun writer(): MemoryWriter = writer +} + +/** + * Allocate the most effective platform-specific memory + */ +actual fun Memory.Companion.allocate(length: Int): Memory { + val buffer = ArrayBuffer(length) + return DataViewMemory(DataView(buffer, 0, length)) +} + +actual fun Memory.Companion.wrap(array: ByteArray): Memory { + @Suppress("CAST_NEVER_SUCCEEDS") val int8Array = array as Int8Array + return DataViewMemory(DataView(int8Array.buffer, int8Array.byteOffset, int8Array.length)) } \ No newline at end of file diff --git a/kmath-memory/src/jvmMain/kotlin/scientifik/memory/ByteBufferMemory.kt b/kmath-memory/src/jvmMain/kotlin/scientifik/memory/ByteBufferMemory.kt index df2e9847a..9ec2b3a09 100644 --- a/kmath-memory/src/jvmMain/kotlin/scientifik/memory/ByteBufferMemory.kt +++ b/kmath-memory/src/jvmMain/kotlin/scientifik/memory/ByteBufferMemory.kt @@ -7,14 +7,6 @@ import java.nio.file.Path import java.nio.file.StandardOpenOption -/** - * Allocate the most effective platform-specific memory - */ -actual fun Memory.Companion.allocate(length: Int): Memory { - val buffer = ByteBuffer.allocate(length) - return ByteBufferMemory(buffer) -} - private class ByteBufferMemory( val buffer: ByteBuffer, val startOffset: Int = 0, @@ -96,6 +88,22 @@ private class ByteBufferMemory( override fun writer(): MemoryWriter = writer } +/** + * Allocate the most effective platform-specific memory + */ +actual fun Memory.Companion.allocate(length: Int): Memory { + val buffer = ByteBuffer.allocate(length) + return ByteBufferMemory(buffer) +} + +actual fun Memory.Companion.wrap(array: ByteArray): Memory { + val buffer = ByteBuffer.wrap(array) + return ByteBufferMemory(buffer) +} + +fun ByteBuffer.asMemory(startOffset: Int = 0, size: Int = limit()): Memory = + ByteBufferMemory(this, startOffset, size) + /** * Use direct memory-mapped buffer from file to read something and close it afterwards. */ diff --git a/kmath-prob/build.gradle.kts b/kmath-prob/build.gradle.kts index 59b25d340..a69d61b73 100644 --- a/kmath-prob/build.gradle.kts +++ b/kmath-prob/build.gradle.kts @@ -8,4 +8,10 @@ kotlin.sourceSets { api(project(":kmath-coroutines")) } } + jvmMain{ + dependencies{ + api("org.apache.commons:commons-rng-sampling:1.3") + api("org.apache.commons:commons-rng-simple:1.3") + } + } } \ No newline at end of file diff --git a/kmath-prob/src/commonMain/kotlin/scientifik/kmath/prob/RandomGenerator.kt b/kmath-prob/src/commonMain/kotlin/scientifik/kmath/prob/RandomGenerator.kt index 13983c6b2..2a225fe47 100644 --- a/kmath-prob/src/commonMain/kotlin/scientifik/kmath/prob/RandomGenerator.kt +++ b/kmath-prob/src/commonMain/kotlin/scientifik/kmath/prob/RandomGenerator.kt @@ -6,10 +6,16 @@ import kotlin.random.Random * A basic generator */ interface RandomGenerator { + fun nextBoolean(): Boolean + fun nextDouble(): Double fun nextInt(): Int + fun nextInt(until: Int): Int fun nextLong(): Long - fun nextBlock(size: Int): ByteArray + fun nextLong(until: Long): Long + + fun fillBytes(array: ByteArray, fromIndex: Int = 0, toIndex: Int = array.size) + fun nextBytes(size: Int): ByteArray = ByteArray(size).also { fillBytes(it) } /** * Create a new generator which is independent from current generator (operations on new generator do not affect this one @@ -21,21 +27,29 @@ interface RandomGenerator { fun fork(): RandomGenerator companion object { - val default by lazy { DefaultGenerator(Random.nextLong()) } + val default by lazy { DefaultGenerator() } + + fun default(seed: Long) = DefaultGenerator(Random(seed)) } } -class DefaultGenerator(seed: Long?) : RandomGenerator { - private val random = seed?.let { Random(it) } ?: Random +inline class DefaultGenerator(val random: Random = Random) : RandomGenerator { + override fun nextBoolean(): Boolean = random.nextBoolean() override fun nextDouble(): Double = random.nextDouble() override fun nextInt(): Int = random.nextInt() + override fun nextInt(until: Int): Int = random.nextInt(until) override fun nextLong(): Long = random.nextLong() - override fun nextBlock(size: Int): ByteArray = random.nextBytes(size) + override fun nextLong(until: Long): Long = random.nextLong(until) - override fun fork(): RandomGenerator = DefaultGenerator(nextLong()) + override fun fillBytes(array: ByteArray, fromIndex: Int, toIndex: Int) { + random.nextBytes(array, fromIndex, toIndex) + } + override fun nextBytes(size: Int): ByteArray = random.nextBytes(size) + + override fun fork(): RandomGenerator = RandomGenerator.default(random.nextLong()) } \ No newline at end of file diff --git a/kmath-prob/src/commonMain/kotlin/scientifik/kmath/prob/UniformDistribution.kt b/kmath-prob/src/commonMain/kotlin/scientifik/kmath/prob/UniformDistribution.kt index ff6e2a551..9d96bff59 100644 --- a/kmath-prob/src/commonMain/kotlin/scientifik/kmath/prob/UniformDistribution.kt +++ b/kmath-prob/src/commonMain/kotlin/scientifik/kmath/prob/UniformDistribution.kt @@ -28,4 +28,7 @@ class UniformDistribution(val range: ClosedFloatingPointRange) : Univari else -> (arg - range.start) / length } } -} \ No newline at end of file +} + +fun Distribution.Companion.uniform(range: ClosedFloatingPointRange): UniformDistribution = + UniformDistribution(range) \ No newline at end of file diff --git a/kmath-prob/src/jvmMain/kotlin/scientifik/kmath/prob/RandomSourceGenerator.kt b/kmath-prob/src/jvmMain/kotlin/scientifik/kmath/prob/RandomSourceGenerator.kt new file mode 100644 index 000000000..f5a73a08b --- /dev/null +++ b/kmath-prob/src/jvmMain/kotlin/scientifik/kmath/prob/RandomSourceGenerator.kt @@ -0,0 +1,67 @@ +package scientifik.kmath.prob + +import org.apache.commons.rng.UniformRandomProvider +import org.apache.commons.rng.simple.RandomSource + +class RandomSourceGenerator(val source: RandomSource, seed: Long?) : RandomGenerator { + internal val random: UniformRandomProvider = seed?.let { + RandomSource.create(source, seed) + } ?: RandomSource.create(source) + + override fun nextBoolean(): Boolean = random.nextBoolean() + + override fun nextDouble(): Double = random.nextDouble() + + override fun nextInt(): Int = random.nextInt() + override fun nextInt(until: Int): Int = random.nextInt(until) + + override fun nextLong(): Long = random.nextLong() + override fun nextLong(until: Long): Long = random.nextLong(until) + + override fun fillBytes(array: ByteArray, fromIndex: Int, toIndex: Int) { + require(toIndex > fromIndex) + random.nextBytes(array, fromIndex, toIndex - fromIndex) + } + + override fun fork(): RandomGenerator = RandomSourceGenerator(source, nextLong()) +} + +inline class RandomGeneratorProvider(val generator: RandomGenerator) : UniformRandomProvider { + override fun nextBoolean(): Boolean = generator.nextBoolean() + + override fun nextFloat(): Float = generator.nextDouble().toFloat() + + override fun nextBytes(bytes: ByteArray) { + generator.fillBytes(bytes) + } + + override fun nextBytes(bytes: ByteArray, start: Int, len: Int) { + generator.fillBytes(bytes, start, start + len) + } + + override fun nextInt(): Int = generator.nextInt() + + override fun nextInt(n: Int): Int = generator.nextInt(n) + + override fun nextDouble(): Double = generator.nextDouble() + + override fun nextLong(): Long = generator.nextLong() + + override fun nextLong(n: Long): Long = generator.nextLong(n) +} + +/** + * Represent this [RandomGenerator] as commons-rng [UniformRandomProvider] preserving and mirroring its current state. + * Getting new value from one of those changes the state of another. + */ +fun RandomGenerator.asUniformRandomProvider(): UniformRandomProvider = if (this is RandomSourceGenerator) { + random +} else { + RandomGeneratorProvider(this) +} + +fun RandomGenerator.Companion.fromSource(source: RandomSource, seed: Long? = null): RandomSourceGenerator = + RandomSourceGenerator(source, seed) + +fun RandomGenerator.Companion.mersenneTwister(seed: Long? = null): RandomSourceGenerator = + fromSource(RandomSource.MT, seed) diff --git a/kmath-prob/src/jvmMain/kotlin/scientifik/kmath/prob/distributions.kt b/kmath-prob/src/jvmMain/kotlin/scientifik/kmath/prob/distributions.kt new file mode 100644 index 000000000..412454994 --- /dev/null +++ b/kmath-prob/src/jvmMain/kotlin/scientifik/kmath/prob/distributions.kt @@ -0,0 +1,109 @@ +package scientifik.kmath.prob + +import org.apache.commons.rng.UniformRandomProvider +import org.apache.commons.rng.sampling.distribution.* +import scientifik.kmath.chains.BlockingIntChain +import scientifik.kmath.chains.BlockingRealChain +import scientifik.kmath.chains.Chain +import java.util.* +import kotlin.math.PI +import kotlin.math.exp +import kotlin.math.pow +import kotlin.math.sqrt + +abstract class ContinuousSamplerDistribution : Distribution { + + private inner class ContinuousSamplerChain(val generator: RandomGenerator) : BlockingRealChain() { + private val sampler = buildCMSampler(generator) + + override fun nextDouble(): Double = sampler.sample() + + override fun fork(): Chain = ContinuousSamplerChain(generator.fork()) + } + + protected abstract fun buildCMSampler(generator: RandomGenerator): ContinuousSampler + + override fun sample(generator: RandomGenerator): BlockingRealChain = ContinuousSamplerChain(generator) +} + +abstract class DiscreteSamplerDistribution : Distribution { + + private inner class ContinuousSamplerChain(val generator: RandomGenerator) : BlockingIntChain() { + private val sampler = buildSampler(generator) + + override fun nextInt(): Int = sampler.sample() + + override fun fork(): Chain = ContinuousSamplerChain(generator.fork()) + } + + protected abstract fun buildSampler(generator: RandomGenerator): DiscreteSampler + + override fun sample(generator: RandomGenerator): BlockingIntChain = ContinuousSamplerChain(generator) +} + +enum class NormalSamplerMethod { + BoxMuller, + Marsaglia, + Ziggurat +} + +private fun normalSampler(method: NormalSamplerMethod, provider: UniformRandomProvider): NormalizedGaussianSampler = + when (method) { + NormalSamplerMethod.BoxMuller -> BoxMullerNormalizedGaussianSampler(provider) + NormalSamplerMethod.Marsaglia -> MarsagliaNormalizedGaussianSampler(provider) + NormalSamplerMethod.Ziggurat -> ZigguratNormalizedGaussianSampler(provider) + } + +fun Distribution.Companion.normal( + method: NormalSamplerMethod = NormalSamplerMethod.Ziggurat +): Distribution = object : ContinuousSamplerDistribution() { + override fun buildCMSampler(generator: RandomGenerator): ContinuousSampler { + val provider: UniformRandomProvider = generator.asUniformRandomProvider() + return normalSampler(method, provider) + } + + override fun probability(arg: Double): Double { + return exp(-arg.pow(2) / 2) / sqrt(PI * 2) + } +} + +fun Distribution.Companion.normal( + mean: Double, + sigma: Double, + method: NormalSamplerMethod = NormalSamplerMethod.Ziggurat +): ContinuousSamplerDistribution = object : ContinuousSamplerDistribution() { + private val sigma2 = sigma.pow(2) + private val norm = sigma * sqrt(PI * 2) + + override fun buildCMSampler(generator: RandomGenerator): ContinuousSampler { + val provider: UniformRandomProvider = generator.asUniformRandomProvider() + val normalizedSampler = normalSampler(method, provider) + return GaussianSampler(normalizedSampler, mean, sigma) + } + + override fun probability(arg: Double): Double { + return exp(-(arg - mean).pow(2) / 2 / sigma2) / norm + } +} + +fun Distribution.Companion.poisson( + lambda: Double +): DiscreteSamplerDistribution = object : DiscreteSamplerDistribution() { + + override fun buildSampler(generator: RandomGenerator): DiscreteSampler { + return PoissonSampler.of(generator.asUniformRandomProvider(), lambda) + } + + private val computedProb: HashMap = hashMapOf(0 to exp(-lambda)) + + override fun probability(arg: Int): Double { + require(arg >= 0) { "The argument must be >= 0" } + return if (arg > 40) { + exp(-(arg - lambda).pow(2) / 2 / lambda) / sqrt(2 * PI * lambda) + } else { + computedProb.getOrPut(arg) { + probability(arg - 1) * lambda / arg + } + } + } +} diff --git a/kmath-prob/src/jvmTest/kotlin/scientifik/kmath/prob/CommonsDistributionsTest.kt b/kmath-prob/src/jvmTest/kotlin/scientifik/kmath/prob/CommonsDistributionsTest.kt new file mode 100644 index 000000000..7638c695e --- /dev/null +++ b/kmath-prob/src/jvmTest/kotlin/scientifik/kmath/prob/CommonsDistributionsTest.kt @@ -0,0 +1,28 @@ +package scientifik.kmath.prob + +import kotlinx.coroutines.flow.take +import kotlinx.coroutines.flow.toList +import kotlinx.coroutines.runBlocking +import org.junit.jupiter.api.Assertions +import org.junit.jupiter.api.Test + +class CommonsDistributionsTest { + @Test + fun testNormalDistributionSuspend() { + val distribution = Distribution.normal(7.0, 2.0) + val generator = RandomGenerator.default(1) + val sample = runBlocking { + distribution.sample(generator).take(1000).toList() + } + Assertions.assertEquals(7.0, sample.average(), 0.1) + } + + @Test + fun testNormalDistributionBlocking() { + val distribution = Distribution.normal(7.0, 2.0) + val generator = RandomGenerator.default(1) + val sample = distribution.sample(generator).nextBlock(1000) + Assertions.assertEquals(7.0, sample.average(), 0.1) + } + +} \ No newline at end of file diff --git a/kmath-prob/src/jvmTest/kotlin/scientifik/kmath/prob/StatisticTest.kt b/kmath-prob/src/jvmTest/kotlin/scientifik/kmath/prob/StatisticTest.kt index af069810f..2613f71d5 100644 --- a/kmath-prob/src/jvmTest/kotlin/scientifik/kmath/prob/StatisticTest.kt +++ b/kmath-prob/src/jvmTest/kotlin/scientifik/kmath/prob/StatisticTest.kt @@ -9,7 +9,7 @@ import kotlin.test.Test class StatisticTest { //create a random number generator. - val generator = DefaultGenerator(1) + val generator = RandomGenerator.default(1) //Create a stateless chain from generator. val data = generator.chain { nextDouble() } //Convert a chaint to Flow and break it into chunks. diff --git a/settings.gradle.kts b/settings.gradle.kts index a08d5f7ee..57173250b 100644 --- a/settings.gradle.kts +++ b/settings.gradle.kts @@ -1,10 +1,12 @@ pluginManagement { + val toolsVersion = "0.5.0" + plugins { - id("scientifik.mpp") version "0.4.1" - id("scientifik.jvm") version "0.4.1" - id("scientifik.atomic") version "0.4.1" - id("scientifik.publish") version "0.4.1" + id("scientifik.mpp") version toolsVersion + id("scientifik.jvm") version toolsVersion + id("scientifik.atomic") version toolsVersion + id("scientifik.publish") version toolsVersion } repositories { @@ -20,7 +22,7 @@ pluginManagement { resolutionStrategy { eachPlugin { when (requested.id.id) { - "scientifik.mpp", "scientifik.jvm", "scientifik.publish" -> useModule("scientifik:gradle-tools:${requested.version}") + "scientifik.mpp", "scientifik.jvm", "scientifik.publish" -> useModule("scientifik:gradle-tools:$toolsVersion") } } }