diff --git a/build.gradle.kts b/build.gradle.kts index aed79909c..7dbe87445 100644 --- a/build.gradle.kts +++ b/build.gradle.kts @@ -15,7 +15,7 @@ allprojects { } group = "space.kscience" - version = "0.3.1" + version = "0.3.2-dev-1" } subprojects { 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 3cb485d7d..72752f855 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 @@ -5,18 +5,13 @@ package space.kscience.kmath.tensors.core +import space.kscience.kmath.PerformancePitfall import space.kscience.kmath.linear.transpose import space.kscience.kmath.nd.* -import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra.div -import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra.dot -import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra.minus -import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra.times -import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra.transposed -import space.kscience.kmath.tensors.core.DoubleTensorAlgebra.Companion.plus +import kotlin.math.abs import kotlin.math.max import kotlin.math.min import kotlin.math.pow -import kotlin.reflect.KFunction3 /** * Type of convergence achieved as a result of executing the Levenberg-Marquardt algorithm. @@ -50,8 +45,8 @@ public enum class TypeOfConvergence { * resultParameters: final parameters. * typeOfConvergence: type of convergence. */ -public data class LMResultInfo ( - var iterations:Int, +public data class LMResultInfo( + var iterations: Int, var funcCalls: Int, var resultChiSq: Double, var resultLambda: Double, @@ -87,7 +82,7 @@ public data class LMResultInfo ( * <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 ( +public data class LMInput( var func: (MutableStructure2D, MutableStructure2D, Int) -> (MutableStructure2D), var startParameters: MutableStructure2D, var independentVariables: MutableStructure2D, @@ -101,7 +96,7 @@ public data class LMInput ( var lambdas: DoubleArray, var updateType: Int, var nargin: Int, - var exampleNumber: Int + var exampleNumber: Int, ) @@ -120,8 +115,10 @@ public data class LMInput ( * @return the 'output'. */ public fun DoubleTensorAlgebra.levenbergMarquardt(inputData: LMInput): LMResultInfo { - val resultInfo = LMResultInfo(0, 0, 0.0, - 0.0, inputData.startParameters, TypeOfConvergence.NoConvergence) + val resultInfo = LMResultInfo( + 0, 0, 0.0, + 0.0, inputData.startParameters, TypeOfConvergence.NoConvergence + ) val eps = 2.2204e-16 @@ -131,18 +128,21 @@ public fun DoubleTensorAlgebra.levenbergMarquardt(inputData: LMInput): LMResultI var p = inputData.startParameters val t = inputData.independentVariables - val Npar = length(p) // number of parameters - 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 - var X2Old = 1e-3 / eps // a really big initial Chi-sq value - var J = zeros(ShapeND(intArrayOf(Npnt, Npar))).as2D() // Jacobian matrix - val DoF = Npnt - Npar // statistical degrees of freedom + val nPar = length(p) // number of parameters + val nPoints = length(inputData.realValues) // number of data points + var pOld = zeros(ShapeND(intArrayOf(nPar, 1))).as2D() // previous set of parameters + var yOld = zeros(ShapeND(intArrayOf(nPoints, 1))).as2D() // previous model, y_old = y_hat(t;p_old) + var x2: Double // a really big initial Chi-sq value + var x2Old: Double // a really big initial Chi-sq value + var jacobian = zeros(ShapeND(intArrayOf(nPoints, nPar))).as2D() // Jacobian matrix + val dof = nPoints - nPar // statistical degrees of freedom 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() + if (inputData.nargin < 5) { + weight = fromArray( + ShapeND(intArrayOf(1, 1)), + doubleArrayOf((inputData.realValues.transpose().dot(inputData.realValues)).as1D()[0]) + ).as2D() } var dp = inputData.pDelta @@ -169,13 +169,13 @@ public fun DoubleTensorAlgebra.levenbergMarquardt(inputData: LMInput): LMResultI var epsilon2 = inputData.epsilons[1] var epsilon3 = inputData.epsilons[2] var epsilon4 = inputData.epsilons[3] - var lambda0 = inputData.lambdas[0] + var lambda0 = inputData.lambdas[0] var lambdaUpFac = inputData.lambdas[1] var lambdaDnFac = inputData.lambdas[2] var updateType = inputData.updateType if (inputData.nargin < 9) { - maxIterations = 10 * Npar + maxIterations = 10 * nPar epsilon1 = 1e-3 epsilon2 = 1e-3 epsilon3 = 1e-1 @@ -190,28 +190,27 @@ public fun DoubleTensorAlgebra.levenbergMarquardt(inputData: LMInput): LMResultI maxParameters = makeColumn(maxParameters) if (length(makeColumn(dp)) == 1) { - dp = ones(ShapeND(intArrayOf(Npar, 1))).div(1 / dp[0, 0]).as2D() + dp = ones(ShapeND(intArrayOf(nPar, 1))).div(1 / dp[0, 0]).as2D() } var stop = false // termination flag if (weight.shape.component1() == 1 || variance(weight) == 0.0) { // identical weights vector - weight = ones(ShapeND(intArrayOf(Npnt, 1))).div(1 / kotlin.math.abs(weight[0, 0])).as2D() - } - else { + weight = ones(ShapeND(intArrayOf(nPoints, 1))).div(1 / abs(weight[0, 0])).as2D() + } else { weight = makeColumn(weight) weight.abs() } // initialize Jacobian with finite difference calculation - 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] + var lmMatxAns = lmMatx(inputData.func, t, pOld, yOld, 1, jacobian, p, inputData.realValues, weight, dp, settings) + var jtWJ = lmMatxAns[0] + var jtWdy = lmMatxAns[1] + x2 = lmMatxAns[2][0, 0] var yHat = lmMatxAns[3] - J = lmMatxAns[4] + jacobian = lmMatxAns[4] - if ( abs(JtWdy).max() < epsilon1 ) { + if (abs(jtWdy).max() < epsilon1) { stop = true } @@ -219,14 +218,13 @@ public fun DoubleTensorAlgebra.levenbergMarquardt(inputData: LMInput): LMResultI var nu = 1 if (updateType == 1) { - lambda = lambda0 // Marquardt: init'l lambda - } - else { - lambda = lambda0 * (makeColumnFromDiagonal(JtWJ)).max() + lambda = lambda0 // Marquardt: init'l lambda + } else { + lambda = lambda0 * (makeColumnFromDiagonal(jtWJ)).max() nu = 2 } - X2Old = X2 // previous value of X2 + x2Old = x2 // previous value of X2 var h: DoubleTensor @@ -235,20 +233,31 @@ public fun DoubleTensorAlgebra.levenbergMarquardt(inputData: LMInput): LMResultI // incremental change in parameters h = if (updateType == 1) { // Marquardt - val solve = solve(JtWJ.plus(makeMatrixWithDiagonal(makeColumnFromDiagonal(JtWJ)).div(1 / lambda)).as2D(), JtWdy) + val solve = + solve((jtWJ + makeMatrixWithDiagonal(makeColumnFromDiagonal(jtWJ)) / (1 / lambda)).as2D(), jtWdy) solve.asDoubleTensor() } else { // Quadratic and Nielsen - val solve = solve(JtWJ.plus(lmEye(Npar).div(1 / lambda)).as2D(), JtWdy) + val solve = solve(jtWJ.plus(lmEye(nPar) * lambda).as2D(), jtWdy) solve.asDoubleTensor() } var pTry = (p + h).as2D() // update the [idx] elements - pTry = smallestElementComparison(largestElementComparison(minParameters, pTry.as2D()), maxParameters) // apply constraints + pTry = smallestElementComparison( + largestElementComparison(minParameters, pTry), + maxParameters + ) // apply constraints - var deltaY = inputData.realValues.minus(evaluateFunction(inputData.func, t, pTry, inputData.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()) { + for (i in 0 until deltaY.shape[0]) { // floating point error; break + for (j in 0 until deltaY.shape[1]) { if (deltaY[i, j] == Double.POSITIVE_INFINITY || deltaY[i, j] == Double.NEGATIVE_INFINITY) { stop = true break @@ -258,49 +267,72 @@ public fun DoubleTensorAlgebra.levenbergMarquardt(inputData: LMInput): LMResultI settings.funcCalls += 1 - val tmp = deltaY.times(weight) - var X2Try = deltaY.as2D().transpose().dot(tmp) // Chi-squared error criteria +// val tmp = deltaY.times(weight) + var X2Try = deltaY.as2D().transpose().dot(deltaY.times(weight)) // Chi-squared error criteria val alpha = 1.0 if (updateType == 2) { // Quadratic // One step of quadratic line update in the h direction for minimum X2 - 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(minParameters, pTry), maxParameters) // apply constraints + val alpha = (jtWdy.transpose() dot h) / ((X2Try - x2) / 2.0 + 2 * (jtWdy.transpose() dot h)) + h = h dot alpha + pTry = (p + h).as2D() // update only [idx] elements + pTry = smallestElementComparison( + largestElementComparison(minParameters, pTry), + maxParameters + ) // apply constraints - deltaY = inputData.realValues.minus(evaluateFunction(inputData.func, t, pTry, inputData.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 + X2Try = deltaY.as2D().transpose() dot deltaY * weight // Chi-squared error criteria } val rho = when (updateType) { // Nielsen 1 -> { val tmp = h.transposed() - .dot(makeMatrixWithDiagonal(makeColumnFromDiagonal(JtWJ)).div(1 / lambda).dot(h).plus(JtWdy)) - X2.minus(X2Try).as2D()[0, 0] / abs(tmp.as2D()).as2D()[0, 0] + .dot((makeMatrixWithDiagonal(makeColumnFromDiagonal(jtWJ)) * lambda dot h) + jtWdy) + (x2 - X2Try)[0, 0] / abs(tmp.as2D())[0, 0] } + else -> { - val tmp = h.transposed().dot(h.div(1 / lambda).plus(JtWdy)) - X2.minus(X2Try).as2D()[0, 0] / abs(tmp.as2D()).as2D()[0, 0] + val tmp = h.transposed().dot((h * lambda) + jtWdy) + x2.minus(X2Try).as2D()[0, 0] / abs(tmp.as2D()).as2D()[0, 0] } } if (rho > epsilon4) { // it IS significantly better - val dX2 = X2.minus(X2Old) - X2Old = X2 + val dX2 = x2.minus(x2Old) + x2Old = x2 pOld = p.copyToTensor().as2D() yOld = yHat.copyToTensor().as2D() p = makeColumn(pTry) // accept p_try - lmMatxAns = lmMatx(inputData.func, t, pOld, yOld, dX2.toInt(), J, p, inputData.realValues, weight, dp, settings) + lmMatxAns = lmMatx( + inputData.func, + t, + pOld, + yOld, + dX2.toInt(), + jacobian, + p, + inputData.realValues, + weight, + dp, + settings + ) // decrease lambda ==> Gauss-Newton method - JtWJ = lmMatxAns[0] - JtWdy = lmMatxAns[1] - X2 = lmMatxAns[2][0, 0] + jtWJ = lmMatxAns[0] + jtWdy = lmMatxAns[1] + x2 = lmMatxAns[2][0, 0] yHat = lmMatxAns[3] - J = lmMatxAns[4] + jacobian = lmMatxAns[4] lambda = when (updateType) { 1 -> { // Levenberg @@ -317,13 +349,14 @@ public fun DoubleTensorAlgebra.levenbergMarquardt(inputData: LMInput): LMResultI } } } 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(inputData.func, t, pOld, yOld, -1, J, p, inputData.realValues, weight, dp, settings) - JtWJ = lmMatxAns[0] - JtWdy = lmMatxAns[1] + x2 = x2Old // do not accept p_try + if (settings.iteration % (2 * nPar) == 0) { // rank-1 update of Jacobian + lmMatxAns = + lmMatx(inputData.func, t, pOld, yOld, -1, jacobian, p, inputData.realValues, weight, dp, settings) + jtWJ = lmMatxAns[0] + jtWdy = lmMatxAns[1] yHat = lmMatxAns[3] - J = lmMatxAns[4] + jacobian = lmMatxAns[4] } // increase lambda ==> gradient descent method @@ -333,7 +366,7 @@ public fun DoubleTensorAlgebra.levenbergMarquardt(inputData: LMInput): LMResultI } 2 -> { // Quadratic - lambda + kotlin.math.abs(((X2Try.as2D()[0, 0] - X2) / 2) / alpha) + lambda + abs(((X2Try.as2D()[0, 0] - x2) / 2) / alpha) } else -> { // Nielsen @@ -343,7 +376,7 @@ public fun DoubleTensorAlgebra.levenbergMarquardt(inputData: LMInput): LMResultI } } - val chiSq = X2 / DoF + val chiSq = x2 / dof resultInfo.iterations = settings.iteration resultInfo.funcCalls = settings.funcCalls resultInfo.resultChiSq = chiSq @@ -351,7 +384,7 @@ public fun DoubleTensorAlgebra.levenbergMarquardt(inputData: LMInput): LMResultI resultInfo.resultParameters = p - if (abs(JtWdy).max() < epsilon1 && settings.iteration > 2) { + if (abs(jtWdy).max() < epsilon1 && settings.iteration > 2) { resultInfo.typeOfConvergence = TypeOfConvergence.InGradient stop = true } @@ -359,7 +392,7 @@ public fun DoubleTensorAlgebra.levenbergMarquardt(inputData: LMInput): LMResultI resultInfo.typeOfConvergence = TypeOfConvergence.InParameters stop = true } - if (X2 / DoF < epsilon3 && settings.iteration > 2) { + if (x2 / dof < epsilon3 && settings.iteration > 2) { resultInfo.typeOfConvergence = TypeOfConvergence.InReducedChiSquare stop = true } @@ -371,10 +404,10 @@ public fun DoubleTensorAlgebra.levenbergMarquardt(inputData: LMInput): LMResultI return resultInfo } -private data class LMSettings ( - var iteration:Int, +private data class LMSettings( + var iteration: Int, var funcCalls: Int, - var exampleNumber:Int + var exampleNumber: Int, ) /* matrix -> column of all elements */ @@ -390,14 +423,14 @@ private fun makeColumn(tensor: MutableStructure2D): MutableStructure2D) : Int { +private fun length(column: MutableStructure2D): Int { return column.shape.component1() } private fun MutableStructure2D.abs() { - for (i in 0 until this.shape.component1()) { - for (j in 0 until this.shape.component2()) { - this[i, j] = kotlin.math.abs(this[i, j]) + for (i in 0 until this.shape[0]) { + for (j in 0 until this.shape[1]) { + this[i, j] = abs(this[i, j]) } } } @@ -413,7 +446,7 @@ private fun abs(input: MutableStructure2D): MutableStructure2D { ).as2D() for (i in 0 until tensor.shape.component1()) { for (j in 0 until tensor.shape.component2()) { - tensor[i, j] = kotlin.math.abs(input[i, j]) + tensor[i, j] = abs(input[i, j]) } } return tensor @@ -441,21 +474,23 @@ private fun lmEye(size: Int): MutableStructure2D { return makeMatrixWithDiagonal(column) } -private fun largestElementComparison(a: MutableStructure2D, b: MutableStructure2D): MutableStructure2D { +private fun largestElementComparison( + a: MutableStructure2D, + b: MutableStructure2D, +): MutableStructure2D { val aSizeX = a.shape.component1() val aSizeY = a.shape.component2() val bSizeX = b.shape.component1() val bSizeY = b.shape.component2() - val tensor = BroadcastDoubleTensorAlgebra.zeros(ShapeND(intArrayOf(max(aSizeX, bSizeX), max(aSizeY, bSizeY)))).as2D() + val tensor = + BroadcastDoubleTensorAlgebra.zeros(ShapeND(intArrayOf(max(aSizeX, bSizeX), max(aSizeY, bSizeY)))).as2D() for (i in 0 until tensor.shape.component1()) { for (j in 0 until tensor.shape.component2()) { if (i < aSizeX && i < bSizeX && j < aSizeY && j < bSizeY) { tensor[i, j] = max(a[i, j], b[i, j]) - } - else if (i < aSizeX && j < aSizeY) { + } else if (i < aSizeX && j < aSizeY) { tensor[i, j] = a[i, j] - } - else { + } else { tensor[i, j] = b[i, j] } } @@ -463,21 +498,23 @@ private fun largestElementComparison(a: MutableStructure2D, b: MutableSt return tensor } -private fun smallestElementComparison(a: MutableStructure2D, b: MutableStructure2D): MutableStructure2D { +private fun smallestElementComparison( + a: MutableStructure2D, + b: MutableStructure2D, +): MutableStructure2D { val aSizeX = a.shape.component1() val aSizeY = a.shape.component2() val bSizeX = b.shape.component1() val bSizeY = b.shape.component2() - val tensor = BroadcastDoubleTensorAlgebra.zeros(ShapeND(intArrayOf(max(aSizeX, bSizeX), max(aSizeY, bSizeY)))).as2D() + val tensor = + BroadcastDoubleTensorAlgebra.zeros(ShapeND(intArrayOf(max(aSizeX, bSizeX), max(aSizeY, bSizeY)))).as2D() for (i in 0 until tensor.shape.component1()) { for (j in 0 until tensor.shape.component2()) { if (i < aSizeX && i < bSizeX && j < aSizeY && j < bSizeY) { tensor[i, j] = min(a[i, j], b[i, j]) - } - else if (i < aSizeX && j < aSizeY) { + } else if (i < aSizeX && j < aSizeY) { tensor[i, j] = a[i, j] - } - else { + } else { tensor[i, j] = b[i, j] } } @@ -485,10 +522,13 @@ private fun smallestElementComparison(a: MutableStructure2D, b: MutableS return tensor } -private fun getZeroIndices(column: MutableStructure2D, epsilon: Double = 0.000001): MutableStructure2D? { +private fun getZeroIndices( + column: MutableStructure2D, + epsilon: Double = 0.000001, +): MutableStructure2D? { var idx = emptyArray() for (i in 0 until column.shape.component1()) { - if (kotlin.math.abs(column[i, 0]) > epsilon) { + if (abs(column[i, 0]) > epsilon) { idx += (i + 1.0) } } @@ -498,18 +538,27 @@ private fun getZeroIndices(column: MutableStructure2D, epsilon: Double = return null } -private fun evaluateFunction(func: (MutableStructure2D, MutableStructure2D, Int) -> MutableStructure2D, - t: MutableStructure2D, p: MutableStructure2D, exampleNumber: Int) - : MutableStructure2D -{ +private fun evaluateFunction( + func: (MutableStructure2D, MutableStructure2D, Int) -> MutableStructure2D, + t: MutableStructure2D, p: MutableStructure2D, exampleNumber: Int, +) + : MutableStructure2D { return func(t, p, exampleNumber) } -private fun lmMatx(func: (MutableStructure2D, MutableStructure2D, Int) -> MutableStructure2D, - t: MutableStructure2D, pOld: MutableStructure2D, yOld: MutableStructure2D, - dX2: Int, JInput: MutableStructure2D, p: MutableStructure2D, - yDat: MutableStructure2D, weight: MutableStructure2D, dp:MutableStructure2D, settings:LMSettings) : Array> -{ +private fun lmMatx( + func: (MutableStructure2D, MutableStructure2D, Int) -> MutableStructure2D, + t: MutableStructure2D, + pOld: MutableStructure2D, + yOld: MutableStructure2D, + dX2: Int, + JInput: MutableStructure2D, + p: MutableStructure2D, + yDat: MutableStructure2D, + weight: MutableStructure2D, + dp: MutableStructure2D, + settings: LMSettings, +): Array> = with(DoubleTensorAlgebra) { // default: dp = 0.001 val Npar = length(p) // number of parameters @@ -520,63 +569,70 @@ private fun lmMatx(func: (MutableStructure2D, MutableStructure2D J = if (settings.iteration % (2 * Npar) == 0 || dX2 > 0) { lmFdJ(func, t, p, yHat, dp, settings).as2D() // finite difference - } - else { + } else { lmBroydenJ(pOld, yOld, J, p, yHat).as2D() // rank-1 update } val deltaY = yDat.minus(yHat) - val chiSq = deltaY.transposed().dot( deltaY.times(weight) ).as2D() - val JtWJ = J.transposed().dot ( J.times( weight.dot(BroadcastDoubleTensorAlgebra.ones(ShapeND(intArrayOf(1, Npar)))) ) ).as2D() - val JtWdy = J.transposed().dot( weight.times(deltaY) ).as2D() + val chiSq = deltaY.transposed().dot(deltaY.times(weight)).as2D() + val JtWJ = + (J.transposed() dot J * (weight dot ones(ShapeND(intArrayOf(1, Npar))))).as2D() + val JtWdy = (J.transposed() dot weight * deltaY).as2D() - return arrayOf(JtWJ,JtWdy,chiSq,yHat,J) + return arrayOf(JtWJ, JtWdy, chiSq, yHat, J) } -private fun lmBroydenJ(pOld: MutableStructure2D, yOld: MutableStructure2D, JInput: MutableStructure2D, - p: MutableStructure2D, y: MutableStructure2D): MutableStructure2D { +private fun lmBroydenJ( + pOld: MutableStructure2D, yOld: MutableStructure2D, JInput: MutableStructure2D, + p: MutableStructure2D, y: MutableStructure2D, +): MutableStructure2D = with(DoubleTensorAlgebra) { var J = JInput.copyToTensor() val h = p.minus(pOld) - val increase = y.minus(yOld).minus( J.dot(h) ).dot(h.transposed()).div( (h.transposed().dot(h)).as2D()[0, 0] ) + val increase = ((y - yOld - (J dot h)) dot h.transposed()) / (h.transposed() dot h)[0, 0] J = J.plus(increase) return J.as2D() } -private fun lmFdJ(func: (MutableStructure2D, MutableStructure2D, exampleNumber: Int) -> MutableStructure2D, - t: MutableStructure2D, p: MutableStructure2D, y: MutableStructure2D, - dp: MutableStructure2D, settings: LMSettings): MutableStructure2D { +@OptIn(PerformancePitfall::class) +private fun lmFdJ( + func: (MutableStructure2D, MutableStructure2D, exampleNumber: Int) -> MutableStructure2D, + t: MutableStructure2D, + p: MutableStructure2D, + y: MutableStructure2D, + dp: MutableStructure2D, + settings: LMSettings, +): MutableStructure2D = with(DoubleTensorAlgebra) { // default: dp = 0.001 * ones(1,n) val m = length(y) // number of data points val n = length(p) // number of parameters - val ps = p.copyToTensor().as2D() - val J = BroadcastDoubleTensorAlgebra.zeros(ShapeND(intArrayOf(m, n))).as2D() // initialize Jacobian to Zero - val del = BroadcastDoubleTensorAlgebra.zeros(ShapeND(intArrayOf(n, 1))).as2D() + val ps = p.copyToTensor() + val J = zero(m, n) // initialize Jacobian to Zero + val del = zero(n, 1) for (j in 0 until n) { - del[j, 0] = dp[j, 0] * (1 + kotlin.math.abs(p[j, 0])) // parameter perturbation + del[j, 0] = dp[j, 0] * (1 + abs(p[j, 0])) // parameter perturbation p[j, 0] = ps[j, 0] + del[j, 0] // perturb parameter p(j) val epsilon = 0.0000001 - if (kotlin.math.abs(del[j, 0]) > epsilon) { + if (abs(del[j, 0]) > epsilon) { val y1 = evaluateFunction(func, t, p, settings.exampleNumber) settings.funcCalls += 1 if (dp[j, 0] < 0) { // backwards difference - for (i in 0 until J.shape.component1()) { - J[i, j] = (y1.as2D().minus(y).as2D())[i, 0] / del[j, 0] + for (i in 0 until J.shape.first()) { + J[i, j] = (y1 - y)[i, 0] / del[j, 0] } - } - else { + } else { // Do tests for it p[j, 0] = ps[j, 0] - del[j, 0] // central difference, additional func call - for (i in 0 until J.shape.component1()) { - J[i, j] = (y1.as2D().minus(evaluateFunction(func, t, p, settings.exampleNumber)).as2D())[i, 0] / (2 * del[j, 0]) + for (i in 0 until J.shape.first()) { + J[i, j] = (y1 - evaluateFunction(func, t, p, settings.exampleNumber))[i, 0] / (2 * del[j, 0]) } settings.funcCalls += 1 } diff --git a/kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestLmAlgorithm.kt b/kmath-tensors/src/jvmTest/kotlin/space/kscience/kmath/tensors/core/TestLmAlgorithm.kt similarity index 52% rename from kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestLmAlgorithm.kt rename to kmath-tensors/src/jvmTest/kotlin/space/kscience/kmath/tensors/core/TestLmAlgorithm.kt index 4b031eb11..4c039190a 100644 --- a/kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestLmAlgorithm.kt +++ b/kmath-tensors/src/jvmTest/kotlin/space/kscience/kmath/tensors/core/TestLmAlgorithm.kt @@ -5,77 +5,83 @@ package space.kscience.kmath.tensors.core -import space.kscience.kmath.nd.MutableStructure2D -import space.kscience.kmath.nd.ShapeND -import space.kscience.kmath.nd.as2D -import space.kscience.kmath.nd.component1 +import space.kscience.kmath.PerformancePitfall +import space.kscience.kmath.nd.* import space.kscience.kmath.operations.invoke -import space.kscience.kmath.tensors.core.DoubleTensorAlgebra.Companion.max -import space.kscience.kmath.tensors.core.DoubleTensorAlgebra.Companion.plus -import space.kscience.kmath.tensors.core.DoubleTensorAlgebra.Companion.pow -import space.kscience.kmath.tensors.core.DoubleTensorAlgebra.Companion.times -import kotlin.math.roundToLong +import space.kscience.kmath.structures.DoubleBuffer import kotlin.test.Test import kotlin.test.assertEquals +@PerformancePitfall class TestLmAlgorithm { companion object { - fun funcEasyForLm(t: MutableStructure2D, p: MutableStructure2D, exampleNumber: Int): MutableStructure2D { + fun funcEasyForLm( + t: MutableStructure2D, + p: MutableStructure2D, + exampleNumber: Int, + ): MutableStructure2D = with(DoubleTensorAlgebra) { val m = t.shape.component1() - var yHat = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf(m, 1))) + val yHat = when (exampleNumber) { + 1 -> exp((t * (-1.0 / p[1, 0]))) * p[0, 0] + (t * p[2, 0]) * exp((t * (-1.0 / p[3, 0]))) - if (exampleNumber == 1) { - yHat = DoubleTensorAlgebra.exp((t.times(-1.0 / p[1, 0]))).times(p[0, 0]) + t.times(p[2, 0]).times( - DoubleTensorAlgebra.exp((t.times(-1.0 / p[3, 0]))) - ) - } - else if (exampleNumber == 2) { - val mt = t.max() - yHat = (t.times(1.0 / mt)).times(p[0, 0]) + - (t.times(1.0 / mt)).pow(2).times(p[1, 0]) + - (t.times(1.0 / mt)).pow(3).times(p[2, 0]) + - (t.times(1.0 / mt)).pow(4).times(p[3, 0]) - } - else if (exampleNumber == 3) { - yHat = DoubleTensorAlgebra.exp((t.times(-1.0 / p[1, 0]))) - .times(p[0, 0]) + DoubleTensorAlgebra.sin((t.times(1.0 / p[3, 0]))).times(p[2, 0]) + 2 -> { + val mt = t.max() + (t * (1.0 / mt)) * p[0, 0] + + (t * (1.0 / mt)).pow(2) * p[1, 0] + + (t * (1.0 / mt)).pow(3) * p[2, 0] + + (t * (1.0 / mt)).pow(4) * p[3, 0] + } + + 3 -> exp(t * (-1.0 / p[1, 0])) * p[0, 0] + + sin((t * (1.0 / p[3, 0]))) * p[2, 0] + + else -> zeros(ShapeND(intArrayOf(m, 1))) } return yHat.as2D() } - fun funcMiddleForLm(t: MutableStructure2D, p: MutableStructure2D, exampleNumber: Int): MutableStructure2D { + fun funcMiddleForLm( + t: MutableStructure2D, + p: MutableStructure2D, + exampleNumber: Int, + ): MutableStructure2D = with(DoubleTensorAlgebra) { val m = t.shape.component1() - var yHat = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf (m, 1))) + var yHat = zeros(ShapeND(intArrayOf(m, 1))) val mt = t.max() - for(i in 0 until p.shape.component1()){ - yHat += (t.times(1.0 / mt)).times(p[i, 0]) + for (i in 0 until p.shape.component1()) { + yHat.plusAssign(t * (1.0 / mt) * p[i, 0]) } - for(i in 0 until 5){ + for (i in 0 until 5) { yHat = funcEasyForLm(yHat.as2D(), p, exampleNumber).asDoubleTensor() } return yHat.as2D() } - fun funcDifficultForLm(t: MutableStructure2D, p: MutableStructure2D, exampleNumber: Int): MutableStructure2D { + fun funcDifficultForLm( + t: MutableStructure2D, + p: MutableStructure2D, + exampleNumber: Int, + ): MutableStructure2D = with(DoubleTensorAlgebra) { val m = t.shape.component1() - var yHat = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf (m, 1))) + var yHat = zeros(ShapeND(intArrayOf(m, 1))) val mt = t.max() - for(i in 0 until p.shape.component1()){ - yHat = yHat.plus( (t.times(1.0 / mt)).times(p[i, 0]) ) + for (i in 0 until p.shape.component1()) { + yHat = yHat + (t * (1.0 / mt)) * p[i, 0] } - for(i in 0 until 4){ + for (i in 0 until 4) { yHat = funcEasyForLm((yHat.as2D() + t).as2D(), p, exampleNumber).asDoubleTensor() } return yHat.as2D() } } + @Test fun testLMEasy() = DoubleTensorAlgebra { val lmMatxYDat = doubleArrayOf( @@ -91,12 +97,12 @@ class TestLmAlgorithm { 14.7665, 13.3718, 15.0587, 13.8320, 14.7873, 13.6824, 14.2579, 14.2154, 13.5818, 13.8157 ) - var exampleNumber = 1 - val p_init = BroadcastDoubleTensorAlgebra.fromArray( + val exampleNumber = 1 + val pInit = fromArray( ShapeND(intArrayOf(4, 1)), doubleArrayOf(5.0, 2.0, 0.2, 10.0) ).as2D() - var t = ones(ShapeND(intArrayOf(100, 1))).as2D() + val t = ones(ShapeND(intArrayOf(100, 1))).as2D() for (i in 0 until 100) { t[i, 0] = t[i, 0] * (i + 1) } @@ -119,22 +125,26 @@ class TestLmAlgorithm { ShapeND(intArrayOf(4, 1)), doubleArrayOf(50.0, 20.0, 2.0, 100.0) ).as2D() - val inputData = LMInput(::funcEasyForLm, p_init, t, yDat, weight, dp, pMin, pMax, 100, - doubleArrayOf(1e-3, 1e-3, 1e-1, 1e-1), doubleArrayOf(1e-2, 11.0, 9.0), 1, 10, exampleNumber) + val inputData = LMInput( + Companion::funcEasyForLm, pInit, t, yDat, weight, dp, pMin, pMax, 100, + doubleArrayOf(1e-3, 1e-3, 1e-1, 1e-1), doubleArrayOf(1e-2, 11.0, 9.0), 1, 10, exampleNumber + ) val result = levenbergMarquardt(inputData) assertEquals(13, result.iterations) assertEquals(31, result.funcCalls) - assertEquals(0.9131368192633, (result.resultChiSq * 1e13).roundToLong() / 1e13) - assertEquals(3.7790980 * 1e-7, (result.resultLambda * 1e13).roundToLong() / 1e13) + assertEquals(0.9131368192633, result.resultChiSq, 1e-13) + assertEquals(3.7790980 * 1e-7, result.resultLambda, 1e-13) assertEquals(result.typeOfConvergence, TypeOfConvergence.InParameters) val expectedParameters = BroadcastDoubleTensorAlgebra.fromArray( ShapeND(intArrayOf(4, 1)), doubleArrayOf(20.527230909086, 9.833627103230, 0.997571256572, 50.174445822506) ).as2D() - result.resultParameters = result.resultParameters.map { x -> (x * 1e12).toLong() / 1e12}.as2D() + result.resultParameters = result.resultParameters.map { x -> (x * 1e12).toLong() / 1e12 }.as2D() val receivedParameters = BroadcastDoubleTensorAlgebra.fromArray( - ShapeND(intArrayOf(4, 1)), doubleArrayOf(result.resultParameters[0, 0], result.resultParameters[1, 0], - result.resultParameters[2, 0], result.resultParameters[3, 0]) + ShapeND(intArrayOf(4, 1)), doubleArrayOf( + result.resultParameters[0, 0], result.resultParameters[1, 0], + result.resultParameters[2, 0], result.resultParameters[3, 0] + ) ).as2D() assertEquals(expectedParameters[0, 0], receivedParameters[0, 0]) assertEquals(expectedParameters[1, 0], receivedParameters[1, 0]) @@ -143,25 +153,25 @@ class TestLmAlgorithm { } @Test - fun TestLMMiddle() = DoubleTensorAlgebra { - val NData = 100 - val tExample = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(NData, 1))).as2D() - for (i in 0 until NData) { + fun testLMMiddle() = DoubleTensorAlgebra { + val nData = 100 + val tExample = one(nData, 1).as2D() + for (i in 0 until nData) { tExample[i, 0] = tExample[i, 0] * (i + 1) } - val Nparams = 20 - val pExample = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1))).as2D() - for (i in 0 until Nparams) { + val nParams = 20 + val pExample = one(nParams, 1).as2D() + for (i in 0 until nParams) { pExample[i, 0] = pExample[i, 0] + i - 25 } val exampleNumber = 1 - val yHat = funcMiddleForLm(tExample, pExample, exampleNumber) + val yHat = funcMiddleForLm(tExample, pExample, exampleNumber) - val pInit = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf(Nparams, 1))).as2D() - for (i in 0 until Nparams) { + val pInit = zeros(ShapeND(intArrayOf(nParams, 1))).as2D() + for (i in 0 until nParams) { pInit[i, 0] = (pExample[i, 0] + 0.9) } @@ -171,13 +181,14 @@ class TestLmAlgorithm { val dp = BroadcastDoubleTensorAlgebra.fromArray( ShapeND(intArrayOf(1, 1)), DoubleArray(1) { -0.01 } ).as2D() - var pMin = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1))) - pMin = pMin.div(1.0 / -50.0) - val pMax = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1))) - pMin = pMin.div(1.0 / 50.0) + var pMin = ones(ShapeND(intArrayOf(nParams, 1))) + pMin = pMin * (-50.0) + val pMax = ones(ShapeND(intArrayOf(nParams, 1))) + pMin = pMin * 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 inputData = LMInput(::funcMiddleForLm, + val inputData = LMInput( + Companion::funcMiddleForLm, pInit.as2D(), t, yDat, @@ -190,63 +201,67 @@ class TestLmAlgorithm { doubleArrayOf(opts[6], opts[7], opts[8]), opts[9].toInt(), 10, - 1) + 1 + ) val result = DoubleTensorAlgebra.levenbergMarquardt(inputData) assertEquals(46, result.iterations) assertEquals(113, result.funcCalls) - assertEquals(0.000005977, (result.resultChiSq * 1e9).roundToLong() / 1e9) - assertEquals(1.0 * 1e-7, (result.resultLambda * 1e13).roundToLong() / 1e13) + assertEquals(0.000005977, result.resultChiSq, 1e-9) + assertEquals(1.0 * 1e-7, result.resultLambda, 1e-13) assertEquals(result.typeOfConvergence, TypeOfConvergence.InReducedChiSquare) - val expectedParameters = BroadcastDoubleTensorAlgebra.fromArray( - ShapeND(intArrayOf(Nparams, 1)), doubleArrayOf( -23.9717, -18.6686, -21.7971, + val expectedParameters = fromArray( + ShapeND(intArrayOf(nParams, 1)), doubleArrayOf( + -23.9717, -18.6686, -21.7971, -20.9681, -22.086, -20.5859, -19.0384, -17.4957, -15.9991, -14.576, -13.2441, - - 12.0201, -10.9256, -9.9878, -9.2309, -8.6589, -8.2365, -7.8783, -7.4598, -6.8511)).as2D() - result.resultParameters = result.resultParameters.map { x -> (x * 1e4).roundToLong() / 1e4}.as2D() - val receivedParameters = zeros(ShapeND(intArrayOf(Nparams, 1))).as2D() - for (i in 0 until Nparams) { + 12.0201, -10.9256, -9.9878, -9.2309, -8.6589, -8.2365, -7.8783, -7.4598, -6.8511 + ) + ) + val receivedParameters = zero(nParams, 1) + for (i in 0 until nParams) { receivedParameters[i, 0] = result.resultParameters[i, 0] - assertEquals(expectedParameters[i, 0], result.resultParameters[i, 0]) + assertEquals(expectedParameters[i, 0], result.resultParameters[i, 0], 1e-2) } } @Test fun TestLMDifficult() = DoubleTensorAlgebra { - val NData = 200 - var tExample = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(NData, 1))).as2D() - for (i in 0 until NData) { + val nData = 200 + val tExample = ones(ShapeND(intArrayOf(nData, 1))).as2D() + for (i in 0 until nData) { tExample[i, 0] = tExample[i, 0] * (i + 1) - 104 } - val Nparams = 15 - var pExample = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1))).as2D() - for (i in 0 until Nparams) { + val nParams = 15 + val pExample = ones(ShapeND(intArrayOf(nParams, 1))).as2D() + for (i in 0 until nParams) { pExample[i, 0] = pExample[i, 0] + i - 25 } val exampleNumber = 1 - var yHat = funcDifficultForLm(tExample, pExample, exampleNumber) + val yHat = funcDifficultForLm(tExample, pExample, exampleNumber) - var pInit = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf(Nparams, 1))).as2D() - for (i in 0 until Nparams) { + val pInit = zeros(ShapeND(intArrayOf(nParams, 1))).as2D() + for (i in 0 until nParams) { pInit[i, 0] = (pExample[i, 0] + 0.9) } - var t = tExample + val t = tExample val yDat = yHat - val weight = 1.0 / Nparams * 1.0 - 0.085 - val dp = BroadcastDoubleTensorAlgebra.fromArray( + val weight = 1.0 / nParams * 1.0 - 0.085 + val dp = fromArray( ShapeND(intArrayOf(1, 1)), DoubleArray(1) { -0.01 } ).as2D() - var pMin = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1))) - pMin = pMin.div(1.0 / -50.0) - val pMax = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1))) - pMin = pMin.div(1.0 / 50.0) + var pMin = ones(ShapeND(intArrayOf(nParams, 1))) + pMin = pMin * (-50.0) + val pMax = ones(ShapeND(intArrayOf(nParams, 1))) + pMin = pMin * (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 inputData = LMInput(::funcDifficultForLm, + val inputData = LMInput( + Companion::funcDifficultForLm, pInit.as2D(), t, yDat, @@ -259,22 +274,37 @@ class TestLmAlgorithm { doubleArrayOf(opts[6], opts[7], opts[8]), opts[9].toInt(), 10, - 1) + 1 + ) val result = DoubleTensorAlgebra.levenbergMarquardt(inputData) assertEquals(2375, result.iterations) assertEquals(4858, result.funcCalls) - assertEquals(5.14347, (result.resultLambda * 1e5).roundToLong() / 1e5) + assertEquals(5.14347, result.resultLambda, 1e-5) assertEquals(result.typeOfConvergence, TypeOfConvergence.InParameters) - val expectedParameters = BroadcastDoubleTensorAlgebra.fromArray( - ShapeND(intArrayOf(Nparams, 1)), doubleArrayOf(-23.6412, -16.7402, -21.5705, -21.0464, - -17.2852, -17.2959, -17.298, 0.9999, -17.2885, -17.3008, -17.2941, -17.2923, -17.2976, -17.3028, -17.2891)).as2D() - result.resultParameters = result.resultParameters.map { x -> (x * 1e4).roundToLong() / 1e4}.as2D() - val receivedParameters = zeros(ShapeND(intArrayOf(Nparams, 1))).as2D() - for (i in 0 until Nparams) { + val expectedParameters = DoubleBuffer( + -23.6412, + -16.7402, + -21.5705, + -21.0464, + -17.2852, + -17.2959, + -17.298, + 0.9999, + -17.2885, + -17.3008, + -17.2941, + -17.2923, + -17.2976, + -17.3028, + -17.2891 + ) + + val receivedParameters = zeros(ShapeND(intArrayOf(nParams, 1))).as2D() + for (i in 0 until nParams) { receivedParameters[i, 0] = result.resultParameters[i, 0] - assertEquals(expectedParameters[i, 0], result.resultParameters[i, 0]) + assertEquals(expectedParameters[i], result.resultParameters[i, 0], 1e-2) } } } \ No newline at end of file