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 8a16a991d..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 @@ -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 @@ -7,11 +10,11 @@ public class DoubleLinearOpsTensorAlgebra : LinearOpsTensorAlgebra, DoubleTensorAlgebra() { - override fun DoubleTensor.inv(): DoubleTensor { - TODO("ANDREI") - } + override fun DoubleTensor.inv(): DoubleTensor = invLU() - override fun DoubleTensor.lu(tol: Double): Pair { + override fun DoubleTensor.det(): DoubleTensor = detLU() + + override fun DoubleTensor.lu(): Pair { checkSquareMatrix(shape) @@ -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) { @@ -138,7 +145,6 @@ public class DoubleLinearOpsTensorAlgebra : TODO("ANDREI") } - override fun DoubleTensor.svd(): Triple { TODO("ALYA") } @@ -147,6 +153,72 @@ 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 + return (0 until m).asSequence().map { lu[it, it] }.fold(sign) { left, right -> left * right } + } + + public fun DoubleTensor.detLU(): DoubleTensor { + val (luTensor, pivotsTensor) = 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) = lu() + val n = shape.size + val invTensor = luTensor.zeroesLike() + + val seq = luTensor.matrixSequence().zip(pivotsTensor.vectorSequence()).zip(invTensor.matrixSequence()) + for ((luP, invMatrix) in seq) { + 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 feb8d11fa..3ac09c8c7 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 @@ -318,7 +321,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 { @@ -350,7 +377,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) @@ -359,10 +386,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) } + } +}