kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/LinearOpsTensorAlgebra.kt

@@ -14,7 +14,7 @@ public interface LinearOpsTensorAlgebra<T, TensorType : TensorStructure<T>, Inde
     public fun TensorType.qr(): TensorType
 
     //https://pytorch.org/docs/stable/generated/torch.lu.html
-    public fun TensorType.lu(tol: T): Pair<TensorType, IndexTensorType>
+    public fun TensorType.lu(): Pair<TensorType, IndexTensorType>
 
     //https://pytorch.org/docs/stable/generated/torch.lu_unpack.html
     public fun luPivot(luTensor: TensorType, pivotsTensor: IndexTensorType): Triple<TensorType, TensorType, TensorType>

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<Double, DoubleTensor, IntTensor>,
     DoubleTensorAlgebra() {
 
-    override fun DoubleTensor.inv(): DoubleTensor {
+    override fun DoubleTensor.inv(): DoubleTensor = invLU()
-        TODO("ANDREI")
-    }
 
-    override fun DoubleTensor.lu(tol: Double): Pair<DoubleTensor, IntTensor> {
+    override fun DoubleTensor.det(): DoubleTensor = detLU()
+
+    override fun DoubleTensor.lu(): Pair<DoubleTensor, IntTensor> {
 
         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<DoubleTensor, DoubleTensor, DoubleTensor> {
         TODO("ALYA")
     }
@@ -147,6 +153,72 @@ public class DoubleLinearOpsTensorAlgebra :
         TODO("ANDREI")
     }
 
+    private fun luMatrixDet(lu: Structure2D<Double>, pivots: Structure1D<Int>): 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<Double>,
+        pivots: Structure1D<Int>,
+        invMatrix : MutableStructure2D<Double>
+    ): 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 <R> DoubleLinearOpsTensorAlgebra(block: DoubleLinearOpsTensorAlgebra.() -> R): R =

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<Double, Dou
     }
 
     override fun DoubleTensor.det(): DoubleTensor {
-        TODO("ANDREI")
+        TODO()
+        /*
+        checkSquareMatrix(shape)
+
+        val n = shape.size
+        val m = shape.last()
+
+        val detTensorShape = IntArray(n - 1) { i -> 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<Double, Dou
     }
 
     public fun DoubleTensor.contentEquals(other: DoubleTensor, eqFunction: (Double, Double) -> 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<Double, Dou
     public fun DoubleTensor.eq(other: DoubleTensor, eqFunction: (Double, Double) -> 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
             }

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())}
         }
     }
 }

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) }
+    }
+}