From 976714475e5a3eb2ff0175b2655ac0d794ecea71 Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Fri, 28 Jul 2023 20:56:31 +0300 Subject: [PATCH] levenbergMarquardt cleanup --- .../core/LevenbergMarquardtAlgorithm.kt | 10 ++-- .../kmath/tensors/core/internal/linUtils.kt | 49 ++++++++++--------- 2 files changed, 31 insertions(+), 28 deletions(-) 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 72752f855..fc87ad1f3 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 @@ -120,8 +120,6 @@ public fun DoubleTensorAlgebra.levenbergMarquardt(inputData: LMInput): LMResultI 0.0, inputData.startParameters, TypeOfConvergence.NoConvergence ) - val eps = 2.2204e-16 - val settings = LMSettings(0, 0, inputData.exampleNumber) settings.funcCalls = 0 // running count of function evaluations @@ -214,7 +212,7 @@ public fun DoubleTensorAlgebra.levenbergMarquardt(inputData: LMInput): LMResultI stop = true } - var lambda = 1.0 + var lambda: Double var nu = 1 if (updateType == 1) { @@ -273,8 +271,8 @@ public fun DoubleTensorAlgebra.levenbergMarquardt(inputData: LMInput): LMResultI 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 - x2) / 2.0 + 2 * (jtWdy.transpose() dot h)) - h = h dot alpha + val alphaTensor = (jtWdy.transpose() dot h) / ((X2Try - x2) / 2.0 + 2 * (jtWdy.transpose() dot h)) + h = h dot alphaTensor pTry = (p + h).as2D() // update only [idx] elements pTry = smallestElementComparison( largestElementComparison(minParameters, pTry), @@ -388,7 +386,7 @@ public fun DoubleTensorAlgebra.levenbergMarquardt(inputData: LMInput): LMResultI resultInfo.typeOfConvergence = TypeOfConvergence.InGradient stop = true } - if ((abs(h.as2D()).div(abs(p) + 1e-12)).max() < epsilon2 && settings.iteration > 2) { + if ((abs(h.as2D()) / (abs(p) + 1e-12)).max() < epsilon2 && settings.iteration > 2) { resultInfo.typeOfConvergence = TypeOfConvergence.InParameters stop = true } 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 086c69e49..7a96001c6 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 @@ -7,7 +7,10 @@ package space.kscience.kmath.tensors.core.internal import space.kscience.kmath.nd.* import space.kscience.kmath.operations.invoke -import space.kscience.kmath.structures.* +import space.kscience.kmath.structures.DoubleBuffer +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 kotlin.math.abs import kotlin.math.max @@ -329,14 +332,16 @@ private fun SIGN(a: Double, b: Double): Double { return -abs(a) } -internal fun MutableStructure2D.svdGolubKahanHelper(u: MutableStructure2D, w: BufferedTensor, - v: MutableStructure2D, iterations: Int, epsilon: Double) { +internal fun MutableStructure2D.svdGolubKahanHelper( + u: MutableStructure2D, w: BufferedTensor, + v: MutableStructure2D, iterations: Int, epsilon: Double, +) { val shape = this.shape val m = shape.component1() val n = shape.component2() - var f = 0.0 + var f: Double val rv1 = DoubleArray(n) - var s = 0.0 + var s: Double var scale = 0.0 var anorm = 0.0 var g = 0.0 @@ -362,10 +367,10 @@ internal fun MutableStructure2D.svdGolubKahanHelper(u: MutableStructure2 s += this[k, i] * this[k, i] } f = this[i, i] - if (f >= 0) { - g = (-1) * abs(sqrt(s)) + g = if (f >= 0) { + -abs(sqrt(s)) } else { - g = abs(sqrt(s)) + abs(sqrt(s)) } val h = f * g - s this[i, i] = f - g @@ -402,10 +407,10 @@ internal fun MutableStructure2D.svdGolubKahanHelper(u: MutableStructure2 s += this[i, k] * this[i, k] } f = this[i, l] - if (f >= 0) { - g = (-1) * abs(sqrt(s)) + g = if (f >= 0) { + -abs(sqrt(s)) } else { - g = abs(sqrt(s)) + abs(sqrt(s)) } val h = f * g - s this[i, l] = f - g @@ -457,7 +462,7 @@ internal fun MutableStructure2D.svdGolubKahanHelper(u: MutableStructure2 for (i in min(n, m) - 1 downTo 0) { l = i + 1 - g = wBuffer[wStart + i] + g = wBuffer[wStart + i] for (j in l until n) { this[i, j] = 0.0 } @@ -484,13 +489,13 @@ internal fun MutableStructure2D.svdGolubKahanHelper(u: MutableStructure2 this[i, i] += 1.0 } - var flag = 0 + var flag: Int var nm = 0 - var c = 0.0 - var h = 0.0 - var y = 0.0 - var z = 0.0 - var x = 0.0 + var c: Double + var h: Double + var y: Double + var z: Double + var x: Double for (k in n - 1 downTo 0) { for (its in 1 until iterations) { flag = 1 @@ -531,7 +536,7 @@ internal fun MutableStructure2D.svdGolubKahanHelper(u: MutableStructure2 } } - z = wBuffer[wStart + k] + z = wBuffer[wStart + k] if (l == k) { if (z < 0.0) { wBuffer[wStart + k] = -z @@ -541,9 +546,9 @@ internal fun MutableStructure2D.svdGolubKahanHelper(u: MutableStructure2 break } - x = wBuffer[wStart + l] + x = wBuffer[wStart + l] nm = k - 1 - y = wBuffer[wStart + nm] + y = wBuffer[wStart + nm] g = rv1[nm] h = rv1[k] f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y) @@ -552,7 +557,7 @@ internal fun MutableStructure2D.svdGolubKahanHelper(u: MutableStructure2 c = 1.0 s = 1.0 - var i = 0 + var i: Int for (j in l until nm + 1) { i = j + 1 g = rv1[i]