diff --git a/CHANGELOG.md b/CHANGELOG.md index 8b4dddf81..4ade9cd9c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,7 @@ - Buffer factories for primitives moved to MutableBuffer.Companion - NDStructure and NDAlgebra to StructureND and AlgebraND respectively - Real -> Double +- DataSets are moved from functions to core ### Deprecated diff --git a/README.md b/README.md index cc9439d27..7080c757e 100644 --- a/README.md +++ b/README.md @@ -1,9 +1,8 @@ [![JetBrains Research](https://jb.gg/badges/research.svg)](https://confluence.jetbrains.com/display/ALL/JetBrains+on+GitHub) [![DOI](https://zenodo.org/badge/129486382.svg)](https://zenodo.org/badge/latestdoi/129486382) - ![Gradle build](https://github.com/mipt-npm/kmath/workflows/Gradle%20build/badge.svg) - -[![Maven Central](https://img.shields.io/maven-central/v/space.kscience/kmath-core.svg?label=Maven%20Central)](https://search.maven.org/search?q=g:%22space.kscience%22%20AND%20a:%22kmath-core%22) +[![Maven Central](https://img.shields.io/maven-central/v/space.kscience/kmath-core.svg?label=Maven%20Central)](https://search.maven.org/search?q=g:%22space.kscience%22) +[![Space](https://img.shields.io/maven-metadata/v?label=Space&metadataUrl=https%3A%2F%2Fmaven.pkg.jetbrains.space%2Fmipt-npm%2Fp%2Fsci%2Fmaven%2Fkscience%2Fkmath%2Fkmath-core%2Fmaven-metadata.xml)](https://maven.pkg.jetbrains.space/mipt-npm/p/sci/maven/space/kscience/) # KMath @@ -147,6 +146,12 @@ performance calculations to code generation. > > > **Maturity**: PROTOTYPE +> +> **Features:** +> - [ejml-vector](kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/EjmlVector.kt) : The Point implementation using SimpleMatrix. +> - [ejml-matrix](kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/EjmlMatrix.kt) : The Matrix implementation using SimpleMatrix. +> - [ejml-linear-space](kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/EjmlLinearSpace.kt) : The LinearSpace implementation using SimpleMatrix. +
* ### [kmath-for-real](kmath-for-real) diff --git a/build.gradle.kts b/build.gradle.kts index 5f1a8b88a..cc863a957 100644 --- a/build.gradle.kts +++ b/build.gradle.kts @@ -1,4 +1,6 @@ +import org.jetbrains.dokka.gradle.DokkaTask import ru.mipt.npm.gradle.KSciencePublishingPlugin +import java.net.URL plugins { id("ru.mipt.npm.gradle.project") @@ -10,21 +12,38 @@ allprojects { maven("https://clojars.org/repo") maven("https://dl.bintray.com/egor-bogomolov/astminer/") maven("https://dl.bintray.com/hotkeytlt/maven") - maven("https://dl.bintray.com/kotlin/kotlin-eap") - maven("https://dl.bintray.com/kotlin/kotlinx") - maven("https://dl.bintray.com/mipt-npm/dev") - maven("https://dl.bintray.com/mipt-npm/kscience") maven("https://jitpack.io") maven("http://logicrunch.research.it.uu.se/maven/") mavenCentral() } group = "space.kscience" - version = "0.3.0-dev-3" + version = "0.3.0-dev-4" } subprojects { if (name.startsWith("kmath")) apply() + + afterEvaluate { + tasks.withType { + dokkaSourceSets.all { + val readmeFile = File(this@subprojects.projectDir, "./README.md") + if (readmeFile.exists()) + includes.setFrom(includes + readmeFile.absolutePath) + + arrayOf( + "http://ejml.org/javadoc/", + "https://commons.apache.org/proper/commons-math/javadocs/api-3.6.1/", + "https://deeplearning4j.org/api/latest/" + ).map { URL("${it}package-list") to URL(it) }.forEach { (a, b) -> + externalDocumentationLink { + packageListUrl.set(a) + url.set(b) + } + } + } + } + } } readme { diff --git a/docs/templates/ARTIFACT-TEMPLATE.md b/docs/templates/ARTIFACT-TEMPLATE.md index cb741bc6f..01d9c51da 100644 --- a/docs/templates/ARTIFACT-TEMPLATE.md +++ b/docs/templates/ARTIFACT-TEMPLATE.md @@ -1,34 +1,28 @@ -> #### Artifact: -> -> This module artifact: `${group}:${name}:${version}`. -> -> Bintray release version: [ ![Download](https://api.bintray.com/packages/mipt-npm/kscience/${name}/images/download.svg) ](https://bintray.com/mipt-npm/kscience/${name}/_latestVersion) -> -> Bintray development version: [ ![Download](https://api.bintray.com/packages/mipt-npm/dev/${name}/images/download.svg) ](https://bintray.com/mipt-npm/dev/${name}/_latestVersion) -> -> **Gradle:** -> -> ```gradle -> repositories { -> maven { url 'https://repo.kotlin.link' } -> maven { url 'https://dl.bintray.com/hotkeytlt/maven' } -> maven { url "https://dl.bintray.com/kotlin/kotlin-eap" } // include for builds based on kotlin-eap -> } -> -> dependencies { -> implementation '${group}:${name}:${version}' -> } -> ``` -> **Gradle Kotlin DSL:** -> -> ```kotlin -> repositories { -> maven("https://repo.kotlin.link") -> maven("https://dl.bintray.com/kotlin/kotlin-eap") // include for builds based on kotlin-eap -> maven("https://dl.bintray.com/hotkeytlt/maven") // required for a -> } -> -> dependencies { -> implementation("${group}:${name}:${version}") -> } -> ``` \ No newline at end of file +## Artifact: + +The Maven coordinates of this project are `${group}:${name}:${version}`. + +**Gradle:** +```gradle +repositories { + maven { url 'https://repo.kotlin.link' } + maven { url 'https://dl.bintray.com/hotkeytlt/maven' } + maven { url "https://dl.bintray.com/kotlin/kotlin-eap" } // include for builds based on kotlin-eap +} + +dependencies { + implementation '${group}:${name}:${version}' +} +``` +**Gradle Kotlin DSL:** +```kotlin +repositories { + maven("https://repo.kotlin.link") + maven("https://dl.bintray.com/kotlin/kotlin-eap") // include for builds based on kotlin-eap + maven("https://dl.bintray.com/hotkeytlt/maven") // required for a +} + +dependencies { + implementation("${group}:${name}:${version}") +} +``` \ No newline at end of file diff --git a/docs/templates/README-TEMPLATE.md b/docs/templates/README-TEMPLATE.md index 4366c8fcd..3502cdccd 100644 --- a/docs/templates/README-TEMPLATE.md +++ b/docs/templates/README-TEMPLATE.md @@ -1,9 +1,8 @@ [![JetBrains Research](https://jb.gg/badges/research.svg)](https://confluence.jetbrains.com/display/ALL/JetBrains+on+GitHub) [![DOI](https://zenodo.org/badge/129486382.svg)](https://zenodo.org/badge/latestdoi/129486382) - ![Gradle build](https://github.com/mipt-npm/kmath/workflows/Gradle%20build/badge.svg) - -[![Maven Central](https://img.shields.io/maven-central/v/space.kscience/kmath-core.svg?label=Maven%20Central)](https://search.maven.org/search?q=g:%22space.kscience%22%20AND%20a:%22kmath-core%22) +[![Maven Central](https://img.shields.io/maven-central/v/space.kscience/kmath-core.svg?label=Maven%20Central)](https://search.maven.org/search?q=g:%22space.kscience%22) +[![Space](https://img.shields.io/maven-metadata/v?label=Space&metadataUrl=https%3A%2F%2Fmaven.pkg.jetbrains.space%2Fmipt-npm%2Fp%2Fsci%2Fmaven%2Fkscience%2Fkmath%2Fkmath-core%2Fmaven-metadata.xml)](https://maven.pkg.jetbrains.space/mipt-npm/p/sci/maven/space/kscience/) # KMath diff --git a/examples/build.gradle.kts b/examples/build.gradle.kts index 5dd40b609..a48b4d0d9 100644 --- a/examples/build.gradle.kts +++ b/examples/build.gradle.kts @@ -106,6 +106,7 @@ kotlin.sourceSets.all { with(languageSettings) { useExperimentalAnnotation("kotlin.contracts.ExperimentalContracts") useExperimentalAnnotation("kotlin.ExperimentalUnsignedTypes") + useExperimentalAnnotation("space.kscience.kmath.misc.UnstableKMathAPI") } } diff --git a/examples/src/benchmarks/kotlin/space/kscience/kmath/benchmarks/ExpressionsInterpretersBenchmark.kt b/examples/src/benchmarks/kotlin/space/kscience/kmath/benchmarks/ExpressionsInterpretersBenchmark.kt index 0899241f9..2438e3979 100644 --- a/examples/src/benchmarks/kotlin/space/kscience/kmath/benchmarks/ExpressionsInterpretersBenchmark.kt +++ b/examples/src/benchmarks/kotlin/space/kscience/kmath/benchmarks/ExpressionsInterpretersBenchmark.kt @@ -9,7 +9,7 @@ import space.kscience.kmath.ast.mstInField import space.kscience.kmath.expressions.Expression import space.kscience.kmath.expressions.expressionInField import space.kscience.kmath.expressions.invoke -import space.kscience.kmath.expressions.symbol +import space.kscience.kmath.misc.symbol import space.kscience.kmath.operations.DoubleField import space.kscience.kmath.operations.bindSymbol import kotlin.random.Random diff --git a/examples/src/main/kotlin/space/kscience/kmath/ast/kotlingradSupport.kt b/examples/src/main/kotlin/space/kscience/kmath/ast/kotlingradSupport.kt index 23c9d5b41..138b3e708 100644 --- a/examples/src/main/kotlin/space/kscience/kmath/ast/kotlingradSupport.kt +++ b/examples/src/main/kotlin/space/kscience/kmath/ast/kotlingradSupport.kt @@ -3,8 +3,8 @@ package space.kscience.kmath.ast import space.kscience.kmath.asm.compile import space.kscience.kmath.expressions.derivative import space.kscience.kmath.expressions.invoke -import space.kscience.kmath.expressions.symbol import space.kscience.kmath.kotlingrad.differentiable +import space.kscience.kmath.misc.symbol import space.kscience.kmath.operations.DoubleField /** diff --git a/examples/src/main/kotlin/space/kscience/kmath/commons/fit/fitWithAutoDiff.kt b/examples/src/main/kotlin/space/kscience/kmath/commons/fit/fitWithAutoDiff.kt index ef0c29a2d..04c55b34c 100644 --- a/examples/src/main/kotlin/space/kscience/kmath/commons/fit/fitWithAutoDiff.kt +++ b/examples/src/main/kotlin/space/kscience/kmath/commons/fit/fitWithAutoDiff.kt @@ -7,11 +7,14 @@ import kscience.plotly.models.ScatterMode import kscience.plotly.models.TraceValues import space.kscience.kmath.commons.optimization.chiSquared import space.kscience.kmath.commons.optimization.minimize -import space.kscience.kmath.expressions.symbol +import space.kscience.kmath.misc.symbol +import space.kscience.kmath.optimization.FunctionOptimization +import space.kscience.kmath.optimization.OptimizationResult import space.kscience.kmath.real.DoubleVector import space.kscience.kmath.real.map import space.kscience.kmath.real.step -import space.kscience.kmath.stat.* +import space.kscience.kmath.stat.RandomGenerator +import space.kscience.kmath.stat.distributions.NormalDistribution import space.kscience.kmath.structures.asIterable import space.kscience.kmath.structures.toList import kotlin.math.pow @@ -33,10 +36,9 @@ operator fun TraceValues.invoke(vector: DoubleVector) { /** * Least squares fie with auto-differentiation. Uses `kmath-commons` and `kmath-for-real` modules. */ -fun main() { - +suspend fun main() { //A generator for a normally distributed values - val generator = Distribution.normal() + val generator = NormalDistribution(2.0, 7.0) //A chain/flow of random values with the given seed val chain = generator.sample(RandomGenerator.default(112667)) @@ -49,7 +51,7 @@ fun main() { //Perform an operation on each x value (much more effective, than numpy) val y = x.map { val value = it.pow(2) + it + 1 - value + chain.nextDouble() * sqrt(value) + value + chain.next() * sqrt(value) } // this will also work, but less effective: // val y = x.pow(2)+ x + 1 + chain.nextDouble() @@ -58,7 +60,7 @@ fun main() { val yErr = y.map { sqrt(it) }//RealVector.same(x.size, sigma) // compute differentiable chi^2 sum for given model ax^2 + bx + c - val chi2 = Fitting.chiSquared(x, y, yErr) { x1 -> + val chi2 = FunctionOptimization.chiSquared(x, y, yErr) { x1 -> //bind variables to autodiff context val a = bind(a) val b = bind(b) @@ -99,4 +101,4 @@ fun main() { } page.makeFile() -} \ No newline at end of file +} diff --git a/examples/src/main/kotlin/space/kscience/kmath/stat/DistributionBenchmark.kt b/examples/src/main/kotlin/space/kscience/kmath/stat/DistributionBenchmark.kt index 1761ed1b5..bfd138502 100644 --- a/examples/src/main/kotlin/space/kscience/kmath/stat/DistributionBenchmark.kt +++ b/examples/src/main/kotlin/space/kscience/kmath/stat/DistributionBenchmark.kt @@ -1,23 +1,24 @@ -package kscience.kmath.commons.prob +package space.kscience.kmath.stat import kotlinx.coroutines.Dispatchers import kotlinx.coroutines.async import kotlinx.coroutines.runBlocking -import org.apache.commons.rng.sampling.distribution.ZigguratNormalizedGaussianSampler +import space.kscience.kmath.stat.samplers.GaussianSampler import org.apache.commons.rng.simple.RandomSource -import space.kscience.kmath.stat.* import java.time.Duration import java.time.Instant +import org.apache.commons.rng.sampling.distribution.GaussianSampler as CMGaussianSampler +import org.apache.commons.rng.sampling.distribution.ZigguratNormalizedGaussianSampler as CMZigguratNormalizedGaussianSampler -private fun runChain(): Duration { +private suspend fun runKMathChained(): Duration { val generator = RandomGenerator.fromSource(RandomSource.MT, 123L) - val normal = Distribution.normal(NormalSamplerMethod.Ziggurat) - val chain = normal.sample(generator) + val normal = GaussianSampler.of(7.0, 2.0) + val chain = normal.sample(generator).blocking() val startTime = Instant.now() var sum = 0.0 repeat(10000001) { counter -> - sum += chain.nextDouble() + sum += chain.next() if (counter % 100000 == 0) { val duration = Duration.between(startTime, Instant.now()) @@ -29,9 +30,15 @@ private fun runChain(): Duration { return Duration.between(startTime, Instant.now()) } -private fun runDirect(): Duration { - val provider = RandomSource.create(RandomSource.MT, 123L) - val sampler = ZigguratNormalizedGaussianSampler(provider) +private fun runApacheDirect(): Duration { + val rng = RandomSource.create(RandomSource.MT, 123L) + + val sampler = CMGaussianSampler.of( + CMZigguratNormalizedGaussianSampler.of(rng), + 7.0, + 2.0 + ) + val startTime = Instant.now() var sum = 0.0 @@ -51,11 +58,9 @@ private fun runDirect(): Duration { /** * 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()}") - } +fun main(): Unit = runBlocking(Dispatchers.Default) { + val chainJob = async { runKMathChained() } + val directJob = async { runApacheDirect() } + println("KMath Chained: ${chainJob.await()}") + println("Apache Direct: ${directJob.await()}") } diff --git a/examples/src/main/kotlin/space/kscience/kmath/stat/DistributionDemo.kt b/examples/src/main/kotlin/space/kscience/kmath/stat/DistributionDemo.kt index 47b8d8111..aac7d51d4 100644 --- a/examples/src/main/kotlin/space/kscience/kmath/stat/DistributionDemo.kt +++ b/examples/src/main/kotlin/space/kscience/kmath/stat/DistributionDemo.kt @@ -3,14 +3,15 @@ package space.kscience.kmath.stat import kotlinx.coroutines.runBlocking import space.kscience.kmath.chains.Chain import space.kscience.kmath.chains.collectWithState +import space.kscience.kmath.stat.distributions.NormalDistribution /** - * The state of distribution averager + * The state of distribution averager. */ private data class AveragingChainState(var num: Int = 0, var value: Double = 0.0) /** - * Averaging + * Averaging. */ private fun Chain.mean(): Chain = collectWithState(AveragingChainState(), { it.copy() }) { chain -> val next = chain.next() @@ -21,7 +22,7 @@ private fun Chain.mean(): Chain = collectWithState(AveragingChai fun main() { - val normal = Distribution.normal() + val normal = NormalDistribution(0.0, 2.0) val chain = normal.sample(RandomGenerator.default).mean() runBlocking { @@ -32,4 +33,4 @@ fun main() { } } } -} \ No newline at end of file +} diff --git a/examples/src/main/kotlin/space/kscience/kmath/structures/StructureReadBenchmark.kt b/examples/src/main/kotlin/space/kscience/kmath/structures/StructureReadBenchmark.kt index 1c8a923c7..8b0a2ed0e 100644 --- a/examples/src/main/kotlin/space/kscience/kmath/structures/StructureReadBenchmark.kt +++ b/examples/src/main/kotlin/space/kscience/kmath/structures/StructureReadBenchmark.kt @@ -34,4 +34,4 @@ fun main() { strides.indices().forEach { res = array[strides.offset(it)] } } println("Array reading finished in $time3 millis") -} \ No newline at end of file +} diff --git a/gradle.properties b/gradle.properties index 7ff50a435..50123b16c 100644 --- a/gradle.properties +++ b/gradle.properties @@ -4,5 +4,5 @@ kotlin.mpp.stability.nowarn=true kotlin.native.enableDependencyPropagation=false kotlin.parallel.tasks.in.project=true org.gradle.configureondemand=true -org.gradle.jvmargs=-XX:MaxMetaspaceSize=2G +org.gradle.jvmargs=-XX:MaxMetaspaceSize=9G org.gradle.parallel=true diff --git a/kmath-ast/README.md b/kmath-ast/README.md index ee14604d2..ff954b914 100644 --- a/kmath-ast/README.md +++ b/kmath-ast/README.md @@ -1,6 +1,6 @@ -# Abstract Syntax Tree Expression Representation and Operations (`kmath-ast`) +# Module kmath-ast -This subproject implements the following features: +Abstract syntax tree expression representation and related optimizations. - [expression-language](src/jvmMain/kotlin/space/kscience/kmath/ast/parser.kt) : Expression language and its parser - [mst](src/commonMain/kotlin/space/kscience/kmath/ast/MST.kt) : MST (Mathematical Syntax Tree) as expression language's syntax intermediate representation @@ -10,40 +10,34 @@ This subproject implements the following features: - [mst-js-codegen](src/jsMain/kotlin/space/kscience/kmath/estree/estree.kt) : Dynamic MST to JS compiler -> #### Artifact: -> -> This module artifact: `space.kscience:kmath-ast:0.3.0-dev-3`. -> -> Bintray release version: [ ![Download](https://api.bintray.com/packages/mipt-npm/kscience/kmath-ast/images/download.svg) ](https://bintray.com/mipt-npm/kscience/kmath-ast/_latestVersion) -> -> Bintray development version: [ ![Download](https://api.bintray.com/packages/mipt-npm/dev/kmath-ast/images/download.svg) ](https://bintray.com/mipt-npm/dev/kmath-ast/_latestVersion) -> -> **Gradle:** -> -> ```gradle -> repositories { -> maven { url 'https://repo.kotlin.link' } -> maven { url 'https://dl.bintray.com/hotkeytlt/maven' } -> maven { url "https://dl.bintray.com/kotlin/kotlin-eap" } // include for builds based on kotlin-eap -> } -> -> dependencies { -> implementation 'space.kscience:kmath-ast:0.3.0-dev-3' -> } -> ``` -> **Gradle Kotlin DSL:** -> -> ```kotlin -> repositories { -> maven("https://repo.kotlin.link") -> maven("https://dl.bintray.com/kotlin/kotlin-eap") // include for builds based on kotlin-eap -> maven("https://dl.bintray.com/hotkeytlt/maven") // required for a -> } -> -> dependencies { -> implementation("space.kscience:kmath-ast:0.3.0-dev-3") -> } -> ``` +## Artifact: + +The Maven coordinates of this project are `space.kscience:kmath-ast:0.3.0-dev-3`. + +**Gradle:** +```gradle +repositories { + maven { url 'https://repo.kotlin.link' } + maven { url 'https://dl.bintray.com/hotkeytlt/maven' } + maven { url "https://dl.bintray.com/kotlin/kotlin-eap" } // include for builds based on kotlin-eap +} + +dependencies { + implementation 'space.kscience:kmath-ast:0.3.0-dev-3' +} +``` +**Gradle Kotlin DSL:** +```kotlin +repositories { + maven("https://repo.kotlin.link") + maven("https://dl.bintray.com/kotlin/kotlin-eap") // include for builds based on kotlin-eap + maven("https://dl.bintray.com/hotkeytlt/maven") // required for a +} + +dependencies { + implementation("space.kscience:kmath-ast:0.3.0-dev-3") +} +``` ## Dynamic expression code generation diff --git a/kmath-ast/docs/README-TEMPLATE.md b/kmath-ast/docs/README-TEMPLATE.md index 80e48008b..db071adb4 100644 --- a/kmath-ast/docs/README-TEMPLATE.md +++ b/kmath-ast/docs/README-TEMPLATE.md @@ -1,6 +1,6 @@ -# Abstract Syntax Tree Expression Representation and Operations (`kmath-ast`) +# Module kmath-ast -This subproject implements the following features: +Abstract syntax tree expression representation and related optimizations. ${features} diff --git a/kmath-ast/src/commonMain/kotlin/space/kscience/kmath/ast/MstExpression.kt b/kmath-ast/src/commonMain/kotlin/space/kscience/kmath/ast/MstExpression.kt index 63dfb38f7..5c43df068 100644 --- a/kmath-ast/src/commonMain/kotlin/space/kscience/kmath/ast/MstExpression.kt +++ b/kmath-ast/src/commonMain/kotlin/space/kscience/kmath/ast/MstExpression.kt @@ -1,6 +1,8 @@ package space.kscience.kmath.ast import space.kscience.kmath.expressions.* +import space.kscience.kmath.misc.StringSymbol +import space.kscience.kmath.misc.Symbol import space.kscience.kmath.operations.* import kotlin.contracts.InvocationKind import kotlin.contracts.contract diff --git a/kmath-ast/src/jsMain/kotlin/space/kscience/kmath/estree/internal/ESTreeBuilder.kt b/kmath-ast/src/jsMain/kotlin/space/kscience/kmath/estree/internal/ESTreeBuilder.kt index b4de9968d..1e966d986 100644 --- a/kmath-ast/src/jsMain/kotlin/space/kscience/kmath/estree/internal/ESTreeBuilder.kt +++ b/kmath-ast/src/jsMain/kotlin/space/kscience/kmath/estree/internal/ESTreeBuilder.kt @@ -3,7 +3,7 @@ package space.kscience.kmath.estree.internal import space.kscience.kmath.estree.internal.astring.generate import space.kscience.kmath.estree.internal.estree.* import space.kscience.kmath.expressions.Expression -import space.kscience.kmath.expressions.Symbol +import space.kscience.kmath.misc.Symbol internal class ESTreeBuilder(val bodyCallback: ESTreeBuilder.() -> BaseExpression) { private class GeneratedExpression(val executable: dynamic, val constants: Array) : Expression { diff --git a/kmath-ast/src/jvmMain/kotlin/space/kscience/kmath/asm/internal/mapIntrinsics.kt b/kmath-ast/src/jvmMain/kotlin/space/kscience/kmath/asm/internal/mapIntrinsics.kt index f54bc070c..0a0c21d8a 100644 --- a/kmath-ast/src/jvmMain/kotlin/space/kscience/kmath/asm/internal/mapIntrinsics.kt +++ b/kmath-ast/src/jvmMain/kotlin/space/kscience/kmath/asm/internal/mapIntrinsics.kt @@ -2,8 +2,8 @@ package space.kscience.kmath.asm.internal -import space.kscience.kmath.expressions.StringSymbol -import space.kscience.kmath.expressions.Symbol +import space.kscience.kmath.misc.StringSymbol +import space.kscience.kmath.misc.Symbol /** * Gets value with given [key] or throws [NoSuchElementException] whenever it is not present. diff --git a/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/expressions/DerivativeStructureExpression.kt b/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/expressions/DerivativeStructureExpression.kt index b74167c3f..58e9687e5 100644 --- a/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/expressions/DerivativeStructureExpression.kt +++ b/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/expressions/DerivativeStructureExpression.kt @@ -2,6 +2,8 @@ package space.kscience.kmath.commons.expressions import org.apache.commons.math3.analysis.differentiation.DerivativeStructure import space.kscience.kmath.expressions.* +import space.kscience.kmath.misc.StringSymbol +import space.kscience.kmath.misc.Symbol import space.kscience.kmath.misc.UnstableKMathAPI import space.kscience.kmath.operations.ExtendedField import space.kscience.kmath.operations.NumbersAddOperations diff --git a/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/linear/CMMatrix.kt b/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/linear/CMMatrix.kt index 89e9649e9..80929e6b9 100644 --- a/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/linear/CMMatrix.kt +++ b/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/linear/CMMatrix.kt @@ -3,6 +3,7 @@ package space.kscience.kmath.commons.linear import org.apache.commons.math3.linear.* import space.kscience.kmath.linear.* import space.kscience.kmath.misc.UnstableKMathAPI +import space.kscience.kmath.nd.StructureFeature import space.kscience.kmath.operations.DoubleField import space.kscience.kmath.structures.DoubleBuffer import kotlin.reflect.KClass @@ -89,7 +90,7 @@ public object CMLinearSpace : LinearSpace { v * this @UnstableKMathAPI - override fun getFeature(structure: Matrix, type: KClass): F? { + override fun getFeature(structure: Matrix, type: KClass): F? { //Return the feature if it is intrinsic to the structure structure.getFeature(type)?.let { return it } diff --git a/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/optimization/CMOptimizationProblem.kt b/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/optimization/CMOptimization.kt similarity index 64% rename from kmath-commons/src/main/kotlin/space/kscience/kmath/commons/optimization/CMOptimizationProblem.kt rename to kmath-commons/src/main/kotlin/space/kscience/kmath/commons/optimization/CMOptimization.kt index 13a10475f..444c505c9 100644 --- a/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/optimization/CMOptimizationProblem.kt +++ b/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/optimization/CMOptimization.kt @@ -9,22 +9,36 @@ import org.apache.commons.math3.optim.nonlinear.scalar.gradient.NonLinearConjuga import org.apache.commons.math3.optim.nonlinear.scalar.noderiv.AbstractSimplex import org.apache.commons.math3.optim.nonlinear.scalar.noderiv.NelderMeadSimplex import org.apache.commons.math3.optim.nonlinear.scalar.noderiv.SimplexOptimizer -import space.kscience.kmath.expressions.* -import space.kscience.kmath.stat.OptimizationFeature -import space.kscience.kmath.stat.OptimizationProblem -import space.kscience.kmath.stat.OptimizationProblemFactory -import space.kscience.kmath.stat.OptimizationResult +import space.kscience.kmath.expressions.DifferentiableExpression +import space.kscience.kmath.expressions.Expression +import space.kscience.kmath.expressions.SymbolIndexer +import space.kscience.kmath.expressions.derivative +import space.kscience.kmath.misc.Symbol +import space.kscience.kmath.misc.UnstableKMathAPI +import space.kscience.kmath.optimization.* import kotlin.reflect.KClass public operator fun PointValuePair.component1(): DoubleArray = point public operator fun PointValuePair.component2(): Double = value -public class CMOptimizationProblem(override val symbols: List) : - OptimizationProblem, SymbolIndexer, OptimizationFeature { +@OptIn(UnstableKMathAPI::class) +public class CMOptimization( + override val symbols: List, +) : FunctionOptimization, NoDerivFunctionOptimization, SymbolIndexer, OptimizationFeature { + private val optimizationData: HashMap, OptimizationData> = HashMap() - private var optimizatorBuilder: (() -> MultivariateOptimizer)? = null - public var convergenceChecker: ConvergenceChecker = SimpleValueChecker(DEFAULT_RELATIVE_TOLERANCE, - DEFAULT_ABSOLUTE_TOLERANCE, DEFAULT_MAX_ITER) + private var optimizerBuilder: (() -> MultivariateOptimizer)? = null + public var convergenceChecker: ConvergenceChecker = SimpleValueChecker( + DEFAULT_RELATIVE_TOLERANCE, + DEFAULT_ABSOLUTE_TOLERANCE, + DEFAULT_MAX_ITER + ) + + override var maximize: Boolean + get() = optimizationData[GoalType::class] == GoalType.MAXIMIZE + set(value) { + optimizationData[GoalType::class] = if (value) GoalType.MAXIMIZE else GoalType.MINIMIZE + } public fun addOptimizationData(data: OptimizationData) { optimizationData[data::class] = data @@ -40,7 +54,7 @@ public class CMOptimizationProblem(override val symbols: List) : addOptimizationData(InitialGuess(map.toDoubleArray())) } - public override fun expression(expression: Expression): Unit { + public override fun function(expression: Expression): Unit { val objectiveFunction = ObjectiveFunction { val args = it.toMap() expression(args) @@ -48,8 +62,8 @@ public class CMOptimizationProblem(override val symbols: List) : addOptimizationData(objectiveFunction) } - public override fun diffExpression(expression: DifferentiableExpression>) { - expression(expression) + public override fun diffFunction(expression: DifferentiableExpression>) { + function(expression) val gradientFunction = ObjectiveFunctionGradient { val args = it.toMap() DoubleArray(symbols.size) { index -> @@ -57,8 +71,8 @@ public class CMOptimizationProblem(override val symbols: List) : } } addOptimizationData(gradientFunction) - if (optimizatorBuilder == null) { - optimizatorBuilder = { + if (optimizerBuilder == null) { + optimizerBuilder = { NonLinearConjugateGradientOptimizer( NonLinearConjugateGradientOptimizer.Formula.FLETCHER_REEVES, convergenceChecker @@ -70,8 +84,8 @@ public class CMOptimizationProblem(override val symbols: List) : public fun simplex(simplex: AbstractSimplex) { addOptimizationData(simplex) //Set optimization builder to simplex if it is not present - if (optimizatorBuilder == null) { - optimizatorBuilder = { SimplexOptimizer(convergenceChecker) } + if (optimizerBuilder == null) { + optimizerBuilder = { SimplexOptimizer(convergenceChecker) } } } @@ -84,7 +98,7 @@ public class CMOptimizationProblem(override val symbols: List) : } public fun optimizer(block: () -> MultivariateOptimizer) { - optimizatorBuilder = block + optimizerBuilder = block } override fun update(result: OptimizationResult) { @@ -92,19 +106,19 @@ public class CMOptimizationProblem(override val symbols: List) : } override fun optimize(): OptimizationResult { - val optimizer = optimizatorBuilder?.invoke() ?: error("Optimizer not defined") + val optimizer = optimizerBuilder?.invoke() ?: error("Optimizer not defined") val (point, value) = optimizer.optimize(*optimizationData.values.toTypedArray()) return OptimizationResult(point.toMap(), value, setOf(this)) } - public companion object : OptimizationProblemFactory { + public companion object : OptimizationProblemFactory { public const val DEFAULT_RELATIVE_TOLERANCE: Double = 1e-4 public const val DEFAULT_ABSOLUTE_TOLERANCE: Double = 1e-4 public const val DEFAULT_MAX_ITER: Int = 1000 - override fun build(symbols: List): CMOptimizationProblem = CMOptimizationProblem(symbols) + override fun build(symbols: List): CMOptimization = CMOptimization(symbols) } } -public fun CMOptimizationProblem.initialGuess(vararg pairs: Pair): Unit = initialGuess(pairs.toMap()) -public fun CMOptimizationProblem.simplexSteps(vararg pairs: Pair): Unit = simplexSteps(pairs.toMap()) +public fun CMOptimization.initialGuess(vararg pairs: Pair): Unit = initialGuess(pairs.toMap()) +public fun CMOptimization.simplexSteps(vararg pairs: Pair): Unit = simplexSteps(pairs.toMap()) diff --git a/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/optimization/cmFit.kt b/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/optimization/cmFit.kt index 5ecd5b756..f84dae693 100644 --- a/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/optimization/cmFit.kt +++ b/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/optimization/cmFit.kt @@ -1,21 +1,21 @@ package space.kscience.kmath.commons.optimization import org.apache.commons.math3.analysis.differentiation.DerivativeStructure -import org.apache.commons.math3.optim.nonlinear.scalar.GoalType import space.kscience.kmath.commons.expressions.DerivativeStructureField import space.kscience.kmath.expressions.DifferentiableExpression import space.kscience.kmath.expressions.Expression -import space.kscience.kmath.expressions.Symbol -import space.kscience.kmath.stat.Fitting -import space.kscience.kmath.stat.OptimizationResult -import space.kscience.kmath.stat.optimizeWith +import space.kscience.kmath.misc.Symbol +import space.kscience.kmath.optimization.FunctionOptimization +import space.kscience.kmath.optimization.OptimizationResult +import space.kscience.kmath.optimization.noDerivOptimizeWith +import space.kscience.kmath.optimization.optimizeWith import space.kscience.kmath.structures.Buffer import space.kscience.kmath.structures.asBuffer /** * Generate a chi squared expression from given x-y-sigma data and inline model. Provides automatic differentiation */ -public fun Fitting.chiSquared( +public fun FunctionOptimization.Companion.chiSquared( x: Buffer, y: Buffer, yErr: Buffer, @@ -25,7 +25,7 @@ public fun Fitting.chiSquared( /** * Generate a chi squared expression from given x-y-sigma data and inline model. Provides automatic differentiation */ -public fun Fitting.chiSquared( +public fun FunctionOptimization.Companion.chiSquared( x: Iterable, y: Iterable, yErr: Iterable, @@ -43,25 +43,26 @@ public fun Fitting.chiSquared( */ public fun Expression.optimize( vararg symbols: Symbol, - configuration: CMOptimizationProblem.() -> Unit, -): OptimizationResult = optimizeWith(CMOptimizationProblem, symbols = symbols, configuration) + configuration: CMOptimization.() -> Unit, +): OptimizationResult = noDerivOptimizeWith(CMOptimization, symbols = symbols, configuration) /** * Optimize differentiable expression */ public fun DifferentiableExpression>.optimize( vararg symbols: Symbol, - configuration: CMOptimizationProblem.() -> Unit, -): OptimizationResult = optimizeWith(CMOptimizationProblem, symbols = symbols, configuration) + configuration: CMOptimization.() -> Unit, +): OptimizationResult = optimizeWith(CMOptimization, symbols = symbols, configuration) public fun DifferentiableExpression>.minimize( vararg startPoint: Pair, - configuration: CMOptimizationProblem.() -> Unit = {}, + configuration: CMOptimization.() -> Unit = {}, ): OptimizationResult { - require(startPoint.isNotEmpty()) { "Must provide a list of symbols for optimization" } - val problem = CMOptimizationProblem(startPoint.map { it.first }).apply(configuration) - problem.diffExpression(this) - problem.initialGuess(startPoint.toMap()) - problem.goal(GoalType.MINIMIZE) - return problem.optimize() + val symbols = startPoint.map { it.first }.toTypedArray() + return optimize(*symbols){ + maximize = false + initialGuess(startPoint.toMap()) + diffFunction(this@minimize) + configuration() + } } \ No newline at end of file diff --git a/kmath-commons/src/test/kotlin/space/kscience/kmath/commons/expressions/DerivativeStructureExpressionTest.kt b/kmath-commons/src/test/kotlin/space/kscience/kmath/commons/expressions/DerivativeStructureExpressionTest.kt index 8d9bab652..b19eb5950 100644 --- a/kmath-commons/src/test/kotlin/space/kscience/kmath/commons/expressions/DerivativeStructureExpressionTest.kt +++ b/kmath-commons/src/test/kotlin/space/kscience/kmath/commons/expressions/DerivativeStructureExpressionTest.kt @@ -1,6 +1,10 @@ package space.kscience.kmath.commons.expressions -import space.kscience.kmath.expressions.* +import space.kscience.kmath.expressions.binding +import space.kscience.kmath.expressions.derivative +import space.kscience.kmath.expressions.invoke +import space.kscience.kmath.misc.Symbol +import space.kscience.kmath.misc.symbol import kotlin.contracts.InvocationKind import kotlin.contracts.contract import kotlin.test.Test diff --git a/kmath-commons/src/test/kotlin/space/kscience/kmath/commons/optimization/OptimizeTest.kt b/kmath-commons/src/test/kotlin/space/kscience/kmath/commons/optimization/OptimizeTest.kt index d29934a4d..36f2639f4 100644 --- a/kmath-commons/src/test/kotlin/space/kscience/kmath/commons/optimization/OptimizeTest.kt +++ b/kmath-commons/src/test/kotlin/space/kscience/kmath/commons/optimization/OptimizeTest.kt @@ -1,13 +1,13 @@ package space.kscience.kmath.commons.optimization -import org.junit.jupiter.api.Test +import kotlinx.coroutines.runBlocking import space.kscience.kmath.commons.expressions.DerivativeStructureExpression -import space.kscience.kmath.expressions.symbol -import space.kscience.kmath.stat.Distribution -import space.kscience.kmath.stat.Fitting +import space.kscience.kmath.misc.symbol +import space.kscience.kmath.optimization.FunctionOptimization import space.kscience.kmath.stat.RandomGenerator -import space.kscience.kmath.stat.normal +import space.kscience.kmath.stat.distributions.NormalDistribution import kotlin.math.pow +import kotlin.test.Test internal class OptimizeTest { val x by symbol @@ -34,28 +34,29 @@ internal class OptimizeTest { simplexSteps(x to 2.0, y to 0.5) //this sets simplex optimizer } + println(result.point) println(result.value) } @Test - fun testCmFit() { + fun testCmFit() = runBlocking { val a by symbol val b by symbol val c by symbol val sigma = 1.0 - val generator = Distribution.normal(0.0, sigma) + val generator = NormalDistribution(0.0, sigma) val chain = generator.sample(RandomGenerator.default(112667)) val x = (1..100).map(Int::toDouble) val y = x.map { - it.pow(2) + it + 1 + chain.nextDouble() + it.pow(2) + it + 1 + chain.next() } val yErr = List(x.size) { sigma } - val chi2 = Fitting.chiSquared(x, y, yErr) { x1 -> + val chi2 = FunctionOptimization.chiSquared(x, y, yErr) { x1 -> val cWithDefault = bindSymbolOrNull(c) ?: one bind(a) * x1.pow(2) + bind(b) * x1 + cWithDefault } @@ -64,5 +65,4 @@ internal class OptimizeTest { println(result) println("Chi2/dof = ${result.value / (x.size - 3)}") } - -} \ No newline at end of file +} diff --git a/kmath-complex/README.md b/kmath-complex/README.md index 9e9cd5b6f..d7b2937fd 100644 --- a/kmath-complex/README.md +++ b/kmath-complex/README.md @@ -1,42 +1,36 @@ -# The Core Module (`kmath-core`) +# Module kmath-complex -Complex and hypercomplex number systems in KMath: +Complex and hypercomplex number systems in KMath. - [complex](src/commonMain/kotlin/space/kscience/kmath/complex/Complex.kt) : Complex Numbers - [quaternion](src/commonMain/kotlin/space/kscience/kmath/complex/Quaternion.kt) : Quaternions -> #### Artifact: -> -> This module artifact: `space.kscience:kmath-complex:0.3.0-dev-3`. -> -> Bintray release version: [ ![Download](https://api.bintray.com/packages/mipt-npm/kscience/kmath-complex/images/download.svg) ](https://bintray.com/mipt-npm/kscience/kmath-complex/_latestVersion) -> -> Bintray development version: [ ![Download](https://api.bintray.com/packages/mipt-npm/dev/kmath-complex/images/download.svg) ](https://bintray.com/mipt-npm/dev/kmath-complex/_latestVersion) -> -> **Gradle:** -> -> ```gradle -> repositories { -> maven { url 'https://repo.kotlin.link' } -> maven { url 'https://dl.bintray.com/hotkeytlt/maven' } -> maven { url "https://dl.bintray.com/kotlin/kotlin-eap" } // include for builds based on kotlin-eap -> } -> -> dependencies { -> implementation 'space.kscience:kmath-complex:0.3.0-dev-3' -> } -> ``` -> **Gradle Kotlin DSL:** -> -> ```kotlin -> repositories { -> maven("https://repo.kotlin.link") -> maven("https://dl.bintray.com/kotlin/kotlin-eap") // include for builds based on kotlin-eap -> maven("https://dl.bintray.com/hotkeytlt/maven") // required for a -> } -> -> dependencies { -> implementation("space.kscience:kmath-complex:0.3.0-dev-3") -> } -> ``` +## Artifact: + +The Maven coordinates of this project are `space.kscience:kmath-complex:0.3.0-dev-3`. + +**Gradle:** +```gradle +repositories { + maven { url 'https://repo.kotlin.link' } + maven { url 'https://dl.bintray.com/hotkeytlt/maven' } + maven { url "https://dl.bintray.com/kotlin/kotlin-eap" } // include for builds based on kotlin-eap +} + +dependencies { + implementation 'space.kscience:kmath-complex:0.3.0-dev-3' +} +``` +**Gradle Kotlin DSL:** +```kotlin +repositories { + maven("https://repo.kotlin.link") + maven("https://dl.bintray.com/kotlin/kotlin-eap") // include for builds based on kotlin-eap + maven("https://dl.bintray.com/hotkeytlt/maven") // required for a +} + +dependencies { + implementation("space.kscience:kmath-complex:0.3.0-dev-3") +} +``` diff --git a/kmath-complex/docs/README-TEMPLATE.md b/kmath-complex/docs/README-TEMPLATE.md index 462fd617e..106d4aff1 100644 --- a/kmath-complex/docs/README-TEMPLATE.md +++ b/kmath-complex/docs/README-TEMPLATE.md @@ -1,6 +1,6 @@ -# The Core Module (`kmath-core`) +# Module kmath-complex -Complex and hypercomplex number systems in KMath: +Complex and hypercomplex number systems in KMath. ${features} diff --git a/kmath-complex/src/commonTest/kotlin/space/kscience/kmath/complex/ExpressionFieldForComplexTest.kt b/kmath-complex/src/commonTest/kotlin/space/kscience/kmath/complex/ExpressionFieldForComplexTest.kt index 81a131318..3837b0d40 100644 --- a/kmath-complex/src/commonTest/kotlin/space/kscience/kmath/complex/ExpressionFieldForComplexTest.kt +++ b/kmath-complex/src/commonTest/kotlin/space/kscience/kmath/complex/ExpressionFieldForComplexTest.kt @@ -3,7 +3,7 @@ package space.kscience.kmath.complex import space.kscience.kmath.expressions.FunctionalExpressionField import space.kscience.kmath.expressions.bindSymbol import space.kscience.kmath.expressions.invoke -import space.kscience.kmath.expressions.symbol +import space.kscience.kmath.misc.symbol import kotlin.test.Test import kotlin.test.assertEquals diff --git a/kmath-core/README.md b/kmath-core/README.md index 4e4b5273d..096c7d833 100644 --- a/kmath-core/README.md +++ b/kmath-core/README.md @@ -1,6 +1,6 @@ -# The Core Module (`kmath-core`) +# Module kmath-core -The core features of KMath: +The core interfaces of KMath. - [algebras](src/commonMain/kotlin/space/kscience/kmath/operations/Algebra.kt) : Algebraic structures like rings, spaces and fields. - [nd](src/commonMain/kotlin/space/kscience/kmath/structures/StructureND.kt) : Many-dimensional structures and operations on them. @@ -13,37 +13,31 @@ performance calculations to code generation. - [autodif](src/commonMain/kotlin/space/kscience/kmath/expressions/SimpleAutoDiff.kt) : Automatic differentiation -> #### Artifact: -> -> This module artifact: `space.kscience:kmath-core:0.3.0-dev-3`. -> -> Bintray release version: [ ![Download](https://api.bintray.com/packages/mipt-npm/kscience/kmath-core/images/download.svg) ](https://bintray.com/mipt-npm/kscience/kmath-core/_latestVersion) -> -> Bintray development version: [ ![Download](https://api.bintray.com/packages/mipt-npm/dev/kmath-core/images/download.svg) ](https://bintray.com/mipt-npm/dev/kmath-core/_latestVersion) -> -> **Gradle:** -> -> ```gradle -> repositories { -> maven { url 'https://repo.kotlin.link' } -> maven { url 'https://dl.bintray.com/hotkeytlt/maven' } -> maven { url "https://dl.bintray.com/kotlin/kotlin-eap" } // include for builds based on kotlin-eap -> } -> -> dependencies { -> implementation 'space.kscience:kmath-core:0.3.0-dev-3' -> } -> ``` -> **Gradle Kotlin DSL:** -> -> ```kotlin -> repositories { -> maven("https://repo.kotlin.link") -> maven("https://dl.bintray.com/kotlin/kotlin-eap") // include for builds based on kotlin-eap -> maven("https://dl.bintray.com/hotkeytlt/maven") // required for a -> } -> -> dependencies { -> implementation("space.kscience:kmath-core:0.3.0-dev-3") -> } -> ``` +## Artifact: + +The Maven coordinates of this project are `space.kscience:kmath-core:0.3.0-dev-3`. + +**Gradle:** +```gradle +repositories { + maven { url 'https://repo.kotlin.link' } + maven { url 'https://dl.bintray.com/hotkeytlt/maven' } + maven { url "https://dl.bintray.com/kotlin/kotlin-eap" } // include for builds based on kotlin-eap +} + +dependencies { + implementation 'space.kscience:kmath-core:0.3.0-dev-3' +} +``` +**Gradle Kotlin DSL:** +```kotlin +repositories { + maven("https://repo.kotlin.link") + maven("https://dl.bintray.com/kotlin/kotlin-eap") // include for builds based on kotlin-eap + maven("https://dl.bintray.com/hotkeytlt/maven") // required for a +} + +dependencies { + implementation("space.kscience:kmath-core:0.3.0-dev-3") +} +``` diff --git a/kmath-core/api/kmath-core.api b/kmath-core/api/kmath-core.api index f0e749d2b..cc6c93c20 100644 --- a/kmath-core/api/kmath-core.api +++ b/kmath-core/api/kmath-core.api @@ -1,3 +1,18 @@ +public final class space/kscience/kmath/data/ColumnarDataKt { +} + +public final class space/kscience/kmath/data/XYColumnarData$DefaultImpls { + public static fun get (Lspace/kscience/kmath/data/XYColumnarData;Lspace/kscience/kmath/misc/Symbol;)Lspace/kscience/kmath/structures/Buffer; +} + +public final class space/kscience/kmath/data/XYColumnarDataKt { + public static synthetic fun asXYData$default (Lspace/kscience/kmath/nd/Structure2D;IIILjava/lang/Object;)Lspace/kscience/kmath/data/XYColumnarData; +} + +public final class space/kscience/kmath/data/XYZColumnarData$DefaultImpls { + public static fun get (Lspace/kscience/kmath/data/XYZColumnarData;Lspace/kscience/kmath/misc/Symbol;)Lspace/kscience/kmath/structures/Buffer; +} + public abstract interface class space/kscience/kmath/domains/Domain { public abstract fun contains (Lspace/kscience/kmath/structures/Buffer;)Z public abstract fun getDimension ()I @@ -14,7 +29,7 @@ public class space/kscience/kmath/expressions/AutoDiffValue { public final class space/kscience/kmath/expressions/DerivationResult { public fun (Ljava/lang/Object;Ljava/util/Map;Lspace/kscience/kmath/operations/Field;)V - public final fun derivative (Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object; + public final fun derivative (Lspace/kscience/kmath/misc/Symbol;)Ljava/lang/Object; public final fun div ()Ljava/lang/Object; public final fun getContext ()Lspace/kscience/kmath/operations/Field; public final fun getValue ()Ljava/lang/Object; @@ -27,7 +42,7 @@ public abstract interface class space/kscience/kmath/expressions/DifferentiableE public final class space/kscience/kmath/expressions/DifferentiableExpressionKt { public static final fun derivative (Lspace/kscience/kmath/expressions/DifferentiableExpression;Ljava/lang/String;)Lspace/kscience/kmath/expressions/Expression; public static final fun derivative (Lspace/kscience/kmath/expressions/DifferentiableExpression;Ljava/util/List;)Lspace/kscience/kmath/expressions/Expression; - public static final fun derivative (Lspace/kscience/kmath/expressions/DifferentiableExpression;[Lspace/kscience/kmath/expressions/Symbol;)Lspace/kscience/kmath/expressions/Expression; + public static final fun derivative (Lspace/kscience/kmath/expressions/DifferentiableExpression;[Lspace/kscience/kmath/misc/Symbol;)Lspace/kscience/kmath/expressions/Expression; } public abstract interface class space/kscience/kmath/expressions/Expression { @@ -36,7 +51,7 @@ public abstract interface class space/kscience/kmath/expressions/Expression { public abstract interface class space/kscience/kmath/expressions/ExpressionAlgebra : space/kscience/kmath/operations/Algebra { public abstract fun bindSymbol (Ljava/lang/String;)Ljava/lang/Object; - public abstract fun bindSymbolOrNull (Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object; + public abstract fun bindSymbolOrNull (Lspace/kscience/kmath/misc/Symbol;)Ljava/lang/Object; public abstract fun const (Ljava/lang/Object;)Ljava/lang/Object; } @@ -56,18 +71,17 @@ public final class space/kscience/kmath/expressions/ExpressionBuildersKt { } public final class space/kscience/kmath/expressions/ExpressionKt { - public static final fun bindSymbol (Lspace/kscience/kmath/expressions/ExpressionAlgebra;Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object; + public static final fun bindSymbol (Lspace/kscience/kmath/expressions/ExpressionAlgebra;Lspace/kscience/kmath/misc/Symbol;)Ljava/lang/Object; public static final fun binding (Lspace/kscience/kmath/expressions/ExpressionAlgebra;)Lkotlin/properties/ReadOnlyProperty; public static final fun callByString (Lspace/kscience/kmath/expressions/Expression;[Lkotlin/Pair;)Ljava/lang/Object; public static final fun callBySymbol (Lspace/kscience/kmath/expressions/Expression;[Lkotlin/Pair;)Ljava/lang/Object; - public static final fun getSymbol ()Lkotlin/properties/ReadOnlyProperty; public static final fun invoke (Lspace/kscience/kmath/expressions/Expression;)Ljava/lang/Object; } public abstract class space/kscience/kmath/expressions/FirstDerivativeExpression : space/kscience/kmath/expressions/DifferentiableExpression { public fun ()V public final fun derivativeOrNull (Ljava/util/List;)Lspace/kscience/kmath/expressions/Expression; - public abstract fun derivativeOrNull (Lspace/kscience/kmath/expressions/Symbol;)Lspace/kscience/kmath/expressions/Expression; + public abstract fun derivativeOrNull (Lspace/kscience/kmath/misc/Symbol;)Lspace/kscience/kmath/expressions/Expression; } public abstract class space/kscience/kmath/expressions/FunctionalExpressionAlgebra : space/kscience/kmath/expressions/ExpressionAlgebra { @@ -77,8 +91,8 @@ public abstract class space/kscience/kmath/expressions/FunctionalExpressionAlgeb public fun binaryOperationFunction (Ljava/lang/String;)Lkotlin/jvm/functions/Function2; public synthetic fun bindSymbol (Ljava/lang/String;)Ljava/lang/Object; public fun bindSymbol (Ljava/lang/String;)Lspace/kscience/kmath/expressions/Expression; - public synthetic fun bindSymbolOrNull (Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object; - public fun bindSymbolOrNull (Lspace/kscience/kmath/expressions/Symbol;)Lspace/kscience/kmath/expressions/Expression; + public synthetic fun bindSymbolOrNull (Lspace/kscience/kmath/misc/Symbol;)Ljava/lang/Object; + public fun bindSymbolOrNull (Lspace/kscience/kmath/misc/Symbol;)Lspace/kscience/kmath/expressions/Expression; public synthetic fun const (Ljava/lang/Object;)Ljava/lang/Object; public fun const (Ljava/lang/Object;)Lspace/kscience/kmath/expressions/Expression; public final fun getAlgebra ()Lspace/kscience/kmath/operations/Algebra; @@ -203,7 +217,7 @@ public class space/kscience/kmath/expressions/FunctionalExpressionRing : space/k public final class space/kscience/kmath/expressions/SimpleAutoDiffExpression : space/kscience/kmath/expressions/FirstDerivativeExpression { public fun (Lspace/kscience/kmath/operations/Field;Lkotlin/jvm/functions/Function1;)V - public fun derivativeOrNull (Lspace/kscience/kmath/expressions/Symbol;)Lspace/kscience/kmath/expressions/Expression; + public fun derivativeOrNull (Lspace/kscience/kmath/misc/Symbol;)Lspace/kscience/kmath/expressions/Expression; public final fun getField ()Lspace/kscience/kmath/operations/Field; public final fun getFunction ()Lkotlin/jvm/functions/Function1; public fun invoke (Ljava/util/Map;)Ljava/lang/Object; @@ -264,8 +278,8 @@ public class space/kscience/kmath/expressions/SimpleAutoDiffField : space/kscien public fun binaryOperationFunction (Ljava/lang/String;)Lkotlin/jvm/functions/Function2; public synthetic fun bindSymbol (Ljava/lang/String;)Ljava/lang/Object; public fun bindSymbol (Ljava/lang/String;)Lspace/kscience/kmath/expressions/AutoDiffValue; - public synthetic fun bindSymbolOrNull (Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object; - public fun bindSymbolOrNull (Lspace/kscience/kmath/expressions/Symbol;)Lspace/kscience/kmath/expressions/AutoDiffValue; + public synthetic fun bindSymbolOrNull (Lspace/kscience/kmath/misc/Symbol;)Ljava/lang/Object; + public fun bindSymbolOrNull (Lspace/kscience/kmath/misc/Symbol;)Lspace/kscience/kmath/expressions/AutoDiffValue; public synthetic fun const (Ljava/lang/Object;)Ljava/lang/Object; public fun const (Ljava/lang/Object;)Lspace/kscience/kmath/expressions/AutoDiffValue; public final fun const (Lkotlin/jvm/functions/Function1;)Lspace/kscience/kmath/expressions/AutoDiffValue; @@ -332,7 +346,7 @@ public final class space/kscience/kmath/expressions/SimpleAutoDiffKt { public static final fun cos (Lspace/kscience/kmath/expressions/SimpleAutoDiffField;Lspace/kscience/kmath/expressions/AutoDiffValue;)Lspace/kscience/kmath/expressions/AutoDiffValue; public static final fun cosh (Lspace/kscience/kmath/expressions/SimpleAutoDiffField;Lspace/kscience/kmath/expressions/AutoDiffValue;)Lspace/kscience/kmath/expressions/AutoDiffValue; public static final fun exp (Lspace/kscience/kmath/expressions/SimpleAutoDiffField;Lspace/kscience/kmath/expressions/AutoDiffValue;)Lspace/kscience/kmath/expressions/AutoDiffValue; - public static final fun grad (Lspace/kscience/kmath/expressions/DerivationResult;[Lspace/kscience/kmath/expressions/Symbol;)Lspace/kscience/kmath/structures/Buffer; + public static final fun grad (Lspace/kscience/kmath/expressions/DerivationResult;[Lspace/kscience/kmath/misc/Symbol;)Lspace/kscience/kmath/structures/Buffer; public static final fun ln (Lspace/kscience/kmath/expressions/SimpleAutoDiffField;Lspace/kscience/kmath/expressions/AutoDiffValue;)Lspace/kscience/kmath/expressions/AutoDiffValue; public static final fun pow (Lspace/kscience/kmath/expressions/SimpleAutoDiffField;Lspace/kscience/kmath/expressions/AutoDiffValue;D)Lspace/kscience/kmath/expressions/AutoDiffValue; public static final fun pow (Lspace/kscience/kmath/expressions/SimpleAutoDiffField;Lspace/kscience/kmath/expressions/AutoDiffValue;I)Lspace/kscience/kmath/expressions/AutoDiffValue; @@ -348,79 +362,13 @@ public final class space/kscience/kmath/expressions/SimpleAutoDiffKt { public static final fun tanh (Lspace/kscience/kmath/expressions/SimpleAutoDiffField;Lspace/kscience/kmath/expressions/AutoDiffValue;)Lspace/kscience/kmath/expressions/AutoDiffValue; } -public final class space/kscience/kmath/expressions/SimpleSymbolIndexer : space/kscience/kmath/expressions/SymbolIndexer { - public static final synthetic fun box-impl (Ljava/util/List;)Lspace/kscience/kmath/expressions/SimpleSymbolIndexer; - public static fun constructor-impl (Ljava/util/List;)Ljava/util/List; - public fun equals (Ljava/lang/Object;)Z - public static fun equals-impl (Ljava/util/List;Ljava/lang/Object;)Z - public static final fun equals-impl0 (Ljava/util/List;Ljava/util/List;)Z - public fun get (Ljava/util/List;Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object; - public fun get (Lspace/kscience/kmath/nd/Structure2D;Lspace/kscience/kmath/expressions/Symbol;Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object; - public fun get (Lspace/kscience/kmath/structures/Buffer;Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object; - public fun get ([DLspace/kscience/kmath/expressions/Symbol;)D - public fun get ([Ljava/lang/Object;Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object; - public static fun get-impl (Ljava/util/List;Ljava/util/List;Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object; - public static fun get-impl (Ljava/util/List;Lspace/kscience/kmath/nd/Structure2D;Lspace/kscience/kmath/expressions/Symbol;Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object; - public static fun get-impl (Ljava/util/List;Lspace/kscience/kmath/structures/Buffer;Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object; - public static fun get-impl (Ljava/util/List;[DLspace/kscience/kmath/expressions/Symbol;)D - public static fun get-impl (Ljava/util/List;[Ljava/lang/Object;Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object; - public fun getSymbols ()Ljava/util/List; - public fun hashCode ()I - public static fun hashCode-impl (Ljava/util/List;)I - public fun indexOf (Lspace/kscience/kmath/expressions/Symbol;)I - public static fun indexOf-impl (Ljava/util/List;Lspace/kscience/kmath/expressions/Symbol;)I - public fun toDoubleArray (Ljava/util/Map;)[D - public static fun toDoubleArray-impl (Ljava/util/List;Ljava/util/Map;)[D - public fun toList (Ljava/util/Map;)Ljava/util/List; - public static fun toList-impl (Ljava/util/List;Ljava/util/Map;)Ljava/util/List; - public fun toMap ([D)Ljava/util/Map; - public static fun toMap-impl (Ljava/util/List;[D)Ljava/util/Map; - public fun toPoint (Ljava/util/Map;Lkotlin/jvm/functions/Function2;)Lspace/kscience/kmath/structures/Buffer; - public static fun toPoint-impl (Ljava/util/List;Ljava/util/Map;Lkotlin/jvm/functions/Function2;)Lspace/kscience/kmath/structures/Buffer; - public fun toString ()Ljava/lang/String; - public static fun toString-impl (Ljava/util/List;)Ljava/lang/String; - public final synthetic fun unbox-impl ()Ljava/util/List; -} - -public final class space/kscience/kmath/expressions/StringSymbol : space/kscience/kmath/expressions/Symbol { - public static final synthetic fun box-impl (Ljava/lang/String;)Lspace/kscience/kmath/expressions/StringSymbol; - public static fun constructor-impl (Ljava/lang/String;)Ljava/lang/String; - public fun equals (Ljava/lang/Object;)Z - public static fun equals-impl (Ljava/lang/String;Ljava/lang/Object;)Z - public static final fun equals-impl0 (Ljava/lang/String;Ljava/lang/String;)Z - public fun getIdentity ()Ljava/lang/String; - public fun hashCode ()I - public static fun hashCode-impl (Ljava/lang/String;)I - public fun toString ()Ljava/lang/String; - public static fun toString-impl (Ljava/lang/String;)Ljava/lang/String; - public final synthetic fun unbox-impl ()Ljava/lang/String; -} - -public abstract interface class space/kscience/kmath/expressions/Symbol { - public abstract fun getIdentity ()Ljava/lang/String; -} - -public abstract interface class space/kscience/kmath/expressions/SymbolIndexer { - public abstract fun get (Ljava/util/List;Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object; - public abstract fun get (Lspace/kscience/kmath/nd/Structure2D;Lspace/kscience/kmath/expressions/Symbol;Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object; - public abstract fun get (Lspace/kscience/kmath/structures/Buffer;Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object; - public abstract fun get ([DLspace/kscience/kmath/expressions/Symbol;)D - public abstract fun get ([Ljava/lang/Object;Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object; - public abstract fun getSymbols ()Ljava/util/List; - public abstract fun indexOf (Lspace/kscience/kmath/expressions/Symbol;)I - public abstract fun toDoubleArray (Ljava/util/Map;)[D - public abstract fun toList (Ljava/util/Map;)Ljava/util/List; - public abstract fun toMap ([D)Ljava/util/Map; - public abstract fun toPoint (Ljava/util/Map;Lkotlin/jvm/functions/Function2;)Lspace/kscience/kmath/structures/Buffer; -} - public final class space/kscience/kmath/expressions/SymbolIndexer$DefaultImpls { - public static fun get (Lspace/kscience/kmath/expressions/SymbolIndexer;Ljava/util/List;Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object; - public static fun get (Lspace/kscience/kmath/expressions/SymbolIndexer;Lspace/kscience/kmath/nd/Structure2D;Lspace/kscience/kmath/expressions/Symbol;Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object; - public static fun get (Lspace/kscience/kmath/expressions/SymbolIndexer;Lspace/kscience/kmath/structures/Buffer;Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object; - public static fun get (Lspace/kscience/kmath/expressions/SymbolIndexer;[DLspace/kscience/kmath/expressions/Symbol;)D - public static fun get (Lspace/kscience/kmath/expressions/SymbolIndexer;[Ljava/lang/Object;Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object; - public static fun indexOf (Lspace/kscience/kmath/expressions/SymbolIndexer;Lspace/kscience/kmath/expressions/Symbol;)I + public static fun get (Lspace/kscience/kmath/expressions/SymbolIndexer;Ljava/util/List;Lspace/kscience/kmath/misc/Symbol;)Ljava/lang/Object; + public static fun get (Lspace/kscience/kmath/expressions/SymbolIndexer;Lspace/kscience/kmath/nd/Structure2D;Lspace/kscience/kmath/misc/Symbol;Lspace/kscience/kmath/misc/Symbol;)Ljava/lang/Object; + public static fun get (Lspace/kscience/kmath/expressions/SymbolIndexer;Lspace/kscience/kmath/structures/Buffer;Lspace/kscience/kmath/misc/Symbol;)Ljava/lang/Object; + public static fun get (Lspace/kscience/kmath/expressions/SymbolIndexer;[DLspace/kscience/kmath/misc/Symbol;)D + public static fun get (Lspace/kscience/kmath/expressions/SymbolIndexer;[Ljava/lang/Object;Lspace/kscience/kmath/misc/Symbol;)Ljava/lang/Object; + public static fun indexOf (Lspace/kscience/kmath/expressions/SymbolIndexer;Lspace/kscience/kmath/misc/Symbol;)I public static fun toDoubleArray (Lspace/kscience/kmath/expressions/SymbolIndexer;Ljava/util/Map;)[D public static fun toList (Lspace/kscience/kmath/expressions/SymbolIndexer;Ljava/util/Map;)Ljava/util/List; public static fun toMap (Lspace/kscience/kmath/expressions/SymbolIndexer;[D)Ljava/util/Map; @@ -428,8 +376,6 @@ public final class space/kscience/kmath/expressions/SymbolIndexer$DefaultImpls { } public final class space/kscience/kmath/expressions/SymbolIndexerKt { - public static final fun withSymbols (Ljava/util/Collection;Lkotlin/jvm/functions/Function1;)Ljava/lang/Object; - public static final fun withSymbols ([Lspace/kscience/kmath/expressions/Symbol;Lkotlin/jvm/functions/Function1;)Ljava/lang/Object; } public final class space/kscience/kmath/linear/BufferedLinearSpace : space/kscience/kmath/linear/LinearSpace { @@ -575,7 +521,7 @@ public final class space/kscience/kmath/linear/MatrixBuilderKt { public static final fun row (Lspace/kscience/kmath/linear/LinearSpace;[Ljava/lang/Object;)Lspace/kscience/kmath/nd/Structure2D; } -public abstract interface class space/kscience/kmath/linear/MatrixFeature { +public abstract interface class space/kscience/kmath/linear/MatrixFeature : space/kscience/kmath/nd/StructureFeature { } public final class space/kscience/kmath/linear/MatrixFeaturesKt { @@ -672,6 +618,35 @@ public final class space/kscience/kmath/misc/CumulativeKt { public static final fun cumulativeSumOfLong (Lkotlin/sequences/Sequence;)Lkotlin/sequences/Sequence; } +public final class space/kscience/kmath/misc/StringSymbol : space/kscience/kmath/misc/Symbol { + public static final synthetic fun box-impl (Ljava/lang/String;)Lspace/kscience/kmath/misc/StringSymbol; + public static fun constructor-impl (Ljava/lang/String;)Ljava/lang/String; + public fun equals (Ljava/lang/Object;)Z + public static fun equals-impl (Ljava/lang/String;Ljava/lang/Object;)Z + public static final fun equals-impl0 (Ljava/lang/String;Ljava/lang/String;)Z + public fun getIdentity ()Ljava/lang/String; + public fun hashCode ()I + public static fun hashCode-impl (Ljava/lang/String;)I + public fun toString ()Ljava/lang/String; + public static fun toString-impl (Ljava/lang/String;)Ljava/lang/String; + public final synthetic fun unbox-impl ()Ljava/lang/String; +} + +public abstract interface class space/kscience/kmath/misc/Symbol { + public static final field Companion Lspace/kscience/kmath/misc/Symbol$Companion; + public abstract fun getIdentity ()Ljava/lang/String; +} + +public final class space/kscience/kmath/misc/Symbol$Companion { + public final fun getX-tWtZOCg ()Ljava/lang/String; + public final fun getY-tWtZOCg ()Ljava/lang/String; + public final fun getZ-tWtZOCg ()Ljava/lang/String; +} + +public final class space/kscience/kmath/misc/SymbolKt { + public static final fun getSymbol ()Lkotlin/properties/ReadOnlyProperty; +} + public abstract interface annotation class space/kscience/kmath/misc/UnstableKMathAPI : java/lang/annotation/Annotation { } @@ -1086,11 +1061,15 @@ public final class space/kscience/kmath/nd/Strides$DefaultImpls { } public abstract interface class space/kscience/kmath/nd/Structure1D : space/kscience/kmath/nd/StructureND, space/kscience/kmath/structures/Buffer { + public static final field Companion Lspace/kscience/kmath/nd/Structure1D$Companion; public abstract fun get ([I)Ljava/lang/Object; public abstract fun getDimension ()I public abstract fun iterator ()Ljava/util/Iterator; } +public final class space/kscience/kmath/nd/Structure1D$Companion { +} + public final class space/kscience/kmath/nd/Structure1D$DefaultImpls { public static fun get (Lspace/kscience/kmath/nd/Structure1D;[I)Ljava/lang/Object; public static fun getDimension (Lspace/kscience/kmath/nd/Structure1D;)I @@ -1132,6 +1111,9 @@ public final class space/kscience/kmath/nd/Structure2DKt { public static final fun as2D (Lspace/kscience/kmath/nd/StructureND;)Lspace/kscience/kmath/nd/Structure2D; } +public abstract interface class space/kscience/kmath/nd/StructureFeature { +} + public abstract interface class space/kscience/kmath/nd/StructureND { public static final field Companion Lspace/kscience/kmath/nd/StructureND$Companion; public abstract fun elements ()Lkotlin/sequences/Sequence; @@ -1201,7 +1183,7 @@ public final class space/kscience/kmath/operations/AlgebraExtensionsKt { } public final class space/kscience/kmath/operations/AlgebraKt { - public static final fun bindSymbol (Lspace/kscience/kmath/operations/Algebra;Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object; + public static final fun bindSymbol (Lspace/kscience/kmath/operations/Algebra;Lspace/kscience/kmath/misc/Symbol;)Ljava/lang/Object; public static final fun invoke (Lspace/kscience/kmath/operations/Algebra;Lkotlin/jvm/functions/Function1;)Ljava/lang/Object; } @@ -2140,6 +2122,7 @@ public final class space/kscience/kmath/structures/BufferKt { public final class space/kscience/kmath/structures/BufferOperationKt { public static final fun asIterable (Lspace/kscience/kmath/structures/Buffer;)Ljava/lang/Iterable; public static final fun asSequence (Lspace/kscience/kmath/structures/Buffer;)Lkotlin/sequences/Sequence; + public static final fun fold (Lspace/kscience/kmath/structures/Buffer;Ljava/lang/Object;Lkotlin/jvm/functions/Function2;)Ljava/lang/Object; public static final fun toList (Lspace/kscience/kmath/structures/Buffer;)Ljava/util/List; } diff --git a/kmath-core/docs/README-TEMPLATE.md b/kmath-core/docs/README-TEMPLATE.md index 83d1ebdce..41cfe1ccb 100644 --- a/kmath-core/docs/README-TEMPLATE.md +++ b/kmath-core/docs/README-TEMPLATE.md @@ -1,6 +1,6 @@ -# The Core Module (`kmath-core`) +# Module kmath-core -The core features of KMath: +The core interfaces of KMath. ${features} diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/data/ColumnarData.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/data/ColumnarData.kt new file mode 100644 index 000000000..761255158 --- /dev/null +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/data/ColumnarData.kt @@ -0,0 +1,34 @@ +package space.kscience.kmath.data + +import space.kscience.kmath.misc.Symbol +import space.kscience.kmath.misc.UnstableKMathAPI +import space.kscience.kmath.nd.Structure2D +import space.kscience.kmath.structures.Buffer + +/** + * A column-based data set with all columns of the same size (not necessary fixed in time). + * The column could be retrieved by a [get] operation. + */ +@UnstableKMathAPI +public interface ColumnarData { + public val size: Int + + public operator fun get(symbol: Symbol): Buffer +} + +/** + * A zero-copy method to represent a [Structure2D] as a two-column x-y data. + * There could more than two columns in the structure. + */ +@UnstableKMathAPI +public fun Structure2D.asColumnarData(mapping: Map): ColumnarData { + require(shape[1] >= mapping.maxOf { it.value }) { "Column index out of bounds" } + return object : ColumnarData { + override val size: Int get() = shape[0] + override fun get(symbol: Symbol): Buffer { + val index = mapping[symbol] ?: error("No column mapping for symbol $symbol") + return columns[index] + } + } +} + diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/data/XYColumnarData.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/data/XYColumnarData.kt new file mode 100644 index 000000000..15239bca1 --- /dev/null +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/data/XYColumnarData.kt @@ -0,0 +1,55 @@ +package space.kscience.kmath.data + +import space.kscience.kmath.misc.Symbol +import space.kscience.kmath.misc.UnstableKMathAPI +import space.kscience.kmath.nd.Structure2D +import space.kscience.kmath.structures.Buffer +import kotlin.math.max + +/** + * The buffer of X values. + */ +@UnstableKMathAPI +public interface XYColumnarData : ColumnarData { + /** + * The buffer of X values + */ + public val x: Buffer + + /** + * The buffer of Y values. + */ + public val y: Buffer + + override fun get(symbol: Symbol): Buffer = when (symbol) { + Symbol.x -> x + Symbol.y -> y + else -> error("A column for symbol $symbol not found") + } +} + +@Suppress("FunctionName") +@UnstableKMathAPI +public fun XYColumnarData(x: Buffer, y: Buffer): XYColumnarData { + require(x.size == y.size) { "Buffer size mismatch. x buffer size is ${x.size}, y buffer size is ${y.size}" } + return object : XYColumnarData { + override val size: Int = x.size + override val x: Buffer = x + override val y: Buffer = y + } +} + + +/** + * A zero-copy method to represent a [Structure2D] as a two-column x-y data. + * There could more than two columns in the structure. + */ +@UnstableKMathAPI +public fun Structure2D.asXYData(xIndex: Int = 0, yIndex: Int = 1): XYColumnarData { + require(shape[1] >= max(xIndex, yIndex)) { "Column index out of bounds" } + return object : XYColumnarData { + override val size: Int get() = this@asXYData.shape[0] + override val x: Buffer get() = columns[xIndex] + override val y: Buffer get() = columns[yIndex] + } +} diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/data/XYZColumnarData.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/data/XYZColumnarData.kt new file mode 100644 index 000000000..f74c6e2d6 --- /dev/null +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/data/XYZColumnarData.kt @@ -0,0 +1,21 @@ +package space.kscience.kmath.data + +import space.kscience.kmath.misc.Symbol +import space.kscience.kmath.misc.UnstableKMathAPI +import space.kscience.kmath.structures.Buffer + +/** + * A [XYColumnarData] with guaranteed [x], [y] and [z] columns designated by corresponding symbols. + * Inherits [XYColumnarData]. + */ +@UnstableKMathAPI +public interface XYZColumnarData : XYColumnarData { + public val z: Buffer + + override fun get(symbol: Symbol): Buffer = when (symbol) { + Symbol.x -> x + Symbol.y -> y + Symbol.z -> z + else -> error("A column for symbol $symbol not found") + } +} \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/domains/DoubleDomain.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/domains/DoubleDomain.kt index 057a4a344..154763159 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/domains/DoubleDomain.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/domains/DoubleDomain.kt @@ -24,7 +24,6 @@ import space.kscience.kmath.misc.UnstableKMathAPI */ @UnstableKMathAPI public interface DoubleDomain : Domain { - /** * Global lower edge * @param num axis number diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/DifferentiableExpression.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/DifferentiableExpression.kt index 5cbc4dbf4..1f0ceaec3 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/DifferentiableExpression.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/DifferentiableExpression.kt @@ -1,5 +1,8 @@ package space.kscience.kmath.expressions +import space.kscience.kmath.misc.StringSymbol +import space.kscience.kmath.misc.Symbol + /** * Represents expression which structure can be differentiated. * diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/Expression.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/Expression.kt index 5ba24aa62..7918f199e 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/Expression.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/Expression.kt @@ -1,26 +1,11 @@ package space.kscience.kmath.expressions +import space.kscience.kmath.misc.StringSymbol +import space.kscience.kmath.misc.Symbol import space.kscience.kmath.operations.Algebra import kotlin.jvm.JvmName import kotlin.properties.ReadOnlyProperty -/** - * A marker interface for a symbol. A symbol mus have an identity - */ -public interface Symbol { - /** - * Identity object for the symbol. Two symbols with the same identity are considered to be the same symbol. - */ - public val identity: String -} - -/** - * A [Symbol] with a [String] identity - */ -public inline class StringSymbol(override val identity: String) : Symbol { - override fun toString(): String = identity -} - /** * An elementary function that could be invoked on a map of arguments. * @@ -92,13 +77,6 @@ public interface ExpressionAlgebra : Algebra { public fun ExpressionAlgebra.bindSymbol(symbol: Symbol): E = bindSymbolOrNull(symbol) ?: error("Symbol $symbol could not be bound to $this") -/** - * A delegate to create a symbol with a string identity in this scope - */ -public val symbol: ReadOnlyProperty = ReadOnlyProperty { _, property -> - StringSymbol(property.name) -} - /** * Bind a symbol by name inside the [ExpressionAlgebra] */ diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/FunctionalExpressionAlgebra.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/FunctionalExpressionAlgebra.kt index cca75754f..ebd9e7f22 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/FunctionalExpressionAlgebra.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/FunctionalExpressionAlgebra.kt @@ -1,5 +1,6 @@ package space.kscience.kmath.expressions +import space.kscience.kmath.misc.Symbol import space.kscience.kmath.operations.* /** diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/SimpleAutoDiff.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/SimpleAutoDiff.kt index 4b0d402ed..d9be4a92e 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/SimpleAutoDiff.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/SimpleAutoDiff.kt @@ -1,6 +1,7 @@ package space.kscience.kmath.expressions import space.kscience.kmath.linear.Point +import space.kscience.kmath.misc.Symbol import space.kscience.kmath.misc.UnstableKMathAPI import space.kscience.kmath.operations.* import space.kscience.kmath.structures.asBuffer diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/SymbolIndexer.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/SymbolIndexer.kt index 580acaafb..4db4b5828 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/SymbolIndexer.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/SymbolIndexer.kt @@ -1,6 +1,8 @@ package space.kscience.kmath.expressions import space.kscience.kmath.linear.Point +import space.kscience.kmath.misc.Symbol +import space.kscience.kmath.misc.UnstableKMathAPI import space.kscience.kmath.nd.Structure2D import space.kscience.kmath.structures.BufferFactory @@ -8,6 +10,7 @@ import space.kscience.kmath.structures.BufferFactory * An environment to easy transform indexed variables to symbols and back. * TODO requires multi-receivers to be beautiful */ +@UnstableKMathAPI public interface SymbolIndexer { public val symbols: List public fun indexOf(symbol: Symbol): Int = symbols.indexOf(symbol) @@ -49,13 +52,16 @@ public interface SymbolIndexer { public fun Map.toDoubleArray(): DoubleArray = DoubleArray(symbols.size) { getValue(symbols[it]) } } +@UnstableKMathAPI public inline class SimpleSymbolIndexer(override val symbols: List) : SymbolIndexer /** * Execute the block with symbol indexer based on given symbol order */ +@UnstableKMathAPI public inline fun withSymbols(vararg symbols: Symbol, block: SymbolIndexer.() -> R): R = with(SimpleSymbolIndexer(symbols.toList()), block) +@UnstableKMathAPI public inline fun withSymbols(symbols: Collection, block: SymbolIndexer.() -> R): R = with(SimpleSymbolIndexer(symbols.toList()), block) \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/LinearSolver.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/LinearSolver.kt index af136c552..3e2dbee3f 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/LinearSolver.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/LinearSolver.kt @@ -3,7 +3,10 @@ package space.kscience.kmath.linear import space.kscience.kmath.nd.as1D /** - * A group of methods to resolve equation A dot X = B, where A and B are matrices or vectors + * A group of methods to solve for *X* in equation *X = A -1 · B*, where *A* and *B* are matrices or + * vectors. + * + * @param T the type of items. */ public interface LinearSolver { /** @@ -23,7 +26,7 @@ public interface LinearSolver { } /** - * Convert matrix to vector if it is possible + * Convert matrix to vector if it is possible. */ public fun Matrix.asVector(): Point = if (this.colNum == 1) @@ -31,4 +34,11 @@ public fun Matrix.asVector(): Point = else error("Can't convert matrix with more than one column to vector") +/** + * Creates an n × 1 [VirtualMatrix], where n is the size of the given buffer. + * + * @param T the type of elements contained in the buffer. + * @receiver a buffer. + * @return the new matrix. + */ public fun Point.asMatrix(): VirtualMatrix = VirtualMatrix(size, 1) { i, _ -> get(i) } diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/LinearSpace.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/LinearSpace.kt index 4f2afc6fa..21094db73 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/LinearSpace.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/LinearSpace.kt @@ -18,6 +18,8 @@ public typealias MutableMatrix = MutableStructure2D /** * Alias or using [Buffer] as a point/vector in a many-dimensional space. + * + * @param T the type of elements contained in the buffer. */ public typealias Point = Buffer @@ -165,7 +167,7 @@ public interface LinearSpace> { * @return a feature object or `null` if it isn't present. */ @UnstableKMathAPI - public fun getFeature(structure: Matrix, type: KClass): F? = structure.getFeature(type) + public fun getFeature(structure: Matrix, type: KClass): F? = structure.getFeature(type) public companion object { @@ -195,7 +197,7 @@ public interface LinearSpace> { * @return a feature object or `null` if it isn't present. */ @UnstableKMathAPI -public inline fun LinearSpace.getFeature(structure: Matrix): F? = +public inline fun LinearSpace.getFeature(structure: Matrix): F? = getFeature(structure, F::class) diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/MatrixFeatures.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/MatrixFeatures.kt index 6b97e89ef..30e3daa7a 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/MatrixFeatures.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/MatrixFeatures.kt @@ -1,10 +1,12 @@ package space.kscience.kmath.linear +import space.kscience.kmath.nd.StructureFeature + /** * A marker interface representing some properties of matrices or additional transformations of them. Features are used * to optimize matrix operations performance in some cases or retrieve the APIs. */ -public interface MatrixFeature +public interface MatrixFeature: StructureFeature /** * Matrices with this feature are considered to have only diagonal non-null elements. diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/MatrixWrapper.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/MatrixWrapper.kt index 868f74cc6..def3b87f7 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/MatrixWrapper.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/MatrixWrapper.kt @@ -1,6 +1,7 @@ package space.kscience.kmath.linear import space.kscience.kmath.misc.UnstableKMathAPI +import space.kscience.kmath.nd.StructureFeature import space.kscience.kmath.nd.getFeature import space.kscience.kmath.operations.Ring import kotlin.reflect.KClass @@ -20,7 +21,7 @@ public class MatrixWrapper internal constructor( */ @UnstableKMathAPI @Suppress("UNCHECKED_CAST") - override fun getFeature(type: KClass): T? = features.singleOrNull { type.isInstance(it) } as? T + override fun getFeature(type: KClass): F? = features.singleOrNull { type.isInstance(it) } as? F ?: origin.getFeature(type) override fun toString(): String { diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/misc/Symbol.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/misc/Symbol.kt new file mode 100644 index 000000000..2c9774b6a --- /dev/null +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/misc/Symbol.kt @@ -0,0 +1,34 @@ +package space.kscience.kmath.misc + +import kotlin.properties.ReadOnlyProperty + +/** + * A marker interface for a symbol. A symbol mus have an identity + */ +public interface Symbol { + /** + * Identity object for the symbol. Two symbols with the same identity are considered to be the same symbol. + */ + public val identity: String + + public companion object{ + public val x: StringSymbol = StringSymbol("x") + public val y: StringSymbol = StringSymbol("y") + public val z: StringSymbol = StringSymbol("z") + } +} + +/** + * A [Symbol] with a [String] identity + */ +public inline class StringSymbol(override val identity: String) : Symbol { + override fun toString(): String = identity +} + + +/** + * A delegate to create a symbol with a string identity in this scope + */ +public val symbol: ReadOnlyProperty = ReadOnlyProperty { _, property -> + StringSymbol(property.name) +} diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/AlgebraND.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/AlgebraND.kt index b23ce947d..2821a6648 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/AlgebraND.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/AlgebraND.kt @@ -67,7 +67,8 @@ public interface AlgebraND> { * @return a feature object or `null` if it isn't present. */ @UnstableKMathAPI - public fun getFeature(structure: StructureND, type: KClass): F? = structure.getFeature(type) + public fun getFeature(structure: StructureND, type: KClass): F? = + structure.getFeature(type) public companion object } @@ -81,7 +82,7 @@ public interface AlgebraND> { * @return a feature object or `null` if it isn't present. */ @UnstableKMathAPI -public inline fun AlgebraND.getFeature(structure: StructureND): F? = +public inline fun AlgebraND.getFeature(structure: StructureND): F? = getFeature(structure, F::class) /** diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/Structure1D.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/Structure1D.kt index 354f3d802..313e8fb31 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/Structure1D.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/Structure1D.kt @@ -17,6 +17,8 @@ public interface Structure1D : StructureND, Buffer { } public override operator fun iterator(): Iterator = (0 until size).asSequence().map(::get).iterator() + + public companion object } /** diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/Structure2D.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/Structure2D.kt index e1a1d37de..5d83ce3dd 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/Structure2D.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/Structure2D.kt @@ -97,7 +97,7 @@ private inline class Structure2DWrapper(val structure: StructureND) : Stru override operator fun get(i: Int, j: Int): T = structure[i, j] @UnstableKMathAPI - override fun getFeature(type: KClass): F? = structure.getFeature(type) + override fun getFeature(type: KClass): F? = structure.getFeature(type) override fun elements(): Sequence> = structure.elements() } diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/StructureND.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/StructureND.kt index 78eac1809..a1aa5e554 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/StructureND.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/StructureND.kt @@ -7,6 +7,8 @@ import kotlin.jvm.JvmName import kotlin.native.concurrent.ThreadLocal import kotlin.reflect.KClass +public interface StructureFeature + /** * Represents n-dimensional structure, i.e. multidimensional container of items of the same type and size. The number * of dimensions and items in an array is defined by its shape, which is a sequence of non-negative integers that @@ -48,7 +50,7 @@ public interface StructureND { * If the feature is not present, null is returned. */ @UnstableKMathAPI - public fun getFeature(type: KClass): F? = null + public fun getFeature(type: KClass): F? = null public companion object { /** @@ -144,7 +146,7 @@ public interface StructureND { public operator fun StructureND.get(vararg index: Int): T = get(index) @UnstableKMathAPI -public inline fun StructureND<*>.getFeature(): T? = getFeature(T::class) +public inline fun StructureND<*>.getFeature(): T? = getFeature(T::class) /** * Represents mutable [StructureND]. diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/Algebra.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/Algebra.kt index 9f57bc4c1..492ec8e88 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/Algebra.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/Algebra.kt @@ -1,6 +1,6 @@ package space.kscience.kmath.operations -import space.kscience.kmath.expressions.Symbol +import space.kscience.kmath.misc.Symbol /** * Stub for DSL the [Algebra] is. @@ -100,8 +100,8 @@ public fun Algebra.bindSymbol(symbol: Symbol): T = bindSymbol(symbo public inline operator fun , R> A.invoke(block: A.() -> R): R = run(block) /** - * Represents linear space without neutral element, i.e. algebraic structure with associative, binary operation [add] - * and scalar multiplication [multiply]. + * Represents group without neutral element (also known as inverse semigroup), i.e. algebraic structure with + * associative, binary operation [add]. * * @param T the type of element of this semispace. */ @@ -177,7 +177,7 @@ public interface GroupOperations : Algebra { } /** - * Represents linear space with neutral element, i.e. algebraic structure with associative, binary operation [add]. + * Represents group, i.e. algebraic structure with associative, binary operation [add]. * * @param T the type of element of this semispace. */ @@ -189,8 +189,8 @@ public interface Group : GroupOperations { } /** - * Represents rng, i.e. algebraic structure with associative, binary, commutative operation [add] and associative, - * operation [multiply] distributive over [add]. + * Represents ring without multiplicative and additive identities, i.e. algebraic structure with + * associative, binary, commutative operation [add] and associative, operation [multiply] distributive over [add]. * * @param T the type of element of this semiring. */ @@ -238,7 +238,7 @@ public interface Ring : Group, RingOperations { } /** - * Represents field without identity elements, i.e. algebraic structure with associative, binary, commutative operations + * Represents field without without multiplicative and additive identities, i.e. algebraic structure with associative, binary, commutative operations * [add] and [multiply]; binary operation [divide] as multiplication of left operand by reciprocal of right one. * * @param T the type of element of this semifield. @@ -276,10 +276,11 @@ public interface FieldOperations : RingOperations { } /** - * Represents field, i.e. algebraic structure with three operations: associative "addition" and "multiplication", - * and "division" and their neutral elements. + * Represents field, i.e. algebraic structure with three operations: associative, commutative addition and + * multiplication, and division. **This interface differs from the eponymous mathematical definition: fields in KMath + * also support associative multiplication by scalar.** * - * @param T the type of element of this semifield. + * @param T the type of element of this field. */ public interface Field : Ring, FieldOperations, ScaleOperations, NumericAlgebra { override fun number(value: Number): T = scale(one, value.toDouble()) diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/BigInt.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/BigInt.kt index 55bb68850..18fbf0fdd 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/BigInt.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/BigInt.kt @@ -12,8 +12,8 @@ import kotlin.math.max import kotlin.math.min import kotlin.math.sign -public typealias Magnitude = UIntArray -public typealias TBase = ULong +private typealias Magnitude = UIntArray +private typealias TBase = ULong /** * Kotlin Multiplatform implementation of Big Integer numbers (KBigInteger). @@ -241,18 +241,18 @@ public class BigInt internal constructor( ) private fun compareMagnitudes(mag1: Magnitude, mag2: Magnitude): Int { - when { - mag1.size > mag2.size -> return 1 - mag1.size < mag2.size -> return -1 + return when { + mag1.size > mag2.size -> 1 + mag1.size < mag2.size -> -1 + else -> { - for (i in mag1.size - 1 downTo 0) { - if (mag1[i] > mag2[i]) { - return 1 - } else if (mag1[i] < mag2[i]) { - return -1 - } + for (i in mag1.size - 1 downTo 0) return when { + mag1[i] > mag2[i] -> 1 + mag1[i] < mag2[i] -> -1 + else -> continue } - return 0 + + 0 } } } @@ -302,10 +302,11 @@ public class BigInt internal constructor( var carry = 0uL for (i in mag.indices) { - val cur: ULong = carry + mag[i].toULong() * x.toULong() + val cur = carry + mag[i].toULong() * x.toULong() result[i] = (cur and BASE).toUInt() carry = cur shr BASE_SIZE } + result[resultLength - 1] = (carry and BASE).toUInt() return stripLeadingZeros(result) @@ -358,6 +359,9 @@ private fun stripLeadingZeros(mag: Magnitude): Magnitude { return mag.sliceArray(IntRange(0, resSize)) } +/** + * Returns the absolute value of the given value [x]. + */ public fun abs(x: BigInt): BigInt = x.abs() /** diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/Buffer.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/Buffer.kt index 72a6e3cae..1b4ccb069 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/Buffer.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/Buffer.kt @@ -40,7 +40,6 @@ public interface Buffer { public operator fun iterator(): Iterator public companion object { - /** * Check the element-by-element match of content of two buffers. */ @@ -110,7 +109,6 @@ public interface MutableBuffer : Buffer { public fun copy(): MutableBuffer public companion object { - /** * Creates a [DoubleBuffer] with the specified [size], where each element is calculated by calling the specified * [initializer] function. diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/bufferOperation.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/bufferOperation.kt index 4355ba17f..5f8bbe21f 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/bufferOperation.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/bufferOperation.kt @@ -70,6 +70,15 @@ public inline fun Buffer.mapIndexed( crossinline block: (index: Int, value: T) -> R, ): Buffer = bufferFactory(size) { block(it, get(it)) } +/** + * Fold given buffer according to [operation] + */ +public inline fun Buffer.fold(initial: R, operation: (acc: R, T) -> R): R { + var accumulator = initial + for (index in this.indices) accumulator = operation(accumulator, get(index)) + return accumulator +} + /** * Zip two buffers using given [transform]. */ diff --git a/kmath-core/src/commonTest/kotlin/space/kscience/kmath/expressions/ExpressionFieldTest.kt b/kmath-core/src/commonTest/kotlin/space/kscience/kmath/expressions/ExpressionFieldTest.kt index b1be0c392..61b1d03f0 100644 --- a/kmath-core/src/commonTest/kotlin/space/kscience/kmath/expressions/ExpressionFieldTest.kt +++ b/kmath-core/src/commonTest/kotlin/space/kscience/kmath/expressions/ExpressionFieldTest.kt @@ -1,5 +1,6 @@ package space.kscience.kmath.expressions +import space.kscience.kmath.misc.symbol import space.kscience.kmath.operations.DoubleField import space.kscience.kmath.operations.invoke import kotlin.test.Test diff --git a/kmath-core/src/commonTest/kotlin/space/kscience/kmath/expressions/SimpleAutoDiffTest.kt b/kmath-core/src/commonTest/kotlin/space/kscience/kmath/expressions/SimpleAutoDiffTest.kt index ca7aca905..666db13d8 100644 --- a/kmath-core/src/commonTest/kotlin/space/kscience/kmath/expressions/SimpleAutoDiffTest.kt +++ b/kmath-core/src/commonTest/kotlin/space/kscience/kmath/expressions/SimpleAutoDiffTest.kt @@ -1,5 +1,7 @@ package space.kscience.kmath.expressions +import space.kscience.kmath.misc.Symbol +import space.kscience.kmath.misc.symbol import space.kscience.kmath.operations.DoubleField import space.kscience.kmath.structures.Buffer import space.kscience.kmath.structures.asBuffer diff --git a/kmath-core/src/commonTest/kotlin/space/kscience/kmath/structures/NumberNDFieldTest.kt b/kmath-core/src/commonTest/kotlin/space/kscience/kmath/structures/NumberNDFieldTest.kt index fb67f0308..376415a56 100644 --- a/kmath-core/src/commonTest/kotlin/space/kscience/kmath/structures/NumberNDFieldTest.kt +++ b/kmath-core/src/commonTest/kotlin/space/kscience/kmath/structures/NumberNDFieldTest.kt @@ -38,12 +38,11 @@ class NumberNDFieldTest { (i * 10 + j).toDouble() } - for (i in 0..2) { + for (i in 0..2) for (j in 0..2) { val expected = (i * 10 + j).toDouble() assertEquals(expected, array[i, j], "Error at index [$i, $j]") } - } } @Test diff --git a/kmath-coroutines/src/commonMain/kotlin/space/kscience/kmath/chains/BlockingDoubleChain.kt b/kmath-coroutines/src/commonMain/kotlin/space/kscience/kmath/chains/BlockingDoubleChain.kt index bc4508994..d024147b4 100644 --- a/kmath-coroutines/src/commonMain/kotlin/space/kscience/kmath/chains/BlockingDoubleChain.kt +++ b/kmath-coroutines/src/commonMain/kotlin/space/kscience/kmath/chains/BlockingDoubleChain.kt @@ -1,12 +1,13 @@ package space.kscience.kmath.chains /** - * Performance optimized chain for real values + * Chunked, specialized chain for real values. */ -public abstract class BlockingDoubleChain : Chain { - public abstract fun nextDouble(): Double +public interface BlockingDoubleChain : Chain { + public override suspend fun next(): Double - override suspend fun next(): Double = nextDouble() - - public fun nextBlock(size: Int): DoubleArray = DoubleArray(size) { nextDouble() } + /** + * Returns an [DoubleArray] chunk of [size] values of [next]. + */ + public suspend fun nextBlock(size: Int): DoubleArray = DoubleArray(size) { next() } } diff --git a/kmath-coroutines/src/commonMain/kotlin/space/kscience/kmath/chains/BlockingIntChain.kt b/kmath-coroutines/src/commonMain/kotlin/space/kscience/kmath/chains/BlockingIntChain.kt index 11da7e503..fb2e453ad 100644 --- a/kmath-coroutines/src/commonMain/kotlin/space/kscience/kmath/chains/BlockingIntChain.kt +++ b/kmath-coroutines/src/commonMain/kotlin/space/kscience/kmath/chains/BlockingIntChain.kt @@ -3,10 +3,7 @@ package space.kscience.kmath.chains /** * Performance optimized chain for integer values */ -public abstract class BlockingIntChain : Chain { - public abstract fun nextInt(): Int - - override suspend fun next(): Int = nextInt() - - public fun nextBlock(size: Int): IntArray = IntArray(size) { nextInt() } +public interface BlockingIntChain : Chain { + public override suspend fun next(): Int + public suspend fun nextBlock(size: Int): IntArray = IntArray(size) { next() } } diff --git a/kmath-coroutines/src/commonMain/kotlin/space/kscience/kmath/chains/Chain.kt b/kmath-coroutines/src/commonMain/kotlin/space/kscience/kmath/chains/Chain.kt index 26d078fcb..a961f2e09 100644 --- a/kmath-coroutines/src/commonMain/kotlin/space/kscience/kmath/chains/Chain.kt +++ b/kmath-coroutines/src/commonMain/kotlin/space/kscience/kmath/chains/Chain.kt @@ -63,12 +63,10 @@ public class MarkovChain(private val seed: suspend () -> R, private public fun value(): R? = value - public override suspend fun next(): R { - mutex.withLock { - val newValue = gen(value ?: seed()) - value = newValue - return newValue - } + public override suspend fun next(): R = mutex.withLock { + val newValue = gen(value ?: seed()) + value = newValue + newValue } public override fun fork(): Chain = MarkovChain(seed = { value ?: seed() }, gen = gen) @@ -90,12 +88,10 @@ public class StatefulChain( public fun value(): R? = value - public override suspend fun next(): R { - mutex.withLock { - val newValue = state.gen(value ?: state.seed()) - value = newValue - return newValue - } + public override suspend fun next(): R = mutex.withLock { + val newValue = state.gen(value ?: state.seed()) + value = newValue + newValue } public override fun fork(): Chain = StatefulChain(forkState(state), seed, forkState, gen) diff --git a/kmath-coroutines/src/commonMain/kotlin/space/kscience/kmath/streaming/BufferFlow.kt b/kmath-coroutines/src/commonMain/kotlin/space/kscience/kmath/streaming/BufferFlow.kt index c271f8889..dc1dd4757 100644 --- a/kmath-coroutines/src/commonMain/kotlin/space/kscience/kmath/streaming/BufferFlow.kt +++ b/kmath-coroutines/src/commonMain/kotlin/space/kscience/kmath/streaming/BufferFlow.kt @@ -28,7 +28,7 @@ public fun Flow.chunked(bufferSize: Int, bufferFactory: BufferFactory) var counter = 0 this@chunked.collect { element -> - list.add(element) + list += element counter++ if (counter == bufferSize) { diff --git a/kmath-coroutines/src/commonMain/kotlin/space/kscience/kmath/streaming/RingBuffer.kt b/kmath-coroutines/src/commonMain/kotlin/space/kscience/kmath/streaming/RingBuffer.kt index f81ad2f0d..e844b765d 100644 --- a/kmath-coroutines/src/commonMain/kotlin/space/kscience/kmath/streaming/RingBuffer.kt +++ b/kmath-coroutines/src/commonMain/kotlin/space/kscience/kmath/streaming/RingBuffer.kt @@ -48,11 +48,9 @@ public class RingBuffer( /** * A safe snapshot operation */ - public suspend fun snapshot(): Buffer { - mutex.withLock { - val copy = buffer.copy() - return VirtualBuffer(size) { i -> copy[startIndex.forward(i)] as T } - } + public suspend fun snapshot(): Buffer = mutex.withLock { + val copy = buffer.copy() + VirtualBuffer(size) { i -> copy[startIndex.forward(i)] as T } } public suspend fun push(element: T) { diff --git a/kmath-dimensions/src/commonTest/kotlin/kscience/dimensions/DMatrixContextTest.kt b/kmath-dimensions/src/commonTest/kotlin/kscience/dimensions/DMatrixContextTest.kt index e2a9628ac..58ed82723 100644 --- a/kmath-dimensions/src/commonTest/kotlin/kscience/dimensions/DMatrixContextTest.kt +++ b/kmath-dimensions/src/commonTest/kotlin/kscience/dimensions/DMatrixContextTest.kt @@ -1,4 +1,4 @@ -package kscience.dimensions +package space.kscience.dimensions import space.kscience.kmath.dimensions.D2 import space.kscience.kmath.dimensions.D3 diff --git a/kmath-ejml/README.md b/kmath-ejml/README.md index 1081b2b7f..2551703a4 100644 --- a/kmath-ejml/README.md +++ b/kmath-ejml/README.md @@ -1,43 +1,37 @@ -# ejml-simple support (`kmath-ejml`) +# Module kmath-ejml -This subproject implements the following features: +EJML based linear algebra implementation. - [ejml-vector](src/main/kotlin/space/kscience/kmath/ejml/EjmlVector.kt) : The Point implementation using SimpleMatrix. - [ejml-matrix](src/main/kotlin/space/kscience/kmath/ejml/EjmlMatrix.kt) : The Matrix implementation using SimpleMatrix. - [ejml-linear-space](src/main/kotlin/space/kscience/kmath/ejml/EjmlLinearSpace.kt) : The LinearSpace implementation using SimpleMatrix. -> #### Artifact: -> -> This module artifact: `space.kscience:kmath-ejml:0.3.0-dev-3`. -> -> Bintray release version: [ ![Download](https://api.bintray.com/packages/mipt-npm/kscience/kmath-ejml/images/download.svg) ](https://bintray.com/mipt-npm/kscience/kmath-ejml/_latestVersion) -> -> Bintray development version: [ ![Download](https://api.bintray.com/packages/mipt-npm/dev/kmath-ejml/images/download.svg) ](https://bintray.com/mipt-npm/dev/kmath-ejml/_latestVersion) -> -> **Gradle:** -> -> ```gradle -> repositories { -> maven { url 'https://repo.kotlin.link' } -> maven { url 'https://dl.bintray.com/hotkeytlt/maven' } -> maven { url "https://dl.bintray.com/kotlin/kotlin-eap" } // include for builds based on kotlin-eap -> } -> -> dependencies { -> implementation 'space.kscience:kmath-ejml:0.3.0-dev-3' -> } -> ``` -> **Gradle Kotlin DSL:** -> -> ```kotlin -> repositories { -> maven("https://repo.kotlin.link") -> maven("https://dl.bintray.com/kotlin/kotlin-eap") // include for builds based on kotlin-eap -> maven("https://dl.bintray.com/hotkeytlt/maven") // required for a -> } -> -> dependencies { -> implementation("space.kscience:kmath-ejml:0.3.0-dev-3") -> } -> ``` +## Artifact: + +The Maven coordinates of this project are `space.kscience:kmath-ejml:0.3.0-dev-3`. + +**Gradle:** +```gradle +repositories { + maven { url 'https://repo.kotlin.link' } + maven { url 'https://dl.bintray.com/hotkeytlt/maven' } + maven { url "https://dl.bintray.com/kotlin/kotlin-eap" } // include for builds based on kotlin-eap +} + +dependencies { + implementation 'space.kscience:kmath-ejml:0.3.0-dev-3' +} +``` +**Gradle Kotlin DSL:** +```kotlin +repositories { + maven("https://repo.kotlin.link") + maven("https://dl.bintray.com/kotlin/kotlin-eap") // include for builds based on kotlin-eap + maven("https://dl.bintray.com/hotkeytlt/maven") // required for a +} + +dependencies { + implementation("space.kscience:kmath-ejml:0.3.0-dev-3") +} +``` diff --git a/kmath-ejml/docs/README-TEMPLATE.md b/kmath-ejml/docs/README-TEMPLATE.md index c53f4a81c..27fcedd65 100644 --- a/kmath-ejml/docs/README-TEMPLATE.md +++ b/kmath-ejml/docs/README-TEMPLATE.md @@ -1,6 +1,6 @@ -# ejml-simple support (`kmath-ejml`) +# Module kmath-ejml -This subproject implements the following features: +EJML based linear algebra implementation. ${features} diff --git a/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/EjmlLinearSpace.kt b/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/EjmlLinearSpace.kt index 6fc0a049c..94429e7b3 100644 --- a/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/EjmlLinearSpace.kt +++ b/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/EjmlLinearSpace.kt @@ -4,6 +4,7 @@ import org.ejml.dense.row.factory.DecompositionFactory_DDRM import org.ejml.simple.SimpleMatrix import space.kscience.kmath.linear.* import space.kscience.kmath.misc.UnstableKMathAPI +import space.kscience.kmath.nd.StructureFeature import space.kscience.kmath.nd.getFeature import space.kscience.kmath.operations.DoubleField import space.kscience.kmath.structures.DoubleBuffer @@ -95,7 +96,7 @@ public object EjmlLinearSpace : LinearSpace { v.toEjml().origin.scale(this).wrapVector() @UnstableKMathAPI - public override fun getFeature(structure: Matrix, type: KClass): F? { + public override fun getFeature(structure: Matrix, type: KClass): F? { //Return the feature if it is intrinsic to the structure structure.getFeature(type)?.let { return it } diff --git a/kmath-for-real/README.md b/kmath-for-real/README.md index 139cd3cef..ad3d33062 100644 --- a/kmath-for-real/README.md +++ b/kmath-for-real/README.md @@ -1,41 +1,37 @@ -# Real number specialization module (`kmath-for-real`) +# Module kmath-for-real + +Specialization of KMath APIs for Double numbers. - [DoubleVector](src/commonMain/kotlin/space/kscience/kmath/real/DoubleVector.kt) : Numpy-like operations for Buffers/Points - [DoubleMatrix](src/commonMain/kotlin/space/kscience/kmath/real/DoubleMatrix.kt) : Numpy-like operations for 2d real structures - [grids](src/commonMain/kotlin/space/kscience/kmath/structures/grids.kt) : Uniform grid generators -> #### Artifact: -> -> This module artifact: `space.kscience:kmath-for-real:0.3.0-dev-3`. -> -> Bintray release version: [ ![Download](https://api.bintray.com/packages/mipt-npm/kscience/kmath-for-real/images/download.svg) ](https://bintray.com/mipt-npm/kscience/kmath-for-real/_latestVersion) -> -> Bintray development version: [ ![Download](https://api.bintray.com/packages/mipt-npm/dev/kmath-for-real/images/download.svg) ](https://bintray.com/mipt-npm/dev/kmath-for-real/_latestVersion) -> -> **Gradle:** -> -> ```gradle -> repositories { -> maven { url 'https://repo.kotlin.link' } -> maven { url 'https://dl.bintray.com/hotkeytlt/maven' } -> maven { url "https://dl.bintray.com/kotlin/kotlin-eap" } // include for builds based on kotlin-eap -> } -> -> dependencies { -> implementation 'space.kscience:kmath-for-real:0.3.0-dev-3' -> } -> ``` -> **Gradle Kotlin DSL:** -> -> ```kotlin -> repositories { -> maven("https://repo.kotlin.link") -> maven("https://dl.bintray.com/kotlin/kotlin-eap") // include for builds based on kotlin-eap -> maven("https://dl.bintray.com/hotkeytlt/maven") // required for a -> } -> -> dependencies { -> implementation("space.kscience:kmath-for-real:0.3.0-dev-3") -> } -> ``` +## Artifact: + +The Maven coordinates of this project are `space.kscience:kmath-for-real:0.3.0-dev-3`. + +**Gradle:** +```gradle +repositories { + maven { url 'https://repo.kotlin.link' } + maven { url 'https://dl.bintray.com/hotkeytlt/maven' } + maven { url "https://dl.bintray.com/kotlin/kotlin-eap" } // include for builds based on kotlin-eap +} + +dependencies { + implementation 'space.kscience:kmath-for-real:0.3.0-dev-3' +} +``` +**Gradle Kotlin DSL:** +```kotlin +repositories { + maven("https://repo.kotlin.link") + maven("https://dl.bintray.com/kotlin/kotlin-eap") // include for builds based on kotlin-eap + maven("https://dl.bintray.com/hotkeytlt/maven") // required for a +} + +dependencies { + implementation("space.kscience:kmath-for-real:0.3.0-dev-3") +} +``` diff --git a/kmath-for-real/docs/README-TEMPLATE.md b/kmath-for-real/docs/README-TEMPLATE.md index 670844bd0..c2ef25aa7 100644 --- a/kmath-for-real/docs/README-TEMPLATE.md +++ b/kmath-for-real/docs/README-TEMPLATE.md @@ -1,4 +1,6 @@ -# Real number specialization module (`kmath-for-real`) +# Module kmath-for-real + +Specialization of KMath APIs for Double numbers. ${features} diff --git a/kmath-for-real/src/commonMain/kotlin/space/kscience/kmath/real/RealVector.kt b/kmath-for-real/src/commonMain/kotlin/space/kscience/kmath/real/RealVector.kt index 433e3c756..6b059ef56 100644 --- a/kmath-for-real/src/commonMain/kotlin/space/kscience/kmath/real/RealVector.kt +++ b/kmath-for-real/src/commonMain/kotlin/space/kscience/kmath/real/RealVector.kt @@ -6,17 +6,13 @@ import space.kscience.kmath.operations.Norm import space.kscience.kmath.structures.Buffer import space.kscience.kmath.structures.MutableBuffer.Companion.double import space.kscience.kmath.structures.asBuffer -import space.kscience.kmath.structures.asIterable +import space.kscience.kmath.structures.fold import space.kscience.kmath.structures.indices import kotlin.math.pow import kotlin.math.sqrt public typealias DoubleVector = Point -public object VectorL2Norm : Norm, Double> { - override fun norm(arg: Point): Double = sqrt(arg.asIterable().sumByDouble(Number::toDouble)) -} - @Suppress("FunctionName") public fun DoubleVector(vararg doubles: Double): DoubleVector = doubles.asBuffer() @@ -102,4 +98,10 @@ public fun DoubleVector.sum(): Double { res += get(i) } return res -} \ No newline at end of file +} + +public object VectorL2Norm : Norm { + override fun norm(arg: DoubleVector): Double = sqrt(arg.fold(0.0) { acc: Double, d: Double -> acc + d.pow(2) }) +} + +public val DoubleVector.norm: Double get() = VectorL2Norm.norm(this) \ No newline at end of file diff --git a/kmath-for-real/src/commonMain/kotlin/space/kscience/kmath/real/grids.kt b/kmath-for-real/src/commonMain/kotlin/space/kscience/kmath/real/grids.kt index 35297a3ac..7f02c7db3 100644 --- a/kmath-for-real/src/commonMain/kotlin/space/kscience/kmath/real/grids.kt +++ b/kmath-for-real/src/commonMain/kotlin/space/kscience/kmath/real/grids.kt @@ -1,7 +1,41 @@ package space.kscience.kmath.real -import space.kscience.kmath.structures.asBuffer -import kotlin.math.abs +import space.kscience.kmath.misc.UnstableKMathAPI +import space.kscience.kmath.structures.Buffer +import space.kscience.kmath.structures.DoubleBuffer +import kotlin.math.floor + +public val ClosedFloatingPointRange.length: Double get() = endInclusive - start + +/** + * Create a Buffer-based grid with equally distributed [numberOfPoints] points. The range could be increasing or decreasing. + * If range has a zero size, then the buffer consisting of [numberOfPoints] equal values is returned. + */ +public fun Buffer.Companion.fromRange(range: ClosedFloatingPointRange, numberOfPoints: Int): DoubleBuffer { + require(numberOfPoints >= 2) { "Number of points in grid must be more than 1" } + val normalizedRange = when { + range.endInclusive > range.start -> range + range.endInclusive < range.start -> range.endInclusive..range.start + else -> return DoubleBuffer(numberOfPoints) { range.start } + } + val step = normalizedRange.length / (numberOfPoints - 1) + return DoubleBuffer(numberOfPoints) { normalizedRange.start + step * it } +} + +/** + * Create a Buffer-based grid with equally distributed points with a fixed [step]. The range could be increasing or decreasing. + * If the step is larger than the range size, single point is returned. + */ +public fun Buffer.Companion.withFixedStep(range: ClosedFloatingPointRange, step: Double): DoubleBuffer { + require(step > 0) { "The grid step must be positive" } + val normalizedRange = when { + range.endInclusive > range.start -> range + range.endInclusive < range.start -> range.endInclusive..range.start + else -> return DoubleBuffer(range.start) + } + val numberOfPoints = floor(normalizedRange.length / step).toInt() + 1 + return DoubleBuffer(numberOfPoints) { normalizedRange.start + step * it } +} /** * Convert double range to sequence. @@ -11,35 +45,5 @@ import kotlin.math.abs * * If step is negative, the same goes from upper boundary downwards */ -public fun ClosedFloatingPointRange.toSequenceWithStep(step: Double): Sequence = when { - step == 0.0 -> error("Zero step in double progression") - - step > 0 -> sequence { - var current = start - - while (current <= endInclusive) { - yield(current) - current += step - } - } - - else -> sequence { - var current = endInclusive - - while (current >= start) { - yield(current) - current += step - } - } -} - -public infix fun ClosedFloatingPointRange.step(step: Double): DoubleVector = - toSequenceWithStep(step).toList().asBuffer() - -/** - * Convert double range to sequence with the fixed number of points - */ -public fun ClosedFloatingPointRange.toSequenceWithPoints(numPoints: Int): Sequence { - require(numPoints > 1) { "The number of points should be more than 2" } - return toSequenceWithStep(abs(endInclusive - start) / (numPoints - 1)) -} +@UnstableKMathAPI +public infix fun ClosedFloatingPointRange.step(step: Double): DoubleBuffer = Buffer.withFixedStep(this, step) \ No newline at end of file diff --git a/kmath-for-real/src/commonTest/kotlin/kaceince/kmath/real/GridTest.kt b/kmath-for-real/src/commonTest/kotlin/kaceince/kmath/real/GridTest.kt index 5a644c8f9..91ee517ab 100644 --- a/kmath-for-real/src/commonTest/kotlin/kaceince/kmath/real/GridTest.kt +++ b/kmath-for-real/src/commonTest/kotlin/kaceince/kmath/real/GridTest.kt @@ -1,13 +1,18 @@ package kaceince.kmath.real +import space.kscience.kmath.real.DoubleVector +import space.kscience.kmath.real.minus +import space.kscience.kmath.real.norm import space.kscience.kmath.real.step import kotlin.test.Test import kotlin.test.assertEquals +import kotlin.test.assertTrue class GridTest { @Test fun testStepGrid() { val grid = 0.0..1.0 step 0.2 assertEquals(6, grid.size) + assertTrue { (grid - DoubleVector(0.0, 0.2, 0.4, 0.6, 0.8, 1.0)).norm < 1e-4 } } } \ No newline at end of file diff --git a/kmath-functions/README.md b/kmath-functions/README.md index d13c4c107..531e97a44 100644 --- a/kmath-functions/README.md +++ b/kmath-functions/README.md @@ -1,6 +1,6 @@ -# Functions (`kmath-functions`) +# Module kmath-functions -Functions and interpolations: +Functions and interpolations. - [piecewise](Piecewise functions.) : src/commonMain/kotlin/space/kscience/kmath/functions/Piecewise.kt - [polynomials](Polynomial functions.) : src/commonMain/kotlin/space/kscience/kmath/functions/Polynomial.kt @@ -8,37 +8,31 @@ Functions and interpolations: - [spline interpolation](Cubic spline XY interpolator.) : src/commonMain/kotlin/space/kscience/kmath/interpolation/SplineInterpolator.kt -> #### Artifact: -> -> This module artifact: `space.kscience:kmath-functions:0.3.0-dev-3`. -> -> Bintray release version: [ ![Download](https://api.bintray.com/packages/mipt-npm/kscience/kmath-functions/images/download.svg) ](https://bintray.com/mipt-npm/kscience/kmath-functions/_latestVersion) -> -> Bintray development version: [ ![Download](https://api.bintray.com/packages/mipt-npm/dev/kmath-functions/images/download.svg) ](https://bintray.com/mipt-npm/dev/kmath-functions/_latestVersion) -> -> **Gradle:** -> -> ```gradle -> repositories { -> maven { url 'https://repo.kotlin.link' } -> maven { url 'https://dl.bintray.com/hotkeytlt/maven' } -> maven { url "https://dl.bintray.com/kotlin/kotlin-eap" } // include for builds based on kotlin-eap -> } -> -> dependencies { -> implementation 'space.kscience:kmath-functions:0.3.0-dev-3' -> } -> ``` -> **Gradle Kotlin DSL:** -> -> ```kotlin -> repositories { -> maven("https://repo.kotlin.link") -> maven("https://dl.bintray.com/kotlin/kotlin-eap") // include for builds based on kotlin-eap -> maven("https://dl.bintray.com/hotkeytlt/maven") // required for a -> } -> -> dependencies { -> implementation("space.kscience:kmath-functions:0.3.0-dev-3") -> } -> ``` +## Artifact: + +The Maven coordinates of this project are `space.kscience:kmath-functions:0.3.0-dev-3`. + +**Gradle:** +```gradle +repositories { + maven { url 'https://repo.kotlin.link' } + maven { url 'https://dl.bintray.com/hotkeytlt/maven' } + maven { url "https://dl.bintray.com/kotlin/kotlin-eap" } // include for builds based on kotlin-eap +} + +dependencies { + implementation 'space.kscience:kmath-functions:0.3.0-dev-3' +} +``` +**Gradle Kotlin DSL:** +```kotlin +repositories { + maven("https://repo.kotlin.link") + maven("https://dl.bintray.com/kotlin/kotlin-eap") // include for builds based on kotlin-eap + maven("https://dl.bintray.com/hotkeytlt/maven") // required for a +} + +dependencies { + implementation("space.kscience:kmath-functions:0.3.0-dev-3") +} +``` diff --git a/kmath-functions/docs/README-TEMPLATE.md b/kmath-functions/docs/README-TEMPLATE.md index 8a34a7cc4..2e163eee5 100644 --- a/kmath-functions/docs/README-TEMPLATE.md +++ b/kmath-functions/docs/README-TEMPLATE.md @@ -1,6 +1,6 @@ -# Functions (`kmath-functions`) +# Module kmath-functions -Functions and interpolations: +Functions and interpolations. ${features} diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Piecewise.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Piecewise.kt index d2470a4b4..2477af618 100644 --- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Piecewise.kt +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Piecewise.kt @@ -2,14 +2,28 @@ package space.kscience.kmath.functions import space.kscience.kmath.operations.Ring +/** + * Represents piecewise-defined function. + * + * @param T the piece key type. + * @param R the sub-function type. + */ public fun interface Piecewise { + /** + * Returns the appropriate sub-function for given piece key. + */ public fun findPiece(arg: T): R? } +/** + * Represents piecewise-defined function where all the sub-functions are polynomials. + */ public fun interface PiecewisePolynomial : Piecewise> /** - * Ordered list of pieces in piecewise function + * Basic [Piecewise] implementation where all the pieces are ordered by the [Comparable] type instances. + * + * @param T the comparable piece key type. */ public class OrderedPiecewisePolynomial>(delimiter: T) : PiecewisePolynomial { @@ -17,22 +31,30 @@ public class OrderedPiecewisePolynomial>(delimiter: T) : private val pieces: MutableList> = arrayListOf() /** - * Dynamically add a piece to the "right" side (beyond maximum argument value of previous piece) - * @param right new rightmost position. If is less then current rightmost position, a error is thrown. + * Dynamically adds a piece to the right side (beyond maximum argument value of previous piece) + * + * @param right new rightmost position. If is less then current rightmost position, an error is thrown. + * @param piece the sub-function. */ public fun putRight(right: T, piece: Polynomial) { require(right > delimiters.last()) { "New delimiter should be to the right of old one" } - delimiters.add(right) - pieces.add(piece) + delimiters += right + pieces += piece } + /** + * Dynamically adds a piece to the left side (beyond maximum argument value of previous piece) + * + * @param left the new leftmost position. If is less then current rightmost position, an error is thrown. + * @param piece the sub-function. + */ public fun putLeft(left: T, piece: Polynomial) { require(left < delimiters.first()) { "New delimiter should be to the left of old one" } delimiters.add(0, left) pieces.add(0, piece) } - override fun findPiece(arg: T): Polynomial? { + public override fun findPiece(arg: T): Polynomial? { if (arg < delimiters.first() || arg >= delimiters.last()) return null else { @@ -46,9 +68,10 @@ public class OrderedPiecewisePolynomial>(delimiter: T) : } /** - * Return a value of polynomial function with given [ring] an given [arg] or null if argument is outside of piecewise definition. + * Return a value of polynomial function with given [ring] an given [arg] or null if argument is outside of piecewise + * definition. */ public fun , C : Ring> PiecewisePolynomial.value(ring: C, arg: T): T? = findPiece(arg)?.value(ring, arg) -public fun , C : Ring> PiecewisePolynomial.asFunction(ring: C): (T) -> T? = { value(ring, it) } \ No newline at end of file +public fun , C : Ring> PiecewisePolynomial.asFunction(ring: C): (T) -> T? = { value(ring, it) } diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Polynomial.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Polynomial.kt index 049909fe1..550785812 100644 --- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Polynomial.kt +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Polynomial.kt @@ -10,16 +10,26 @@ import kotlin.math.max import kotlin.math.pow /** - * Polynomial coefficients without fixation on specific context they are applied to - * @param coefficients constant is the leftmost coefficient + * Polynomial coefficients model without fixation on specific context they are applied to. + * + * @param coefficients constant is the leftmost coefficient. */ public inline class Polynomial(public val coefficients: List) +/** + * Returns a [Polynomial] instance with given [coefficients]. + */ @Suppress("FunctionName") public fun Polynomial(vararg coefficients: T): Polynomial = Polynomial(coefficients.toList()) +/** + * Evaluates the value of the given double polynomial for given double argument. + */ public fun Polynomial.value(): Double = coefficients.reduceIndexed { index, acc, d -> acc + d.pow(index) } +/** + * Evaluates the value of the given polynomial for given argument. + */ public fun > Polynomial.value(ring: C, arg: T): T = ring { if (coefficients.isEmpty()) return@ring zero var res = coefficients.first() @@ -35,19 +45,23 @@ public fun > Polynomial.value(ring: C, arg: T): T = ring } /** - * Represent the polynomial as a regular context-less function + * Represent the polynomial as a regular context-less function. */ public fun > Polynomial.asFunction(ring: C): (T) -> T = { value(ring, it) } /** - * An algebra for polynomials + * Space of polynomials. + * + * @param T the type of operated polynomials. + * @param C the intersection of [Ring] of [T] and [ScaleOperations] of [T]. + * @param ring the [C] instance. */ public class PolynomialSpace( private val ring: C, ) : Group>, ScaleOperations> where C : Ring, C : ScaleOperations { public override val zero: Polynomial = Polynomial(emptyList()) - override fun Polynomial.unaryMinus(): Polynomial = with(ring) { + override fun Polynomial.unaryMinus(): Polynomial = ring { Polynomial(coefficients.map { -it }) } @@ -64,6 +78,9 @@ public class PolynomialSpace( public override fun scale(a: Polynomial, value: Double): Polynomial = ring { Polynomial(List(a.coefficients.size) { index -> a.coefficients[index] * value }) } + /** + * Evaluates the polynomial for the given value [arg]. + */ public operator fun Polynomial.invoke(arg: T): T = value(ring, arg) } diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/Integrator.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/Integrator.kt index ebc53ad2e..e421fb680 100644 --- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/Integrator.kt +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/Integrator.kt @@ -1,11 +1,11 @@ package space.kscience.kmath.integration /** - * A general interface for all integrators + * A general interface for all integrators. */ public interface Integrator { /** - * Run one integration pass and return a new [Integrand] with a new set of features + * Runs one integration pass and return a new [Integrand] with a new set of features. */ public fun integrate(integrand: I): I -} \ No newline at end of file +} diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/UnivariateIntegrand.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/UnivariateIntegrand.kt index 9389318e8..ca4bbf6b8 100644 --- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/UnivariateIntegrand.kt +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/UnivariateIntegrand.kt @@ -61,4 +61,4 @@ public fun UnivariateIntegrator.integrate( return integrate( UnivariateIntegrand(function, *features.toTypedArray()) ).value ?: error("Unexpected: no value after integration.") -} \ No newline at end of file +} diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/Interpolator.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/Interpolator.kt index fa978a9bc..9fad30abb 100644 --- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/Interpolator.kt +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/Interpolator.kt @@ -1,45 +1,50 @@ +@file:OptIn(UnstableKMathAPI::class) + package space.kscience.kmath.interpolation +import space.kscience.kmath.data.XYColumnarData import space.kscience.kmath.functions.PiecewisePolynomial import space.kscience.kmath.functions.value +import space.kscience.kmath.misc.UnstableKMathAPI import space.kscience.kmath.operations.Ring import space.kscience.kmath.structures.Buffer import space.kscience.kmath.structures.asBuffer -public fun interface Interpolator { - public fun interpolate(points: XYPointSet): (X) -> Y +public fun interface Interpolator { + public fun interpolate(points: XYColumnarData): (X) -> Y } -public interface PolynomialInterpolator> : Interpolator { +public interface PolynomialInterpolator> : Interpolator { public val algebra: Ring public fun getDefaultValue(): T = error("Out of bounds") - public fun interpolatePolynomials(points: XYPointSet): PiecewisePolynomial + public fun interpolatePolynomials(points: XYColumnarData): PiecewisePolynomial - override fun interpolate(points: XYPointSet): (T) -> T = { x -> + override fun interpolate(points: XYColumnarData): (T) -> T = { x -> interpolatePolynomials(points).value(algebra, x) ?: getDefaultValue() } } + public fun > PolynomialInterpolator.interpolatePolynomials( x: Buffer, y: Buffer, ): PiecewisePolynomial { - val pointSet = BufferXYPointSet(x, y) + val pointSet = XYColumnarData(x, y) return interpolatePolynomials(pointSet) } public fun > PolynomialInterpolator.interpolatePolynomials( data: Map, ): PiecewisePolynomial { - val pointSet = BufferXYPointSet(data.keys.toList().asBuffer(), data.values.toList().asBuffer()) + val pointSet = XYColumnarData(data.keys.toList().asBuffer(), data.values.toList().asBuffer()) return interpolatePolynomials(pointSet) } public fun > PolynomialInterpolator.interpolatePolynomials( data: List>, ): PiecewisePolynomial { - val pointSet = BufferXYPointSet(data.map { it.first }.asBuffer(), data.map { it.second }.asBuffer()) + val pointSet = XYColumnarData(data.map { it.first }.asBuffer(), data.map { it.second }.asBuffer()) return interpolatePolynomials(pointSet) } diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/LinearInterpolator.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/LinearInterpolator.kt index c939384e3..37d378ad0 100644 --- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/LinearInterpolator.kt +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/LinearInterpolator.kt @@ -1,16 +1,25 @@ package space.kscience.kmath.interpolation +import space.kscience.kmath.data.XYColumnarData import space.kscience.kmath.functions.OrderedPiecewisePolynomial import space.kscience.kmath.functions.PiecewisePolynomial import space.kscience.kmath.functions.Polynomial +import space.kscience.kmath.misc.UnstableKMathAPI import space.kscience.kmath.operations.Field import space.kscience.kmath.operations.invoke +@OptIn(UnstableKMathAPI::class) +internal fun > insureSorted(points: XYColumnarData<*, T, *>) { + for (i in 0 until points.size - 1) + require(points.x[i + 1] > points.x[i]) { "Input data is not sorted at index $i" } +} + /** * Reference JVM implementation: https://github.com/apache/commons-math/blob/master/src/main/java/org/apache/commons/math4/analysis/interpolation/LinearInterpolator.java */ public class LinearInterpolator>(public override val algebra: Field) : PolynomialInterpolator { - public override fun interpolatePolynomials(points: XYPointSet): PiecewisePolynomial = algebra { + @OptIn(UnstableKMathAPI::class) + public override fun interpolatePolynomials(points: XYColumnarData): PiecewisePolynomial = algebra { require(points.size > 0) { "Point array should not be empty" } insureSorted(points) diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/SplineInterpolator.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/SplineInterpolator.kt index ddbe743f0..3a3dfab59 100644 --- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/SplineInterpolator.kt +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/SplineInterpolator.kt @@ -1,15 +1,20 @@ package space.kscience.kmath.interpolation +import space.kscience.kmath.data.XYColumnarData import space.kscience.kmath.functions.OrderedPiecewisePolynomial import space.kscience.kmath.functions.PiecewisePolynomial import space.kscience.kmath.functions.Polynomial +import space.kscience.kmath.misc.UnstableKMathAPI import space.kscience.kmath.operations.Field import space.kscience.kmath.operations.invoke import space.kscience.kmath.structures.MutableBufferFactory /** - * Generic spline interpolator. Not recommended for performance critical places, use platform-specific and type specific ones. - * Based on https://github.com/apache/commons-math/blob/eb57d6d457002a0bb5336d789a3381a24599affe/src/main/java/org/apache/commons/math4/analysis/interpolation/SplineInterpolator.java + * Generic spline interpolator. Not recommended for performance critical places, use platform-specific and type + * specific ones. + * + * Based on + * https://github.com/apache/commons-math/blob/eb57d6d457002a0bb5336d789a3381a24599affe/src/main/java/org/apache/commons/math4/analysis/interpolation/SplineInterpolator.java */ public class SplineInterpolator>( public override val algebra: Field, @@ -17,7 +22,8 @@ public class SplineInterpolator>( ) : PolynomialInterpolator { //TODO possibly optimize zeroed buffers - public override fun interpolatePolynomials(points: XYPointSet): PiecewisePolynomial = algebra { + @OptIn(UnstableKMathAPI::class) + public override fun interpolatePolynomials(points: XYColumnarData): PiecewisePolynomial = algebra { require(points.size >= 3) { "Can't use spline interpolator with less than 3 points" } insureSorted(points) // Number of intervals. The number of data points is n + 1. diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/XYPointSet.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/XYPointSet.kt deleted file mode 100644 index 1ff7b6351..000000000 --- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/XYPointSet.kt +++ /dev/null @@ -1,53 +0,0 @@ -package space.kscience.kmath.interpolation - -import space.kscience.kmath.nd.Structure2D -import space.kscience.kmath.structures.Buffer - -public interface XYPointSet { - public val size: Int - public val x: Buffer - public val y: Buffer -} - -public interface XYZPointSet : XYPointSet { - public val z: Buffer -} - -internal fun > insureSorted(points: XYPointSet) { - for (i in 0 until points.size - 1) - require(points.x[i + 1] > points.x[i]) { "Input data is not sorted at index $i" } -} - -public class NDStructureColumn(public val structure: Structure2D, public val column: Int) : Buffer { - public override val size: Int - get() = structure.rowNum - - init { - require(column < structure.colNum) { "Column index is outside of structure column range" } - } - - public override operator fun get(index: Int): T = structure[index, column] - public override operator fun iterator(): Iterator = sequence { repeat(size) { yield(get(it)) } }.iterator() -} - -public class BufferXYPointSet( - public override val x: Buffer, - public override val y: Buffer, -) : XYPointSet { - public override val size: Int - get() = x.size - - init { - require(x.size == y.size) { "Sizes of x and y buffers should be the same" } - } -} - -public fun Structure2D.asXYPointSet(): XYPointSet { - require(shape[1] == 2) { "Structure second dimension should be of size 2" } - - return object : XYPointSet { - override val size: Int get() = this@asXYPointSet.shape[0] - override val x: Buffer get() = NDStructureColumn(this@asXYPointSet, 0) - override val y: Buffer get() = NDStructureColumn(this@asXYPointSet, 1) - } -} diff --git a/kmath-kotlingrad/src/main/kotlin/space/kscience/kmath/kotlingrad/DifferentiableMstExpression.kt b/kmath-kotlingrad/src/main/kotlin/space/kscience/kmath/kotlingrad/DifferentiableMstExpression.kt index 39a7248b4..1275b0c90 100644 --- a/kmath-kotlingrad/src/main/kotlin/space/kscience/kmath/kotlingrad/DifferentiableMstExpression.kt +++ b/kmath-kotlingrad/src/main/kotlin/space/kscience/kmath/kotlingrad/DifferentiableMstExpression.kt @@ -5,7 +5,7 @@ import space.kscience.kmath.ast.MST import space.kscience.kmath.ast.MstAlgebra import space.kscience.kmath.ast.MstExpression import space.kscience.kmath.expressions.DifferentiableExpression -import space.kscience.kmath.expressions.Symbol +import space.kscience.kmath.misc.Symbol import space.kscience.kmath.operations.NumericAlgebra /** @@ -18,8 +18,10 @@ import space.kscience.kmath.operations.NumericAlgebra * @param A the [NumericAlgebra] of [T]. * @property expr the underlying [MstExpression]. */ -public inline class DifferentiableMstExpression(public val expr: MstExpression) : - DifferentiableExpression> where A : NumericAlgebra, T : Number { +public inline class DifferentiableMstExpression( + public val expr: MstExpression, +) : DifferentiableExpression> where A : NumericAlgebra { + public constructor(algebra: A, mst: MST) : this(MstExpression(algebra, mst)) /** diff --git a/kmath-nd4j/README.md b/kmath-nd4j/README.md index 2771722eb..938d05c33 100644 --- a/kmath-nd4j/README.md +++ b/kmath-nd4j/README.md @@ -1,46 +1,40 @@ -# ND4J NDStructure implementation (`kmath-nd4j`) +# Module kmath-nd4j -This subproject implements the following features: +ND4J based implementations of KMath abstractions. - [nd4jarraystructure](#) : NDStructure wrapper for INDArray - [nd4jarrayrings](#) : Rings over Nd4jArrayStructure of Int and Long - [nd4jarrayfields](#) : Fields over Nd4jArrayStructure of Float and Double -> #### Artifact: -> -> This module artifact: `space.kscience:kmath-nd4j:0.3.0-dev-3`. -> -> Bintray release version: [ ![Download](https://api.bintray.com/packages/mipt-npm/kscience/kmath-nd4j/images/download.svg) ](https://bintray.com/mipt-npm/kscience/kmath-nd4j/_latestVersion) -> -> Bintray development version: [ ![Download](https://api.bintray.com/packages/mipt-npm/dev/kmath-nd4j/images/download.svg) ](https://bintray.com/mipt-npm/dev/kmath-nd4j/_latestVersion) -> -> **Gradle:** -> -> ```gradle -> repositories { -> maven { url 'https://repo.kotlin.link' } -> maven { url 'https://dl.bintray.com/hotkeytlt/maven' } -> maven { url "https://dl.bintray.com/kotlin/kotlin-eap" } // include for builds based on kotlin-eap -> } -> -> dependencies { -> implementation 'space.kscience:kmath-nd4j:0.3.0-dev-3' -> } -> ``` -> **Gradle Kotlin DSL:** -> -> ```kotlin -> repositories { -> maven("https://repo.kotlin.link") -> maven("https://dl.bintray.com/kotlin/kotlin-eap") // include for builds based on kotlin-eap -> maven("https://dl.bintray.com/hotkeytlt/maven") // required for a -> } -> -> dependencies { -> implementation("space.kscience:kmath-nd4j:0.3.0-dev-3") -> } -> ``` +## Artifact: + +The Maven coordinates of this project are `space.kscience:kmath-nd4j:0.3.0-dev-3`. + +**Gradle:** +```gradle +repositories { + maven { url 'https://repo.kotlin.link' } + maven { url 'https://dl.bintray.com/hotkeytlt/maven' } + maven { url "https://dl.bintray.com/kotlin/kotlin-eap" } // include for builds based on kotlin-eap +} + +dependencies { + implementation 'space.kscience:kmath-nd4j:0.3.0-dev-3' +} +``` +**Gradle Kotlin DSL:** +```kotlin +repositories { + maven("https://repo.kotlin.link") + maven("https://dl.bintray.com/kotlin/kotlin-eap") // include for builds based on kotlin-eap + maven("https://dl.bintray.com/hotkeytlt/maven") // required for a +} + +dependencies { + implementation("space.kscience:kmath-nd4j:0.3.0-dev-3") +} +``` ## Examples diff --git a/kmath-nd4j/docs/README-TEMPLATE.md b/kmath-nd4j/docs/README-TEMPLATE.md index 9783e74be..5f325cab5 100644 --- a/kmath-nd4j/docs/README-TEMPLATE.md +++ b/kmath-nd4j/docs/README-TEMPLATE.md @@ -1,6 +1,6 @@ -# ND4J NDStructure implementation (`kmath-nd4j`) +# Module kmath-nd4j -This subproject implements the following features: +ND4J based implementations of KMath abstractions. ${features} diff --git a/kmath-stat/build.gradle.kts b/kmath-stat/build.gradle.kts index 5b29a9e64..bc3890b1e 100644 --- a/kmath-stat/build.gradle.kts +++ b/kmath-stat/build.gradle.kts @@ -3,14 +3,6 @@ plugins { } kotlin.sourceSets { - all { - languageSettings.apply { - useExperimentalAnnotation("kotlinx.coroutines.FlowPreview") - useExperimentalAnnotation("kotlinx.coroutines.ExperimentalCoroutinesApi") - useExperimentalAnnotation("kotlinx.coroutines.ObsoleteCoroutinesApi") - } - } - commonMain { dependencies { api(project(":kmath-coroutines")) diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/optimization/FunctionOptimization.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/optimization/FunctionOptimization.kt new file mode 100644 index 000000000..528a5744e --- /dev/null +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/optimization/FunctionOptimization.kt @@ -0,0 +1,89 @@ +package space.kscience.kmath.optimization + +import space.kscience.kmath.expressions.AutoDiffProcessor +import space.kscience.kmath.expressions.DifferentiableExpression +import space.kscience.kmath.expressions.Expression +import space.kscience.kmath.expressions.ExpressionAlgebra +import space.kscience.kmath.misc.Symbol +import space.kscience.kmath.operations.ExtendedField +import space.kscience.kmath.structures.Buffer +import space.kscience.kmath.structures.indices + +/** + * A likelihood function optimization problem with provided derivatives + */ +public interface FunctionOptimization : Optimization { + /** + * The optimization direction. If true search for function maximum, if false, search for the minimum + */ + public var maximize: Boolean + + /** + * Define the initial guess for the optimization problem + */ + public fun initialGuess(map: Map) + + /** + * Set a differentiable expression as objective function as function and gradient provider + */ + public fun diffFunction(expression: DifferentiableExpression>) + + public companion object { + /** + * Generate a chi squared expression from given x-y-sigma data and inline model. Provides automatic differentiation + */ + public fun chiSquared( + autoDiff: AutoDiffProcessor>, + x: Buffer, + y: Buffer, + yErr: Buffer, + model: A.(I) -> I, + ): DifferentiableExpression> where A : ExtendedField, A : ExpressionAlgebra { + require(x.size == y.size) { "X and y buffers should be of the same size" } + require(y.size == yErr.size) { "Y and yErr buffer should of the same size" } + + return autoDiff.process { + var sum = zero + + x.indices.forEach { + val xValue = const(x[it]) + val yValue = const(y[it]) + val yErrValue = const(yErr[it]) + val modelValue = model(xValue) + sum += ((yValue - modelValue) / yErrValue).pow(2) + } + + sum + } + } + } +} + +/** + * Define a chi-squared-based objective function + */ +public fun FunctionOptimization.chiSquared( + autoDiff: AutoDiffProcessor>, + x: Buffer, + y: Buffer, + yErr: Buffer, + model: A.(I) -> I, +) where A : ExtendedField, A : ExpressionAlgebra { + val chiSquared = FunctionOptimization.chiSquared(autoDiff, x, y, yErr, model) + diffFunction(chiSquared) + maximize = false +} + +/** + * Optimize differentiable expression using specific [OptimizationProblemFactory] + */ +public fun > DifferentiableExpression>.optimizeWith( + factory: OptimizationProblemFactory, + vararg symbols: Symbol, + configuration: F.() -> Unit, +): OptimizationResult { + require(symbols.isNotEmpty()) { "Must provide a list of symbols for optimization" } + val problem = factory(symbols.toList(), configuration) + problem.diffFunction(this) + return problem.optimize() +} diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/optimization/NoDerivFunctionOptimization.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/optimization/NoDerivFunctionOptimization.kt new file mode 100644 index 000000000..b8785dd8c --- /dev/null +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/optimization/NoDerivFunctionOptimization.kt @@ -0,0 +1,69 @@ +package space.kscience.kmath.optimization + +import space.kscience.kmath.expressions.Expression +import space.kscience.kmath.misc.Symbol +import space.kscience.kmath.structures.Buffer +import space.kscience.kmath.structures.indices +import kotlin.math.pow + +/** + * A likelihood function optimization problem + */ +public interface NoDerivFunctionOptimization : Optimization { + /** + * The optimization direction. If true search for function maximum, if false, search for the minimum + */ + public var maximize: Boolean + + /** + * Define the initial guess for the optimization problem + */ + public fun initialGuess(map: Map) + + /** + * Set an objective function expression + */ + public fun function(expression: Expression) + + public companion object { + /** + * Generate a chi squared expression from given x-y-sigma model represented by an expression. Does not provide derivatives + */ + public fun chiSquared( + x: Buffer, + y: Buffer, + yErr: Buffer, + model: Expression, + xSymbol: Symbol = Symbol.x, + ): Expression { + require(x.size == y.size) { "X and y buffers should be of the same size" } + require(y.size == yErr.size) { "Y and yErr buffer should of the same size" } + + return Expression { arguments -> + x.indices.sumByDouble { + val xValue = x[it] + val yValue = y[it] + val yErrValue = yErr[it] + val modifiedArgs = arguments + (xSymbol to xValue) + val modelValue = model(modifiedArgs) + ((yValue - modelValue) / yErrValue).pow(2) + } + } + } + } +} + + +/** + * Optimize expression without derivatives using specific [OptimizationProblemFactory] + */ +public fun > Expression.noDerivOptimizeWith( + factory: OptimizationProblemFactory, + vararg symbols: Symbol, + configuration: F.() -> Unit, +): OptimizationResult { + require(symbols.isNotEmpty()) { "Must provide a list of symbols for optimization" } + val problem = factory(symbols.toList(), configuration) + problem.function(this) + return problem.optimize() +} diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/optimization/Optimization.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/optimization/Optimization.kt new file mode 100644 index 000000000..0b13e07ba --- /dev/null +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/optimization/Optimization.kt @@ -0,0 +1,44 @@ +package space.kscience.kmath.optimization + +import space.kscience.kmath.misc.Symbol + +public interface OptimizationFeature + +public class OptimizationResult( + public val point: Map, + public val value: T, + public val features: Set = emptySet(), +) { + override fun toString(): String { + return "OptimizationResult(point=$point, value=$value)" + } +} + +public operator fun OptimizationResult.plus( + feature: OptimizationFeature, +): OptimizationResult = OptimizationResult(point, value, features + feature) + +/** + * An optimization problem builder over [T] variables + */ +public interface Optimization { + + /** + * Update the problem from previous optimization run + */ + public fun update(result: OptimizationResult) + + /** + * Make an optimization run + */ + public fun optimize(): OptimizationResult +} + +public fun interface OptimizationProblemFactory> { + public fun build(symbols: List): P +} + +public operator fun > OptimizationProblemFactory.invoke( + symbols: List, + block: P.() -> Unit, +): P = build(symbols).apply(block) diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/optimization/XYFit.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/optimization/XYFit.kt new file mode 100644 index 000000000..c3106c819 --- /dev/null +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/optimization/XYFit.kt @@ -0,0 +1,40 @@ +package space.kscience.kmath.optimization + +import space.kscience.kmath.data.ColumnarData +import space.kscience.kmath.expressions.AutoDiffProcessor +import space.kscience.kmath.expressions.DifferentiableExpression +import space.kscience.kmath.expressions.Expression +import space.kscience.kmath.expressions.ExpressionAlgebra +import space.kscience.kmath.misc.Symbol +import space.kscience.kmath.misc.UnstableKMathAPI +import space.kscience.kmath.operations.ExtendedField +import space.kscience.kmath.operations.Field + +@UnstableKMathAPI +public interface XYFit : Optimization { + + public val algebra: Field + + /** + * Set X-Y data for this fit optionally including x and y errors + */ + public fun data( + dataSet: ColumnarData, + xSymbol: Symbol, + ySymbol: Symbol, + xErrSymbol: Symbol? = null, + yErrSymbol: Symbol? = null, + ) + + public fun model(model: (T) -> DifferentiableExpression) + + /** + * Set the differentiable model for this fit + */ + public fun model( + autoDiff: AutoDiffProcessor>, + modelFunction: A.(I) -> I, + ): Unit where A : ExtendedField, A : ExpressionAlgebra = model { arg -> + autoDiff.process { modelFunction(const(arg)) } + } +} \ No newline at end of file diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/Distribution.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/Distribution.kt index 25f6a87b5..095182160 100644 --- a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/Distribution.kt +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/Distribution.kt @@ -1,17 +1,29 @@ package space.kscience.kmath.stat +import kotlinx.coroutines.flow.first import space.kscience.kmath.chains.Chain import space.kscience.kmath.chains.collect import space.kscience.kmath.structures.Buffer import space.kscience.kmath.structures.BufferFactory -import space.kscience.kmath.structures.DoubleBuffer +import space.kscience.kmath.structures.IntBuffer +import space.kscience.kmath.structures.MutableBuffer +import kotlin.jvm.JvmName -public interface Sampler { +/** + * Sampler that generates chains of values of type [T]. + */ +public fun interface Sampler { + /** + * Generates a chain of samples. + * + * @param generator the randomness provider. + * @return the new chain. + */ public fun sample(generator: RandomGenerator): Chain } /** - * A distribution of typed objects + * A distribution of typed objects. */ public interface Distribution : Sampler { /** @@ -20,11 +32,7 @@ public interface Distribution : Sampler { */ public fun probability(arg: T): Double - /** - * Create a chain of samples from this distribution. - * The chain is not guaranteed to be stateless, but different sample chains should be independent. - */ - override fun sample(generator: RandomGenerator): Chain + public override fun sample(generator: RandomGenerator): Chain /** * An empty companion. Distribution factories should be written as its extensions @@ -63,16 +71,27 @@ public fun Sampler.sampleBuffer( //clear list from previous run tmp.clear() //Fill list - repeat(size) { - tmp.add(chain.next()) - } + repeat(size) { tmp += chain.next() } //return new buffer with elements from tmp bufferFactory(size) { tmp[it] } } } /** - * Generate a bunch of samples from real distributions + * Samples one value from this [Sampler]. */ +public suspend fun Sampler.next(generator: RandomGenerator): T = sample(generator).first() + +/** + * Generates [size] real samples and chunks them into some buffers. + */ +@JvmName("sampleRealBuffer") public fun Sampler.sampleBuffer(generator: RandomGenerator, size: Int): Chain> = - sampleBuffer(generator, size, ::DoubleBuffer) + sampleBuffer(generator, size, MutableBuffer.Companion::double) + +/** + * Generates [size] integer samples and chunks them into some buffers. + */ +@JvmName("sampleIntBuffer") +public fun Sampler.sampleBuffer(generator: RandomGenerator, size: Int): Chain> = + sampleBuffer(generator, size, ::IntBuffer) diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/FactorizedDistribution.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/FactorizedDistribution.kt index ff7a13652..3dd506b67 100644 --- a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/FactorizedDistribution.kt +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/FactorizedDistribution.kt @@ -14,7 +14,7 @@ public interface NamedDistribution : Distribution> public class FactorizedDistribution(public val distributions: Collection>) : NamedDistribution { override fun probability(arg: Map): Double = - distributions.fold(1.0) { acc, distr -> acc * distr.probability(arg) } + distributions.fold(1.0) { acc, dist -> acc * dist.probability(arg) } override fun sample(generator: RandomGenerator): Chain> { val chains = distributions.map { it.sample(generator) } @@ -38,6 +38,6 @@ public class DistributionBuilder { private val distributions = ArrayList>() public infix fun String.to(distribution: Distribution) { - distributions.add(NamedDistributionWrapper(this, distribution)) + distributions += NamedDistributionWrapper(this, distribution) } } diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/Fitting.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/Fitting.kt deleted file mode 100644 index b006c8ba2..000000000 --- a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/Fitting.kt +++ /dev/null @@ -1,63 +0,0 @@ -package space.kscience.kmath.stat - -import space.kscience.kmath.expressions.* -import space.kscience.kmath.operations.ExtendedField -import space.kscience.kmath.structures.Buffer -import space.kscience.kmath.structures.indices -import kotlin.math.pow - -public object Fitting { - - /** - * Generate a chi squared expression from given x-y-sigma data and inline model. Provides automatic differentiation - */ - public fun chiSquared( - autoDiff: AutoDiffProcessor>, - x: Buffer, - y: Buffer, - yErr: Buffer, - model: A.(I) -> I, - ): DifferentiableExpression> where A : ExtendedField, A : ExpressionAlgebra { - require(x.size == y.size) { "X and y buffers should be of the same size" } - require(y.size == yErr.size) { "Y and yErr buffer should of the same size" } - - return autoDiff.process { - var sum = zero - - x.indices.forEach { - val xValue = const(x[it]) - val yValue = const(y[it]) - val yErrValue = const(yErr[it]) - val modelValue = model(xValue) - sum += ((yValue - modelValue) / yErrValue).pow(2) - } - - sum - } - } - - /** - * Generate a chi squared expression from given x-y-sigma model represented by an expression. Does not provide derivatives - */ - public fun chiSquared( - x: Buffer, - y: Buffer, - yErr: Buffer, - model: Expression, - xSymbol: Symbol = StringSymbol("x"), - ): Expression { - require(x.size == y.size) { "X and y buffers should be of the same size" } - require(y.size == yErr.size) { "Y and yErr buffer should of the same size" } - - return Expression { arguments -> - x.indices.sumByDouble { - val xValue = x[it] - val yValue = y[it] - val yErrValue = yErr[it] - val modifiedArgs = arguments + (xSymbol to xValue) - val modelValue = model(modifiedArgs) - ((yValue - modelValue) / yErrValue).pow(2) - } - } - } -} diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/OptimizationProblem.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/OptimizationProblem.kt deleted file mode 100644 index ffa24fa98..000000000 --- a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/OptimizationProblem.kt +++ /dev/null @@ -1,88 +0,0 @@ -package space.kscience.kmath.stat - -import space.kscience.kmath.expressions.DifferentiableExpression -import space.kscience.kmath.expressions.Expression -import space.kscience.kmath.expressions.Symbol - -public interface OptimizationFeature - -public class OptimizationResult( - public val point: Map, - public val value: T, - public val features: Set = emptySet(), -) { - override fun toString(): String { - return "OptimizationResult(point=$point, value=$value)" - } -} - -public operator fun OptimizationResult.plus( - feature: OptimizationFeature, -): OptimizationResult = OptimizationResult(point, value, features + feature) - -/** - * A configuration builder for optimization problem - */ -public interface OptimizationProblem { - /** - * Define the initial guess for the optimization problem - */ - public fun initialGuess(map: Map) - - /** - * Set an objective function expression - */ - public fun expression(expression: Expression) - - /** - * Set a differentiable expression as objective function as function and gradient provider - */ - public fun diffExpression(expression: DifferentiableExpression>) - - /** - * Update the problem from previous optimization run - */ - public fun update(result: OptimizationResult) - - /** - * Make an optimization run - */ - public fun optimize(): OptimizationResult -} - -public fun interface OptimizationProblemFactory> { - public fun build(symbols: List): P -} - -public operator fun > OptimizationProblemFactory.invoke( - symbols: List, - block: P.() -> Unit, -): P = build(symbols).apply(block) - -/** - * Optimize expression without derivatives using specific [OptimizationProblemFactory] - */ -public fun > Expression.optimizeWith( - factory: OptimizationProblemFactory, - vararg symbols: Symbol, - configuration: F.() -> Unit, -): OptimizationResult { - require(symbols.isNotEmpty()) { "Must provide a list of symbols for optimization" } - val problem = factory(symbols.toList(), configuration) - problem.expression(this) - return problem.optimize() -} - -/** - * Optimize differentiable expression using specific [OptimizationProblemFactory] - */ -public fun > DifferentiableExpression>.optimizeWith( - factory: OptimizationProblemFactory, - vararg symbols: Symbol, - configuration: F.() -> Unit, -): OptimizationResult { - require(symbols.isNotEmpty()) { "Must provide a list of symbols for optimization" } - val problem = factory(symbols.toList(), configuration) - problem.diffExpression(this) - return problem.optimize() -} diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/RandomChain.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/RandomChain.kt index 6e1f36c8a..2f117a035 100644 --- a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/RandomChain.kt +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/RandomChain.kt @@ -1,17 +1,22 @@ package space.kscience.kmath.stat +import space.kscience.kmath.chains.BlockingDoubleChain +import space.kscience.kmath.chains.BlockingIntChain import space.kscience.kmath.chains.Chain /** * A possibly stateful chain producing random values. + * + * @property generator the underlying [RandomGenerator] instance. */ public class RandomChain( public val generator: RandomGenerator, - private val gen: suspend RandomGenerator.() -> R, + private val gen: suspend RandomGenerator.() -> R ) : Chain { override suspend fun next(): R = generator.gen() - override fun fork(): Chain = RandomChain(generator.fork(), gen) } public fun RandomGenerator.chain(gen: suspend RandomGenerator.() -> R): RandomChain = RandomChain(this, gen) +public fun Chain.blocking(): BlockingDoubleChain = object : Chain by this, BlockingDoubleChain {} +public fun Chain.blocking(): BlockingIntChain = object : Chain by this, BlockingIntChain {} diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/RandomGenerator.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/RandomGenerator.kt index 1a4f3b75d..bad2334e9 100644 --- a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/RandomGenerator.kt +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/RandomGenerator.kt @@ -82,6 +82,8 @@ public interface RandomGenerator { /** * Implements [RandomGenerator] by delegating all operations to [Random]. + * + * @property random the underlying [Random] object. */ public class DefaultGenerator(public val random: Random = Random) : RandomGenerator { public override fun nextBoolean(): Boolean = random.nextBoolean() diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/SamplerAlgebra.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/SamplerAlgebra.kt index c5ec99dae..25ec7eca6 100644 --- a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/SamplerAlgebra.kt +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/SamplerAlgebra.kt @@ -8,16 +8,28 @@ import space.kscience.kmath.operations.Group import space.kscience.kmath.operations.ScaleOperations import space.kscience.kmath.operations.invoke -public class BasicSampler(public val chainBuilder: (RandomGenerator) -> Chain) : Sampler { - public override fun sample(generator: RandomGenerator): Chain = chainBuilder(generator) -} - +/** + * Implements [Sampler] by sampling only certain [value]. + * + * @property value the value to sample. + */ public class ConstantSampler(public val value: T) : Sampler { public override fun sample(generator: RandomGenerator): Chain = ConstantChain(value) } /** - * A space for samplers. Allows to perform simple operations on distributions + * Implements [Sampler] by delegating sampling to value of [chainBuilder]. + * + * @property chainBuilder the provider of [Chain]. + */ +public class BasicSampler(public val chainBuilder: (RandomGenerator) -> Chain) : Sampler { + public override fun sample(generator: RandomGenerator): Chain = chainBuilder(generator) +} + +/** + * A space of samplers. Allows to perform simple operations on distributions. + * + * @property algebra the space to provide addition and scalar multiplication for [T]. */ public class SamplerSpace(public val algebra: S) : Group>, ScaleOperations> where S : Group, S : ScaleOperations { @@ -29,8 +41,10 @@ public class SamplerSpace(public val algebra: S) : Group> } public override fun scale(a: Sampler, value: Double): Sampler = BasicSampler { generator -> - a.sample(generator).map { algebra { it * value } } + a.sample(generator).map { a -> + algebra { a * value } + } } - override fun Sampler.unaryMinus(): Sampler = scale(this, -1.0) + public override fun Sampler.unaryMinus(): Sampler = scale(this, -1.0) } diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/distributions/NormalDistribution.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/distributions/NormalDistribution.kt new file mode 100644 index 000000000..6515cbaa7 --- /dev/null +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/distributions/NormalDistribution.kt @@ -0,0 +1,41 @@ +package space.kscience.kmath.stat.distributions + +import space.kscience.kmath.chains.Chain +import space.kscience.kmath.stat.RandomGenerator +import space.kscience.kmath.stat.UnivariateDistribution +import space.kscience.kmath.stat.internal.InternalErf +import space.kscience.kmath.stat.samplers.GaussianSampler +import space.kscience.kmath.stat.samplers.NormalizedGaussianSampler +import space.kscience.kmath.stat.samplers.ZigguratNormalizedGaussianSampler +import kotlin.math.* + +/** + * Implements [UnivariateDistribution] for the normal (gaussian) distribution. + */ +public inline class NormalDistribution(public val sampler: GaussianSampler) : UnivariateDistribution { + public constructor( + mean: Double, + standardDeviation: Double, + normalized: NormalizedGaussianSampler = ZigguratNormalizedGaussianSampler.of(), + ) : this(GaussianSampler.of(mean, standardDeviation, normalized)) + + public override fun probability(arg: Double): Double { + val x1 = (arg - sampler.mean) / sampler.standardDeviation + return exp(-0.5 * x1 * x1 - (ln(sampler.standardDeviation) + 0.5 * ln(2 * PI))) + } + + public override fun sample(generator: RandomGenerator): Chain = sampler.sample(generator) + + public override fun cumulative(arg: Double): Double { + val dev = arg - sampler.mean + + return when { + abs(dev) > 40 * sampler.standardDeviation -> if (dev < 0) 0.0 else 1.0 + else -> 0.5 * InternalErf.erfc(-dev / (sampler.standardDeviation * SQRT2)) + } + } + + private companion object { + private val SQRT2 = sqrt(2.0) + } +} diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/internal/InternalErf.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/internal/InternalErf.kt new file mode 100644 index 000000000..4e1623867 --- /dev/null +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/internal/InternalErf.kt @@ -0,0 +1,15 @@ +package space.kscience.kmath.stat.internal + +import kotlin.math.abs + +/** + * Based on Commons Math implementation. + * See [https://commons.apache.org/proper/commons-math/javadocs/api-3.3/org/apache/commons/math3/special/Erf.html]. + */ +internal object InternalErf { + fun erfc(x: Double): Double { + if (abs(x) > 40) return if (x > 0) 0.0 else 2.0 + val ret = InternalGamma.regularizedGammaQ(0.5, x * x, 10000) + return if (x < 0) 2 - ret else ret + } +} \ No newline at end of file diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/internal/InternalGamma.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/internal/InternalGamma.kt new file mode 100644 index 000000000..4f5adbe97 --- /dev/null +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/internal/InternalGamma.kt @@ -0,0 +1,238 @@ +package space.kscience.kmath.stat.internal + +import kotlin.math.* + +private abstract class ContinuedFraction protected constructor() { + protected abstract fun getA(n: Int, x: Double): Double + protected abstract fun getB(n: Int, x: Double): Double + + fun evaluate(x: Double, maxIterations: Int): Double { + val small = 1e-50 + var hPrev = getA(0, x) + if (hPrev == 0.0 || abs(0.0 - hPrev) <= small) hPrev = small + var n = 1 + var dPrev = 0.0 + var cPrev = hPrev + var hN = hPrev + + while (n < maxIterations) { + val a = getA(n, x) + val b = getB(n, x) + var dN = a + b * dPrev + if (dN == 0.0 || abs(0.0 - dN) <= small) dN = small + var cN = a + b / cPrev + if (cN == 0.0 || abs(0.0 - cN) <= small) cN = small + dN = 1 / dN + val deltaN = cN * dN + hN = hPrev * deltaN + check(!hN.isInfinite()) { "hN is infinite" } + check(!hN.isNaN()) { "hN is NaN" } + if (abs(deltaN - 1.0) < 10e-9) break + dPrev = dN + cPrev = cN + hPrev = hN + n++ + } + + check(n < maxIterations) { "n is more than maxIterations" } + return hN + } +} + +internal object InternalGamma { + const val LANCZOS_G = 607.0 / 128.0 + + private val LANCZOS = doubleArrayOf( + 0.99999999999999709182, + 57.156235665862923517, + -59.597960355475491248, + 14.136097974741747174, + -0.49191381609762019978, + .33994649984811888699e-4, + .46523628927048575665e-4, + -.98374475304879564677e-4, + .15808870322491248884e-3, + -.21026444172410488319e-3, + .21743961811521264320e-3, + -.16431810653676389022e-3, + .84418223983852743293e-4, + -.26190838401581408670e-4, + .36899182659531622704e-5 + ) + + private val HALF_LOG_2_PI = 0.5 * ln(2.0 * PI) + private const val INV_GAMMA1P_M1_A0 = .611609510448141581788E-08 + private const val INV_GAMMA1P_M1_A1 = .624730830116465516210E-08 + private const val INV_GAMMA1P_M1_B1 = .203610414066806987300E+00 + private const val INV_GAMMA1P_M1_B2 = .266205348428949217746E-01 + private const val INV_GAMMA1P_M1_B3 = .493944979382446875238E-03 + private const val INV_GAMMA1P_M1_B4 = -.851419432440314906588E-05 + private const val INV_GAMMA1P_M1_B5 = -.643045481779353022248E-05 + private const val INV_GAMMA1P_M1_B6 = .992641840672773722196E-06 + private const val INV_GAMMA1P_M1_B7 = -.607761895722825260739E-07 + private const val INV_GAMMA1P_M1_B8 = .195755836614639731882E-09 + private const val INV_GAMMA1P_M1_P0 = .6116095104481415817861E-08 + private const val INV_GAMMA1P_M1_P1 = .6871674113067198736152E-08 + private const val INV_GAMMA1P_M1_P2 = .6820161668496170657918E-09 + private const val INV_GAMMA1P_M1_P3 = .4686843322948848031080E-10 + private const val INV_GAMMA1P_M1_P4 = .1572833027710446286995E-11 + private const val INV_GAMMA1P_M1_P5 = -.1249441572276366213222E-12 + private const val INV_GAMMA1P_M1_P6 = .4343529937408594255178E-14 + private const val INV_GAMMA1P_M1_Q1 = .3056961078365221025009E+00 + private const val INV_GAMMA1P_M1_Q2 = .5464213086042296536016E-01 + private const val INV_GAMMA1P_M1_Q3 = .4956830093825887312020E-02 + private const val INV_GAMMA1P_M1_Q4 = .2692369466186361192876E-03 + private const val INV_GAMMA1P_M1_C = -.422784335098467139393487909917598E+00 + private const val INV_GAMMA1P_M1_C0 = .577215664901532860606512090082402E+00 + private const val INV_GAMMA1P_M1_C1 = -.655878071520253881077019515145390E+00 + private const val INV_GAMMA1P_M1_C2 = -.420026350340952355290039348754298E-01 + private const val INV_GAMMA1P_M1_C3 = .166538611382291489501700795102105E+00 + private const val INV_GAMMA1P_M1_C4 = -.421977345555443367482083012891874E-01 + private const val INV_GAMMA1P_M1_C5 = -.962197152787697356211492167234820E-02 + private const val INV_GAMMA1P_M1_C6 = .721894324666309954239501034044657E-02 + private const val INV_GAMMA1P_M1_C7 = -.116516759185906511211397108401839E-02 + private const val INV_GAMMA1P_M1_C8 = -.215241674114950972815729963053648E-03 + private const val INV_GAMMA1P_M1_C9 = .128050282388116186153198626328164E-03 + private const val INV_GAMMA1P_M1_C10 = -.201348547807882386556893914210218E-04 + private const val INV_GAMMA1P_M1_C11 = -.125049348214267065734535947383309E-05 + private const val INV_GAMMA1P_M1_C12 = .113302723198169588237412962033074E-05 + private const val INV_GAMMA1P_M1_C13 = -.205633841697760710345015413002057E-06 + + fun logGamma(x: Double): Double = when { + x.isNaN() || x <= 0.0 -> Double.NaN + x < 0.5 -> logGamma1p(x) - ln(x) + x <= 2.5 -> logGamma1p(x - 0.5 - 0.5) + + x <= 8.0 -> { + val n = floor(x - 1.5).toInt() + val prod = (1..n).fold(1.0, { prod, i -> prod * (x - i) }) + logGamma1p(x - (n + 1)) + ln(prod) + } + + else -> { + val tmp = x + LANCZOS_G + .5 + (x + .5) * ln(tmp) - tmp + HALF_LOG_2_PI + ln(lanczos(x) / x) + } + } + + private fun regularizedGammaP( + a: Double, + x: Double, + maxIterations: Int = Int.MAX_VALUE + ): Double = when { + a.isNaN() || x.isNaN() || a <= 0.0 || x < 0.0 -> Double.NaN + x == 0.0 -> 0.0 + x >= a + 1 -> 1.0 - regularizedGammaQ(a, x, maxIterations) + + else -> { + // calculate series + var n = 0.0 // current element index + var an = 1.0 / a // n-th element in the series + var sum = an // partial sum + + while (abs(an / sum) > 10e-15 && n < maxIterations && sum < Double.POSITIVE_INFINITY) { + // compute next element in the series + n += 1.0 + an *= x / (a + n) + + // update partial sum + sum += an + } + + when { + n >= maxIterations -> throw error("Maximal iterations is exceeded $maxIterations") + sum.isInfinite() -> 1.0 + else -> exp(-x + a * ln(x) - logGamma(a)) * sum + } + } + } + + fun regularizedGammaQ( + a: Double, + x: Double, + maxIterations: Int = Int.MAX_VALUE + ): Double = when { + a.isNaN() || x.isNaN() || a <= 0.0 || x < 0.0 -> Double.NaN + x == 0.0 -> 1.0 + x < a + 1.0 -> 1.0 - regularizedGammaP(a, x, maxIterations) + + else -> 1.0 / object : ContinuedFraction() { + override fun getA(n: Int, x: Double): Double = 2.0 * n + 1.0 - a + x + override fun getB(n: Int, x: Double): Double = n * (a - n) + }.evaluate(x, maxIterations) * exp(-x + a * ln(x) - logGamma(a)) + } + + private fun lanczos(x: Double): Double = + (LANCZOS.size - 1 downTo 1).sumByDouble { LANCZOS[it] / (x + it) } + LANCZOS[0] + + private fun invGamma1pm1(x: Double): Double { + require(x >= -0.5) + require(x <= 1.5) + val ret: Double + val t = if (x <= 0.5) x else x - 0.5 - 0.5 + + if (t < 0.0) { + val a = INV_GAMMA1P_M1_A0 + t * INV_GAMMA1P_M1_A1 + var b = INV_GAMMA1P_M1_B8 + b = INV_GAMMA1P_M1_B7 + t * b + b = INV_GAMMA1P_M1_B6 + t * b + b = INV_GAMMA1P_M1_B5 + t * b + b = INV_GAMMA1P_M1_B4 + t * b + b = INV_GAMMA1P_M1_B3 + t * b + b = INV_GAMMA1P_M1_B2 + t * b + b = INV_GAMMA1P_M1_B1 + t * b + b = 1.0 + t * b + var c = INV_GAMMA1P_M1_C13 + t * (a / b) + c = INV_GAMMA1P_M1_C12 + t * c + c = INV_GAMMA1P_M1_C11 + t * c + c = INV_GAMMA1P_M1_C10 + t * c + c = INV_GAMMA1P_M1_C9 + t * c + c = INV_GAMMA1P_M1_C8 + t * c + c = INV_GAMMA1P_M1_C7 + t * c + c = INV_GAMMA1P_M1_C6 + t * c + c = INV_GAMMA1P_M1_C5 + t * c + c = INV_GAMMA1P_M1_C4 + t * c + c = INV_GAMMA1P_M1_C3 + t * c + c = INV_GAMMA1P_M1_C2 + t * c + c = INV_GAMMA1P_M1_C1 + t * c + c = INV_GAMMA1P_M1_C + t * c + ret = (if (x > 0.5) t * c / x else x * (c + 0.5 + 0.5)) + } else { + var p = INV_GAMMA1P_M1_P6 + p = INV_GAMMA1P_M1_P5 + t * p + p = INV_GAMMA1P_M1_P4 + t * p + p = INV_GAMMA1P_M1_P3 + t * p + p = INV_GAMMA1P_M1_P2 + t * p + p = INV_GAMMA1P_M1_P1 + t * p + p = INV_GAMMA1P_M1_P0 + t * p + var q = INV_GAMMA1P_M1_Q4 + q = INV_GAMMA1P_M1_Q3 + t * q + q = INV_GAMMA1P_M1_Q2 + t * q + q = INV_GAMMA1P_M1_Q1 + t * q + q = 1.0 + t * q + var c = INV_GAMMA1P_M1_C13 + p / q * t + c = INV_GAMMA1P_M1_C12 + t * c + c = INV_GAMMA1P_M1_C11 + t * c + c = INV_GAMMA1P_M1_C10 + t * c + c = INV_GAMMA1P_M1_C9 + t * c + c = INV_GAMMA1P_M1_C8 + t * c + c = INV_GAMMA1P_M1_C7 + t * c + c = INV_GAMMA1P_M1_C6 + t * c + c = INV_GAMMA1P_M1_C5 + t * c + c = INV_GAMMA1P_M1_C4 + t * c + c = INV_GAMMA1P_M1_C3 + t * c + c = INV_GAMMA1P_M1_C2 + t * c + c = INV_GAMMA1P_M1_C1 + t * c + c = INV_GAMMA1P_M1_C0 + t * c + ret = (if (x > 0.5) t / x * (c - 0.5 - 0.5) else x * c) + } + + return ret + } + + private fun logGamma1p(x: Double): Double { + require(x >= -0.5) + require(x <= 1.5) + return -ln1p(invGamma1pm1(x)) + } +} diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/internal/InternalUtils.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/internal/InternalUtils.kt new file mode 100644 index 000000000..722eee946 --- /dev/null +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/internal/InternalUtils.kt @@ -0,0 +1,70 @@ +package space.kscience.kmath.stat.internal + +import kotlin.math.ln +import kotlin.math.min + +internal object InternalUtils { + private val FACTORIALS = longArrayOf( + 1L, 1L, 2L, + 6L, 24L, 120L, + 720L, 5040L, 40320L, + 362880L, 3628800L, 39916800L, + 479001600L, 6227020800L, 87178291200L, + 1307674368000L, 20922789888000L, 355687428096000L, + 6402373705728000L, 121645100408832000L, 2432902008176640000L + ) + + private const val BEGIN_LOG_FACTORIALS = 2 + + fun factorial(n: Int): Long = FACTORIALS[n] + + fun validateProbabilities(probabilities: DoubleArray?): Double { + require(!(probabilities == null || probabilities.isEmpty())) { "Probabilities must not be empty." } + + val sumProb = probabilities.sumByDouble { prob -> + require(!(prob < 0 || prob.isInfinite() || prob.isNaN())) { "Invalid probability: $prob" } + prob + } + + require(!(sumProb.isInfinite() || sumProb <= 0)) { "Invalid sum of probabilities: $sumProb" } + return sumProb + } + + class FactorialLog private constructor(numValues: Int, cache: DoubleArray?) { + private val logFactorials: DoubleArray = DoubleArray(numValues) + + init { + val endCopy: Int + + if (cache != null && cache.size > BEGIN_LOG_FACTORIALS) { + // Copy available values. + endCopy = min(cache.size, numValues) + + cache.copyInto( + logFactorials, + BEGIN_LOG_FACTORIALS, + BEGIN_LOG_FACTORIALS, endCopy + ) + } else + // All values to be computed + endCopy = BEGIN_LOG_FACTORIALS + + // Compute remaining values. + (endCopy until numValues).forEach { i -> + if (i < FACTORIALS.size) + logFactorials[i] = ln(FACTORIALS[i].toDouble()) + else + logFactorials[i] = logFactorials[i - 1] + ln(i.toDouble()) + } + } + + fun value(n: Int): Double { + if (n < logFactorials.size) return logFactorials[n] + return if (n < FACTORIALS.size) ln(FACTORIALS[n].toDouble()) else InternalGamma.logGamma(n + 1.0) + } + + companion object { + fun create(): FactorialLog = FactorialLog(0, null) + } + } +} \ No newline at end of file diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/samplers/AhrensDieterExponentialSampler.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/samplers/AhrensDieterExponentialSampler.kt new file mode 100644 index 000000000..504c6b881 --- /dev/null +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/samplers/AhrensDieterExponentialSampler.kt @@ -0,0 +1,73 @@ +package space.kscience.kmath.stat.samplers + +import space.kscience.kmath.chains.Chain +import space.kscience.kmath.stat.RandomGenerator +import space.kscience.kmath.stat.Sampler +import space.kscience.kmath.stat.chain +import space.kscience.kmath.stat.internal.InternalUtils +import kotlin.math.ln +import kotlin.math.pow + +/** + * Sampling from an [exponential distribution](http://mathworld.wolfram.com/ExponentialDistribution.html). + * + * Based on Commons RNG implementation. + * See [https://commons.apache.org/proper/commons-rng/commons-rng-sampling/apidocs/org/apache/commons/rng/sampling/distribution/AhrensDieterExponentialSampler.html]. + */ +public class AhrensDieterExponentialSampler private constructor(public val mean: Double) : Sampler { + public override fun sample(generator: RandomGenerator): Chain = generator.chain { + // Step 1: + var a = 0.0 + var u = nextDouble() + + // Step 2 and 3: + while (u < 0.5) { + a += EXPONENTIAL_SA_QI[0] + u *= 2.0 + } + + // Step 4 (now u >= 0.5): + u += u - 1 + // Step 5: + if (u <= EXPONENTIAL_SA_QI[0]) return@chain mean * (a + u) + // Step 6: + var i = 0 // Should be 1, be we iterate before it in while using 0. + var u2 = nextDouble() + var umin = u2 + + // Step 7 and 8: + do { + ++i + u2 = nextDouble() + if (u2 < umin) umin = u2 + // Step 8: + } while (u > EXPONENTIAL_SA_QI[i]) // Ensured to exit since EXPONENTIAL_SA_QI[MAX] = 1. + + mean * (a + umin * EXPONENTIAL_SA_QI[0]) + } + + override fun toString(): String = "Ahrens-Dieter Exponential deviate" + + public companion object { + private val EXPONENTIAL_SA_QI by lazy { DoubleArray(16) } + + init { + /** + * Filling EXPONENTIAL_SA_QI table. + * Note that we don't want qi = 0 in the table. + */ + val ln2 = ln(2.0) + var qi = 0.0 + + EXPONENTIAL_SA_QI.indices.forEach { i -> + qi += ln2.pow(i + 1.0) / InternalUtils.factorial(i + 1) + EXPONENTIAL_SA_QI[i] = qi + } + } + + public fun of(mean: Double): AhrensDieterExponentialSampler { + require(mean > 0) { "mean is not strictly positive: $mean" } + return AhrensDieterExponentialSampler(mean) + } + } +} diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/samplers/AhrensDieterMarsagliaTsangGammaSampler.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/samplers/AhrensDieterMarsagliaTsangGammaSampler.kt new file mode 100644 index 000000000..81182f6cd --- /dev/null +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/samplers/AhrensDieterMarsagliaTsangGammaSampler.kt @@ -0,0 +1,120 @@ +package space.kscience.kmath.stat.samplers + +import space.kscience.kmath.chains.Chain +import space.kscience.kmath.stat.RandomGenerator +import space.kscience.kmath.stat.Sampler +import space.kscience.kmath.stat.chain +import space.kscience.kmath.stat.next +import kotlin.math.* + +/** + * Sampling from the [gamma distribution](http://mathworld.wolfram.com/GammaDistribution.html). + * - For 0 < alpha < 1: + * Ahrens, J. H. and Dieter, U., Computer methods for sampling from gamma, beta, Poisson and binomial distributions, Computing, 12, 223-246, 1974. + * - For alpha >= 1: + * Marsaglia and Tsang, A Simple Method for Generating Gamma Variables. ACM Transactions on Mathematical Software, Volume 26 Issue 3, September, 2000. + * + * Based on Commons RNG implementation. + * + * See [https://commons.apache.org/proper/commons-rng/commons-rng-sampling/apidocs/org/apache/commons/rng/sampling/distribution/AhrensDieterMarsagliaTsangGammaSampler.html]. + */ +public class AhrensDieterMarsagliaTsangGammaSampler private constructor( + alpha: Double, + theta: Double +) : Sampler { + private val delegate: BaseGammaSampler = + if (alpha < 1) AhrensDieterGammaSampler(alpha, theta) else MarsagliaTsangGammaSampler(alpha, theta) + + private abstract class BaseGammaSampler internal constructor( + protected val alpha: Double, + protected val theta: Double + ) : Sampler { + init { + require(alpha > 0) { "alpha is not strictly positive: $alpha" } + require(theta > 0) { "theta is not strictly positive: $theta" } + } + + override fun toString(): String = "Ahrens-Dieter-Marsaglia-Tsang Gamma deviate" + } + + private class AhrensDieterGammaSampler(alpha: Double, theta: Double) : + BaseGammaSampler(alpha, theta) { + private val oneOverAlpha: Double = 1.0 / alpha + private val bGSOptim: Double = 1.0 + alpha / E + + override fun sample(generator: RandomGenerator): Chain = generator.chain { + var x: Double + + // [1]: p. 228, Algorithm GS. + while (true) { + // Step 1: + val u = generator.nextDouble() + val p = bGSOptim * u + + if (p <= 1) { + // Step 2: + x = p.pow(oneOverAlpha) + val u2 = generator.nextDouble() + + if (u2 > exp(-x)) // Reject. + continue + + break + } + + // Step 3: + x = -ln((bGSOptim - p) * oneOverAlpha) + val u2: Double = generator.nextDouble() + if (u2 <= x.pow(alpha - 1.0)) break + // Reject and continue. + } + + x * theta + } + } + + private class MarsagliaTsangGammaSampler(alpha: Double, theta: Double) : + BaseGammaSampler(alpha, theta) { + private val dOptim: Double + private val cOptim: Double + private val gaussian: NormalizedGaussianSampler + + init { + gaussian = ZigguratNormalizedGaussianSampler.of() + dOptim = alpha - ONE_THIRD + cOptim = ONE_THIRD / sqrt(dOptim) + } + + override fun sample(generator: RandomGenerator): Chain = generator.chain { + var v: Double + + while (true) { + val x = gaussian.next(generator) + val oPcTx = 1 + cOptim * x + v = oPcTx * oPcTx * oPcTx + if (v <= 0) continue + val x2 = x * x + val u = generator.nextDouble() + // Squeeze. + if (u < 1 - 0.0331 * x2 * x2) break + if (ln(u) < 0.5 * x2 + dOptim * (1 - v + ln(v))) break + } + + theta * dOptim * v + } + + companion object { + private const val ONE_THIRD = 1.0 / 3.0 + } + } + + public override fun sample(generator: RandomGenerator): Chain = delegate.sample(generator) + public override fun toString(): String = delegate.toString() + + public companion object { + public fun of( + alpha: Double, + theta: Double + ): Sampler = AhrensDieterMarsagliaTsangGammaSampler(alpha, theta) + } +} \ No newline at end of file diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/samplers/AliasMethodDiscreteSampler.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/samplers/AliasMethodDiscreteSampler.kt new file mode 100644 index 000000000..cae97db65 --- /dev/null +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/samplers/AliasMethodDiscreteSampler.kt @@ -0,0 +1,286 @@ +package space.kscience.kmath.stat.samplers + +import space.kscience.kmath.chains.Chain +import space.kscience.kmath.stat.RandomGenerator +import space.kscience.kmath.stat.Sampler +import space.kscience.kmath.stat.chain +import space.kscience.kmath.stat.internal.InternalUtils +import kotlin.math.ceil +import kotlin.math.max +import kotlin.math.min + +/** + * Distribution sampler that uses the Alias method. It can be used to sample from n values each with an associated + * probability. This implementation is based on the detailed explanation of the alias method by Keith Schartz and + * implements Vose's algorithm. + * + * Vose, M.D., A linear algorithm for generating random numbers with a given distribution, IEEE Transactions on + * Software Engineering, 17, 972-975, 1991. he algorithm will sample values in O(1) time after a pre-processing step + * of O(n) time. + * + * The alias tables are constructed using fraction probabilities with an assumed denominator of 253. In the generic + * case sampling uses UniformRandomProvider.nextInt(int) and the upper 53-bits from UniformRandomProvider.nextLong(). + * + * Zero padding the input probabilities can be used to make more sampling more efficient. Any zero entry will always be + * aliased removing the requirement to compute a long. Increased sampling speed comes at the cost of increased storage + * space. The algorithm requires approximately 12 bytes of storage per input probability, that is n * 12 for size n. + * Zero-padding only requires 4 bytes of storage per padded value as the probability is known to be zero. + * + * An optimisation is performed for small table sizes that are a power of 2. In this case the sampling uses 1 or 2 + * calls from UniformRandomProvider.nextInt() to generate up to 64-bits for creation of an 11-bit index and 53-bits + * for the long. This optimisation requires a generator with a high cycle length for the lower order bits. + * + * Larger table sizes that are a power of 2 will benefit from fast algorithms for UniformRandomProvider.nextInt(int) + * that exploit the power of 2. + * + * Based on Commons RNG implementation. + * See [https://commons.apache.org/proper/commons-rng/commons-rng-sampling/apidocs/org/apache/commons/rng/sampling/distribution/AliasMethodDiscreteSampler.html]. + */ +public open class AliasMethodDiscreteSampler private constructor( + // Deliberate direct storage of input arrays + protected val probability: LongArray, + protected val alias: IntArray +) : Sampler { + + private class SmallTableAliasMethodDiscreteSampler( + probability: LongArray, + alias: IntArray + ) : AliasMethodDiscreteSampler(probability, alias) { + // Assume the table size is a power of 2 and create the mask + private val mask: Int = alias.size - 1 + + override fun sample(generator: RandomGenerator): Chain = generator.chain { + val bits = generator.nextInt() + // Isolate lower bits + val j = bits and mask + + // Optimisation for zero-padded input tables + if (j >= probability.size) + // No probability must use the alias + return@chain alias[j] + + // Create a uniform random deviate as a long. + // This replicates functionality from the o.a.c.rng.core.utils.NumberFactory.makeLong + val longBits = generator.nextInt().toLong() shl 32 or (bits.toLong() and hex_ffffffff) + // Choose between the two. Use a 53-bit long for the probability. + if (longBits ushr 11 < probability[j]) j else alias[j] + } + + private companion object { + private const val hex_ffffffff = 4294967295L + } + } + + public override fun sample(generator: RandomGenerator): Chain = generator.chain { + // This implements the algorithm as per Vose (1991): + // v = uniform() in [0, 1) + // j = uniform(n) in [0, n) + // if v < prob[j] then + // return j + // else + // return alias[j] + val j = generator.nextInt(alias.size) + + // Optimisation for zero-padded input tables + // No probability must use the alias + if (j >= probability.size) return@chain alias[j] + + // Note: We could check the probability before computing a deviate. + // p(j) == 0 => alias[j] + // p(j) == 1 => j + // However it is assumed these edge cases are rare: + // + // The probability table will be 1 for approximately 1/n samples, i.e. only the + // last unpaired probability. This is only worth checking for when the table size (n) + // is small. But in that case the user should zero-pad the table for performance. + // + // The probability table will be 0 when an input probability was zero. We + // will assume this is also rare if modelling a discrete distribution where + // all samples are possible. The edge case for zero-padded tables is handled above. + + // Choose between the two. Use a 53-bit long for the probability. + if (generator.nextLong() ushr 11 < probability[j]) j else alias[j] + } + + public override fun toString(): String = "Alias method" + + public companion object { + private const val DEFAULT_ALPHA = 0 + private const val ZERO = 0.0 + private const val ONE_AS_NUMERATOR = 1L shl 53 + private const val CONVERT_TO_NUMERATOR: Double = ONE_AS_NUMERATOR.toDouble() + private const val MAX_SMALL_POWER_2_SIZE = 1 shl 11 + + public fun of( + probabilities: DoubleArray, + alpha: Int = DEFAULT_ALPHA + ): Sampler { + // The Alias method balances N categories with counts around the mean into N sections, + // each allocated 'mean' observations. + // + // Consider 4 categories with counts 6,3,2,1. The histogram can be balanced into a + // 2D array as 4 sections with a height of the mean: + // + // 6 + // 6 + // 6 + // 63 => 6366 -- + // 632 6326 |-- mean + // 6321 6321 -- + // + // section abcd + // + // Each section is divided as: + // a: 6=1/1 + // b: 3=1/1 + // c: 2=2/3; 6=1/3 (6 is the alias) + // d: 1=1/3; 6=2/3 (6 is the alias) + // + // The sample is obtained by randomly selecting a section, then choosing which category + // from the pair based on a uniform random deviate. + val sumProb = InternalUtils.validateProbabilities(probabilities) + // Allow zero-padding + val n = computeSize(probabilities.size, alpha) + // Partition into small and large by splitting on the average. + val mean = sumProb / n + // The cardinality of smallSize + largeSize = n. + // So fill the same array from either end. + val indices = IntArray(n) + var large = n + var small = 0 + + probabilities.indices.forEach { i -> + if (probabilities[i] >= mean) indices[--large] = i else indices[small++] = i + } + + small = fillRemainingIndices(probabilities.size, indices, small) + // This may be smaller than the input length if the probabilities were already padded. + val nonZeroIndex = findLastNonZeroIndex(probabilities) + // The probabilities are modified so use a copy. + // Note: probabilities are required only up to last nonZeroIndex + val remainingProbabilities = probabilities.copyOf(nonZeroIndex + 1) + // Allocate the final tables. + // Probability table may be truncated (when zero padded). + // The alias table is full length. + val probability = LongArray(remainingProbabilities.size) + val alias = IntArray(n) + + // This loop uses each large in turn to fill the alias table for small probabilities that + // do not reach the requirement to fill an entire section alone (i.e. p < mean). + // Since the sum of the small should be less than the sum of the large it should use up + // all the small first. However floating point round-off can result in + // misclassification of items as small or large. The Vose algorithm handles this using + // a while loop conditioned on the size of both sets and a subsequent loop to use + // unpaired items. + while (large != n && small != 0) { + // Index of the small and the large probabilities. + val j = indices[--small] + val k = indices[large++] + + // Optimisation for zero-padded input: + // p(j) = 0 above the last nonZeroIndex + if (j > nonZeroIndex) + // The entire amount for the section is taken from the alias. + remainingProbabilities[k] -= mean + else { + val pj = remainingProbabilities[j] + // Item j is a small probability that is below the mean. + // Compute the weight of the section for item j: pj / mean. + // This is scaled by 2^53 and the ceiling function used to round-up + // the probability to a numerator of a fraction in the range [1,2^53]. + // Ceiling ensures non-zero values. + probability[j] = ceil(CONVERT_TO_NUMERATOR * (pj / mean)).toLong() + // The remaining amount for the section is taken from the alias. + // Effectively: probabilities[k] -= (mean - pj) + remainingProbabilities[k] += pj - mean + } + + // If not j then the alias is k + alias[j] = k + + // Add the remaining probability from large to the appropriate list. + if (remainingProbabilities[k] >= mean) indices[--large] = k else indices[small++] = k + } + + // Final loop conditions to consume unpaired items. + // Note: The large set should never be non-empty but this can occur due to round-off + // error so consume from both. + fillTable(probability, alias, indices, 0, small) + fillTable(probability, alias, indices, large, n) + + // Change the algorithm for small power of 2 sized tables + return if (isSmallPowerOf2(n)) + SmallTableAliasMethodDiscreteSampler(probability, alias) + else + AliasMethodDiscreteSampler(probability, alias) + } + + private fun fillRemainingIndices(length: Int, indices: IntArray, small: Int): Int { + var updatedSmall = small + (length until indices.size).forEach { i -> indices[updatedSmall++] = i } + return updatedSmall + } + + private fun findLastNonZeroIndex(probabilities: DoubleArray): Int { + // No bounds check is performed when decrementing as the array contains at least one + // value above zero. + var nonZeroIndex = probabilities.size - 1 + while (probabilities[nonZeroIndex] == ZERO) nonZeroIndex-- + return nonZeroIndex + } + + private fun computeSize(length: Int, alpha: Int): Int { + // If No padding + if (alpha < 0) return length + // Use the number of leading zeros function to find the next power of 2, + // i.e. ceil(log2(x)) + var pow2 = 32 - numberOfLeadingZeros(length - 1) + // Increase by the alpha. Clip this to limit to a positive integer (2^30) + pow2 = min(30, pow2 + alpha) + // Use max to handle a length above the highest possible power of 2 + return max(length, 1 shl pow2) + } + + private fun fillTable( + probability: LongArray, + alias: IntArray, + indices: IntArray, + start: Int, + end: Int + ) = (start until end).forEach { i -> + val index = indices[i] + probability[index] = ONE_AS_NUMERATOR + alias[index] = index + } + + private fun isSmallPowerOf2(n: Int): Boolean = n <= MAX_SMALL_POWER_2_SIZE && n and n - 1 == 0 + + private fun numberOfLeadingZeros(i: Int): Int { + var mutI = i + if (mutI <= 0) return if (mutI == 0) 32 else 0 + var n = 31 + + if (mutI >= 1 shl 16) { + n -= 16 + mutI = mutI ushr 16 + } + + if (mutI >= 1 shl 8) { + n -= 8 + mutI = mutI ushr 8 + } + + if (mutI >= 1 shl 4) { + n -= 4 + mutI = mutI ushr 4 + } + + if (mutI >= 1 shl 2) { + n -= 2 + mutI = mutI ushr 2 + } + + return n - (mutI ushr 1) + } + } +} diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/samplers/BoxMullerNormalizedGaussianSampler.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/samplers/BoxMullerNormalizedGaussianSampler.kt new file mode 100644 index 000000000..04beb448d --- /dev/null +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/samplers/BoxMullerNormalizedGaussianSampler.kt @@ -0,0 +1,48 @@ +package space.kscience.kmath.stat.samplers + +import space.kscience.kmath.chains.Chain +import space.kscience.kmath.stat.RandomGenerator +import space.kscience.kmath.stat.Sampler +import space.kscience.kmath.stat.chain +import kotlin.math.* + +/** + * [Box-Muller algorithm](https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform) for sampling from a Gaussian + * distribution. + * + * Based on Commons RNG implementation. + * See [https://commons.apache.org/proper/commons-rng/commons-rng-sampling/apidocs/org/apache/commons/rng/sampling/distribution/BoxMullerNormalizedGaussianSampler.html]. + */ +public class BoxMullerNormalizedGaussianSampler private constructor() : NormalizedGaussianSampler, Sampler { + private var nextGaussian: Double = Double.NaN + + public override fun sample(generator: RandomGenerator): Chain = generator.chain { + val random: Double + + if (nextGaussian.isNaN()) { + // Generate a pair of Gaussian numbers. + val x = nextDouble() + val y = nextDouble() + val alpha = 2 * PI * x + val r = sqrt(-2 * ln(y)) + // Return the first element of the generated pair. + random = r * cos(alpha) + // Keep second element of the pair for next invocation. + nextGaussian = r * sin(alpha) + } else { + // Use the second element of the pair (generated at the + // previous invocation). + random = nextGaussian + // Both elements of the pair have been used. + nextGaussian = Double.NaN + } + + random + } + + public override fun toString(): String = "Box-Muller normalized Gaussian deviate" + + public companion object { + public fun of(): BoxMullerNormalizedGaussianSampler = BoxMullerNormalizedGaussianSampler() + } +} diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/samplers/GaussianSampler.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/samplers/GaussianSampler.kt new file mode 100644 index 000000000..eba26cfb5 --- /dev/null +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/samplers/GaussianSampler.kt @@ -0,0 +1,43 @@ +package space.kscience.kmath.stat.samplers + +import space.kscience.kmath.chains.Chain +import space.kscience.kmath.chains.map +import space.kscience.kmath.stat.RandomGenerator +import space.kscience.kmath.stat.Sampler + +/** + * Sampling from a Gaussian distribution with given mean and standard deviation. + * + * Based on Commons RNG implementation. + * See [https://commons.apache.org/proper/commons-rng/commons-rng-sampling/apidocs/org/apache/commons/rng/sampling/distribution/GaussianSampler.html]. + * + * @property mean the mean of the distribution. + * @property standardDeviation the variance of the distribution. + */ +public class GaussianSampler private constructor( + public val mean: Double, + public val standardDeviation: Double, + private val normalized: NormalizedGaussianSampler +) : Sampler { + public override fun sample(generator: RandomGenerator): Chain = normalized + .sample(generator) + .map { standardDeviation * it + mean } + + override fun toString(): String = "Gaussian deviate [$normalized]" + + public companion object { + public fun of( + mean: Double, + standardDeviation: Double, + normalized: NormalizedGaussianSampler = ZigguratNormalizedGaussianSampler.of() + ): GaussianSampler { + require(standardDeviation > 0.0) { "standard deviation is not strictly positive: $standardDeviation" } + + return GaussianSampler( + mean, + standardDeviation, + normalized + ) + } + } +} diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/samplers/KempSmallMeanPoissonSampler.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/samplers/KempSmallMeanPoissonSampler.kt new file mode 100644 index 000000000..1d7f90023 --- /dev/null +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/samplers/KempSmallMeanPoissonSampler.kt @@ -0,0 +1,63 @@ +package space.kscience.kmath.stat.samplers + +import space.kscience.kmath.chains.Chain +import space.kscience.kmath.stat.RandomGenerator +import space.kscience.kmath.stat.Sampler +import space.kscience.kmath.stat.chain +import kotlin.math.exp + +/** + * Sampler for the Poisson distribution. + * - Kemp, A, W, (1981) Efficient Generation of Logarithmically Distributed Pseudo-Random Variables. Journal of the Royal Statistical Society. Vol. 30, No. 3, pp. 249-253. + * This sampler is suitable for mean < 40. For large means, LargeMeanPoissonSampler should be used instead. + * + * Note: The algorithm uses a recurrence relation to compute the Poisson probability and a rolling summation for the cumulative probability. When the mean is large the initial probability (Math.exp(-mean)) is zero and an exception is raised by the constructor. + * + * Sampling uses 1 call to UniformRandomProvider.nextDouble(). This method provides an alternative to the SmallMeanPoissonSampler for slow generators of double. + * + * Based on Commons RNG implementation. + * See [https://commons.apache.org/proper/commons-rng/commons-rng-sampling/apidocs/org/apache/commons/rng/sampling/distribution/KempSmallMeanPoissonSampler.html]. + */ +public class KempSmallMeanPoissonSampler private constructor( + private val p0: Double, + private val mean: Double +) : Sampler { + public override fun sample(generator: RandomGenerator): Chain = generator.chain { + // Note on the algorithm: + // - X is the unknown sample deviate (the output of the algorithm) + // - x is the current value from the distribution + // - p is the probability of the current value x, p(X=x) + // - u is effectively the cumulative probability that the sample X + // is equal or above the current value x, p(X>=x) + // So if p(X>=x) > p(X=x) the sample must be above x, otherwise it is x + var u = nextDouble() + var x = 0 + var p = p0 + + while (u > p) { + u -= p + // Compute the next probability using a recurrence relation. + // p(x+1) = p(x) * mean / (x+1) + p *= mean / ++x + // The algorithm listed in Kemp (1981) does not check that the rolling probability + // is positive. This check is added to ensure no errors when the limit of the summation + // 1 - sum(p(x)) is above 0 due to cumulative error in floating point arithmetic. + if (p == 0.0) return@chain x + } + + x + } + + public override fun toString(): String = "Kemp Small Mean Poisson deviate" + + public companion object { + public fun of(mean: Double): KempSmallMeanPoissonSampler { + require(mean > 0) { "Mean is not strictly positive: $mean" } + val p0 = exp(-mean) + // Probability must be positive. As mean increases then p(0) decreases. + require(p0 > 0) { "No probability for mean: $mean" } + return KempSmallMeanPoissonSampler(p0, mean) + } + } +} + diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/samplers/LargeMeanPoissonSampler.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/samplers/LargeMeanPoissonSampler.kt new file mode 100644 index 000000000..de1e7cc89 --- /dev/null +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/samplers/LargeMeanPoissonSampler.kt @@ -0,0 +1,130 @@ +package space.kscience.kmath.stat.samplers + +import space.kscience.kmath.chains.Chain +import space.kscience.kmath.chains.ConstantChain +import space.kscience.kmath.stat.RandomGenerator +import space.kscience.kmath.stat.Sampler +import space.kscience.kmath.stat.chain +import space.kscience.kmath.stat.internal.InternalUtils +import space.kscience.kmath.stat.next +import kotlin.math.* + +/** + * Sampler for the Poisson distribution. + * - For large means, we use the rejection algorithm described in + * Devroye, Luc. (1981).The Computer Generation of Poisson Random Variables + * Computing vol. 26 pp. 197-207. + * + * This sampler is suitable for mean >= 40. + * + * Based on Commons RNG implementation. + * See [https://commons.apache.org/proper/commons-rng/commons-rng-sampling/apidocs/org/apache/commons/rng/sampling/distribution/LargeMeanPoissonSampler.html]. + */ +public class LargeMeanPoissonSampler private constructor(public val mean: Double) : Sampler { + private val exponential: Sampler = AhrensDieterExponentialSampler.of(1.0) + private val gaussian: Sampler = ZigguratNormalizedGaussianSampler.of() + private val factorialLog: InternalUtils.FactorialLog = NO_CACHE_FACTORIAL_LOG + private val lambda: Double = floor(mean) + private val logLambda: Double = ln(lambda) + private val logLambdaFactorial: Double = getFactorialLog(lambda.toInt()) + private val delta: Double = sqrt(lambda * ln(32 * lambda / PI + 1)) + private val halfDelta: Double = delta / 2 + private val twolpd: Double = 2 * lambda + delta + private val c1: Double = 1 / (8 * lambda) + private val a1: Double = sqrt(PI * twolpd) * exp(c1) + private val a2: Double = twolpd / delta * exp(-delta * (1 + delta) / twolpd) + private val aSum: Double = a1 + a2 + 1 + private val p1: Double = a1 / aSum + private val p2: Double = a2 / aSum + + private val smallMeanPoissonSampler: Sampler = if (mean - lambda < Double.MIN_VALUE) + NO_SMALL_MEAN_POISSON_SAMPLER + else // Not used. + KempSmallMeanPoissonSampler.of(mean - lambda) + + public override fun sample(generator: RandomGenerator): Chain = generator.chain { + // This will never be null. It may be a no-op delegate that returns zero. + val y2 = smallMeanPoissonSampler.next(generator) + var x: Double + var y: Double + var v: Double + var a: Int + var t: Double + var qr: Double + var qa: Double + + while (true) { + // Step 1: + val u = generator.nextDouble() + + if (u <= p1) { + // Step 2: + val n = gaussian.next(generator) + x = n * sqrt(lambda + halfDelta) - 0.5 + if (x > delta || x < -lambda) continue + y = if (x < 0) floor(x) else ceil(x) + val e = exponential.next(generator) + v = -e - 0.5 * n * n + c1 + } else { + // Step 3: + if (u > p1 + p2) { + y = lambda + break + } + + x = delta + twolpd / delta * exponential.next(generator) + y = ceil(x) + v = -exponential.next(generator) - delta * (x + 1) / twolpd + } + + // The Squeeze Principle + // Step 4.1: + a = if (x < 0) 1 else 0 + t = y * (y + 1) / (2 * lambda) + + // Step 4.2 + if (v < -t && a == 0) { + y += lambda + break + } + + // Step 4.3: + qr = t * ((2 * y + 1) / (6 * lambda) - 1) + qa = qr - t * t / (3 * (lambda + a * (y + 1))) + + // Step 4.4: + if (v < qa) { + y += lambda + break + } + + // Step 4.5: + if (v > qr) continue + + // Step 4.6: + if (v < y * logLambda - getFactorialLog((y + lambda).toInt()) + logLambdaFactorial) { + y += lambda + break + } + } + + min(y2 + y.toLong(), Int.MAX_VALUE.toLong()).toInt() + } + + private fun getFactorialLog(n: Int): Double = factorialLog.value(n) + public override fun toString(): String = "Large Mean Poisson deviate" + + public companion object { + private const val MAX_MEAN: Double = 0.5 * Int.MAX_VALUE + private val NO_CACHE_FACTORIAL_LOG: InternalUtils.FactorialLog = InternalUtils.FactorialLog.create() + + private val NO_SMALL_MEAN_POISSON_SAMPLER: Sampler = Sampler { ConstantChain(0) } + + public fun of(mean: Double): LargeMeanPoissonSampler { + require(mean >= 1) { "mean is not >= 1: $mean" } + // The algorithm is not valid if Math.floor(mean) is not an integer. + require(mean <= MAX_MEAN) { "mean $mean > $MAX_MEAN" } + return LargeMeanPoissonSampler(mean) + } + } +} diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/samplers/MarsagliaNormalizedGaussianSampler.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/samplers/MarsagliaNormalizedGaussianSampler.kt new file mode 100644 index 000000000..8a659642f --- /dev/null +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/samplers/MarsagliaNormalizedGaussianSampler.kt @@ -0,0 +1,61 @@ +package space.kscience.kmath.stat.samplers + +import space.kscience.kmath.chains.Chain +import space.kscience.kmath.stat.RandomGenerator +import space.kscience.kmath.stat.Sampler +import space.kscience.kmath.stat.chain +import kotlin.math.ln +import kotlin.math.sqrt + +/** + * [Marsaglia polar method](https://en.wikipedia.org/wiki/Marsaglia_polar_method) for sampling from a Gaussian + * distribution with mean 0 and standard deviation 1. This is a variation of the algorithm implemented in + * [BoxMullerNormalizedGaussianSampler]. + * + * Based on Commons RNG implementation. + * See [https://commons.apache.org/proper/commons-rng/commons-rng-sampling/apidocs/org/apache/commons/rng/sampling/distribution/MarsagliaNormalizedGaussianSampler.html] + */ +public class MarsagliaNormalizedGaussianSampler private constructor() : NormalizedGaussianSampler, Sampler { + private var nextGaussian = Double.NaN + + public override fun sample(generator: RandomGenerator): Chain = generator.chain { + if (nextGaussian.isNaN()) { + val alpha: Double + var x: Double + + // Rejection scheme for selecting a pair that lies within the unit circle. + while (true) { + // Generate a pair of numbers within [-1 , 1). + x = 2.0 * generator.nextDouble() - 1.0 + val y = 2.0 * generator.nextDouble() - 1.0 + val r2 = x * x + y * y + + if (r2 < 1 && r2 > 0) { + // Pair (x, y) is within unit circle. + alpha = sqrt(-2 * ln(r2) / r2) + // Keep second element of the pair for next invocation. + nextGaussian = alpha * y + // Return the first element of the generated pair. + break + } + // Pair is not within the unit circle: Generate another one. + } + + // Return the first element of the generated pair. + alpha * x + } else { + // Use the second element of the pair (generated at the + // previous invocation). + val r = nextGaussian + // Both elements of the pair have been used. + nextGaussian = Double.NaN + r + } + } + + public override fun toString(): String = "Box-Muller (with rejection) normalized Gaussian deviate" + + public companion object { + public fun of(): MarsagliaNormalizedGaussianSampler = MarsagliaNormalizedGaussianSampler() + } +} diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/samplers/NormalizedGaussianSampler.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/samplers/NormalizedGaussianSampler.kt new file mode 100644 index 000000000..4eb3d60e0 --- /dev/null +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/samplers/NormalizedGaussianSampler.kt @@ -0,0 +1,9 @@ +package space.kscience.kmath.stat.samplers + +import space.kscience.kmath.stat.Sampler + +/** + * Marker interface for a sampler that generates values from an N(0,1) + * [Gaussian distribution](https://en.wikipedia.org/wiki/Normal_distribution). + */ +public interface NormalizedGaussianSampler : Sampler diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/samplers/PoissonSampler.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/samplers/PoissonSampler.kt new file mode 100644 index 000000000..0c0234892 --- /dev/null +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/samplers/PoissonSampler.kt @@ -0,0 +1,30 @@ +package space.kscience.kmath.stat.samplers + +import space.kscience.kmath.chains.Chain +import space.kscience.kmath.stat.RandomGenerator +import space.kscience.kmath.stat.Sampler + +/** + * Sampler for the Poisson distribution. + * - For small means, a Poisson process is simulated using uniform deviates, as described in + * Knuth (1969). Seminumerical Algorithms. The Art of Computer Programming, Volume 2. Chapter 3.4.1.F.3 + * Important integer-valued distributions: The Poisson distribution. Addison Wesley. + * The Poisson process (and hence, the returned value) is bounded by 1000 * mean. + * - For large means, we use the rejection algorithm described in + * Devroye, Luc. (1981). The Computer Generation of Poisson Random Variables Computing vol. 26 pp. 197-207. + * + * Based on Commons RNG implementation. + * See [https://commons.apache.org/proper/commons-rng/commons-rng-sampling/apidocs/org/apache/commons/rng/sampling/distribution/PoissonSampler.html]. + */ +public class PoissonSampler private constructor(mean: Double) : Sampler { + private val poissonSamplerDelegate: Sampler = of(mean) + public override fun sample(generator: RandomGenerator): Chain = poissonSamplerDelegate.sample(generator) + public override fun toString(): String = poissonSamplerDelegate.toString() + + public companion object { + private const val PIVOT = 40.0 + + public fun of(mean: Double): Sampler =// Each sampler should check the input arguments. + if (mean < PIVOT) SmallMeanPoissonSampler.of(mean) else LargeMeanPoissonSampler.of(mean) + } +} diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/samplers/SmallMeanPoissonSampler.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/samplers/SmallMeanPoissonSampler.kt new file mode 100644 index 000000000..0fe7ff161 --- /dev/null +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/samplers/SmallMeanPoissonSampler.kt @@ -0,0 +1,50 @@ +package space.kscience.kmath.stat.samplers + +import space.kscience.kmath.chains.Chain +import space.kscience.kmath.stat.RandomGenerator +import space.kscience.kmath.stat.Sampler +import space.kscience.kmath.stat.chain +import kotlin.math.ceil +import kotlin.math.exp + +/** + * Sampler for the Poisson distribution. + * - For small means, a Poisson process is simulated using uniform deviates, as described in + * Knuth (1969). Seminumerical Algorithms. The Art of Computer Programming, Volume 2. Chapter 3.4.1.F.3 Important + * integer-valued distributions: The Poisson distribution. Addison Wesley. + * - The Poisson process (and hence, the returned value) is bounded by 1000 * mean. + * This sampler is suitable for mean < 40. For large means, [LargeMeanPoissonSampler] should be used instead. + * + * Based on Commons RNG implementation. + * + * See [https://commons.apache.org/proper/commons-rng/commons-rng-sampling/apidocs/org/apache/commons/rng/sampling/distribution/SmallMeanPoissonSampler.html]. + */ +public class SmallMeanPoissonSampler private constructor(mean: Double) : Sampler { + private val p0: Double = exp(-mean) + + private val limit: Int = (if (p0 > 0) + ceil(1000 * mean) + else + throw IllegalArgumentException("No p(x=0) probability for mean: $mean")).toInt() + + public override fun sample(generator: RandomGenerator): Chain = generator.chain { + var n = 0 + var r = 1.0 + + while (n < limit) { + r *= nextDouble() + if (r >= p0) n++ else break + } + + n + } + + public override fun toString(): String = "Small Mean Poisson deviate" + + public companion object { + public fun of(mean: Double): SmallMeanPoissonSampler { + require(mean > 0) { "mean is not strictly positive: $mean" } + return SmallMeanPoissonSampler(mean) + } + } +} diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/samplers/ZigguratNormalizedGaussianSampler.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/samplers/ZigguratNormalizedGaussianSampler.kt new file mode 100644 index 000000000..90815209f --- /dev/null +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/samplers/ZigguratNormalizedGaussianSampler.kt @@ -0,0 +1,88 @@ +package space.kscience.kmath.stat.samplers + +import space.kscience.kmath.chains.Chain +import space.kscience.kmath.stat.RandomGenerator +import space.kscience.kmath.stat.Sampler +import space.kscience.kmath.stat.chain +import kotlin.math.* + +/** + * [Marsaglia and Tsang "Ziggurat"](https://en.wikipedia.org/wiki/Ziggurat_algorithm) method for sampling from a + * Gaussian distribution with mean 0 and standard deviation 1. The algorithm is explained in this paper and this + * implementation has been adapted from the C code provided therein. + * + * Based on Commons RNG implementation. + * See [https://commons.apache.org/proper/commons-rng/commons-rng-sampling/apidocs/org/apache/commons/rng/sampling/distribution/ZigguratNormalizedGaussianSampler.html]. + */ +public class ZigguratNormalizedGaussianSampler private constructor() : + NormalizedGaussianSampler, Sampler { + + private fun sampleOne(generator: RandomGenerator): Double { + val j = generator.nextLong() + val i = (j and LAST.toLong()).toInt() + return if (abs(j) < K[i]) j * W[i] else fix(generator, j, i) + } + + public override fun sample(generator: RandomGenerator): Chain = generator.chain { sampleOne(this) } + public override fun toString(): String = "Ziggurat normalized Gaussian deviate" + + private fun fix(generator: RandomGenerator, hz: Long, iz: Int): Double { + var x = hz * W[iz] + + return when { + iz == 0 -> { + var y: Double + + do { + y = -ln(generator.nextDouble()) + x = -ln(generator.nextDouble()) * ONE_OVER_R + } while (y + y < x * x) + + val out = R + x + if (hz > 0) out else -out + } + + F[iz] + generator.nextDouble() * (F[iz - 1] - F[iz]) < gauss(x) -> x + else -> sampleOne(generator) + } + } + + public companion object { + private const val R: Double = 3.442619855899 + private const val ONE_OVER_R: Double = 1 / R + private const val V: Double = 9.91256303526217e-3 + private val MAX: Double = 2.0.pow(63.0) + private val ONE_OVER_MAX: Double = 1.0 / MAX + private const val LEN: Int = 128 + private const val LAST: Int = LEN - 1 + private val K: LongArray = LongArray(LEN) + private val W: DoubleArray = DoubleArray(LEN) + private val F: DoubleArray = DoubleArray(LEN) + + init { + // Filling the tables. + var d = R + var t = d + var fd = gauss(d) + val q = V / fd + K[0] = (d / q * MAX).toLong() + K[1] = 0 + W[0] = q * ONE_OVER_MAX + W[LAST] = d * ONE_OVER_MAX + F[0] = 1.0 + F[LAST] = fd + + (LAST - 1 downTo 1).forEach { i -> + d = sqrt(-2 * ln(V / d + fd)) + fd = gauss(d) + K[i + 1] = (d / t * MAX).toLong() + t = d + F[i] = fd + W[i] = d * ONE_OVER_MAX + } + } + + public fun of(): ZigguratNormalizedGaussianSampler = ZigguratNormalizedGaussianSampler() + private fun gauss(x: Double): Double = exp(-0.5 * x * x) + } +} diff --git a/kmath-stat/src/jvmMain/kotlin/space/kscience/kmath/stat/RandomSourceGenerator.kt b/kmath-stat/src/jvmMain/kotlin/space/kscience/kmath/stat/RandomSourceGenerator.kt index 9e752d571..f6a5d6605 100644 --- a/kmath-stat/src/jvmMain/kotlin/space/kscience/kmath/stat/RandomSourceGenerator.kt +++ b/kmath-stat/src/jvmMain/kotlin/space/kscience/kmath/stat/RandomSourceGenerator.kt @@ -3,10 +3,14 @@ package space.kscience.kmath.stat import org.apache.commons.rng.UniformRandomProvider import org.apache.commons.rng.simple.RandomSource -public class RandomSourceGenerator(public val source: RandomSource, seed: Long?) : RandomGenerator { - internal val random: UniformRandomProvider = seed?.let { - RandomSource.create(source, seed) - } ?: RandomSource.create(source) +/** + * Implements [RandomGenerator] by delegating all operations to [RandomSource]. + * + * @property source the underlying [RandomSource] object. + */ +public class RandomSourceGenerator internal constructor(public val source: RandomSource, seed: Long?) : RandomGenerator { + internal val random: UniformRandomProvider = seed?.let { RandomSource.create(source, seed) } + ?: RandomSource.create(source) public override fun nextBoolean(): Boolean = random.nextBoolean() public override fun nextDouble(): Double = random.nextDouble() @@ -23,22 +27,84 @@ public class RandomSourceGenerator(public val source: RandomSource, seed: Long?) public override fun fork(): RandomGenerator = RandomSourceGenerator(source, nextLong()) } +/** + * Implements [UniformRandomProvider] by delegating all operations to [RandomGenerator]. + * + * @property generator the underlying [RandomGenerator] object. + */ public inline class RandomGeneratorProvider(public val generator: RandomGenerator) : UniformRandomProvider { + /** + * Generates a [Boolean] value. + * + * @return the next random value. + */ public override fun nextBoolean(): Boolean = generator.nextBoolean() + + /** + * Generates a [Float] value between 0 and 1. + * + * @return the next random value between 0 and 1. + */ public override fun nextFloat(): Float = generator.nextDouble().toFloat() - public override fun nextBytes(bytes: ByteArray) { - generator.fillBytes(bytes) - } + /** + * Generates [Byte] values and places them into a user-supplied array. + * + * The number of random bytes produced is equal to the length of the the byte array. + * + * @param bytes byte array in which to put the random bytes. + */ + public override fun nextBytes(bytes: ByteArray): Unit = generator.fillBytes(bytes) + /** + * Generates [Byte] values and places them into a user-supplied array. + * + * The array is filled with bytes extracted from random integers. This implies that the number of random bytes + * generated may be larger than the length of the byte array. + * + * @param bytes the array in which to put the generated bytes. + * @param start the index at which to start inserting the generated bytes. + * @param len the number of bytes to insert. + */ public override fun nextBytes(bytes: ByteArray, start: Int, len: Int) { generator.fillBytes(bytes, start, start + len) } + /** + * Generates an [Int] value. + * + * @return the next random value. + */ public override fun nextInt(): Int = generator.nextInt() + + /** + * Generates an [Int] value between 0 (inclusive) and the specified value (exclusive). + * + * @param n the bound on the random number to be returned. Must be positive. + * @return a random integer between 0 (inclusive) and [n] (exclusive). + */ public override fun nextInt(n: Int): Int = generator.nextInt(n) + + /** + * Generates a [Double] value between 0 and 1. + * + * @return the next random value between 0 and 1. + */ public override fun nextDouble(): Double = generator.nextDouble() + + /** + * Generates a [Long] value. + * + * @return the next random value. + */ public override fun nextLong(): Long = generator.nextLong() + + /** + * Generates a [Long] value between 0 (inclusive) and the specified value (exclusive). + * + * @param n Bound on the random number to be returned. Must be positive. + * @return a random long value between 0 (inclusive) and [n] (exclusive). + */ public override fun nextLong(n: Long): Long = generator.nextLong(n) } @@ -51,8 +117,14 @@ public fun RandomGenerator.asUniformRandomProvider(): UniformRandomProvider = if else RandomGeneratorProvider(this) +/** + * Returns [RandomSourceGenerator] with given [RandomSource] and [seed]. + */ public fun RandomGenerator.Companion.fromSource(source: RandomSource, seed: Long? = null): RandomSourceGenerator = RandomSourceGenerator(source, seed) +/** + * Returns [RandomSourceGenerator] with [RandomSource.MT] algorithm and given [seed]. + */ public fun RandomGenerator.Companion.mersenneTwister(seed: Long? = null): RandomSourceGenerator = fromSource(RandomSource.MT, seed) diff --git a/kmath-stat/src/jvmMain/kotlin/space/kscience/kmath/stat/distributions.kt b/kmath-stat/src/jvmMain/kotlin/space/kscience/kmath/stat/distributions.kt deleted file mode 100644 index 8b5551a16..000000000 --- a/kmath-stat/src/jvmMain/kotlin/space/kscience/kmath/stat/distributions.kt +++ /dev/null @@ -1,98 +0,0 @@ -package space.kscience.kmath.stat - -import org.apache.commons.rng.UniformRandomProvider -import org.apache.commons.rng.sampling.distribution.* -import space.kscience.kmath.chains.BlockingDoubleChain -import space.kscience.kmath.chains.BlockingIntChain -import space.kscience.kmath.chains.Chain -import kotlin.math.PI -import kotlin.math.exp -import kotlin.math.pow -import kotlin.math.sqrt - -public abstract class ContinuousSamplerDistribution : Distribution { - private inner class ContinuousSamplerChain(val generator: RandomGenerator) : BlockingDoubleChain() { - 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 - - public override fun sample(generator: RandomGenerator): BlockingDoubleChain = ContinuousSamplerChain(generator) -} - -public 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 - - public override fun sample(generator: RandomGenerator): BlockingIntChain = ContinuousSamplerChain(generator) -} - -public 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) - } - -public fun Distribution.Companion.normal( - method: NormalSamplerMethod = NormalSamplerMethod.Ziggurat, -): ContinuousSamplerDistribution = object : ContinuousSamplerDistribution() { - override fun buildCMSampler(generator: RandomGenerator): ContinuousSampler { - val provider = generator.asUniformRandomProvider() - return normalSampler(method, provider) - } - - override fun probability(arg: Double): Double = exp(-arg.pow(2) / 2) / sqrt(PI * 2) -} - -/** - * A univariate normal distribution with given [mean] and [sigma]. [method] defines commons-rng generation method - */ -public 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 = generator.asUniformRandomProvider() - val normalizedSampler = normalSampler(method, provider) - return GaussianSampler(normalizedSampler, mean, sigma) - } - - override fun probability(arg: Double): Double = exp(-(arg - mean).pow(2) / 2 / sigma2) / norm -} - -public fun Distribution.Companion.poisson(lambda: Double): DiscreteSamplerDistribution = - object : DiscreteSamplerDistribution() { - private val computedProb: MutableMap = hashMapOf(0 to exp(-lambda)) - - override fun buildSampler(generator: RandomGenerator): DiscreteSampler = - PoissonSampler.of(generator.asUniformRandomProvider(), 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-stat/src/jvmTest/kotlin/space/kscience/kmath/stat/CommonsDistributionsTest.kt b/kmath-stat/src/jvmTest/kotlin/space/kscience/kmath/stat/CommonsDistributionsTest.kt index 70708a5c8..76aac65c4 100644 --- a/kmath-stat/src/jvmTest/kotlin/space/kscience/kmath/stat/CommonsDistributionsTest.kt +++ b/kmath-stat/src/jvmTest/kotlin/space/kscience/kmath/stat/CommonsDistributionsTest.kt @@ -5,23 +5,22 @@ import kotlinx.coroutines.flow.toList import kotlinx.coroutines.runBlocking import org.junit.jupiter.api.Assertions import org.junit.jupiter.api.Test +import space.kscience.kmath.stat.samplers.GaussianSampler internal class CommonsDistributionsTest { @Test fun testNormalDistributionSuspend() { - val distribution = Distribution.normal(7.0, 2.0) + val distribution = GaussianSampler.of(7.0, 2.0) val generator = RandomGenerator.default(1) - val sample = runBlocking { - distribution.sample(generator).take(1000).toList() - } + 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 distribution = GaussianSampler.of(7.0, 2.0) val generator = RandomGenerator.default(1) - val sample = distribution.sample(generator).nextBlock(1000) + val sample = runBlocking { distribution.sample(generator).blocking().nextBlock(1000) } Assertions.assertEquals(7.0, sample.average(), 0.1) } } diff --git a/kmath-stat/src/jvmTest/kotlin/space/kscience/kmath/stat/SamplerTest.kt b/kmath-stat/src/jvmTest/kotlin/space/kscience/kmath/stat/SamplerTest.kt index 244b5107f..497a843c0 100644 --- a/kmath-stat/src/jvmTest/kotlin/space/kscience/kmath/stat/SamplerTest.kt +++ b/kmath-stat/src/jvmTest/kotlin/space/kscience/kmath/stat/SamplerTest.kt @@ -7,11 +7,8 @@ class SamplerTest { @Test fun bufferSamplerTest() { - val sampler: Sampler = - BasicSampler { it.chain { nextDouble() } } + val sampler = Sampler { it.chain { nextDouble() } } val data = sampler.sampleBuffer(RandomGenerator.default, 100) - runBlocking { - println(data.next()) - } + runBlocking { println(data.next()) } } } \ No newline at end of file diff --git a/settings.gradle.kts b/settings.gradle.kts index b4d7b3049..4467d5ed6 100644 --- a/settings.gradle.kts +++ b/settings.gradle.kts @@ -4,12 +4,11 @@ pluginManagement { mavenLocal() gradlePluginPortal() jcenter() - maven("https://dl.bintray.com/kotlin/kotlin-eap") maven("https://dl.bintray.com/kotlin/kotlinx") } - val toolsVersion = "0.9.1" - val kotlinVersion = "1.4.31" + val toolsVersion = "0.9.3" + val kotlinVersion = "1.4.32" plugins { id("kotlinx.benchmark") version "0.2.0-dev-20"