forked from kscience/kmath
Merge pull request #513 from margarita0303/dev
Added Levenberg-Marquardt algorithm and svd Golub-Kahan
This commit is contained in:
SPC-code
GitHub
commit
7e46c7de4e
@ -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)
|
||||
}
|
||||
@ -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)
|
||||
}
|
||||
@ -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)
|
||||
}
|
||||
@ -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
|
||||
}
|
||||
@ -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()}")
|
||||
}
|
||||
@ -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)
|
||||
}
|
||||
@ -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
|
||||
|
||||
@ -4,7 +4,13 @@ plugins {
|
||||
|
||||
kscience{
|
||||
jvm()
|
||||
js()
|
||||
js {
|
||||
browser {
|
||||
testTask {
|
||||
useMocha().timeout = "0"
|
||||
}
|
||||
}
|
||||
}
|
||||
native()
|
||||
|
||||
dependencies {
|
||||
|
||||
@ -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>
|
||||
}
|
||||
|
||||
@ -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
|
||||
|
||||
|
||||
public val DoubleField.tensorAlgebra: DoubleTensorAlgebra get() = DoubleTensorAlgebra
|
||||
@ -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()
|
||||
}
|
||||
@ -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]
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -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`.
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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])
|
||||
}
|
||||
}
|
||||
}
|
||||
Loading…
Reference in New Issue
Block a user