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

v0.3.0-dev-9 #324

Merged
altavir merged 265 commits from dev into master 2021-05-08 17:16:29 +03:00
Conversation 0 Commits 265 Files Changed 182 +94 -89
5 changed files with 94 additions and 89 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files
Showing only changes of commit f0cdb9b657 - Show all commits

17
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 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) { for (linearIndex in 0 until linearSize) {
val totalMultiIndex = resTensor.linearStructure.index(linearIndex) val totalMultiIndex = resTensor.linearStructure.index(linearIndex)
val curMultiIndex = tensor.shape.copyOf() 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 var totalDim = 0
for (shape in shapes) { for (shape in shapes) {
totalDim = max(totalDim, shape.size) totalDim = max(totalDim, shape.size)
@ -51,7 +51,7 @@ internal inline fun broadcastShapes(vararg shapes: IntArray): IntArray {
return totalShape 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) { if (tensor.shape.size > newShape.size) {
throw RuntimeException("Tensor is not compatible with the new shape") 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 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 totalShape = broadcastShapes(*(tensors.map { it.shape }).toTypedArray())
val n = totalShape.reduce { acc, i -> acc * i } val n = totalShape.reduce { acc, i -> acc * i }
@ -85,7 +85,7 @@ internal inline fun broadcastTensors(vararg tensors: DoubleTensor): List<DoubleT
return res return res
} }
internal inline fun broadcastOuterTensors(vararg tensors: DoubleTensor): List<DoubleTensor> { internal fun broadcastOuterTensors(vararg tensors: DoubleTensor): List<DoubleTensor> {
val onlyTwoDims = tensors.asSequence().onEach { val onlyTwoDims = tensors.asSequence().onEach {
require(it.shape.size >= 2) { require(it.shape.size >= 2) {
throw RuntimeException("Tensors must have at least 2 dimensions") throw RuntimeException("Tensors must have at least 2 dimensions")
@ -99,7 +99,7 @@ internal inline fun broadcastOuterTensors(vararg tensors: DoubleTensor): List<Do
val totalShape = broadcastShapes(*(tensors.map { it.shape.sliceArray(0..it.shape.size - 3) }).toTypedArray()) val totalShape = broadcastShapes(*(tensors.map { it.shape.sliceArray(0..it.shape.size - 3) }).toTypedArray())
val n = totalShape.reduce { acc, i -> acc * i } val n = totalShape.reduce { acc, i -> acc * i }
val res = ArrayList<DoubleTensor>(0) return buildList {
for (tensor in tensors) { for (tensor in tensors) {
val matrixShape = tensor.shape.sliceArray(tensor.shape.size - 2 until tensor.shape.size).copyOf() val matrixShape = tensor.shape.sliceArray(tensor.shape.size - 2 until tensor.shape.size).copyOf()
val matrixSize = matrixShape[0] * matrixShape[1] val matrixSize = matrixShape[0] * matrixShape[1]
@ -137,8 +137,7 @@ internal inline fun broadcastOuterTensors(vararg tensors: DoubleTensor): List<Do
newTensor.mutableBuffer.array()[newTensor.bufferStart + curLinearIndex] newTensor.mutableBuffer.array()[newTensor.bufferStart + curLinearIndex]
} }
} }
res += resTensor add(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 import space.kscience.kmath.tensors.core.algebras.DoubleTensorAlgebra
internal inline fun checkEmptyShape(shape: IntArray): Unit = internal fun checkEmptyShape(shape: IntArray): Unit =
check(shape.isNotEmpty()) { check(shape.isNotEmpty()) {
"Illegal empty shape provided" "Illegal empty shape provided"
} }
internal inline fun checkEmptyDoubleBuffer(buffer: DoubleArray): Unit = internal fun checkEmptyDoubleBuffer(buffer: DoubleArray): Unit =
check(buffer.isNotEmpty()) { check(buffer.isNotEmpty()) {
"Illegal empty buffer provided" "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)) { check(buffer.size == shape.reduce(Int::times)) {
"Inconsistent shape ${shape.toList()} for buffer of size ${buffer.size} provided" "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) { check(a.shape contentEquals b.shape) {
"Incompatible shapes ${a.shape.toList()} and ${b.shape.toList()} " "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)) { check((i < dim) and (j < dim)) {
"Cannot transpose $i to $j for a tensor of dim $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)) 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 val n = shape.size
check(n >= 2) { check(n >= 2) {
"Expected tensor with 2 or more dimensions, got size $n instead" "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 tensor: TensorStructure<Double>, epsilon: Double = 1e-6
): Unit = ): Unit =
check(tensor.eq(tensor.transpose(), epsilon)) { check(tensor.eq(tensor.transpose(), epsilon)) {
"Tensor is not symmetric about the last 2 dimensions at precision $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 tensor: DoubleTensor, epsilon: Double = 1e-6
): Unit { ): Unit {
checkSymmetric(tensor, epsilon) 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 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 n = shape.size
val vectorOffset = shape[n - 1] val vectorOffset = shape[n - 1]
val vectorShape = intArrayOf(shape.last()) 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 { internal fun <T> BufferedTensor<T>.matrixSequence(): Sequence<BufferedTensor<T>> = sequence {
check(shape.size >= 2) { "todo" }
val n = shape.size 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 matrixOffset = shape[n - 1] * shape[n - 2]
val matrixShape = intArrayOf(shape[n - 2], shape[n - 1]) val matrixShape = intArrayOf(shape[n - 2], shape[n - 1])
for (offset in 0 until numElements step matrixOffset) { for (offset in 0 until numElements step matrixOffset) {
@ -46,8 +46,7 @@ internal inline fun <T> BufferedTensor<T>.forEachMatrix(matrixAction: (BufferedT
} }
} }
internal fun dotHelper(
internal inline fun dotHelper(
a: MutableStructure2D<Double>, a: MutableStructure2D<Double>,
b: MutableStructure2D<Double>, b: MutableStructure2D<Double>,
res: MutableStructure2D<Double>, res: MutableStructure2D<Double>,
@ -64,10 +63,11 @@ internal inline fun dotHelper(
} }
} }
internal inline fun luHelper( internal fun luHelper(
lu: MutableStructure2D<Double>, lu: MutableStructure2D<Double>,
pivots: MutableStructure1D<Int>, pivots: MutableStructure1D<Int>,
epsilon: Double): Boolean { epsilon: Double
): Boolean {
val m = lu.rowNum val m = lu.rowNum
@ -114,7 +114,7 @@ internal inline fun luHelper(
return false return false
} }
internal inline fun <T> BufferedTensor<T>.setUpPivots(): IntTensor { internal fun <T> BufferedTensor<T>.setUpPivots(): IntTensor {
val n = this.shape.size val n = this.shape.size
val m = this.shape.last() val m = this.shape.last()
val pivotsShape = IntArray(n - 1) { i -> this.shape[i] } 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, tensor: DoubleTensor,
epsilon: Double): Pair<DoubleTensor, IntTensor>? { epsilon: Double
): Pair<DoubleTensor, IntTensor>? {
checkSquareMatrix(tensor.shape) checkSquareMatrix(tensor.shape)
val luTensor = tensor.copy() val luTensor = tensor.copy()
val pivotsTensor = tensor.setUpPivots() val pivotsTensor = tensor.setUpPivots()
for ((lu, pivots) in luTensor.matrixSequence().zip(pivotsTensor.vectorSequence())) 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 null
return Pair(luTensor, pivotsTensor) return Pair(luTensor, pivotsTensor)
} }
internal inline fun pivInit( internal fun pivInit(
p: MutableStructure2D<Double>, p: MutableStructure2D<Double>,
pivot: MutableStructure1D<Int>, pivot: MutableStructure1D<Int>,
n: Int n: Int
@ -151,7 +152,7 @@ internal inline fun pivInit(
} }
} }
internal inline fun luPivotHelper( internal fun luPivotHelper(
l: MutableStructure2D<Double>, l: MutableStructure2D<Double>,
u: MutableStructure2D<Double>, u: MutableStructure2D<Double>,
lu: MutableStructure2D<Double>, lu: MutableStructure2D<Double>,
@ -172,7 +173,7 @@ internal inline fun luPivotHelper(
} }
} }
internal inline fun choleskyHelper( internal fun choleskyHelper(
a: MutableStructure2D<Double>, a: MutableStructure2D<Double>,
l: MutableStructure2D<Double>, l: MutableStructure2D<Double>,
n: Int 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) { if (lu[0, 0] == 0.0) {
return 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 } 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>, lu: MutableStructure2D<Double>,
pivots: MutableStructure1D<Int>, pivots: MutableStructure1D<Int>,
invMatrix: MutableStructure2D<Double> invMatrix: MutableStructure2D<Double>
@ -229,7 +230,7 @@ internal inline fun luMatrixInv(
} }
} }
internal inline fun DoubleLinearOpsTensorAlgebra.qrHelper( internal fun DoubleLinearOpsTensorAlgebra.qrHelper(
matrix: DoubleTensor, matrix: DoubleTensor,
q: DoubleTensor, q: DoubleTensor,
r: MutableStructure2D<Double> 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 val (n, m) = a.shape
var v: DoubleTensor var v: DoubleTensor
val b: 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, matrix: DoubleTensor,
USV: Pair<BufferedTensor<Double>, Pair<BufferedTensor<Double>, BufferedTensor<Double>>>, USV: Pair<BufferedTensor<Double>, Pair<BufferedTensor<Double>, BufferedTensor<Double>>>,
m: Int, n: Int, epsilon: 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 (i in 0 until n)
for (j in 0 until n) { for (j in 0 until n) {
if (i == j) { 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() 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 distribution = GaussianSampler(0.0, 1.0)
val generator = RandomGenerator.default(seed) val generator = RandomGenerator.default(seed)
return distribution.sample(generator).nextBufferBlocking(n).toDoubleArray() 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 unnorm = getRandomNormals(n, seed)
val norm = sqrt(unnorm.map { it * it }.sum()) val norm = sqrt(unnorm.map { it * it }.sum())
return unnorm.map { it / norm }.toDoubleArray() 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 val ii = n + i
check(ii >= 0) { check(ii >= 0) {
"Out of bound index $i for tensor of dim $n" "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 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 ten = 10.0
val approxOrder = if(value == 0.0) 0 else ceil(log10(abs(value))).toInt() val approxOrder = if (value == 0.0) 0 else ceil(log10(abs(value))).toInt()
val order = if( val order = if (
((value % ten) == 0.0) or ((value % ten) == 0.0) or
(value == 1.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 lead = value / ten.pow(order)
val leadDisplay = round(lead*ten.pow(digits)) / ten.pow(digits) 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 orderDisplay = if (order == 0) "" else if (order > 0) "E+$order" else "E$order"
val valueDisplay = "$leadDisplay$orderDisplay" 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 fLength = digits + 6
val endSpace = " ".repeat(fLength - res.length) val endSpace = " ".repeat(fLength - res.length)
return "$res$endSpace" return "$res$endSpace"
} }
internal inline fun DoubleTensor.toPrettyString(): String = buildString { internal fun DoubleTensor.toPrettyString(): String = buildString {
var offset = 0 var offset = 0
val shape = this@toPrettyString.shape val shape = this@toPrettyString.shape
val linearStructure = this@toPrettyString.linearStructure val linearStructure = this@toPrettyString.linearStructure
@ -72,32 +73,36 @@ internal inline fun DoubleTensor.toPrettyString(): String = buildString {
append(initString) append(initString)
var charOffset = 3 var charOffset = 3
for (vector in vectorSequence()) { for (vector in vectorSequence()) {
append(" ".repeat(charOffset)) repeat(charOffset) { append(' ') }
val index = linearStructure.index(offset) val index = linearStructure.index(offset)
for (ind in index.reversed()) { for (ind in index.reversed()) {
if (ind != 0) { if (ind != 0) {
break break
} }
append("[") append('[')
charOffset += 1 charOffset += 1
} }
val values = vector.as1D().toMutableList().map(::format) val values = vector.as1D().toMutableList().map(::format)
append(values.joinToString(", ")) values.joinTo(this, separator = ", ")
append("]")
append(']')
charOffset -= 1 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) { if (ind != maxInd - 1) {
break return@forEach
} }
append("]") append(']')
charOffset -=1 charOffset -= 1
} }
offset += vectorSize offset += vectorSize
if (this@toPrettyString.numElements == offset) { if (this@toPrettyString.numElements == offset) {
break break
} }
append(",\n") append(",\n")
} }
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 svd = tensor.svd()
val tensorSVD = svd.first val tensorSVD = svd.first