From e8dafad6c582646f649f8d94484e9c2643ac0de3 Mon Sep 17 00:00:00 2001 From: Margarita Lashina Date: Wed, 7 Jun 2023 05:25:32 +0300 Subject: [PATCH] the input data is placed in a separate class, to which the documentation is written --- .../StaticLm/staticDifficultTest.kt | 20 +- .../StaticLm/staticEasyTest.kt | 17 +- .../StaticLm/staticMiddleTest.kt | 20 +- .../StreamingLm/streamLm.kt | 32 +-- .../LevenbergMarquardt/functionsToOptimize.kt | 14 +- .../core/LevenbergMarquardtAlgorithm.kt | 193 +++++++++++------- .../kmath/tensors/core/TestLmAlgorithm.kt | 45 ++-- 7 files changed, 197 insertions(+), 144 deletions(-) diff --git a/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StaticLm/staticDifficultTest.kt b/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StaticLm/staticDifficultTest.kt index 95a62e21e..e6f575262 100644 --- a/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StaticLm/staticDifficultTest.kt +++ b/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StaticLm/staticDifficultTest.kt @@ -12,7 +12,8 @@ import space.kscience.kmath.tensors.LevenbergMarquardt.funcDifficultForLm import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra.div import space.kscience.kmath.tensors.core.DoubleTensorAlgebra -import space.kscience.kmath.tensors.core.lm +import space.kscience.kmath.tensors.core.LMInput +import space.kscience.kmath.tensors.core.levenbergMarquardt import kotlin.math.roundToInt fun main() { @@ -39,9 +40,7 @@ fun main() { var t = t_example val y_dat = y_hat - val weight = BroadcastDoubleTensorAlgebra.fromArray( - ShapeND(intArrayOf(1, 1)), DoubleArray(1) { 1.0 / Nparams * 1.0 - 0.085 } - ).as2D() + val weight = 1.0 / Nparams * 1.0 - 0.085 val dp = BroadcastDoubleTensorAlgebra.fromArray( ShapeND(intArrayOf(1, 1)), DoubleArray(1) { -0.01 } ).as2D() @@ -52,8 +51,7 @@ fun main() { val opts = doubleArrayOf(3.0, 10000.0, 1e-6, 1e-6, 1e-6, 1e-6, 1e-2, 11.0, 9.0, 1.0) // val opts = doubleArrayOf(3.0, 10000.0, 1e-6, 1e-6, 1e-6, 1e-6, 1e-3, 11.0, 9.0, 1.0) - val result = DoubleTensorAlgebra.lm( - ::funcDifficultForLm, + val inputData = LMInput(::funcDifficultForLm, p_init.as2D(), t, y_dat, @@ -61,10 +59,14 @@ fun main() { dp, p_min.as2D(), p_max.as2D(), - opts, + opts[1].toInt(), + doubleArrayOf(opts[2], opts[3], opts[4], opts[5]), + doubleArrayOf(opts[6], opts[7], opts[8]), + opts[9].toInt(), 10, - 1 - ) + 1) + + val result = DoubleTensorAlgebra.levenbergMarquardt(inputData) println("Parameters:") for (i in 0 until result.resultParameters.shape.component1()) { diff --git a/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StaticLm/staticEasyTest.kt b/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StaticLm/staticEasyTest.kt index 0b26838a0..507943031 100644 --- a/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StaticLm/staticEasyTest.kt +++ b/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StaticLm/staticEasyTest.kt @@ -12,14 +12,13 @@ import space.kscience.kmath.tensors.LevenbergMarquardt.funcDifficultForLm import space.kscience.kmath.tensors.LevenbergMarquardt.funcEasyForLm import space.kscience.kmath.tensors.LevenbergMarquardt.getStartDataForFuncEasy import space.kscience.kmath.tensors.core.DoubleTensorAlgebra -import space.kscience.kmath.tensors.core.lm +import space.kscience.kmath.tensors.core.LMInput +import space.kscience.kmath.tensors.core.levenbergMarquardt import kotlin.math.roundToInt fun main() { val startedData = getStartDataForFuncEasy() - - val result = DoubleTensorAlgebra.lm( - ::funcEasyForLm, + val inputData = LMInput(::funcEasyForLm, DoubleTensorAlgebra.ones(ShapeND(intArrayOf(4, 1))).as2D(), startedData.t, startedData.y_dat, @@ -27,10 +26,14 @@ fun main() { startedData.dp, startedData.p_min, startedData.p_max, - startedData.opts, + startedData.opts[1].toInt(), + doubleArrayOf(startedData.opts[2], startedData.opts[3], startedData.opts[4], startedData.opts[5]), + doubleArrayOf(startedData.opts[6], startedData.opts[7], startedData.opts[8]), + startedData.opts[9].toInt(), 10, - startedData.example_number - ) + startedData.example_number) + + val result = DoubleTensorAlgebra.levenbergMarquardt(inputData) println("Parameters:") for (i in 0 until result.resultParameters.shape.component1()) { diff --git a/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StaticLm/staticMiddleTest.kt b/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StaticLm/staticMiddleTest.kt index b60ea1897..0659db103 100644 --- a/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StaticLm/staticMiddleTest.kt +++ b/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StaticLm/staticMiddleTest.kt @@ -12,7 +12,8 @@ import space.kscience.kmath.tensors.LevenbergMarquardt.funcMiddleForLm import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra.div import space.kscience.kmath.tensors.core.DoubleTensorAlgebra -import space.kscience.kmath.tensors.core.lm +import space.kscience.kmath.tensors.core.LMInput +import space.kscience.kmath.tensors.core.levenbergMarquardt import kotlin.math.roundToInt fun main() { val NData = 100 @@ -38,9 +39,7 @@ fun main() { var t = t_example val y_dat = y_hat - val weight = BroadcastDoubleTensorAlgebra.fromArray( - ShapeND(intArrayOf(1, 1)), DoubleArray(1) { 1.0 } - ).as2D() + val weight = 1.0 val dp = BroadcastDoubleTensorAlgebra.fromArray( ShapeND(intArrayOf(1, 1)), DoubleArray(1) { -0.01 } ).as2D() @@ -50,8 +49,7 @@ fun main() { p_min = p_min.div(1.0 / 50.0) val opts = doubleArrayOf(3.0, 7000.0, 1e-5, 1e-5, 1e-5, 1e-5, 1e-5, 11.0, 9.0, 1.0) - val result = DoubleTensorAlgebra.lm( - ::funcMiddleForLm, + val inputData = LMInput(::funcMiddleForLm, p_init.as2D(), t, y_dat, @@ -59,10 +57,14 @@ fun main() { dp, p_min.as2D(), p_max.as2D(), - opts, + opts[1].toInt(), + doubleArrayOf(opts[2], opts[3], opts[4], opts[5]), + doubleArrayOf(opts[6], opts[7], opts[8]), + opts[9].toInt(), 10, - 1 - ) + 1) + + val result = DoubleTensorAlgebra.levenbergMarquardt(inputData) println("Parameters:") for (i in 0 until result.resultParameters.shape.component1()) { diff --git a/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StreamingLm/streamLm.kt b/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StreamingLm/streamLm.kt index 99fb97923..fe96b2fe9 100644 --- a/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StreamingLm/streamLm.kt +++ b/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StreamingLm/streamLm.kt @@ -11,7 +11,8 @@ import space.kscience.kmath.nd.* import space.kscience.kmath.tensors.LevenbergMarquardt.StartDataLm import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra.zeros import space.kscience.kmath.tensors.core.DoubleTensorAlgebra -import space.kscience.kmath.tensors.core.lm +import space.kscience.kmath.tensors.core.LMInput +import space.kscience.kmath.tensors.core.levenbergMarquardt import kotlin.random.Random import kotlin.reflect.KFunction3 @@ -31,20 +32,23 @@ fun streamLm(lm_func: KFunction3, MutableStructure2D< var steps = numberOfLaunches val isEndless = (steps <= 0) + val inputData = LMInput(lm_func, + p_init, + t, + y_dat, + weight, + dp, + p_min, + p_max, + opts[1].toInt(), + doubleArrayOf(opts[2], opts[3], opts[4], opts[5]), + doubleArrayOf(opts[6], opts[7], opts[8]), + opts[9].toInt(), + 10, + example_number) + while (isEndless || steps > 0) { - val result = DoubleTensorAlgebra.lm( - lm_func, - p_init, - t, - y_dat, - weight, - dp, - p_min, - p_max, - opts, - 10, - example_number - ) + val result = DoubleTensorAlgebra.levenbergMarquardt(inputData) emit(result.resultParameters) delay(launchFrequencyInMs) p_init = result.resultParameters diff --git a/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/functionsToOptimize.kt b/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/functionsToOptimize.kt index 537b86da3..7ccb37ed0 100644 --- a/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/functionsToOptimize.kt +++ b/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/functionsToOptimize.kt @@ -24,7 +24,7 @@ public data class StartDataLm ( var p_init: MutableStructure2D, var t: MutableStructure2D, var y_dat: MutableStructure2D, - var weight: MutableStructure2D, + var weight: Double, var dp: MutableStructure2D, var p_min: MutableStructure2D, var p_max: MutableStructure2D, @@ -113,9 +113,7 @@ fun getStartDataForFuncDifficult(): StartDataLm { var t = t_example val y_dat = y_hat - val weight = BroadcastDoubleTensorAlgebra.fromArray( - ShapeND(intArrayOf(1, 1)), DoubleArray(1) { 1.0 / Nparams * 1.0 - 0.085 } - ).as2D() + val weight = 1.0 / Nparams * 1.0 - 0.085 val dp = BroadcastDoubleTensorAlgebra.fromArray( ShapeND(intArrayOf(1, 1)), DoubleArray(1) { -0.01 } ).as2D() @@ -154,9 +152,7 @@ fun getStartDataForFuncMiddle(): StartDataLm { } var t = t_example val y_dat = y_hat - val weight = BroadcastDoubleTensorAlgebra.fromArray( - ShapeND(intArrayOf(1, 1)), DoubleArray(1) { 1.0 } - ).as2D() + val weight = 1.0 val dp = BroadcastDoubleTensorAlgebra.fromArray( ShapeND(intArrayOf(1, 1)), DoubleArray(1) { -0.01 } ).as2D() @@ -202,9 +198,7 @@ fun getStartDataForFuncEasy(): StartDataLm { ShapeND(intArrayOf(100, 1)), lm_matx_y_dat ).as2D() - val weight = BroadcastDoubleTensorAlgebra.fromArray( - ShapeND(intArrayOf(1, 1)), DoubleArray(1) { 4.0 } - ).as2D() + val weight = 4.0 val dp = BroadcastDoubleTensorAlgebra.fromArray( ShapeND(intArrayOf(1, 1)), DoubleArray(1) { -0.01 } diff --git a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/LevenbergMarquardtAlgorithm.kt b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/LevenbergMarquardtAlgorithm.kt index af41a7a00..ec9d92c88 100644 --- a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/LevenbergMarquardtAlgorithm.kt +++ b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/LevenbergMarquardtAlgorithm.kt @@ -19,21 +19,21 @@ import kotlin.math.pow import kotlin.reflect.KFunction3 /** - * Type of convergence achieved as a result of executing the Levenberg-Marquardt algorithm + * Type of convergence achieved as a result of executing the Levenberg-Marquardt algorithm. * * InGradient: gradient convergence achieved - * (max(J^T W dy) < epsilon1 = opts[2], + * (max(J^T W dy) < epsilon1, * where J - Jacobi matrix (dy^/dp) for the current approximation y^, - * W - weight matrix from input, dy = (y - y^(p))) + * W - weight matrix from input, dy = (y - y^(p))). * InParameters: convergence in parameters achieved - * (max(h_i / p_i) < epsilon2 = opts[3], - * where h_i - offset for parameter p_i on the current iteration) + * (max(h_i / p_i) < epsilon2, + * where h_i - offset for parameter p_i on the current iteration). * InReducedChiSquare: chi-squared convergence achieved - * (chi squared value divided by (m - n + 1) < epsilon2 = opts[4], - * where n - number of parameters, m - amount of points - * NoConvergence: the maximum number of iterations has been reached without reaching any convergence + * (chi squared value divided by (m - n + 1) < epsilon2, + * where n - number of parameters, m - amount of points). + * NoConvergence: the maximum number of iterations has been reached without reaching any convergence. */ -public enum class TypeOfConvergence{ +public enum class TypeOfConvergence { InGradient, InParameters, InReducedChiSquare, @@ -41,14 +41,14 @@ public enum class TypeOfConvergence{ } /** - * Class for the data obtained as a result of the execution of the Levenberg-Marquardt algorithm + * The data obtained as a result of the execution of the Levenberg-Marquardt algorithm. * - * iterations: number of completed iterations - * funcCalls: the number of evaluations of the input function during execution - * resultChiSq: chi squared value on final parameters - * resultLambda: final lambda parameter used to calculate the offset - * resultParameters: final parameters - * typeOfConvergence: type of convergence + * iterations: number of completed iterations. + * funcCalls: the number of evaluations of the input function during execution. + * resultChiSq: chi squared value on final parameters. + * resultLambda: final lambda parameter used to calculate the offset. + * resultParameters: final parameters. + * typeOfConvergence: type of convergence. */ public data class LMResultInfo ( var iterations:Int, @@ -59,26 +59,65 @@ public data class LMResultInfo ( var typeOfConvergence: TypeOfConvergence, ) -public fun DoubleTensorAlgebra.lm( - func: KFunction3, MutableStructure2D, Int, MutableStructure2D>, - pInput: MutableStructure2D, tInput: MutableStructure2D, yDatInput: MutableStructure2D, - weightInput: MutableStructure2D, dpInput: MutableStructure2D, pMinInput: MutableStructure2D, - pMaxInput: MutableStructure2D, optsInput: DoubleArray, nargin: Int, exampleNumber: Int): LMResultInfo { - +/** + * Input data for the Levenberg-Marquardt function. + * + * func: function of n independent variables x, m parameters an example number, + * rotating a vector of n values y, in which each of the y_i is calculated at its x_i with the given parameters. + * startParameters: starting parameters. + * independentVariables: independent variables, for each of which the real value is known. + * realValues: real values obtained with given independent variables but unknown parameters. + * weight: measurement error for realValues (denominator in each term of sum of weighted squared errors). + * pDelta: delta when calculating the derivative with respect to parameters. + * minParameters: the lower bound of parameter values. + * maxParameters: upper limit of parameter values. + * maxIterations: maximum allowable number of iterations. + * epsilons: epsilon1 - convergence tolerance for gradient, + * epsilon2 - convergence tolerance for parameters, + * epsilon3 - convergence tolerance for reduced chi-square, + * epsilon4 - determines acceptance of a step. + * lambdas: lambda0 - starting lambda value for parameter offset count, + * lambdaUp - factor for increasing lambda, + * lambdaDown - factor for decreasing lambda. + * updateType: 1: Levenberg-Marquardt lambda update, + * 2: Quadratic update, + * 3: Nielsen's lambda update equations. + * nargin: a value that determines which options to use by default + * (<5 - use weight by default, <6 - use pDelta by default, <7 - use minParameters by default, + * <8 - use maxParameters by default, <9 - use updateType by default). + * exampleNumber: a parameter for a function with which you can choose its behavior. + */ +public data class LMInput ( + var func: KFunction3, MutableStructure2D, Int, MutableStructure2D>, + var startParameters: MutableStructure2D, + var independentVariables: MutableStructure2D, + var realValues: MutableStructure2D, + var weight: Double, + var pDelta: MutableStructure2D, + var minParameters: MutableStructure2D, + var maxParameters: MutableStructure2D, + var maxIterations: Int, + var epsilons: DoubleArray, + var lambdas: DoubleArray, + var updateType: Int, + var nargin: Int, + var exampleNumber: Int +) +public fun DoubleTensorAlgebra.levenbergMarquardt(inputData: LMInput): LMResultInfo { val resultInfo = LMResultInfo(0, 0, 0.0, - 0.0, pInput, TypeOfConvergence.NoConvergence) + 0.0, inputData.startParameters, TypeOfConvergence.NoConvergence) val eps = 2.2204e-16 - val settings = LMSettings(0, 0, exampleNumber) + val settings = LMSettings(0, 0, inputData.exampleNumber) settings.funcCalls = 0 // running count of function evaluations - var p = pInput - val t = tInput + var p = inputData.startParameters + val t = inputData.independentVariables val Npar = length(p) // number of parameters - val Npnt = length(yDatInput) // number of data points + val Npnt = length(inputData.realValues) // number of data points var pOld = zeros(ShapeND(intArrayOf(Npar, 1))).as2D() // previous set of parameters var yOld = zeros(ShapeND(intArrayOf(Npnt, 1))).as2D() // previous model, y_old = y_hat(t;p_old) var X2 = 1e-3 / eps // a really big initial Chi-sq value @@ -86,50 +125,55 @@ public fun DoubleTensorAlgebra.lm( var J = zeros(ShapeND(intArrayOf(Npnt, Npar))).as2D() // Jacobian matrix val DoF = Npnt - Npar // statistical degrees of freedom - var weight = weightInput - if (nargin < 5) { - weight = fromArray(ShapeND(intArrayOf(1, 1)), doubleArrayOf((yDatInput.transpose().dot(yDatInput)).as1D()[0])).as2D() + var weight = fromArray(ShapeND(intArrayOf(1, 1)), doubleArrayOf(inputData.weight)).as2D() + if (inputData.nargin < 5) { + weight = fromArray(ShapeND(intArrayOf(1, 1)), doubleArrayOf((inputData.realValues.transpose().dot(inputData.realValues)).as1D()[0])).as2D() } - var dp = dpInput - if (nargin < 6) { + var dp = inputData.pDelta + if (inputData.nargin < 6) { dp = fromArray(ShapeND(intArrayOf(1, 1)), doubleArrayOf(0.001)).as2D() } - var pMin = pMinInput - if (nargin < 7) { - pMin = p - pMin.abs() - pMin = pMin.div(-100.0).as2D() + var minParameters = inputData.minParameters + if (inputData.nargin < 7) { + minParameters = p + minParameters.abs() + minParameters = minParameters.div(-100.0).as2D() } - var pMax = pMaxInput - if (nargin < 8) { - pMax = p - pMax.abs() - pMax = pMax.div(100.0).as2D() + var maxParameters = inputData.maxParameters + if (inputData.nargin < 8) { + maxParameters = p + maxParameters.abs() + maxParameters = maxParameters.div(100.0).as2D() } - var opts = optsInput - if (nargin < 10) { - opts = doubleArrayOf(3.0, 10.0 * Npar, 1e-3, 1e-3, 1e-1, 1e-1, 1e-2, 11.0, 9.0, 1.0) + var maxIterations = inputData.maxIterations + var epsilon1 = inputData.epsilons[0] // convergence tolerance for gradient + var epsilon2 = inputData.epsilons[1] // convergence tolerance for parameters + var epsilon3 = inputData.epsilons[2] // convergence tolerance for Chi-square + var epsilon4 = inputData.epsilons[3] // determines acceptance of a L-M step + var lambda0 = inputData.lambdas[0] // initial value of damping paramter, lambda + var lambdaUpFac = inputData.lambdas[1] // factor for increasing lambda + var lambdaDnFac = inputData.lambdas[2] // factor for decreasing lambda + var updateType = inputData.updateType // 1: Levenberg-Marquardt lambda update + // 2: Quadratic update + // 3: Nielsen's lambda update equations + if (inputData.nargin < 9) { + maxIterations = 10 * Npar + epsilon1 = 1e-3 + epsilon2 = 1e-3 + epsilon3 = 1e-1 + epsilon4 = 1e-1 + lambda0 = 1e-2 + lambdaUpFac = 11.0 + lambdaDnFac = 9.0 + updateType = 1 } - val prnt = opts[0] // >1 intermediate results; >2 plots - val maxIterations = opts[1].toInt() // maximum number of iterations - val epsilon1 = opts[2] // convergence tolerance for gradient - val epsilon2 = opts[3] // convergence tolerance for parameters - val epsilon3 = opts[4] // convergence tolerance for Chi-square - val epsilon4 = opts[5] // determines acceptance of a L-M step - val lambda0 = opts[6] // initial value of damping paramter, lambda - val lambdaUpFac = opts[7] // factor for increasing lambda - val lambdaDnFac = opts[8] // factor for decreasing lambda - val updateType = opts[9].toInt() // 1: Levenberg-Marquardt lambda update - // 2: Quadratic update - // 3: Nielsen's lambda update equations - - pMin = makeColumn(pMin) - pMax = makeColumn(pMax) + minParameters = makeColumn(minParameters) + maxParameters = makeColumn(maxParameters) if (length(makeColumn(dp)) == 1) { dp = ones(ShapeND(intArrayOf(Npar, 1))).div(1 / dp[0, 0]).as2D() @@ -146,7 +190,7 @@ public fun DoubleTensorAlgebra.lm( } // initialize Jacobian with finite difference calculation - var lmMatxAns = lmMatx(func, t, pOld, yOld, 1, J, p, yDatInput, weight, dp, settings) + var lmMatxAns = lmMatx(inputData.func, t, pOld, yOld, 1, J, p, inputData.realValues, weight, dp, settings) var JtWJ = lmMatxAns[0] var JtWdy = lmMatxAns[1] X2 = lmMatxAns[2][0, 0] @@ -189,9 +233,9 @@ public fun DoubleTensorAlgebra.lm( } var pTry = (p + h).as2D() // update the [idx] elements - pTry = smallestElementComparison(largestElementComparison(pMin, pTry.as2D()), pMax) // apply constraints + pTry = smallestElementComparison(largestElementComparison(minParameters, pTry.as2D()), maxParameters) // apply constraints - var deltaY = yDatInput.minus(evaluateFunction(func, t, pTry, exampleNumber)) // residual error using p_try + var deltaY = inputData.realValues.minus(evaluateFunction(inputData.func, t, pTry, inputData.exampleNumber)) // residual error using p_try for (i in 0 until deltaY.shape.component1()) { // floating point error; break for (j in 0 until deltaY.shape.component2()) { @@ -214,9 +258,9 @@ public fun DoubleTensorAlgebra.lm( val alpha = JtWdy.transpose().dot(h) / ((X2Try.minus(X2)).div(2.0).plus(2 * JtWdy.transpose().dot(h))) h = h.dot(alpha) pTry = p.plus(h).as2D() // update only [idx] elements - pTry = smallestElementComparison(largestElementComparison(pMin, pTry), pMax) // apply constraints + pTry = smallestElementComparison(largestElementComparison(minParameters, pTry), maxParameters) // apply constraints - deltaY = yDatInput.minus(evaluateFunction(func, t, pTry, exampleNumber)) // residual error using p_try + deltaY = inputData.realValues.minus(evaluateFunction(inputData.func, t, pTry, inputData.exampleNumber)) // residual error using p_try settings.funcCalls += 1 X2Try = deltaY.as2D().transpose().dot(deltaY.times(weight)) // Chi-squared error criteria @@ -242,7 +286,7 @@ public fun DoubleTensorAlgebra.lm( yOld = yHat.copyToTensor().as2D() p = makeColumn(pTry) // accept p_try - lmMatxAns = lmMatx(func, t, pOld, yOld, dX2.toInt(), J, p, yDatInput, weight, dp, settings) + lmMatxAns = lmMatx(inputData.func, t, pOld, yOld, dX2.toInt(), J, p, inputData.realValues, weight, dp, settings) // decrease lambda ==> Gauss-Newton method JtWJ = lmMatxAns[0] @@ -268,7 +312,7 @@ public fun DoubleTensorAlgebra.lm( } else { // it IS NOT better X2 = X2Old // do not accept p_try if (settings.iteration % (2 * Npar) == 0) { // rank-1 update of Jacobian - lmMatxAns = lmMatx(func, t, pOld, yOld, -1, J, p, yDatInput, weight, dp, settings) + lmMatxAns = lmMatx(inputData.func, t, pOld, yOld, -1, J, p, inputData.realValues, weight, dp, settings) JtWJ = lmMatxAns[0] JtWdy = lmMatxAns[1] yHat = lmMatxAns[3] @@ -292,14 +336,13 @@ public fun DoubleTensorAlgebra.lm( } } - if (prnt > 1) { - val chiSq = X2 / DoF - resultInfo.iterations = settings.iteration - resultInfo.funcCalls = settings.funcCalls - resultInfo.resultChiSq = chiSq - resultInfo.resultLambda = lambda - resultInfo.resultParameters = p - } + val chiSq = X2 / DoF + resultInfo.iterations = settings.iteration + resultInfo.funcCalls = settings.funcCalls + resultInfo.resultChiSq = chiSq + resultInfo.resultLambda = lambda + resultInfo.resultParameters = p + if (abs(JtWdy).max() < epsilon1 && settings.iteration > 2) { resultInfo.typeOfConvergence = TypeOfConvergence.InGradient diff --git a/kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestLmAlgorithm.kt b/kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestLmAlgorithm.kt index 114e5f879..4e1df3b08 100644 --- a/kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestLmAlgorithm.kt +++ b/kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestLmAlgorithm.kt @@ -105,9 +105,7 @@ class TestLmAlgorithm { ShapeND(intArrayOf(100, 1)), lm_matx_y_dat ).as2D() - val weight = BroadcastDoubleTensorAlgebra.fromArray( - ShapeND(intArrayOf(1, 1)), DoubleArray(1) { 4.0 } - ).as2D() + val weight = 4.0 val dp = BroadcastDoubleTensorAlgebra.fromArray( ShapeND(intArrayOf(1, 1)), DoubleArray(1) { -0.01 } @@ -123,7 +121,12 @@ class TestLmAlgorithm { val opts = doubleArrayOf(3.0, 100.0, 1e-3, 1e-3, 1e-1, 1e-1, 1e-2, 11.0, 9.0, 1.0) - val result = lm(::funcEasyForLm, p_init, t, y_dat, weight, dp, p_min, p_max, opts, 10, example_number) + val inputData = LMInput(::funcEasyForLm, p_init, t, y_dat, weight, dp, p_min, p_max, opts[1].toInt(), + doubleArrayOf(opts[2], opts[3], opts[4], opts[5]), + doubleArrayOf(opts[6], opts[7], opts[8]), + opts[9].toInt(), 10, example_number) + + val result = levenbergMarquardt(inputData) assertEquals(13, result.iterations) assertEquals(31, result.funcCalls) assertEquals(0.9131368192633, (result.resultChiSq * 1e13).roundToLong() / 1e13) @@ -168,9 +171,7 @@ class TestLmAlgorithm { var t = t_example val y_dat = y_hat - val weight = BroadcastDoubleTensorAlgebra.fromArray( - ShapeND(intArrayOf(1, 1)), DoubleArray(1) { 1.0 } - ).as2D() + val weight = 1.0 val dp = BroadcastDoubleTensorAlgebra.fromArray( ShapeND(intArrayOf(1, 1)), DoubleArray(1) { -0.01 } ).as2D() @@ -180,8 +181,7 @@ class TestLmAlgorithm { p_min = p_min.div(1.0 / 50.0) val opts = doubleArrayOf(3.0, 7000.0, 1e-5, 1e-5, 1e-5, 1e-5, 1e-5, 11.0, 9.0, 1.0) - val result = DoubleTensorAlgebra.lm( - ::funcMiddleForLm, + val inputData = LMInput(::funcMiddleForLm, p_init.as2D(), t, y_dat, @@ -189,10 +189,14 @@ class TestLmAlgorithm { dp, p_min.as2D(), p_max.as2D(), - opts, + opts[1].toInt(), + doubleArrayOf(opts[2], opts[3], opts[4], opts[5]), + doubleArrayOf(opts[6], opts[7], opts[8]), + opts[9].toInt(), 10, - 1 - ) + 1) + + val result = DoubleTensorAlgebra.levenbergMarquardt(inputData) } @Test @@ -220,9 +224,7 @@ class TestLmAlgorithm { var t = t_example val y_dat = y_hat - val weight = BroadcastDoubleTensorAlgebra.fromArray( - ShapeND(intArrayOf(1, 1)), DoubleArray(1) { 1.0 / Nparams * 1.0 - 0.085 } - ).as2D() + val weight = 1.0 / Nparams * 1.0 - 0.085 val dp = BroadcastDoubleTensorAlgebra.fromArray( ShapeND(intArrayOf(1, 1)), DoubleArray(1) { -0.01 } ).as2D() @@ -232,8 +234,7 @@ class TestLmAlgorithm { p_min = p_min.div(1.0 / 50.0) val opts = doubleArrayOf(3.0, 7000.0, 1e-2, 1e-3, 1e-2, 1e-2, 1e-2, 11.0, 9.0, 1.0) - val result = DoubleTensorAlgebra.lm( - ::funcDifficultForLm, + val inputData = LMInput(::funcDifficultForLm, p_init.as2D(), t, y_dat, @@ -241,9 +242,13 @@ class TestLmAlgorithm { dp, p_min.as2D(), p_max.as2D(), - opts, + opts[1].toInt(), + doubleArrayOf(opts[2], opts[3], opts[4], opts[5]), + doubleArrayOf(opts[6], opts[7], opts[8]), + opts[9].toInt(), 10, - 1 - ) + 1) + + val result = DoubleTensorAlgebra.levenbergMarquardt(inputData) } } \ No newline at end of file