diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUDecomposition.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUDecomposition.kt index 7b645c74b..62e9ba903 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUDecomposition.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUDecomposition.kt @@ -188,57 +188,57 @@ class LUPDecompositionBuilder, F : Field>(val context: F, v } -//class LUSolver, F : Field>(val singularityCheck: (T) -> Boolean) : LinearSolver { -// -// -// override fun solve(a: Matrix, b: Matrix): Matrix { -// val decomposition = LUPDecompositionBuilder(ring, singularityCheck).decompose(a) -// -// if (b.rowNum != a.colNum) { -// error("Matrix dimension mismatch expected ${a.rowNum}, but got ${b.colNum}") -// } -// -// -//// val bp = Array(a.rowNum) { Array(b.colNum){ring.zero} } -//// for (row in 0 until a.rowNum) { -//// val bpRow = bp[row] -//// val pRow = decomposition.pivot[row] -//// for (col in 0 until b.colNum) { -//// 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 -// } +class LUSolver, F : Field>(val singularityCheck: (T) -> Boolean) : LinearSolver { + + + override fun solve(a: Matrix, b: Matrix): Matrix { + val decomposition = LUPDecompositionBuilder(ring, singularityCheck).decompose(a) + + if (b.rowNum != a.colNum) { + error("Matrix dimension mismatch expected ${a.rowNum}, but got ${b.colNum}") + } + + +// val bp = Array(a.rowNum) { Array(b.colNum){ring.zero} } +// for (row in 0 until a.rowNum) { +// val bpRow = bp[row] +// val pRow = decomposition.pivot[row] +// for (col in 0 until b.colNum) { +// bpRow[col] = b[pRow, col] // } // } -// -// // 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] } + } +} diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt index a8a9d7ed3..14fcfc2ff 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt @@ -10,15 +10,32 @@ typealias MutableBufferFactory = (Int, (Int) -> T) -> MutableBuffer */ interface Buffer { + /** + * 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 + /** + * 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 + companion object { /** @@ -61,6 +78,8 @@ interface MutableBuffer : Buffer { */ fun copy(): MutableBuffer + 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(private val list: List) : Buffer { override fun get(index: Int): T = list[index] override fun iterator(): Iterator = list.iterator() + + override fun transform(transformation: (index: Int, value: T) -> T): ListBuffer = + list.mapIndexed(transformation).asBuffer() } fun List.asBuffer() = ListBuffer(this) @@ -110,11 +132,20 @@ inline class MutableListBuffer(private val list: MutableList) : MutableBuf override fun iterator(): Iterator = list.iterator() + override fun transform(transformation: (index: Int, value: T) -> T): ListBuffer = + list.mapIndexed(transformation).asBuffer() + override fun copy(): MutableBuffer = MutableListBuffer(ArrayList(list)) + + override fun transformInPlace(transformation: (index: Int, value: T) -> T) = list.forEachIndexed { index, value -> + list[index] = transformation(index, value) + } } +fun MutableList.asBuffer() = MutableListBuffer(this) + class ArrayBuffer(private val array: Array) : MutableBuffer { - //Can't inline because array invariant + //Can't inline because array is invariant override val size: Int get() = array.size @@ -127,8 +158,18 @@ class ArrayBuffer(private val array: Array) : MutableBuffer { override fun iterator(): Iterator = array.iterator() override fun copy(): MutableBuffer = ArrayBuffer(array.copyOf()) + + override fun transform(transformation: (index: Int, value: T) -> T): ArrayBuffer = + 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 Array.asBuffer() = ArrayBuffer(this) + inline class DoubleBuffer(private val array: DoubleArray) : MutableBuffer { override val size: Int get() = array.size @@ -141,8 +182,18 @@ inline class DoubleBuffer(private val array: DoubleArray) : MutableBuffer = array.iterator() override fun copy(): MutableBuffer = 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 { override val size: Int get() = array.size @@ -155,8 +206,18 @@ inline class ShortBuffer(private val array: ShortArray) : MutableBuffer { override fun iterator(): Iterator = array.iterator() override fun copy(): MutableBuffer = 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 { override val size: Int get() = array.size @@ -169,8 +230,18 @@ inline class IntBuffer(private val array: IntArray) : MutableBuffer { override fun iterator(): Iterator = array.iterator() override fun copy(): MutableBuffer = 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 { override val size: Int get() = array.size @@ -183,14 +254,26 @@ inline class LongBuffer(private val array: LongArray) : MutableBuffer { override fun iterator(): Iterator = array.iterator() override fun copy(): MutableBuffer = 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(private val buffer: MutableBuffer) : Buffer { override val size: Int get() = buffer.size override fun get(index: Int): T = buffer.get(index) override fun iterator(): Iterator = buffer.iterator() + + override fun transform(transformation: (index: Int, value: T) -> T): Buffer = buffer.transform(transformation) } /** @@ -210,6 +293,9 @@ class VirtualBuffer(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 = + VirtualBuffer(size) { index -> transformation(index, generator(index)) } } /**