kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUDecomposition.kt

@@ -188,57 +188,57 @@ class LUPDecompositionBuilder<T : Comparable<T>, F : Field<T>>(val context: F, v
 }
 
 
-//class LUSolver<T : Comparable<T>, F : Field<T>>(val singularityCheck: (T) -> Boolean) : LinearSolver<T, F> {
+class LUSolver<T : Comparable<T>, F : Field<T>>(val singularityCheck: (T) -> Boolean) : LinearSolver<T, F> {
-//
+
-//
+
-//    override fun solve(a: Matrix<T>, b: Matrix<T>): Matrix<T> {
+    override fun solve(a: Matrix<T>, b: Matrix<T>): Matrix<T> {
-//        val decomposition = LUPDecompositionBuilder(ring, singularityCheck).decompose(a)
+        val decomposition = LUPDecompositionBuilder(ring, singularityCheck).decompose(a)
-//
+
-//        if (b.rowNum != a.colNum) {
+        if (b.rowNum != a.colNum) {
-//            error("Matrix dimension mismatch expected ${a.rowNum}, but got ${b.colNum}")
+            error("Matrix dimension mismatch expected ${a.rowNum}, but got ${b.colNum}")
-//        }
+        }
-//
+
-//
+
-////        val bp = Array(a.rowNum) { Array<T>(b.colNum){ring.zero} }
+//        val bp = Array(a.rowNum) { Array<T>(b.colNum){ring.zero} }
-////        for (row in 0 until a.rowNum) {
+//        for (row in 0 until a.rowNum) {
-////            val bpRow = bp[row]
+//            val bpRow = bp[row]
-////            val pRow = decomposition.pivot[row]
+//            val pRow = decomposition.pivot[row]
-////            for (col in 0 until b.colNum) {
+//            for (col in 0 until b.colNum) {
-////                bpRow[col] = b[pRow, col]
+//                bpRow[col] = b[pRow, col]
-////            }
-////        }
-//
-//        // Apply permutations to b
-//        val bp = produce(a.rowNum, a.colNum) { i, j -> b[decomposition.pivot[i], j] }
-//
-//        // Solve LY = b
-//        for (col in 0 until a.rowNum) {
-//            val bpCol = bp[col]
-//            for (i in col + 1 until a.rowNum) {
-//                val bpI = bp[i]
-//                val luICol = decomposition.lu[i, col]
-//                for (j in 0 until b.colNum) {
-//                    bpI[j] -= bpCol[j] * luICol
-//                }
-//            }
-//        }
-//
-//        // Solve UX = Y
-//        for (col in a.rowNum - 1 downTo 0) {
-//            val bpCol = bp[col]
-//            val luDiag = decomposition.lu[col, col]
-//            for (j in 0 until b.colNum) {
-//                bpCol[j] /= luDiag
-//            }
-//            for (i in 0 until col) {
-//                val bpI = bp[i]
-//                val luICol = decomposition.lu[i, col]
-//                for (j in 0 until b.colNum) {
-//                    bpI[j] -= bpCol[j] * luICol
-//                }
-//            }
-//        }
-//
-//        return produce(a.rowNum, a.colNum) { i, j -> bp[i][j] }
 //            }
 //        }
+
+        // Apply permutations to b
+        val bp = produce(a.rowNum, a.colNum) { i, j -> b[decomposition.pivot[i], j] }
+
+        // Solve LY = b
+        for (col in 0 until a.rowNum) {
+            val bpCol = bp[col]
+            for (i in col + 1 until a.rowNum) {
+                val bpI = bp[i]
+                val luICol = decomposition.lu[i, col]
+                for (j in 0 until b.colNum) {
+                    bpI[j] -= bpCol[j] * luICol
+                }
+            }
+        }
+
+        // Solve UX = Y
+        for (col in a.rowNum - 1 downTo 0) {
+            val bpCol = bp[col]
+            val luDiag = decomposition.lu[col, col]
+            for (j in 0 until b.colNum) {
+                bpCol[j] /= luDiag
+            }
+            for (i in 0 until col) {
+                val bpI = bp[i]
+                val luICol = decomposition.lu[i, col]
+                for (j in 0 until b.colNum) {
+                    bpI[j] -= bpCol[j] * luICol
+                }
+            }
+        }
+
+        return produce(a.rowNum, a.colNum) { i, j -> bp[i][j] }
+    }
+}

kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt

@@ -10,15 +10,32 @@ typealias MutableBufferFactory<T> = (Int, (Int) -> T) -> MutableBuffer<T>
  */
 interface Buffer<T> {
 
+    /**
+     * The size of the buffer
+     */
     val size: Int
 
+    /**
+     * Get element at given index
+     */
     operator fun get(index: Int): T
 
+    /**
+     * Iterate over all elements
+     */
     operator fun iterator(): Iterator<T>
 
+    /**
+     * Check content eqiality with another buffer
+     */
     fun contentEquals(other: Buffer<*>): Boolean =
         asSequence().mapIndexed { index, value -> value == other[index] }.all { it }
 
+    /**
+     * Map the contents of the buffer to new buffer of the same type
+     */
+    fun transform(transformation: (index: Int, value: T) -> T): Buffer<T>
+
     companion object {
 
         /**
@@ -61,6 +78,8 @@ interface MutableBuffer<T> : Buffer<T> {
      */
     fun copy(): MutableBuffer<T>
 
+    fun transformInPlace(transformation: (index: Int, value: T) -> T)
+
     companion object {
         /**
          * Create a boxing mutable buffer of given type
@@ -93,6 +112,9 @@ inline class ListBuffer<T>(private val list: List<T>) : Buffer<T> {
     override fun get(index: Int): T = list[index]
 
     override fun iterator(): Iterator<T> = list.iterator()
+
+    override fun transform(transformation: (index: Int, value: T) -> T): ListBuffer<T> =
+        list.mapIndexed(transformation).asBuffer()
 }
 
 fun <T> List<T>.asBuffer() = ListBuffer(this)
@@ -110,11 +132,20 @@ inline class MutableListBuffer<T>(private val list: MutableList<T>) : MutableBuf
 
     override fun iterator(): Iterator<T> = list.iterator()
 
+    override fun transform(transformation: (index: Int, value: T) -> T): ListBuffer<T> =
+        list.mapIndexed(transformation).asBuffer()
+
     override fun copy(): MutableBuffer<T> = MutableListBuffer(ArrayList(list))
+
+    override fun transformInPlace(transformation: (index: Int, value: T) -> T) = list.forEachIndexed { index, value ->
+        list[index] = transformation(index, value)
+    }
 }
 
+fun <T> MutableList<T>.asBuffer() = MutableListBuffer(this)
+
 class ArrayBuffer<T>(private val array: Array<T>) : MutableBuffer<T> {
-    //Can't inline because array invariant
+    //Can't inline because array is invariant
     override val size: Int
         get() = array.size
 
@@ -127,7 +158,17 @@ class ArrayBuffer<T>(private val array: Array<T>) : MutableBuffer<T> {
     override fun iterator(): Iterator<T> = array.iterator()
 
     override fun copy(): MutableBuffer<T> = ArrayBuffer(array.copyOf())
+
+    override fun transform(transformation: (index: Int, value: T) -> T): ArrayBuffer<T> =
+        Array(size) { index -> transformation(index, get(index)) }.asBuffer()
+
+
+    override fun transformInPlace(transformation: (index: Int, value: T) -> T) = array.forEachIndexed { index, value ->
+        array[index] = transformation(index, value)
     }
+}
+
+fun <T> Array<T>.asBuffer() = ArrayBuffer(this)
 
 inline class DoubleBuffer(private val array: DoubleArray) : MutableBuffer<Double> {
     override val size: Int get() = array.size
@@ -141,7 +182,17 @@ inline class DoubleBuffer(private val array: DoubleArray) : MutableBuffer<Double
     override fun iterator(): Iterator<Double> = array.iterator()
 
     override fun copy(): MutableBuffer<Double> = DoubleBuffer(array.copyOf())
+
+    override fun transform(transformation: (index: Int, value: Double) -> Double): DoubleBuffer =
+        DoubleArray(size) { index -> transformation(index, get(index)) }.asBuffer()
+
+    override fun transformInPlace(transformation: (index: Int, value: Double) -> Double) =
+        array.forEachIndexed { index, value ->
+            array[index] = transformation(index, value)
         }
+}
+
+fun DoubleArray.asBuffer() = DoubleBuffer(this)
 
 inline class ShortBuffer(private val array: ShortArray) : MutableBuffer<Short> {
     override val size: Int get() = array.size
@@ -155,7 +206,17 @@ inline class ShortBuffer(private val array: ShortArray) : MutableBuffer<Short> {
     override fun iterator(): Iterator<Short> = array.iterator()
 
     override fun copy(): MutableBuffer<Short> = ShortBuffer(array.copyOf())
+
+    override fun transform(transformation: (index: Int, value: Short) -> Short): ShortBuffer =
+        ShortArray(size) { index -> transformation(index, get(index)) }.asBuffer()
+
+    override fun transformInPlace(transformation: (index: Int, value: Short) -> Short) =
+        array.forEachIndexed { index, value ->
+            array[index] = transformation(index, value)
         }
+}
+
+fun ShortArray.asBuffer() = ShortBuffer(this)
 
 inline class IntBuffer(private val array: IntArray) : MutableBuffer<Int> {
     override val size: Int get() = array.size
@@ -169,7 +230,17 @@ inline class IntBuffer(private val array: IntArray) : MutableBuffer<Int> {
     override fun iterator(): Iterator<Int> = array.iterator()
 
     override fun copy(): MutableBuffer<Int> = IntBuffer(array.copyOf())
+
+    override fun transform(transformation: (index: Int, value: Int) -> Int): IntBuffer =
+        IntArray(size) { index -> transformation(index, get(index)) }.asBuffer()
+
+    override fun transformInPlace(transformation: (index: Int, value: Int) -> Int) =
+        array.forEachIndexed { index, value ->
+            array[index] = transformation(index, value)
         }
+}
+
+fun IntArray.asBuffer() = IntBuffer(this)
 
 inline class LongBuffer(private val array: LongArray) : MutableBuffer<Long> {
     override val size: Int get() = array.size
@@ -183,7 +254,17 @@ inline class LongBuffer(private val array: LongArray) : MutableBuffer<Long> {
     override fun iterator(): Iterator<Long> = array.iterator()
 
     override fun copy(): MutableBuffer<Long> = LongBuffer(array.copyOf())
+
+    override fun transform(transformation: (index: Int, value: Long) -> Long): LongBuffer =
+        LongArray(size) { index -> transformation(index, get(index)) }.asBuffer()
+
+    override fun transformInPlace(transformation: (index: Int, value: Long) -> Long) =
+        array.forEachIndexed { index, value ->
+            array[index] = transformation(index, value)
         }
+}
+
+fun LongArray.asBuffer() = LongBuffer(this)
 
 inline class ReadOnlyBuffer<T>(private val buffer: MutableBuffer<T>) : Buffer<T> {
     override val size: Int get() = buffer.size
@@ -191,6 +272,8 @@ inline class ReadOnlyBuffer<T>(private val buffer: MutableBuffer<T>) : Buffer<T>
     override fun get(index: Int): T = buffer.get(index)
 
     override fun iterator(): Iterator<T> = buffer.iterator()
+
+    override fun transform(transformation: (index: Int, value: T) -> T): Buffer<T> = buffer.transform(transformation)
 }
 
 /**
@@ -210,6 +293,9 @@ class VirtualBuffer<T>(override val size: Int, private val generator: (Int) -> T
         }
     }
 
+    // TODO this composition could become very compex very fast, replace it by boxed generator?
+    override fun transform(transformation: (index: Int, value: T) -> T): Buffer<T> =
+        VirtualBuffer(size) { index -> transformation(index, generator(index)) }
 }
 
 /**