diff --git a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/broadcastUtils.kt b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/broadcastUtils.kt index 2d4d6f0f7..58d8654af 100644 --- a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/broadcastUtils.kt +++ b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/broadcastUtils.kt @@ -2,7 +2,7 @@ package space.kscience.kmath.tensors.core import kotlin.math.max -internal inline fun multiIndexBroadCasting(tensor: DoubleTensor, resTensor: DoubleTensor, linearSize: Int) { +internal fun multiIndexBroadCasting(tensor: DoubleTensor, resTensor: DoubleTensor, linearSize: Int) { for (linearIndex in 0 until linearSize) { val totalMultiIndex = resTensor.linearStructure.index(linearIndex) val curMultiIndex = tensor.shape.copyOf() @@ -23,7 +23,7 @@ internal inline fun multiIndexBroadCasting(tensor: DoubleTensor, resTensor: Doub } } -internal inline fun broadcastShapes(vararg shapes: IntArray): IntArray { +internal fun broadcastShapes(vararg shapes: IntArray): IntArray { var totalDim = 0 for (shape in shapes) { totalDim = max(totalDim, shape.size) @@ -51,7 +51,7 @@ internal inline fun broadcastShapes(vararg shapes: IntArray): IntArray { return totalShape } -internal inline fun broadcastTo(tensor: DoubleTensor, newShape: IntArray): DoubleTensor { +internal fun broadcastTo(tensor: DoubleTensor, newShape: IntArray): DoubleTensor { if (tensor.shape.size > newShape.size) { throw RuntimeException("Tensor is not compatible with the new shape") } @@ -71,7 +71,7 @@ internal inline fun broadcastTo(tensor: DoubleTensor, newShape: IntArray): Doubl return resTensor } -internal inline fun broadcastTensors(vararg tensors: DoubleTensor): List { +internal fun broadcastTensors(vararg tensors: DoubleTensor): List { val totalShape = broadcastShapes(*(tensors.map { it.shape }).toTypedArray()) val n = totalShape.reduce { acc, i -> acc * i } @@ -85,7 +85,7 @@ internal inline fun broadcastTensors(vararg tensors: DoubleTensor): List { +internal fun broadcastOuterTensors(vararg tensors: DoubleTensor): List { val onlyTwoDims = tensors.asSequence().onEach { require(it.shape.size >= 2) { throw RuntimeException("Tensors must have at least 2 dimensions") @@ -99,46 +99,45 @@ internal inline fun broadcastOuterTensors(vararg tensors: DoubleTensor): List acc * i } - val res = ArrayList(0) - for (tensor in tensors) { - val matrixShape = tensor.shape.sliceArray(tensor.shape.size - 2 until tensor.shape.size).copyOf() - val matrixSize = matrixShape[0] * matrixShape[1] - val matrix = DoubleTensor(matrixShape, DoubleArray(matrixSize)) + return buildList { + for (tensor in tensors) { + val matrixShape = tensor.shape.sliceArray(tensor.shape.size - 2 until tensor.shape.size).copyOf() + val matrixSize = matrixShape[0] * matrixShape[1] + val matrix = DoubleTensor(matrixShape, DoubleArray(matrixSize)) - val outerTensor = DoubleTensor(totalShape, DoubleArray(n)) - val resTensor = DoubleTensor(totalShape + matrixShape, DoubleArray(n * matrixSize)) + val outerTensor = DoubleTensor(totalShape, DoubleArray(n)) + val resTensor = DoubleTensor(totalShape + matrixShape, DoubleArray(n * matrixSize)) - for (linearIndex in 0 until n) { - val totalMultiIndex = outerTensor.linearStructure.index(linearIndex) - var curMultiIndex = tensor.shape.sliceArray(0..tensor.shape.size - 3).copyOf() - curMultiIndex = IntArray(totalMultiIndex.size - curMultiIndex.size) { 1 } + curMultiIndex + for (linearIndex in 0 until n) { + val totalMultiIndex = outerTensor.linearStructure.index(linearIndex) + var curMultiIndex = tensor.shape.sliceArray(0..tensor.shape.size - 3).copyOf() + curMultiIndex = IntArray(totalMultiIndex.size - curMultiIndex.size) { 1 } + curMultiIndex - val newTensor = DoubleTensor(curMultiIndex + matrixShape, tensor.mutableBuffer.array()) + val newTensor = DoubleTensor(curMultiIndex + matrixShape, tensor.mutableBuffer.array()) - for (i in curMultiIndex.indices) { - if (curMultiIndex[i] != 1) { - curMultiIndex[i] = totalMultiIndex[i] - } else { - curMultiIndex[i] = 0 + for (i in curMultiIndex.indices) { + if (curMultiIndex[i] != 1) { + curMultiIndex[i] = totalMultiIndex[i] + } else { + curMultiIndex[i] = 0 + } + } + + for (i in 0 until matrixSize) { + val curLinearIndex = newTensor.linearStructure.offset( + curMultiIndex + + matrix.linearStructure.index(i) + ) + val newLinearIndex = resTensor.linearStructure.offset( + totalMultiIndex + + matrix.linearStructure.index(i) + ) + + resTensor.mutableBuffer.array()[resTensor.bufferStart + newLinearIndex] = + newTensor.mutableBuffer.array()[newTensor.bufferStart + curLinearIndex] } } - - for (i in 0 until matrixSize) { - val curLinearIndex = newTensor.linearStructure.offset( - curMultiIndex + - matrix.linearStructure.index(i) - ) - val newLinearIndex = resTensor.linearStructure.offset( - totalMultiIndex + - matrix.linearStructure.index(i) - ) - - resTensor.mutableBuffer.array()[resTensor.bufferStart + newLinearIndex] = - newTensor.mutableBuffer.array()[newTensor.bufferStart + curLinearIndex] - } + add(resTensor) } - res += resTensor } - - return res } \ No newline at end of file diff --git a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/checks.kt b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/checks.kt index fd5c8413d..fd98be8b2 100644 --- a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/checks.kt +++ b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/checks.kt @@ -5,38 +5,38 @@ import space.kscience.kmath.tensors.core.algebras.DoubleLinearOpsTensorAlgebra import space.kscience.kmath.tensors.core.algebras.DoubleTensorAlgebra -internal inline fun checkEmptyShape(shape: IntArray): Unit = +internal fun checkEmptyShape(shape: IntArray): Unit = check(shape.isNotEmpty()) { "Illegal empty shape provided" } -internal inline fun checkEmptyDoubleBuffer(buffer: DoubleArray): Unit = +internal fun checkEmptyDoubleBuffer(buffer: DoubleArray): Unit = check(buffer.isNotEmpty()) { "Illegal empty buffer provided" } -internal inline fun checkBufferShapeConsistency(shape: IntArray, buffer: DoubleArray): Unit = +internal fun checkBufferShapeConsistency(shape: IntArray, buffer: DoubleArray): Unit = check(buffer.size == shape.reduce(Int::times)) { "Inconsistent shape ${shape.toList()} for buffer of size ${buffer.size} provided" } -internal inline fun checkShapesCompatible(a: TensorStructure, b: TensorStructure): Unit = +internal fun checkShapesCompatible(a: TensorStructure, b: TensorStructure): Unit = check(a.shape contentEquals b.shape) { "Incompatible shapes ${a.shape.toList()} and ${b.shape.toList()} " } -internal inline fun checkTranspose(dim: Int, i: Int, j: Int): Unit = +internal fun checkTranspose(dim: Int, i: Int, j: Int): Unit = check((i < dim) and (j < dim)) { "Cannot transpose $i to $j for a tensor of dim $dim" } -internal inline fun checkView(a: TensorStructure, shape: IntArray): Unit = +internal fun checkView(a: TensorStructure, shape: IntArray): Unit = check(a.shape.reduce(Int::times) == shape.reduce(Int::times)) -internal inline fun checkSquareMatrix(shape: IntArray): Unit { +internal fun checkSquareMatrix(shape: IntArray): Unit { val n = shape.size check(n >= 2) { "Expected tensor with 2 or more dimensions, got size $n instead" @@ -46,14 +46,14 @@ internal inline fun checkSquareMatrix(shape: IntArray): Unit { } } -internal inline fun DoubleTensorAlgebra.checkSymmetric( +internal fun DoubleTensorAlgebra.checkSymmetric( tensor: TensorStructure, epsilon: Double = 1e-6 ): Unit = check(tensor.eq(tensor.transpose(), epsilon)) { "Tensor is not symmetric about the last 2 dimensions at precision $epsilon" } -internal inline fun DoubleLinearOpsTensorAlgebra.checkPositiveDefinite( +internal fun DoubleLinearOpsTensorAlgebra.checkPositiveDefinite( tensor: DoubleTensor, epsilon: Double = 1e-6 ): Unit { checkSymmetric(tensor, epsilon) diff --git a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/linUtils.kt b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/linUtils.kt index 776a5e4d5..a152b3a17 100644 --- a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/linUtils.kt +++ b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/linUtils.kt @@ -13,7 +13,7 @@ import kotlin.math.sign import kotlin.math.sqrt -internal inline fun BufferedTensor.vectorSequence(): Sequence> = sequence { +internal fun BufferedTensor.vectorSequence(): Sequence> = sequence { val n = shape.size val vectorOffset = shape[n - 1] val vectorShape = intArrayOf(shape.last()) @@ -23,9 +23,9 @@ internal inline fun BufferedTensor.vectorSequence(): Sequence BufferedTensor.matrixSequence(): Sequence> = sequence { - check(shape.size >= 2) { "todo" } +internal fun BufferedTensor.matrixSequence(): Sequence> = sequence { val n = shape.size + check(n >= 2) { "Expected tensor with 2 or more dimensions, got size $n" } val matrixOffset = shape[n - 1] * shape[n - 2] val matrixShape = intArrayOf(shape[n - 2], shape[n - 1]) for (offset in 0 until numElements step matrixOffset) { @@ -46,8 +46,7 @@ internal inline fun BufferedTensor.forEachMatrix(matrixAction: (BufferedT } } - -internal inline fun dotHelper( +internal fun dotHelper( a: MutableStructure2D, b: MutableStructure2D, res: MutableStructure2D, @@ -64,10 +63,11 @@ internal inline fun dotHelper( } } -internal inline fun luHelper( +internal fun luHelper( lu: MutableStructure2D, pivots: MutableStructure1D, - epsilon: Double): Boolean { + epsilon: Double +): Boolean { val m = lu.rowNum @@ -114,7 +114,7 @@ internal inline fun luHelper( return false } -internal inline fun BufferedTensor.setUpPivots(): IntTensor { +internal fun BufferedTensor.setUpPivots(): IntTensor { val n = this.shape.size val m = this.shape.last() val pivotsShape = IntArray(n - 1) { i -> this.shape[i] } @@ -126,22 +126,23 @@ internal inline fun BufferedTensor.setUpPivots(): IntTensor { ) } -internal inline fun DoubleLinearOpsTensorAlgebra.computeLU( +internal fun DoubleLinearOpsTensorAlgebra.computeLU( tensor: DoubleTensor, - epsilon: Double): Pair? { + epsilon: Double +): Pair? { checkSquareMatrix(tensor.shape) val luTensor = tensor.copy() val pivotsTensor = tensor.setUpPivots() for ((lu, pivots) in luTensor.matrixSequence().zip(pivotsTensor.vectorSequence())) - if(luHelper(lu.as2D(), pivots.as1D(), epsilon)) + if (luHelper(lu.as2D(), pivots.as1D(), epsilon)) return null return Pair(luTensor, pivotsTensor) } -internal inline fun pivInit( +internal fun pivInit( p: MutableStructure2D, pivot: MutableStructure1D, n: Int @@ -151,7 +152,7 @@ internal inline fun pivInit( } } -internal inline fun luPivotHelper( +internal fun luPivotHelper( l: MutableStructure2D, u: MutableStructure2D, lu: MutableStructure2D, @@ -172,7 +173,7 @@ internal inline fun luPivotHelper( } } -internal inline fun choleskyHelper( +internal fun choleskyHelper( a: MutableStructure2D, l: MutableStructure2D, n: Int @@ -193,7 +194,7 @@ internal inline fun choleskyHelper( } } -internal inline fun luMatrixDet(lu: MutableStructure2D, pivots: MutableStructure1D): Double { +internal fun luMatrixDet(lu: MutableStructure2D, pivots: MutableStructure1D): Double { if (lu[0, 0] == 0.0) { return 0.0 } @@ -202,7 +203,7 @@ internal inline fun luMatrixDet(lu: MutableStructure2D, pivots: MutableS return (0 until m).asSequence().map { lu[it, it] }.fold(sign) { left, right -> left * right } } -internal inline fun luMatrixInv( +internal fun luMatrixInv( lu: MutableStructure2D, pivots: MutableStructure1D, invMatrix: MutableStructure2D @@ -229,7 +230,7 @@ internal inline fun luMatrixInv( } } -internal inline fun DoubleLinearOpsTensorAlgebra.qrHelper( +internal fun DoubleLinearOpsTensorAlgebra.qrHelper( matrix: DoubleTensor, q: DoubleTensor, r: MutableStructure2D @@ -259,7 +260,7 @@ internal inline fun DoubleLinearOpsTensorAlgebra.qrHelper( } } -internal inline fun DoubleLinearOpsTensorAlgebra.svd1d(a: DoubleTensor, epsilon: Double = 1e-10): DoubleTensor { +internal fun DoubleLinearOpsTensorAlgebra.svd1d(a: DoubleTensor, epsilon: Double = 1e-10): DoubleTensor { val (n, m) = a.shape var v: DoubleTensor val b: DoubleTensor @@ -283,7 +284,7 @@ internal inline fun DoubleLinearOpsTensorAlgebra.svd1d(a: DoubleTensor, epsilon: } } -internal inline fun DoubleLinearOpsTensorAlgebra.svdHelper( +internal fun DoubleLinearOpsTensorAlgebra.svdHelper( matrix: DoubleTensor, USV: Pair, Pair, BufferedTensor>>, m: Int, n: Int, epsilon: Double @@ -335,7 +336,7 @@ internal inline fun DoubleLinearOpsTensorAlgebra.svdHelper( } } -internal inline fun cleanSymHelper(matrix: MutableStructure2D, n: Int) { +internal fun cleanSymHelper(matrix: MutableStructure2D, n: Int) { for (i in 0 until n) for (j in 0 until n) { if (i == j) { diff --git a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/utils.kt b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/utils.kt index cb23dbdd6..58d280307 100644 --- a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/utils.kt +++ b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/utils.kt @@ -23,19 +23,19 @@ internal fun Buffer.array(): DoubleArray = when (this) { else -> this.toDoubleArray() } -internal inline fun getRandomNormals(n: Int, seed: Long): DoubleArray { +internal fun getRandomNormals(n: Int, seed: Long): DoubleArray { val distribution = GaussianSampler(0.0, 1.0) val generator = RandomGenerator.default(seed) return distribution.sample(generator).nextBufferBlocking(n).toDoubleArray() } -internal inline fun getRandomUnitVector(n: Int, seed: Long): DoubleArray { +internal fun getRandomUnitVector(n: Int, seed: Long): DoubleArray { val unnorm = getRandomNormals(n, seed) val norm = sqrt(unnorm.map { it * it }.sum()) return unnorm.map { it / norm }.toDoubleArray() } -internal inline fun minusIndexFrom(n: Int, i: Int) : Int = if (i >= 0) i else { +internal fun minusIndexFrom(n: Int, i: Int): Int = if (i >= 0) i else { val ii = n + i check(ii >= 0) { "Out of bound index $i for tensor of dim $n" @@ -43,27 +43,28 @@ internal inline fun minusIndexFrom(n: Int, i: Int) : Int = if (i >= 0) i else { ii } -internal inline fun BufferedTensor.minusIndex(i: Int): Int = minusIndexFrom(this.dimension, i) +internal fun BufferedTensor.minusIndex(i: Int): Int = minusIndexFrom(this.dimension, i) -internal inline fun format(value: Double, digits: Int = 4): String { +internal fun format(value: Double, digits: Int = 4): String { val ten = 10.0 - val approxOrder = if(value == 0.0) 0 else ceil(log10(abs(value))).toInt() - val order = if( + val approxOrder = if (value == 0.0) 0 else ceil(log10(abs(value))).toInt() + val order = if ( ((value % ten) == 0.0) or (value == 1.0) or - ((1/value) % ten == 0.0)) approxOrder else approxOrder - 1 + ((1 / value) % ten == 0.0) + ) approxOrder else approxOrder - 1 val lead = value / ten.pow(order) - val leadDisplay = round(lead*ten.pow(digits)) / ten.pow(digits) - val orderDisplay = if(order == 0) "" else if(order > 0) "E+$order" else "E$order" + val leadDisplay = round(lead * ten.pow(digits)) / ten.pow(digits) + val orderDisplay = if (order == 0) "" else if (order > 0) "E+$order" else "E$order" val valueDisplay = "$leadDisplay$orderDisplay" - val res = if(value < 0.0) valueDisplay else " $valueDisplay" + val res = if (value < 0.0) valueDisplay else " $valueDisplay" val fLength = digits + 6 val endSpace = " ".repeat(fLength - res.length) return "$res$endSpace" } -internal inline fun DoubleTensor.toPrettyString(): String = buildString { +internal fun DoubleTensor.toPrettyString(): String = buildString { var offset = 0 val shape = this@toPrettyString.shape val linearStructure = this@toPrettyString.linearStructure @@ -72,32 +73,36 @@ internal inline fun DoubleTensor.toPrettyString(): String = buildString { append(initString) var charOffset = 3 for (vector in vectorSequence()) { - append(" ".repeat(charOffset)) + repeat(charOffset) { append(' ') } val index = linearStructure.index(offset) for (ind in index.reversed()) { if (ind != 0) { break } - append("[") + append('[') charOffset += 1 } val values = vector.as1D().toMutableList().map(::format) - append(values.joinToString(", ")) - append("]") + values.joinTo(this, separator = ", ") + + append(']') charOffset -= 1 - for ((ind, maxInd) in index.reversed().zip(shape.reversed()).drop(1)){ + + index.reversed().zip(shape.reversed()).drop(1).forEach { (ind, maxInd) -> if (ind != maxInd - 1) { - break + return@forEach } - append("]") - charOffset -=1 + append(']') + charOffset -= 1 } + offset += vectorSize if (this@toPrettyString.numElements == offset) { break } + append(",\n") } append("\n)") diff --git a/kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleLinearOpsAlgebra.kt b/kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleLinearOpsAlgebra.kt index 1148c0aad..2282d7fcb 100644 --- a/kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleLinearOpsAlgebra.kt +++ b/kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleLinearOpsAlgebra.kt @@ -182,7 +182,7 @@ internal class TestDoubleLinearOpsTensorAlgebra { } -private inline fun DoubleLinearOpsTensorAlgebra.testSVDFor(tensor: DoubleTensor, epsilon: Double = 1e-10): Unit { +private fun DoubleLinearOpsTensorAlgebra.testSVDFor(tensor: DoubleTensor, epsilon: Double = 1e-10): Unit { val svd = tensor.svd() val tensorSVD = svd.first