From 048a1ceaed926614d5307ae3ca399cba1fcfa951 Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Fri, 22 May 2020 21:24:20 +0300 Subject: [PATCH 1/3] Moved probability distributions to commons-rng and to prob module --- build.gradle.kts | 2 +- .../commons/prob/CMRandomGeneratorWrapper.kt | 32 ------ .../kmath/commons/prob/CommonsDistribution.kt | 82 -------------- .../kotlin/scientifik/kmath/chains/Chain.kt | 8 +- .../kotlin/scientifik/memory/Memory.kt | 6 + .../scientifik/memory/DataViewMemory.kt | 22 ++-- .../scientifik/memory/ByteBufferMemory.kt | 24 ++-- kmath-prob/build.gradle.kts | 6 + .../scientifik/kmath/prob/RandomGenerator.kt | 26 ++++- .../kmath/prob/UniformDistribution.kt | 5 +- .../kmath/prob/RandomSourceGenerator.kt | 61 ++++++++++ .../scientifik/kmath/prob/distributions.kt | 107 ++++++++++++++++++ .../kmath/prob/CommonsDistributionsTest.kt | 19 ++++ .../scientifik/kmath/prob/StatisticTest.kt | 2 +- 14 files changed, 261 insertions(+), 141 deletions(-) delete mode 100644 kmath-commons/src/main/kotlin/scientifik/kmath/commons/prob/CMRandomGeneratorWrapper.kt delete mode 100644 kmath-commons/src/main/kotlin/scientifik/kmath/commons/prob/CommonsDistribution.kt create mode 100644 kmath-prob/src/jvmMain/kotlin/scientifik/kmath/prob/RandomSourceGenerator.kt create mode 100644 kmath-prob/src/jvmMain/kotlin/scientifik/kmath/prob/distributions.kt create mode 100644 kmath-prob/src/jvmTest/kotlin/scientifik/kmath/prob/CommonsDistributionsTest.kt diff --git a/build.gradle.kts b/build.gradle.kts index 4a7e99d90..6d102a77a 100644 --- a/build.gradle.kts +++ b/build.gradle.kts @@ -2,7 +2,7 @@ plugins { id("scientifik.publish") apply false } -val kmathVersion by extra("0.1.4-dev-6") +val kmathVersion by extra("0.1.4-dev-7") val bintrayRepo by extra("scientifik") val githubProject by extra("kmath") 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-coroutines/src/commonMain/kotlin/scientifik/kmath/chains/Chain.kt b/kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/chains/Chain.kt index fd68636f9..85db55480 100644 --- a/kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/chains/Chain.kt +++ b/kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/chains/Chain.kt @@ -38,9 +38,13 @@ 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 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..2a97fa4be --- /dev/null +++ b/kmath-prob/src/jvmMain/kotlin/scientifik/kmath/prob/RandomSourceGenerator.kt @@ -0,0 +1,61 @@ +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, null) + } ?: 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) +} 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..59588d0ef --- /dev/null +++ b/kmath-prob/src/jvmMain/kotlin/scientifik/kmath/prob/distributions.kt @@ -0,0 +1,107 @@ +package scientifik.kmath.prob + +import org.apache.commons.rng.UniformRandomProvider +import org.apache.commons.rng.sampling.distribution.* +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) : Chain { + private val sampler = buildSampler(generator) + + override suspend fun next(): Double = sampler.sample() + + override fun fork(): Chain = ContinuousSamplerChain(generator.fork()) + } + + protected abstract fun buildSampler(generator: RandomGenerator): ContinuousSampler + + override fun sample(generator: RandomGenerator): Chain = ContinuousSamplerChain(generator) +} + +abstract class DiscreteSamplerDistribution : Distribution { + + private inner class ContinuousSamplerChain(val generator: RandomGenerator) : Chain { + private val sampler = buildSampler(generator) + + override suspend fun next(): Int = sampler.sample() + + override fun fork(): Chain = ContinuousSamplerChain(generator.fork()) + } + + protected abstract fun buildSampler(generator: RandomGenerator): DiscreteSampler + + override fun sample(generator: RandomGenerator): Chain = 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 buildSampler(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 +): Distribution = object : ContinuousSamplerDistribution() { + private val sigma2 = sigma.pow(2) + private val norm = sigma * sqrt(PI * 2) + + override fun buildSampler(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 +): Distribution = 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..d6fd629b7 --- /dev/null +++ b/kmath-prob/src/jvmTest/kotlin/scientifik/kmath/prob/CommonsDistributionsTest.kt @@ -0,0 +1,19 @@ +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 testNormalDistribution(){ + 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) + } +} \ 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. From 3f68c0c34e548c8476e6d23da22e026d4bd206b6 Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Fri, 22 May 2020 21:28:14 +0300 Subject: [PATCH 2/3] Fix demo for distributions --- examples/build.gradle.kts | 1 + .../kotlin/scientifik/kmath/commons/prob/DistributionDemo.kt | 3 ++- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/examples/build.gradle.kts b/examples/build.gradle.kts index 8ab0dd1eb..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")) 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 From 8533a31677b2b78cfa4d22b214cd4490dcfe0204 Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Sat, 23 May 2020 10:45:02 +0300 Subject: [PATCH 3/3] Optimized blocking chains for primitive numbers generation. --- .../commons/prob/DistributionBenchmark.kt | 71 +++++++++++++++++++ .../kmath/chains/BlockingIntChain.kt | 12 ++++ .../kmath/chains/BlockingRealChain.kt | 12 ++++ .../kotlin/scientifik/kmath/chains/Chain.kt | 1 - .../scientifik/kmath/streaming/BufferFlow.kt | 34 +++++---- .../kmath/prob/RandomSourceGenerator.kt | 8 ++- .../scientifik/kmath/prob/distributions.kt | 26 +++---- .../kmath/prob/CommonsDistributionsTest.kt | 13 +++- 8 files changed, 149 insertions(+), 28 deletions(-) create mode 100644 examples/src/main/kotlin/scientifik/kmath/commons/prob/DistributionBenchmark.kt create mode 100644 kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/chains/BlockingIntChain.kt create mode 100644 kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/chains/BlockingRealChain.kt 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/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 85db55480..5635499e5 100644 --- a/kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/chains/Chain.kt +++ b/kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/chains/Chain.kt @@ -48,7 +48,6 @@ interface Chain: Flow { } companion object - } 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-prob/src/jvmMain/kotlin/scientifik/kmath/prob/RandomSourceGenerator.kt b/kmath-prob/src/jvmMain/kotlin/scientifik/kmath/prob/RandomSourceGenerator.kt index 2a97fa4be..f5a73a08b 100644 --- a/kmath-prob/src/jvmMain/kotlin/scientifik/kmath/prob/RandomSourceGenerator.kt +++ b/kmath-prob/src/jvmMain/kotlin/scientifik/kmath/prob/RandomSourceGenerator.kt @@ -5,7 +5,7 @@ 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, null) + RandomSource.create(source, seed) } ?: RandomSource.create(source) override fun nextBoolean(): Boolean = random.nextBoolean() @@ -59,3 +59,9 @@ fun RandomGenerator.asUniformRandomProvider(): UniformRandomProvider = if (this } 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 index 59588d0ef..412454994 100644 --- a/kmath-prob/src/jvmMain/kotlin/scientifik/kmath/prob/distributions.kt +++ b/kmath-prob/src/jvmMain/kotlin/scientifik/kmath/prob/distributions.kt @@ -2,6 +2,8 @@ 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 @@ -11,32 +13,32 @@ import kotlin.math.sqrt abstract class ContinuousSamplerDistribution : Distribution { - private inner class ContinuousSamplerChain(val generator: RandomGenerator) : Chain { - private val sampler = buildSampler(generator) + private inner class ContinuousSamplerChain(val generator: RandomGenerator) : BlockingRealChain() { + private val sampler = buildCMSampler(generator) - override suspend fun next(): Double = sampler.sample() + override fun nextDouble(): Double = sampler.sample() override fun fork(): Chain = ContinuousSamplerChain(generator.fork()) } - protected abstract fun buildSampler(generator: RandomGenerator): ContinuousSampler + protected abstract fun buildCMSampler(generator: RandomGenerator): ContinuousSampler - override fun sample(generator: RandomGenerator): Chain = ContinuousSamplerChain(generator) + override fun sample(generator: RandomGenerator): BlockingRealChain = ContinuousSamplerChain(generator) } abstract class DiscreteSamplerDistribution : Distribution { - private inner class ContinuousSamplerChain(val generator: RandomGenerator) : Chain { + private inner class ContinuousSamplerChain(val generator: RandomGenerator) : BlockingIntChain() { private val sampler = buildSampler(generator) - override suspend fun next(): Int = sampler.sample() + 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): Chain = ContinuousSamplerChain(generator) + override fun sample(generator: RandomGenerator): BlockingIntChain = ContinuousSamplerChain(generator) } enum class NormalSamplerMethod { @@ -55,7 +57,7 @@ private fun normalSampler(method: NormalSamplerMethod, provider: UniformRandomPr fun Distribution.Companion.normal( method: NormalSamplerMethod = NormalSamplerMethod.Ziggurat ): Distribution = object : ContinuousSamplerDistribution() { - override fun buildSampler(generator: RandomGenerator): ContinuousSampler { + override fun buildCMSampler(generator: RandomGenerator): ContinuousSampler { val provider: UniformRandomProvider = generator.asUniformRandomProvider() return normalSampler(method, provider) } @@ -69,11 +71,11 @@ fun Distribution.Companion.normal( mean: Double, sigma: Double, method: NormalSamplerMethod = NormalSamplerMethod.Ziggurat -): Distribution = object : ContinuousSamplerDistribution() { +): ContinuousSamplerDistribution = object : ContinuousSamplerDistribution() { private val sigma2 = sigma.pow(2) private val norm = sigma * sqrt(PI * 2) - override fun buildSampler(generator: RandomGenerator): ContinuousSampler { + override fun buildCMSampler(generator: RandomGenerator): ContinuousSampler { val provider: UniformRandomProvider = generator.asUniformRandomProvider() val normalizedSampler = normalSampler(method, provider) return GaussianSampler(normalizedSampler, mean, sigma) @@ -86,7 +88,7 @@ fun Distribution.Companion.normal( fun Distribution.Companion.poisson( lambda: Double -): Distribution = object : DiscreteSamplerDistribution() { +): DiscreteSamplerDistribution = object : DiscreteSamplerDistribution() { override fun buildSampler(generator: RandomGenerator): DiscreteSampler { return PoissonSampler.of(generator.asUniformRandomProvider(), lambda) diff --git a/kmath-prob/src/jvmTest/kotlin/scientifik/kmath/prob/CommonsDistributionsTest.kt b/kmath-prob/src/jvmTest/kotlin/scientifik/kmath/prob/CommonsDistributionsTest.kt index d6fd629b7..7638c695e 100644 --- a/kmath-prob/src/jvmTest/kotlin/scientifik/kmath/prob/CommonsDistributionsTest.kt +++ b/kmath-prob/src/jvmTest/kotlin/scientifik/kmath/prob/CommonsDistributionsTest.kt @@ -8,12 +8,21 @@ import org.junit.jupiter.api.Test class CommonsDistributionsTest { @Test - fun testNormalDistribution(){ - val distribution = Distribution.normal(7.0,2.0) + 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