From 7f4f4c7703275e80cb51e38af3a28a7672dbbf5f Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Tue, 24 Jan 2023 21:01:26 +0300 Subject: [PATCH] Update of symbolic operations --- CHANGELOG.md | 2 + build.gradle.kts | 2 +- .../space/kscience/kmath/fit/chiSquared.kt | 2 +- .../kotlin/space/kscience/kmath/fit/qowFit.kt | 21 ++-- .../expressions/ExpressionWithDefault.kt | 31 +++++ .../kscience/kmath/expressions/NamedMatrix.kt | 37 ++++++ .../kmath/optimization/OptimizationProblem.kt | 17 ++- .../kmath/optimization/QowOptimizer.kt | 116 ++++++++++-------- .../kscience/kmath/optimization/XYFit.kt | 7 +- .../kmath/stat/chiSquaredExpression.kt | 7 +- 10 files changed, 173 insertions(+), 69 deletions(-) create mode 100644 kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/ExpressionWithDefault.kt create mode 100644 kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/NamedMatrix.kt rename kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/specialExpressions.kt => kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/chiSquaredExpression.kt (87%) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1b3ab75fd..b8cb0d976 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,8 @@ ## [Unreleased] ### Added +- `NamedMatrix` - matrix with symbol-based indexing +- `Expression` with default arguments - Type-aliases for numbers like `Float64` - 2D optimal trajectory computation in a separate module `kmath-trajectory` - Autodiff for generic algebra elements in core! diff --git a/build.gradle.kts b/build.gradle.kts index 4efb6c802..8731dbc70 100644 --- a/build.gradle.kts +++ b/build.gradle.kts @@ -15,7 +15,7 @@ allprojects { } group = "space.kscience" - version = "0.3.1-dev-8" + version = "0.3.1-dev-9" } subprojects { diff --git a/examples/src/main/kotlin/space/kscience/kmath/fit/chiSquared.kt b/examples/src/main/kotlin/space/kscience/kmath/fit/chiSquared.kt index c5f694ce0..258ed0c84 100644 --- a/examples/src/main/kotlin/space/kscience/kmath/fit/chiSquared.kt +++ b/examples/src/main/kotlin/space/kscience/kmath/fit/chiSquared.kt @@ -10,7 +10,6 @@ import kotlinx.html.h3 import space.kscience.kmath.commons.optimization.CMOptimizer import space.kscience.kmath.distributions.NormalDistribution import space.kscience.kmath.expressions.autodiff -import space.kscience.kmath.expressions.chiSquaredExpression import space.kscience.kmath.expressions.symbol import space.kscience.kmath.operations.asIterable import space.kscience.kmath.operations.toList @@ -22,6 +21,7 @@ import space.kscience.kmath.random.RandomGenerator import space.kscience.kmath.real.DoubleVector import space.kscience.kmath.real.map import space.kscience.kmath.real.step +import space.kscience.kmath.stat.chiSquaredExpression import space.kscience.plotly.* import space.kscience.plotly.models.ScatterMode import space.kscience.plotly.models.TraceValues diff --git a/examples/src/main/kotlin/space/kscience/kmath/fit/qowFit.kt b/examples/src/main/kotlin/space/kscience/kmath/fit/qowFit.kt index 9650fa630..fe7f48b72 100644 --- a/examples/src/main/kotlin/space/kscience/kmath/fit/qowFit.kt +++ b/examples/src/main/kotlin/space/kscience/kmath/fit/qowFit.kt @@ -15,10 +15,7 @@ import space.kscience.kmath.expressions.binding import space.kscience.kmath.expressions.symbol import space.kscience.kmath.operations.asIterable import space.kscience.kmath.operations.toList -import space.kscience.kmath.optimization.QowOptimizer -import space.kscience.kmath.optimization.chiSquaredOrNull -import space.kscience.kmath.optimization.fitWith -import space.kscience.kmath.optimization.resultPoint +import space.kscience.kmath.optimization.* import space.kscience.kmath.random.RandomGenerator import space.kscience.kmath.real.map import space.kscience.kmath.real.step @@ -32,6 +29,8 @@ import kotlin.math.sqrt private val a by symbol private val b by symbol private val c by symbol +private val d by symbol +private val e by symbol /** @@ -64,16 +63,22 @@ suspend fun main() { val result = XYErrorColumnarData.of(x, y, yErr).fitWith( QowOptimizer, Double.autodiff, - mapOf(a to 0.9, b to 1.2, c to 2.0) + mapOf(a to 0.9, b to 1.2, c to 2.0, e to 1.0, d to 1.0, e to 0.0), + OptimizationParameters(a, b, c, d) ) { 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 + val d by binding + val e by binding + + a * arg.pow(2) + b * arg + c + d * arg.pow(3) + e / arg } + println("Resulting chi2/dof: ${result.chiSquaredOrNull}/${result.dof}") + //display a page with plot and numerical results val page = Plotly.page { plot { @@ -89,7 +94,7 @@ suspend fun main() { scatter { mode = ScatterMode.lines x(x) - y(x.map { result.model(result.resultPoint + (Symbol.x to it)) }) + y(x.map { result.model(result.startPoint + result.resultPoint + (Symbol.x to it)) }) name = "fit" } } @@ -98,7 +103,7 @@ suspend fun main() { +"Fit result: ${result.resultPoint}" } h3 { - +"Chi2/dof = ${result.chiSquaredOrNull!! / (x.size - 3)}" + +"Chi2/dof = ${result.chiSquaredOrNull!! / result.dof}" } } diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/ExpressionWithDefault.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/ExpressionWithDefault.kt new file mode 100644 index 000000000..c802fe04c --- /dev/null +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/ExpressionWithDefault.kt @@ -0,0 +1,31 @@ +/* + * Copyright 2018-2023 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.expressions + +public class ExpressionWithDefault( + private val origin: Expression, + private val defaultArgs: Map, +) : Expression { + override fun invoke(arguments: Map): T = origin.invoke(defaultArgs + arguments) +} + +public fun Expression.withDefaultArgs(defaultArgs: Map): ExpressionWithDefault = + ExpressionWithDefault(this, defaultArgs) + + +public class DiffExpressionWithDefault( + private val origin: DifferentiableExpression, + private val defaultArgs: Map, +) : DifferentiableExpression { + + override fun invoke(arguments: Map): T = origin.invoke(defaultArgs + arguments) + + override fun derivativeOrNull(symbols: List): Expression? = + origin.derivativeOrNull(symbols)?.withDefaultArgs(defaultArgs) +} + +public fun DifferentiableExpression.withDefaultArgs(defaultArgs: Map): DiffExpressionWithDefault = + DiffExpressionWithDefault(this, defaultArgs) \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/NamedMatrix.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/NamedMatrix.kt new file mode 100644 index 000000000..f0b27ffa2 --- /dev/null +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/NamedMatrix.kt @@ -0,0 +1,37 @@ +/* + * Copyright 2018-2023 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.expressions + +import space.kscience.kmath.linear.Matrix +import space.kscience.kmath.misc.UnstableKMathAPI +import space.kscience.kmath.structures.getOrNull + +public class NamedMatrix(public val values: Matrix, public val indexer: SymbolIndexer) : Matrix by values { + public operator fun get(i: Symbol, j: Symbol): T = get(indexer.indexOf(i), indexer.indexOf(j)) + + public companion object { + + public fun toStringWithSymbols(values: Matrix<*>, indexer: SymbolIndexer): String = buildString { + appendLine(indexer.symbols.joinToString(separator = "\t", prefix = "\t\t")) + indexer.symbols.forEach { i -> + append(i.identity + "\t") + values.rows.getOrNull(indexer.indexOf(i))?.let { row -> + indexer.symbols.forEach { j -> + append(row.getOrNull(indexer.indexOf(j)).toString()) + append("\t") + } + appendLine() + } + } + } + } +} + +public fun Matrix.named(indexer: SymbolIndexer): NamedMatrix = NamedMatrix(this, indexer) + +public fun Matrix.named(symbols: List): NamedMatrix = named(SimpleSymbolIndexer(symbols)) \ No newline at end of file diff --git a/kmath-optimization/src/commonMain/kotlin/space/kscience/kmath/optimization/OptimizationProblem.kt b/kmath-optimization/src/commonMain/kotlin/space/kscience/kmath/optimization/OptimizationProblem.kt index db154be42..9fdcfc53d 100644 --- a/kmath-optimization/src/commonMain/kotlin/space/kscience/kmath/optimization/OptimizationProblem.kt +++ b/kmath-optimization/src/commonMain/kotlin/space/kscience/kmath/optimization/OptimizationProblem.kt @@ -6,8 +6,8 @@ package space.kscience.kmath.optimization import space.kscience.kmath.expressions.DifferentiableExpression +import space.kscience.kmath.expressions.NamedMatrix import space.kscience.kmath.expressions.Symbol -import space.kscience.kmath.linear.Matrix import space.kscience.kmath.misc.* import kotlin.reflect.KClass @@ -32,7 +32,10 @@ public interface OptimizationPrior : OptimizationFeature, DifferentiableExpre override val key: FeatureKey get() = OptimizationPrior::class } -public class OptimizationCovariance(public val covariance: Matrix) : OptimizationFeature { +/** + * Covariance matrix for + */ +public class OptimizationCovariance(public val covariance: NamedMatrix) : OptimizationFeature { override fun toString(): String = "Covariance($covariance)" } @@ -57,10 +60,20 @@ public class OptimizationLog(private val loggable: Loggable) : Loggable by logga override fun toString(): String = "Log($loggable)" } +/** + * Free parameters of the optimization + */ public class OptimizationParameters(public val symbols: List) : OptimizationFeature { public constructor(vararg symbols: Symbol) : this(listOf(*symbols)) override fun toString(): String = "Parameters($symbols)" } +/** + * Maximum allowed number of iterations + */ +public class OptimizationIterations(public val maxIterations: Int) : OptimizationFeature { + override fun toString(): String = "Iterations($maxIterations)" +} + diff --git a/kmath-optimization/src/commonMain/kotlin/space/kscience/kmath/optimization/QowOptimizer.kt b/kmath-optimization/src/commonMain/kotlin/space/kscience/kmath/optimization/QowOptimizer.kt index 669860934..bbebec6af 100644 --- a/kmath-optimization/src/commonMain/kotlin/space/kscience/kmath/optimization/QowOptimizer.kt +++ b/kmath-optimization/src/commonMain/kotlin/space/kscience/kmath/optimization/QowOptimizer.kt @@ -5,10 +5,7 @@ package space.kscience.kmath.optimization -import space.kscience.kmath.expressions.DifferentiableExpression -import space.kscience.kmath.expressions.Symbol -import space.kscience.kmath.expressions.SymbolIndexer -import space.kscience.kmath.expressions.derivative +import space.kscience.kmath.expressions.* import space.kscience.kmath.linear.* import space.kscience.kmath.misc.UnstableKMathAPI import space.kscience.kmath.misc.log @@ -16,6 +13,7 @@ import space.kscience.kmath.operations.DoubleField import space.kscience.kmath.operations.DoubleL2Norm import space.kscience.kmath.operations.algebra import space.kscience.kmath.structures.DoubleBuffer +import kotlin.math.abs public class QowRuns(public val runs: Int) : OptimizationFeature { @@ -40,18 +38,24 @@ public object QowOptimizer : Optimizer { @OptIn(UnstableKMathAPI::class) private class QoWeight( val problem: XYFit, - val parameters: Map, - ) : Map by parameters, SymbolIndexer { - override val symbols: List = parameters.keys.toList() + val freeParameters: Map, + ) : SymbolIndexer { + val size get() = freeParameters.size + + override val symbols: List = freeParameters.keys.toList() val data get() = problem.data + val allParameters by lazy { + problem.startPoint + freeParameters + } + /** * 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) { d, s -> - problem.distance(d).derivative(symbols[s])(parameters) + problem.distance(d).derivative(symbols[s]).invoke(allParameters) } } @@ -60,29 +64,31 @@ public object QowOptimizer : Optimizer { */ val dispersion: Point by lazy { DoubleBuffer(problem.data.size) { d -> - 1.0/problem.weight(d).invoke(parameters) + 1.0 / problem.weight(d).invoke(allParameters) } } - val prior: DifferentiableExpression? get() = problem.getFeature>() + val prior: DifferentiableExpression? + get() = problem.getFeature>()?.withDefaultArgs(allParameters) - override fun toString(): String = parameters.toString() + override fun toString(): String = freeParameters.toString() } /** * The signed distance from the model to the [d]-th point of data. */ - private fun QoWeight.distance(d: Int, parameters: Map): Double = problem.distance(d)(parameters) + private fun QoWeight.distance(d: Int, parameters: Map): Double = + problem.distance(d)(allParameters + parameters) /** * The derivative of [distance] */ private fun QoWeight.distanceDerivative(symbol: Symbol, d: Int, parameters: Map): Double = - problem.distance(d).derivative(symbol)(parameters) + problem.distance(d).derivative(symbol).invoke(allParameters + parameters) /** - * Теоретическая ковариация весовых функций. + * Theoretical covariance of weight functions * * D(\phi)=E(\phi_k(\theta_0) \phi_l(\theta_0))= disDeriv_k * disDeriv_l /sigma^2 */ @@ -92,7 +98,7 @@ public object QowOptimizer : Optimizer { } /** - * Экспериментальная ковариация весов. Формула (22) из + * Experimental covariance Eq (22) from * http://arxiv.org/abs/physics/0604127 */ private fun QoWeight.covarFExp(theta: Map): Matrix = @@ -115,10 +121,9 @@ public object QowOptimizer : Optimizer { * Equation derivatives for Newton run */ private fun QoWeight.getEqDerivValues( - theta: Map = parameters, + theta: Map = freeParameters, ): Matrix = with(linearSpace) { - //Возвращает производную k-того Eq по l-тому параметру - //val res = Array(fitDim) { DoubleArray(fitDim) } + //Derivative of k Eq over l parameter val sderiv = buildMatrix(data.size, size) { d, s -> distanceDerivative(symbols[s], d, theta) } @@ -140,16 +145,15 @@ public object QowOptimizer : Optimizer { /** - * Значения уравнений метода квазиоптимальных весов + * Quasi optimal weights equations values */ - private fun QoWeight.getEqValues(theta: Map = this): Point { + private fun QoWeight.getEqValues(theta: Map): Point { val distances = DoubleBuffer(data.size) { d -> distance(d, theta) } - return DoubleBuffer(size) { s -> val base = (0 until data.size).sumOf { d -> distances[d] * derivs[d, s] / dispersion[d] } - //Поправка на априорную вероятность + //Prior probability correction prior?.let { prior -> - base - prior.derivative(symbols[s])(theta) / prior(theta) + base - prior.derivative(symbols[s]).invoke(theta) / prior(theta) } ?: base } } @@ -157,15 +161,13 @@ public object QowOptimizer : Optimizer { private fun QoWeight.newtonianStep( theta: Map, - eqvalues: Point, + eqValues: Point, ): QoWeight = linearSpace { - with(this@newtonianStep) { - val start = theta.toPoint() - val invJacob = solver.inverse(this@newtonianStep.getEqDerivValues(theta)) + val start = theta.toPoint() + val invJacob = solver.inverse(getEqDerivValues(theta)) - val step = invJacob.dot(eqvalues) - return QoWeight(problem, theta + (start - step).toMap()) - } + val step = invJacob.dot(eqValues) + return QoWeight(problem, theta + (start - step).toMap()) } private fun QoWeight.newtonianRun( @@ -177,10 +179,10 @@ public object QowOptimizer : Optimizer { val logger = problem.getFeature() 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 par = freeParameters + + logger?.log { "Starting newtonian iteration from: \n\t$allParameters" } var eqvalues = getEqValues(par) //Values of the weight functions @@ -193,48 +195,48 @@ public object QowOptimizer : Optimizer { logger?.log { "Starting step number $i" } val currentSolution = if (fast) { - //Берет значения матрицы в той точке, где считается вес - newtonianStep(this, eqvalues) + //Matrix values in the point of weight computation + newtonianStep(freeParameters, eqvalues) } else { - //Берет значения матрицы в точке par + //Matrix values in a current point newtonianStep(par, eqvalues) } // здесь должен стоять учет границ параметров logger?.log { "Parameter values after step are: \n\t$currentSolution" } - eqvalues = getEqValues(currentSolution) - val currentDis = DoubleL2Norm.norm(eqvalues)// невязка после шага + eqvalues = getEqValues(currentSolution.freeParameters) + val currentDis = DoubleL2Norm.norm(eqvalues)// discrepancy after the step logger?.log { "The discrepancy after step is: $currentDis." } if (currentDis >= dis && i > 1) { - //дополнительно проверяем, чтобы был сделан хотя бы один шаг + //Check if one step is made flag = true logger?.log { "The discrepancy does not decrease. Stopping iteration." } + } else if (abs(dis - currentDis) <= tolerance) { + flag = true + par = currentSolution.freeParameters + logger?.log { "Relative discrepancy tolerance threshold is reached. Stopping iteration." } } else { - par = currentSolution + par = currentSolution.freeParameters dis = currentDis } if (i >= maxSteps) { flag = true logger?.log { "Maximum number of iterations reached. Stopping iteration." } } - if (dis <= tolerance) { - flag = true - logger?.log { "Tolerance threshold is reached. Stopping iteration." } - } } return QoWeight(problem, par) } - private fun QoWeight.covariance(): Matrix { + private fun QoWeight.covariance(): NamedMatrix { val logger = problem.getFeature() logger?.log { """ - Starting errors estimation using quasioptimal weights method. The starting weight is: - ${problem.startPoint} + Starting errors estimation using quasi-optimal weights method. The starting weight is: + $allParameters """.trimIndent() } @@ -248,19 +250,27 @@ public object QowOptimizer : Optimizer { // valid = false // } // } - return covar + logger?.log { + "Covariance matrix:" + "\n" + NamedMatrix.toStringWithSymbols(covar, this) + } + return covar.named(symbols) } override suspend fun optimize(problem: XYFit): XYFit { val qowRuns = problem.getFeature()?.runs ?: 2 + val iterations = problem.getFeature()?.maxIterations ?: 50 + val freeParameters: Map = problem.getFeature()?.let { op -> + problem.startPoint.filterKeys { it in op.symbols } + } ?: problem.startPoint - var qow = QoWeight(problem, problem.startPoint) - var res = qow.newtonianRun() + var qow = QoWeight(problem, freeParameters) + var res = qow.newtonianRun(maxSteps = iterations) repeat(qowRuns - 1) { - qow = QoWeight(problem, res.parameters) - res = qow.newtonianRun() + qow = QoWeight(problem, res.freeParameters) + res = qow.newtonianRun(maxSteps = iterations) } - return res.problem.withFeature(OptimizationResult(res.parameters)) + val covariance = res.covariance() + return res.problem.withFeature(OptimizationResult(res.freeParameters), OptimizationCovariance(covariance)) } } \ No newline at end of file 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 index bba7a7113..5c28826ee 100644 --- a/kmath-optimization/src/commonMain/kotlin/space/kscience/kmath/optimization/XYFit.kt +++ b/kmath-optimization/src/commonMain/kotlin/space/kscience/kmath/optimization/XYFit.kt @@ -152,7 +152,7 @@ public suspend fun XYColumnarData.fitWith( */ public val XYFit.chiSquaredOrNull: Double? get() { - val result = resultPointOrNull ?: return null + val result = startPoint + (resultPointOrNull ?: return null) return data.indices.sumOf { index -> @@ -164,4 +164,7 @@ public val XYFit.chiSquaredOrNull: Double? ((y - mu) / yErr).pow(2) } - } \ No newline at end of file + } + +public val XYFit.dof: Int + get() = data.size - (getFeature()?.symbols?.size ?: startPoint.size) \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/specialExpressions.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/chiSquaredExpression.kt similarity index 87% rename from kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/specialExpressions.kt rename to kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/chiSquaredExpression.kt index 59dfeb8ea..ca9755ad5 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/specialExpressions.kt +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/chiSquaredExpression.kt @@ -1,10 +1,13 @@ /* - * Copyright 2018-2022 KMath contributors. + * Copyright 2018-2023 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.expressions +package space.kscience.kmath.stat +import space.kscience.kmath.expressions.AutoDiffProcessor +import space.kscience.kmath.expressions.DifferentiableExpression +import space.kscience.kmath.expressions.ExpressionAlgebra import space.kscience.kmath.operations.ExtendedField import space.kscience.kmath.operations.asIterable import space.kscience.kmath.structures.Buffer