kscience/kmath
9
0
Fork 1
Code Issues Pull Requests 6 Packages Projects Releases 16 Wiki Activity

levenbergMarquardt cleanup

Browse Source
This commit is contained in:
Alexander Nozik 2023-07-28 20:39:05 +03:00
parent 14f0fa1a6f
commit 1e2a8a40e5
3 changed files with 313 additions and 227 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
build.gradle.kts
View File

@ -15,7 +15,7 @@ allprojects {
} }
group = "space.kscience" group = "space.kscience"
version = "0.3.1" version = "0.3.2-dev-1"
} }
subprojects { subprojects {

308
kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/LevenbergMarquardtAlgorithm.kt
View File

@ -5,18 +5,13 @@
package space.kscience.kmath.tensors.core package space.kscience.kmath.tensors.core
import space.kscience.kmath.PerformancePitfall
import space.kscience.kmath.linear.transpose import space.kscience.kmath.linear.transpose
import space.kscience.kmath.nd.* import space.kscience.kmath.nd.*
import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra.div import kotlin.math.abs
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.max
import kotlin.math.min import kotlin.math.min
import kotlin.math.pow import kotlin.math.pow
import kotlin.reflect.KFunction3
/** /**
* Type of convergence achieved as a result of executing the Levenberg-Marquardt algorithm. * Type of convergence achieved as a result of executing the Levenberg-Marquardt algorithm.
@ -50,8 +45,8 @@ public enum class TypeOfConvergence {
* resultParameters: final parameters. * resultParameters: final parameters.
* typeOfConvergence: type of convergence. * typeOfConvergence: type of convergence.
*/ */
public data class LMResultInfo ( public data class LMResultInfo(
var iterations:Int, var iterations: Int,
var funcCalls: Int, var funcCalls: Int,
var resultChiSq: Double, var resultChiSq: Double,
var resultLambda: Double, var resultLambda: Double,
@ -87,7 +82,7 @@ public data class LMResultInfo (
* <8 - use maxParameters by default, <9 - use updateType 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. * exampleNumber: a parameter for a function with which you can choose its behavior.
*/ */
public data class LMInput ( public data class LMInput(
var func: (MutableStructure2D<Double>, MutableStructure2D<Double>, Int) -> (MutableStructure2D<Double>), var func: (MutableStructure2D<Double>, MutableStructure2D<Double>, Int) -> (MutableStructure2D<Double>),
var startParameters: MutableStructure2D<Double>, var startParameters: MutableStructure2D<Double>,
var independentVariables: MutableStructure2D<Double>, var independentVariables: MutableStructure2D<Double>,
@ -101,7 +96,7 @@ public data class LMInput (
var lambdas: DoubleArray, var lambdas: DoubleArray,
var updateType: Int, var updateType: Int,
var nargin: Int, var nargin: Int,
var exampleNumber: Int var exampleNumber: Int,
) )
@ -120,8 +115,10 @@ public data class LMInput (
* @return the 'output'. * @return the 'output'.
*/ */
public fun DoubleTensorAlgebra.levenbergMarquardt(inputData: LMInput): LMResultInfo { public fun DoubleTensorAlgebra.levenbergMarquardt(inputData: LMInput): LMResultInfo {
val resultInfo = LMResultInfo(0, 0, 0.0, val resultInfo = LMResultInfo(
0.0, inputData.startParameters, TypeOfConvergence.NoConvergence) 0, 0, 0.0,
0.0, inputData.startParameters, TypeOfConvergence.NoConvergence
)
val eps = 2.2204e-16 val eps = 2.2204e-16
@ -131,18 +128,21 @@ public fun DoubleTensorAlgebra.levenbergMarquardt(inputData: LMInput): LMResultI
var p = inputData.startParameters var p = inputData.startParameters
val t = inputData.independentVariables val t = inputData.independentVariables
val Npar = length(p) // number of parameters val nPar = length(p) // number of parameters
val Npnt = length(inputData.realValues) // number of data points val nPoints = length(inputData.realValues) // number of data points
var pOld = zeros(ShapeND(intArrayOf(Npar, 1))).as2D() // previous set of parameters 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 yOld = zeros(ShapeND(intArrayOf(nPoints, 1))).as2D() // previous model, y_old = y_hat(t;p_old)
var X2 = 1e-3 / eps // a really big initial Chi-sq value var x2: Double // a really big initial Chi-sq value
var X2Old = 1e-3 / eps // a really big initial Chi-sq value var x2Old: Double // a really big initial Chi-sq value
var J = zeros(ShapeND(intArrayOf(Npnt, Npar))).as2D() // Jacobian matrix var jacobian = zeros(ShapeND(intArrayOf(nPoints, nPar))).as2D() // Jacobian matrix
val DoF = Npnt - Npar // statistical degrees of freedom val dof = nPoints - nPar // statistical degrees of freedom
var weight = fromArray(ShapeND(intArrayOf(1, 1)), doubleArrayOf(inputData.weight)).as2D() var weight = fromArray(ShapeND(intArrayOf(1, 1)), doubleArrayOf(inputData.weight)).as2D()
if (inputData.nargin < 5) { if (inputData.nargin < 5) {
weight = fromArray(ShapeND(intArrayOf(1, 1)), doubleArrayOf((inputData.realValues.transpose().dot(inputData.realValues)).as1D()[0])).as2D() weight = fromArray(
ShapeND(intArrayOf(1, 1)),
doubleArrayOf((inputData.realValues.transpose().dot(inputData.realValues)).as1D()[0])
).as2D()
} }
var dp = inputData.pDelta var dp = inputData.pDelta
@ -175,7 +175,7 @@ public fun DoubleTensorAlgebra.levenbergMarquardt(inputData: LMInput): LMResultI
var updateType = inputData.updateType var updateType = inputData.updateType
if (inputData.nargin < 9) { if (inputData.nargin < 9) {
maxIterations = 10 * Npar maxIterations = 10 * nPar
epsilon1 = 1e-3 epsilon1 = 1e-3
epsilon2 = 1e-3 epsilon2 = 1e-3
epsilon3 = 1e-1 epsilon3 = 1e-1
@ -190,28 +190,27 @@ public fun DoubleTensorAlgebra.levenbergMarquardt(inputData: LMInput): LMResultI
maxParameters = makeColumn(maxParameters) maxParameters = makeColumn(maxParameters)
if (length(makeColumn(dp)) == 1) { if (length(makeColumn(dp)) == 1) {
dp = ones(ShapeND(intArrayOf(Npar, 1))).div(1 / dp[0, 0]).as2D() dp = ones(ShapeND(intArrayOf(nPar, 1))).div(1 / dp[0, 0]).as2D()
} }
var stop = false // termination flag var stop = false // termination flag
if (weight.shape.component1() == 1 || variance(weight) == 0.0) { // identical weights vector 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() weight = ones(ShapeND(intArrayOf(nPoints, 1))).div(1 / abs(weight[0, 0])).as2D()
} } else {
else {
weight = makeColumn(weight) weight = makeColumn(weight)
weight.abs() weight.abs()
} }
// initialize Jacobian with finite difference calculation // initialize Jacobian with finite difference calculation
var lmMatxAns = lmMatx(inputData.func, t, pOld, yOld, 1, J, p, inputData.realValues, weight, dp, settings) var lmMatxAns = lmMatx(inputData.func, t, pOld, yOld, 1, jacobian, p, inputData.realValues, weight, dp, settings)
var JtWJ = lmMatxAns[0] var jtWJ = lmMatxAns[0]
var JtWdy = lmMatxAns[1] var jtWdy = lmMatxAns[1]
X2 = lmMatxAns[2][0, 0] x2 = lmMatxAns[2][0, 0]
var yHat = lmMatxAns[3] var yHat = lmMatxAns[3]
J = lmMatxAns[4] jacobian = lmMatxAns[4]
if ( abs(JtWdy).max() < epsilon1 ) { if (abs(jtWdy).max() < epsilon1) {
stop = true stop = true
} }
@ -220,13 +219,12 @@ public fun DoubleTensorAlgebra.levenbergMarquardt(inputData: LMInput): LMResultI
if (updateType == 1) { if (updateType == 1) {
lambda = lambda0 // Marquardt: init'l lambda lambda = lambda0 // Marquardt: init'l lambda
} } else {
else { lambda = lambda0 * (makeColumnFromDiagonal(jtWJ)).max()
lambda = lambda0 * (makeColumnFromDiagonal(JtWJ)).max()
nu = 2 nu = 2
} }
X2Old = X2 // previous value of X2 x2Old = x2 // previous value of X2
var h: DoubleTensor var h: DoubleTensor
@ -235,20 +233,31 @@ public fun DoubleTensorAlgebra.levenbergMarquardt(inputData: LMInput): LMResultI
// incremental change in parameters // incremental change in parameters
h = if (updateType == 1) { // Marquardt h = if (updateType == 1) { // Marquardt
val solve = solve(JtWJ.plus(makeMatrixWithDiagonal(makeColumnFromDiagonal(JtWJ)).div(1 / lambda)).as2D(), JtWdy) val solve =
solve((jtWJ + makeMatrixWithDiagonal(makeColumnFromDiagonal(jtWJ)) / (1 / lambda)).as2D(), jtWdy)
solve.asDoubleTensor() solve.asDoubleTensor()
} else { // Quadratic and Nielsen } else { // Quadratic and Nielsen
val solve = solve(JtWJ.plus(lmEye(Npar).div(1 / lambda)).as2D(), JtWdy) val solve = solve(jtWJ.plus(lmEye(nPar) * lambda).as2D(), jtWdy)
solve.asDoubleTensor() solve.asDoubleTensor()
} }
var pTry = (p + h).as2D() // update the [idx] elements var pTry = (p + h).as2D() // update the [idx] elements
pTry = smallestElementComparison(largestElementComparison(minParameters, pTry.as2D()), maxParameters) // apply constraints pTry = smallestElementComparison(
largestElementComparison(minParameters, pTry),
maxParameters
) // apply constraints
var deltaY = inputData.realValues.minus(evaluateFunction(inputData.func, t, pTry, inputData.exampleNumber)) // residual error using p_try 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 (i in 0 until deltaY.shape[0]) { // floating point error; break
for (j in 0 until deltaY.shape.component2()) { for (j in 0 until deltaY.shape[1]) {
if (deltaY[i, j] == Double.POSITIVE_INFINITY || deltaY[i, j] == Double.NEGATIVE_INFINITY) { if (deltaY[i, j] == Double.POSITIVE_INFINITY || deltaY[i, j] == Double.NEGATIVE_INFINITY) {
stop = true stop = true
break break
@ -258,49 +267,72 @@ public fun DoubleTensorAlgebra.levenbergMarquardt(inputData: LMInput): LMResultI
settings.funcCalls += 1 settings.funcCalls += 1
val tmp = deltaY.times(weight) // val tmp = deltaY.times(weight)
var X2Try = deltaY.as2D().transpose().dot(tmp) // Chi-squared error criteria var X2Try = deltaY.as2D().transpose().dot(deltaY.times(weight)) // Chi-squared error criteria
val alpha = 1.0 val alpha = 1.0
if (updateType == 2) { // Quadratic if (updateType == 2) { // Quadratic
// One step of quadratic line update in the h direction for minimum X2 // 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))) val alpha = (jtWdy.transpose() dot h) / ((X2Try - x2) / 2.0 + 2 * (jtWdy.transpose() dot h))
h = h.dot(alpha) h = h dot alpha
pTry = p.plus(h).as2D() // update only [idx] elements pTry = (p + h).as2D() // update only [idx] elements
pTry = smallestElementComparison(largestElementComparison(minParameters, pTry), maxParameters) // apply constraints pTry = smallestElementComparison(
largestElementComparison(minParameters, pTry),
maxParameters
) // apply constraints
deltaY = inputData.realValues.minus(evaluateFunction(inputData.func, t, pTry, inputData.exampleNumber)) // residual error using p_try deltaY = inputData.realValues.minus(
evaluateFunction(
inputData.func,
t,
pTry,
inputData.exampleNumber
)
) // residual error using p_try
settings.funcCalls += 1 settings.funcCalls += 1
X2Try = deltaY.as2D().transpose().dot(deltaY.times(weight)) // Chi-squared error criteria X2Try = deltaY.as2D().transpose() dot deltaY * weight // Chi-squared error criteria
} }
val rho = when (updateType) { // Nielsen val rho = when (updateType) { // Nielsen
1 -> { 1 -> {
val tmp = h.transposed() val tmp = h.transposed()
.dot(makeMatrixWithDiagonal(makeColumnFromDiagonal(JtWJ)).div(1 / lambda).dot(h).plus(JtWdy)) .dot((makeMatrixWithDiagonal(makeColumnFromDiagonal(jtWJ)) * lambda dot h) + jtWdy)
X2.minus(X2Try).as2D()[0, 0] / abs(tmp.as2D()).as2D()[0, 0] (x2 - X2Try)[0, 0] / abs(tmp.as2D())[0, 0]
} }
else -> { else -> {
val tmp = h.transposed().dot(h.div(1 / lambda).plus(JtWdy)) val tmp = h.transposed().dot((h * lambda) + jtWdy)
X2.minus(X2Try).as2D()[0, 0] / abs(tmp.as2D()).as2D()[0, 0] x2.minus(X2Try).as2D()[0, 0] / abs(tmp.as2D()).as2D()[0, 0]
} }
} }
if (rho > epsilon4) { // it IS significantly better if (rho > epsilon4) { // it IS significantly better
val dX2 = X2.minus(X2Old) val dX2 = x2.minus(x2Old)
X2Old = X2 x2Old = x2
pOld = p.copyToTensor().as2D() pOld = p.copyToTensor().as2D()
yOld = yHat.copyToTensor().as2D() yOld = yHat.copyToTensor().as2D()
p = makeColumn(pTry) // accept p_try p = makeColumn(pTry) // accept p_try
lmMatxAns = lmMatx(inputData.func, t, pOld, yOld, dX2.toInt(), J, p, inputData.realValues, weight, dp, settings) lmMatxAns = lmMatx(
inputData.func,
t,
pOld,
yOld,
dX2.toInt(),
jacobian,
p,
inputData.realValues,
weight,
dp,
settings
)
// decrease lambda ==> Gauss-Newton method // decrease lambda ==> Gauss-Newton method
JtWJ = lmMatxAns[0] jtWJ = lmMatxAns[0]
JtWdy = lmMatxAns[1] jtWdy = lmMatxAns[1]
X2 = lmMatxAns[2][0, 0] x2 = lmMatxAns[2][0, 0]
yHat = lmMatxAns[3] yHat = lmMatxAns[3]
J = lmMatxAns[4] jacobian = lmMatxAns[4]
lambda = when (updateType) { lambda = when (updateType) {
1 -> { // Levenberg 1 -> { // Levenberg
@ -317,13 +349,14 @@ public fun DoubleTensorAlgebra.levenbergMarquardt(inputData: LMInput): LMResultI
} }
} }
} else { // it IS NOT better } else { // it IS NOT better
X2 = X2Old // do not accept p_try x2 = x2Old // do not accept p_try
if (settings.iteration % (2 * Npar) == 0) { // rank-1 update of Jacobian 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) lmMatxAns =
JtWJ = lmMatxAns[0] lmMatx(inputData.func, t, pOld, yOld, -1, jacobian, p, inputData.realValues, weight, dp, settings)
JtWdy = lmMatxAns[1] jtWJ = lmMatxAns[0]
jtWdy = lmMatxAns[1]
yHat = lmMatxAns[3] yHat = lmMatxAns[3]
J = lmMatxAns[4] jacobian = lmMatxAns[4]
} }
// increase lambda ==> gradient descent method // increase lambda ==> gradient descent method
@ -333,7 +366,7 @@ public fun DoubleTensorAlgebra.levenbergMarquardt(inputData: LMInput): LMResultI
} }
2 -> { // Quadratic 2 -> { // Quadratic
lambda + kotlin.math.abs(((X2Try.as2D()[0, 0] - X2) / 2) / alpha) lambda + abs(((X2Try.as2D()[0, 0] - x2) / 2) / alpha)
} }
else -> { // Nielsen else -> { // Nielsen
@ -343,7 +376,7 @@ public fun DoubleTensorAlgebra.levenbergMarquardt(inputData: LMInput): LMResultI
} }
} }
val chiSq = X2 / DoF val chiSq = x2 / dof
resultInfo.iterations = settings.iteration resultInfo.iterations = settings.iteration
resultInfo.funcCalls = settings.funcCalls resultInfo.funcCalls = settings.funcCalls
resultInfo.resultChiSq = chiSq resultInfo.resultChiSq = chiSq
@ -351,7 +384,7 @@ public fun DoubleTensorAlgebra.levenbergMarquardt(inputData: LMInput): LMResultI
resultInfo.resultParameters = p resultInfo.resultParameters = p
if (abs(JtWdy).max() < epsilon1 && settings.iteration > 2) { if (abs(jtWdy).max() < epsilon1 && settings.iteration > 2) {
resultInfo.typeOfConvergence = TypeOfConvergence.InGradient resultInfo.typeOfConvergence = TypeOfConvergence.InGradient
stop = true stop = true
} }
@ -359,7 +392,7 @@ public fun DoubleTensorAlgebra.levenbergMarquardt(inputData: LMInput): LMResultI
resultInfo.typeOfConvergence = TypeOfConvergence.InParameters resultInfo.typeOfConvergence = TypeOfConvergence.InParameters
stop = true stop = true
} }
if (X2 / DoF < epsilon3 && settings.iteration > 2) { if (x2 / dof < epsilon3 && settings.iteration > 2) {
resultInfo.typeOfConvergence = TypeOfConvergence.InReducedChiSquare resultInfo.typeOfConvergence = TypeOfConvergence.InReducedChiSquare
stop = true stop = true
} }
@ -371,10 +404,10 @@ public fun DoubleTensorAlgebra.levenbergMarquardt(inputData: LMInput): LMResultI
return resultInfo return resultInfo
} }
private data class LMSettings ( private data class LMSettings(
var iteration:Int, var iteration: Int,
var funcCalls: Int, var funcCalls: Int,
var exampleNumber:Int var exampleNumber: Int,
) )
/* matrix -> column of all elements */ /* matrix -> column of all elements */
@ -390,14 +423,14 @@ private fun makeColumn(tensor: MutableStructure2D<Double>): MutableStructure2D<D
} }
/* column length */ /* column length */
private fun length(column: MutableStructure2D<Double>) : Int { private fun length(column: MutableStructure2D<Double>): Int {
return column.shape.component1() return column.shape.component1()
} }
private fun MutableStructure2D<Double>.abs() { private fun MutableStructure2D<Double>.abs() {
for (i in 0 until this.shape.component1()) { for (i in 0 until this.shape[0]) {
for (j in 0 until this.shape.component2()) { for (j in 0 until this.shape[1]) {
this[i, j] = kotlin.math.abs(this[i, j]) this[i, j] = abs(this[i, j])
} }
} }
} }
@ -413,7 +446,7 @@ private fun abs(input: MutableStructure2D<Double>): MutableStructure2D<Double> {
).as2D() ).as2D()
for (i in 0 until tensor.shape.component1()) { for (i in 0 until tensor.shape.component1()) {
for (j in 0 until tensor.shape.component2()) { for (j in 0 until tensor.shape.component2()) {
tensor[i, j] = kotlin.math.abs(input[i, j]) tensor[i, j] = abs(input[i, j])
} }
} }
return tensor return tensor
@ -441,21 +474,23 @@ private fun lmEye(size: Int): MutableStructure2D<Double> {
return makeMatrixWithDiagonal(column) return makeMatrixWithDiagonal(column)
} }
private fun largestElementComparison(a: MutableStructure2D<Double>, b: MutableStructure2D<Double>): MutableStructure2D<Double> { private fun largestElementComparison(
a: MutableStructure2D<Double>,
b: MutableStructure2D<Double>,
): MutableStructure2D<Double> {
val aSizeX = a.shape.component1() val aSizeX = a.shape.component1()
val aSizeY = a.shape.component2() val aSizeY = a.shape.component2()
val bSizeX = b.shape.component1() val bSizeX = b.shape.component1()
val bSizeY = b.shape.component2() val bSizeY = b.shape.component2()
val tensor = BroadcastDoubleTensorAlgebra.zeros(ShapeND(intArrayOf(max(aSizeX, bSizeX), max(aSizeY, bSizeY)))).as2D() val tensor =
BroadcastDoubleTensorAlgebra.zeros(ShapeND(intArrayOf(max(aSizeX, bSizeX), max(aSizeY, bSizeY)))).as2D()
for (i in 0 until tensor.shape.component1()) { for (i in 0 until tensor.shape.component1()) {
for (j in 0 until tensor.shape.component2()) { for (j in 0 until tensor.shape.component2()) {
if (i < aSizeX && i < bSizeX && j < aSizeY && j < bSizeY) { if (i < aSizeX && i < bSizeX && j < aSizeY && j < bSizeY) {
tensor[i, j] = max(a[i, j], b[i, j]) tensor[i, j] = max(a[i, j], b[i, j])
} } else if (i < aSizeX && j < aSizeY) {
else if (i < aSizeX && j < aSizeY) {
tensor[i, j] = a[i, j] tensor[i, j] = a[i, j]
} } else {
else {
tensor[i, j] = b[i, j] tensor[i, j] = b[i, j]
} }
} }
@ -463,21 +498,23 @@ private fun largestElementComparison(a: MutableStructure2D<Double>, b: MutableSt
return tensor return tensor
} }
private fun smallestElementComparison(a: MutableStructure2D<Double>, b: MutableStructure2D<Double>): MutableStructure2D<Double> { private fun smallestElementComparison(
a: MutableStructure2D<Double>,
b: MutableStructure2D<Double>,
): MutableStructure2D<Double> {
val aSizeX = a.shape.component1() val aSizeX = a.shape.component1()
val aSizeY = a.shape.component2() val aSizeY = a.shape.component2()
val bSizeX = b.shape.component1() val bSizeX = b.shape.component1()
val bSizeY = b.shape.component2() val bSizeY = b.shape.component2()
val tensor = BroadcastDoubleTensorAlgebra.zeros(ShapeND(intArrayOf(max(aSizeX, bSizeX), max(aSizeY, bSizeY)))).as2D() val tensor =
BroadcastDoubleTensorAlgebra.zeros(ShapeND(intArrayOf(max(aSizeX, bSizeX), max(aSizeY, bSizeY)))).as2D()
for (i in 0 until tensor.shape.component1()) { for (i in 0 until tensor.shape.component1()) {
for (j in 0 until tensor.shape.component2()) { for (j in 0 until tensor.shape.component2()) {
if (i < aSizeX && i < bSizeX && j < aSizeY && j < bSizeY) { if (i < aSizeX && i < bSizeX && j < aSizeY && j < bSizeY) {
tensor[i, j] = min(a[i, j], b[i, j]) tensor[i, j] = min(a[i, j], b[i, j])
} } else if (i < aSizeX && j < aSizeY) {
else if (i < aSizeX && j < aSizeY) {
tensor[i, j] = a[i, j] tensor[i, j] = a[i, j]
} } else {
else {
tensor[i, j] = b[i, j] tensor[i, j] = b[i, j]
} }
} }
@ -485,10 +522,13 @@ private fun smallestElementComparison(a: MutableStructure2D<Double>, b: MutableS
return tensor return tensor
} }
private fun getZeroIndices(column: MutableStructure2D<Double>, epsilon: Double = 0.000001): MutableStructure2D<Double>? { private fun getZeroIndices(
column: MutableStructure2D<Double>,
epsilon: Double = 0.000001,
): MutableStructure2D<Double>? {
var idx = emptyArray<Double>() var idx = emptyArray<Double>()
for (i in 0 until column.shape.component1()) { for (i in 0 until column.shape.component1()) {
if (kotlin.math.abs(column[i, 0]) > epsilon) { if (abs(column[i, 0]) > epsilon) {
idx += (i + 1.0) idx += (i + 1.0)
} }
} }
@ -498,18 +538,27 @@ private fun getZeroIndices(column: MutableStructure2D<Double>, epsilon: Double =
return null return null
} }
private fun evaluateFunction(func: (MutableStructure2D<Double>, MutableStructure2D<Double>, Int) -> MutableStructure2D<Double>, private fun evaluateFunction(
t: MutableStructure2D<Double>, p: MutableStructure2D<Double>, exampleNumber: Int) func: (MutableStructure2D<Double>, MutableStructure2D<Double>, Int) -> MutableStructure2D<Double>,
: MutableStructure2D<Double> t: MutableStructure2D<Double>, p: MutableStructure2D<Double>, exampleNumber: Int,
{ )
: MutableStructure2D<Double> {
return func(t, p, exampleNumber) return func(t, p, exampleNumber)
} }
private fun lmMatx(func: (MutableStructure2D<Double>, MutableStructure2D<Double>, Int) -> MutableStructure2D<Double>, private fun lmMatx(
t: MutableStructure2D<Double>, pOld: MutableStructure2D<Double>, yOld: MutableStructure2D<Double>, func: (MutableStructure2D<Double>, MutableStructure2D<Double>, Int) -> MutableStructure2D<Double>,
dX2: Int, JInput: MutableStructure2D<Double>, p: MutableStructure2D<Double>, t: MutableStructure2D<Double>,
yDat: MutableStructure2D<Double>, weight: MutableStructure2D<Double>, dp:MutableStructure2D<Double>, settings:LMSettings) : Array<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>> = with(DoubleTensorAlgebra) {
// default: dp = 0.001 // default: dp = 0.001
val Npar = length(p) // number of parameters val Npar = length(p) // number of parameters
@ -520,63 +569,70 @@ private fun lmMatx(func: (MutableStructure2D<Double>, MutableStructure2D<Double>
J = if (settings.iteration % (2 * Npar) == 0 || dX2 > 0) { J = if (settings.iteration % (2 * Npar) == 0 || dX2 > 0) {
lmFdJ(func, t, p, yHat, dp, settings).as2D() // finite difference lmFdJ(func, t, p, yHat, dp, settings).as2D() // finite difference
} } else {
else {
lmBroydenJ(pOld, yOld, J, p, yHat).as2D() // rank-1 update lmBroydenJ(pOld, yOld, J, p, yHat).as2D() // rank-1 update
} }
val deltaY = yDat.minus(yHat) val deltaY = yDat.minus(yHat)
val chiSq = deltaY.transposed().dot( deltaY.times(weight) ).as2D() 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 JtWJ =
val JtWdy = J.transposed().dot( weight.times(deltaY) ).as2D() (J.transposed() dot J * (weight dot ones(ShapeND(intArrayOf(1, Npar))))).as2D()
val JtWdy = (J.transposed() dot weight * deltaY).as2D()
return arrayOf(JtWJ,JtWdy,chiSq,yHat,J) return arrayOf(JtWJ, JtWdy, chiSq, yHat, J)
} }
private fun lmBroydenJ(pOld: MutableStructure2D<Double>, yOld: MutableStructure2D<Double>, JInput: MutableStructure2D<Double>, private fun lmBroydenJ(
p: MutableStructure2D<Double>, y: MutableStructure2D<Double>): MutableStructure2D<Double> { pOld: MutableStructure2D<Double>, yOld: MutableStructure2D<Double>, JInput: MutableStructure2D<Double>,
p: MutableStructure2D<Double>, y: MutableStructure2D<Double>,
): MutableStructure2D<Double> = with(DoubleTensorAlgebra) {
var J = JInput.copyToTensor() var J = JInput.copyToTensor()
val h = p.minus(pOld) 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] ) val increase = ((y - yOld - (J dot h)) dot h.transposed()) / (h.transposed() dot h)[0, 0]
J = J.plus(increase) J = J.plus(increase)
return J.as2D() return J.as2D()
} }
private fun lmFdJ(func: (MutableStructure2D<Double>, MutableStructure2D<Double>, exampleNumber: Int) -> MutableStructure2D<Double>, @OptIn(PerformancePitfall::class)
t: MutableStructure2D<Double>, p: MutableStructure2D<Double>, y: MutableStructure2D<Double>, private fun lmFdJ(
dp: MutableStructure2D<Double>, settings: LMSettings): MutableStructure2D<Double> { 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> = with(DoubleTensorAlgebra) {
// default: dp = 0.001 * ones(1,n) // default: dp = 0.001 * ones(1,n)
val m = length(y) // number of data points val m = length(y) // number of data points
val n = length(p) // number of parameters val n = length(p) // number of parameters
val ps = p.copyToTensor().as2D() val ps = p.copyToTensor()
val J = BroadcastDoubleTensorAlgebra.zeros(ShapeND(intArrayOf(m, n))).as2D() // initialize Jacobian to Zero val J = zero(m, n) // initialize Jacobian to Zero
val del = BroadcastDoubleTensorAlgebra.zeros(ShapeND(intArrayOf(n, 1))).as2D() val del = zero(n, 1)
for (j in 0 until n) { for (j in 0 until n) {
del[j, 0] = dp[j, 0] * (1 + kotlin.math.abs(p[j, 0])) // parameter perturbation del[j, 0] = dp[j, 0] * (1 + abs(p[j, 0])) // parameter perturbation
p[j, 0] = ps[j, 0] + del[j, 0] // perturb parameter p(j) p[j, 0] = ps[j, 0] + del[j, 0] // perturb parameter p(j)
val epsilon = 0.0000001 val epsilon = 0.0000001
if (kotlin.math.abs(del[j, 0]) > epsilon) { if (abs(del[j, 0]) > epsilon) {
val y1 = evaluateFunction(func, t, p, settings.exampleNumber) val y1 = evaluateFunction(func, t, p, settings.exampleNumber)
settings.funcCalls += 1 settings.funcCalls += 1
if (dp[j, 0] < 0) { // backwards difference if (dp[j, 0] < 0) { // backwards difference
for (i in 0 until J.shape.component1()) { for (i in 0 until J.shape.first()) {
J[i, j] = (y1.as2D().minus(y).as2D())[i, 0] / del[j, 0] J[i, j] = (y1 - y)[i, 0] / del[j, 0]
} }
} } else {
else {
// Do tests for it // Do tests for it
p[j, 0] = ps[j, 0] - del[j, 0] // central difference, additional func call p[j, 0] = ps[j, 0] - del[j, 0] // central difference, additional func call
for (i in 0 until J.shape.component1()) { for (i in 0 until J.shape.first()) {
J[i, j] = (y1.as2D().minus(evaluateFunction(func, t, p, settings.exampleNumber)).as2D())[i, 0] / (2 * del[j, 0]) J[i, j] = (y1 - evaluateFunction(func, t, p, settings.exampleNumber))[i, 0] / (2 * del[j, 0])
} }
settings.funcCalls += 1 settings.funcCalls += 1
} }

218
kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestLmAlgorithm.kt → kmath-tensors/src/jvmTest/kotlin/space/kscience/kmath/tensors/core/TestLmAlgorithm.kt
View File

@ -5,77 +5,83 @@
package space.kscience.kmath.tensors.core package space.kscience.kmath.tensors.core
import space.kscience.kmath.nd.MutableStructure2D import space.kscience.kmath.PerformancePitfall
import space.kscience.kmath.nd.ShapeND import space.kscience.kmath.nd.*
import space.kscience.kmath.nd.as2D
import space.kscience.kmath.nd.component1
import space.kscience.kmath.operations.invoke import space.kscience.kmath.operations.invoke
import space.kscience.kmath.tensors.core.DoubleTensorAlgebra.Companion.max import space.kscience.kmath.structures.DoubleBuffer
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.Test
import kotlin.test.assertEquals import kotlin.test.assertEquals
@PerformancePitfall
class TestLmAlgorithm { class TestLmAlgorithm {
companion object { companion object {
fun funcEasyForLm(t: MutableStructure2D<Double>, p: MutableStructure2D<Double>, exampleNumber: Int): MutableStructure2D<Double> { fun funcEasyForLm(
t: MutableStructure2D<Double>,
p: MutableStructure2D<Double>,
exampleNumber: Int,
): MutableStructure2D<Double> = with(DoubleTensorAlgebra) {
val m = t.shape.component1() val m = t.shape.component1()
var yHat = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf(m, 1))) val yHat = when (exampleNumber) {
1 -> exp((t * (-1.0 / p[1, 0]))) * p[0, 0] + (t * p[2, 0]) * exp((t * (-1.0 / p[3, 0])))
if (exampleNumber == 1) { 2 -> {
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() val mt = t.max()
yHat = (t.times(1.0 / mt)).times(p[0, 0]) + (t * (1.0 / mt)) * p[0, 0] +
(t.times(1.0 / mt)).pow(2).times(p[1, 0]) + (t * (1.0 / mt)).pow(2) * p[1, 0] +
(t.times(1.0 / mt)).pow(3).times(p[2, 0]) + (t * (1.0 / mt)).pow(3) * p[2, 0] +
(t.times(1.0 / mt)).pow(4).times(p[3, 0]) (t * (1.0 / mt)).pow(4) * p[3, 0]
} }
else if (exampleNumber == 3) {
yHat = DoubleTensorAlgebra.exp((t.times(-1.0 / p[1, 0]))) 3 -> exp(t * (-1.0 / p[1, 0])) * p[0, 0] +
.times(p[0, 0]) + DoubleTensorAlgebra.sin((t.times(1.0 / p[3, 0]))).times(p[2, 0]) sin((t * (1.0 / p[3, 0]))) * p[2, 0]
else -> zeros(ShapeND(intArrayOf(m, 1)))
} }
return yHat.as2D() return yHat.as2D()
} }
fun funcMiddleForLm(t: MutableStructure2D<Double>, p: MutableStructure2D<Double>, exampleNumber: Int): MutableStructure2D<Double> { fun funcMiddleForLm(
t: MutableStructure2D<Double>,
p: MutableStructure2D<Double>,
exampleNumber: Int,
): MutableStructure2D<Double> = with(DoubleTensorAlgebra) {
val m = t.shape.component1() val m = t.shape.component1()
var yHat = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf (m, 1))) var yHat = zeros(ShapeND(intArrayOf(m, 1)))
val mt = t.max() val mt = t.max()
for(i in 0 until p.shape.component1()){ for (i in 0 until p.shape.component1()) {
yHat += (t.times(1.0 / mt)).times(p[i, 0]) yHat.plusAssign(t * (1.0 / mt) * p[i, 0])
} }
for(i in 0 until 5){ for (i in 0 until 5) {
yHat = funcEasyForLm(yHat.as2D(), p, exampleNumber).asDoubleTensor() yHat = funcEasyForLm(yHat.as2D(), p, exampleNumber).asDoubleTensor()
} }
return yHat.as2D() return yHat.as2D()
} }
fun funcDifficultForLm(t: MutableStructure2D<Double>, p: MutableStructure2D<Double>, exampleNumber: Int): MutableStructure2D<Double> { fun funcDifficultForLm(
t: MutableStructure2D<Double>,
p: MutableStructure2D<Double>,
exampleNumber: Int,
): MutableStructure2D<Double> = with(DoubleTensorAlgebra) {
val m = t.shape.component1() val m = t.shape.component1()
var yHat = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf (m, 1))) var yHat = zeros(ShapeND(intArrayOf(m, 1)))
val mt = t.max() val mt = t.max()
for(i in 0 until p.shape.component1()){ for (i in 0 until p.shape.component1()) {
yHat = yHat.plus( (t.times(1.0 / mt)).times(p[i, 0]) ) yHat = yHat + (t * (1.0 / mt)) * p[i, 0]
} }
for(i in 0 until 4){ for (i in 0 until 4) {
yHat = funcEasyForLm((yHat.as2D() + t).as2D(), p, exampleNumber).asDoubleTensor() yHat = funcEasyForLm((yHat.as2D() + t).as2D(), p, exampleNumber).asDoubleTensor()
} }
return yHat.as2D() return yHat.as2D()
} }
} }
@Test @Test
fun testLMEasy() = DoubleTensorAlgebra { fun testLMEasy() = DoubleTensorAlgebra {
val lmMatxYDat = doubleArrayOf( val lmMatxYDat = doubleArrayOf(
@ -91,12 +97,12 @@ class TestLmAlgorithm {
14.7665, 13.3718, 15.0587, 13.8320, 14.7873, 13.6824, 14.2579, 14.2154, 13.5818, 13.8157 14.7665, 13.3718, 15.0587, 13.8320, 14.7873, 13.6824, 14.2579, 14.2154, 13.5818, 13.8157
) )
var exampleNumber = 1 val exampleNumber = 1
val p_init = BroadcastDoubleTensorAlgebra.fromArray( val pInit = fromArray(
ShapeND(intArrayOf(4, 1)), doubleArrayOf(5.0, 2.0, 0.2, 10.0) ShapeND(intArrayOf(4, 1)), doubleArrayOf(5.0, 2.0, 0.2, 10.0)
).as2D() ).as2D()
var t = ones(ShapeND(intArrayOf(100, 1))).as2D() val t = ones(ShapeND(intArrayOf(100, 1))).as2D()
for (i in 0 until 100) { for (i in 0 until 100) {
t[i, 0] = t[i, 0] * (i + 1) t[i, 0] = t[i, 0] * (i + 1)
} }
@ -119,22 +125,26 @@ class TestLmAlgorithm {
ShapeND(intArrayOf(4, 1)), doubleArrayOf(50.0, 20.0, 2.0, 100.0) ShapeND(intArrayOf(4, 1)), doubleArrayOf(50.0, 20.0, 2.0, 100.0)
).as2D() ).as2D()
val inputData = LMInput(::funcEasyForLm, p_init, t, yDat, weight, dp, pMin, pMax, 100, val inputData = LMInput(
doubleArrayOf(1e-3, 1e-3, 1e-1, 1e-1), doubleArrayOf(1e-2, 11.0, 9.0), 1, 10, exampleNumber) Companion::funcEasyForLm, pInit, 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) val result = levenbergMarquardt(inputData)
assertEquals(13, result.iterations) assertEquals(13, result.iterations)
assertEquals(31, result.funcCalls) assertEquals(31, result.funcCalls)
assertEquals(0.9131368192633, (result.resultChiSq * 1e13).roundToLong() / 1e13) assertEquals(0.9131368192633, result.resultChiSq, 1e-13)
assertEquals(3.7790980 * 1e-7, (result.resultLambda * 1e13).roundToLong() / 1e13) assertEquals(3.7790980 * 1e-7, result.resultLambda, 1e-13)
assertEquals(result.typeOfConvergence, TypeOfConvergence.InParameters) assertEquals(result.typeOfConvergence, TypeOfConvergence.InParameters)
val expectedParameters = BroadcastDoubleTensorAlgebra.fromArray( val expectedParameters = BroadcastDoubleTensorAlgebra.fromArray(
ShapeND(intArrayOf(4, 1)), doubleArrayOf(20.527230909086, 9.833627103230, 0.997571256572, 50.174445822506) ShapeND(intArrayOf(4, 1)), doubleArrayOf(20.527230909086, 9.833627103230, 0.997571256572, 50.174445822506)
).as2D() ).as2D()
result.resultParameters = result.resultParameters.map { x -> (x * 1e12).toLong() / 1e12}.as2D() result.resultParameters = result.resultParameters.map { x -> (x * 1e12).toLong() / 1e12 }.as2D()
val receivedParameters = BroadcastDoubleTensorAlgebra.fromArray( val receivedParameters = BroadcastDoubleTensorAlgebra.fromArray(
ShapeND(intArrayOf(4, 1)), doubleArrayOf(result.resultParameters[0, 0], result.resultParameters[1, 0], ShapeND(intArrayOf(4, 1)), doubleArrayOf(
result.resultParameters[2, 0], result.resultParameters[3, 0]) result.resultParameters[0, 0], result.resultParameters[1, 0],
result.resultParameters[2, 0], result.resultParameters[3, 0]
)
).as2D() ).as2D()
assertEquals(expectedParameters[0, 0], receivedParameters[0, 0]) assertEquals(expectedParameters[0, 0], receivedParameters[0, 0])
assertEquals(expectedParameters[1, 0], receivedParameters[1, 0]) assertEquals(expectedParameters[1, 0], receivedParameters[1, 0])
@ -143,16 +153,16 @@ class TestLmAlgorithm {
} }
@Test @Test
fun TestLMMiddle() = DoubleTensorAlgebra { fun testLMMiddle() = DoubleTensorAlgebra {
val NData = 100 val nData = 100
val tExample = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(NData, 1))).as2D() val tExample = one(nData, 1).as2D()
for (i in 0 until NData) { for (i in 0 until nData) {
tExample[i, 0] = tExample[i, 0] * (i + 1) tExample[i, 0] = tExample[i, 0] * (i + 1)
} }
val Nparams = 20 val nParams = 20
val pExample = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1))).as2D() val pExample = one(nParams, 1).as2D()
for (i in 0 until Nparams) { for (i in 0 until nParams) {
pExample[i, 0] = pExample[i, 0] + i - 25 pExample[i, 0] = pExample[i, 0] + i - 25
} }
@ -160,8 +170,8 @@ class TestLmAlgorithm {
val yHat = funcMiddleForLm(tExample, pExample, exampleNumber) val yHat = funcMiddleForLm(tExample, pExample, exampleNumber)
val pInit = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf(Nparams, 1))).as2D() val pInit = zeros(ShapeND(intArrayOf(nParams, 1))).as2D()
for (i in 0 until Nparams) { for (i in 0 until nParams) {
pInit[i, 0] = (pExample[i, 0] + 0.9) pInit[i, 0] = (pExample[i, 0] + 0.9)
} }
@ -171,13 +181,14 @@ class TestLmAlgorithm {
val dp = BroadcastDoubleTensorAlgebra.fromArray( val dp = BroadcastDoubleTensorAlgebra.fromArray(
ShapeND(intArrayOf(1, 1)), DoubleArray(1) { -0.01 } ShapeND(intArrayOf(1, 1)), DoubleArray(1) { -0.01 }
).as2D() ).as2D()
var pMin = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1))) var pMin = ones(ShapeND(intArrayOf(nParams, 1)))
pMin = pMin.div(1.0 / -50.0) pMin = pMin * (-50.0)
val pMax = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1))) val pMax = ones(ShapeND(intArrayOf(nParams, 1)))
pMin = pMin.div(1.0 / 50.0) pMin = pMin * 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 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, val inputData = LMInput(
Companion::funcMiddleForLm,
pInit.as2D(), pInit.as2D(),
t, t,
yDat, yDat,
@ -190,63 +201,67 @@ class TestLmAlgorithm {
doubleArrayOf(opts[6], opts[7], opts[8]), doubleArrayOf(opts[6], opts[7], opts[8]),
opts[9].toInt(), opts[9].toInt(),
10, 10,
1) 1
)
val result = DoubleTensorAlgebra.levenbergMarquardt(inputData) val result = DoubleTensorAlgebra.levenbergMarquardt(inputData)
assertEquals(46, result.iterations) assertEquals(46, result.iterations)
assertEquals(113, result.funcCalls) assertEquals(113, result.funcCalls)
assertEquals(0.000005977, (result.resultChiSq * 1e9).roundToLong() / 1e9) assertEquals(0.000005977, result.resultChiSq, 1e-9)
assertEquals(1.0 * 1e-7, (result.resultLambda * 1e13).roundToLong() / 1e13) assertEquals(1.0 * 1e-7, result.resultLambda, 1e-13)
assertEquals(result.typeOfConvergence, TypeOfConvergence.InReducedChiSquare) assertEquals(result.typeOfConvergence, TypeOfConvergence.InReducedChiSquare)
val expectedParameters = BroadcastDoubleTensorAlgebra.fromArray( val expectedParameters = fromArray(
ShapeND(intArrayOf(Nparams, 1)), doubleArrayOf( -23.9717, -18.6686, -21.7971, 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, - -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() 12.0201, -10.9256, -9.9878, -9.2309, -8.6589, -8.2365, -7.8783, -7.4598, -6.8511
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) { val receivedParameters = zero(nParams, 1)
for (i in 0 until nParams) {
receivedParameters[i, 0] = result.resultParameters[i, 0] receivedParameters[i, 0] = result.resultParameters[i, 0]
assertEquals(expectedParameters[i, 0], result.resultParameters[i, 0]) assertEquals(expectedParameters[i, 0], result.resultParameters[i, 0], 1e-2)
} }
} }
@Test @Test
fun TestLMDifficult() = DoubleTensorAlgebra { fun TestLMDifficult() = DoubleTensorAlgebra {
val NData = 200 val nData = 200
var tExample = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(NData, 1))).as2D() val tExample = ones(ShapeND(intArrayOf(nData, 1))).as2D()
for (i in 0 until NData) { for (i in 0 until nData) {
tExample[i, 0] = tExample[i, 0] * (i + 1) - 104 tExample[i, 0] = tExample[i, 0] * (i + 1) - 104
} }
val Nparams = 15 val nParams = 15
var pExample = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1))).as2D() val pExample = ones(ShapeND(intArrayOf(nParams, 1))).as2D()
for (i in 0 until Nparams) { for (i in 0 until nParams) {
pExample[i, 0] = pExample[i, 0] + i - 25 pExample[i, 0] = pExample[i, 0] + i - 25
} }
val exampleNumber = 1 val exampleNumber = 1
var yHat = funcDifficultForLm(tExample, pExample, exampleNumber) val yHat = funcDifficultForLm(tExample, pExample, exampleNumber)
var pInit = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf(Nparams, 1))).as2D() val pInit = zeros(ShapeND(intArrayOf(nParams, 1))).as2D()
for (i in 0 until Nparams) { for (i in 0 until nParams) {
pInit[i, 0] = (pExample[i, 0] + 0.9) pInit[i, 0] = (pExample[i, 0] + 0.9)
} }
var t = tExample val t = tExample
val yDat = yHat val yDat = yHat
val weight = 1.0 / Nparams * 1.0 - 0.085 val weight = 1.0 / nParams * 1.0 - 0.085
val dp = BroadcastDoubleTensorAlgebra.fromArray( val dp = fromArray(
ShapeND(intArrayOf(1, 1)), DoubleArray(1) { -0.01 } ShapeND(intArrayOf(1, 1)), DoubleArray(1) { -0.01 }
).as2D() ).as2D()
var pMin = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1))) var pMin = ones(ShapeND(intArrayOf(nParams, 1)))
pMin = pMin.div(1.0 / -50.0) pMin = pMin * (-50.0)
val pMax = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1))) val pMax = ones(ShapeND(intArrayOf(nParams, 1)))
pMin = pMin.div(1.0 / 50.0) pMin = pMin * (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 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, val inputData = LMInput(
Companion::funcDifficultForLm,
pInit.as2D(), pInit.as2D(),
t, t,
yDat, yDat,
@ -259,22 +274,37 @@ class TestLmAlgorithm {
doubleArrayOf(opts[6], opts[7], opts[8]), doubleArrayOf(opts[6], opts[7], opts[8]),
opts[9].toInt(), opts[9].toInt(),
10, 10,
1) 1
)
val result = DoubleTensorAlgebra.levenbergMarquardt(inputData) val result = DoubleTensorAlgebra.levenbergMarquardt(inputData)
assertEquals(2375, result.iterations) assertEquals(2375, result.iterations)
assertEquals(4858, result.funcCalls) assertEquals(4858, result.funcCalls)
assertEquals(5.14347, (result.resultLambda * 1e5).roundToLong() / 1e5) assertEquals(5.14347, result.resultLambda, 1e-5)
assertEquals(result.typeOfConvergence, TypeOfConvergence.InParameters) assertEquals(result.typeOfConvergence, TypeOfConvergence.InParameters)
val expectedParameters = BroadcastDoubleTensorAlgebra.fromArray( val expectedParameters = DoubleBuffer(
ShapeND(intArrayOf(Nparams, 1)), doubleArrayOf(-23.6412, -16.7402, -21.5705, -21.0464, -23.6412,
-17.2852, -17.2959, -17.298, 0.9999, -17.2885, -17.3008, -17.2941, -17.2923, -17.2976, -17.3028, -17.2891)).as2D() -16.7402,
result.resultParameters = result.resultParameters.map { x -> (x * 1e4).roundToLong() / 1e4}.as2D() -21.5705,
val receivedParameters = zeros(ShapeND(intArrayOf(Nparams, 1))).as2D() -21.0464,
for (i in 0 until Nparams) { -17.2852,
-17.2959,
-17.298,
0.9999,
-17.2885,
-17.3008,
-17.2941,
-17.2923,
-17.2976,
-17.3028,
-17.2891
)
val receivedParameters = zeros(ShapeND(intArrayOf(nParams, 1))).as2D()
for (i in 0 until nParams) {
receivedParameters[i, 0] = result.resultParameters[i, 0] receivedParameters[i, 0] = result.resultParameters[i, 0]
assertEquals(expectedParameters[i, 0], result.resultParameters[i, 0]) assertEquals(expectedParameters[i], result.resultParameters[i, 0], 1e-2)
} }
} }
} }