From 8d81d2d8d5a39d0309dbd2c4b70b516092f5bfea Mon Sep 17 00:00:00 2001 From: Margarita Lashina Date: Tue, 6 Jun 2023 01:41:08 +0300 Subject: [PATCH] move lm-algorithm from DoubleTensorAlgebra as extension --- .../tensors/api/LinearOpsTensorAlgebra.kt | 14 - .../kmath/tensors/core/DoubleTensorAlgebra.kt | 363 ------------ .../core/LevenbergMarquardtAlgorithm.kt | 553 ++++++++++++++++++ .../kmath/tensors/core/internal/linUtils.kt | 229 +------- .../tensors/core/TestDoubleTensorAlgebra.kt | 100 ---- .../kmath/tensors/core/TestLmAlgorithm.kt | 1 - 6 files changed, 554 insertions(+), 706 deletions(-) create mode 100644 kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/LevenbergMarquardtAlgorithm.kt diff --git a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/api/LinearOpsTensorAlgebra.kt b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/api/LinearOpsTensorAlgebra.kt index 7964656f1..6b3859316 100644 --- a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/api/LinearOpsTensorAlgebra.kt +++ b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/api/LinearOpsTensorAlgebra.kt @@ -7,15 +7,7 @@ package space.kscience.kmath.tensors.api import space.kscience.kmath.nd.MutableStructure2D import space.kscience.kmath.nd.StructureND -import space.kscience.kmath.nd.as2D import space.kscience.kmath.operations.Field -import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra -import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra.dot -import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra.map -import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra.transposed -import space.kscience.kmath.tensors.core.DoubleTensorAlgebra -import space.kscience.kmath.tensors.core.internal.LMSettings -import kotlin.reflect.KFunction3 /** * Common linear algebra operations. Operates on [Tensor]. @@ -137,10 +129,4 @@ public interface LinearOpsTensorAlgebra> : TensorPartialDivision var typeOfConvergence: TypeOfConvergence, var epsilon: Double ) - - public fun lm( - func: KFunction3, MutableStructure2D, LMSettings, 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 } diff --git a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleTensorAlgebra.kt b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleTensorAlgebra.kt index 2b8bf84c0..c8cf56888 100644 --- a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleTensorAlgebra.kt +++ b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleTensorAlgebra.kt @@ -9,7 +9,6 @@ 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.operations.DoubleBufferOps import space.kscience.kmath.operations.DoubleField @@ -19,7 +18,6 @@ import space.kscience.kmath.tensors.api.LinearOpsTensorAlgebra import space.kscience.kmath.tensors.api.Tensor import space.kscience.kmath.tensors.core.internal.* import kotlin.math.* -import kotlin.reflect.KFunction3 /** * Implementation of basic operations over double tensors and basic algebra operations on them. @@ -716,367 +714,6 @@ public open class DoubleTensorAlgebra : val aInverse = aSvd.third.dot(s).dot(aSvd.first.transposed()) return aInverse.dot(b).as2D() } - - override fun lm( - func: KFunction3, MutableStructure2D, LMSettings, 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): LinearOpsTensorAlgebra.LMResultInfo { - - var resultInfo = LinearOpsTensorAlgebra.LMResultInfo(0, 0, example_number, 0.0, - 0.0, p_input, LinearOpsTensorAlgebra.TypeOfConvergence.noConvergence, 0.0) - - val eps:Double = 2.2204e-16 - - var settings = LMSettings(0, 0, example_number) - settings.func_calls = 0 // running count of function evaluations - - var p = p_input - val y_dat = y_dat_input - val t = t_input - - val Npar = length(p) // number of parameters - val Npnt = length(y_dat) // number of data points - var p_old = zeros(ShapeND(intArrayOf(Npar, 1))).as2D() // previous set of parameters - var y_old = 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 X2_old = 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 - - var corr_p = 0 - var sigma_p = 0 - var sigma_y = 0 - var R_sq = 0 - var cvg_hist = 0 - - if (length(t) != length(y_dat)) { - // println("lm.m error: the length of t must equal the length of y_dat") - val length_t = length(t) - val length_y_dat = length(y_dat) - X2 = 0.0 - - corr_p = 0 - sigma_p = 0 - sigma_y = 0 - R_sq = 0 - cvg_hist = 0 - } - - var weight = weight_input - if (nargin < 5) { - fromArray(ShapeND(intArrayOf(1, 1)), doubleArrayOf(1.0)).as2D() - } - - var dp = dp_input - if (nargin < 6) { - dp = fromArray(ShapeND(intArrayOf(1, 1)), doubleArrayOf(-0.001)).as2D() - } - - var p_min = p_min_input - if (nargin < 7) { - p_min = p - p_min.abs() - p_min = p_min.div(-100.0).as2D() - } - - var p_max = p_max_input - if (nargin < 8) { - p_max = p - p_max.abs() - p_max = p_max.div(100.0).as2D() - } - - var c = c_input - if (nargin < 9) { - c = fromArray(ShapeND(intArrayOf(1, 1)), doubleArrayOf(1.0)).as2D() - } - - var opts = opts_input - 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) - } - - val prnt = opts[0] // >1 intermediate results; >2 plots - val MaxIter = opts[1].toInt() // maximum number of iterations - val epsilon_1 = opts[2] // convergence tolerance for gradient - val epsilon_2 = opts[3] // convergence tolerance for parameters - val epsilon_3 = opts[4] // convergence tolerance for Chi-square - val epsilon_4 = opts[5] // determines acceptance of a L-M step - val lambda_0 = opts[6] // initial value of damping paramter, lambda - val lambda_UP_fac = opts[7] // factor for increasing lambda - val lambda_DN_fac = opts[8] // factor for decreasing lambda - val Update_Type = opts[9].toInt() // 1: Levenberg-Marquardt lambda update - // 2: Quadratic update - // 3: Nielsen's lambda update equations - - p_min = make_column(p_min) - p_max = make_column(p_max) - - if (length(make_column(dp)) == 1) { - dp = ones(ShapeND(intArrayOf(Npar, 1))).div(1 / dp[0, 0]).as2D() - } - - 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 - - if (weight.shape.component1() == 1 || variance(weight) == 0.0) { // identical weights vector - weight = ones(ShapeND(intArrayOf(Npnt, 1))).div(1 / abs(weight[0, 0])).as2D() - // println("using uniform weights for error analysis") - } - else { - weight = make_column(weight) - weight.abs() - } - - // initialize Jacobian with finite difference calculation - var lm_matx_ans = lm_matx(func, t, p_old, y_old,1, J, p, y_dat, weight, dp, settings) - var JtWJ = lm_matx_ans[0] - var JtWdy = lm_matx_ans[1] - X2 = lm_matx_ans[2][0, 0] - var y_hat = lm_matx_ans[3] - J = lm_matx_ans[4] - - if ( abs(JtWdy).max()!! < epsilon_1 ) { -// println(" *** Your Initial Guess is Extremely Close to Optimal ***\n") -// println(" *** epsilon_1 = %e\n$epsilon_1") - stop = true - } - - var lambda = 1.0 - var nu = 1 - when (Update_Type) { - 1 -> lambda = lambda_0 // Marquardt: init'l lambda - else -> { // Quadratic and Nielsen - lambda = lambda_0 * (diag(JtWJ)).max()!! - nu = 2 - } - } - - X2_old = X2 // previous value of X2 - var cvg_hst = ones(ShapeND(intArrayOf(MaxIter, Npar + 3))) // initialize convergence history - - var h: DoubleTensor - var dX2 = X2 - while (!stop && settings.iteration <= MaxIter) { //--- Start Main Loop - settings.iteration += 1 - - // incremental change in parameters - h = when (Update_Type) { - 1 -> { // Marquardt - val solve = solve(JtWJ.plus(make_matrx_with_diagonal(diag(JtWJ)).div(1 / lambda)).as2D(), JtWdy) - solve.asDoubleTensor() - } - - else -> { // Quadratic and Nielsen - val solve = solve(JtWJ.plus(lm_eye(Npar).div(1 / lambda)).as2D(), JtWdy) - solve.asDoubleTensor() - } - } - - 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 - - for (i in 0 until delta_y.shape.component1()) { // floating point error; break - for (j in 0 until delta_y.shape.component2()) { - if (delta_y[i, j] == Double.POSITIVE_INFINITY || delta_y[i, j] == Double.NEGATIVE_INFINITY) { - stop = true - break - } - } - } - - settings.func_calls += 1 - - val tmp = delta_y.times(weight) - var X2_try = delta_y.as2D().transpose().dot(tmp) // Chi-squared error criteria - - val alpha = 1.0 - if (Update_Type == 2) { // Quadratic - // One step of quadratic line update in the h direction for minimum X2 - - val alpha = JtWdy.transpose().dot(h) / ( (X2_try.minus(X2)).div(2.0).plus(2 * JtWdy.transpose().dot(h)) ) - h = h.dot(alpha) - 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 - - val tmp = delta_y.times(weight) - X2_try = delta_y.as2D().transpose().dot(tmp) // Chi-squared error criteria - } - - val rho = when (Update_Type) { // Nielsen - 1 -> { - val tmp = h.transposed().dot(make_matrx_with_diagonal(diag(JtWJ)).div(1 / lambda).dot(h).plus(JtWdy)) - X2.minus(X2_try).as2D()[0, 0] / abs(tmp.as2D()).as2D()[0, 0] - } - else -> { - val tmp = h.transposed().dot(h.div(1 / lambda).plus(JtWdy)) - X2.minus(X2_try).as2D()[0, 0] / abs(tmp.as2D()).as2D()[0, 0] - } - } - - if (rho > epsilon_4) { // it IS significantly better - val dX2 = X2.minus(X2_old) - X2_old = X2 - p_old = p.copyToTensor().as2D() - y_old = y_hat.copyToTensor().as2D() - p = make_column(p_try) // accept p_try - - lm_matx_ans = lm_matx(func, t, p_old, y_old, dX2.toInt(), J, p, y_dat, weight, dp, settings) - // decrease lambda ==> Gauss-Newton method - - JtWJ = lm_matx_ans[0] - JtWdy = lm_matx_ans[1] - X2 = lm_matx_ans[2][0, 0] - y_hat = lm_matx_ans[3] - J = lm_matx_ans[4] - - lambda = when (Update_Type) { - 1 -> { // Levenberg - max(lambda / lambda_DN_fac, 1e-7); - } - 2 -> { // Quadratic - max( lambda / (1 + alpha) , 1e-7 ); - } - else -> { // Nielsen - nu = 2 - lambda * max( 1.0 / 3, 1 - (2 * rho - 1).pow(3) ) - } - } - } - else { // it IS NOT better - X2 = X2_old // do not accept p_try - if (settings.iteration % (2 * Npar) == 0 ) { // rank-1 update of Jacobian - lm_matx_ans = lm_matx(func, t, p_old, y_old,-1, J, p, y_dat, weight, dp, settings) - JtWJ = lm_matx_ans[0] - JtWdy = lm_matx_ans[1] - dX2 = lm_matx_ans[2][0, 0] - y_hat = lm_matx_ans[3] - J = lm_matx_ans[4] - } - - // increase lambda ==> gradient descent method - lambda = when (Update_Type) { - 1 -> { // Levenberg - min(lambda * lambda_UP_fac, 1e7) - } - 2 -> { // Quadratic - lambda + abs(((X2_try.as2D()[0, 0] - X2) / 2) / alpha) - } - else -> { // Nielsen - nu *= 2 - lambda * (nu / 2) - } - } - } - - if (prnt > 1) { - val chi_sq = X2 / DoF -// println("Iteration $settings | chi_sq=$chi_sq | lambda=$lambda") -// print("param: ") -// for (pn in 0 until Npar) { -// print(p[pn, 0].toString() + " ") -// } -// print("\ndp/p: ") -// for (pn in 0 until Npar) { -// print((h.as2D()[pn, 0] / p[pn, 0]).toString() + " ") -// } - resultInfo.iterations = settings.iteration - resultInfo.func_calls = settings.func_calls - resultInfo.result_chi_sq = chi_sq - resultInfo.result_lambda = lambda - resultInfo.result_parameters = p - } - - // update convergence history ... save _reduced_ Chi-square - // cvg_hst(iteration,:) = [ func_calls p' X2/DoF lambda ]; - - if (abs(JtWdy).max()!! < epsilon_1 && settings.iteration > 2) { -// println(" **** Convergence in r.h.s. (\"JtWdy\") ****") -// println(" **** epsilon_1 = $epsilon_1") - resultInfo.typeOfConvergence = LinearOpsTensorAlgebra.TypeOfConvergence.inRHS_JtWdy - resultInfo.epsilon = epsilon_1 - stop = true - } - if ((abs(h.as2D()).div(abs(p) + 1e-12)).max() < epsilon_2 && settings.iteration > 2) { -// println(" **** Convergence in Parameters ****") -// println(" **** epsilon_2 = $epsilon_2") - resultInfo.typeOfConvergence = LinearOpsTensorAlgebra.TypeOfConvergence.inParameters - resultInfo.epsilon = epsilon_2 - stop = true - } - if (X2 / DoF < epsilon_3 && settings.iteration > 2) { -// println(" **** Convergence in reduced Chi-square **** ") -// println(" **** epsilon_3 = $epsilon_3") - resultInfo.typeOfConvergence = LinearOpsTensorAlgebra.TypeOfConvergence.inReducedChi_square - resultInfo.epsilon = epsilon_3 - stop = true - } - if (settings.iteration == MaxIter) { -// println(" !! Maximum Number of Iterations Reached Without Convergence !!") - resultInfo.typeOfConvergence = LinearOpsTensorAlgebra.TypeOfConvergence.noConvergence - resultInfo.epsilon = 0.0 - print("noConvergence, MaxIter = ") - println(MaxIter) - stop = true - } - } // --- End of Main Loop - - // --- convergence achieved, find covariance and confidence intervals - - // ---- Error Analysis ---- - -// if (weight.shape.component1() == 1 || weight.variance() == 0.0) { -// weight = DoF / (delta_y.transpose().dot(delta_y)) * ones(intArrayOf(Npt, 1)) -// } - -// if (nargout > 1) { -// val redX2 = X2 / DoF -// } -// -// lm_matx_ans = lm_matx(func, t, p_old, y_old, -1, J, p, y_dat, weight, dp) -// JtWJ = lm_matx_ans[0] -// JtWdy = lm_matx_ans[1] -// X2 = lm_matx_ans[2][0, 0] -// y_hat = lm_matx_ans[3] -// J = lm_matx_ans[4] -// -// if (nargout > 2) { // standard error of parameters -// covar_p = inv(JtWJ); -// siif nagma_p = sqrt(diag(covar_p)); -// } -// -// if (nargout > 3) { // standard error of the fit -// /// sigma_y = sqrt(diag(J * covar_p * J')); // slower version of below -// sigma_y = zeros(Npnt,1); -// for i=1:Npnt -// sigma_y(i) = J(i,:) * covar_p * J(i,:)'; -// end -// sigma_y = sqrt(sigma_y); -// } -// -// if (nargout > 4) { // parameter correlation matrix -// corr_p = covar_p ./ [sigma_p*sigma_p']; -// } -// -// if (nargout > 5) { // coefficient of multiple determination -// R_sq = corr([y_dat y_hat]); -// R_sq = R_sq(1,2).^2; -// } -// -// if (nargout > 6) { // convergence history -// cvg_hst = cvg_hst(1:iteration,:); -// } - - return resultInfo - } } public val Double.Companion.tensorAlgebra: DoubleTensorAlgebra get() = DoubleTensorAlgebra 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 new file mode 100644 index 000000000..f4b50626a --- /dev/null +++ b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/LevenbergMarquardtAlgorithm.kt @@ -0,0 +1,553 @@ +/* + * Copyright 2018-2023 KMath contributors. + * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file. + */ + +package space.kscience.kmath.tensors.core + +import space.kscience.kmath.linear.transpose +import space.kscience.kmath.nd.* +import space.kscience.kmath.tensors.api.LinearOpsTensorAlgebra +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.max +import kotlin.math.min +import kotlin.math.pow +import kotlin.reflect.KFunction3 + +public fun DoubleTensorAlgebra.lm( + func: KFunction3, MutableStructure2D, LMSettings, 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): LinearOpsTensorAlgebra.LMResultInfo { + + val resultInfo = LinearOpsTensorAlgebra.LMResultInfo(0, 0, example_number, 0.0, + 0.0, p_input, LinearOpsTensorAlgebra.TypeOfConvergence.noConvergence, 0.0) + + val eps:Double = 2.2204e-16 + + val settings = LMSettings(0, 0, example_number) + settings.func_calls = 0 // running count of function evaluations + + var p = p_input + val y_dat = y_dat_input + val t = t_input + + val Npar = length(p) // number of parameters + val Npnt = length(y_dat) // number of data points + var p_old = zeros(ShapeND(intArrayOf(Npar, 1))).as2D() // previous set of parameters + var y_old = 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 X2_old = 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 + + var corr_p = 0 + var sigma_p = 0 + var sigma_y = 0 + var R_sq = 0 + var cvg_hist = 0 + + if (length(t) != length(y_dat)) { + // println("lm.m error: the length of t must equal the length of y_dat") + val length_t = length(t) + val length_y_dat = length(y_dat) + X2 = 0.0 + + corr_p = 0 + sigma_p = 0 + sigma_y = 0 + R_sq = 0 + cvg_hist = 0 + } + + var weight = weight_input + if (nargin < 5) { + weight = fromArray(ShapeND(intArrayOf(1, 1)), doubleArrayOf((y_dat.transpose().dot(y_dat)).as1D()[0])).as2D() + } + + var dp = dp_input + if (nargin < 6) { + dp = fromArray(ShapeND(intArrayOf(1, 1)), doubleArrayOf(0.001)).as2D() + } + + var p_min = p_min_input + if (nargin < 7) { + p_min = p + p_min.abs() + p_min = p_min.div(-100.0).as2D() + } + + var p_max = p_max_input + if (nargin < 8) { + p_max = p + p_max.abs() + p_max = p_max.div(100.0).as2D() + } + + var c = c_input + if (nargin < 9) { + c = fromArray(ShapeND(intArrayOf(1, 1)), doubleArrayOf(1.0)).as2D() + } + + var opts = opts_input + 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) + } + + val prnt = opts[0] // >1 intermediate results; >2 plots + val MaxIter = opts[1].toInt() // maximum number of iterations + val epsilon_1 = opts[2] // convergence tolerance for gradient + val epsilon_2 = opts[3] // convergence tolerance for parameters + val epsilon_3 = opts[4] // convergence tolerance for Chi-square + val epsilon_4 = opts[5] // determines acceptance of a L-M step + val lambda_0 = opts[6] // initial value of damping paramter, lambda + val lambda_UP_fac = opts[7] // factor for increasing lambda + val lambda_DN_fac = opts[8] // factor for decreasing lambda + val Update_Type = opts[9].toInt() // 1: Levenberg-Marquardt lambda update + // 2: Quadratic update + // 3: Nielsen's lambda update equations + + p_min = make_column(p_min) + p_max = make_column(p_max) + + if (length(make_column(dp)) == 1) { + dp = ones(ShapeND(intArrayOf(Npar, 1))).div(1 / dp[0, 0]).as2D() + } + + 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 + + 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() + // println("using uniform weights for error analysis") + } + else { + weight = make_column(weight) + weight.abs() + } + + // initialize Jacobian with finite difference calculation + var lm_matx_ans = lm_matx(func, t, p_old, y_old,1, J, p, y_dat, weight, dp, settings) + var JtWJ = lm_matx_ans[0] + var JtWdy = lm_matx_ans[1] + X2 = lm_matx_ans[2][0, 0] + var y_hat = lm_matx_ans[3] + J = lm_matx_ans[4] + + if ( abs(JtWdy).max()!! < epsilon_1 ) { +// println(" *** Your Initial Guess is Extremely Close to Optimal ***\n") +// println(" *** epsilon_1 = %e\n$epsilon_1") + stop = true + } + + var lambda = 1.0 + var nu = 1 + when (Update_Type) { + 1 -> lambda = lambda_0 // Marquardt: init'l lambda + else -> { // Quadratic and Nielsen + lambda = lambda_0 * (diag(JtWJ)).max()!! + nu = 2 + } + } + + X2_old = X2 // previous value of X2 + var cvg_hst = ones(ShapeND(intArrayOf(MaxIter, Npar + 3))) // initialize convergence history + + var h: DoubleTensor + var dX2 = X2 + while (!stop && settings.iteration <= MaxIter) { //--- Start Main Loop + settings.iteration += 1 + + // incremental change in parameters + h = when (Update_Type) { + 1 -> { // Marquardt + val solve = solve(JtWJ.plus(make_matrx_with_diagonal(diag(JtWJ)).div(1 / lambda)).as2D(), JtWdy) + solve.asDoubleTensor() + } + + else -> { // Quadratic and Nielsen + val solve = solve(JtWJ.plus(lm_eye(Npar).div(1 / lambda)).as2D(), JtWdy) + solve.asDoubleTensor() + } + } + + 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 + + for (i in 0 until delta_y.shape.component1()) { // floating point error; break + for (j in 0 until delta_y.shape.component2()) { + if (delta_y[i, j] == Double.POSITIVE_INFINITY || delta_y[i, j] == Double.NEGATIVE_INFINITY) { + stop = true + break + } + } + } + + settings.func_calls += 1 + + val tmp = delta_y.times(weight) + var X2_try = delta_y.as2D().transpose().dot(tmp) // Chi-squared error criteria + + val alpha = 1.0 + if (Update_Type == 2) { // Quadratic + // One step of quadratic line update in the h direction for minimum X2 + + val alpha = JtWdy.transpose().dot(h) / ( (X2_try.minus(X2)).div(2.0).plus(2 * JtWdy.transpose().dot(h)) ) + h = h.dot(alpha) + 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 + + val tmp = delta_y.times(weight) + X2_try = delta_y.as2D().transpose().dot(tmp) // Chi-squared error criteria + } + + val rho = when (Update_Type) { // Nielsen + 1 -> { + val tmp = h.transposed().dot(make_matrx_with_diagonal(diag(JtWJ)).div(1 / lambda).dot(h).plus(JtWdy)) + X2.minus(X2_try).as2D()[0, 0] / abs(tmp.as2D()).as2D()[0, 0] + } + else -> { + val tmp = h.transposed().dot(h.div(1 / lambda).plus(JtWdy)) + X2.minus(X2_try).as2D()[0, 0] / abs(tmp.as2D()).as2D()[0, 0] + } + } + + if (rho > epsilon_4) { // it IS significantly better + val dX2 = X2.minus(X2_old) + X2_old = X2 + p_old = p.copyToTensor().as2D() + y_old = y_hat.copyToTensor().as2D() + p = make_column(p_try) // accept p_try + + lm_matx_ans = lm_matx(func, t, p_old, y_old, dX2.toInt(), J, p, y_dat, weight, dp, settings) + // decrease lambda ==> Gauss-Newton method + + JtWJ = lm_matx_ans[0] + JtWdy = lm_matx_ans[1] + X2 = lm_matx_ans[2][0, 0] + y_hat = lm_matx_ans[3] + J = lm_matx_ans[4] + + lambda = when (Update_Type) { + 1 -> { // Levenberg + max(lambda / lambda_DN_fac, 1e-7); + } + 2 -> { // Quadratic + max( lambda / (1 + alpha) , 1e-7 ); + } + else -> { // Nielsen + nu = 2 + lambda * max( 1.0 / 3, 1 - (2 * rho - 1).pow(3) ) + } + } + } + else { // it IS NOT better + X2 = X2_old // do not accept p_try + if (settings.iteration % (2 * Npar) == 0 ) { // rank-1 update of Jacobian + lm_matx_ans = lm_matx(func, t, p_old, y_old,-1, J, p, y_dat, weight, dp, settings) + JtWJ = lm_matx_ans[0] + JtWdy = lm_matx_ans[1] + dX2 = lm_matx_ans[2][0, 0] + y_hat = lm_matx_ans[3] + J = lm_matx_ans[4] + } + + // increase lambda ==> gradient descent method + lambda = when (Update_Type) { + 1 -> { // Levenberg + min(lambda * lambda_UP_fac, 1e7) + } + 2 -> { // Quadratic + lambda + kotlin.math.abs(((X2_try.as2D()[0, 0] - X2) / 2) / alpha) + } + else -> { // Nielsen + nu *= 2 + lambda * (nu / 2) + } + } + } + + if (prnt > 1) { + val chi_sq = X2 / DoF +// println("Iteration $settings | chi_sq=$chi_sq | lambda=$lambda") +// print("param: ") +// for (pn in 0 until Npar) { +// print(p[pn, 0].toString() + " ") +// } +// print("\ndp/p: ") +// for (pn in 0 until Npar) { +// print((h.as2D()[pn, 0] / p[pn, 0]).toString() + " ") +// } + resultInfo.iterations = settings.iteration + resultInfo.func_calls = settings.func_calls + resultInfo.result_chi_sq = chi_sq + resultInfo.result_lambda = lambda + resultInfo.result_parameters = p + } + + // update convergence history ... save _reduced_ Chi-square + // cvg_hst(iteration,:) = [ func_calls p' X2/DoF lambda ]; + + if (abs(JtWdy).max()!! < epsilon_1 && settings.iteration > 2) { +// println(" **** Convergence in r.h.s. (\"JtWdy\") ****") +// println(" **** epsilon_1 = $epsilon_1") + resultInfo.typeOfConvergence = LinearOpsTensorAlgebra.TypeOfConvergence.inRHS_JtWdy + resultInfo.epsilon = epsilon_1 + stop = true + } + if ((abs(h.as2D()).div(abs(p) + 1e-12)).max() < epsilon_2 && settings.iteration > 2) { +// println(" **** Convergence in Parameters ****") +// println(" **** epsilon_2 = $epsilon_2") + resultInfo.typeOfConvergence = LinearOpsTensorAlgebra.TypeOfConvergence.inParameters + resultInfo.epsilon = epsilon_2 + stop = true + } + if (X2 / DoF < epsilon_3 && settings.iteration > 2) { +// println(" **** Convergence in reduced Chi-square **** ") +// println(" **** epsilon_3 = $epsilon_3") + resultInfo.typeOfConvergence = LinearOpsTensorAlgebra.TypeOfConvergence.inReducedChi_square + resultInfo.epsilon = epsilon_3 + stop = true + } + if (settings.iteration == MaxIter) { +// println(" !! Maximum Number of Iterations Reached Without Convergence !!") + resultInfo.typeOfConvergence = LinearOpsTensorAlgebra.TypeOfConvergence.noConvergence + resultInfo.epsilon = 0.0 + stop = true + } + } // --- End of Main Loop + return resultInfo +} + +public data class LMSettings ( + var iteration:Int, + var func_calls: Int, + var example_number:Int +) + +/* matrix -> column of all elemnets */ +public 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()) { + for (j in 0 until tensor.shape.component2()) { + buffer[i * tensor.shape.component2() + j] = tensor[i, j] + } + } + val column = BroadcastDoubleTensorAlgebra.fromArray(ShapeND(shape), buffer).as2D() + return column +} + +/* column length */ +public fun length(column: MutableStructure2D) : Int { + return column.shape.component1() +} + +public 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]) + } + } +} + +public fun abs(input: MutableStructure2D): MutableStructure2D { + val tensor = BroadcastDoubleTensorAlgebra.ones( + ShapeND( + intArrayOf( + input.shape.component1(), + input.shape.component2() + ) + ) + ).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]) + } + } + return tensor +} + +public 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] + } + return tensor +} + +public 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) { + tensor[i, i] = column[i, 0] + } + return tensor +} + +public 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 { + val a_sizeX = a.shape.component1() + val a_sizeY = a.shape.component2() + val b_sizeX = b.shape.component1() + val b_sizeY = b.shape.component2() + val tensor = BroadcastDoubleTensorAlgebra.zeros(ShapeND(intArrayOf(max(a_sizeX, b_sizeX), max(a_sizeY, b_sizeY)))).as2D() + for (i in 0 until tensor.shape.component1()) { + for (j in 0 until tensor.shape.component2()) { + if (i < a_sizeX && i < b_sizeX && j < a_sizeY && j < b_sizeY) { + tensor[i, j] = max(a[i, j], b[i, j]) + } + else if (i < a_sizeX && j < a_sizeY) { + tensor[i, j] = a[i, j] + } + else { + tensor[i, j] = b[i, j] + } + } + } + return tensor +} + +public 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() + val b_sizeY = b.shape.component2() + val tensor = BroadcastDoubleTensorAlgebra.zeros(ShapeND(intArrayOf(max(a_sizeX, b_sizeX), max(a_sizeY, b_sizeY)))).as2D() + for (i in 0 until tensor.shape.component1()) { + for (j in 0 until tensor.shape.component2()) { + if (i < a_sizeX && i < b_sizeX && j < a_sizeY && j < b_sizeY) { + tensor[i, j] = min(a[i, j], b[i, j]) + } + else if (i < a_sizeX && j < a_sizeY) { + tensor[i, j] = a[i, j] + } + else { + tensor[i, j] = b[i, j] + } + } + } + return tensor +} + +public 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) { + idx += (i + 1.0) + } + } + if (idx.size > 0) { + return BroadcastDoubleTensorAlgebra.fromArray(ShapeND(intArrayOf(idx.size, 1)), idx.toDoubleArray()).as2D() + } + return null +} + +public fun feval(func: (MutableStructure2D, MutableStructure2D, LMSettings) -> MutableStructure2D, + t: MutableStructure2D, p: MutableStructure2D, settings: LMSettings) + : MutableStructure2D +{ + return func(t, p, settings) +} + +public fun lm_matx(func: (MutableStructure2D, MutableStructure2D, LMSettings) -> 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> +{ + // default: dp = 0.001 + + 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 + + var J = J_input + + if (settings.iteration % (2 * Npar) == 0 || dX2 > 0) { + J = lm_FD_J(func, t, p, y_hat, dp, settings).as2D() // finite difference + } + else { + J = lm_Broyden_J(p_old, y_old, J, p, y_hat).as2D() // rank-1 update + } + + val delta_y = y_dat.minus(y_hat) + + val Chi_sq = delta_y.transposed().dot( delta_y.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(delta_y) ).as2D() + + return arrayOf(JtWJ,JtWdy,Chi_sq,y_hat,J) +} + +public fun lm_Broyden_J(p_old: MutableStructure2D, y_old: MutableStructure2D, J_input: MutableStructure2D, + p: MutableStructure2D, y: MutableStructure2D): MutableStructure2D { + var J = J_input.copyToTensor() + + val h = p.minus(p_old) + val increase = y.minus(y_old).minus( J.dot(h) ).dot(h.transposed()).div( (h.transposed().dot(h)).as2D()[0, 0] ) + J = J.plus(increase) + + return J.as2D() +} + +public fun lm_FD_J(func: (MutableStructure2D, MutableStructure2D, settings: LMSettings) -> MutableStructure2D, + t: MutableStructure2D, p: MutableStructure2D, y: MutableStructure2D, + dp: MutableStructure2D, settings: LMSettings): MutableStructure2D { + // 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() + + for (j in 0 until n) { + + del[j, 0] = dp[j, 0] * (1 + kotlin.math.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) { + val y1 = feval(func, t, p, settings) + settings.func_calls += 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] + } + } + else { + // Do tests for it + 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]) + } + settings.func_calls += 1 + } + } + + p[j, 0] = ps[j, 0] // restore p(j) + } + + return J.as2D() +} diff --git a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/linUtils.kt b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/linUtils.kt index 6e5456c62..086c69e49 100644 --- a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/linUtils.kt +++ b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/linUtils.kt @@ -9,12 +9,6 @@ import space.kscience.kmath.nd.* import space.kscience.kmath.operations.invoke import space.kscience.kmath.structures.* import space.kscience.kmath.tensors.core.* -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 @@ -607,225 +601,4 @@ internal fun MutableStructure2D.svdGolubKahanHelper(u: MutableStructure2 u[j, i] = this[j, i] } } -} - -data class LMSettings ( - var iteration:Int, - var func_calls: Int, - var example_number:Int -) - -/* matrix -> column of all elemnets */ -fun make_column(tensor: MutableStructure2D) : MutableStructure2D { - val shape = intArrayOf(tensor.shape.component1() * tensor.shape.component2(), 1) - var buffer = DoubleArray(tensor.shape.component1() * tensor.shape.component2()) - for (i in 0 until tensor.shape.component1()) { - for (j in 0 until tensor.shape.component2()) { - buffer[i * tensor.shape.component2() + j] = tensor[i, j] - } - } - var column = BroadcastDoubleTensorAlgebra.fromArray(ShapeND(shape), buffer).as2D() - return column -} - -/* column length */ -fun length(column: MutableStructure2D) : Int { - return column.shape.component1() -} - -fun MutableStructure2D.abs() { - for (i in 0 until this.shape.component1()) { - for (j in 0 until this.shape.component2()) { - this[i, j] = abs(this[i, j]) - } - } -} - -fun abs(input: MutableStructure2D): MutableStructure2D { - val tensor = BroadcastDoubleTensorAlgebra.ones( - ShapeND( - intArrayOf( - input.shape.component1(), - input.shape.component2() - ) - ) - ).as2D() - for (i in 0 until tensor.shape.component1()) { - for (j in 0 until tensor.shape.component2()) { - tensor[i, j] = abs(input[i, j]) - } - } - return tensor -} - -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] - } - return tensor -} - -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) { - tensor[i, i] = column[i, 0] - } - return tensor -} - -fun lm_eye(size: Int): MutableStructure2D { - val column = BroadcastDoubleTensorAlgebra.ones(ShapeND(intArrayOf(size, 1))).as2D() - return make_matrx_with_diagonal(column) -} - -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() - val b_sizeY = b.shape.component2() - val tensor = BroadcastDoubleTensorAlgebra.zeros(ShapeND(intArrayOf(max(a_sizeX, b_sizeX), max(a_sizeY, b_sizeY)))).as2D() - for (i in 0 until tensor.shape.component1()) { - for (j in 0 until tensor.shape.component2()) { - if (i < a_sizeX && i < b_sizeX && j < a_sizeY && j < b_sizeY) { - tensor[i, j] = max(a[i, j], b[i, j]) - } - else if (i < a_sizeX && j < a_sizeY) { - tensor[i, j] = a[i, j] - } - else { - tensor[i, j] = b[i, j] - } - } - } - return tensor -} - -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() - val b_sizeY = b.shape.component2() - val tensor = BroadcastDoubleTensorAlgebra.zeros(ShapeND(intArrayOf(max(a_sizeX, b_sizeX), max(a_sizeY, b_sizeY)))).as2D() - for (i in 0 until tensor.shape.component1()) { - for (j in 0 until tensor.shape.component2()) { - if (i < a_sizeX && i < b_sizeX && j < a_sizeY && j < b_sizeY) { - tensor[i, j] = min(a[i, j], b[i, j]) - } - else if (i < a_sizeX && j < a_sizeY) { - tensor[i, j] = a[i, j] - } - else { - tensor[i, j] = b[i, j] - } - } - } - return tensor -} - -fun get_zero_indices(column: MutableStructure2D, epsilon: Double = 0.000001): MutableStructure2D? { - var idx = emptyArray() - for (i in 0 until column.shape.component1()) { - if (abs(column[i, 0]) > epsilon) { - idx += (i + 1.0) - } - } - if (idx.size > 0) { - return BroadcastDoubleTensorAlgebra.fromArray(ShapeND(intArrayOf(idx.size, 1)), idx.toDoubleArray()).as2D() - } - return null -} - -fun feval(func: (MutableStructure2D, MutableStructure2D, LMSettings) -> MutableStructure2D, - t: MutableStructure2D, p: MutableStructure2D, settings: LMSettings) - : MutableStructure2D -{ - return func(t, p, settings) -} - -fun lm_matx(func: (MutableStructure2D, MutableStructure2D, LMSettings) -> 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> -{ - // default: dp = 0.001 - - 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 - - var J = J_input - - if (settings.iteration % (2 * Npar) == 0 || dX2 > 0) { - J = lm_FD_J(func, t, p, y_hat, dp, settings).as2D() // finite difference - } - else { - J = lm_Broyden_J(p_old, y_old, J, p, y_hat).as2D() // rank-1 update - } - - val delta_y = y_dat.minus(y_hat) - - val Chi_sq = delta_y.transposed().dot( delta_y.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(delta_y) ).as2D() - - return arrayOf(JtWJ,JtWdy,Chi_sq,y_hat,J) -} - -fun lm_Broyden_J(p_old: MutableStructure2D, y_old: MutableStructure2D, J_input: MutableStructure2D, - p: MutableStructure2D, y: MutableStructure2D): MutableStructure2D { - var J = J_input.copyToTensor() - - val h = p.minus(p_old) - val increase = y.minus(y_old).minus( J.dot(h) ).dot(h.transposed()).div( (h.transposed().dot(h)).as2D()[0, 0] ) - J = J.plus(increase) - - return J.as2D() -} - -fun lm_FD_J(func: (MutableStructure2D, MutableStructure2D, settings: LMSettings) -> MutableStructure2D, - t: MutableStructure2D, p: MutableStructure2D, y: MutableStructure2D, - dp: MutableStructure2D, settings: LMSettings): MutableStructure2D { - // 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() - - for (j in 0 until n) { - - 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) { - val y1 = feval(func, t, p, settings) - settings.func_calls += 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] - } - } - else { - // Do tests for it - 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]) - } - settings.func_calls += 1 - } - } - - p[j, 0] = ps[j, 0] // restore p(j) - } - - return J.as2D() -} +} \ No newline at end of file diff --git a/kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleTensorAlgebra.kt b/kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleTensorAlgebra.kt index 5aea9b879..7222fc7a6 100644 --- a/kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleTensorAlgebra.kt +++ b/kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleTensorAlgebra.kt @@ -8,11 +8,7 @@ package space.kscience.kmath.tensors.core import space.kscience.kmath.nd.* import space.kscience.kmath.operations.invoke -import space.kscience.kmath.tensors.api.LinearOpsTensorAlgebra -import space.kscience.kmath.tensors.core.internal.LMSettings import space.kscience.kmath.testutils.assertBufferEquals -import kotlin.math.roundToInt -import kotlin.math.roundToLong import kotlin.test.Test import kotlin.test.assertEquals import kotlin.test.assertFalse @@ -209,100 +205,4 @@ internal class TestDoubleTensorAlgebra { assertTrue { ShapeND(5, 5) contentEquals res.shape } assertEquals(2.0, res[4, 4]) } - - @Test - fun testLM() = DoubleTensorAlgebra { - fun lm_func(t: MutableStructure2D, p: MutableStructure2D, settings: LMSettings): 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]))) - ) - } - 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]) - } - - return y_hat.as2D() - } - - val lm_matx_y_dat = doubleArrayOf( - 19.6594, 18.6096, 17.6792, 17.2747, 16.3065, 17.1458, 16.0467, 16.7023, 15.7809, 15.9807, - 14.7620, 15.1128, 16.0973, 15.1934, 15.8636, 15.4763, 15.6860, 15.1895, 15.3495, 16.6054, - 16.2247, 15.9854, 16.1421, 17.0960, 16.7769, 17.1997, 17.2767, 17.5882, 17.5378, 16.7894, - 17.7648, 18.2512, 18.1581, 16.7037, 17.8475, 17.9081, 18.3067, 17.9632, 18.2817, 19.1427, - 18.8130, 18.5658, 18.0056, 18.4607, 18.5918, 18.2544, 18.3731, 18.7511, 19.3181, 17.3066, - 17.9632, 19.0513, 18.7528, 18.2928, 18.5967, 17.8567, 17.7859, 18.4016, 18.9423, 18.4959, - 17.8000, 18.4251, 17.7829, 17.4645, 17.5221, 17.3517, 17.4637, 17.7563, 16.8471, 17.4558, - 17.7447, 17.1487, 17.3183, 16.8312, 17.7551, 17.0942, 15.6093, 16.4163, 15.3755, 16.6725, - 16.2332, 16.2316, 16.2236, 16.5361, 15.3721, 15.3347, 15.5815, 15.6319, 14.4538, 14.6044, - 14.7665, 13.3718, 15.0587, 13.8320, 14.7873, 13.6824, 14.2579, 14.2154, 13.5818, 13.8157 - ) - - var example_number = 1 - val p_init = BroadcastDoubleTensorAlgebra.fromArray( - ShapeND(intArrayOf(4, 1)), doubleArrayOf(5.0, 2.0, 0.2, 10.0) - ).as2D() - - var t = ones(ShapeND(intArrayOf(100, 1))).as2D() - for (i in 0 until 100) { - t[i, 0] = t[i, 0] * (i + 1) - } - - val y_dat = BroadcastDoubleTensorAlgebra.fromArray( - ShapeND(intArrayOf(100, 1)), lm_matx_y_dat - ).as2D() - - val weight = BroadcastDoubleTensorAlgebra.fromArray( - ShapeND(intArrayOf(1, 1)), DoubleArray(1) { 4.0 } - ).as2D() - - val dp = BroadcastDoubleTensorAlgebra.fromArray( - ShapeND(intArrayOf(1, 1)), DoubleArray(1) { -0.01 } - ).as2D() - - val p_min = BroadcastDoubleTensorAlgebra.fromArray( - ShapeND(intArrayOf(4, 1)), doubleArrayOf(-50.0, -20.0, -2.0, -100.0) - ).as2D() - - val p_max = BroadcastDoubleTensorAlgebra.fromArray( - ShapeND(intArrayOf(4, 1)), doubleArrayOf(50.0, 20.0, 2.0, 100.0) - ).as2D() - - val consts = BroadcastDoubleTensorAlgebra.fromArray( - ShapeND(intArrayOf(1, 1)), doubleArrayOf(0.0) - ).as2D() - - 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(::lm_func, 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(result.typeOfConvergence, LinearOpsTensorAlgebra.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() - 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]) - ).as2D() - assertEquals(expectedParameters[0, 0], receivedParameters[0, 0]) - assertEquals(expectedParameters[1, 0], receivedParameters[1, 0]) - assertEquals(expectedParameters[2, 0], receivedParameters[2, 0]) - assertEquals(expectedParameters[3, 0], receivedParameters[3, 0]) - } } 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 f554a742c..db8bc9d20 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 @@ -15,7 +15,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.internal.LMSettings import kotlin.math.roundToLong import kotlin.test.Test import kotlin.test.assertEquals