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

Inv and det #262

Merged
AndreiKingsley merged 5 commits from andrew into feature/tensor-algebra 2021-03-24 18:52:12 +03:00
Conversation 0 Commits 5 Files Changed 5 +185 -14
5 changed files with 185 additions and 14 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/LinearOpsTensorAlgebra.kt
View File

@ -14,7 +14,7 @@ public interface LinearOpsTensorAlgebra<T, TensorType : TensorStructure<T>, Inde
public fun TensorType.qr(): TensorType
//https://pytorch.org/docs/stable/generated/torch.lu.html
public fun TensorType.lu(tol: T): Pair<TensorType, IndexTensorType>
public fun TensorType.lu(): Pair<TensorType, IndexTensorType>
//https://pytorch.org/docs/stable/generated/torch.lu_unpack.html
public fun luPivot(luTensor: TensorType, pivotsTensor: IndexTensorType): Triple<TensorType, TensorType, TensorType>

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

@ -1,5 +1,8 @@
package space.kscience.kmath.tensors.core
import space.kscience.kmath.nd.MutableStructure2D
import space.kscience.kmath.nd.Structure1D
import space.kscience.kmath.nd.Structure2D
import space.kscience.kmath.tensors.LinearOpsTensorAlgebra
import kotlin.math.sqrt
@ -7,11 +10,11 @@ public class DoubleLinearOpsTensorAlgebra :
LinearOpsTensorAlgebra<Double, DoubleTensor, IntTensor>,
DoubleTensorAlgebra() {
override fun DoubleTensor.inv(): DoubleTensor {
TODO("ANDREI")
}
override fun DoubleTensor.inv(): DoubleTensor = invLU()
override fun DoubleTensor.lu(tol: Double): Pair<DoubleTensor, IntTensor> {
override fun DoubleTensor.det(): DoubleTensor = detLU()
override fun DoubleTensor.lu(): Pair<DoubleTensor, IntTensor> {
checkSquareMatrix(shape)
@ -20,6 +23,8 @@ public class DoubleLinearOpsTensorAlgebra :
val n = shape.size
val m = shape.last()
val pivotsShape = IntArray(n - 1) { i -> shape[i] }
pivotsShape[n - 2] = m + 1
val pivotsTensor = IntTensor(
pivotsShape,
IntArray(pivotsShape.reduce(Int::times)) { 0 }
@ -54,6 +59,8 @@ public class DoubleLinearOpsTensorAlgebra :
lu[maxInd, k] = tmp
}
pivots[m] += 1
}
for (j in i + 1 until m) {
@ -138,7 +145,6 @@ public class DoubleLinearOpsTensorAlgebra :
TODO("ANDREI")
}
override fun DoubleTensor.svd(): Triple<DoubleTensor, DoubleTensor, DoubleTensor> {
TODO("ALYA")
}
@ -147,6 +153,72 @@ public class DoubleLinearOpsTensorAlgebra :
TODO("ANDREI")
}
private fun luMatrixDet(lu: Structure2D<Double>, pivots: Structure1D<Int>): Double {
val m = lu.shape[0]
val sign = if((pivots[m] - m) % 2 == 0) 1.0 else -1.0
return (0 until m).asSequence().map { lu[it, it] }.fold(sign) { left, right -> left * right }
}
public fun DoubleTensor.detLU(): DoubleTensor {
val (luTensor, pivotsTensor) = lu()
val n = shape.size
val detTensorShape = IntArray(n - 1) { i -> shape[i] }
detTensorShape[n - 2] = 1
val resBuffer = DoubleArray(detTensorShape.reduce(Int::times)) { 0.0 }
val detTensor = DoubleTensor(
detTensorShape,
resBuffer
)
luTensor.matrixSequence().zip(pivotsTensor.vectorSequence()).forEachIndexed { index, (luMatrix, pivots) ->
resBuffer[index] = luMatrixDet(luMatrix, pivots)
}
return detTensor
}
private fun luMatrixInv(
lu: Structure2D<Double>,
pivots: Structure1D<Int>,
invMatrix : MutableStructure2D<Double>
): Unit {
val m = lu.shape[0]
for (j in 0 until m) {
for (i in 0 until m) {
if (pivots[i] == j){
invMatrix[i, j] = 1.0
}
for (k in 0 until i){
invMatrix[i, j] -= lu[i, k] * invMatrix[k, j]
}
}
for (i in m - 1 downTo 0) {
for (k in i + 1 until m) {
invMatrix[i, j] -= lu[i, k] * invMatrix[k, j]
}
invMatrix[i, j] /= lu[i, i]
}
}
}
public fun DoubleTensor.invLU(): DoubleTensor {
val (luTensor, pivotsTensor) = lu()
val n = shape.size
val invTensor = luTensor.zeroesLike()
val seq = luTensor.matrixSequence().zip(pivotsTensor.vectorSequence()).zip(invTensor.matrixSequence())
for ((luP, invMatrix) in seq) {
val (lu, pivots) = luP
luMatrixInv(lu, pivots, invMatrix)
}
return invTensor
}
}
public inline fun <R> DoubleLinearOpsTensorAlgebra(block: DoubleLinearOpsTensorAlgebra.() -> R): R =

29
kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleTensorAlgebra.kt
View File

@ -1,5 +1,8 @@
package space.kscience.kmath.tensors.core
import space.kscience.kmath.linear.Matrix
import space.kscience.kmath.nd.MutableStructure2D
import space.kscience.kmath.nd.Structure2D
import space.kscience.kmath.tensors.TensorPartialDivisionAlgebra
import kotlin.math.abs
@ -318,7 +321,31 @@ public open class DoubleTensorAlgebra : TensorPartialDivisionAlgebra<Double, Dou
}
override fun DoubleTensor.det(): DoubleTensor {
TODO("ANDREI")
TODO()
/*
checkSquareMatrix(shape)
val n = shape.size
val m = shape.last()
val detTensorShape = IntArray(n - 1) { i -> shape[i] }
detTensorShape[n - 1] = 1
val resBuffer = DoubleArray(detTensorShape.reduce(Int::times)) { 0.0 }
val detTensor = DoubleTensor(
detTensorShape,
resBuffer
)
this.matrixSequence().forEachIndexed{i, matrix ->
// todo need Matrix determinant algo
// todo resBuffer[i] = matrix.det()
}
return detTensor
*/
}
override fun DoubleTensor.square(): DoubleTensor {

7
kmath-core/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleAnalyticTensorAlgebra.kt
View File

@ -3,7 +3,6 @@ package space.kscience.kmath.tensors.core
import kotlin.math.abs
import kotlin.math.exp
import kotlin.test.Test
import kotlin.test.assertEquals
import kotlin.test.assertTrue
class TestDoubleAnalyticTensorAlgebra {
@ -16,9 +15,9 @@ class TestDoubleAnalyticTensorAlgebra {
return this.map(transform).toDoubleArray()
}
fun DoubleArray.deltaEqual(other: DoubleArray, delta: Double = 1e-5): Boolean {
fun DoubleArray.epsEqual(other: DoubleArray, eps: Double = 1e-5): Boolean {
for ((elem1, elem2) in this.asSequence().zip(other.asSequence())) {
if (abs(elem1 - elem2) > delta) {
if (abs(elem1 - elem2) > eps) {
return false
}
}
@ -29,7 +28,7 @@ class TestDoubleAnalyticTensorAlgebra {
fun testExp() = DoubleAnalyticTensorAlgebra {
tensor.exp().let {
assertTrue { shape contentEquals it.shape }
assertTrue { buffer.fmap(::exp).deltaEqual(it.buffer.array())}
assertTrue { buffer.fmap(::exp).epsEqual(it.buffer.array())}
}
}
}

73
kmath-core/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleLinearOpsAlgebra.kt Normal file
View File

@ -0,0 +1,73 @@
package space.kscience.kmath.tensors.core
import kotlin.math.abs
import kotlin.math.exp
import kotlin.test.Test
import kotlin.test.assertEquals
import kotlin.test.assertTrue
class TestDoubleLinearOpsTensorAlgebra {
private val eps = 1e-5
private fun Double.epsEqual(other: Double): Boolean {
return abs(this - other) < eps
}
fun DoubleArray.epsEqual(other: DoubleArray, eps: Double = 1e-5): Boolean {
for ((elem1, elem2) in this.asSequence().zip(other.asSequence())) {
if (abs(elem1 - elem2) > eps) {
return false
}
}
return true
}
@Test
fun testDetLU() = DoubleLinearOpsTensorAlgebra {
val tensor = fromArray(
intArrayOf(2, 2, 2),
doubleArrayOf(
1.0, 3.0,
1.0, 2.0,
1.5, 1.0,
10.0, 2.0
)
)
val expectedShape = intArrayOf(2, 1)
val expectedBuffer = doubleArrayOf(
-1.0,
-7.0
)
val detTensor = tensor.detLU()
assertTrue { detTensor.shape contentEquals expectedShape }
assertTrue { detTensor.buffer.array().epsEqual(expectedBuffer) }
}
@Test
fun testInvLU() = DoubleLinearOpsTensorAlgebra {
val tensor = fromArray(
intArrayOf(2, 2, 2),
doubleArrayOf(
1.0, 0.0,
0.0, 2.0,
1.0, 1.0,
1.0, 0.0
)
)
val expectedShape = intArrayOf(2, 2, 2)
val expectedBuffer = doubleArrayOf(
1.0, 0.0,
0.0, 0.5,
0.0, 1.0,
1.0, -1.0
)
val invTensor = tensor.invLU()
assertTrue { invTensor.shape contentEquals expectedShape }
assertTrue { invTensor.buffer.array().epsEqual(expectedBuffer) }
}
}