From fea53af0eeaa3522a5754bfe83688f010a37a5f4 Mon Sep 17 00:00:00 2001 From: Andrei Kislitsyn Date: Wed, 24 Mar 2021 14:00:47 +0300 Subject: [PATCH 1/3] fix --- .../space/kscience/kmath/tensors/LinearOpsTensorAlgebra.kt | 2 +- .../kmath/tensors/core/DoubleLinearOpsTensorAlgebra.kt | 3 +-- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/LinearOpsTensorAlgebra.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/LinearOpsTensorAlgebra.kt index c4fa3e0a8..e412af7c6 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/LinearOpsTensorAlgebra.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/LinearOpsTensorAlgebra.kt @@ -14,7 +14,7 @@ public interface LinearOpsTensorAlgebra, Inde public fun TensorType.qr(): TensorType //https://pytorch.org/docs/stable/generated/torch.lu.html - public fun TensorType.lu(tol: T): Pair + public fun TensorType.lu(): Pair //https://pytorch.org/docs/stable/generated/torch.lu_unpack.html public fun luPivot(luTensor: TensorType, pivotsTensor: IndexTensorType): Triple diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleLinearOpsTensorAlgebra.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleLinearOpsTensorAlgebra.kt index e0abc49b7..d090ce79c 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleLinearOpsTensorAlgebra.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleLinearOpsTensorAlgebra.kt @@ -11,7 +11,7 @@ public class DoubleLinearOpsTensorAlgebra : TODO("Not yet implemented") } - override fun DoubleTensor.lu(tol: Double): Pair { + override fun DoubleTensor.lu(): Pair { checkSquareMatrix(shape) @@ -138,7 +138,6 @@ public class DoubleLinearOpsTensorAlgebra : TODO("Not yet implemented") } - override fun DoubleTensor.svd(): Triple { TODO("Not yet implemented") } -- 2.34.1 From 206bcfc909d56a859ec783a30cc75a8ade178786 Mon Sep 17 00:00:00 2001 From: Andrei Kislitsyn Date: Wed, 24 Mar 2021 18:08:36 +0300 Subject: [PATCH 2/3] lu inv and det complete + tests --- .../kmath/tensors/LinearOpsTensorAlgebra.kt | 2 +- .../core/DoubleLinearOpsTensorAlgebra.kt | 79 +++++++++++++++++++ .../kmath/tensors/core/DoubleTensorAlgebra.kt | 35 +++++++- .../core/TestDoubleAnalyticTensorAlgebra.kt | 7 +- .../core/TestDoubleLinearOpsAlgebra.kt | 73 +++++++++++++++++ 5 files changed, 187 insertions(+), 9 deletions(-) create mode 100644 kmath-core/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleLinearOpsAlgebra.kt diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/LinearOpsTensorAlgebra.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/LinearOpsTensorAlgebra.kt index e412af7c6..f551d524a 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/LinearOpsTensorAlgebra.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/LinearOpsTensorAlgebra.kt @@ -13,7 +13,7 @@ public interface LinearOpsTensorAlgebra, Inde //https://pytorch.org/docs/stable/linalg.html#torch.linalg.qr public fun TensorType.qr(): TensorType - //https://pytorch.org/docs/stable/generated/torch.lu.html + //htt ps://pytorch.org/docs/stable/generated/torch.lu.html public fun TensorType.lu(): Pair //https://pytorch.org/docs/stable/generated/torch.lu_unpack.html diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleLinearOpsTensorAlgebra.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleLinearOpsTensorAlgebra.kt index 7779ad029..c8c4c413f 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleLinearOpsTensorAlgebra.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleLinearOpsTensorAlgebra.kt @@ -1,5 +1,8 @@ package space.kscience.kmath.tensors.core +import space.kscience.kmath.nd.MutableStructure2D +import space.kscience.kmath.nd.Structure1D +import space.kscience.kmath.nd.Structure2D import space.kscience.kmath.tensors.LinearOpsTensorAlgebra import kotlin.math.sqrt @@ -20,6 +23,8 @@ public class DoubleLinearOpsTensorAlgebra : val n = shape.size val m = shape.last() val pivotsShape = IntArray(n - 1) { i -> shape[i] } + pivotsShape[n - 2] = m + 1 + val pivotsTensor = IntTensor( pivotsShape, IntArray(pivotsShape.reduce(Int::times)) { 0 } @@ -54,6 +59,8 @@ public class DoubleLinearOpsTensorAlgebra : lu[maxInd, k] = tmp } + pivots[m] += 1 + } for (j in i + 1 until m) { @@ -146,6 +153,78 @@ public class DoubleLinearOpsTensorAlgebra : TODO("ANDREI") } + private fun luMatrixDet(lu: Structure2D, pivots: Structure1D): Double { + val m = lu.shape[0] + val sign = if((pivots[m] - m) % 2 == 0) 1.0 else -1.0 + var det = sign + for (i in 0 until m){ + det *= lu[i, i] + } + return det + } + + public fun DoubleTensor.detLU(): DoubleTensor { + val (luTensor, pivotsTensor) = this.lu() + val n = shape.size + + val detTensorShape = IntArray(n - 1) { i -> shape[i] } + detTensorShape[n - 2] = 1 + val resBuffer = DoubleArray(detTensorShape.reduce(Int::times)) { 0.0 } + + val detTensor = DoubleTensor( + detTensorShape, + resBuffer + ) + + luTensor.matrixSequence().zip(pivotsTensor.vectorSequence()).forEachIndexed { index, (luMatrix, pivots) -> + resBuffer[index] = luMatrixDet(luMatrix, pivots) + } + + return detTensor + } + + private fun luMatrixInv( + lu: Structure2D, + pivots: Structure1D, + invMatrix : MutableStructure2D + ): Unit { + val m = lu.shape[0] + + for (j in 0 until m) { + for (i in 0 until m) { + if (pivots[i] == j){ + invMatrix[i, j] = 1.0 + } + + for (k in 0 until i){ + invMatrix[i, j] -= lu[i, k] * invMatrix[k, j] + } + } + + for (i in m - 1 downTo 0) { + for (k in i + 1 until m) { + invMatrix[i, j] -= lu[i, k] * invMatrix[k, j] + } + invMatrix[i, j] /= lu[i, i] + } + } + } + + public fun DoubleTensor.invLU(): DoubleTensor { + val (luTensor, pivotsTensor) = this.lu() + val n = shape.size + val invTensor = luTensor.zeroesLike() + + for ( + (luP, invMatrix) in + luTensor.matrixSequence().zip(pivotsTensor.vectorSequence()).zip(invTensor.matrixSequence()) + ) { + val (lu, pivots) = luP + luMatrixInv(lu, pivots, invMatrix) + } + + return invTensor + } } public inline fun DoubleLinearOpsTensorAlgebra(block: DoubleLinearOpsTensorAlgebra.() -> R): R = diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleTensorAlgebra.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleTensorAlgebra.kt index c67687a09..ff65ca027 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleTensorAlgebra.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleTensorAlgebra.kt @@ -1,5 +1,8 @@ package space.kscience.kmath.tensors.core +import space.kscience.kmath.linear.Matrix +import space.kscience.kmath.nd.MutableStructure2D +import space.kscience.kmath.nd.Structure2D import space.kscience.kmath.tensors.TensorPartialDivisionAlgebra import kotlin.math.abs @@ -262,7 +265,31 @@ public open class DoubleTensorAlgebra : TensorPartialDivisionAlgebra shape[i] } + detTensorShape[n - 1] = 1 + val resBuffer = DoubleArray(detTensorShape.reduce(Int::times)) { 0.0 } + + val detTensor = DoubleTensor( + detTensorShape, + resBuffer + ) + + this.matrixSequence().forEachIndexed{i, matrix -> + // todo need Matrix determinant algo + // todo resBuffer[i] = matrix.det() + } + + + return detTensor + + */ } override fun DoubleTensor.square(): DoubleTensor { @@ -294,7 +321,7 @@ public open class DoubleTensorAlgebra : TensorPartialDivisionAlgebra Boolean): Boolean { - if (!(this.shape contentEquals other.shape)){ + if (!(this.shape contentEquals other.shape)) { return false } return this.eq(other, eqFunction) @@ -303,10 +330,10 @@ public open class DoubleTensorAlgebra : TensorPartialDivisionAlgebra Boolean): Boolean { // todo broadcasting checking val n = this.linearStructure.size - if (n != other.linearStructure.size){ + if (n != other.linearStructure.size) { return false } - for (i in 0 until n){ + for (i in 0 until n) { if (!eqFunction(this.buffer[this.bufferStart + i], other.buffer[other.bufferStart + i])) { return false } diff --git a/kmath-core/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleAnalyticTensorAlgebra.kt b/kmath-core/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleAnalyticTensorAlgebra.kt index 8028ce175..4dcdb7848 100644 --- a/kmath-core/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleAnalyticTensorAlgebra.kt +++ b/kmath-core/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleAnalyticTensorAlgebra.kt @@ -3,7 +3,6 @@ package space.kscience.kmath.tensors.core import kotlin.math.abs import kotlin.math.exp import kotlin.test.Test -import kotlin.test.assertEquals import kotlin.test.assertTrue class TestDoubleAnalyticTensorAlgebra { @@ -16,9 +15,9 @@ class TestDoubleAnalyticTensorAlgebra { return this.map(transform).toDoubleArray() } - fun DoubleArray.deltaEqual(other: DoubleArray, delta: Double = 1e-5): Boolean { + fun DoubleArray.epsEqual(other: DoubleArray, eps: Double = 1e-5): Boolean { for ((elem1, elem2) in this.asSequence().zip(other.asSequence())) { - if (abs(elem1 - elem2) > delta) { + if (abs(elem1 - elem2) > eps) { return false } } @@ -29,7 +28,7 @@ class TestDoubleAnalyticTensorAlgebra { fun testExp() = DoubleAnalyticTensorAlgebra { tensor.exp().let { assertTrue { shape contentEquals it.shape } - assertTrue { buffer.fmap(::exp).deltaEqual(it.buffer.array())} + assertTrue { buffer.fmap(::exp).epsEqual(it.buffer.array())} } } } \ No newline at end of file diff --git a/kmath-core/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleLinearOpsAlgebra.kt b/kmath-core/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleLinearOpsAlgebra.kt new file mode 100644 index 000000000..2da79382f --- /dev/null +++ b/kmath-core/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleLinearOpsAlgebra.kt @@ -0,0 +1,73 @@ +package space.kscience.kmath.tensors.core + +import kotlin.math.abs +import kotlin.math.exp +import kotlin.test.Test +import kotlin.test.assertEquals +import kotlin.test.assertTrue + +class TestDoubleLinearOpsTensorAlgebra { + + private val eps = 1e-5 + + private fun Double.epsEqual(other: Double): Boolean { + return abs(this - other) < eps + } + + fun DoubleArray.epsEqual(other: DoubleArray, eps: Double = 1e-5): Boolean { + for ((elem1, elem2) in this.asSequence().zip(other.asSequence())) { + if (abs(elem1 - elem2) > eps) { + return false + } + } + return true + } + + @Test + fun testDetLU() = DoubleLinearOpsTensorAlgebra { + val tensor = fromArray( + intArrayOf(2, 2, 2), + doubleArrayOf( + 1.0, 3.0, + 1.0, 2.0, + 1.5, 1.0, + 10.0, 2.0 + ) + ) + + val expectedShape = intArrayOf(2, 1) + val expectedBuffer = doubleArrayOf( + -1.0, + -7.0 + ) + val detTensor = tensor.detLU() + + assertTrue { detTensor.shape contentEquals expectedShape } + assertTrue { detTensor.buffer.array().epsEqual(expectedBuffer) } + } + + @Test + fun testInvLU() = DoubleLinearOpsTensorAlgebra { + val tensor = fromArray( + intArrayOf(2, 2, 2), + doubleArrayOf( + 1.0, 0.0, + 0.0, 2.0, + 1.0, 1.0, + 1.0, 0.0 + ) + ) + + val expectedShape = intArrayOf(2, 2, 2) + val expectedBuffer = doubleArrayOf( + 1.0, 0.0, + 0.0, 0.5, + 0.0, 1.0, + 1.0, -1.0 + ) + + val invTensor = tensor.invLU() + assertTrue { invTensor.shape contentEquals expectedShape } + assertTrue { invTensor.buffer.array().epsEqual(expectedBuffer) } + } +} -- 2.34.1 From ab8137000146c52f10b3213c960b03917a0ada4f Mon Sep 17 00:00:00 2001 From: Andrei Kislitsyn Date: Wed, 24 Mar 2021 18:42:41 +0300 Subject: [PATCH 3/3] fixes --- .../kmath/tensors/LinearOpsTensorAlgebra.kt | 2 +- .../core/DoubleLinearOpsTensorAlgebra.kt | 22 +++++++------------ 2 files changed, 9 insertions(+), 15 deletions(-) diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/LinearOpsTensorAlgebra.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/LinearOpsTensorAlgebra.kt index f551d524a..e412af7c6 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/LinearOpsTensorAlgebra.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/LinearOpsTensorAlgebra.kt @@ -13,7 +13,7 @@ public interface LinearOpsTensorAlgebra, Inde //https://pytorch.org/docs/stable/linalg.html#torch.linalg.qr public fun TensorType.qr(): TensorType - //htt ps://pytorch.org/docs/stable/generated/torch.lu.html + //https://pytorch.org/docs/stable/generated/torch.lu.html public fun TensorType.lu(): Pair //https://pytorch.org/docs/stable/generated/torch.lu_unpack.html diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleLinearOpsTensorAlgebra.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleLinearOpsTensorAlgebra.kt index c8c4c413f..0b9d824f0 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleLinearOpsTensorAlgebra.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleLinearOpsTensorAlgebra.kt @@ -10,9 +10,9 @@ public class DoubleLinearOpsTensorAlgebra : LinearOpsTensorAlgebra, DoubleTensorAlgebra() { - override fun DoubleTensor.inv(): DoubleTensor { - TODO("ANDREI") - } + override fun DoubleTensor.inv(): DoubleTensor = invLU() + + override fun DoubleTensor.det(): DoubleTensor = detLU() override fun DoubleTensor.lu(): Pair { @@ -156,15 +156,11 @@ public class DoubleLinearOpsTensorAlgebra : private fun luMatrixDet(lu: Structure2D, pivots: Structure1D): Double { val m = lu.shape[0] val sign = if((pivots[m] - m) % 2 == 0) 1.0 else -1.0 - var det = sign - for (i in 0 until m){ - det *= lu[i, i] - } - return det + return (0 until m).asSequence().map { lu[it, it] }.fold(sign) { left, right -> left * right } } public fun DoubleTensor.detLU(): DoubleTensor { - val (luTensor, pivotsTensor) = this.lu() + val (luTensor, pivotsTensor) = lu() val n = shape.size val detTensorShape = IntArray(n - 1) { i -> shape[i] } @@ -211,14 +207,12 @@ public class DoubleLinearOpsTensorAlgebra : } public fun DoubleTensor.invLU(): DoubleTensor { - val (luTensor, pivotsTensor) = this.lu() + val (luTensor, pivotsTensor) = lu() val n = shape.size val invTensor = luTensor.zeroesLike() - for ( - (luP, invMatrix) in - luTensor.matrixSequence().zip(pivotsTensor.vectorSequence()).zip(invTensor.matrixSequence()) - ) { + val seq = luTensor.matrixSequence().zip(pivotsTensor.vectorSequence()).zip(invTensor.matrixSequence()) + for ((luP, invMatrix) in seq) { val (lu, pivots) = luP luMatrixInv(lu, pivots, invMatrix) } -- 2.34.1