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 new file mode 100644 index 000000000..e6f575262 --- /dev/null +++ b/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StaticLm/staticDifficultTest.kt @@ -0,0 +1,90 @@ +/* + * 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.LevenbergMarquardt.StaticLm + +import space.kscience.kmath.nd.ShapeND +import space.kscience.kmath.nd.as2D +import space.kscience.kmath.nd.component1 +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.LMInput +import space.kscience.kmath.tensors.core.levenbergMarquardt +import kotlin.math.roundToInt + +fun main() { + val NData = 200 + var t_example = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(NData, 1))).as2D() + for (i in 0 until NData) { + t_example[i, 0] = t_example[i, 0] * (i + 1) - 104 + } + + val Nparams = 15 + var p_example = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1))).as2D() + for (i in 0 until Nparams) { + p_example[i, 0] = p_example[i, 0] + i - 25 + } + + val exampleNumber = 1 + + 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) { + p_init[i, 0] = (p_example[i, 0] + 0.9) + } + + var t = t_example + val y_dat = y_hat + val weight = 1.0 / Nparams * 1.0 - 0.085 + val dp = BroadcastDoubleTensorAlgebra.fromArray( + ShapeND(intArrayOf(1, 1)), DoubleArray(1) { -0.01 } + ).as2D() + var p_min = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1))) + p_min = p_min.div(1.0 / -50.0) + val p_max = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1))) + p_min = p_min.div(1.0 / 50.0) + val opts = doubleArrayOf(3.0, 10000.0, 1e-6, 1e-6, 1e-6, 1e-6, 1e-2, 11.0, 9.0, 1.0) +// val opts = doubleArrayOf(3.0, 10000.0, 1e-6, 1e-6, 1e-6, 1e-6, 1e-3, 11.0, 9.0, 1.0) + + val inputData = LMInput(::funcDifficultForLm, + p_init.as2D(), + t, + y_dat, + weight, + dp, + p_min.as2D(), + p_max.as2D(), + opts[1].toInt(), + doubleArrayOf(opts[2], opts[3], opts[4], opts[5]), + doubleArrayOf(opts[6], opts[7], opts[8]), + opts[9].toInt(), + 10, + 1) + + val result = DoubleTensorAlgebra.levenbergMarquardt(inputData) + + println("Parameters:") + 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.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 + println("$x $y") + } + + println("Сhi_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 new file mode 100644 index 000000000..507943031 --- /dev/null +++ b/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StaticLm/staticEasyTest.kt @@ -0,0 +1,57 @@ +/* + * 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.LevenbergMarquardt.StaticLm + +import space.kscience.kmath.nd.ShapeND +import space.kscience.kmath.nd.as2D +import space.kscience.kmath.nd.component1 +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.LMInput +import space.kscience.kmath.tensors.core.levenbergMarquardt +import kotlin.math.roundToInt + +fun main() { + val startedData = getStartDataForFuncEasy() + val inputData = LMInput(::funcEasyForLm, + DoubleTensorAlgebra.ones(ShapeND(intArrayOf(4, 1))).as2D(), + startedData.t, + startedData.y_dat, + startedData.weight, + startedData.dp, + startedData.p_min, + startedData.p_max, + startedData.opts[1].toInt(), + doubleArrayOf(startedData.opts[2], startedData.opts[3], startedData.opts[4], startedData.opts[5]), + doubleArrayOf(startedData.opts[6], startedData.opts[7], startedData.opts[8]), + startedData.opts[9].toInt(), + 10, + startedData.example_number) + + val result = DoubleTensorAlgebra.levenbergMarquardt(inputData) + + println("Parameters:") + 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.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 + println("$x $y") + } + + println("Сhi_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 new file mode 100644 index 000000000..0659db103 --- /dev/null +++ b/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StaticLm/staticMiddleTest.kt @@ -0,0 +1,88 @@ +/* + * 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.LevenbergMarquardt.StaticLm + +import space.kscience.kmath.nd.ShapeND +import space.kscience.kmath.nd.as2D +import space.kscience.kmath.nd.component1 +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.LMInput +import space.kscience.kmath.tensors.core.levenbergMarquardt +import kotlin.math.roundToInt +fun main() { + val NData = 100 + var t_example = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(NData, 1))).as2D() + for (i in 0 until NData) { + t_example[i, 0] = t_example[i, 0] * (i + 1) + } + + val Nparams = 20 + var p_example = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1))).as2D() + for (i in 0 until Nparams) { + p_example[i, 0] = p_example[i, 0] + i - 25 + } + + val exampleNumber = 1 + + 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) + } + + var t = t_example + val y_dat = y_hat + val weight = 1.0 + val dp = BroadcastDoubleTensorAlgebra.fromArray( + ShapeND(intArrayOf(1, 1)), DoubleArray(1) { -0.01 } + ).as2D() + var p_min = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1))) + p_min = p_min.div(1.0 / -50.0) + val p_max = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1))) + p_min = p_min.div(1.0 / 50.0) + val opts = doubleArrayOf(3.0, 7000.0, 1e-5, 1e-5, 1e-5, 1e-5, 1e-5, 11.0, 9.0, 1.0) + + val inputData = LMInput(::funcMiddleForLm, + p_init.as2D(), + t, + y_dat, + weight, + dp, + p_min.as2D(), + p_max.as2D(), + opts[1].toInt(), + doubleArrayOf(opts[2], opts[3], opts[4], opts[5]), + doubleArrayOf(opts[6], opts[7], opts[8]), + opts[9].toInt(), + 10, + 1) + + val result = DoubleTensorAlgebra.levenbergMarquardt(inputData) + + println("Parameters:") + 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.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 + println("$x $y") + } + + println("Сhi_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 new file mode 100644 index 000000000..b2818ef2a --- /dev/null +++ b/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StreamingLm/streamLm.kt @@ -0,0 +1,68 @@ +/* + * 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.LevenbergMarquardt.StreamingLm + +import kotlinx.coroutines.delay +import kotlinx.coroutines.flow.* +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.LMInput +import space.kscience.kmath.tensors.core.levenbergMarquardt +import kotlin.random.Random +import kotlin.reflect.KFunction3 + +fun streamLm(lm_func: (MutableStructure2D, MutableStructure2D, Int) -> (MutableStructure2D), + startData: StartDataLm, launchFrequencyInMs: Long, numberOfLaunches: Int): Flow> = flow{ + + var example_number = startData.example_number + var p_init = startData.p_init + var t = startData.t + var y_dat = startData.y_dat + val weight = startData.weight + val dp = startData.dp + val p_min = startData.p_min + val p_max = startData.p_max + val opts = startData.opts + + var steps = numberOfLaunches + val isEndless = (steps <= 0) + + val inputData = LMInput(lm_func, + p_init, + t, + y_dat, + weight, + dp, + p_min, + p_max, + opts[1].toInt(), + doubleArrayOf(opts[2], opts[3], opts[4], opts[5]), + doubleArrayOf(opts[6], opts[7], opts[8]), + opts[9].toInt(), + 10, + example_number) + + while (isEndless || steps > 0) { + val result = DoubleTensorAlgebra.levenbergMarquardt(inputData) + emit(result.resultParameters) + delay(launchFrequencyInMs) + inputData.realValues = generateNewYDat(y_dat, 0.1) + inputData.startParameters = result.resultParameters + if (!isEndless) steps -= 1 + } +} + +fun generateNewYDat(y_dat: MutableStructure2D, delta: Double): MutableStructure2D{ + val n = y_dat.shape.component1() + val y_dat_new = zeros(ShapeND(intArrayOf(n, 1))).as2D() + for (i in 0 until n) { + val randomEps = Random.nextDouble(delta + delta) - delta + y_dat_new[i, 0] = y_dat[i, 0] + randomEps + } + return y_dat_new +} \ No newline at end of file diff --git a/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StreamingLm/streamingLmTest.kt b/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StreamingLm/streamingLmTest.kt new file mode 100644 index 000000000..c9dd5029e --- /dev/null +++ b/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StreamingLm/streamingLmTest.kt @@ -0,0 +1,32 @@ +/* + * 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.LevenbergMarquardt.StreamingLm + +import space.kscience.kmath.nd.* +import space.kscience.kmath.tensors.LevenbergMarquardt.* +import kotlin.math.roundToInt + +suspend fun main(){ + val startData = getStartDataForFuncDifficult() + // Создание потока: + val lmFlow = streamLm(::funcDifficultForLm, startData, 0, 100) + var initialTime = System.currentTimeMillis() + var lastTime: Long + val launches = mutableListOf() + // Запуск потока + lmFlow.collect { parameters -> + lastTime = System.currentTimeMillis() + launches.add(lastTime - initialTime) + initialTime = lastTime + for (i in 0 until parameters.shape.component1()) { + val x = (parameters[i, 0] * 10000).roundToInt() / 10000.0 + print("$x ") + if (i == parameters.shape.component1() - 1) println() + } + } + + println("Average without first is: ${launches.subList(1, launches.size - 1).average()}") +} \ No newline at end of file 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 new file mode 100644 index 000000000..7ccb37ed0 --- /dev/null +++ b/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/functionsToOptimize.kt @@ -0,0 +1,222 @@ +/* + * 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.LevenbergMarquardt + +import space.kscience.kmath.nd.MutableStructure2D +import space.kscience.kmath.nd.ShapeND +import space.kscience.kmath.nd.as2D +import space.kscience.kmath.nd.component1 +import space.kscience.kmath.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.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.asDoubleTensor + +public data class StartDataLm ( + var lm_matx_y_dat: MutableStructure2D, + var example_number: Int, + var p_init: MutableStructure2D, + var t: MutableStructure2D, + var y_dat: MutableStructure2D, + var weight: Double, + var dp: MutableStructure2D, + var p_min: MutableStructure2D, + var p_max: MutableStructure2D, + var consts: MutableStructure2D, + var opts: DoubleArray +) + +fun funcEasyForLm(t: MutableStructure2D, p: MutableStructure2D, exampleNumber: Int): MutableStructure2D { + val m = t.shape.component1() + var y_hat = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf(m, 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 (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, exampleNumber: Int): MutableStructure2D { + val m = t.shape.component1() + var y_hat = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf (m, 1))) + + val mt = t.max() + for(i in 0 until p.shape.component1()){ + y_hat += (t.times(1.0 / mt)).times(p[i, 0]) + } + + for(i in 0 until 5){ + y_hat = funcEasyForLm(y_hat.as2D(), p, exampleNumber).asDoubleTensor() + } + + return y_hat.as2D() +} + +fun funcDifficultForLm(t: MutableStructure2D, p: MutableStructure2D, exampleNumber: Int): MutableStructure2D { + val m = t.shape.component1() + 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]) ) + } + + 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() + for (i in 0 until NData) { + t_example[i, 0] = t_example[i, 0] * (i + 1) - 104 + } + + val Nparams = 15 + var p_example = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1))).as2D() + for (i in 0 until Nparams) { + p_example[i, 0] = p_example[i, 0] + i - 25 + } + + val exampleNumber = 1 + + 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) { + p_init[i, 0] = (p_example[i, 0] + 0.9) + } + + var t = t_example + val y_dat = y_hat + val weight = 1.0 / Nparams * 1.0 - 0.085 + val dp = BroadcastDoubleTensorAlgebra.fromArray( + ShapeND(intArrayOf(1, 1)), DoubleArray(1) { -0.01 } + ).as2D() + var p_min = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1))) + p_min = p_min.div(1.0 / -50.0) + val p_max = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1))) + p_min = p_min.div(1.0 / 50.0) + val consts = BroadcastDoubleTensorAlgebra.fromArray( + ShapeND(intArrayOf(1, 1)), doubleArrayOf(0.0) + ).as2D() + val opts = doubleArrayOf(3.0, 10000.0, 1e-2, 1e-3, 1e-2, 1e-2, 1e-2, 11.0, 9.0, 1.0) + + return StartDataLm(y_dat, 1, p_init, t, y_dat, weight, dp, p_min.as2D(), p_max.as2D(), consts, opts) +} + +fun getStartDataForFuncMiddle(): StartDataLm { + val NData = 100 + var t_example = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(NData, 1))).as2D() + for (i in 0 until NData) { + t_example[i, 0] = t_example[i, 0] * (i + 1) + } + + val Nparams = 20 + var p_example = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1))).as2D() + for (i in 0 until Nparams) { + p_example[i, 0] = p_example[i, 0] + i - 25 + } + + val exampleNumber = 1 + + 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] + 10.0) + } + var t = t_example + val y_dat = y_hat + val weight = 1.0 + val dp = BroadcastDoubleTensorAlgebra.fromArray( + ShapeND(intArrayOf(1, 1)), DoubleArray(1) { -0.01 } + ).as2D() + var p_min = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1))) + p_min = p_min.div(1.0 / -50.0) + val p_max = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1))) + p_min = p_min.div(1.0 / 50.0) + val consts = BroadcastDoubleTensorAlgebra.fromArray( + ShapeND(intArrayOf(1, 1)), doubleArrayOf(0.0) + ).as2D() + val opts = doubleArrayOf(3.0, 10000.0, 1e-5, 1e-5, 1e-5, 1e-5, 1e-5, 11.0, 9.0, 1.0) + + var example_number = 1 + + return StartDataLm(y_dat, example_number, p_init, t, y_dat, weight, dp, p_min.as2D(), p_max.as2D(), consts, opts) +} + +fun getStartDataForFuncEasy(): StartDataLm { + 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 = DoubleTensorAlgebra.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 = 4.0 + + 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) + + return StartDataLm(y_dat, example_number, p_init, t, y_dat, weight, dp, p_min, p_max, consts, opts) +} \ No newline at end of file diff --git a/kmath-tensors/build.gradle.kts b/kmath-tensors/build.gradle.kts index 79c39bae7..2497314f0 100644 --- a/kmath-tensors/build.gradle.kts +++ b/kmath-tensors/build.gradle.kts @@ -4,7 +4,13 @@ plugins { kscience{ jvm() - js() + js { + browser { + testTask { + useMocha().timeout = "0" + } + } + } native() dependencies { 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 faff2eb80..f2c7f1821 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 @@ -5,6 +5,7 @@ package space.kscience.kmath.tensors.api +import space.kscience.kmath.nd.MutableStructure2D import space.kscience.kmath.nd.StructureND import space.kscience.kmath.operations.Field @@ -103,4 +104,11 @@ public interface LinearOpsTensorAlgebra> : TensorPartialDivision */ public fun symEig(structureND: StructureND): Pair, StructureND> + /** Returns the solution to the equation Ax = B for the square matrix A as `input1` and + * for the square matrix B as `input2`. + * + * @receiver the `input1` and the `input2`. + * @return the square matrix x which is the solution of the equation. + */ + public fun solve(a: MutableStructure2D, b: MutableStructure2D): MutableStructure2D } 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 7a076b8b2..0a293783c 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 @@ -17,10 +17,7 @@ import space.kscience.kmath.tensors.api.AnalyticTensorAlgebra import space.kscience.kmath.tensors.api.LinearOpsTensorAlgebra import space.kscience.kmath.tensors.api.Tensor import space.kscience.kmath.tensors.core.internal.* -import kotlin.math.abs -import kotlin.math.ceil -import kotlin.math.floor -import kotlin.math.sqrt +import kotlin.math.* /** * Implementation of basic operations over double tensors and basic algebra operations on them. @@ -706,14 +703,18 @@ public open class DoubleTensorAlgebra : override fun svd( structureND: StructureND, ): Triple, StructureND, StructureND> = - svd(structureND = structureND, epsilon = 1e-10) + svdGolubKahan(structureND = structureND, epsilon = 1e-10) override fun symEig(structureND: StructureND): Pair = symEigJacobi(structureND = structureND, maxIteration = 50, epsilon = 1e-15) + override fun solve(a: MutableStructure2D, b: MutableStructure2D): MutableStructure2D { + val aSvd = DoubleTensorAlgebra.svd(a) + val s = BroadcastDoubleTensorAlgebra.diagonalEmbedding(aSvd.second.map {1.0 / it}) + val aInverse = aSvd.third.dot(s).dot(aSvd.first.transposed()) + return aInverse.dot(b).as2D() + } } public val Double.Companion.tensorAlgebra: DoubleTensorAlgebra get() = DoubleTensorAlgebra -public val DoubleField.tensorAlgebra: DoubleTensorAlgebra get() = DoubleTensorAlgebra - - +public val DoubleField.tensorAlgebra: DoubleTensorAlgebra get() = DoubleTensorAlgebra \ No newline at end of file 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..3cb485d7d --- /dev/null +++ b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/LevenbergMarquardtAlgorithm.kt @@ -0,0 +1,589 @@ +/* + * 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.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 + +/** + * Type of convergence achieved as a result of executing the Levenberg-Marquardt algorithm. + * + * InGradient: gradient convergence achieved + * (max(J^T W dy) < epsilon1, + * where J - Jacobi matrix (dy^/dp) for the current approximation y^, + * W - weight matrix from input, dy = (y - y^(p))). + * InParameters: convergence in parameters achieved + * (max(h_i / p_i) < epsilon2, + * where h_i - offset for parameter p_i on the current iteration). + * InReducedChiSquare: chi-squared convergence achieved + * (chi squared value divided by (m - n + 1) < epsilon2, + * where n - number of parameters, m - amount of points). + * NoConvergence: the maximum number of iterations has been reached without reaching any convergence. + */ +public enum class TypeOfConvergence { + InGradient, + InParameters, + InReducedChiSquare, + NoConvergence +} + +/** + * 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 funcCalls: Int, + var resultChiSq: Double, + var resultLambda: Double, + var resultParameters: MutableStructure2D, + var typeOfConvergence: TypeOfConvergence, +) + +/** + * Input data for the Levenberg-Marquardt function. + * + * func: function of n independent variables x, m parameters an example number, + * rotating a vector of n values y, in which each of the y_i is calculated at its x_i with the given parameters. + * startParameters: starting parameters. + * independentVariables: independent variables, for each of which the real value is known. + * realValues: real values obtained with given independent variables but unknown parameters. + * weight: measurement error for realValues (denominator in each term of sum of weighted squared errors). + * pDelta: delta when calculating the derivative with respect to parameters. + * minParameters: the lower bound of parameter values. + * maxParameters: upper limit of parameter values. + * maxIterations: maximum allowable number of iterations. + * epsilons: epsilon1 - convergence tolerance for gradient, + * epsilon2 - convergence tolerance for parameters, + * epsilon3 - convergence tolerance for reduced chi-square, + * epsilon4 - determines acceptance of a step. + * lambdas: lambda0 - starting lambda value for parameter offset count, + * lambdaUp - factor for increasing lambda, + * lambdaDown - factor for decreasing lambda. + * updateType: 1: Levenberg-Marquardt lambda update, + * 2: Quadratic update, + * 3: Nielsen's lambda update equations. + * nargin: a value that determines which options to use by default + * (<5 - use weight by default, <6 - use pDelta by default, <7 - use minParameters by default, + * <8 - use maxParameters by default, <9 - use updateType by default). + * exampleNumber: a parameter for a function with which you can choose its behavior. + */ +public data class LMInput ( + var func: (MutableStructure2D, MutableStructure2D, Int) -> (MutableStructure2D), + var startParameters: MutableStructure2D, + var independentVariables: MutableStructure2D, + var realValues: MutableStructure2D, + var weight: Double, + var pDelta: MutableStructure2D, + var minParameters: MutableStructure2D, + var maxParameters: MutableStructure2D, + var maxIterations: Int, + var epsilons: DoubleArray, + var lambdas: DoubleArray, + var updateType: Int, + var nargin: Int, + var exampleNumber: Int +) + + +/** + * Levenberg-Marquardt optimization. + * + * An optimization method that iteratively searches for the optimal function parameters + * that best describe the dataset. The 'input' is the function being optimized, a set of real data + * (calculated with independent variables, but with an unknown set of parameters), a set of + * independent variables, and variables for adjusting the algorithm, described in the documentation for the LMInput class. + * The function returns number of completed iterations, the number of evaluations of the input function during execution, + * chi squared value on final parameters, final lambda parameter used to calculate the offset, final parameters + * and type of convergence in the 'output'. + * + * @receiver the `input`. + * @return the 'output'. + */ +public fun DoubleTensorAlgebra.levenbergMarquardt(inputData: LMInput): LMResultInfo { + val resultInfo = LMResultInfo(0, 0, 0.0, + 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 + + var p = inputData.startParameters + val t = inputData.independentVariables + + val Npar = length(p) // number of parameters + val Npnt = length(inputData.realValues) // number of data points + var pOld = zeros(ShapeND(intArrayOf(Npar, 1))).as2D() // previous set of parameters + var yOld = zeros(ShapeND(intArrayOf(Npnt, 1))).as2D() // previous model, y_old = y_hat(t;p_old) + var X2 = 1e-3 / eps // a really big initial Chi-sq value + var X2Old = 1e-3 / eps // a really big initial Chi-sq value + var J = zeros(ShapeND(intArrayOf(Npnt, Npar))).as2D() // Jacobian matrix + val DoF = Npnt - Npar // statistical degrees of freedom + + var weight = fromArray(ShapeND(intArrayOf(1, 1)), doubleArrayOf(inputData.weight)).as2D() + if (inputData.nargin < 5) { + weight = fromArray(ShapeND(intArrayOf(1, 1)), doubleArrayOf((inputData.realValues.transpose().dot(inputData.realValues)).as1D()[0])).as2D() + } + + var dp = inputData.pDelta + if (inputData.nargin < 6) { + dp = fromArray(ShapeND(intArrayOf(1, 1)), doubleArrayOf(0.001)).as2D() + } + + var minParameters = inputData.minParameters + if (inputData.nargin < 7) { + minParameters = p + minParameters.abs() + minParameters = minParameters.div(-100.0).as2D() + } + + var maxParameters = inputData.maxParameters + if (inputData.nargin < 8) { + maxParameters = p + maxParameters.abs() + maxParameters = maxParameters.div(100.0).as2D() + } + + var maxIterations = inputData.maxIterations + var epsilon1 = inputData.epsilons[0] + var epsilon2 = inputData.epsilons[1] + var epsilon3 = inputData.epsilons[2] + var epsilon4 = inputData.epsilons[3] + var lambda0 = inputData.lambdas[0] + var lambdaUpFac = inputData.lambdas[1] + var lambdaDnFac = inputData.lambdas[2] + var updateType = inputData.updateType + + if (inputData.nargin < 9) { + maxIterations = 10 * Npar + epsilon1 = 1e-3 + epsilon2 = 1e-3 + epsilon3 = 1e-1 + epsilon4 = 1e-1 + lambda0 = 1e-2 + lambdaUpFac = 11.0 + lambdaDnFac = 9.0 + updateType = 1 + } + + minParameters = makeColumn(minParameters) + maxParameters = makeColumn(maxParameters) + + if (length(makeColumn(dp)) == 1) { + dp = ones(ShapeND(intArrayOf(Npar, 1))).div(1 / dp[0, 0]).as2D() + } + + var stop = false // termination flag + + if (weight.shape.component1() == 1 || variance(weight) == 0.0) { // identical weights vector + weight = ones(ShapeND(intArrayOf(Npnt, 1))).div(1 / kotlin.math.abs(weight[0, 0])).as2D() + } + else { + weight = makeColumn(weight) + weight.abs() + } + + // initialize Jacobian with finite difference calculation + var lmMatxAns = lmMatx(inputData.func, t, pOld, yOld, 1, J, p, inputData.realValues, weight, dp, settings) + var JtWJ = lmMatxAns[0] + var JtWdy = lmMatxAns[1] + X2 = lmMatxAns[2][0, 0] + var yHat = lmMatxAns[3] + J = lmMatxAns[4] + + if ( abs(JtWdy).max() < epsilon1 ) { + stop = true + } + + var lambda = 1.0 + var nu = 1 + + if (updateType == 1) { + lambda = lambda0 // Marquardt: init'l lambda + } + else { + lambda = lambda0 * (makeColumnFromDiagonal(JtWJ)).max() + nu = 2 + } + + X2Old = X2 // previous value of X2 + + var h: DoubleTensor + + while (!stop && settings.iteration <= maxIterations) { + settings.iteration += 1 + + // incremental change in parameters + h = if (updateType == 1) { // Marquardt + val solve = solve(JtWJ.plus(makeMatrixWithDiagonal(makeColumnFromDiagonal(JtWJ)).div(1 / lambda)).as2D(), JtWdy) + solve.asDoubleTensor() + } else { // Quadratic and Nielsen + val solve = solve(JtWJ.plus(lmEye(Npar).div(1 / lambda)).as2D(), JtWdy) + solve.asDoubleTensor() + } + + var pTry = (p + h).as2D() // update the [idx] elements + pTry = smallestElementComparison(largestElementComparison(minParameters, pTry.as2D()), maxParameters) // apply constraints + + var deltaY = inputData.realValues.minus(evaluateFunction(inputData.func, t, pTry, inputData.exampleNumber)) // residual error using p_try + + for (i in 0 until deltaY.shape.component1()) { // floating point error; break + for (j in 0 until deltaY.shape.component2()) { + if (deltaY[i, j] == Double.POSITIVE_INFINITY || deltaY[i, j] == Double.NEGATIVE_INFINITY) { + stop = true + break + } + } + } + + settings.funcCalls += 1 + + val tmp = deltaY.times(weight) + var X2Try = deltaY.as2D().transpose().dot(tmp) // Chi-squared error criteria + + val alpha = 1.0 + if (updateType == 2) { // Quadratic + // One step of quadratic line update in the h direction for minimum X2 + val alpha = JtWdy.transpose().dot(h) / ((X2Try.minus(X2)).div(2.0).plus(2 * JtWdy.transpose().dot(h))) + h = h.dot(alpha) + pTry = p.plus(h).as2D() // update only [idx] elements + pTry = smallestElementComparison(largestElementComparison(minParameters, pTry), maxParameters) // apply constraints + + deltaY = inputData.realValues.minus(evaluateFunction(inputData.func, t, pTry, inputData.exampleNumber)) // residual error using p_try + settings.funcCalls += 1 + + X2Try = deltaY.as2D().transpose().dot(deltaY.times(weight)) // Chi-squared error criteria + } + + val rho = when (updateType) { // Nielsen + 1 -> { + val tmp = h.transposed() + .dot(makeMatrixWithDiagonal(makeColumnFromDiagonal(JtWJ)).div(1 / lambda).dot(h).plus(JtWdy)) + X2.minus(X2Try).as2D()[0, 0] / abs(tmp.as2D()).as2D()[0, 0] + } + else -> { + val tmp = h.transposed().dot(h.div(1 / lambda).plus(JtWdy)) + X2.minus(X2Try).as2D()[0, 0] / abs(tmp.as2D()).as2D()[0, 0] + } + } + + if (rho > epsilon4) { // it IS significantly better + val dX2 = X2.minus(X2Old) + X2Old = X2 + pOld = p.copyToTensor().as2D() + yOld = yHat.copyToTensor().as2D() + p = makeColumn(pTry) // accept p_try + + lmMatxAns = lmMatx(inputData.func, t, pOld, yOld, dX2.toInt(), J, p, inputData.realValues, weight, dp, settings) + // decrease lambda ==> Gauss-Newton method + JtWJ = lmMatxAns[0] + JtWdy = lmMatxAns[1] + X2 = lmMatxAns[2][0, 0] + yHat = lmMatxAns[3] + J = lmMatxAns[4] + + lambda = when (updateType) { + 1 -> { // Levenberg + max(lambda / lambdaDnFac, 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 = X2Old // do not accept p_try + if (settings.iteration % (2 * Npar) == 0) { // rank-1 update of Jacobian + lmMatxAns = lmMatx(inputData.func, t, pOld, yOld, -1, J, p, inputData.realValues, weight, dp, settings) + JtWJ = lmMatxAns[0] + JtWdy = lmMatxAns[1] + yHat = lmMatxAns[3] + J = lmMatxAns[4] + } + + // increase lambda ==> gradient descent method + lambda = when (updateType) { + 1 -> { // Levenberg + min(lambda * lambdaUpFac, 1e7) + } + + 2 -> { // Quadratic + lambda + kotlin.math.abs(((X2Try.as2D()[0, 0] - X2) / 2) / alpha) + } + + else -> { // Nielsen + nu *= 2 + lambda * (nu / 2) + } + } + } + + val chiSq = X2 / DoF + resultInfo.iterations = settings.iteration + resultInfo.funcCalls = settings.funcCalls + resultInfo.resultChiSq = chiSq + resultInfo.resultLambda = lambda + resultInfo.resultParameters = p + + + if (abs(JtWdy).max() < epsilon1 && settings.iteration > 2) { + resultInfo.typeOfConvergence = TypeOfConvergence.InGradient + stop = true + } + if ((abs(h.as2D()).div(abs(p) + 1e-12)).max() < epsilon2 && settings.iteration > 2) { + resultInfo.typeOfConvergence = TypeOfConvergence.InParameters + stop = true + } + if (X2 / DoF < epsilon3 && settings.iteration > 2) { + resultInfo.typeOfConvergence = TypeOfConvergence.InReducedChiSquare + stop = true + } + if (settings.iteration == maxIterations) { + resultInfo.typeOfConvergence = TypeOfConvergence.NoConvergence + stop = true + } + } + return resultInfo +} + +private data class LMSettings ( + var iteration:Int, + var funcCalls: Int, + var exampleNumber:Int +) + +/* matrix -> column of all elements */ +private fun makeColumn(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] + } + } + return BroadcastDoubleTensorAlgebra.fromArray(ShapeND(shape), buffer).as2D() +} + +/* column length */ +private fun length(column: MutableStructure2D) : Int { + return column.shape.component1() +} + +private fun MutableStructure2D.abs() { + for (i in 0 until this.shape.component1()) { + for (j in 0 until this.shape.component2()) { + this[i, j] = kotlin.math.abs(this[i, j]) + } + } +} + +private 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 +} + +private fun makeColumnFromDiagonal(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 +} + +private fun makeMatrixWithDiagonal(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 +} + +private fun lmEye(size: Int): MutableStructure2D { + val column = BroadcastDoubleTensorAlgebra.ones(ShapeND(intArrayOf(size, 1))).as2D() + return makeMatrixWithDiagonal(column) +} + +private fun largestElementComparison(a: MutableStructure2D, b: MutableStructure2D): MutableStructure2D { + val aSizeX = a.shape.component1() + val aSizeY = a.shape.component2() + val bSizeX = b.shape.component1() + val bSizeY = b.shape.component2() + val tensor = BroadcastDoubleTensorAlgebra.zeros(ShapeND(intArrayOf(max(aSizeX, bSizeX), max(aSizeY, bSizeY)))).as2D() + for (i in 0 until tensor.shape.component1()) { + for (j in 0 until tensor.shape.component2()) { + if (i < aSizeX && i < bSizeX && j < aSizeY && j < bSizeY) { + tensor[i, j] = max(a[i, j], b[i, j]) + } + else if (i < aSizeX && j < aSizeY) { + tensor[i, j] = a[i, j] + } + else { + tensor[i, j] = b[i, j] + } + } + } + return tensor +} + +private fun smallestElementComparison(a: MutableStructure2D, b: MutableStructure2D): MutableStructure2D { + val aSizeX = a.shape.component1() + val aSizeY = a.shape.component2() + val bSizeX = b.shape.component1() + val bSizeY = b.shape.component2() + val tensor = BroadcastDoubleTensorAlgebra.zeros(ShapeND(intArrayOf(max(aSizeX, bSizeX), max(aSizeY, bSizeY)))).as2D() + for (i in 0 until tensor.shape.component1()) { + for (j in 0 until tensor.shape.component2()) { + if (i < aSizeX && i < bSizeX && j < aSizeY && j < bSizeY) { + tensor[i, j] = min(a[i, j], b[i, j]) + } + else if (i < aSizeX && j < aSizeY) { + tensor[i, j] = a[i, j] + } + else { + tensor[i, j] = b[i, j] + } + } + } + return tensor +} + +private fun getZeroIndices(column: MutableStructure2D, epsilon: Double = 0.000001): MutableStructure2D? { + var idx = emptyArray() + for (i in 0 until column.shape.component1()) { + if (kotlin.math.abs(column[i, 0]) > epsilon) { + idx += (i + 1.0) + } + } + if (idx.isNotEmpty()) { + return BroadcastDoubleTensorAlgebra.fromArray(ShapeND(intArrayOf(idx.size, 1)), idx.toDoubleArray()).as2D() + } + return null +} + +private fun evaluateFunction(func: (MutableStructure2D, MutableStructure2D, Int) -> MutableStructure2D, + t: MutableStructure2D, p: MutableStructure2D, exampleNumber: Int) + : MutableStructure2D +{ + return func(t, p, exampleNumber) +} + +private fun lmMatx(func: (MutableStructure2D, MutableStructure2D, Int) -> MutableStructure2D, + t: MutableStructure2D, pOld: MutableStructure2D, yOld: MutableStructure2D, + dX2: Int, JInput: MutableStructure2D, p: MutableStructure2D, + yDat: MutableStructure2D, weight: MutableStructure2D, dp:MutableStructure2D, settings:LMSettings) : Array> +{ + // default: dp = 0.001 + val Npar = length(p) // number of parameters + + val yHat = evaluateFunction(func, t, p, settings.exampleNumber) // evaluate model using parameters 'p' + settings.funcCalls += 1 + + var J = JInput + + J = if (settings.iteration % (2 * Npar) == 0 || dX2 > 0) { + lmFdJ(func, t, p, yHat, dp, settings).as2D() // finite difference + } + else { + lmBroydenJ(pOld, yOld, J, p, yHat).as2D() // rank-1 update + } + + val deltaY = yDat.minus(yHat) + + val chiSq = deltaY.transposed().dot( deltaY.times(weight) ).as2D() + val JtWJ = J.transposed().dot ( J.times( weight.dot(BroadcastDoubleTensorAlgebra.ones(ShapeND(intArrayOf(1, Npar)))) ) ).as2D() + val JtWdy = J.transposed().dot( weight.times(deltaY) ).as2D() + + return arrayOf(JtWJ,JtWdy,chiSq,yHat,J) +} + +private fun lmBroydenJ(pOld: MutableStructure2D, yOld: MutableStructure2D, JInput: MutableStructure2D, + p: MutableStructure2D, y: MutableStructure2D): MutableStructure2D { + var J = JInput.copyToTensor() + + val h = p.minus(pOld) + val increase = y.minus(yOld).minus( J.dot(h) ).dot(h.transposed()).div( (h.transposed().dot(h)).as2D()[0, 0] ) + J = J.plus(increase) + + return J.as2D() +} + +private fun lmFdJ(func: (MutableStructure2D, MutableStructure2D, exampleNumber: Int) -> MutableStructure2D, + t: MutableStructure2D, p: MutableStructure2D, y: MutableStructure2D, + dp: MutableStructure2D, settings: LMSettings): MutableStructure2D { + // 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 = evaluateFunction(func, t, p, settings.exampleNumber) + settings.funcCalls += 1 + + if (dp[j, 0] < 0) { // backwards difference + for (i in 0 until J.shape.component1()) { + J[i, j] = (y1.as2D().minus(y).as2D())[i, 0] / del[j, 0] + } + } + else { + // Do tests for it + p[j, 0] = ps[j, 0] - del[j, 0] // central difference, additional func call + for (i in 0 until J.shape.component1()) { + J[i, j] = (y1.as2D().minus(evaluateFunction(func, t, p, settings.exampleNumber)).as2D())[i, 0] / (2 * del[j, 0]) + } + settings.funcCalls += 1 + } + } + + p[j, 0] = ps[j, 0] + } + + 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 cf3697e76..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 @@ -7,12 +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.DoubleBuffer -import space.kscience.kmath.structures.IntBuffer -import space.kscience.kmath.structures.asBuffer -import space.kscience.kmath.structures.indices +import space.kscience.kmath.structures.* import space.kscience.kmath.tensors.core.* import kotlin.math.abs +import kotlin.math.max import kotlin.math.min import kotlin.math.sqrt @@ -308,3 +306,299 @@ internal fun DoubleTensorAlgebra.svdHelper( matrixV.source[i] = vBuffer[i] } } + +private fun pythag(a: Double, b: Double): Double { + val at: Double = abs(a) + val bt: Double = abs(b) + val ct: Double + val result: Double + if (at > bt) { + ct = bt / at + result = at * sqrt(1.0 + ct * ct) + } else if (bt > 0.0) { + ct = at / bt + result = bt * sqrt(1.0 + ct * ct) + } else result = 0.0 + return result +} + +private fun SIGN(a: Double, b: Double): Double { + if (b >= 0.0) + return abs(a) + else + return -abs(a) +} + +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 + val rv1 = DoubleArray(n) + var s = 0.0 + var scale = 0.0 + var anorm = 0.0 + var g = 0.0 + var l = 0 + + val wStart = 0 + val wBuffer = w.source + + for (i in 0 until n) { + /* left-hand reduction */ + l = i + 1 + rv1[i] = scale * g + g = 0.0 + s = 0.0 + scale = 0.0 + if (i < m) { + for (k in i until m) { + scale += abs(this[k, i]); + } + if (abs(scale) > epsilon) { + for (k in i until m) { + this[k, i] = (this[k, i] / scale) + s += this[k, i] * this[k, i] + } + f = this[i, i] + if (f >= 0) { + g = (-1) * abs(sqrt(s)) + } else { + g = abs(sqrt(s)) + } + val h = f * g - s + this[i, i] = f - g + if (i != n - 1) { + for (j in l until n) { + s = 0.0 + for (k in i until m) { + s += this[k, i] * this[k, j] + } + f = s / h + for (k in i until m) { + this[k, j] += f * this[k, i] + } + } + } + for (k in i until m) { + this[k, i] = this[k, i] * scale + } + } + } + + wBuffer[wStart + i] = scale * g + /* right-hand reduction */ + g = 0.0 + s = 0.0 + scale = 0.0 + if (i < m && i != n - 1) { + for (k in l until n) { + scale += abs(this[i, k]) + } + if (abs(scale) > epsilon) { + for (k in l until n) { + this[i, k] = this[i, k] / scale + s += this[i, k] * this[i, k] + } + f = this[i, l] + if (f >= 0) { + g = (-1) * abs(sqrt(s)) + } else { + g = abs(sqrt(s)) + } + val h = f * g - s + this[i, l] = f - g + for (k in l until n) { + rv1[k] = this[i, k] / h + } + if (i != m - 1) { + for (j in l until m) { + s = 0.0 + for (k in l until n) { + s += this[j, k] * this[i, k] + } + for (k in l until n) { + this[j, k] += s * rv1[k] + } + } + } + for (k in l until n) { + this[i, k] = this[i, k] * scale + } + } + } + anorm = max(anorm, (abs(wBuffer[wStart + i]) + abs(rv1[i]))); + } + + for (i in n - 1 downTo 0) { + if (i < n - 1) { + if (abs(g) > epsilon) { + for (j in l until n) { + v[j, i] = (this[i, j] / this[i, l]) / g + } + for (j in l until n) { + s = 0.0 + for (k in l until n) + s += this[i, k] * v[k, j] + for (k in l until n) + v[k, j] += s * v[k, i] + } + } + for (j in l until n) { + v[i, j] = 0.0 + v[j, i] = 0.0 + } + } + v[i, i] = 1.0 + g = rv1[i] + l = i + } + + for (i in min(n, m) - 1 downTo 0) { + l = i + 1 + g = wBuffer[wStart + i] + for (j in l until n) { + this[i, j] = 0.0 + } + if (abs(g) > epsilon) { + g = 1.0 / g + for (j in l until n) { + s = 0.0 + for (k in l until m) { + s += this[k, i] * this[k, j] + } + f = (s / this[i, i]) * g + for (k in i until m) { + this[k, j] += f * this[k, i] + } + } + for (j in i until m) { + this[j, i] *= g + } + } else { + for (j in i until m) { + this[j, i] = 0.0 + } + } + this[i, i] += 1.0 + } + + var flag = 0 + var nm = 0 + var c = 0.0 + var h = 0.0 + var y = 0.0 + var z = 0.0 + var x = 0.0 + for (k in n - 1 downTo 0) { + for (its in 1 until iterations) { + flag = 1 + for (newl in k downTo 0) { + nm = newl - 1 + if (abs(rv1[newl]) + anorm == anorm) { + flag = 0 + l = newl + break + } + if (abs(wBuffer[wStart + nm]) + anorm == anorm) { + l = newl + break + } + } + + if (flag != 0) { + c = 0.0 + s = 1.0 + for (i in l until k + 1) { + f = s * rv1[i] + rv1[i] = c * rv1[i] + if (abs(f) + anorm == anorm) { + break + } + g = wBuffer[wStart + i] + h = pythag(f, g) + wBuffer[wStart + i] = h + h = 1.0 / h + c = g * h + s = (-f) * h + for (j in 0 until m) { + y = this[j, nm] + z = this[j, i] + this[j, nm] = y * c + z * s + this[j, i] = z * c - y * s + } + } + } + + z = wBuffer[wStart + k] + if (l == k) { + if (z < 0.0) { + wBuffer[wStart + k] = -z + for (j in 0 until n) + v[j, k] = -v[j, k] + } + break + } + + x = wBuffer[wStart + l] + nm = k - 1 + y = wBuffer[wStart + nm] + g = rv1[nm] + h = rv1[k] + f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y) + g = pythag(f, 1.0) + f = ((x - z) * (x + z) + h * ((y / (f + SIGN(g, f))) - h)) / x + c = 1.0 + s = 1.0 + + var i = 0 + for (j in l until nm + 1) { + i = j + 1 + g = rv1[i] + y = wBuffer[wStart + i] + h = s * g + g = c * g + z = pythag(f, h) + rv1[j] = z + c = f / z + s = h / z + f = x * c + g * s + g = g * c - x * s + h = y * s + y *= c + + for (jj in 0 until n) { + x = v[jj, j]; + z = v[jj, i]; + v[jj, j] = x * c + z * s; + v[jj, i] = z * c - x * s; + } + z = pythag(f, h) + wBuffer[wStart + j] = z + if (abs(z) > epsilon) { + z = 1.0 / z + c = f * z + s = h * z + } + f = c * g + s * y + x = c * y - s * g + for (jj in 0 until m) { + y = this[jj, j] + z = this[jj, i] + this[jj, j] = y * c + z * s + this[jj, i] = z * c - y * s + } + } + rv1[l] = 0.0 + rv1[k] = f + wBuffer[wStart + k] = x + } + } + + for (i in 0 until n) { + for (j in 0 until m) { + u[j, i] = this[j, i] + } + } +} \ No newline at end of file diff --git a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/tensorOps.kt b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/tensorOps.kt index e5dc55f68..88ffb0bfe 100644 --- a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/tensorOps.kt +++ b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/tensorOps.kt @@ -212,6 +212,36 @@ public fun DoubleTensorAlgebra.svd( return Triple(uTensor.transposed(), sTensor, vTensor.transposed()) } +public fun DoubleTensorAlgebra.svdGolubKahan( + structureND: StructureND, + iterations: Int = 30, epsilon: Double = 1e-10 +): Triple { + val size = structureND.dimension + val commonShape = structureND.shape.slice(0 until size - 2) + val (n, m) = structureND.shape.slice(size - 2 until size) + val uTensor = zeros(commonShape + intArrayOf(n, m)) + val sTensor = zeros(commonShape + intArrayOf(m)) + val vTensor = zeros(commonShape + intArrayOf(m, m)) + + val matrices = structureND.asDoubleTensor().matrices + val uTensors = uTensor.matrices + val sTensorVectors = sTensor.vectors + val vTensors = vTensor.matrices + + for (index in matrices.indices) { + val matrix = matrices[index] + val matrixSize = matrix.shape.linearSize + val curMatrix = DoubleTensor( + matrix.shape, + matrix.source.view(0, matrixSize).copy() + ) + curMatrix.as2D().svdGolubKahanHelper(uTensors[index].as2D(), sTensorVectors[index], vTensors[index].as2D(), + iterations, epsilon) + } + + return Triple(uTensor, sTensor, vTensor) +} + /** * Returns eigenvalues and eigenvectors of a real symmetric matrix input or a batch of real symmetric matrices, * represented by a pair `eigenvalues to eigenvectors`. 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 cae01bed8..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 @@ -6,9 +6,7 @@ package space.kscience.kmath.tensors.core -import space.kscience.kmath.nd.ShapeND -import space.kscience.kmath.nd.contentEquals -import space.kscience.kmath.nd.get +import space.kscience.kmath.nd.* import space.kscience.kmath.operations.invoke import space.kscience.kmath.testutils.assertBufferEquals import kotlin.test.Test 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 new file mode 100644 index 000000000..4b031eb11 --- /dev/null +++ b/kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestLmAlgorithm.kt @@ -0,0 +1,280 @@ +/* + * 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.nd.MutableStructure2D +import space.kscience.kmath.nd.ShapeND +import space.kscience.kmath.nd.as2D +import space.kscience.kmath.nd.component1 +import space.kscience.kmath.operations.invoke +import space.kscience.kmath.tensors.core.DoubleTensorAlgebra.Companion.max +import space.kscience.kmath.tensors.core.DoubleTensorAlgebra.Companion.plus +import space.kscience.kmath.tensors.core.DoubleTensorAlgebra.Companion.pow +import space.kscience.kmath.tensors.core.DoubleTensorAlgebra.Companion.times +import kotlin.math.roundToLong +import kotlin.test.Test +import kotlin.test.assertEquals + +class TestLmAlgorithm { + companion object { + fun funcEasyForLm(t: MutableStructure2D, p: MutableStructure2D, exampleNumber: Int): MutableStructure2D { + val m = t.shape.component1() + var yHat = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf(m, 1))) + + if (exampleNumber == 1) { + yHat = DoubleTensorAlgebra.exp((t.times(-1.0 / p[1, 0]))).times(p[0, 0]) + t.times(p[2, 0]).times( + DoubleTensorAlgebra.exp((t.times(-1.0 / p[3, 0]))) + ) + } + else if (exampleNumber == 2) { + val mt = t.max() + yHat = (t.times(1.0 / mt)).times(p[0, 0]) + + (t.times(1.0 / mt)).pow(2).times(p[1, 0]) + + (t.times(1.0 / mt)).pow(3).times(p[2, 0]) + + (t.times(1.0 / mt)).pow(4).times(p[3, 0]) + } + else if (exampleNumber == 3) { + yHat = DoubleTensorAlgebra.exp((t.times(-1.0 / p[1, 0]))) + .times(p[0, 0]) + DoubleTensorAlgebra.sin((t.times(1.0 / p[3, 0]))).times(p[2, 0]) + } + + return yHat.as2D() + } + + fun funcMiddleForLm(t: MutableStructure2D, p: MutableStructure2D, exampleNumber: Int): MutableStructure2D { + val m = t.shape.component1() + var yHat = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf (m, 1))) + + val mt = t.max() + for(i in 0 until p.shape.component1()){ + yHat += (t.times(1.0 / mt)).times(p[i, 0]) + } + + for(i in 0 until 5){ + yHat = funcEasyForLm(yHat.as2D(), p, exampleNumber).asDoubleTensor() + } + + return yHat.as2D() + } + + fun funcDifficultForLm(t: MutableStructure2D, p: MutableStructure2D, exampleNumber: Int): MutableStructure2D { + val m = t.shape.component1() + var yHat = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf (m, 1))) + + val mt = t.max() + for(i in 0 until p.shape.component1()){ + yHat = yHat.plus( (t.times(1.0 / mt)).times(p[i, 0]) ) + } + + for(i in 0 until 4){ + yHat = funcEasyForLm((yHat.as2D() + t).as2D(), p, exampleNumber).asDoubleTensor() + } + + return yHat.as2D() + } + } + @Test + fun testLMEasy() = DoubleTensorAlgebra { + val lmMatxYDat = 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 exampleNumber = 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 yDat = BroadcastDoubleTensorAlgebra.fromArray( + ShapeND(intArrayOf(100, 1)), lmMatxYDat + ).as2D() + + val weight = 4.0 + + val dp = BroadcastDoubleTensorAlgebra.fromArray( + ShapeND(intArrayOf(1, 1)), DoubleArray(1) { -0.01 } + ).as2D() + + val pMin = BroadcastDoubleTensorAlgebra.fromArray( + ShapeND(intArrayOf(4, 1)), doubleArrayOf(-50.0, -20.0, -2.0, -100.0) + ).as2D() + + val pMax = BroadcastDoubleTensorAlgebra.fromArray( + ShapeND(intArrayOf(4, 1)), doubleArrayOf(50.0, 20.0, 2.0, 100.0) + ).as2D() + + val inputData = LMInput(::funcEasyForLm, p_init, t, yDat, weight, dp, pMin, pMax, 100, + doubleArrayOf(1e-3, 1e-3, 1e-1, 1e-1), doubleArrayOf(1e-2, 11.0, 9.0), 1, 10, exampleNumber) + + val result = levenbergMarquardt(inputData) + assertEquals(13, result.iterations) + assertEquals(31, result.funcCalls) + assertEquals(0.9131368192633, (result.resultChiSq * 1e13).roundToLong() / 1e13) + assertEquals(3.7790980 * 1e-7, (result.resultLambda * 1e13).roundToLong() / 1e13) + assertEquals(result.typeOfConvergence, TypeOfConvergence.InParameters) + val expectedParameters = BroadcastDoubleTensorAlgebra.fromArray( + ShapeND(intArrayOf(4, 1)), doubleArrayOf(20.527230909086, 9.833627103230, 0.997571256572, 50.174445822506) + ).as2D() + result.resultParameters = result.resultParameters.map { x -> (x * 1e12).toLong() / 1e12}.as2D() + val receivedParameters = BroadcastDoubleTensorAlgebra.fromArray( + ShapeND(intArrayOf(4, 1)), doubleArrayOf(result.resultParameters[0, 0], result.resultParameters[1, 0], + result.resultParameters[2, 0], result.resultParameters[3, 0]) + ).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]) + } + + @Test + fun TestLMMiddle() = DoubleTensorAlgebra { + val NData = 100 + val tExample = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(NData, 1))).as2D() + for (i in 0 until NData) { + tExample[i, 0] = tExample[i, 0] * (i + 1) + } + + val Nparams = 20 + val pExample = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1))).as2D() + for (i in 0 until Nparams) { + pExample[i, 0] = pExample[i, 0] + i - 25 + } + + val exampleNumber = 1 + + val yHat = funcMiddleForLm(tExample, pExample, exampleNumber) + + val pInit = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf(Nparams, 1))).as2D() + for (i in 0 until Nparams) { + pInit[i, 0] = (pExample[i, 0] + 0.9) + } + + val t = tExample + val yDat = yHat + val weight = 1.0 + val dp = BroadcastDoubleTensorAlgebra.fromArray( + ShapeND(intArrayOf(1, 1)), DoubleArray(1) { -0.01 } + ).as2D() + var pMin = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1))) + pMin = pMin.div(1.0 / -50.0) + val pMax = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1))) + pMin = pMin.div(1.0 / 50.0) + val opts = doubleArrayOf(3.0, 7000.0, 1e-5, 1e-5, 1e-5, 1e-5, 1e-5, 11.0, 9.0, 1.0) + + val inputData = LMInput(::funcMiddleForLm, + pInit.as2D(), + t, + yDat, + weight, + dp, + pMin.as2D(), + pMax.as2D(), + opts[1].toInt(), + doubleArrayOf(opts[2], opts[3], opts[4], opts[5]), + doubleArrayOf(opts[6], opts[7], opts[8]), + opts[9].toInt(), + 10, + 1) + + val result = DoubleTensorAlgebra.levenbergMarquardt(inputData) + + assertEquals(46, result.iterations) + assertEquals(113, result.funcCalls) + assertEquals(0.000005977, (result.resultChiSq * 1e9).roundToLong() / 1e9) + assertEquals(1.0 * 1e-7, (result.resultLambda * 1e13).roundToLong() / 1e13) + assertEquals(result.typeOfConvergence, TypeOfConvergence.InReducedChiSquare) + val expectedParameters = BroadcastDoubleTensorAlgebra.fromArray( + ShapeND(intArrayOf(Nparams, 1)), doubleArrayOf( -23.9717, -18.6686, -21.7971, + -20.9681, -22.086, -20.5859, -19.0384, -17.4957, -15.9991, -14.576, -13.2441, - + 12.0201, -10.9256, -9.9878, -9.2309, -8.6589, -8.2365, -7.8783, -7.4598, -6.8511)).as2D() + result.resultParameters = result.resultParameters.map { x -> (x * 1e4).roundToLong() / 1e4}.as2D() + val receivedParameters = zeros(ShapeND(intArrayOf(Nparams, 1))).as2D() + for (i in 0 until Nparams) { + receivedParameters[i, 0] = result.resultParameters[i, 0] + assertEquals(expectedParameters[i, 0], result.resultParameters[i, 0]) + } + } + + @Test + fun TestLMDifficult() = DoubleTensorAlgebra { + val NData = 200 + var tExample = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(NData, 1))).as2D() + for (i in 0 until NData) { + tExample[i, 0] = tExample[i, 0] * (i + 1) - 104 + } + + val Nparams = 15 + var pExample = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1))).as2D() + for (i in 0 until Nparams) { + pExample[i, 0] = pExample[i, 0] + i - 25 + } + + val exampleNumber = 1 + + var yHat = funcDifficultForLm(tExample, pExample, exampleNumber) + + var pInit = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf(Nparams, 1))).as2D() + for (i in 0 until Nparams) { + pInit[i, 0] = (pExample[i, 0] + 0.9) + } + + var t = tExample + val yDat = yHat + val weight = 1.0 / Nparams * 1.0 - 0.085 + val dp = BroadcastDoubleTensorAlgebra.fromArray( + ShapeND(intArrayOf(1, 1)), DoubleArray(1) { -0.01 } + ).as2D() + var pMin = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1))) + pMin = pMin.div(1.0 / -50.0) + val pMax = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1))) + pMin = pMin.div(1.0 / 50.0) + val opts = doubleArrayOf(3.0, 7000.0, 1e-2, 1e-3, 1e-2, 1e-2, 1e-2, 11.0, 9.0, 1.0) + + val inputData = LMInput(::funcDifficultForLm, + pInit.as2D(), + t, + yDat, + weight, + dp, + pMin.as2D(), + pMax.as2D(), + opts[1].toInt(), + doubleArrayOf(opts[2], opts[3], opts[4], opts[5]), + doubleArrayOf(opts[6], opts[7], opts[8]), + opts[9].toInt(), + 10, + 1) + + val result = DoubleTensorAlgebra.levenbergMarquardt(inputData) + + assertEquals(2375, result.iterations) + assertEquals(4858, result.funcCalls) + assertEquals(5.14347, (result.resultLambda * 1e5).roundToLong() / 1e5) + assertEquals(result.typeOfConvergence, TypeOfConvergence.InParameters) + val expectedParameters = BroadcastDoubleTensorAlgebra.fromArray( + ShapeND(intArrayOf(Nparams, 1)), doubleArrayOf(-23.6412, -16.7402, -21.5705, -21.0464, + -17.2852, -17.2959, -17.298, 0.9999, -17.2885, -17.3008, -17.2941, -17.2923, -17.2976, -17.3028, -17.2891)).as2D() + result.resultParameters = result.resultParameters.map { x -> (x * 1e4).roundToLong() / 1e4}.as2D() + val receivedParameters = zeros(ShapeND(intArrayOf(Nparams, 1))).as2D() + for (i in 0 until Nparams) { + receivedParameters[i, 0] = result.resultParameters[i, 0] + assertEquals(expectedParameters[i, 0], result.resultParameters[i, 0]) + } + } +} \ No newline at end of file