From 1fa0da2810d69aa6d9dac669b4ebc353dd5265d1 Mon Sep 17 00:00:00 2001 From: Andrei Kislitsyn Date: Wed, 17 Mar 2021 17:53:14 +0300 Subject: [PATCH] complete lu and matrix mapping --- .../kscience/kmath/tensors/BufferedTensor.kt | 95 +++++++++---------- .../tensors/DoubleLinearOpsTensorAlgebra.kt | 83 ++++++++-------- .../kscience/kmath/tensors/TensorAlgebra.kt | 2 +- 3 files changed, 85 insertions(+), 95 deletions(-) diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/BufferedTensor.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/BufferedTensor.kt index c9a401ad5..c48e47f4c 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/BufferedTensor.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/BufferedTensor.kt @@ -1,5 +1,7 @@ package space.kscience.kmath.tensors +import space.kscience.kmath.linear.Matrix +import space.kscience.kmath.nd.* import space.kscience.kmath.structures.* @@ -29,59 +31,50 @@ public open class BufferedTensor( override fun hashCode(): Int = 0 + // todo rename to vector + public inline fun forEachVector(vectorAction : (MutableStructure1D) -> Unit): Unit { + check(shape.size >= 1) {"todo"} + val vectorOffset = strides.strides[0] + val vectorShape = intArrayOf(shape.last()) + for (offset in 0 until numel step vectorOffset) { + val vector = BufferedTensor(vectorShape, buffer, offset).as1D() + vectorAction(vector) + } + } + + public inline fun forEachMatrix(matrixAction : (MutableStructure2D) -> Unit): Unit { + check(shape.size >= 2) {"todo"} + val matrixOffset = strides.strides[1] + val matrixShape = intArrayOf(shape[shape.size - 2], shape.last()) //todo better way? + for (offset in 0 until numel step matrixOffset) { + val matrix = BufferedTensor(matrixShape, buffer, offset).as2D() + matrixAction(matrix) + } + } + // todo remove code copy-pasting + + public fun vectorSequence(): Sequence> = sequence { + check(shape.size >= 1) {"todo"} + val vectorOffset = strides.strides[0] + val vectorShape = intArrayOf(shape.last()) + for (offset in 0 until numel step vectorOffset) { + val vector = BufferedTensor(vectorShape, buffer, offset).as1D() + yield(vector) + } + } + + public fun matrixSequence(): Sequence> = sequence { + check(shape.size >= 2) {"todo"} + val matrixOffset = strides.strides[1] + val matrixShape = intArrayOf(shape[shape.size - 2], shape.last()) //todo better way? + for (offset in 0 until numel step matrixOffset) { + val matrix = BufferedTensor(matrixShape, buffer, offset).as2D() + yield(matrix) + } + } + } -/* -//todo make generator mb nextMatrixIndex? -public class InnerMatrix(private val tensor: BufferedTensor){ - private var offset: Int = 0 - private val n : Int = tensor.shape.size - //stride? - private val step = tensor.shape[n - 1] * tensor.shape[n - 2] - - public operator fun get(i: Int, j: Int): T { - val index = tensor.strides.index(offset) - index[n - 2] = i - index[n - 1] = j - return tensor[index] - } - - public operator fun set(i: Int, j: Int, value: T): Unit { - val index = tensor.strides.index(offset) - index[n - 2] = i - index[n - 1] = j - tensor[index] = value - } - - public fun makeStep(){ - offset += step - } -} - -public class InnerVector(private val tensor: BufferedTensor){ - private var offset: Int = 0 - private val n : Int = tensor.shape.size - //stride? - private val step = tensor.shape[n - 1] - - public operator fun get(i: Int): T { - val index = tensor.strides.index(offset) - index[n - 1] = i - return tensor[index] - } - - public operator fun set(i: Int, value: T): Unit { - val index = tensor.strides.index(offset) - index[n - 1] = i - tensor[index] = value - } - - public fun makeStep(){ - offset += step - } -} -//todo default buffer = arrayOf(0)??? - */ public class IntTensor( shape: IntArray, diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/DoubleLinearOpsTensorAlgebra.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/DoubleLinearOpsTensorAlgebra.kt index 689eda9e0..3f44305b1 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/DoubleLinearOpsTensorAlgebra.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/DoubleLinearOpsTensorAlgebra.kt @@ -9,27 +9,21 @@ public class DoubleLinearOpsTensorAlgebra : } override fun DoubleTensor.lu(): Pair { - /* + // todo checks + val luTensor = this.copy() - val lu = InnerMatrix(luTensor) - //stride TODO!!! move to generator? - var matCnt = 1 - for (i in 0 until this.shape.size - 2) { - matCnt *= this.shape[i] - } + val n = this.shape.size - val m = this.shape[n - 1] - val pivotsShape = IntArray(n - 1) { i -> - this.shape[i] - } + val m = this.shape.last() + val pivotsShape = IntArray(n - 1) { i -> this.shape[i] } val pivotsTensor = IntTensor( pivotsShape, - IntArray(matCnt * m) { 0 } + IntArray(pivotsShape.reduce(Int::times)) { 0 } //todo default??? ) - val pivot = InnerVector(pivotsTensor) - for (i in 0 until matCnt) { - for (row in 0 until m) pivot[row] = row + + for ((lu, pivots) in luTensor.matrixSequence().zip(pivotsTensor.vectorSequence())){ + for (row in 0 until m) pivots[row] = row for (i in 0 until m) { var maxA = -1.0 @@ -47,9 +41,9 @@ public class DoubleLinearOpsTensorAlgebra : if (iMax != i) { - val j = pivot[i] - pivot[i] = pivot[iMax] - pivot[iMax] = j + val j = pivots[i] + pivots[i] = pivots[iMax] + pivots[iMax] = j for (k in 0 until m) { val tmp = lu[i, k] @@ -66,42 +60,45 @@ public class DoubleLinearOpsTensorAlgebra : } } } - lu.makeStep() - pivot.makeStep() } - return Pair(luTensor, pivotsTensor)*/ - TODO("Andrei, use view, get, as2D, as1D") + return Pair(luTensor, pivotsTensor) + } - override fun luPivot(lu: DoubleTensor, pivots: IntTensor): Triple { - /* - // todo checks - val n = lu.shape[0] - val p = lu.zeroesLike() - pivots.buffer.unsafeToIntArray().forEachIndexed { i, pivot -> - p[i, pivot] = 1.0 + override fun luPivot(luTensor: DoubleTensor, pivotsTensor: IntTensor): Triple { + //todo checks + val n = luTensor.shape.last() + val pTensor = luTensor.zeroesLike() + for ((p, pivot) in pTensor.matrixSequence().zip(pivotsTensor.vectorSequence())){ + for (i in 0 until n){ + p[i, pivot[i]] = 1.0 + } } - val l = lu.zeroesLike() - val u = lu.zeroesLike() - for (i in 0 until n) { - for (j in 0 until n) { - if (i == j) { - l[i, j] = 1.0 - } - if (j < i) { - l[i, j] = lu[i, j] - } - if (j >= i) { - u[i, j] = lu[i, j] + val lTensor = luTensor.zeroesLike() + val uTensor = luTensor.zeroesLike() + + for ((pairLU, lu) in lTensor.matrixSequence().zip(uTensor.matrixSequence()).zip(luTensor.matrixSequence())){ + val (l, u) = pairLU + for (i in 0 until n) { + for (j in 0 until n) { + if (i == j) { + l[i, j] = 1.0 + } + if (j < i) { + l[i, j] = lu[i, j] + } + if (j >= i) { + u[i, j] = lu[i, j] + } } } } - return Triple(p, l, u)*/ - TODO("Andrei, first we need implement get(Int)") + return Triple(pTensor, lTensor, uTensor) + } override fun DoubleTensor.cholesky(): DoubleTensor { diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/TensorAlgebra.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/TensorAlgebra.kt index b1adf2962..7ec920d88 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/TensorAlgebra.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/TensorAlgebra.kt @@ -5,7 +5,7 @@ public interface TensorAlgebra> { public fun zeros(shape: IntArray): TensorType - public fun TensorType.zeroesLike(): TensorType + public fun TensorType.zeroesLike(): TensorType // mb it shouldn't be tensor but algebra method (like in numpy/torch) ? public fun ones(shape: IntArray): TensorType public fun TensorType.onesLike(): TensorType