kscience/kmath
9
0
Fork 1
Code Issues Pull Requests 6 Packages Projects Releases 16 Wiki Activity

refactoring #303

Merged
AndreiKingsley merged 2 commits from andrew into feature/tensor-algebra 2021-04-30 21:09:54 +03:00
Conversation 0 Commits 2 Files Changed 6 +123 -111
6 changed files with 123 additions and 111 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

51
kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/algebras/DoubleLinearOpsTensorAlgebra.kt
View File

@ -53,14 +53,15 @@ public object DoubleLinearOpsTensorAlgebra :
val lTensor = luTensor.zeroesLike()
val uTensor = luTensor.zeroesLike()
for ((pairLU, lu) in lTensor.matrixSequence().zip(uTensor.matrixSequence())
.zip(luTensor.tensor.matrixSequence())) {
val (l, u) = pairLU
luPivotHelper(l.as2D(), u.as2D(), lu.as2D(), n)
}
lTensor.matrixSequence()
.zip(uTensor.matrixSequence())
.zip(luTensor.tensor.matrixSequence())
.forEach { (pairLU, lu) ->
val (l, u) = pairLU
luPivotHelper(l.as2D(), u.as2D(), lu.as2D(), n)
}
return Triple(pTensor, lTensor, uTensor)
}
public fun TensorStructure<Double>.cholesky(epsilon: Double): DoubleTensor {
@ -82,11 +83,13 @@ public object DoubleLinearOpsTensorAlgebra :
checkSquareMatrix(shape)
val qTensor = zeroesLike()
val rTensor = zeroesLike()
val seq = tensor.matrixSequence().zip((qTensor.matrixSequence().zip(rTensor.matrixSequence())))
for ((matrix, qr) in seq) {
tensor.matrixSequence()
.zip((qTensor.matrixSequence()
.zip(rTensor.matrixSequence()))).forEach { (matrix, qr) ->
val (q, r) = qr
qrHelper(matrix.asTensor(), q.asTensor(), r.as2D())
}
return qTensor to rTensor
}
@ -97,20 +100,24 @@ public object DoubleLinearOpsTensorAlgebra :
val size = tensor.linearStructure.dim
val commonShape = tensor.shape.sliceArray(0 until size - 2)
val (n, m) = tensor.shape.sliceArray(size - 2 until size)
val resU = zeros(commonShape + intArrayOf(min(n, m), n))
val resS = zeros(commonShape + intArrayOf(min(n, m)))
val resV = zeros(commonShape + intArrayOf(min(n, m), m))
val uTensor = zeros(commonShape + intArrayOf(min(n, m), n))
val sTensor = zeros(commonShape + intArrayOf(min(n, m)))
val vTensor = zeros(commonShape + intArrayOf(min(n, m), m))
for ((matrix, USV) in tensor.matrixSequence()
.zip(resU.matrixSequence().zip(resS.vectorSequence().zip(resV.matrixSequence())))) {
val matrixSize = matrix.shape.reduce { acc, i -> acc * i }
val curMatrix = DoubleTensor(
matrix.shape,
matrix.mutableBuffer.array().slice(matrix.bufferStart until matrix.bufferStart + matrixSize).toDoubleArray()
)
svdHelper(curMatrix, USV, m, n, epsilon)
}
return Triple(resU.transpose(), resS, resV.transpose())
tensor.matrixSequence()
.zip(uTensor.matrixSequence()
.zip(sTensor.vectorSequence()
.zip(vTensor.matrixSequence()))).forEach { (matrix, USV) ->
val matrixSize = matrix.shape.reduce { acc, i -> acc * i }
val curMatrix = DoubleTensor(
matrix.shape,
matrix.mutableBuffer.array().slice(matrix.bufferStart until matrix.bufferStart + matrixSize)
.toDoubleArray()
)
svdHelper(curMatrix, USV, m, n, epsilon)
}
return Triple(uTensor.transpose(), sTensor, vTensor.transpose())
}
override fun TensorStructure<Double>.symEig(): Pair<DoubleTensor, DoubleTensor> =
@ -127,7 +134,7 @@ public object DoubleLinearOpsTensorAlgebra :
cleanSymHelper(matrix.as2D(), n)
val eig = (utv dot s.view(shp)).view(s.shape)
return Pair(eig, v)
return eig to v
}
public fun TensorStructure<Double>.detLU(epsilon: Double = 1e-9): DoubleTensor {

77
kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/broadcastUtils.kt
View File

@ -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<DoubleTensor> {
internal fun broadcastTensors(vararg tensors: DoubleTensor): List<DoubleTensor> {
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<DoubleT
return res
}
internal inline fun broadcastOuterTensors(vararg tensors: DoubleTensor): List<DoubleTensor> {
internal fun broadcastOuterTensors(vararg tensors: DoubleTensor): List<DoubleTensor> {
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<Do
val totalShape = broadcastShapes(*(tensors.map { it.shape.sliceArray(0..it.shape.size - 3) }).toTypedArray())
val n = totalShape.reduce { acc, i -> acc * i }
val res = ArrayList<DoubleTensor>(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
}

18
kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/checks.kt
View File

@ -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 <T> checkShapesCompatible(a: TensorStructure<T>, b: TensorStructure<T>): Unit =
internal fun <T> checkShapesCompatible(a: TensorStructure<T>, b: TensorStructure<T>): 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 <T> checkView(a: TensorStructure<T>, shape: IntArray): Unit =
internal fun <T> checkView(a: TensorStructure<T>, 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<Double>, 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)

41
kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/linUtils.kt
View File

@ -13,7 +13,7 @@ import kotlin.math.sign
import kotlin.math.sqrt
internal inline fun <T> BufferedTensor<T>.vectorSequence(): Sequence<BufferedTensor<T>> = sequence {
internal fun <T> BufferedTensor<T>.vectorSequence(): Sequence<BufferedTensor<T>> = sequence {
val n = shape.size
val vectorOffset = shape[n - 1]
val vectorShape = intArrayOf(shape.last())
@ -23,9 +23,9 @@ internal inline fun <T> BufferedTensor<T>.vectorSequence(): Sequence<BufferedTen
}
}
internal inline fun <T> BufferedTensor<T>.matrixSequence(): Sequence<BufferedTensor<T>> = sequence {
check(shape.size >= 2) { "todo" }
internal fun <T> BufferedTensor<T>.matrixSequence(): Sequence<BufferedTensor<T>> = 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 <T> BufferedTensor<T>.forEachMatrix(matrixAction: (BufferedT
}
}
internal inline fun dotHelper(
internal fun dotHelper(
a: MutableStructure2D<Double>,
b: MutableStructure2D<Double>,
res: MutableStructure2D<Double>,
@ -64,10 +63,11 @@ internal inline fun dotHelper(
}
}
internal inline fun luHelper(
internal fun luHelper(
lu: MutableStructure2D<Double>,
pivots: MutableStructure1D<Int>,
epsilon: Double): Boolean {
epsilon: Double
): Boolean {
val m = lu.rowNum
@ -114,7 +114,7 @@ internal inline fun luHelper(
return false
}
internal inline fun <T> BufferedTensor<T>.setUpPivots(): IntTensor {
internal fun <T> BufferedTensor<T>.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 <T> BufferedTensor<T>.setUpPivots(): IntTensor {
)
}
internal inline fun DoubleLinearOpsTensorAlgebra.computeLU(
internal fun DoubleLinearOpsTensorAlgebra.computeLU(
tensor: DoubleTensor,
epsilon: Double): Pair<DoubleTensor, IntTensor>? {
epsilon: Double
): Pair<DoubleTensor, IntTensor>? {
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<Double>,
pivot: MutableStructure1D<Int>,
n: Int
@ -151,7 +152,7 @@ internal inline fun pivInit(
}
}
internal inline fun luPivotHelper(
internal fun luPivotHelper(
l: MutableStructure2D<Double>,
u: MutableStructure2D<Double>,
lu: MutableStructure2D<Double>,
@ -172,7 +173,7 @@ internal inline fun luPivotHelper(
}
}
internal inline fun choleskyHelper(
internal fun choleskyHelper(
a: MutableStructure2D<Double>,
l: MutableStructure2D<Double>,
n: Int
@ -193,7 +194,7 @@ internal inline fun choleskyHelper(
}
}
internal inline fun luMatrixDet(lu: MutableStructure2D<Double>, pivots: MutableStructure1D<Int>): Double {
internal fun luMatrixDet(lu: MutableStructure2D<Double>, pivots: MutableStructure1D<Int>): Double {
if (lu[0, 0] == 0.0) {
return 0.0
}
@ -202,7 +203,7 @@ internal inline fun luMatrixDet(lu: MutableStructure2D<Double>, 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<Double>,
pivots: MutableStructure1D<Int>,
invMatrix: MutableStructure2D<Double>
@ -229,7 +230,7 @@ internal inline fun luMatrixInv(
}
}
internal inline fun DoubleLinearOpsTensorAlgebra.qrHelper(
internal fun DoubleLinearOpsTensorAlgebra.qrHelper(
matrix: DoubleTensor,
q: DoubleTensor,
r: MutableStructure2D<Double>
@ -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<BufferedTensor<Double>, Pair<BufferedTensor<Double>, BufferedTensor<Double>>>,
m: Int, n: Int, epsilon: Double
@ -335,7 +336,7 @@ internal inline fun DoubleLinearOpsTensorAlgebra.svdHelper(
}
}
internal inline fun cleanSymHelper(matrix: MutableStructure2D<Double>, n: Int) {
internal fun cleanSymHelper(matrix: MutableStructure2D<Double>, n: Int) {
for (i in 0 until n)
for (j in 0 until n) {
if (i == j) {

45
kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/utils.kt
View File

@ -23,19 +23,19 @@ internal fun Buffer<Double>.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 <T> BufferedTensor<T>.minusIndex(i: Int): Int = minusIndexFrom(this.dimension, i)
internal fun <T> BufferedTensor<T>.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)")

2
kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleLinearOpsAlgebra.kt
View File

@ -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