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 e996b7f7e..24cfa955a 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,6 @@ 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.LMSettings import space.kscience.kmath.tensors.core.lm import kotlin.math.roundToInt @@ -29,9 +28,9 @@ fun main() { p_example[i, 0] = p_example[i, 0] + i - 25 } - val settings = LMSettings(0, 0, 1) + val exampleNumber = 1 - var y_hat = funcDifficultForLm(t_example, p_example, settings) + var y_hat = funcDifficultForLm(t_example, p_example, exampleNumber) var p_init = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf(Nparams, 1))).as2D() for (i in 0 until Nparams) { @@ -72,14 +71,14 @@ fun main() { ) println("Parameters:") - for (i in 0 until result.result_parameters.shape.component1()) { - val x = (result.result_parameters[i, 0] * 10000).roundToInt() / 10000.0 + for (i in 0 until result.resultParameters.shape.component1()) { + val x = (result.resultParameters[i, 0] * 10000).roundToInt() / 10000.0 print("$x ") } println() println("Y true and y received:") - var y_hat_after = funcDifficultForLm(t_example, result.result_parameters, settings) + var y_hat_after = funcDifficultForLm(t_example, result.resultParameters, exampleNumber) for (i in 0 until y_hat.shape.component1()) { val x = (y_hat[i, 0] * 10000).roundToInt() / 10000.0 val y = (y_hat_after[i, 0] * 10000).roundToInt() / 10000.0 @@ -87,7 +86,7 @@ fun main() { } println("Сhi_sq:") - println(result.result_chi_sq) + println(result.resultChiSq) println("Number of iterations:") println(result.iterations) } \ No newline at end of file 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 d9e5f350e..c7b2b0def 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,7 +12,6 @@ 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.LMSettings import space.kscience.kmath.tensors.core.lm import kotlin.math.roundToInt @@ -35,14 +34,14 @@ fun main() { ) println("Parameters:") - for (i in 0 until result.result_parameters.shape.component1()) { - val x = (result.result_parameters[i, 0] * 10000).roundToInt() / 10000.0 + for (i in 0 until result.resultParameters.shape.component1()) { + val x = (result.resultParameters[i, 0] * 10000).roundToInt() / 10000.0 print("$x ") } println() println("Y true and y received:") - var y_hat_after = funcDifficultForLm(startedData.t, result.result_parameters, LMSettings(0, 0, startedData.example_number)) + var y_hat_after = funcDifficultForLm(startedData.t, result.resultParameters, startedData.example_number) for (i in 0 until startedData.y_dat.shape.component1()) { val x = (startedData.y_dat[i, 0] * 10000).roundToInt() / 10000.0 val y = (y_hat_after[i, 0] * 10000).roundToInt() / 10000.0 @@ -50,7 +49,7 @@ fun main() { } println("Сhi_sq:") - println(result.result_chi_sq) + println(result.resultChiSq) println("Number of iterations:") println(result.iterations) } \ No newline at end of file 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 b5553a930..471143102 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,8 +12,6 @@ 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.DoubleTensorAlgebra.Companion.times -import space.kscience.kmath.tensors.core.LMSettings import space.kscience.kmath.tensors.core.lm import kotlin.math.roundToInt fun main() { @@ -29,16 +27,15 @@ fun main() { p_example[i, 0] = p_example[i, 0] + i - 25 } - val settings = LMSettings(0, 0, 1) + val exampleNumber = 1 - var y_hat = funcMiddleForLm(t_example, p_example, settings) + var y_hat = funcMiddleForLm(t_example, p_example, exampleNumber) var p_init = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf(Nparams, 1))).as2D() for (i in 0 until Nparams) { p_init[i, 0] = (p_example[i, 0] + 0.9) } -// val p_init = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1))) -// val p_init = p_example + var t = t_example val y_dat = y_hat val weight = BroadcastDoubleTensorAlgebra.fromArray( @@ -54,7 +51,7 @@ fun main() { val consts = BroadcastDoubleTensorAlgebra.fromArray( ShapeND(intArrayOf(1, 1)), doubleArrayOf(0.0) ).as2D() - val opts = doubleArrayOf(3.0, 10000.0, 1e-3, 1e-3, 1e-3, 1e-3, 1e-15, 11.0, 9.0, 1.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, @@ -72,14 +69,14 @@ fun main() { ) println("Parameters:") - for (i in 0 until result.result_parameters.shape.component1()) { - val x = (result.result_parameters[i, 0] * 10000).roundToInt() / 10000.0 + for (i in 0 until result.resultParameters.shape.component1()) { + val x = (result.resultParameters[i, 0] * 10000).roundToInt() / 10000.0 print("$x ") } println() - var y_hat_after = funcMiddleForLm(t_example, result.result_parameters, settings) + var y_hat_after = funcMiddleForLm(t_example, result.resultParameters, exampleNumber) for (i in 0 until y_hat.shape.component1()) { val x = (y_hat[i, 0] * 10000).roundToInt() / 10000.0 val y = (y_hat_after[i, 0] * 10000).roundToInt() / 10000.0 @@ -87,7 +84,7 @@ fun main() { } println("Сhi_sq:") - println(result.result_chi_sq) + println(result.resultChiSq) println("Number of iterations:") println(result.iterations) } \ No newline at end of file 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 46bda40c6..f031d82bf 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,12 +11,11 @@ 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.LMSettings import space.kscience.kmath.tensors.core.lm import kotlin.random.Random import kotlin.reflect.KFunction3 -fun streamLm(lm_func: KFunction3, MutableStructure2D, LMSettings, MutableStructure2D>, +fun streamLm(lm_func: KFunction3, MutableStructure2D, Int, MutableStructure2D>, startData: StartDataLm, launchFrequencyInMs: Long, numberOfLaunches: Int): Flow> = flow{ var example_number = startData.example_number @@ -48,9 +47,9 @@ fun streamLm(lm_func: KFunction3, MutableStructure2D< 10, example_number ) - emit(result.result_parameters) + emit(result.resultParameters) delay(launchFrequencyInMs) - p_init = result.result_parameters + p_init = result.resultParameters y_dat = generateNewYDat(y_dat, 0.1) if (!isEndless) steps -= 1 } 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 a22de71d8..537b86da3 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 @@ -16,7 +16,6 @@ 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 space.kscience.kmath.tensors.core.LMSettings import space.kscience.kmath.tensors.core.asDoubleTensor public data class StartDataLm ( @@ -33,23 +32,31 @@ public data class StartDataLm ( var opts: DoubleArray ) -fun funcDifficultForLm(t: MutableStructure2D, p: MutableStructure2D, settings: LMSettings): MutableStructure2D { +fun funcEasyForLm(t: MutableStructure2D, p: MutableStructure2D, exampleNumber: Int): MutableStructure2D { val m = t.shape.component1() - var y_hat = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf (m, 1))) + var y_hat = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf(m, 1))) - val mt = t.max() - for(i in 0 until p.shape.component1()){ - y_hat = y_hat.plus( (t.times(1.0 / mt)).times(p[i, 0]) ) + if (exampleNumber == 1) { + y_hat = 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]))) + ) } - - for(i in 0 until 4){ - y_hat = funcEasyForLm((y_hat.as2D() + t).as2D(), p, settings).asDoubleTensor() + else if (exampleNumber == 2) { + val mt = t.max() + y_hat = (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) { + y_hat = 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]) } return y_hat.as2D() } -fun funcMiddleForLm(t: MutableStructure2D, p: MutableStructure2D, settings: LMSettings): MutableStructure2D { +fun funcMiddleForLm(t: MutableStructure2D, p: MutableStructure2D, exampleNumber: Int): MutableStructure2D { val m = t.shape.component1() var y_hat = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf (m, 1))) @@ -59,36 +66,29 @@ fun funcMiddleForLm(t: MutableStructure2D, p: MutableStructure2D } for(i in 0 until 5){ - y_hat = funcEasyForLm(y_hat.as2D(), p, settings).asDoubleTensor() + y_hat = funcEasyForLm(y_hat.as2D(), p, exampleNumber).asDoubleTensor() } return y_hat.as2D() } -fun funcEasyForLm(t: MutableStructure2D, p: MutableStructure2D, settings: LMSettings): MutableStructure2D { +fun funcDifficultForLm(t: MutableStructure2D, p: MutableStructure2D, exampleNumber: Int): MutableStructure2D { val m = t.shape.component1() var y_hat = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf (m, 1))) - if (settings.example_number == 1) { - y_hat = 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]))) - ) + val mt = t.max() + for(i in 0 until p.shape.component1()){ + y_hat = y_hat.plus( (t.times(1.0 / mt)).times(p[i, 0]) ) } - else if (settings.example_number == 2) { - val mt = t.max() - y_hat = (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 (settings.example_number == 3) { - y_hat = 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]) + + for(i in 0 until 4){ + y_hat = funcEasyForLm((y_hat.as2D() + t).as2D(), p, exampleNumber).asDoubleTensor() } return y_hat.as2D() } + fun getStartDataForFuncDifficult(): StartDataLm { val NData = 200 var t_example = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(NData, 1))).as2D() @@ -102,9 +102,9 @@ fun getStartDataForFuncDifficult(): StartDataLm { p_example[i, 0] = p_example[i, 0] + i - 25 } - val settings = LMSettings(0, 0, 1) + val exampleNumber = 1 - var y_hat = funcDifficultForLm(t_example, p_example, settings) + var y_hat = funcDifficultForLm(t_example, p_example, exampleNumber) var p_init = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf(Nparams, 1))).as2D() for (i in 0 until Nparams) { @@ -144,9 +144,9 @@ fun getStartDataForFuncMiddle(): StartDataLm { p_example[i, 0] = p_example[i, 0] + i - 25 } - val settings = LMSettings(0, 0, 1) + val exampleNumber = 1 - var y_hat = funcMiddleForLm(t_example, p_example, settings) + var y_hat = funcMiddleForLm(t_example, p_example, exampleNumber) var p_init = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf(Nparams, 1))).as2D() for (i in 0 until Nparams) { 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 deb7ee300..7fbf6ecd8 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 @@ -40,36 +40,39 @@ public enum class TypeOfConvergence{ NoConvergence } +/** + * Class for 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 + */ public data class LMResultInfo ( var iterations:Int, - var func_calls: Int, - var example_number: Int, - var result_chi_sq: Double, - var result_lambda: Double, - var result_parameters: MutableStructure2D, + var funcCalls: Int, + var resultChiSq: Double, + var resultLambda: Double, + var resultParameters: MutableStructure2D, var typeOfConvergence: TypeOfConvergence, - var epsilon: Double -) - -public data class LMSettings ( - var iteration:Int, - var func_calls: Int, - var example_number:Int ) public fun DoubleTensorAlgebra.lm( - func: KFunction3, MutableStructure2D, LMSettings, MutableStructure2D>, + func: KFunction3, MutableStructure2D, Int, MutableStructure2D>, p_input: MutableStructure2D, t_input: MutableStructure2D, y_dat_input: MutableStructure2D, weight_input: MutableStructure2D, dp_input: MutableStructure2D, p_min_input: MutableStructure2D, p_max_input: MutableStructure2D, c_input: MutableStructure2D, opts_input: DoubleArray, nargin: Int, example_number: Int): LMResultInfo { - val resultInfo = LMResultInfo(0, 0, example_number, 0.0, - 0.0, p_input, TypeOfConvergence.NoConvergence, 0.0) + + val resultInfo = LMResultInfo(0, 0, 0.0, + 0.0, p_input, TypeOfConvergence.NoConvergence) val eps:Double = 2.2204e-16 val settings = LMSettings(0, 0, example_number) - settings.func_calls = 0 // running count of function evaluations + settings.funcCalls = 0 // running count of function evaluations var p = p_input val y_dat = y_dat_input @@ -160,7 +163,7 @@ public fun DoubleTensorAlgebra.lm( val idx = get_zero_indices(dp) // indices of the parameters to be fit val Nfit = idx?.shape?.component1() // number of parameters to fit var stop = false // termination flag - val y_init = feval(func, t, p, settings) // residual error using p_try + val y_init = feval(func, t, p, example_number) // residual error using p_try 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() @@ -219,7 +222,7 @@ public fun DoubleTensorAlgebra.lm( var p_try = (p + h).as2D() // update the [idx] elements p_try = smallest_element_comparison(largest_element_comparison(p_min, p_try.as2D()), p_max) // apply constraints - var delta_y = y_dat.minus(feval(func, t, p_try, settings)) // residual error using p_try + var delta_y = y_dat.minus(feval(func, t, p_try, example_number)) // residual error using p_try for (i in 0 until delta_y.shape.component1()) { // floating point error; break for (j in 0 until delta_y.shape.component2()) { @@ -230,7 +233,7 @@ public fun DoubleTensorAlgebra.lm( } } - settings.func_calls += 1 + settings.funcCalls += 1 val tmp = delta_y.times(weight) var X2_try = delta_y.as2D().transpose().dot(tmp) // Chi-squared error criteria @@ -244,8 +247,8 @@ public fun DoubleTensorAlgebra.lm( p_try = p.plus(h).as2D() // update only [idx] elements p_try = smallest_element_comparison(largest_element_comparison(p_min, p_try), p_max) // apply constraints - var delta_y = y_dat.minus(feval(func, t, p_try, settings)) // residual error using p_try - settings.func_calls += 1 + var delta_y = y_dat.minus(feval(func, t, p_try, example_number)) // residual error using p_try + settings.funcCalls += 1 val tmp = delta_y.times(weight) X2_try = delta_y.as2D().transpose().dot(tmp) // Chi-squared error criteria @@ -320,10 +323,10 @@ public fun DoubleTensorAlgebra.lm( if (prnt > 1) { val chi_sq = X2 / DoF resultInfo.iterations = settings.iteration - resultInfo.func_calls = settings.func_calls - resultInfo.result_chi_sq = chi_sq - resultInfo.result_lambda = lambda - resultInfo.result_parameters = p + resultInfo.funcCalls = settings.funcCalls + resultInfo.resultChiSq = chi_sq + resultInfo.resultLambda = lambda + resultInfo.resultParameters = p } // update convergence history ... save _reduced_ Chi-square @@ -331,30 +334,32 @@ public fun DoubleTensorAlgebra.lm( if (abs(JtWdy).max()!! < epsilon_1 && settings.iteration > 2) { resultInfo.typeOfConvergence = TypeOfConvergence.InGradient - resultInfo.epsilon = epsilon_1 stop = true } if ((abs(h.as2D()).div(abs(p) + 1e-12)).max() < epsilon_2 && settings.iteration > 2) { resultInfo.typeOfConvergence = TypeOfConvergence.InParameters - resultInfo.epsilon = epsilon_2 stop = true } if (X2 / DoF < epsilon_3 && settings.iteration > 2) { resultInfo.typeOfConvergence = TypeOfConvergence.InReducedChiSquare - resultInfo.epsilon = epsilon_3 stop = true } if (settings.iteration == MaxIter) { resultInfo.typeOfConvergence = TypeOfConvergence.NoConvergence - resultInfo.epsilon = 0.0 stop = true } } return resultInfo } +private data class LMSettings ( + var iteration:Int, + var funcCalls: Int, + var exampleNumber:Int +) + /* matrix -> column of all elemnets */ -public fun make_column(tensor: MutableStructure2D) : MutableStructure2D { +private fun make_column(tensor: MutableStructure2D) : MutableStructure2D { val shape = intArrayOf(tensor.shape.component1() * tensor.shape.component2(), 1) val buffer = DoubleArray(tensor.shape.component1() * tensor.shape.component2()) for (i in 0 until tensor.shape.component1()) { @@ -367,11 +372,11 @@ public fun make_column(tensor: MutableStructure2D) : MutableStructure2D< } /* column length */ -public fun length(column: MutableStructure2D) : Int { +private fun length(column: MutableStructure2D) : Int { return column.shape.component1() } -public fun MutableStructure2D.abs() { +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]) @@ -379,7 +384,7 @@ public fun MutableStructure2D.abs() { } } -public fun abs(input: MutableStructure2D): MutableStructure2D { +private fun abs(input: MutableStructure2D): MutableStructure2D { val tensor = BroadcastDoubleTensorAlgebra.ones( ShapeND( intArrayOf( @@ -396,7 +401,7 @@ public fun abs(input: MutableStructure2D): MutableStructure2D { return tensor } -public fun diag(input: MutableStructure2D): MutableStructure2D { +private fun diag(input: MutableStructure2D): MutableStructure2D { val tensor = BroadcastDoubleTensorAlgebra.ones(ShapeND(intArrayOf(input.shape.component1(), 1))).as2D() for (i in 0 until tensor.shape.component1()) { tensor[i, 0] = input[i, i] @@ -404,7 +409,7 @@ public fun diag(input: MutableStructure2D): MutableStructure2D { return tensor } -public fun make_matrx_with_diagonal(column: MutableStructure2D): MutableStructure2D { +private fun make_matrx_with_diagonal(column: MutableStructure2D): MutableStructure2D { val size = column.shape.component1() val tensor = BroadcastDoubleTensorAlgebra.zeros(ShapeND(intArrayOf(size, size))).as2D() for (i in 0 until size) { @@ -413,12 +418,12 @@ public fun make_matrx_with_diagonal(column: MutableStructure2D): Mutable return tensor } -public fun lm_eye(size: Int): MutableStructure2D { +private fun lm_eye(size: Int): MutableStructure2D { val column = BroadcastDoubleTensorAlgebra.ones(ShapeND(intArrayOf(size, 1))).as2D() return make_matrx_with_diagonal(column) } -public fun largest_element_comparison(a: MutableStructure2D, b: MutableStructure2D): MutableStructure2D { +private fun largest_element_comparison(a: MutableStructure2D, b: MutableStructure2D): MutableStructure2D { val a_sizeX = a.shape.component1() val a_sizeY = a.shape.component2() val b_sizeX = b.shape.component1() @@ -440,7 +445,7 @@ public fun largest_element_comparison(a: MutableStructure2D, b: MutableS return tensor } -public fun smallest_element_comparison(a: MutableStructure2D, b: MutableStructure2D): MutableStructure2D { +private fun smallest_element_comparison(a: MutableStructure2D, b: MutableStructure2D): MutableStructure2D { val a_sizeX = a.shape.component1() val a_sizeY = a.shape.component2() val b_sizeX = b.shape.component1() @@ -462,7 +467,7 @@ public fun smallest_element_comparison(a: MutableStructure2D, b: Mutable return tensor } -public fun get_zero_indices(column: MutableStructure2D, epsilon: Double = 0.000001): MutableStructure2D? { +private fun get_zero_indices(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) { @@ -475,14 +480,14 @@ public fun get_zero_indices(column: MutableStructure2D, epsilon: Double return null } -public fun feval(func: (MutableStructure2D, MutableStructure2D, LMSettings) -> MutableStructure2D, - t: MutableStructure2D, p: MutableStructure2D, settings: LMSettings) +private fun feval(func: (MutableStructure2D, MutableStructure2D, Int) -> MutableStructure2D, + t: MutableStructure2D, p: MutableStructure2D, exampleNumber: Int) : MutableStructure2D { - return func(t, p, settings) + return func(t, p, exampleNumber) } -public fun lm_matx(func: (MutableStructure2D, MutableStructure2D, LMSettings) -> MutableStructure2D, +private fun lm_matx(func: (MutableStructure2D, MutableStructure2D, Int) -> MutableStructure2D, t: MutableStructure2D, p_old: MutableStructure2D, y_old: MutableStructure2D, dX2: Int, J_input: MutableStructure2D, p: MutableStructure2D, y_dat: MutableStructure2D, weight: MutableStructure2D, dp:MutableStructure2D, settings:LMSettings) : Array> @@ -492,8 +497,8 @@ public fun lm_matx(func: (MutableStructure2D, MutableStructure2D val Npnt = length(y_dat) // number of data points val Npar = length(p) // number of parameters - val y_hat = feval(func, t, p, settings) // evaluate model using parameters 'p' - settings.func_calls += 1 + val y_hat = feval(func, t, p, settings.exampleNumber) // evaluate model using parameters 'p' + settings.funcCalls += 1 var J = J_input @@ -513,7 +518,7 @@ public fun lm_matx(func: (MutableStructure2D, MutableStructure2D return arrayOf(JtWJ,JtWdy,Chi_sq,y_hat,J) } -public fun lm_Broyden_J(p_old: MutableStructure2D, y_old: MutableStructure2D, J_input: MutableStructure2D, +private fun lm_Broyden_J(p_old: MutableStructure2D, y_old: MutableStructure2D, J_input: MutableStructure2D, p: MutableStructure2D, y: MutableStructure2D): MutableStructure2D { var J = J_input.copyToTensor() @@ -524,7 +529,7 @@ public fun lm_Broyden_J(p_old: MutableStructure2D, y_old: MutableStructu return J.as2D() } -public fun lm_FD_J(func: (MutableStructure2D, MutableStructure2D, settings: LMSettings) -> MutableStructure2D, +private fun lm_FD_J(func: (MutableStructure2D, MutableStructure2D, exampleNumber: Int) -> MutableStructure2D, t: MutableStructure2D, p: MutableStructure2D, y: MutableStructure2D, dp: MutableStructure2D, settings: LMSettings): MutableStructure2D { // default: dp = 0.001 * ones(1,n) @@ -543,8 +548,8 @@ public fun lm_FD_J(func: (MutableStructure2D, MutableStructure2D val epsilon = 0.0000001 if (kotlin.math.abs(del[j, 0]) > epsilon) { - val y1 = feval(func, t, p, settings) - settings.func_calls += 1 + val y1 = feval(func, t, p, settings.exampleNumber) + settings.funcCalls += 1 if (dp[j, 0] < 0) { // backwards difference for (i in 0 until J.shape.component1()) { @@ -556,9 +561,9 @@ public fun lm_FD_J(func: (MutableStructure2D, MutableStructure2D println("Potential mistake") 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(feval(func, t, p, settings)).as2D())[i, 0] / (2 * del[j, 0]) + J[i, j] = (y1.as2D().minus(feval(func, t, p, settings.exampleNumber)).as2D())[i, 0] / (2 * del[j, 0]) } - settings.func_calls += 1 + settings.funcCalls += 1 } } 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 1112043fc..c3e5fa16f 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 @@ -20,23 +20,23 @@ import kotlin.test.assertEquals class TestLmAlgorithm { companion object { - fun funcEasyForLm(t: MutableStructure2D, p: MutableStructure2D, settings: LMSettings): MutableStructure2D { + fun funcEasyForLm(t: MutableStructure2D, p: MutableStructure2D, exampleNumber: Int): MutableStructure2D { val m = t.shape.component1() var y_hat = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf(m, 1))) - if (settings.example_number == 1) { + if (exampleNumber == 1) { y_hat = 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 (settings.example_number == 2) { + else if (exampleNumber == 2) { val mt = t.max() y_hat = (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 (settings.example_number == 3) { + else if (exampleNumber == 3) { y_hat = 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]) } @@ -44,7 +44,7 @@ class TestLmAlgorithm { return y_hat.as2D() } - fun funcMiddleForLm(t: MutableStructure2D, p: MutableStructure2D, settings: LMSettings): MutableStructure2D { + fun funcMiddleForLm(t: MutableStructure2D, p: MutableStructure2D, exampleNumber: Int): MutableStructure2D { val m = t.shape.component1() var y_hat = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf (m, 1))) @@ -54,13 +54,13 @@ class TestLmAlgorithm { } for(i in 0 until 5){ - y_hat = funcEasyForLm(y_hat.as2D(), p, settings).asDoubleTensor() + y_hat = funcEasyForLm(y_hat.as2D(), p, exampleNumber).asDoubleTensor() } return y_hat.as2D() } - fun funcDifficultForLm(t: MutableStructure2D, p: MutableStructure2D, settings: LMSettings): MutableStructure2D { + fun funcDifficultForLm(t: MutableStructure2D, p: MutableStructure2D, exampleNumber: Int): MutableStructure2D { val m = t.shape.component1() var y_hat = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf (m, 1))) @@ -70,12 +70,11 @@ class TestLmAlgorithm { } for(i in 0 until 4){ - y_hat = funcEasyForLm((y_hat.as2D() + t).as2D(), p, settings).asDoubleTensor() + y_hat = funcEasyForLm((y_hat.as2D() + t).as2D(), p, exampleNumber).asDoubleTensor() } return y_hat.as2D() } - } @Test fun testLMEasy() = DoubleTensorAlgebra { @@ -130,18 +129,17 @@ class TestLmAlgorithm { val result = lm(::funcEasyForLm, p_init, t, y_dat, weight, dp, p_min, p_max, consts, opts, 10, example_number) assertEquals(13, result.iterations) - assertEquals(31, result.func_calls) - assertEquals(1, result.example_number) - assertEquals(0.9131368192633, (result.result_chi_sq * 1e13).roundToLong() / 1e13) - assertEquals(3.7790980 * 1e-7, (result.result_lambda * 1e13).roundToLong() / 1e13) + assertEquals(31, result.funcCalls) + assertEquals(0.9131368192633, (result.resultChiSq * 1e13).roundToLong() / 1e13) + assertEquals(3.7790980 * 1e-7, (result.resultLambda * 1e13).roundToLong() / 1e13) assertEquals(result.typeOfConvergence, TypeOfConvergence.InParameters) val expectedParameters = BroadcastDoubleTensorAlgebra.fromArray( ShapeND(intArrayOf(4, 1)), doubleArrayOf(20.527230909086, 9.833627103230, 0.997571256572, 50.174445822506) ).as2D() - result.result_parameters = result.result_parameters.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.result_parameters[0, 0], result.result_parameters[1, 0], - result.result_parameters[2, 0], result.result_parameters[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]) @@ -163,9 +161,9 @@ class TestLmAlgorithm { p_example[i, 0] = p_example[i, 0] + i - 25 } - val settings = LMSettings(0, 0, 1) + val exampleNumber = 1 - var y_hat = funcMiddleForLm(t_example, p_example, settings) + var y_hat = funcMiddleForLm(t_example, p_example, exampleNumber) var p_init = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf(Nparams, 1))).as2D() for (i in 0 until Nparams) { @@ -219,9 +217,9 @@ class TestLmAlgorithm { p_example[i, 0] = p_example[i, 0] + i - 25 } - val settings = LMSettings(0, 0, 1) + val exampleNumber = 1 - var y_hat = funcDifficultForLm(t_example, p_example, settings) + var y_hat = funcDifficultForLm(t_example, p_example, exampleNumber) var p_init = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf(Nparams, 1))).as2D() for (i in 0 until Nparams) {