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..c376d6d0f --- /dev/null +++ b/examples/src/main/kotlin/space/kscience/kmath/tensors/margarita.kt @@ -0,0 +1,65 @@ +/* + * 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.tensorAlgebra +import kotlin.math.* + +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 + ) + 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() + 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 new file mode 100644 index 000000000..ec4c8cc32 --- /dev/null +++ b/examples/src/main/kotlin/space/kscience/kmath/tensors/svdcmp.kt @@ -0,0 +1,325 @@ +package space.kscience.kmath.tensors + +import space.kscience.kmath.nd.* +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) +} + +// matrix v is not transposed at the output + +internal fun MutableStructure2D.svdGolabKahan(v: MutableStructure2D, w: MutableStructure2D) { + 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 + 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 + } + + // до этого момента все правильно считается + // дальше - нет + + for (i in min(n, m) - 1 downTo 0) { + l = i + 1 + g = w[i, 0] + for (j in l until n) { + this[i, j] = 0.0 + } + if (g != 0.0) { + // !!!!! вот тут деление на почти ноль + 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 + } + +// 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 + + 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 + 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 + } + } +} \ No newline at end of file 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