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

added the rest of the algorithm

Browse Source
This commit is contained in:
Margarita 2022-05-24 23:23:35 +03:00
parent 86efe48217
commit 37922365b6
2 changed files with 131 additions and 136 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
examples/src/main/kotlin/space/kscience/kmath/tensors/margarita.kt
View File

@ -57,17 +57,14 @@ fun main(): Unit = Double.tensorAlgebra.withBroadcast {
5.000000, 6.000000, 7.000000 5.000000, 6.000000, 7.000000
) )
val buffer2 = doubleArrayOf( 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, 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 tensor = fromArray(shape, buffer).as2D()
val v = fromArray(shape, buffer2).as2D() val v = fromArray(intArrayOf(3, 3), buffer2).as2D()
tensor.print() tensor.print()
tensor.svdcmp(v) tensor.svdcmp(v)
// tensor.print()

262
examples/src/main/kotlin/space/kscience/kmath/tensors/svdcmp.kt
View File

@ -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. * 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 { fun pythag(a: Double, b: Double): Double {
val at: Double = abs(a) val at: Double = abs(a)
val bt: Double = abs(b) val bt: Double = abs(b)
@ -46,9 +45,9 @@ internal fun MutableStructure2D<Double>.svdcmp(v: MutableStructure2D<Double>) {
var anorm = 0.0 var anorm = 0.0
var g = 0.0 var g = 0.0
var l = 0 var l = 0
val w_shape = intArrayOf(m, 1) val w_shape = intArrayOf(n, 1)
var w_buffer = doubleArrayOf(0.000000) 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) w_buffer += doubleArrayOf(0.000000)
} }
val w = BroadcastDoubleTensorAlgebra.fromArray(w_shape, w_buffer).as2D() val w = BroadcastDoubleTensorAlgebra.fromArray(w_shape, w_buffer).as2D()
@ -198,8 +197,8 @@ internal fun MutableStructure2D<Double>.svdcmp(v: MutableStructure2D<Double>) {
this[i, i] += 1.0 this[i, i] += 1.0
} }
println("matrix") // println("matrix")
this.print() // this.print()
// тут матрица должна выглядеть так: // тут матрица должна выглядеть так:
// 0.134840 -0.762770 0.522117 // 0.134840 -0.762770 0.522117
// -0.269680 -0.476731 -0.245388 // -0.269680 -0.476731 -0.245388
@ -207,132 +206,131 @@ internal fun MutableStructure2D<Double>.svdcmp(v: MutableStructure2D<Double>) {
// -0.539360 0.095346 -0.297540 // -0.539360 0.095346 -0.297540
// -0.674200 0.381385 0.548193 // -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 flag = 0
// var nm = 0 var nm = 0
// var c = 0.0 var c = 0.0
// var h = 0.0 var h = 0.0
// var y = 0.0 var y = 0.0
// var z = 0.0 var z = 0.0
// var x = 0.0 var x = 0.0
// println("L = " + l) for (k in n - 1 downTo 0) {
// for (k in n - 1 downTo 0) { for (its in 1 until 30) {
// for (its in 1 until 30) { flag = 1
// flag = 1 for (newl in k downTo 0) {
// nm = newl - 1
// for (newl in k downTo 0) { if (abs(rv1[newl]) + anorm == anorm) {
// nm = newl - 1 flag = 0
// if (abs(rv1[newl]) + anorm == anorm) { l = newl
// flag = 0 break
// l = newl }
//// println("newl before break1 = " + newl) if (abs(w[nm, 0]) + anorm == anorm) {
// println("break1") l = newl
// break break
// } }
// if (abs(w[nm, 0] + anorm) == anorm) { }
// l = newl
// println("break2") if (flag != 0) {
// break c = 0.0
// } s = 1.0
// } for (i in l until k) {
// f = s * rv1[i]
//// println("NEWL = " + l) rv1[i] = c * rv1[i]
// if (abs(f) + anorm == anorm) {
//// l = 0 break
// }
// if (flag != 0) { h = pythag(f, g)
// c = 0.0 w[i, 0] = h
// s = 1.0 h = 1.0 / h
// for (i in l until k) { c = g * h
// f = s * rv1[i] s = (-f) * h
// rv1[i] = c * rv1[i] for (j in 0 until m) {
// if (abs(f) + anorm == anorm) { y = this[j, nm]
// println("break3") z = this[j, i]
// break 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 z = w[k, 0]
// for (j in 0 until m) { if (l == k) {
// y = this[j, nm] if (z < 0.0) {
// z = this[j, i] w[k, 0] = -z
// this[j, nm] = y * c + z * s for (j in 0 until n)
// this[j, i] = z * c - y * s v[j, k] = -v[j, k]
// } }
// } break
// } }
//
// z = w[k, 0] if (its == 30) {
//// println("l = " + l) return
//// println("k = " + k) }
// if (l == k) {
// if (z < 0.0) { x = w[l, 0]
// w[k, 0] = -z nm = k - 1
// for (j in 0 until n) y = w[nm, 0]
// v[j, k] = -v[j, k] g = rv1[nm]
// } h = rv1[k]
// println("break4") f = ((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y)
// break g = pythag(f,1.0)
// } f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x
// c = 1.0
// if (its == 30) { s = 1.0
// return
// } var i = 0
// for (j in l until nm + 1) {
// x = w[l, 0] i = j + 1
// nm = k - 1 g = rv1[i]
//// println("nm = " + nm) y = w[i, 0]
// y = w[nm, 0] h = s * g
// g = rv1[nm] g = c * g
// h = rv1[k] z = pythag(f,h)
// f = ((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y) rv1[j] = z
// g = pythag(f,1.0) c = f / z
// f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x s = h / z
// c = 1.0 f = x * c + g * s
// s = 1.0 g = g * c - x * s
// h = y * s
// var i = 0 // непонятно, где она должна быть объявлена y *= c
// for (j in l until nm) {
// i = j + 1 for (jj in 0 until n) {
// g = rv1[i] x=v[jj, j];
// y = w[i, 0] z=v[jj, i];
// h = s * g v[jj, j] = x * c + z * s;
// g = c * g v[jj, i] = z * c - x * s;
// z = pythag(f,h) }
// rv1[j] = z z = pythag(f,h)
// c = f / z w[j, 0] = z
// s = h / z if (z != 0.0) {
// f = x * c + g * s z = 1.0 / z
// g = g * c - x * s c = f * z
// h = y * s s = h * z
// y *= c }
// for (jj in 0 until n) { f = c * g + s * y
// x=v[jj, j]; x = c * y - s * g
// z=v[jj, i]; for (jj in 0 until m) {
// v[jj, j] = x * c + z * s; y = this[jj, j]
// v[jj, i] = z * c - x * s; z = this[jj, i]
// } this[jj, j] = y * c + z * s
// z = pythag(f,h) this[jj, i] = z * c - y * s
// w[j, 0] = z }
// if (z != 0.0) { }
// z = 1.0 / z rv1[l] = 0.0
// c = f * z rv1[k] = f
// s = h * z w[k, 0] = x
// } }
// f = c * g + s * y }
// x = c * y - s * g
// for (jj in 0 until m) { println("u")
// y = this[jj, j] this.print()
// z = this[jj, i] println("w")
// this[jj, j] = y * c + z * s w.print()
// this[jj, i] = z * c - y * s println("v")
// } v.print()
// }
// rv1[l] = 0.0
// rv1[k] = f
// w[k, 0] = x
// }
// }
} }