Fitting refactor

This commit is contained in:
Alexander Nozik 2020-10-28 09:08:37 +03:00
parent 1c1580c8e6
commit f8c3d1793c
4 changed files with 105 additions and 109 deletions

View File

@ -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<Double>,
y: Buffer<Double>,
yErr: Buffer<Double>,
model: Expression<Double>,
xSymbol: Symbol = StringSymbol("x"),
): Expression<Double> {
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<Double>,
y: Buffer<Double>,
yErr: Buffer<Double>,
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<Double>.optimize(
vararg symbols: Symbol,
configuration: CMOptimizationProblem.() -> Unit,
): OptimizationResult<Double> = optimizeWith(CMOptimizationProblem, symbols = symbols, configuration)
/**
* Optimize differentiable expression
*/
public fun DifferentiableExpression<Double>.optimize(
vararg symbols: Symbol,
configuration: CMOptimizationProblem.() -> Unit,
): OptimizationResult<Double> = optimizeWith(CMOptimizationProblem, symbols = symbols, configuration)
public fun DifferentiableExpression<Double>.minimize(
vararg startPoint: Pair<Symbol, Double>,
configuration: CMOptimizationProblem.() -> Unit = {},
): OptimizationResult<Double> {
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()
}

View File

@ -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<Double>,
y: Buffer<Double>,
yErr: Buffer<Double>,
model: DerivativeStructureField.(x: DerivativeStructure) -> DerivativeStructure,
): DifferentiableExpression<Double> = 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<Double>,
y: Iterable<Double>,
yErr: Iterable<Double>,
model: DerivativeStructureField.(x: DerivativeStructure) -> DerivativeStructure,
): DifferentiableExpression<Double> = chiSquared(
DerivativeStructureField,
x.toList().asBuffer(),
y.toList().asBuffer(),
yErr.toList().asBuffer(),
model
)
/**
* Optimize expression without derivatives
*/
public fun Expression<Double>.optimize(
vararg symbols: Symbol,
configuration: CMOptimizationProblem.() -> Unit,
): OptimizationResult<Double> = optimizeWith(CMOptimizationProblem, symbols = symbols, configuration)
/**
* Optimize differentiable expression
*/
public fun DifferentiableExpression<Double>.optimize(
vararg symbols: Symbol,
configuration: CMOptimizationProblem.() -> Unit,
): OptimizationResult<Double> = optimizeWith(CMOptimizationProblem, symbols = symbols, configuration)
public fun DifferentiableExpression<Double>.minimize(
vararg startPoint: Pair<Symbol, Double>,
configuration: CMOptimizationProblem.() -> Unit = {},
): OptimizationResult<Double> {
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()
}

View File

@ -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)}")
}
}

View File

@ -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<Double>,
y: Buffer<Double>,
yErr: Buffer<Double>,
model: Expression<Double>,
xSymbol: Symbol = StringSymbol("x"),
): Expression<Double> {
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)
}
}
}
}