From f8c3d1793c80f7f6fcdbd53b39e9c004cd535a3a Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Wed, 28 Oct 2020 09:08:37 +0300 Subject: [PATCH] Fitting refactor --- .../kmath/commons/optimization/CMFit.kt | 95 ------------------- .../kmath/commons/optimization/cmFit.kt | 68 +++++++++++++ .../commons/optimization/OptimizeTest.kt | 20 ++-- .../kmath/prob/{Fit.kt => Fitting.kt} | 31 +++++- 4 files changed, 105 insertions(+), 109 deletions(-) delete mode 100644 kmath-commons/src/main/kotlin/kscience/kmath/commons/optimization/CMFit.kt create mode 100644 kmath-commons/src/main/kotlin/kscience/kmath/commons/optimization/cmFit.kt rename kmath-prob/src/commonMain/kotlin/kscience/kmath/prob/{Fit.kt => Fitting.kt} (52%) diff --git a/kmath-commons/src/main/kotlin/kscience/kmath/commons/optimization/CMFit.kt b/kmath-commons/src/main/kotlin/kscience/kmath/commons/optimization/CMFit.kt deleted file mode 100644 index 3143dcca5..000000000 --- a/kmath-commons/src/main/kotlin/kscience/kmath/commons/optimization/CMFit.kt +++ /dev/null @@ -1,95 +0,0 @@ -package kscience.kmath.commons.optimization - -import kscience.kmath.commons.expressions.DerivativeStructureExpression -import kscience.kmath.commons.expressions.DerivativeStructureField -import kscience.kmath.expressions.DifferentiableExpression -import kscience.kmath.expressions.Expression -import kscience.kmath.expressions.StringSymbol -import kscience.kmath.expressions.Symbol -import kscience.kmath.structures.Buffer -import kscience.kmath.structures.indices -import org.apache.commons.math3.analysis.differentiation.DerivativeStructure -import org.apache.commons.math3.optim.nonlinear.scalar.GoalType -import kotlin.math.pow - - -public object CMFit { - - /** - * Generate a chi squared expression from given x-y-sigma model represented by an expression. Does not provide derivatives - * TODO move to prob/stat - */ - 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) - } - } - } - - /** - * Generate a chi squared expression from given x-y-sigma data and inline model. Provides automatic differentiation - */ - public fun chiSquared( - x: Buffer, - y: Buffer, - yErr: Buffer, - model: DerivativeStructureField.(x: DerivativeStructure) -> DerivativeStructure, - ): DerivativeStructureExpression { - 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 DerivativeStructureExpression { - var sum = zero - x.indices.forEach { - val xValue = x[it] - val yValue = y[it] - val yErrValue = yErr[it] - val modelValue = model(const(xValue)) - sum += ((yValue - modelValue) / yErrValue).pow(2) - } - sum - } - } -} - -/** - * Optimize expression without derivatives - */ -public fun Expression.optimize( - vararg symbols: Symbol, - configuration: CMOptimizationProblem.() -> Unit, -): OptimizationResult = optimizeWith(CMOptimizationProblem, symbols = symbols, configuration) - - -/** - * Optimize differentiable expression - */ -public fun DifferentiableExpression.optimize( - vararg symbols: Symbol, - configuration: CMOptimizationProblem.() -> Unit, -): OptimizationResult = optimizeWith(CMOptimizationProblem, symbols = symbols, configuration) - -public fun DifferentiableExpression.minimize( - vararg startPoint: Pair, - configuration: CMOptimizationProblem.() -> Unit = {}, -): OptimizationResult { - require(startPoint.isNotEmpty()) { "Must provide a list of symbols for optimization" } - val problem = CMOptimizationProblem(startPoint.map { it.first }).apply(configuration) - problem.diffExpression(this) - problem.initialGuess(startPoint.toMap()) - problem.goal(GoalType.MINIMIZE) - return problem.optimize() -} \ No newline at end of file diff --git a/kmath-commons/src/main/kotlin/kscience/kmath/commons/optimization/cmFit.kt b/kmath-commons/src/main/kotlin/kscience/kmath/commons/optimization/cmFit.kt new file mode 100644 index 000000000..24df3177d --- /dev/null +++ b/kmath-commons/src/main/kotlin/kscience/kmath/commons/optimization/cmFit.kt @@ -0,0 +1,68 @@ +package kscience.kmath.commons.optimization + +import kscience.kmath.commons.expressions.DerivativeStructureField +import kscience.kmath.expressions.DifferentiableExpression +import kscience.kmath.expressions.Expression +import kscience.kmath.expressions.Symbol +import kscience.kmath.prob.Fitting +import kscience.kmath.structures.Buffer +import kscience.kmath.structures.asBuffer +import org.apache.commons.math3.analysis.differentiation.DerivativeStructure +import org.apache.commons.math3.optim.nonlinear.scalar.GoalType + + +/** + * Generate a chi squared expression from given x-y-sigma data and inline model. Provides automatic differentiation + */ +public fun Fitting.chiSquared( + x: Buffer, + y: Buffer, + yErr: Buffer, + model: DerivativeStructureField.(x: DerivativeStructure) -> DerivativeStructure, +): DifferentiableExpression = chiSquared(DerivativeStructureField, x, y, yErr, model) + +/** + * Generate a chi squared expression from given x-y-sigma data and inline model. Provides automatic differentiation + */ +public fun Fitting.chiSquared( + x: Iterable, + y: Iterable, + yErr: Iterable, + model: DerivativeStructureField.(x: DerivativeStructure) -> DerivativeStructure, +): DifferentiableExpression = chiSquared( + DerivativeStructureField, + x.toList().asBuffer(), + y.toList().asBuffer(), + yErr.toList().asBuffer(), + model +) + + +/** + * Optimize expression without derivatives + */ +public fun Expression.optimize( + vararg symbols: Symbol, + configuration: CMOptimizationProblem.() -> Unit, +): OptimizationResult = optimizeWith(CMOptimizationProblem, symbols = symbols, configuration) + + +/** + * Optimize differentiable expression + */ +public fun DifferentiableExpression.optimize( + vararg symbols: Symbol, + configuration: CMOptimizationProblem.() -> Unit, +): OptimizationResult = optimizeWith(CMOptimizationProblem, symbols = symbols, configuration) + +public fun DifferentiableExpression.minimize( + vararg startPoint: Pair, + configuration: CMOptimizationProblem.() -> Unit = {}, +): OptimizationResult { + require(startPoint.isNotEmpty()) { "Must provide a list of symbols for optimization" } + val problem = CMOptimizationProblem(startPoint.map { it.first }).apply(configuration) + problem.diffExpression(this) + problem.initialGuess(startPoint.toMap()) + problem.goal(GoalType.MINIMIZE) + return problem.optimize() +} \ No newline at end of file diff --git a/kmath-commons/src/test/kotlin/kscience/kmath/commons/optimization/OptimizeTest.kt b/kmath-commons/src/test/kotlin/kscience/kmath/commons/optimization/OptimizeTest.kt index d9fc5ebef..502ed40f8 100644 --- a/kmath-commons/src/test/kotlin/kscience/kmath/commons/optimization/OptimizeTest.kt +++ b/kmath-commons/src/test/kotlin/kscience/kmath/commons/optimization/OptimizeTest.kt @@ -3,6 +3,7 @@ package kscience.kmath.commons.optimization import kscience.kmath.commons.expressions.DerivativeStructureExpression import kscience.kmath.expressions.symbol import kscience.kmath.prob.Distribution +import kscience.kmath.prob.Fitting import kscience.kmath.prob.RandomGenerator import kscience.kmath.prob.normal import kscience.kmath.structures.asBuffer @@ -39,7 +40,7 @@ internal class OptimizeTest { } @Test - fun testFit() { + fun testCmFit() { val a by symbol val b by symbol val c by symbol @@ -52,15 +53,14 @@ internal class OptimizeTest { it.pow(2) + it + 1 + chain.nextDouble() } val yErr = x.map { sigma } - with(CMFit) { - val chi2 = chiSquared(x.asBuffer(), y.asBuffer(), yErr.asBuffer()) { x -> - val cWithDefault = bindOrNull(c)?: one - bind(a) * x.pow(2) + bind(b) * x + cWithDefault - } - - val result = chi2.minimize(a to 1.5, b to 0.9, c to 1.0) - println(result) - println("Chi2/dof = ${result.value / (x.size - 3)}") + val chi2 = Fitting.chiSquared(x.asBuffer(), y.asBuffer(), yErr.asBuffer()) { x -> + val cWithDefault = bindOrNull(c) ?: one + bind(a) * x.pow(2) + bind(b) * x + cWithDefault } + + val result = chi2.minimize(a to 1.5, b to 0.9, c to 1.0) + println(result) + println("Chi2/dof = ${result.value / (x.size - 3)}") } + } \ No newline at end of file diff --git a/kmath-prob/src/commonMain/kotlin/kscience/kmath/prob/Fit.kt b/kmath-prob/src/commonMain/kotlin/kscience/kmath/prob/Fitting.kt similarity index 52% rename from kmath-prob/src/commonMain/kotlin/kscience/kmath/prob/Fit.kt rename to kmath-prob/src/commonMain/kotlin/kscience/kmath/prob/Fitting.kt index efe582212..97548d676 100644 --- a/kmath-prob/src/commonMain/kotlin/kscience/kmath/prob/Fit.kt +++ b/kmath-prob/src/commonMain/kotlin/kscience/kmath/prob/Fitting.kt @@ -1,13 +1,12 @@ package kscience.kmath.prob -import kscience.kmath.expressions.AutoDiffProcessor -import kscience.kmath.expressions.DifferentiableExpression -import kscience.kmath.expressions.ExpressionAlgebra +import kscience.kmath.expressions.* import kscience.kmath.operations.ExtendedField import kscience.kmath.structures.Buffer import kscience.kmath.structures.indices +import kotlin.math.pow -public object Fit { +public object Fitting { /** * Generate a chi squared expression from given x-y-sigma data and inline model. Provides automatic differentiation @@ -33,4 +32,28 @@ public object Fit { 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) + } + } + } } \ No newline at end of file