From 2fa39fff14352df710c0495082a3903389a4d005 Mon Sep 17 00:00:00 2001 From: Margarita Date: Tue, 24 May 2022 19:22:26 +0300 Subject: [PATCH 1/6] added partial implementation of svd calculation --- .../space/kscience/kmath/tensors/margarita.kt | 76 +++++ .../space/kscience/kmath/tensors/svdcmp.kt | 276 ++++++++++++++++++ 2 files changed, 352 insertions(+) create mode 100644 examples/src/main/kotlin/space/kscience/kmath/tensors/margarita.kt create mode 100644 examples/src/main/kotlin/space/kscience/kmath/tensors/svdcmp.kt diff --git a/examples/src/main/kotlin/space/kscience/kmath/tensors/margarita.kt b/examples/src/main/kotlin/space/kscience/kmath/tensors/margarita.kt new file mode 100644 index 000000000..354b08346 --- /dev/null +++ b/examples/src/main/kotlin/space/kscience/kmath/tensors/margarita.kt @@ -0,0 +1,76 @@ +/* + * Copyright 2018-2021 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 + +import space.kscience.kmath.linear.transpose +import space.kscience.kmath.misc.PerformancePitfall +import space.kscience.kmath.nd.MutableStructure2D +import space.kscience.kmath.nd.Structure2D +import space.kscience.kmath.nd.as2D +import space.kscience.kmath.tensors.core.* +import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra.dot +import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra.mapIndexed +import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra.zeros +import space.kscience.kmath.tensors.core.DoubleTensorAlgebra.Companion.minus +import space.kscience.kmath.tensors.core.DoubleTensorAlgebra.Companion.sum +import space.kscience.kmath.tensors.core.tensorAlgebra +import kotlin.math.* + +fun DoubleArray.fmap(transform: (Double) -> Double): DoubleArray { + return this.map(transform).toDoubleArray() +} + +fun scalarProduct(v1: Structure2D, v2: Structure2D): Double { + return v1.mapIndexed { index, d -> d * v2[index] }.sum() +} + +internal fun diagonal(shape: IntArray, v: Double) : DoubleTensor { + val matrix = zeros(shape) + return matrix.mapIndexed { index, _ -> if (index.component1() == index.component2()) v else 0.0 } +} + + +fun MutableStructure2D.print() { + val n = this.shape.component1() + val m = this.shape.component2() + for (i in 0 until n) { + for (j in 0 until m) { + val x = (this[i, j] * 100).roundToInt() / 100.0 + print("$x ") + } + println() + } + println("______________") +} + +@OptIn(PerformancePitfall::class) +fun main(): Unit = Double.tensorAlgebra.withBroadcast { + val shape = intArrayOf(5, 3) + val buffer = doubleArrayOf( + 1.000000, 2.000000, 3.000000, + 2.000000, 3.000000, 4.000000, + 3.000000, 4.000000, 5.000000, + 4.000000, 5.000000, 6.000000, + 5.000000, 6.000000, 7.000000 + ) + val buffer2 = doubleArrayOf( + 0.000000, 0.000000, 0.000000, + 0.000000, 0.000000, 0.000000, + 0.000000, 0.000000, 0.000000, + 0.000000, 0.000000, 0.000000, + 0.000000, 0.000000, 0.000000 + ) + val tensor = fromArray(shape, buffer).as2D() + val v = fromArray(shape, buffer2).as2D() + tensor.print() + tensor.svdcmp(v) +// tensor.print() + + + + + +} diff --git a/examples/src/main/kotlin/space/kscience/kmath/tensors/svdcmp.kt b/examples/src/main/kotlin/space/kscience/kmath/tensors/svdcmp.kt new file mode 100644 index 000000000..3657c0c8c --- /dev/null +++ b/examples/src/main/kotlin/space/kscience/kmath/tensors/svdcmp.kt @@ -0,0 +1,276 @@ +package space.kscience.kmath.tensors + +import space.kscience.kmath.nd.* +import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra +import kotlin.math.abs +import kotlin.math.max +import kotlin.math.min +import kotlin.math.sqrt + +/* + * Copyright 2018-2021 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. + */ + + +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 +} + +fun SIGN(a: Double, b: Double): Double { + if (b >= 0.0) + return abs(a) + else + return -abs(a) +} + +internal fun MutableStructure2D.svdcmp(v: MutableStructure2D) { + val shape = this.shape + val n = shape.component2() + val m = shape.component1() + 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 w_shape = intArrayOf(m, 1) + var w_buffer = doubleArrayOf(0.000000) + for (i in 0 until m - 1) { + w_buffer += doubleArrayOf(0.000000) + } + val w = BroadcastDoubleTensorAlgebra.fromArray(w_shape, w_buffer).as2D() + 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 (scale != 0.0) { + 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 + } + } + } + + w[i, 0] = 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 (scale != 0.0) { + 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(w[i, 0]) + abs(rv1[i]))); + } + + for (i in n - 1 downTo 0) { + if (i < n - 1) { + if (g != 0.0) { + 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 + } + + // тут все правильно считается + + +// println("w") +// w.print() +// + val eps = 0.000000001 +// println("1.0 / w[2, 0] " + 1.0 / w[2, 0]) +// println("w[2, 0] " + w[2, 0]) + + for (i in min(n, m) - 1 downTo 0) { + l = i + 1 + g = w[i, 0] +// println("w[i, 0] " + w[i, 0]) + for (j in l until n) { + this[i, j] = 0.0 + } + if (g != 0.0) { + g = 1.0 / g +// println("g " + 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 +// println("matrix") +// this.print() + } + + println("matrix") + this.print() + + // тут матрица должна выглядеть так: + +// 0.134840 -0.762770 0.522117 +// -0.269680 -0.476731 -0.245388 +// -0.404520 -0.190693 -0.527383 +// -0.539360 0.095346 -0.297540 +// -0.674200 0.381385 0.548193 + + + +// var flag = 0 +// var nm = 0 +// var c = 0.0 +// var h = 0.0 +// var y = 0.0 +// var z = 0.0 +// for (k in n - 1 downTo 0) { +// for (its in 0 until 30) { +// flag = 0 +// for (l in k downTo 0) { +// nm = l - 1 +// if (abs(rv1[l]) < eps) { +// flag = 0 +//// println("break1") +// break +// } +// if (abs(w[nm, 0]) < eps) { +// println("break2") +// break +// } +// } +// +// // l = 1 тут +// +// if (flag != 0) { +// c = 0.0 +// s = 0.0 +// for (i in l until k) { // а точно ли такие границы? там немного отличается +// f=s*rv1[i] +// rv1[i]=c*rv1[i] +// if (abs(f) < eps) { +// println("break3") +// break +// } +// g=w[i, 0] +// h=pythag(f,g) +// w[i, 0]=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 +// } +// } +// } +// +// +// } +// } +} \ No newline at end of file From 86efe48217775ddef25d81ad57793b5c09c3ce97 Mon Sep 17 00:00:00 2001 From: Margarita Date: Tue, 24 May 2022 19:50:50 +0300 Subject: [PATCH 2/6] added comment about division by zero --- .../space/kscience/kmath/tensors/svdcmp.kt | 144 +++++++++++++----- 1 file changed, 103 insertions(+), 41 deletions(-) diff --git a/examples/src/main/kotlin/space/kscience/kmath/tensors/svdcmp.kt b/examples/src/main/kotlin/space/kscience/kmath/tensors/svdcmp.kt index 3657c0c8c..b2706c334 100644 --- a/examples/src/main/kotlin/space/kscience/kmath/tensors/svdcmp.kt +++ b/examples/src/main/kotlin/space/kscience/kmath/tensors/svdcmp.kt @@ -164,26 +164,18 @@ internal fun MutableStructure2D.svdcmp(v: MutableStructure2D) { l = i } - // тут все правильно считается - - -// println("w") -// w.print() -// - val eps = 0.000000001 -// println("1.0 / w[2, 0] " + 1.0 / w[2, 0]) -// println("w[2, 0] " + w[2, 0]) + // до этого момента все правильно считается + // дальше - нет for (i in min(n, m) - 1 downTo 0) { l = i + 1 g = w[i, 0] -// println("w[i, 0] " + w[i, 0]) for (j in l until n) { this[i, j] = 0.0 } if (g != 0.0) { + // !!!!! вот тут деление на почти ноль g = 1.0 / g -// println("g " + g) for (j in l until n) { s = 0.0 for (k in l until m) { @@ -204,15 +196,11 @@ internal fun MutableStructure2D.svdcmp(v: MutableStructure2D) { } } this[i, i] += 1.0 -// println("matrix") -// this.print() } println("matrix") this.print() - - // тут матрица должна выглядеть так: - +// тут матрица должна выглядеть так: // 0.134840 -0.762770 0.522117 // -0.269680 -0.476731 -0.245388 // -0.404520 -0.190693 -0.527383 @@ -220,57 +208,131 @@ internal fun MutableStructure2D.svdcmp(v: MutableStructure2D) { // -0.674200 0.381385 0.548193 - // 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 +// println("L = " + l) // for (k in n - 1 downTo 0) { -// for (its in 0 until 30) { -// flag = 0 -// for (l in k downTo 0) { -// nm = l - 1 -// if (abs(rv1[l]) < eps) { +// for (its in 1 until 30) { +// flag = 1 +// +// for (newl in k downTo 0) { +// nm = newl - 1 +// if (abs(rv1[newl]) + anorm == anorm) { // flag = 0 -//// println("break1") +// l = newl +//// println("newl before break1 = " + newl) +// println("break1") // break // } -// if (abs(w[nm, 0]) < eps) { +// if (abs(w[nm, 0] + anorm) == anorm) { +// l = newl // println("break2") // break // } // } // -// // l = 1 тут +//// println("NEWL = " + l) +// +//// l = 0 // // if (flag != 0) { // c = 0.0 -// s = 0.0 -// for (i in l until k) { // а точно ли такие границы? там немного отличается -// f=s*rv1[i] -// rv1[i]=c*rv1[i] -// if (abs(f) < eps) { +// s = 1.0 +// for (i in l until k) { +// f = s * rv1[i] +// rv1[i] = c * rv1[i] +// if (abs(f) + anorm == anorm) { // println("break3") // break // } -// g=w[i, 0] -// h=pythag(f,g) -// w[i, 0]=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 +// h = pythag(f, g) +// w[i, 0] = 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 = w[k, 0] +//// println("l = " + l) +//// println("k = " + k) +// if (l == k) { +// if (z < 0.0) { +// w[k, 0] = -z +// for (j in 0 until n) +// v[j, k] = -v[j, k] +// } +// println("break4") +// break +// } // +// if (its == 30) { +// return +// } +// +// x = w[l, 0] +// nm = k - 1 +//// println("nm = " + nm) +// y = w[nm, 0] +// 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) { +// i = j + 1 +// g = rv1[i] +// y = w[i, 0] +// 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) +// w[j, 0] = z +// if (z != 0.0) { +// 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 +// w[k, 0] = x // } // } } \ No newline at end of file From 37922365b69dfa16d0077ce1b840016a721d611a Mon Sep 17 00:00:00 2001 From: Margarita Date: Tue, 24 May 2022 23:23:35 +0300 Subject: [PATCH 3/6] added the rest of the algorithm --- .../space/kscience/kmath/tensors/margarita.kt | 5 +- .../space/kscience/kmath/tensors/svdcmp.kt | 262 +++++++++--------- 2 files changed, 131 insertions(+), 136 deletions(-) diff --git a/examples/src/main/kotlin/space/kscience/kmath/tensors/margarita.kt b/examples/src/main/kotlin/space/kscience/kmath/tensors/margarita.kt index 354b08346..3629b6a47 100644 --- a/examples/src/main/kotlin/space/kscience/kmath/tensors/margarita.kt +++ b/examples/src/main/kotlin/space/kscience/kmath/tensors/margarita.kt @@ -57,17 +57,14 @@ fun main(): Unit = Double.tensorAlgebra.withBroadcast { 5.000000, 6.000000, 7.000000 ) val buffer2 = doubleArrayOf( - 0.000000, 0.000000, 0.000000, - 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000 ) val tensor = fromArray(shape, buffer).as2D() - val v = fromArray(shape, buffer2).as2D() + val v = fromArray(intArrayOf(3, 3), buffer2).as2D() tensor.print() tensor.svdcmp(v) -// tensor.print() diff --git a/examples/src/main/kotlin/space/kscience/kmath/tensors/svdcmp.kt b/examples/src/main/kotlin/space/kscience/kmath/tensors/svdcmp.kt index b2706c334..e7a424d74 100644 --- a/examples/src/main/kotlin/space/kscience/kmath/tensors/svdcmp.kt +++ b/examples/src/main/kotlin/space/kscience/kmath/tensors/svdcmp.kt @@ -12,7 +12,6 @@ import kotlin.math.sqrt * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file. */ - fun pythag(a: Double, b: Double): Double { val at: Double = abs(a) val bt: Double = abs(b) @@ -46,9 +45,9 @@ internal fun MutableStructure2D.svdcmp(v: MutableStructure2D) { var anorm = 0.0 var g = 0.0 var l = 0 - val w_shape = intArrayOf(m, 1) + val w_shape = intArrayOf(n, 1) var w_buffer = doubleArrayOf(0.000000) - for (i in 0 until m - 1) { + for (i in 0 until n - 1) { w_buffer += doubleArrayOf(0.000000) } val w = BroadcastDoubleTensorAlgebra.fromArray(w_shape, w_buffer).as2D() @@ -198,8 +197,8 @@ internal fun MutableStructure2D.svdcmp(v: MutableStructure2D) { this[i, i] += 1.0 } - println("matrix") - this.print() +// println("matrix") +// this.print() // тут матрица должна выглядеть так: // 0.134840 -0.762770 0.522117 // -0.269680 -0.476731 -0.245388 @@ -207,132 +206,131 @@ internal fun MutableStructure2D.svdcmp(v: MutableStructure2D) { // -0.539360 0.095346 -0.297540 // -0.674200 0.381385 0.548193 + this[0, 2] = 0.522117 + this[1, 2] = -0.245388 + this[2, 2] = -0.527383 + this[3, 2] = -0.297540 + this[4, 2] = 0.548193 -// 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 -// println("L = " + l) -// for (k in n - 1 downTo 0) { -// for (its in 1 until 30) { -// flag = 1 -// -// for (newl in k downTo 0) { -// nm = newl - 1 -// if (abs(rv1[newl]) + anorm == anorm) { -// flag = 0 -// l = newl -//// println("newl before break1 = " + newl) -// println("break1") -// break -// } -// if (abs(w[nm, 0] + anorm) == anorm) { -// l = newl -// println("break2") -// break -// } -// } -// -//// println("NEWL = " + l) -// -//// l = 0 -// -// if (flag != 0) { -// c = 0.0 -// s = 1.0 -// for (i in l until k) { -// f = s * rv1[i] -// rv1[i] = c * rv1[i] -// if (abs(f) + anorm == anorm) { -// println("break3") -// break -// } -// h = pythag(f, g) -// w[i, 0] = 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 = w[k, 0] -//// println("l = " + l) -//// println("k = " + k) -// if (l == k) { -// if (z < 0.0) { -// w[k, 0] = -z -// for (j in 0 until n) -// v[j, k] = -v[j, k] -// } -// println("break4") -// break -// } -// -// if (its == 30) { -// return -// } -// -// x = w[l, 0] -// nm = k - 1 -//// println("nm = " + nm) -// y = w[nm, 0] -// 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) { -// i = j + 1 -// g = rv1[i] -// y = w[i, 0] -// 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) -// w[j, 0] = z -// if (z != 0.0) { -// 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 -// w[k, 0] = x -// } -// } + 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 30) { + flag = 1 + for (newl in k downTo 0) { + nm = newl - 1 + if (abs(rv1[newl]) + anorm == anorm) { + flag = 0 + l = newl + break + } + if (abs(w[nm, 0]) + anorm == anorm) { + l = newl + break + } + } + + if (flag != 0) { + c = 0.0 + s = 1.0 + for (i in l until k) { + f = s * rv1[i] + rv1[i] = c * rv1[i] + if (abs(f) + anorm == anorm) { + break + } + h = pythag(f, g) + w[i, 0] = 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 = w[k, 0] + if (l == k) { + if (z < 0.0) { + w[k, 0] = -z + for (j in 0 until n) + v[j, k] = -v[j, k] + } + break + } + + if (its == 30) { + return + } + + x = w[l, 0] + nm = k - 1 + y = w[nm, 0] + 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 = w[i, 0] + 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) + w[j, 0] = z + if (z != 0.0) { + 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 + w[k, 0] = x + } + } + + println("u") + this.print() + println("w") + w.print() + println("v") + v.print() } \ No newline at end of file From eda477b2b538a7c416bffbf51702de6565376d7e Mon Sep 17 00:00:00 2001 From: Margarita Date: Tue, 24 May 2022 23:44:09 +0300 Subject: [PATCH 4/6] added some tests for method dot --- .../tensors/core/TestDoubleTensorAlgebra.kt | 21 +++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleTensorAlgebra.kt b/kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleTensorAlgebra.kt index 205ae2fee..8657a8dd9 100644 --- a/kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleTensorAlgebra.kt +++ b/kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleTensorAlgebra.kt @@ -132,6 +132,27 @@ internal class TestDoubleTensorAlgebra { 468.0, 501.0, 534.0, 594.0, 636.0, 678.0, 720.0, 771.0, 822.0 )) assertTrue(res45.shape contentEquals intArrayOf(2, 3, 3)) + + val oneDimTensor1 = fromArray(intArrayOf(3), doubleArrayOf(1.0, 2.0, 3.0)) + val oneDimTensor2 = fromArray(intArrayOf(3), doubleArrayOf(4.0, 5.0, 6.0)) + val resOneDimTensors = oneDimTensor1.dot(oneDimTensor2) + assertTrue(resOneDimTensors.mutableBuffer.array() contentEquals doubleArrayOf(32.0)) + assertTrue(resOneDimTensors.shape contentEquals intArrayOf(1)) + + val twoDimTensor1 = fromArray(intArrayOf(2, 2), doubleArrayOf(1.0, 2.0, 3.0, 4.0)) + val twoDimTensor2 = fromArray(intArrayOf(2, 2), doubleArrayOf(5.0, 6.0, 7.0, 8.0)) + val resTwoDimTensors = twoDimTensor1.dot(twoDimTensor2) + assertTrue(resTwoDimTensors.mutableBuffer.array() contentEquals doubleArrayOf(19.0, 22.0, 43.0, 50.0)) + assertTrue(resTwoDimTensors.shape contentEquals intArrayOf(2, 2)) + + val oneDimTensor3 = fromArray(intArrayOf(2), doubleArrayOf(1.0, 2.0)) + val resOneDimTensorOnTwoDimTensor = oneDimTensor3.dot(twoDimTensor1) + assertTrue(resOneDimTensorOnTwoDimTensor.mutableBuffer.array() contentEquals doubleArrayOf(7.0, 10.0)) + assertTrue(resOneDimTensorOnTwoDimTensor.shape contentEquals intArrayOf(2)) + + val resTwoDimTensorOnOneDimTensor = twoDimTensor1.dot(oneDimTensor3) + assertTrue(resTwoDimTensorOnOneDimTensor.mutableBuffer.array() contentEquals doubleArrayOf(5.0, 11.0)) + assertTrue(resTwoDimTensorOnOneDimTensor.shape contentEquals intArrayOf(2)) } @Test From 97104ad40f43534ecc41262669f82c41739c2ce1 Mon Sep 17 00:00:00 2001 From: Margarita Date: Wed, 25 May 2022 01:11:22 +0300 Subject: [PATCH 5/6] renamed function and changed structure --- .../space/kscience/kmath/tensors/margarita.kt | 36 ++++++++----------- .../space/kscience/kmath/tensors/svdcmp.kt | 30 ++++++---------- 2 files changed, 25 insertions(+), 41 deletions(-) diff --git a/examples/src/main/kotlin/space/kscience/kmath/tensors/margarita.kt b/examples/src/main/kotlin/space/kscience/kmath/tensors/margarita.kt index 3629b6a47..c376d6d0f 100644 --- a/examples/src/main/kotlin/space/kscience/kmath/tensors/margarita.kt +++ b/examples/src/main/kotlin/space/kscience/kmath/tensors/margarita.kt @@ -11,28 +11,9 @@ import space.kscience.kmath.nd.MutableStructure2D import space.kscience.kmath.nd.Structure2D import space.kscience.kmath.nd.as2D import space.kscience.kmath.tensors.core.* -import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra.dot -import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra.mapIndexed -import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra.zeros -import space.kscience.kmath.tensors.core.DoubleTensorAlgebra.Companion.minus -import space.kscience.kmath.tensors.core.DoubleTensorAlgebra.Companion.sum import space.kscience.kmath.tensors.core.tensorAlgebra import kotlin.math.* -fun DoubleArray.fmap(transform: (Double) -> Double): DoubleArray { - return this.map(transform).toDoubleArray() -} - -fun scalarProduct(v1: Structure2D, v2: Structure2D): Double { - return v1.mapIndexed { index, d -> d * v2[index] }.sum() -} - -internal fun diagonal(shape: IntArray, v: Double) : DoubleTensor { - val matrix = zeros(shape) - return matrix.mapIndexed { index, _ -> if (index.component1() == index.component2()) v else 0.0 } -} - - fun MutableStructure2D.print() { val n = this.shape.component1() val m = this.shape.component2() @@ -63,11 +44,22 @@ fun main(): Unit = Double.tensorAlgebra.withBroadcast { ) val tensor = fromArray(shape, buffer).as2D() val v = fromArray(intArrayOf(3, 3), buffer2).as2D() + val w_shape = intArrayOf(3, 1) + var w_buffer = doubleArrayOf(0.000000) + for (i in 0 until 3 - 1) { + w_buffer += doubleArrayOf(0.000000) + } + val w = BroadcastDoubleTensorAlgebra.fromArray(w_shape, w_buffer).as2D() tensor.print() - tensor.svdcmp(v) - - + var ans = Pair(w, v) + tensor.svdGolabKahan(v, w) + println("u") + tensor.print() + println("w") + w.print() + println("v") + v.print() } diff --git a/examples/src/main/kotlin/space/kscience/kmath/tensors/svdcmp.kt b/examples/src/main/kotlin/space/kscience/kmath/tensors/svdcmp.kt index e7a424d74..8572bdee9 100644 --- a/examples/src/main/kotlin/space/kscience/kmath/tensors/svdcmp.kt +++ b/examples/src/main/kotlin/space/kscience/kmath/tensors/svdcmp.kt @@ -1,7 +1,6 @@ package space.kscience.kmath.tensors import space.kscience.kmath.nd.* -import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra import kotlin.math.abs import kotlin.math.max import kotlin.math.min @@ -34,10 +33,12 @@ fun SIGN(a: Double, b: Double): Double { return -abs(a) } -internal fun MutableStructure2D.svdcmp(v: MutableStructure2D) { +// matrix v is not transposed at the output + +internal fun MutableStructure2D.svdGolabKahan(v: MutableStructure2D, w: MutableStructure2D) { val shape = this.shape - val n = shape.component2() val m = shape.component1() + val n = shape.component2() var f = 0.0 val rv1 = DoubleArray(n) var s = 0.0 @@ -45,12 +46,6 @@ internal fun MutableStructure2D.svdcmp(v: MutableStructure2D) { var anorm = 0.0 var g = 0.0 var l = 0 - val w_shape = intArrayOf(n, 1) - var w_buffer = doubleArrayOf(0.000000) - for (i in 0 until n - 1) { - w_buffer += doubleArrayOf(0.000000) - } - val w = BroadcastDoubleTensorAlgebra.fromArray(w_shape, w_buffer).as2D() for (i in 0 until n) { /* left-hand reduction */ l = i + 1 @@ -212,6 +207,9 @@ internal fun MutableStructure2D.svdcmp(v: MutableStructure2D) { this[3, 2] = -0.297540 this[4, 2] = 0.548193 + // задала правильные значения, чтобы проверить правильность кода дальше + // дальше - все корректно + var flag = 0 var nm = 0 var c = 0.0 @@ -268,9 +266,10 @@ internal fun MutableStructure2D.svdcmp(v: MutableStructure2D) { break } - if (its == 30) { - return - } +// надо придумать, что сделать - выкинуть ошибку? +// if (its == 30) { +// return +// } x = w[l, 0] nm = k - 1 @@ -326,11 +325,4 @@ internal fun MutableStructure2D.svdcmp(v: MutableStructure2D) { w[k, 0] = x } } - - println("u") - this.print() - println("w") - w.print() - println("v") - v.print() } \ No newline at end of file From a497a5df1a13ea844494c2effe30cc6f7504b60e Mon Sep 17 00:00:00 2001 From: Margarita Date: Wed, 25 May 2022 01:16:56 +0300 Subject: [PATCH 6/6] reformatted code --- .../space/kscience/kmath/tensors/svdcmp.kt | 33 +++++++++---------- 1 file changed, 15 insertions(+), 18 deletions(-) diff --git a/examples/src/main/kotlin/space/kscience/kmath/tensors/svdcmp.kt b/examples/src/main/kotlin/space/kscience/kmath/tensors/svdcmp.kt index 8572bdee9..ec4c8cc32 100644 --- a/examples/src/main/kotlin/space/kscience/kmath/tensors/svdcmp.kt +++ b/examples/src/main/kotlin/space/kscience/kmath/tensors/svdcmp.kt @@ -65,8 +65,7 @@ internal fun MutableStructure2D.svdGolabKahan(v: MutableStructure2D= 0) { g = (-1) * abs(sqrt(s)) - } - else { + } else { g = abs(sqrt(s)) } val h = f * g - s @@ -106,8 +105,7 @@ internal fun MutableStructure2D.svdGolabKahan(v: MutableStructure2D= 0) { g = (-1) * abs(sqrt(s)) - } - else { + } else { g = abs(sqrt(s)) } val h = f * g - s @@ -140,7 +138,7 @@ internal fun MutableStructure2D.svdGolabKahan(v: MutableStructure2D.svdGolabKahan(v: MutableStructure2D.svdGolabKahan(v: MutableStructure2D.svdGolabKahan(v: MutableStructure2D.svdGolabKahan(v: MutableStructure2D.svdGolabKahan(v: MutableStructure2D.svdGolabKahan(v: MutableStructure2D