diff --git a/examples/build.gradle.kts b/examples/build.gradle.kts index 4cc6fecc0..d06005321 100644 --- a/examples/build.gradle.kts +++ b/examples/build.gradle.kts @@ -20,6 +20,7 @@ dependencies { implementation(project(":kmath-coroutines")) implementation(project(":kmath-commons")) implementation(project(":kmath-complex")) + implementation(project(":kmath-optimization")) implementation(project(":kmath-stat")) implementation(project(":kmath-viktor")) implementation(project(":kmath-dimensions")) diff --git a/examples/src/main/kotlin/space/kscience/kmath/commons/fit/fitWithAutoDiff.kt b/examples/src/main/kotlin/space/kscience/kmath/fit/chiSquared.kt similarity index 91% rename from examples/src/main/kotlin/space/kscience/kmath/commons/fit/fitWithAutoDiff.kt rename to examples/src/main/kotlin/space/kscience/kmath/fit/chiSquared.kt index 8d95ebb4a..c1070de52 100644 --- a/examples/src/main/kotlin/space/kscience/kmath/commons/fit/fitWithAutoDiff.kt +++ b/examples/src/main/kotlin/space/kscience/kmath/fit/chiSquared.kt @@ -3,19 +3,17 @@ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file. */ -package space.kscience.kmath.commons.fit +package space.kscience.kmath.fit import kotlinx.html.br import kotlinx.html.h3 import space.kscience.kmath.commons.expressions.DSProcessor import space.kscience.kmath.commons.optimization.CMOptimizer import space.kscience.kmath.distributions.NormalDistribution +import space.kscience.kmath.expressions.binding import space.kscience.kmath.expressions.chiSquaredExpression import space.kscience.kmath.expressions.symbol -import space.kscience.kmath.optimization.FunctionOptimizationTarget -import space.kscience.kmath.optimization.optimizeWith -import space.kscience.kmath.optimization.resultPoint -import space.kscience.kmath.optimization.resultValue +import space.kscience.kmath.optimization.* import space.kscience.kmath.real.DoubleVector import space.kscience.kmath.real.map import space.kscience.kmath.real.step @@ -25,6 +23,7 @@ import space.kscience.kmath.structures.toList import space.kscience.plotly.* import space.kscience.plotly.models.ScatterMode import space.kscience.plotly.models.TraceValues +import kotlin.math.abs import kotlin.math.pow import kotlin.math.sqrt @@ -45,7 +44,7 @@ operator fun TraceValues.invoke(vector: DoubleVector) { */ suspend fun main() { //A generator for a normally distributed values - val generator = NormalDistribution(2.0, 7.0) + val generator = NormalDistribution(0.0, 1.0) //A chain/flow of random values with the given seed val chain = generator.sample(RandomGenerator.default(112667)) @@ -56,7 +55,7 @@ suspend fun main() { //Perform an operation on each x value (much more effective, than numpy) - val y = x.map { + val y = x.map { it -> val value = it.pow(2) + it + 1 value + chain.next() * sqrt(value) } diff --git a/examples/src/main/kotlin/space/kscience/kmath/fit/qowFit.kt b/examples/src/main/kotlin/space/kscience/kmath/fit/qowFit.kt new file mode 100644 index 000000000..944f80697 --- /dev/null +++ b/examples/src/main/kotlin/space/kscience/kmath/fit/qowFit.kt @@ -0,0 +1,105 @@ +/* + * Copyright 2018-2021 KMath contributors. + * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file. + */ + +package space.kscience.kmath.fit + +import kotlinx.html.br +import kotlinx.html.h3 +import space.kscience.kmath.commons.expressions.DSProcessor +import space.kscience.kmath.data.XYErrorColumnarData +import space.kscience.kmath.distributions.NormalDistribution +import space.kscience.kmath.expressions.Symbol +import space.kscience.kmath.expressions.binding +import space.kscience.kmath.expressions.symbol +import space.kscience.kmath.optimization.QowOptimizer +import space.kscience.kmath.optimization.fitWith +import space.kscience.kmath.optimization.resultPoint +import space.kscience.kmath.real.map +import space.kscience.kmath.real.step +import space.kscience.kmath.stat.RandomGenerator +import space.kscience.kmath.structures.asIterable +import space.kscience.kmath.structures.toList +import space.kscience.plotly.* +import space.kscience.plotly.models.ScatterMode +import kotlin.math.abs +import kotlin.math.pow +import kotlin.math.sqrt + +// Forward declaration of symbols that will be used in expressions. +private val a by symbol +private val b by symbol +private val c by symbol + + +/** + * Least squares fie with auto-differentiation. Uses `kmath-commons` and `kmath-for-real` modules. + */ +suspend fun main() { + //A generator for a normally distributed values + val generator = NormalDistribution(0.0, 1.0) + + //A chain/flow of random values with the given seed + val chain = generator.sample(RandomGenerator.default(112667)) + + + //Create a uniformly distributed x values like numpy.arrange + val x = 1.0..100.0 step 1.0 + + + //Perform an operation on each x value (much more effective, than numpy) + val y = x.map { it -> + val value = it.pow(2) + it + 100 + value + chain.next() * sqrt(value) + } + // this will also work, but less effective: + // val y = x.pow(2)+ x + 1 + chain.nextDouble() + + // create same errors for all xs + val yErr = y.map { sqrt(abs(it)) } + require(yErr.asIterable().all { it > 0 }) { "All errors must be strictly positive" } + + val result = XYErrorColumnarData.of(x, y, yErr).fitWith( + QowOptimizer, + DSProcessor, + mapOf(a to 1.0, b to 1.2, c to 99.0) + ) { arg -> + //bind variables to autodiff context + val a by binding + val b by binding + //Include default value for c if it is not provided as a parameter + val c = bindSymbolOrNull(c) ?: one + a * arg.pow(2) + b * arg + c + } + + //display a page with plot and numerical results + val page = Plotly.page { + plot { + scatter { + mode = ScatterMode.markers + x(x) + y(y) + error_y { + array = yErr.toList() + } + name = "data" + } + scatter { + mode = ScatterMode.lines + x(x) + y(x.map { result.model(result.resultPoint + (Symbol.x to it)) }) + name = "fit" + } + } + br() + h3 { + +"Fit result: ${result.resultPoint}" + } +// h3 { +// +"Chi2/dof = ${result.resultValue / (x.size - 3)}" +// } + } + + page.makeFile() +} diff --git a/examples/src/main/kotlin/space/kscience/kmath/functions/matrixIntegration.kt b/examples/src/main/kotlin/space/kscience/kmath/functions/matrixIntegration.kt index 2619d3d74..e5ba29d73 100644 --- a/examples/src/main/kotlin/space/kscience/kmath/functions/matrixIntegration.kt +++ b/examples/src/main/kotlin/space/kscience/kmath/functions/matrixIntegration.kt @@ -22,7 +22,7 @@ fun main(): Unit = DoubleField { } //Define a function in a nd space - val function: (Double) -> StructureND = { x: Double -> 3 * number(x).pow(2) + 2 * diagonal(x) + 1 } + val function: (Double) -> StructureND = { x: Double -> 3 * x.pow(2) + 2 * diagonal(x) + 1 } //get the result of the integration val result = gaussIntegrator.integrate(0.0..10.0, function = function) diff --git a/kmath-commons/build.gradle.kts b/kmath-commons/build.gradle.kts index a208c956c..96c17a215 100644 --- a/kmath-commons/build.gradle.kts +++ b/kmath-commons/build.gradle.kts @@ -9,6 +9,7 @@ dependencies { api(project(":kmath-core")) api(project(":kmath-complex")) api(project(":kmath-coroutines")) + api(project(":kmath-optimization")) api(project(":kmath-stat")) api(project(":kmath-functions")) api("org.apache.commons:commons-math3:3.6.1") diff --git a/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/optimization/CMOptimizer.kt b/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/optimization/CMOptimizer.kt index abf95daf6..11eb6fba8 100644 --- a/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/optimization/CMOptimizer.kt +++ b/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/optimization/CMOptimizer.kt @@ -18,6 +18,7 @@ import space.kscience.kmath.expressions.SymbolIndexer import space.kscience.kmath.expressions.derivative import space.kscience.kmath.expressions.withSymbols import space.kscience.kmath.misc.UnstableKMathAPI +import space.kscience.kmath.misc.log import space.kscience.kmath.optimization.* import kotlin.collections.set import kotlin.reflect.KClass @@ -108,15 +109,17 @@ public object CMOptimizer : Optimizer> { val objectiveFunction = ObjectiveFunction { val args = startPoint + it.toMap() - problem.expression(args) + val res = problem.expression(args) + res } addOptimizationData(objectiveFunction) val gradientFunction = ObjectiveFunctionGradient { val args = startPoint + it.toMap() - DoubleArray(symbols.size) { index -> + val res = DoubleArray(symbols.size) { index -> problem.expression.derivative(symbols[index])(args) } + res } addOptimizationData(gradientFunction) 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 3c57f5467..56252ab34 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 @@ -42,8 +42,8 @@ internal class AutoDiffTest { @Test fun autoDifTest() { val f = DerivativeStructureExpression { - val x by binding() - val y by binding() + val x by binding + val y by binding x.pow(2) + 2 * x * y + y.pow(2) + 1 } 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 170fceceb..facbf3a9a 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 @@ -15,6 +15,7 @@ import space.kscience.kmath.expressions.chiSquaredExpression import space.kscience.kmath.expressions.symbol import space.kscience.kmath.optimization.* import space.kscience.kmath.stat.RandomGenerator +import space.kscience.kmath.structures.DoubleBuffer import space.kscience.kmath.structures.asBuffer import space.kscience.kmath.structures.map import kotlin.math.pow @@ -58,7 +59,7 @@ internal class OptimizeTest { it.pow(2) + it + 1 + chain.next() } - val yErr = List(x.size) { sigma }.asBuffer() + val yErr = DoubleBuffer(x.size) { sigma } val chi2 = DSProcessor.chiSquaredExpression( x, y, yErr 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 index bd0e13b9d..2fce772cc 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/data/XYColumnarData.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/data/XYColumnarData.kt @@ -32,19 +32,21 @@ public interface XYColumnarData : ColumnarData { Symbol.y -> y else -> null } -} -@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 + public companion object{ + @UnstableKMathAPI + public fun of(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 + } + } } } + /** * Represent a [ColumnarData] as an [XYColumnarData]. The presence or respective columns is checked on creation. */ diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/data/XYErrorColumnarData.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/data/XYErrorColumnarData.kt index 7199de888..8ddd6406f 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/data/XYErrorColumnarData.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/data/XYErrorColumnarData.kt @@ -5,15 +5,13 @@ package space.kscience.kmath.data -import space.kscience.kmath.data.XYErrorColumnarData.Companion import space.kscience.kmath.expressions.Symbol -import space.kscience.kmath.expressions.symbol import space.kscience.kmath.misc.UnstableKMathAPI import space.kscience.kmath.structures.Buffer /** - * A [ColumnarData] with additional [Companion.yErr] column for an [Symbol.y] error + * A [ColumnarData] with additional [Symbol.yError] column for an [Symbol.y] error * Inherits [XYColumnarData]. */ @UnstableKMathAPI @@ -23,11 +21,24 @@ public interface XYErrorColumnarData : XYColumnarData = when (symbol) { Symbol.x -> x Symbol.y -> y - Companion.yErr -> yErr + Symbol.yError -> yErr else -> error("A column for symbol $symbol not found") } - public companion object{ - public val yErr: Symbol by symbol + public companion object { + public fun of( + x: Buffer, y: Buffer, yErr: Buffer + ): XYErrorColumnarData { + require(x.size == y.size) { "Buffer size mismatch. x buffer size is ${x.size}, y buffer size is ${y.size}" } + require(y.size == yErr.size) { "Buffer size mismatch. y buffer size is ${x.size}, yErr buffer size is ${y.size}" } + + return object : XYErrorColumnarData { + override val size: Int = x.size + override val x: Buffer = x + override val y: Buffer = y + override val yErr: Buffer = yErr + } + } } -} \ No newline at end of file +} + 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 5105c2bec..e94cb98eb 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 @@ -68,6 +68,7 @@ public interface ExpressionAlgebra : Algebra { /** * Bind a symbol by name inside the [ExpressionAlgebra] */ -public fun ExpressionAlgebra.binding(): ReadOnlyProperty = ReadOnlyProperty { _, property -> - bindSymbol(property.name) ?: error("A variable with name ${property.name} does not exist") -} +public val ExpressionAlgebra.binding: ReadOnlyProperty + get() = ReadOnlyProperty { _, property -> + bindSymbol(property.name) ?: error("A variable with name ${property.name} does not exist") + } diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/specialExpressions.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/specialExpressions.kt index ede4a779c..6b17dfca5 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/specialExpressions.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/specialExpressions.kt @@ -7,13 +7,18 @@ package space.kscience.kmath.expressions import space.kscience.kmath.operations.ExtendedField import space.kscience.kmath.structures.Buffer +import space.kscience.kmath.structures.asIterable import space.kscience.kmath.structures.indices +import kotlin.jvm.JvmName /** * Generate a chi squared expression from given x-y-sigma data and inline model. Provides automatic * differentiation. + * + * **WARNING** All elements of [yErr] must be positive. */ -public fun AutoDiffProcessor.chiSquaredExpression( +@JvmName("genericChiSquaredExpression") +public fun , I : Any, A> AutoDiffProcessor.chiSquaredExpression( x: Buffer, y: Buffer, yErr: Buffer, @@ -35,4 +40,14 @@ public fun AutoDiffProcessor.chiSquaredExpression sum } +} + +public fun AutoDiffProcessor.chiSquaredExpression( + x: Buffer, + y: Buffer, + yErr: Buffer, + model: A.(I) -> I, +): DifferentiableExpression where A : ExtendedField, A : ExpressionAlgebra { + require(yErr.asIterable().all { it > 0.0 }) { "All errors must be strictly positive" } + return chiSquaredExpression(x, y, yErr, model) } \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/misc/Featured.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/misc/Featured.kt index be1c8380c..29b7caec6 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/misc/Featured.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/misc/Featured.kt @@ -5,6 +5,7 @@ package space.kscience.kmath.misc +import kotlin.jvm.JvmInline import kotlin.reflect.KClass /** @@ -29,7 +30,8 @@ public interface Feature> { /** * A container for a set of features */ -public class FeatureSet> private constructor(public val features: Map, F>) : Featured { +@JvmInline +public value class FeatureSet> private constructor(public val features: Map, F>) : Featured { @Suppress("UNCHECKED_CAST") override fun getFeature(type: FeatureKey): T? = features[type]?.let { it as T } @@ -48,6 +50,9 @@ public class FeatureSet> private constructor(public val features: public operator fun iterator(): Iterator = features.values.iterator() + override fun toString(): String = features.values.joinToString(prefix = "[ ", postfix = " ]") + + public companion object { public fun > of(vararg features: F): FeatureSet = FeatureSet(features.associateBy { it.key }) public fun > of(features: Iterable): FeatureSet = diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/misc/logging.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/misc/logging.kt index d13840841..9dfc564c3 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/misc/logging.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/misc/logging.kt @@ -5,10 +5,18 @@ package space.kscience.kmath.misc -public interface Loggable { - public fun log(tag: String = INFO, block: () -> String) +import space.kscience.kmath.misc.Loggable.Companion.INFO + +public fun interface Loggable { + public fun log(tag: String, block: () -> String) public companion object { public const val INFO: String = "INFO" + + public val console: Loggable = Loggable { tag, block -> + println("[$tag] ${block()}") + } } -} \ No newline at end of file +} + +public fun Loggable.log(block: () -> String): Unit = log(INFO, block) \ No newline at end of file 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 4d1b00b3d..80c5943cf 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 @@ -16,7 +16,7 @@ class ExpressionFieldTest { @Test fun testExpression() { val expression = with(FunctionalExpressionField(DoubleField)) { - val x by binding() + val x by binding x * x + 2 * x + one } @@ -27,7 +27,7 @@ class ExpressionFieldTest { @Test fun separateContext() { fun FunctionalExpressionField.expression(): Expression { - val x by binding() + val x by binding return x * x + 2 * x + one } @@ -38,7 +38,7 @@ class ExpressionFieldTest { @Test fun valueExpression() { val expressionBuilder: FunctionalExpressionField.() -> Expression = { - val x by binding() + val x by binding x * x + 2 * x + one } diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/GaussIntegrator.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/GaussIntegrator.kt index 9b938394a..9785d7744 100644 --- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/GaussIntegrator.kt +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/GaussIntegrator.kt @@ -73,7 +73,8 @@ public class GaussIntegrator( } /** - * Create a Gauss-Legendre integrator for this field. + * Create a Gauss integrator for this field. By default, uses Legendre rule to compute points and weights. + * Custom rules could be provided by [GaussIntegratorRuleFactory] feature. * @see [GaussIntegrator] */ public val Field.gaussIntegrator: GaussIntegrator get() = GaussIntegrator(this) 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 f4a0abd5d..b13adefa5 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 @@ -42,20 +42,20 @@ public fun > PolynomialInterpolator.interpolatePolynomials( x: Buffer, y: Buffer, ): PiecewisePolynomial { - val pointSet = XYColumnarData(x, y) + val pointSet = XYColumnarData.of(x, y) return interpolatePolynomials(pointSet) } public fun > PolynomialInterpolator.interpolatePolynomials( data: Map, ): PiecewisePolynomial { - val pointSet = XYColumnarData(data.keys.toList().asBuffer(), data.values.toList().asBuffer()) + val pointSet = XYColumnarData.of(data.keys.toList().asBuffer(), data.values.toList().asBuffer()) return interpolatePolynomials(pointSet) } public fun > PolynomialInterpolator.interpolatePolynomials( data: List>, ): PiecewisePolynomial { - val pointSet = XYColumnarData(data.map { it.first }.asBuffer(), data.map { it.second }.asBuffer()) + val pointSet = XYColumnarData.of(data.map { it.first }.asBuffer(), data.map { it.second }.asBuffer()) return interpolatePolynomials(pointSet) } diff --git a/kmath-optimization/build.gradle.kts b/kmath-optimization/build.gradle.kts new file mode 100644 index 000000000..ca013b2f4 --- /dev/null +++ b/kmath-optimization/build.gradle.kts @@ -0,0 +1,20 @@ +plugins { + id("ru.mipt.npm.gradle.mpp") + id("ru.mipt.npm.gradle.native") +} + +kscience { + useAtomic() +} + +kotlin.sourceSets { + commonMain { + dependencies { + api(project(":kmath-coroutines")) + } + } +} + +readme { + maturity = ru.mipt.npm.gradle.Maturity.EXPERIMENTAL +} \ No newline at end of file diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/optimization/FunctionOptimization.kt b/kmath-optimization/src/commonMain/kotlin/space/kscience/kmath/optimization/FunctionOptimization.kt similarity index 65% rename from kmath-stat/src/commonMain/kotlin/space/kscience/kmath/optimization/FunctionOptimization.kt rename to kmath-optimization/src/commonMain/kotlin/space/kscience/kmath/optimization/FunctionOptimization.kt index 1af6c4bda..02602b068 100644 --- a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/optimization/FunctionOptimization.kt +++ b/kmath-optimization/src/commonMain/kotlin/space/kscience/kmath/optimization/FunctionOptimization.kt @@ -5,13 +5,11 @@ package space.kscience.kmath.optimization -import space.kscience.kmath.expressions.* +import space.kscience.kmath.expressions.DifferentiableExpression +import space.kscience.kmath.expressions.Symbol import space.kscience.kmath.misc.FeatureSet -import space.kscience.kmath.operations.ExtendedField -import space.kscience.kmath.structures.Buffer -import space.kscience.kmath.structures.indices -public class OptimizationValue(public val value: T) : OptimizationFeature{ +public class OptimizationValue(public val value: T) : OptimizationFeature { override fun toString(): String = "Value($value)" } @@ -23,9 +21,28 @@ public enum class FunctionOptimizationTarget : OptimizationFeature { public class FunctionOptimization( override val features: FeatureSet, public val expression: DifferentiableExpression, -) : OptimizationProblem{ +) : OptimizationProblem { - public companion object + + override fun equals(other: Any?): Boolean { + if (this === other) return true + if (other == null || this::class != other::class) return false + + other as FunctionOptimization<*> + + if (features != other.features) return false + if (expression != other.expression) return false + + return true + } + + override fun hashCode(): Int { + var result = features.hashCode() + result = 31 * result + expression.hashCode() + return result + } + + override fun toString(): String = "FunctionOptimization(features=$features)" } public fun FunctionOptimization.withFeatures( @@ -47,7 +64,7 @@ public suspend fun DifferentiableExpression.optimizeWith( return optimizer.optimize(problem) } -public val FunctionOptimization.resultValueOrNull:T? +public val FunctionOptimization.resultValueOrNull: T? get() = getFeature>()?.point?.let { expression(it) } public val FunctionOptimization.resultValue: T diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/optimization/OptimizationBuilder.kt b/kmath-optimization/src/commonMain/kotlin/space/kscience/kmath/optimization/OptimizationBuilder.kt similarity index 94% rename from kmath-stat/src/commonMain/kotlin/space/kscience/kmath/optimization/OptimizationBuilder.kt rename to kmath-optimization/src/commonMain/kotlin/space/kscience/kmath/optimization/OptimizationBuilder.kt index 7d52ae26e..10a25c029 100644 --- a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/optimization/OptimizationBuilder.kt +++ b/kmath-optimization/src/commonMain/kotlin/space/kscience/kmath/optimization/OptimizationBuilder.kt @@ -72,15 +72,15 @@ public suspend fun DifferentiableExpression.optimizeWith( public class XYOptimizationBuilder( public val data: XYColumnarData, public val model: DifferentiableExpression, -) : OptimizationBuilder() { +) : OptimizationBuilder() { public var pointToCurveDistance: PointToCurveDistance = PointToCurveDistance.byY public var pointWeight: PointWeight = PointWeight.byYSigma - override fun build(): XYOptimization = XYOptimization( - FeatureSet.of(features), + override fun build(): XYFit = XYFit( data, model, + FeatureSet.of(features), pointToCurveDistance, pointWeight ) @@ -90,4 +90,4 @@ public fun XYOptimization( data: XYColumnarData, model: DifferentiableExpression, builder: XYOptimizationBuilder.() -> Unit, -): XYOptimization = XYOptimizationBuilder(data, model).apply(builder).build() \ No newline at end of file +): XYFit = XYOptimizationBuilder(data, model).apply(builder).build() \ No newline at end of file diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/optimization/OptimizationProblem.kt b/kmath-optimization/src/commonMain/kotlin/space/kscience/kmath/optimization/OptimizationProblem.kt similarity index 100% rename from kmath-stat/src/commonMain/kotlin/space/kscience/kmath/optimization/OptimizationProblem.kt rename to kmath-optimization/src/commonMain/kotlin/space/kscience/kmath/optimization/OptimizationProblem.kt diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/optimization/Optimizer.kt b/kmath-optimization/src/commonMain/kotlin/space/kscience/kmath/optimization/Optimizer.kt similarity index 100% rename from kmath-stat/src/commonMain/kotlin/space/kscience/kmath/optimization/Optimizer.kt rename to kmath-optimization/src/commonMain/kotlin/space/kscience/kmath/optimization/Optimizer.kt diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/optimization/QowOptimizer.kt b/kmath-optimization/src/commonMain/kotlin/space/kscience/kmath/optimization/QowOptimizer.kt similarity index 77% rename from kmath-stat/src/commonMain/kotlin/space/kscience/kmath/optimization/QowOptimizer.kt rename to kmath-optimization/src/commonMain/kotlin/space/kscience/kmath/optimization/QowOptimizer.kt index 365c58952..f0969d2f6 100644 --- a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/optimization/QowOptimizer.kt +++ b/kmath-optimization/src/commonMain/kotlin/space/kscience/kmath/optimization/QowOptimizer.kt @@ -11,6 +11,7 @@ import space.kscience.kmath.expressions.SymbolIndexer import space.kscience.kmath.expressions.derivative import space.kscience.kmath.linear.* import space.kscience.kmath.misc.UnstableKMathAPI +import space.kscience.kmath.misc.log import space.kscience.kmath.operations.DoubleField import space.kscience.kmath.structures.DoubleBuffer import space.kscience.kmath.structures.DoubleL2Norm @@ -21,14 +22,14 @@ import space.kscience.kmath.structures.DoubleL2Norm * See [the article](http://arxiv.org/abs/physics/0604127). */ @UnstableKMathAPI -public class QowOptimizer : Optimizer { +public object QowOptimizer : Optimizer { private val linearSpace: LinearSpace = LinearSpace.double private val solver: LinearSolver = linearSpace.lupSolver() @OptIn(UnstableKMathAPI::class) - private inner class QoWeight( - val problem: XYOptimization, + private class QoWeight( + val problem: XYFit, val parameters: Map, ) : Map by parameters, SymbolIndexer { override val symbols: List = parameters.keys.toList() @@ -39,8 +40,8 @@ public class QowOptimizer : Optimizer { * Derivatives of the spectrum over parameters. First index in the point number, second one - index of parameter */ val derivs: Matrix by lazy { - linearSpace.buildMatrix(problem.data.size, symbols.size) { i, k -> - problem.distance(i).derivative(symbols[k])(parameters) + linearSpace.buildMatrix(problem.data.size, symbols.size) { d, s -> + problem.distance(d).derivative(symbols[s])(parameters) } } @@ -48,25 +49,27 @@ public class QowOptimizer : Optimizer { * Array of dispersions in each point */ val dispersion: Point by lazy { - DoubleBuffer(problem.data.size) { i -> - problem.weight(i).invoke(parameters) + DoubleBuffer(problem.data.size) { d -> + problem.weight(d).invoke(parameters) } } val prior: DifferentiableExpression? get() = problem.getFeature>() + + override fun toString(): String = parameters.toString() } /** - * The signed distance from the model to the [i]-th point of data. + * The signed distance from the model to the [d]-th point of data. */ - private fun QoWeight.distance(i: Int, parameters: Map): Double = problem.distance(i)(parameters) + private fun QoWeight.distance(d: Int, parameters: Map): Double = problem.distance(d)(parameters) /** * The derivative of [distance] */ - private fun QoWeight.distanceDerivative(symbol: Symbol, i: Int, parameters: Map): Double = - problem.distance(i).derivative(symbol)(parameters) + private fun QoWeight.distanceDerivative(symbol: Symbol, d: Int, parameters: Map): Double = + problem.distance(d).derivative(symbol)(parameters) /** * Теоретическая ковариация весовых функций. @@ -74,8 +77,8 @@ public class QowOptimizer : Optimizer { * D(\phi)=E(\phi_k(\theta_0) \phi_l(\theta_0))= disDeriv_k * disDeriv_l /sigma^2 */ private fun QoWeight.covarF(): Matrix = - linearSpace.matrix(size, size).symmetric { k, l -> - (0 until data.size).sumOf { i -> derivs[k, i] * derivs[l, i] / dispersion[i] } + linearSpace.matrix(size, size).symmetric { s1, s2 -> + (0 until data.size).sumOf { d -> derivs[d, s1] * derivs[d, s2] / dispersion[d] } } /** @@ -89,12 +92,12 @@ public class QowOptimizer : Optimizer { * количество вызывов функции будет dim^2 вместо dim Первый индекс - * номер точки, второй - номер переменной, по которой берется производная */ - val eqvalues = linearSpace.buildMatrix(data.size, size) { i, l -> - distance(i, theta) * derivs[l, i] / dispersion[i] + val eqvalues = linearSpace.buildMatrix(data.size, size) { d, s -> + distance(d, theta) * derivs[d, s] / dispersion[d] } - buildMatrix(size, size) { k, l -> - (0 until data.size).sumOf { i -> eqvalues[i, l] * eqvalues[i, k] } + buildMatrix(size, size) { s1, s2 -> + (0 until data.size).sumOf { d -> eqvalues[d, s2] * eqvalues[d, s1] } } } @@ -106,20 +109,20 @@ public class QowOptimizer : Optimizer { ): Matrix = with(linearSpace) { //Возвращает производную k-того Eq по l-тому параметру //val res = Array(fitDim) { DoubleArray(fitDim) } - val sderiv = buildMatrix(data.size, size) { i, l -> - distanceDerivative(symbols[l], i, theta) + val sderiv = buildMatrix(data.size, size) { d, s -> + distanceDerivative(symbols[s], d, theta) } - buildMatrix(size, size) { k, l -> - val base = (0 until data.size).sumOf { i -> - require(dispersion[i] > 0) - sderiv[i, l] * derivs[k, i] / dispersion[i] + buildMatrix(size, size) { s1, s2 -> + val base = (0 until data.size).sumOf { d -> + require(dispersion[d] > 0) + sderiv[d, s2] * derivs[d, s1] / dispersion[d] } prior?.let { prior -> //Check if this one is correct val pi = prior(theta) - val deriv1 = prior.derivative(symbols[k])(theta) - val deriv2 = prior.derivative(symbols[l])(theta) + val deriv1 = prior.derivative(symbols[s1])(theta) + val deriv2 = prior.derivative(symbols[s2])(theta) base + deriv1 * deriv2 / pi / pi } ?: base } @@ -130,13 +133,13 @@ public class QowOptimizer : Optimizer { * Значения уравнений метода квазиоптимальных весов */ private fun QoWeight.getEqValues(theta: Map = this): Point { - val distances = DoubleBuffer(data.size) { i -> distance(i, theta) } + val distances = DoubleBuffer(data.size) { d -> distance(d, theta) } - return DoubleBuffer(size) { k -> - val base = (0 until data.size).sumOf { i -> distances[i] * derivs[k, i] / dispersion[i] } + return DoubleBuffer(size) { s -> + val base = (0 until data.size).sumOf { d -> distances[d] * derivs[d, s] / dispersion[d] } //Поправка на априорную вероятность prior?.let { prior -> - base - prior.derivative(symbols[k])(theta) / prior(theta) + base - prior.derivative(symbols[s])(theta) / prior(theta) } ?: base } } @@ -163,15 +166,15 @@ public class QowOptimizer : Optimizer { val logger = problem.getFeature() - var dis: Double//норма невязки - // Для удобства работаем всегда с полным набором параметров + var dis: Double //discrepancy value + // Working with the full set of parameters var par = problem.startPoint logger?.log { "Starting newtonian iteration from: \n\t$par" } - var eqvalues = getEqValues(par)//значения функций + var eqvalues = getEqValues(par) //Values of the weight functions - dis = DoubleL2Norm.norm(eqvalues)// невязка + dis = DoubleL2Norm.norm(eqvalues) // discrepancy logger?.log { "Starting discrepancy is $dis" } var i = 0 var flag = false @@ -238,7 +241,8 @@ public class QowOptimizer : Optimizer { return covar } - override suspend fun optimize(problem: XYOptimization): XYOptimization { + override suspend fun optimize(problem: XYFit): XYFit { + val qowSteps = 2 val initialWeight = QoWeight(problem, problem.startPoint) val res = initialWeight.newtonianRun() return res.problem.withFeature(OptimizationResult(res.parameters)) diff --git a/kmath-optimization/src/commonMain/kotlin/space/kscience/kmath/optimization/XYFit.kt b/kmath-optimization/src/commonMain/kotlin/space/kscience/kmath/optimization/XYFit.kt new file mode 100644 index 000000000..cbd765923 --- /dev/null +++ b/kmath-optimization/src/commonMain/kotlin/space/kscience/kmath/optimization/XYFit.kt @@ -0,0 +1,125 @@ +/* + * Copyright 2018-2021 KMath contributors. + * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file. + */ +@file:OptIn(UnstableKMathAPI::class) + +package space.kscience.kmath.optimization + +import space.kscience.kmath.data.XYColumnarData +import space.kscience.kmath.expressions.* +import space.kscience.kmath.misc.FeatureSet +import space.kscience.kmath.misc.Loggable +import space.kscience.kmath.misc.UnstableKMathAPI +import space.kscience.kmath.operations.ExtendedField +import space.kscience.kmath.operations.bindSymbol +import kotlin.math.pow + +/** + * Specify the way to compute distance from point to the curve as DifferentiableExpression + */ +public interface PointToCurveDistance : OptimizationFeature { + public fun distance(problem: XYFit, index: Int): DifferentiableExpression + + public companion object { + public val byY: PointToCurveDistance = object : PointToCurveDistance { + override fun distance(problem: XYFit, index: Int): DifferentiableExpression { + val x = problem.data.x[index] + val y = problem.data.y[index] + + return object : DifferentiableExpression { + override fun derivativeOrNull( + symbols: List + ): Expression? = problem.model.derivativeOrNull(symbols)?.let { derivExpression -> + Expression { arguments -> + derivExpression.invoke(arguments + (Symbol.x to x)) + } + } + + override fun invoke(arguments: Map): Double = + problem.model(arguments + (Symbol.x to x)) - y + } + } + + override fun toString(): String = "PointToCurveDistanceByY" + } + } +} + +/** + * Compute a wight of the point. The more the weight, the more impact this point will have on the fit. + * By default, uses Dispersion^-1 + */ +public interface PointWeight : OptimizationFeature { + public fun weight(problem: XYFit, index: Int): DifferentiableExpression + + public companion object { + public fun bySigma(sigmaSymbol: Symbol): PointWeight = object : PointWeight { + override fun weight(problem: XYFit, index: Int): DifferentiableExpression = + object : DifferentiableExpression { + override fun invoke(arguments: Map): Double { + return problem.data[sigmaSymbol]?.get(index)?.pow(-2) ?: 1.0 + } + + override fun derivativeOrNull(symbols: List): Expression = Expression { 0.0 } + } + + override fun toString(): String = "PointWeightBySigma($sigmaSymbol)" + + } + + public val byYSigma: PointWeight = bySigma(Symbol.yError) + } +} + +/** + * A fit problem for X-Y-Yerr data. Also known as "least-squares" problem. + */ +public class XYFit( + public val data: XYColumnarData, + public val model: DifferentiableExpression, + override val features: FeatureSet, + internal val pointToCurveDistance: PointToCurveDistance = PointToCurveDistance.byY, + internal val pointWeight: PointWeight = PointWeight.byYSigma, +) : OptimizationProblem { + public fun distance(index: Int): DifferentiableExpression = pointToCurveDistance.distance(this, index) + + public fun weight(index: Int): DifferentiableExpression = pointWeight.weight(this, index) +} + +public fun XYFit.withFeature(vararg features: OptimizationFeature): XYFit { + return XYFit(data, model, this.features.with(*features), pointToCurveDistance, pointWeight) +} + +/** + * Fit given dta with + */ +public suspend fun XYColumnarData.fitWith( + optimizer: Optimizer, + processor: AutoDiffProcessor, + startingPoint: Map, + vararg features: OptimizationFeature = emptyArray(), + xSymbol: Symbol = Symbol.x, + pointToCurveDistance: PointToCurveDistance = PointToCurveDistance.byY, + pointWeight: PointWeight = PointWeight.byYSigma, + model: A.(I) -> I +): XYFit where A : ExtendedField, A : ExpressionAlgebra { + val modelExpression = processor.differentiate { + val x = bindSymbol(xSymbol) + model(x) + } + + var actualFeatures = FeatureSet.of(*features, OptimizationStartPoint(startingPoint)) + + if (actualFeatures.getFeature() == null) { + actualFeatures = actualFeatures.with(OptimizationLog(Loggable.console)) + } + val problem = XYFit( + this, + modelExpression, + actualFeatures, + pointToCurveDistance, + pointWeight + ) + return optimizer.optimize(problem) +} \ No newline at end of file diff --git a/kmath-optimization/src/commonMain/kotlin/space/kscience/kmath/optimization/logLikelihood.kt b/kmath-optimization/src/commonMain/kotlin/space/kscience/kmath/optimization/logLikelihood.kt new file mode 100644 index 000000000..6701887a2 --- /dev/null +++ b/kmath-optimization/src/commonMain/kotlin/space/kscience/kmath/optimization/logLikelihood.kt @@ -0,0 +1,66 @@ +/* + * Copyright 2018-2021 KMath contributors. + * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file. + */ + +package space.kscience.kmath.optimization + +import space.kscience.kmath.data.XYColumnarData +import space.kscience.kmath.data.indices +import space.kscience.kmath.expressions.DifferentiableExpression +import space.kscience.kmath.expressions.Expression +import space.kscience.kmath.expressions.Symbol +import space.kscience.kmath.expressions.derivative +import space.kscience.kmath.misc.UnstableKMathAPI +import kotlin.math.PI +import kotlin.math.ln +import kotlin.math.pow +import kotlin.math.sqrt + + +private val oneOver2Pi = 1.0 / sqrt(2 * PI) + +@UnstableKMathAPI +internal fun XYFit.logLikelihood(): DifferentiableExpression = object : DifferentiableExpression { + override fun derivativeOrNull(symbols: List): Expression = Expression { arguments -> + data.indices.sumOf { index -> + val d = distance(index)(arguments) + val weight = weight(index)(arguments) + val weightDerivative = weight(index)(arguments) + + // -1 / (sqrt(2 PI) * sigma) + 2 (x-mu)/ 2 sigma^2 * d mu/ d theta - (x-mu)^2 / 2 * d w/ d theta + return@sumOf -oneOver2Pi * sqrt(weight) + //offset derivative + d * model.derivative(symbols)(arguments) * weight - //model derivative + d.pow(2) * weightDerivative / 2 //weight derivative + } + } + + override fun invoke(arguments: Map): Double { + return data.indices.sumOf { index -> + val d = distance(index)(arguments) + val weight = weight(index)(arguments) + //1/sqrt(2 PI sigma^2) - (x-mu)^2/ (2 * sigma^2) + oneOver2Pi * ln(weight) - d.pow(2) * weight + } / 2 + } + +} + +/** + * Optimize given XY (least squares) [problem] using this function [Optimizer]. + * The problem is treated as maximum likelihood problem and is done via maximizing logarithmic likelihood, respecting + * possible weight dependency on the model and parameters. + */ +@UnstableKMathAPI +public suspend fun Optimizer>.maximumLogLikelihood(problem: XYFit): XYFit { + val functionOptimization = FunctionOptimization(problem.features, problem.logLikelihood()) + val result = optimize(functionOptimization.withFeatures(FunctionOptimizationTarget.MAXIMIZE)) + return XYFit(problem.data, problem.model, result.features) +} + +@UnstableKMathAPI +public suspend fun Optimizer>.maximumLogLikelihood( + data: XYColumnarData, + model: DifferentiableExpression, + builder: XYOptimizationBuilder.() -> Unit, +): XYFit = maximumLogLikelihood(XYOptimization(data, model, builder)) diff --git a/kmath-stat/build.gradle.kts b/kmath-stat/build.gradle.kts index e3e396b6f..41a1666f8 100644 --- a/kmath-stat/build.gradle.kts +++ b/kmath-stat/build.gradle.kts @@ -1,6 +1,5 @@ plugins { - kotlin("multiplatform") - id("ru.mipt.npm.gradle.common") + id("ru.mipt.npm.gradle.mpp") id("ru.mipt.npm.gradle.native") } diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/optimization/XYOptimization.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/optimization/XYOptimization.kt deleted file mode 100644 index 68c0c77eb..000000000 --- a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/optimization/XYOptimization.kt +++ /dev/null @@ -1,189 +0,0 @@ -/* - * Copyright 2018-2021 KMath contributors. - * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file. - */ -@file:OptIn(UnstableKMathAPI::class) - -package space.kscience.kmath.optimization - -import space.kscience.kmath.data.XYColumnarData -import space.kscience.kmath.data.indices -import space.kscience.kmath.expressions.DifferentiableExpression -import space.kscience.kmath.expressions.Expression -import space.kscience.kmath.expressions.Symbol -import space.kscience.kmath.expressions.derivative -import space.kscience.kmath.misc.FeatureSet -import space.kscience.kmath.misc.UnstableKMathAPI -import kotlin.math.PI -import kotlin.math.ln -import kotlin.math.pow -import kotlin.math.sqrt - -/** - * Specify the way to compute distance from point to the curve as DifferentiableExpression - */ -public interface PointToCurveDistance : OptimizationFeature { - public fun distance(problem: XYOptimization, index: Int): DifferentiableExpression - - public companion object { - public val byY: PointToCurveDistance = object : PointToCurveDistance { - override fun distance(problem: XYOptimization, index: Int): DifferentiableExpression { - - val x = problem.data.x[index] - val y = problem.data.y[index] - return object : DifferentiableExpression { - override fun derivativeOrNull(symbols: List): Expression? = - problem.model.derivativeOrNull(symbols) - - override fun invoke(arguments: Map): Double = - problem.model(arguments + (Symbol.x to x)) - y - } - } - - override fun toString(): String = "PointToCurveDistanceByY" - } - } -} - -/** - * Compute a wight of the point. The more the weight, the more impact this point will have on the fit. - * By default uses Dispersion^-1 - */ -public interface PointWeight : OptimizationFeature { - public fun weight(problem: XYOptimization, index: Int): DifferentiableExpression - - public companion object { - public fun bySigma(sigmaSymbol: Symbol): PointWeight = object : PointWeight { - override fun weight(problem: XYOptimization, index: Int): DifferentiableExpression = - object : DifferentiableExpression { - override fun invoke(arguments: Map): Double { - return problem.data[sigmaSymbol]?.get(index)?.pow(-2) ?: 1.0 - } - - override fun derivativeOrNull(symbols: List): Expression = Expression { 0.0 } - } - - override fun toString(): String = "PointWeightBySigma($sigmaSymbol)" - - } - - public val byYSigma: PointWeight = bySigma(Symbol.yError) - } -} - -/** - * An optimization for XY data. - */ -public class XYOptimization( - override val features: FeatureSet, - public val data: XYColumnarData, - public val model: DifferentiableExpression, - internal val pointToCurveDistance: PointToCurveDistance = PointToCurveDistance.byY, - internal val pointWeight: PointWeight = PointWeight.byYSigma, -) : OptimizationProblem { - public fun distance(index: Int): DifferentiableExpression = pointToCurveDistance.distance(this, index) - - public fun weight(index: Int): DifferentiableExpression = pointWeight.weight(this, index) -} - -public fun XYOptimization.withFeature(vararg features: OptimizationFeature): XYOptimization { - return XYOptimization(this.features.with(*features), data, model, pointToCurveDistance, pointWeight) -} - -private val oneOver2Pi = 1.0 / sqrt(2 * PI) - -internal fun XYOptimization.likelihood(): DifferentiableExpression = object : DifferentiableExpression { - override fun derivativeOrNull(symbols: List): Expression = Expression { arguments -> - data.indices.sumOf { index -> - - val d = distance(index)(arguments) - val weight = weight(index)(arguments) - val weightDerivative = weight(index)(arguments) - - // -1 / (sqrt(2 PI) * sigma) + 2 (x-mu)/ 2 sigma^2 * d mu/ d theta - (x-mu)^2 / 2 * d w/ d theta - return@sumOf -oneOver2Pi * sqrt(weight) + //offset derivative - d * model.derivative(symbols)(arguments) * weight - //model derivative - d.pow(2) * weightDerivative / 2 //weight derivative - } - } - - override fun invoke(arguments: Map): Double { - return data.indices.sumOf { index -> - val d = distance(index)(arguments) - val weight = weight(index)(arguments) - //1/sqrt(2 PI sigma^2) - (x-mu)^2/ (2 * sigma^2) - oneOver2Pi * ln(weight) - d.pow(2) * weight - } / 2 - } - -} - -/** - * Optimize given XY (least squares) [problem] using this function [Optimizer]. - * The problem is treated as maximum likelihood problem and is done via maximizing logarithmic likelihood, respecting - * possible weight dependency on the model and parameters. - */ -public suspend fun Optimizer>.maximumLogLikelihood(problem: XYOptimization): XYOptimization { - val functionOptimization = FunctionOptimization(problem.features, problem.likelihood()) - val result = optimize(functionOptimization.withFeatures(FunctionOptimizationTarget.MAXIMIZE)) - return XYOptimization(result.features, problem.data, problem.model) -} - -public suspend fun Optimizer>.maximumLogLikelihood( - data: XYColumnarData, - model: DifferentiableExpression, - builder: XYOptimizationBuilder.() -> Unit, -): XYOptimization = maximumLogLikelihood(XYOptimization(data, model, builder)) - -//public suspend fun XYColumnarData.fitWith( -// optimizer: XYOptimization, -// problemBuilder: XYOptimizationBuilder.() -> Unit = {}, -// -//) - - -// -//@UnstableKMathAPI -//public interface XYFit : OptimizationProblem { -// -// 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)) } -// } -//} - -// -///** -// * 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) -// function(chiSquared) -// maximize = false -//} \ No newline at end of file diff --git a/settings.gradle.kts b/settings.gradle.kts index f05092bb1..d1cbbe74c 100644 --- a/settings.gradle.kts +++ b/settings.gradle.kts @@ -26,6 +26,7 @@ include( ":kmath-histograms", ":kmath-commons", ":kmath-viktor", + ":kmath-optimization", ":kmath-stat", ":kmath-nd4j", ":kmath-dimensions",