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 new file mode 100644 index 000000000..9ffe2db61 --- /dev/null +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/BufferedTensor.kt @@ -0,0 +1,33 @@ +package space.kscience.kmath.tensors + +import space.kscience.kmath.linear.BufferMatrix +import space.kscience.kmath.nd.MutableNDBuffer +import space.kscience.kmath.structures.* +import space.kscience.kmath.structures.BufferAccessor2D + +public open class BufferedTensor( + override val shape: IntArray, + buffer: MutableBuffer +) : + TensorStructure, + MutableNDBuffer( + TensorStrides(shape), + buffer + ) { + + public operator fun get(i: Int, j: Int): T{ + check(this.dimension == 2) {"Not matrix"} + return this[intArrayOf(i, j)] + } + + public operator fun set(i: Int, j: Int, value: T): Unit{ + check(this.dimension == 2) {"Not matrix"} + this[intArrayOf(i, j)] = value + } + +} + +public class IntTensor( + shape: IntArray, + buffer: IntArray +) : BufferedTensor(shape, IntBuffer(buffer)) diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/RealTensorAlgebra.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/RealTensorAlgebra.kt index 81a290c38..f91f44519 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/RealTensorAlgebra.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/RealTensorAlgebra.kt @@ -1,25 +1,20 @@ package space.kscience.kmath.tensors -import space.kscience.kmath.nd.MutableNDBuffer import space.kscience.kmath.structures.RealBuffer import space.kscience.kmath.structures.array +import kotlin.math.abs import kotlin.math.max public class RealTensor( - override val shape: IntArray, + shape: IntArray, buffer: DoubleArray -) : - TensorStructure, - MutableNDBuffer( - TensorStrides(shape), - RealBuffer(buffer) - ) +) : BufferedTensor(shape, RealBuffer(buffer)) public class RealTensorAlgebra : TensorPartialDivisionAlgebra { override fun RealTensor.value(): Double { - check(this.shape contentEquals intArrayOf(1)) { + check(this.shape contentEquals intArrayOf(1)) { "Inconsistent value for tensor of shape ${shape.toList()}" } return this.buffer.array[0] @@ -40,7 +35,10 @@ public class RealTensorAlgebra : TensorPartialDivisionAlgebra { - TODO() + override fun RealTensor.lu(): Pair { + // todo checks + val lu = this.copy() + val m = this.shape[0] + val pivot = IntArray(m) + + + // Initialize permutation array and parity + for (row in 0 until m) pivot[row] = row + var even = true + + for (i in 0 until m) { + var maxA = -1.0 + var iMax = i + + for (k in i until m) { + val absA = abs(lu[k, i]) + if (absA > maxA) { + maxA = absA + iMax = k + } + } + + //todo check singularity + + if (iMax != i) { + + val j = pivot[i] + pivot[i] = pivot[iMax] + pivot[iMax] = j + even != even + + for (k in 0 until m) { + val tmp = lu[i, k] + lu[i, k] = lu[iMax, k] + lu[iMax, k] = tmp + } + + } + + for (j in i + 1 until m) { + lu[j, i] /= lu[i, i] + for (k in i + 1 until m) { + lu[j, k] -= lu[j, i] * lu[i, k] + } + } + } + return Pair(lu, IntTensor(intArrayOf(m), pivot)) } - override fun luUnpack(A_LU: RealTensor, pivots: RealTensor): Triple { - TODO("Not yet implemented") + override fun luUnpack(lu: RealTensor, pivots: IntTensor): Triple { + // todo checks + val n = lu.shape[0] + val p = zeroesLike(lu) + pivots.buffer.array.forEachIndexed { i, pivot -> + p[i, pivot] = 1.0 + } + val l = zeroesLike(lu) + val u = zeroesLike(lu) + + 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) } override fun RealTensor.svd(): Triple { 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 ffc0d28cc..ec257e45c 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 @@ -74,11 +74,12 @@ public interface TensorPartialDivisionAlgebra //https://pytorch.org/docs/stable/generated/torch.log.html public fun TensorType.log(): TensorType + // todo change type of pivots //https://pytorch.org/docs/stable/generated/torch.lu.html - public fun TensorType.lu(): Pair + public fun TensorType.lu(): Pair //https://pytorch.org/docs/stable/generated/torch.lu_unpack.html - public fun luUnpack(A_LU: TensorType, pivots: TensorType): Triple + public fun luUnpack(A_LU: TensorType, pivots: IntTensor): Triple //https://pytorch.org/docs/stable/generated/torch.svd.html public fun TensorType.svd(): Triple