From 88d0c19a741f32145877c65d1cf7f70ff20b9ade Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Fri, 19 Mar 2021 11:07:27 +0300 Subject: [PATCH] Refactor structure features. Basic curve fitting --- build.gradle.kts | 6 +- .../kmath/commons/fit/fitWithAutoDiff.kt | 8 +- .../kscience/kmath/commons/linear/CMMatrix.kt | 3 +- ...timizationProblem.kt => CMOptimization.kt} | 42 +++--- .../kmath/commons/optimization/cmFit.kt | 22 ++-- .../commons/optimization/OptimizeTest.kt | 4 +- kmath-core/api/kmath-core.api | 9 +- .../kscience/kmath/linear/LinearSpace.kt | 4 +- .../kscience/kmath/linear/MatrixFeatures.kt | 4 +- .../kscience/kmath/linear/MatrixWrapper.kt | 3 +- .../space/kscience/kmath/nd/AlgebraND.kt | 5 +- .../space/kscience/kmath/nd/Structure1D.kt | 2 + .../space/kscience/kmath/nd/Structure2D.kt | 2 +- .../space/kscience/kmath/nd/StructureND.kt | 6 +- .../kmath/chains/BlockingDoubleChain.kt | 2 +- .../kscience/kmath/ejml/EjmlLinearSpace.kt | 3 +- .../kotlin/space/kscience/kmath/real/grids.kt | 74 ++++++----- kmath-stat/build.gradle.kts | 8 -- .../kscience/kmath/optimization/DataFit.kt | 17 +++ .../optimization/FunctionOptimization.kt | 122 ++++++++++++++++++ .../kmath/optimization/Optimization.kt | 44 +++++++ .../space/kscience/kmath/stat/Fitting.kt | 63 --------- .../kmath/stat/OptimizationProblem.kt | 88 ------------- .../space/kscience/kmath/stat/RandomChain.kt | 2 +- .../kscience/kmath/stat/distributions.kt | 25 ++-- 25 files changed, 311 insertions(+), 257 deletions(-) rename kmath-commons/src/main/kotlin/space/kscience/kmath/commons/optimization/{CMOptimizationProblem.kt => CMOptimization.kt} (75%) create mode 100644 kmath-stat/src/commonMain/kotlin/space/kscience/kmath/optimization/DataFit.kt create mode 100644 kmath-stat/src/commonMain/kotlin/space/kscience/kmath/optimization/FunctionOptimization.kt create mode 100644 kmath-stat/src/commonMain/kotlin/space/kscience/kmath/optimization/Optimization.kt delete mode 100644 kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/Fitting.kt delete mode 100644 kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/OptimizationProblem.kt diff --git a/build.gradle.kts b/build.gradle.kts index 5f1a8b88a..9570d7744 100644 --- a/build.gradle.kts +++ b/build.gradle.kts @@ -1,5 +1,3 @@ -import ru.mipt.npm.gradle.KSciencePublishingPlugin - plugins { id("ru.mipt.npm.gradle.project") } @@ -20,11 +18,11 @@ allprojects { } group = "space.kscience" - version = "0.3.0-dev-3" + version = "0.3.0-dev-4" } subprojects { - if (name.startsWith("kmath")) apply() + if (name.startsWith("kmath")) apply(plugin = "maven-publish") } readme { 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..6141a1058 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 @@ -8,10 +8,14 @@ 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.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.Distribution +import space.kscience.kmath.stat.RandomGenerator +import space.kscience.kmath.stat.normal import space.kscience.kmath.structures.asIterable import space.kscience.kmath.structures.toList import kotlin.math.pow @@ -58,7 +62,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) 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 75% 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..6200b61a9 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 @@ -10,21 +10,25 @@ 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.optimization.FunctionOptimization +import space.kscience.kmath.optimization.OptimizationFeature +import space.kscience.kmath.optimization.OptimizationProblemFactory +import space.kscience.kmath.optimization.OptimizationResult 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 { +public class CMOptimization( + override val symbols: List, +) : FunctionOptimization, 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 + ) public fun addOptimizationData(data: OptimizationData) { optimizationData[data::class] = data @@ -57,8 +61,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 +74,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 +88,7 @@ public class CMOptimizationProblem(override val symbols: List) : } public fun optimizer(block: () -> MultivariateOptimizer) { - optimizatorBuilder = block + optimizerBuilder = block } override fun update(result: OptimizationResult) { @@ -92,19 +96,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..384414e6d 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 @@ -6,16 +6,16 @@ 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.optimization.FunctionOptimization +import space.kscience.kmath.optimization.OptimizationResult +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,23 +43,23 @@ public fun Fitting.chiSquared( */ public fun Expression.optimize( vararg symbols: Symbol, - configuration: CMOptimizationProblem.() -> Unit, -): OptimizationResult = optimizeWith(CMOptimizationProblem, symbols = symbols, configuration) + configuration: CMOptimization.() -> Unit, +): OptimizationResult = optimizeWith(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) + val problem = CMOptimization(startPoint.map { it.first }).apply(configuration) problem.diffExpression(this) problem.initialGuess(startPoint.toMap()) problem.goal(GoalType.MINIMIZE) 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..cbbe1457e 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 @@ -3,8 +3,8 @@ package space.kscience.kmath.commons.optimization import org.junit.jupiter.api.Test import space.kscience.kmath.commons.expressions.DerivativeStructureExpression import space.kscience.kmath.expressions.symbol +import space.kscience.kmath.optimization.FunctionOptimization import space.kscience.kmath.stat.Distribution -import space.kscience.kmath.stat.Fitting import space.kscience.kmath.stat.RandomGenerator import space.kscience.kmath.stat.normal import kotlin.math.pow @@ -55,7 +55,7 @@ internal class OptimizeTest { 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 } diff --git a/kmath-core/api/kmath-core.api b/kmath-core/api/kmath-core.api index 04325379e..8c62198cd 100644 --- a/kmath-core/api/kmath-core.api +++ b/kmath-core/api/kmath-core.api @@ -575,7 +575,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 { @@ -1060,11 +1060,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 @@ -1104,6 +1108,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; 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 6a587270b..dfc4c7c9b 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 @@ -164,7 +164,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 { @@ -194,7 +194,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/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 8cf5d18cb..5483ed28f 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 @@ -15,6 +15,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 3834cd97c..5dfdd028a 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 @@ -69,7 +69,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-coroutines/src/commonMain/kotlin/space/kscience/kmath/chains/BlockingDoubleChain.kt b/kmath-coroutines/src/commonMain/kotlin/space/kscience/kmath/chains/BlockingDoubleChain.kt index bc4508994..ba6adf35b 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 @@ -8,5 +8,5 @@ public abstract class BlockingDoubleChain : Chain { override suspend fun next(): Double = nextDouble() - public fun nextBlock(size: Int): DoubleArray = DoubleArray(size) { nextDouble() } + public open fun nextBlock(size: Int): DoubleArray = DoubleArray(size) { nextDouble() } } 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 a82fe933e..f84a57c9d 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 @@ -89,7 +90,7 @@ public object EjmlLinearSpace : LinearSpace { v.toEjml().origin.scale(this).wrapVector() @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-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..446a3d9cb 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,43 @@ 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. + */ +@UnstableKMathAPI +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 / (numberOfPoints - 1) } +} + +/** + * 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. + */ +@UnstableKMathAPI +public fun Buffer.Companion.fromRange(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() + return DoubleBuffer(numberOfPoints) { normalizedRange.start + step * it / (numberOfPoints - 1) } +} /** * Convert double range to sequence. @@ -11,35 +47,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.fromRange(this, step) \ No newline at end of file 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/DataFit.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/optimization/DataFit.kt new file mode 100644 index 000000000..e9a5dea86 --- /dev/null +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/optimization/DataFit.kt @@ -0,0 +1,17 @@ +package space.kscience.kmath.optimization + +import space.kscience.kmath.expressions.DifferentiableExpression +import space.kscience.kmath.expressions.StringSymbol +import space.kscience.kmath.expressions.Symbol +import space.kscience.kmath.structures.Buffer + +public interface DataFit : Optimization { + + public fun modelAndData( + x: Buffer, + y: Buffer, + yErr: Buffer, + model: DifferentiableExpression, + xSymbol: Symbol = StringSymbol("x"), + ) +} \ No newline at end of file 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..bf37f5c64 --- /dev/null +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/optimization/FunctionOptimization.kt @@ -0,0 +1,122 @@ +package space.kscience.kmath.optimization + +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 + +/** + * A likelihood function optimization problem + */ +public interface FunctionOptimization: Optimization, DataFit { + /** + * 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>) + + override fun modelAndData( + x: Buffer, + y: Buffer, + yErr: Buffer, + model: DifferentiableExpression, + xSymbol: Symbol, + ) { + 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" } + + } + + 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 + } + } + + /** + * 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) + } + } + } + } +} + +/** + * 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/optimization/Optimization.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/optimization/Optimization.kt new file mode 100644 index 000000000..370274b41 --- /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.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) + +/** + * 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/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..978094ffd 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 @@ -14,4 +14,4 @@ public class RandomChain( override fun fork(): Chain = RandomChain(generator.fork(), gen) } -public fun RandomGenerator.chain(gen: suspend RandomGenerator.() -> R): RandomChain = RandomChain(this, gen) +public fun RandomGenerator.chain(gen: suspend RandomGenerator.() -> R): RandomChain = RandomChain(this, gen) \ No newline at end of file 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 index 8b5551a16..c3d711789 100644 --- a/kmath-stat/src/jvmMain/kotlin/space/kscience/kmath/stat/distributions.kt +++ b/kmath-stat/src/jvmMain/kotlin/space/kscience/kmath/stat/distributions.kt @@ -80,19 +80,20 @@ public fun Distribution.Companion.normal( 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)) +public fun Distribution.Companion.poisson( + lambda: Double, +): DiscreteSamplerDistribution = object : DiscreteSamplerDistribution() { + private val computedProb: HashMap = hashMapOf(0 to exp(-lambda)) - override fun buildSampler(generator: RandomGenerator): DiscreteSampler = - PoissonSampler.of(generator.asUniformRandomProvider(), 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" } + 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 } - } + return if (arg > 40) + exp(-(arg - lambda).pow(2) / 2 / lambda) / sqrt(2 * PI * lambda) + else + computedProb.getOrPut(arg) { probability(arg - 1) * lambda / arg } } +}