diff --git a/kmath-commons/src/main/kotlin/kscience/kmath/commons/expressions/DerivativeStructureExpression.kt b/kmath-commons/src/main/kotlin/kscience/kmath/commons/expressions/DerivativeStructureExpression.kt index 376fea7a3..272501729 100644 --- a/kmath-commons/src/main/kotlin/kscience/kmath/commons/expressions/DerivativeStructureExpression.kt +++ b/kmath-commons/src/main/kotlin/kscience/kmath/commons/expressions/DerivativeStructureExpression.kt @@ -38,13 +38,13 @@ public class DerivativeStructureField( key.identity to DerivativeStructureSymbol(key, value) } - override fun const(value: Double): DerivativeStructure = DerivativeStructure(order, bindings.size, value) + override fun const(value: Double): DerivativeStructure = DerivativeStructure(bindings.size, order, value) public override fun bindOrNull(symbol: Symbol): DerivativeStructureSymbol? = variables[symbol.identity] public fun bind(symbol: Symbol): DerivativeStructureSymbol = variables.getValue(symbol.identity) - public fun Number.const(): DerivativeStructure = const(toDouble()) + //public fun Number.const(): DerivativeStructure = const(toDouble()) public fun DerivativeStructure.derivative(parameter: Symbol, order: Int = 1): Double { return derivative(mapOf(parameter to order)) 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 index 4ffd0559d..a62630ed3 100644 --- a/kmath-commons/src/main/kotlin/kscience/kmath/commons/optimization/CMFit.kt +++ b/kmath-commons/src/main/kotlin/kscience/kmath/commons/optimization/CMFit.kt @@ -17,9 +17,9 @@ 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 core/separate module + * TODO move to prob/stat */ - public fun chiSquaredExpression( + public fun chiSquared( x: Buffer, y: Buffer, yErr: Buffer, @@ -35,7 +35,7 @@ public object CMFit { val yErrValue = yErr[it] val modifiedArgs = arguments + (xSymbol to xValue) val modelValue = model(modifiedArgs) - ((yValue - modelValue) / yErrValue).pow(2) / 2 + ((yValue - modelValue) / yErrValue).pow(2) } } } @@ -43,7 +43,7 @@ public object CMFit { /** * Generate a chi squared expression from given x-y-sigma data and inline model. Provides automatic differentiation */ - public fun chiSquaredExpression( + public fun chiSquared( x: Buffer, y: Buffer, yErr: Buffer, @@ -58,7 +58,7 @@ public object CMFit { val yValue = y[it] val yErrValue = yErr[it] val modelValue = model(const(xValue)) - sum += ((yValue - modelValue) / yErrValue).pow(2) / 2 + sum += ((yValue - modelValue) / yErrValue).pow(2) } sum } @@ -92,12 +92,13 @@ public fun DifferentiableExpression.optimize( } public fun DifferentiableExpression.minimize( - vararg symbols: Symbol, - configuration: CMOptimizationProblem.() -> Unit, + vararg startPoint: Pair, + configuration: CMOptimizationProblem.() -> Unit = {}, ): OptimizationResult { - require(symbols.isNotEmpty()) { "Must provide a list of symbols for optimization" } - val problem = CMOptimizationProblem(symbols.toList()).apply(configuration) + 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/CMOptimizationProblem.kt b/kmath-commons/src/main/kotlin/kscience/kmath/commons/optimization/CMOptimizationProblem.kt index b5ea59d6b..2ca907d05 100644 --- a/kmath-commons/src/main/kotlin/kscience/kmath/commons/optimization/CMOptimizationProblem.kt +++ b/kmath-commons/src/main/kotlin/kscience/kmath/commons/optimization/CMOptimizationProblem.kt @@ -17,14 +17,13 @@ public operator fun PointValuePair.component2(): Double = value public class CMOptimizationProblem( override val symbols: List, -) : OptimizationProblem, SymbolIndexer { - protected val optimizationData: HashMap, OptimizationData> = HashMap() +) : OptimizationProblem, 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 fun addOptimizationData(data: OptimizationData) { + public fun addOptimizationData(data: OptimizationData) { optimizationData[data::class] = data } @@ -32,6 +31,8 @@ public class CMOptimizationProblem( addOptimizationData(MaxEval.unlimited()) } + public fun exportOptimizationData(): List = optimizationData.values.toList() + public fun initialGuess(map: Map): Unit { addOptimizationData(InitialGuess(map.toDoubleArray())) } @@ -90,7 +91,7 @@ public class CMOptimizationProblem( override fun optimize(): OptimizationResult { val optimizer = optimizatorBuilder?.invoke() ?: error("Optimizer not defined") val (point, value) = optimizer.optimize(*optimizationData.values.toTypedArray()) - return OptimizationResult(point.toMap(), value) + return OptimizationResult(point.toMap(), value, setOf(this)) } public companion object { diff --git a/kmath-commons/src/main/kotlin/kscience/kmath/commons/optimization/OptimizationProblem.kt b/kmath-commons/src/main/kotlin/kscience/kmath/commons/optimization/OptimizationProblem.kt index e52450be1..a246a817b 100644 --- a/kmath-commons/src/main/kotlin/kscience/kmath/commons/optimization/OptimizationProblem.kt +++ b/kmath-commons/src/main/kotlin/kscience/kmath/commons/optimization/OptimizationProblem.kt @@ -4,14 +4,19 @@ import kscience.kmath.expressions.DifferentiableExpression import kscience.kmath.expressions.Expression import kscience.kmath.expressions.Symbol -public interface OptimizationResultFeature +public interface OptimizationFeature +//TODO move to prob/stat public class OptimizationResult( public val point: Map, public val value: T, - public val features: Set = emptySet(), -) + public val features: Set = emptySet(), +){ + override fun toString(): String { + return "OptimizationResult(point=$point, value=$value)" + } +} /** * A configuration builder for optimization problem diff --git a/kmath-commons/src/main/kotlin/kscience/kmath/commons/random/CMRandomGeneratorWrapper.kt b/kmath-commons/src/main/kotlin/kscience/kmath/commons/random/CMRandomGeneratorWrapper.kt index 58609deae..9600f6901 100644 --- a/kmath-commons/src/main/kotlin/kscience/kmath/commons/random/CMRandomGeneratorWrapper.kt +++ b/kmath-commons/src/main/kotlin/kscience/kmath/commons/random/CMRandomGeneratorWrapper.kt @@ -2,8 +2,9 @@ package kscience.kmath.commons.random import kscience.kmath.prob.RandomGenerator -public class CMRandomGeneratorWrapper(public val factory: (IntArray) -> RandomGenerator) : - org.apache.commons.math3.random.RandomGenerator { +public class CMRandomGeneratorWrapper( + public val factory: (IntArray) -> RandomGenerator, +) : org.apache.commons.math3.random.RandomGenerator { private var generator: RandomGenerator = factory(intArrayOf()) public override fun nextBoolean(): Boolean = generator.nextBoolean() 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 07bda2aa4..ff5542235 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 @@ -2,14 +2,19 @@ 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.RandomGenerator +import kscience.kmath.prob.normal +import kscience.kmath.structures.asBuffer import org.junit.jupiter.api.Test +import kotlin.math.pow internal class OptimizeTest { val x by symbol val y by symbol val normal = DerivativeStructureExpression { - exp(-bind(x).pow(2) / 2) + exp(- bind(y).pow(2) / 2) + exp(-bind(x).pow(2) / 2) + exp(-bind(y).pow(2) / 2) } @Test @@ -32,4 +37,29 @@ internal class OptimizeTest { println(result.point) println(result.value) } + + @Test + fun testFit() { + val a by symbol + val b by symbol + val c by symbol + + val sigma = 1.0 + val generator = Distribution.normal(0.0, sigma) + val chain = generator.sample(RandomGenerator.default(1126)) + val x = (1..100).map { it.toDouble() } + val y = x.map { it -> + it.pow(2) + it + 1 + chain.nextDouble() + } + val yErr = x.map { sigma } + with(CMFit) { + val chi2 = chiSquared(x.asBuffer(), y.asBuffer(), yErr.asBuffer()) { x -> + bind(a) * x.pow(2) + bind(b) * x + bind(c) + } + + 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