From 19c1af18743b899e566c254a121cd4a961d0ca59 Mon Sep 17 00:00:00 2001 From: Margarita Lashina Date: Wed, 3 May 2023 21:25:30 +0300 Subject: [PATCH] added helper functions for levenberg-marquardt algorithm --- .../kmath/tensors/core/internal/linUtils.kt | 228 ++++++++++++++++++ 1 file changed, 228 insertions(+) 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 cf3697e76..c559803ec 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 @@ -12,7 +12,14 @@ import space.kscience.kmath.structures.IntBuffer import space.kscience.kmath.structures.asBuffer import space.kscience.kmath.structures.indices 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 import kotlin.math.sqrt @@ -308,3 +315,224 @@ internal fun DoubleTensorAlgebra.svdHelper( matrixV.source[i] = vBuffer[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() +}