teldufalsari/kmath
2
0
Fork 0
forked from kscience/kmath
Code Pull Requests 1 Activity

Merge remote-tracking branch 'space/dev' into dev

Browse Source
This commit is contained in:
Alexander Nozik 2023-07-18 11:15:27 +03:00
commit 14f0fa1a6f
16 changed files with 1781 additions and 18 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

90
examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StaticLm/staticDifficultTest.kt Normal file
View File

@ -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)
}

57
examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StaticLm/staticEasyTest.kt Normal file
View File

@ -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)
}

88
examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StaticLm/staticMiddleTest.kt Normal file
View File

@ -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)
}

68
examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StreamingLm/streamLm.kt Normal file
View File

@ -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<Double>, MutableStructure2D<Double>, Int) -> (MutableStructure2D<Double>),
startData: StartDataLm, launchFrequencyInMs: Long, numberOfLaunches: Int): Flow<MutableStructure2D<Double>> = 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<Double>, delta: Double): MutableStructure2D<Double>{
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
}

32
examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StreamingLm/streamingLmTest.kt Normal file
View File

@ -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<Long>()
// Запуск потока
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()}")
}

222
examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/functionsToOptimize.kt Normal file
View File

@ -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<Double>,
var example_number: Int,
var p_init: MutableStructure2D<Double>,
var t: MutableStructure2D<Double>,
var y_dat: MutableStructure2D<Double>,
var weight: Double,
var dp: MutableStructure2D<Double>,
var p_min: MutableStructure2D<Double>,
var p_max: MutableStructure2D<Double>,
var consts: MutableStructure2D<Double>,
var opts: DoubleArray
)
fun funcEasyForLm(t: MutableStructure2D<Double>, p: MutableStructure2D<Double>, exampleNumber: Int): MutableStructure2D<Double> {
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<Double>, p: MutableStructure2D<Double>, exampleNumber: Int): MutableStructure2D<Double> {
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<Double>, p: MutableStructure2D<Double>, exampleNumber: Int): MutableStructure2D<Double> {
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)
}

2
kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/_generated.kt
View File

@ -19,9 +19,9 @@ import org.ejml.sparse.csc.factory.DecompositionFactory_DSCC
import org.ejml.sparse.csc.factory.DecompositionFactory_FSCC
import org.ejml.sparse.csc.factory.LinearSolverFactory_DSCC
import org.ejml.sparse.csc.factory.LinearSolverFactory_FSCC
import space.kscience.kmath.UnstableKMathAPI
import space.kscience.kmath.linear.*
import space.kscience.kmath.linear.Matrix
import space.kscience.kmath.UnstableKMathAPI
import space.kscience.kmath.nd.StructureFeature
import space.kscience.kmath.operations.DoubleField
import space.kscience.kmath.operations.FloatField

2
kmath-symja/build.gradle.kts
View File

@ -10,7 +10,7 @@ plugins {
description = "Symja integration module"
dependencies {
api("org.matheclipse:matheclipse-core:2.0.0-SNAPSHOT") {
api("org.matheclipse:matheclipse-core:2.0.0") {
// Incorrect transitive dependencies
exclude("org.apfloat", "apfloat")
exclude("org.hipparchus", "hipparchus-clustering")

8
kmath-tensors/build.gradle.kts
View File

@ -4,7 +4,13 @@ plugins {
kscience{
jvm()
js()
js {
browser {
testTask {
useMocha().timeout = "0"
}
}
}
native()
dependencies {

8
kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/api/LinearOpsTensorAlgebra.kt
View File

@ -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<T, A : Field<T>> : TensorPartialDivision
*/
public fun symEig(structureND: StructureND<T>): Pair<StructureND<T>, StructureND<T>>
/** 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<Double>, b: MutableStructure2D<Double>): MutableStructure2D<Double>
}

15
kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleTensorAlgebra.kt
View File

@ -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<Double>,
): Triple<StructureND<Double>, StructureND<Double>, StructureND<Double>> =
svd(structureND = structureND, epsilon = 1e-10)
svdGolubKahan(structureND = structureND, epsilon = 1e-10)
override fun symEig(structureND: StructureND<Double>): Pair<DoubleTensor, DoubleTensor> =
symEigJacobi(structureND = structureND, maxIteration = 50, epsilon = 1e-15)
override fun solve(a: MutableStructure2D<Double>, b: MutableStructure2D<Double>): MutableStructure2D<Double> {
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

589
kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/LevenbergMarquardtAlgorithm.kt Normal file
View File

@ -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<Double>,
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<Double>, MutableStructure2D<Double>, Int) -> (MutableStructure2D<Double>),
var startParameters: MutableStructure2D<Double>,
var independentVariables: MutableStructure2D<Double>,
var realValues: MutableStructure2D<Double>,
var weight: Double,
var pDelta: MutableStructure2D<Double>,
var minParameters: MutableStructure2D<Double>,
var maxParameters: MutableStructure2D<Double>,
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<Double>): MutableStructure2D<Double> {
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<Double>) : Int {
return column.shape.component1()
}
private fun MutableStructure2D<Double>.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<Double>): MutableStructure2D<Double> {
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<Double>): MutableStructure2D<Double> {
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<Double>): MutableStructure2D<Double> {
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<Double> {
val column = BroadcastDoubleTensorAlgebra.ones(ShapeND(intArrayOf(size, 1))).as2D()
return makeMatrixWithDiagonal(column)
}
private fun largestElementComparison(a: MutableStructure2D<Double>, b: MutableStructure2D<Double>): MutableStructure2D<Double> {
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<Double>, b: MutableStructure2D<Double>): MutableStructure2D<Double> {
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<Double>, epsilon: Double = 0.000001): MutableStructure2D<Double>? {
var idx = emptyArray<Double>()
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<Double>, MutableStructure2D<Double>, Int) -> MutableStructure2D<Double>,
t: MutableStructure2D<Double>, p: MutableStructure2D<Double>, exampleNumber: Int)
: MutableStructure2D<Double>
{
return func(t, p, exampleNumber)
}
private fun lmMatx(func: (MutableStructure2D<Double>, MutableStructure2D<Double>, Int) -> MutableStructure2D<Double>,
t: MutableStructure2D<Double>, pOld: MutableStructure2D<Double>, yOld: MutableStructure2D<Double>,
dX2: Int, JInput: MutableStructure2D<Double>, p: MutableStructure2D<Double>,
yDat: MutableStructure2D<Double>, weight: MutableStructure2D<Double>, dp:MutableStructure2D<Double>, settings:LMSettings) : Array<MutableStructure2D<Double>>
{
// 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<Double>, yOld: MutableStructure2D<Double>, JInput: MutableStructure2D<Double>,
p: MutableStructure2D<Double>, y: MutableStructure2D<Double>): MutableStructure2D<Double> {
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<Double>, MutableStructure2D<Double>, exampleNumber: Int) -> MutableStructure2D<Double>,
t: MutableStructure2D<Double>, p: MutableStructure2D<Double>, y: MutableStructure2D<Double>,
dp: MutableStructure2D<Double>, settings: LMSettings): MutableStructure2D<Double> {
// 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()
}

302
kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/linUtils.kt
View File

@ -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<Double>.svdGolubKahanHelper(u: MutableStructure2D<Double>, w: BufferedTensor<Double>,
v: MutableStructure2D<Double>, 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]
}
}
}

30
kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/tensorOps.kt
View File

@ -212,6 +212,36 @@ public fun DoubleTensorAlgebra.svd(
return Triple(uTensor.transposed(), sTensor, vTensor.transposed())
}
public fun DoubleTensorAlgebra.svdGolubKahan(
structureND: StructureND<Double>,
iterations: Int = 30, epsilon: Double = 1e-10
): Triple<DoubleTensor, DoubleTensor, DoubleTensor> {
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`.

4
kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleTensorAlgebra.kt
View File

@ -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

280
kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestLmAlgorithm.kt Normal file
View File

@ -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<Double>, p: MutableStructure2D<Double>, exampleNumber: Int): MutableStructure2D<Double> {
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<Double>, p: MutableStructure2D<Double>, exampleNumber: Int): MutableStructure2D<Double> {
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<Double>, p: MutableStructure2D<Double>, exampleNumber: Int): MutableStructure2D<Double> {
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])
}
}
}