Adding more systematic checks

This commit is contained in:
rgrit91 2021-01-16 09:54:03 +00:00
commit ca2082405a
11 changed files with 234 additions and 68 deletions

View File

@ -7,10 +7,16 @@ import kscience.kmath.structures.asBuffer
import kotlin.math.sqrt import kotlin.math.sqrt
/** /**
* A 2d structure plus optional matrix-specific features * A [Matrix] that holds [MatrixFeature] objects.
*
* @param T the type of items.
*/ */
public interface FeaturedMatrix<T : Any> : Matrix<T> { public interface FeaturedMatrix<T : Any> : Matrix<T> {
override val shape: IntArray get() = intArrayOf(rowNum, colNum) public override val shape: IntArray get() = intArrayOf(rowNum, colNum)
/**
* The set of features this matrix possesses.
*/
public val features: Set<MatrixFeature> public val features: Set<MatrixFeature>
/** /**

View File

@ -4,16 +4,15 @@ import kscience.kmath.operations.*
import kscience.kmath.structures.* import kscience.kmath.structures.*
/** /**
* Common implementation of [LUPDecompositionFeature] * Common implementation of [LupDecompositionFeature].
*/ */
public class LUPDecomposition<T : Any>( public class LUPDecomposition<T : Any>(
public val context: MatrixContext<T, FeaturedMatrix<T>>, public val context: MatrixContext<T, FeaturedMatrix<T>>,
public val elementContext: Field<T>, public val elementContext: Field<T>,
public val lu: Structure2D<T>, public val lu: Matrix<T>,
public val pivot: IntArray, public val pivot: IntArray,
private val even: Boolean, private val even: Boolean,
) : LUPDecompositionFeature<T>, DeterminantFeature<T> { ) : LupDecompositionFeature<T>, DeterminantFeature<T> {
/** /**
* Returns the matrix L of the decomposition. * Returns the matrix L of the decomposition.
* *
@ -151,7 +150,10 @@ public inline fun <reified T : Comparable<T>, F : Field<T>> GenericMatrixContext
public fun MatrixContext<Double, FeaturedMatrix<Double>>.lup(matrix: Matrix<Double>): LUPDecomposition<Double> = public fun MatrixContext<Double, FeaturedMatrix<Double>>.lup(matrix: Matrix<Double>): LUPDecomposition<Double> =
lup(Buffer.Companion::real, RealField, matrix) { it < 1e-11 } lup(Buffer.Companion::real, RealField, matrix) { it < 1e-11 }
public fun <T : Any> LUPDecomposition<T>.solveWithLUP(factory: MutableBufferFactory<T>, matrix: Matrix<T>): FeaturedMatrix<T> { public fun <T : Any> LUPDecomposition<T>.solveWithLUP(
factory: MutableBufferFactory<T>,
matrix: Matrix<T>
): FeaturedMatrix<T> {
require(matrix.rowNum == pivot.size) { "Matrix dimension mismatch. Expected ${pivot.size}, but got ${matrix.colNum}" } require(matrix.rowNum == pivot.size) { "Matrix dimension mismatch. Expected ${pivot.size}, but got ${matrix.colNum}" }
BufferAccessor2D(matrix.rowNum, matrix.colNum, factory).run { BufferAccessor2D(matrix.rowNum, matrix.colNum, factory).run {

View File

@ -1,62 +1,154 @@
package kscience.kmath.linear package kscience.kmath.linear
/** /**
* A marker interface representing some matrix feature like diagonal, sparse, zero, etc. Features used to optimize matrix * A marker interface representing some properties of matrices or additional transformations of them. Features are used
* operations performance in some cases. * to optimize matrix operations performance in some cases or retrieve the APIs.
*/ */
public interface MatrixFeature public interface MatrixFeature
/** /**
* The matrix with this feature is considered to have only diagonal non-null elements * Matrices with this feature are considered to have only diagonal non-null elements.
*/ */
public object DiagonalFeature : MatrixFeature public object DiagonalFeature : MatrixFeature
/** /**
* Matrix with this feature has all zero elements * Matrices with this feature have all zero elements.
*/ */
public object ZeroFeature : MatrixFeature public object ZeroFeature : MatrixFeature
/** /**
* Matrix with this feature have unit elements on diagonal and zero elements in all other places * Matrices with this feature have unit elements on diagonal and zero elements in all other places.
*/ */
public object UnitFeature : MatrixFeature public object UnitFeature : MatrixFeature
/** /**
* Inverted matrix feature * Matrices with this feature can be inverted: [inverse] = `a`<sup>-1</sup> where `a` is the owning matrix.
*
* @param T the type of matrices' items.
*/ */
public interface InverseMatrixFeature<T : Any> : MatrixFeature { public interface InverseMatrixFeature<T : Any> : MatrixFeature {
/**
* The inverse matrix of the matrix that owns this feature.
*/
public val inverse: FeaturedMatrix<T> public val inverse: FeaturedMatrix<T>
} }
/** /**
* A determinant container * Matrices with this feature can compute their determinant.
*/ */
public interface DeterminantFeature<T : Any> : MatrixFeature { public interface DeterminantFeature<T : Any> : MatrixFeature {
/**
* The determinant of the matrix that owns this feature.
*/
public val determinant: T public val determinant: T
} }
/**
* Produces a [DeterminantFeature] where the [DeterminantFeature.determinant] is [determinant].
*
* @param determinant the value of determinant.
* @return a new [DeterminantFeature].
*/
@Suppress("FunctionName") @Suppress("FunctionName")
public fun <T : Any> DeterminantFeature(determinant: T): DeterminantFeature<T> = object : DeterminantFeature<T> { public fun <T : Any> DeterminantFeature(determinant: T): DeterminantFeature<T> = object : DeterminantFeature<T> {
override val determinant: T = determinant override val determinant: T = determinant
} }
/** /**
* Lower triangular matrix * Matrices with this feature are lower triangular ones.
*/ */
public object LFeature : MatrixFeature public object LFeature : MatrixFeature
/** /**
* Upper triangular feature * Matrices with this feature are upper triangular ones.
*/ */
public object UFeature : MatrixFeature public object UFeature : MatrixFeature
/** /**
* TODO add documentation * Matrices with this feature support LU factorization with partial pivoting: *[p] &middot; a = [l] &middot; [u]* where
* *a* is the owning matrix.
*
* @param T the type of matrices' items.
*/ */
public interface LUPDecompositionFeature<T : Any> : MatrixFeature { public interface LupDecompositionFeature<T : Any> : MatrixFeature {
/**
* The lower triangular matrix in this decomposition. It may have [LFeature].
*/
public val l: FeaturedMatrix<T> public val l: FeaturedMatrix<T>
/**
* The upper triangular matrix in this decomposition. It may have [UFeature].
*/
public val u: FeaturedMatrix<T> public val u: FeaturedMatrix<T>
/**
* The permutation matrix in this decomposition.
*/
public val p: FeaturedMatrix<T> public val p: FeaturedMatrix<T>
} }
/**
* Matrices with this feature are orthogonal ones: *a &middot; a<sup>T</sup> = u* where *a* is the owning matrix, *u*
* is the unit matrix ([UnitFeature]).
*/
public object OrthogonalFeature : MatrixFeature
/**
* Matrices with this feature support QR factorization: *a = [q] &middot; [r]* where *a* is the owning matrix.
*
* @param T the type of matrices' items.
*/
public interface QRDecompositionFeature<T : Any> : MatrixFeature {
/**
* The orthogonal matrix in this decomposition. It may have [OrthogonalFeature].
*/
public val q: FeaturedMatrix<T>
/**
* The upper triangular matrix in this decomposition. It may have [UFeature].
*/
public val r: FeaturedMatrix<T>
}
/**
* Matrices with this feature support Cholesky factorization: *a = [l] &middot; [l]<sup>H</sup>* where *a* is the
* owning matrix.
*
* @param T the type of matrices' items.
*/
public interface CholeskyDecompositionFeature<T : Any> : MatrixFeature {
/**
* The triangular matrix in this decomposition. It may have either [UFeature] or [LFeature].
*/
public val l: FeaturedMatrix<T>
}
/**
* Matrices with this feature support SVD: *a = [u] &middot; [s] &middot; [v]<sup>H</sup>* where *a* is the owning
* matrix.
*
* @param T the type of matrices' items.
*/
public interface SingularValueDecompositionFeature<T : Any> : MatrixFeature {
/**
* The matrix in this decomposition. It is unitary, and it consists from left singular vectors.
*/
public val u: FeaturedMatrix<T>
/**
* The matrix in this decomposition. Its main diagonal elements are singular values.
*/
public val s: FeaturedMatrix<T>
/**
* The matrix in this decomposition. It is unitary, and it consists from right singular vectors.
*/
public val v: FeaturedMatrix<T>
/**
* The buffer of singular values of this SVD.
*/
public val singularValues: Point<T>
}
//TODO add sparse matrix feature //TODO add sparse matrix feature

View File

@ -1,12 +1,40 @@
package kscience.kmath.structures package kscience.kmath.structures
/** /**
* A structure that is guaranteed to be two-dimensional * A structure that is guaranteed to be two-dimensional.
*
* @param T the type of items.
*/ */
public interface Structure2D<T> : NDStructure<T> { public interface Structure2D<T> : NDStructure<T> {
/**
* The number of rows in this structure.
*/
public val rowNum: Int get() = shape[0] public val rowNum: Int get() = shape[0]
/**
* The number of columns in this structure.
*/
public val colNum: Int get() = shape[1] public val colNum: Int get() = shape[1]
/**
* The buffer of rows of this structure. It gets elements from the structure dynamically.
*/
public val rows: Buffer<Buffer<T>>
get() = VirtualBuffer(rowNum) { i -> VirtualBuffer(colNum) { j -> get(i, j) } }
/**
* The buffer of columns of this structure. It gets elements from the structure dynamically.
*/
public val columns: Buffer<Buffer<T>>
get() = VirtualBuffer(colNum) { j -> VirtualBuffer(rowNum) { i -> get(i, j) } }
/**
* Retrieves an element from the structure by two indices.
*
* @param i the first index.
* @param j the second index.
* @return an element.
*/
public operator fun get(i: Int, j: Int): T public operator fun get(i: Int, j: Int): T
override operator fun get(index: IntArray): T { override operator fun get(index: IntArray): T {
@ -14,15 +42,9 @@ public interface Structure2D<T> : NDStructure<T> {
return get(index[0], index[1]) return get(index[0], index[1])
} }
public val rows: Buffer<Buffer<T>>
get() = VirtualBuffer(rowNum) { i -> VirtualBuffer(colNum) { j -> get(i, j) } }
public val columns: Buffer<Buffer<T>>
get() = VirtualBuffer(colNum) { j -> VirtualBuffer(rowNum) { i -> get(i, j) } }
override fun elements(): Sequence<Pair<IntArray, T>> = sequence { override fun elements(): Sequence<Pair<IntArray, T>> = sequence {
for (i in (0 until rowNum)) for (i in 0 until rowNum)
for (j in (0 until colNum)) yield(intArrayOf(i, j) to get(i, j)) for (j in 0 until colNum) yield(intArrayOf(i, j) to get(i, j))
} }
public companion object public companion object
@ -47,4 +69,9 @@ public fun <T> NDStructure<T>.as2D(): Structure2D<T> = if (shape.size == 2)
else else
error("Can't create 2d-structure from ${shape.size}d-structure") error("Can't create 2d-structure from ${shape.size}d-structure")
/**
* Alias for [Structure2D] with more familiar name.
*
* @param T the type of items.
*/
public typealias Matrix<T> = Structure2D<T> public typealias Matrix<T> = Structure2D<T>

View File

@ -1,12 +1,10 @@
package kscience.kmath.ejml package kscience.kmath.ejml
import kscience.kmath.linear.*
import kscience.kmath.structures.NDStructure
import kscience.kmath.structures.RealBuffer
import org.ejml.dense.row.factory.DecompositionFactory_DDRM import org.ejml.dense.row.factory.DecompositionFactory_DDRM
import org.ejml.simple.SimpleMatrix import org.ejml.simple.SimpleMatrix
import kscience.kmath.linear.DeterminantFeature
import kscience.kmath.linear.FeaturedMatrix
import kscience.kmath.linear.LUPDecompositionFeature
import kscience.kmath.linear.MatrixFeature
import kscience.kmath.structures.NDStructure
/** /**
* Represents featured matrix over EJML [SimpleMatrix]. * Represents featured matrix over EJML [SimpleMatrix].
@ -14,42 +12,71 @@ import kscience.kmath.structures.NDStructure
* @property origin the underlying [SimpleMatrix]. * @property origin the underlying [SimpleMatrix].
* @author Iaroslav Postovalov * @author Iaroslav Postovalov
*/ */
public class EjmlMatrix(public val origin: SimpleMatrix, features: Set<MatrixFeature>? = null) : FeaturedMatrix<Double> { public class EjmlMatrix(public val origin: SimpleMatrix, features: Set<MatrixFeature> = emptySet()) :
FeaturedMatrix<Double> {
public override val rowNum: Int public override val rowNum: Int
get() = origin.numRows() get() = origin.numRows()
public override val colNum: Int public override val colNum: Int
get() = origin.numCols() get() = origin.numCols()
public override val shape: IntArray public override val shape: IntArray by lazy { intArrayOf(rowNum, colNum) }
get() = intArrayOf(origin.numRows(), origin.numCols())
public override val features: Set<MatrixFeature> = setOf( public override val features: Set<MatrixFeature> = hashSetOf(
object : LUPDecompositionFeature<Double>, DeterminantFeature<Double> { object : InverseMatrixFeature<Double> {
override val determinant: Double override val inverse: FeaturedMatrix<Double> by lazy { EjmlMatrix(origin.invert()) }
get() = origin.determinant() },
private val lup by lazy { object : DeterminantFeature<Double> {
val ludecompositionF64 = DecompositionFactory_DDRM.lu(origin.numRows(), origin.numCols()) override val determinant: Double by lazy(origin::determinant)
.also { it.decompose(origin.ddrm.copy()) } },
Triple( object : SingularValueDecompositionFeature<Double> {
EjmlMatrix(SimpleMatrix(ludecompositionF64.getRowPivot(null))), private val svd by lazy {
EjmlMatrix(SimpleMatrix(ludecompositionF64.getLower(null))), DecompositionFactory_DDRM.svd(origin.numRows(), origin.numCols(), true, true, false)
EjmlMatrix(SimpleMatrix(ludecompositionF64.getUpper(null))), .apply { decompose(origin.ddrm.copy()) }
)
} }
override val l: FeaturedMatrix<Double> override val u: FeaturedMatrix<Double> by lazy { EjmlMatrix(SimpleMatrix(svd.getU(null, false))) }
get() = lup.second override val s: FeaturedMatrix<Double> by lazy { EjmlMatrix(SimpleMatrix(svd.getW(null))) }
override val v: FeaturedMatrix<Double> by lazy { EjmlMatrix(SimpleMatrix(svd.getV(null, false))) }
override val singularValues: Point<Double> by lazy { RealBuffer(svd.singularValues) }
},
override val u: FeaturedMatrix<Double> object : QRDecompositionFeature<Double> {
get() = lup.third private val qr by lazy {
DecompositionFactory_DDRM.qr().apply { decompose(origin.ddrm.copy()) }
}
override val p: FeaturedMatrix<Double> override val q: FeaturedMatrix<Double> by lazy { EjmlMatrix(SimpleMatrix(qr.getQ(null, false))) }
get() = lup.first override val r: FeaturedMatrix<Double> by lazy { EjmlMatrix(SimpleMatrix(qr.getR(null, false))) }
} },
) union features.orEmpty()
object : CholeskyDecompositionFeature<Double> {
override val l: FeaturedMatrix<Double> by lazy {
val cholesky =
DecompositionFactory_DDRM.chol(rowNum, true).apply { decompose(origin.ddrm.copy()) }
EjmlMatrix(SimpleMatrix(cholesky.getT(null)), setOf(LFeature))
}
},
object : LupDecompositionFeature<Double> {
private val lup by lazy {
DecompositionFactory_DDRM.lu(origin.numRows(), origin.numCols()).apply { decompose(origin.ddrm.copy()) }
}
override val l: FeaturedMatrix<Double> by lazy {
EjmlMatrix(SimpleMatrix(lup.getLower(null)), setOf(LFeature))
}
override val u: FeaturedMatrix<Double> by lazy {
EjmlMatrix(SimpleMatrix(lup.getUpper(null)), setOf(UFeature))
}
override val p: FeaturedMatrix<Double> by lazy { EjmlMatrix(SimpleMatrix(lup.getRowPivot(null))) }
},
) union features
public override fun suggestFeature(vararg features: MatrixFeature): EjmlMatrix = public override fun suggestFeature(vararg features: MatrixFeature): EjmlMatrix =
EjmlMatrix(origin, this.features + features) EjmlMatrix(origin, this.features + features)

View File

@ -17,7 +17,6 @@ public fun Matrix<Double>.toEjml(): EjmlMatrix =
* @author Iaroslav Postovalov * @author Iaroslav Postovalov
*/ */
public object EjmlMatrixContext : MatrixContext<Double, EjmlMatrix> { public object EjmlMatrixContext : MatrixContext<Double, EjmlMatrix> {
/** /**
* Converts this vector to EJML one. * Converts this vector to EJML one.
*/ */
@ -33,6 +32,11 @@ public object EjmlMatrixContext : MatrixContext<Double, EjmlMatrix> {
} }
}) })
override fun point(size: Int, initializer: (Int) -> Double): Point<Double> =
EjmlVector(SimpleMatrix(size, 1).also {
(0 until it.numRows()).forEach { row -> it[row, 0] = initializer(row) }
})
public override fun Matrix<Double>.dot(other: Matrix<Double>): EjmlMatrix = public override fun Matrix<Double>.dot(other: Matrix<Double>): EjmlMatrix =
EjmlMatrix(toEjml().origin.mult(other.toEjml().origin)) EjmlMatrix(toEjml().origin.mult(other.toEjml().origin))
@ -73,12 +77,3 @@ public fun EjmlMatrixContext.solve(a: Matrix<Double>, b: Matrix<Double>): EjmlMa
*/ */
public fun EjmlMatrixContext.solve(a: Matrix<Double>, b: Point<Double>): EjmlVector = public fun EjmlMatrixContext.solve(a: Matrix<Double>, b: Point<Double>): EjmlVector =
EjmlVector(a.toEjml().origin.solve(b.toEjml().origin)) EjmlVector(a.toEjml().origin.solve(b.toEjml().origin))
/**
* Returns the inverse of given matrix: b = a^(-1).
*
* @param a the matrix.
* @return the inverse of this matrix.
* @author Iaroslav Postovalov
*/
public fun EjmlMatrixContext.inverse(a: Matrix<Double>): EjmlMatrix = EjmlMatrix(a.toEjml().origin.invert())

View File

@ -1,7 +1,7 @@
package kscience.kmath.ejml package kscience.kmath.ejml
import kscience.kmath.linear.DeterminantFeature import kscience.kmath.linear.DeterminantFeature
import kscience.kmath.linear.LUPDecompositionFeature import kscience.kmath.linear.LupDecompositionFeature
import kscience.kmath.linear.MatrixFeature import kscience.kmath.linear.MatrixFeature
import kscience.kmath.linear.getFeature import kscience.kmath.linear.getFeature
import org.ejml.dense.row.factory.DecompositionFactory_DDRM import org.ejml.dense.row.factory.DecompositionFactory_DDRM
@ -44,7 +44,7 @@ internal class EjmlMatrixTest {
val w = EjmlMatrix(m) val w = EjmlMatrix(m)
val det = w.getFeature<DeterminantFeature<Double>>() ?: fail() val det = w.getFeature<DeterminantFeature<Double>>() ?: fail()
assertEquals(m.determinant(), det.determinant) assertEquals(m.determinant(), det.determinant)
val lup = w.getFeature<LUPDecompositionFeature<Double>>() ?: fail() val lup = w.getFeature<LupDecompositionFeature<Double>>() ?: fail()
val ludecompositionF64 = DecompositionFactory_DDRM.lu(m.numRows(), m.numCols()) val ludecompositionF64 = DecompositionFactory_DDRM.lu(m.numRows(), m.numCols())
.also { it.decompose(m.ddrm.copy()) } .also { it.decompose(m.ddrm.copy()) }

View File

@ -31,6 +31,7 @@ extern "C"
TorchTensorHandle copy_to_long(TorchTensorHandle tensor_handle); TorchTensorHandle copy_to_long(TorchTensorHandle tensor_handle);
TorchTensorHandle copy_to_int(TorchTensorHandle tensor_handle); TorchTensorHandle copy_to_int(TorchTensorHandle tensor_handle);
void swap_tensors(TorchTensorHandle lhs_handle, TorchTensorHandle rhs_handle); void swap_tensors(TorchTensorHandle lhs_handle, TorchTensorHandle rhs_handle);
TorchTensorHandle view_tensor(TorchTensorHandle tensor_handle, int *shape, int dim);
char *tensor_to_string(TorchTensorHandle tensor_handle); char *tensor_to_string(TorchTensorHandle tensor_handle);
void dispose_char(char *ptr); void dispose_char(char *ptr);

View File

@ -74,6 +74,10 @@ void swap_tensors(TorchTensorHandle lhs_handle, TorchTensorHandle rhs_handle)
{ {
std::swap(ctorch::cast(lhs_handle), ctorch::cast(rhs_handle)); std::swap(ctorch::cast(lhs_handle), ctorch::cast(rhs_handle));
} }
TorchTensorHandle view_tensor(TorchTensorHandle tensor_handle, int *shape, int dim)
{
return new torch::Tensor(ctorch::cast(tensor_handle).view(ctorch::to_vec_int(shape, dim)));
}
char *tensor_to_string(TorchTensorHandle tensor_handle) char *tensor_to_string(TorchTensorHandle tensor_handle)
{ {

View File

@ -102,6 +102,9 @@ public sealed class TorchTensorAlgebra<
public infix fun TorchTensorType.swap(otherTensor: TorchTensorType): Unit { public infix fun TorchTensorType.swap(otherTensor: TorchTensorType): Unit {
swap_tensors(this.tensorHandle, otherTensor.tensorHandle) swap_tensors(this.tensorHandle, otherTensor.tensorHandle)
} }
public fun TorchTensorType.view(shape: IntArray): TorchTensorType =
wrap(view_tensor(this.tensorHandle, shape.toCValues(), shape.size)!!)
} }
public sealed class TorchTensorFieldAlgebra<T, TVar : CPrimitiveVar, public sealed class TorchTensorFieldAlgebra<T, TVar : CPrimitiveVar,

View File

@ -60,4 +60,13 @@ class TestTorchTensor {
} }
assertTrue(tensorInt.copyToFloat().copyToArray() contentEquals floatArrayOf(4f,5f,6f)) assertTrue(tensorInt.copyToFloat().copyToArray() contentEquals floatArrayOf(4f,5f,6f))
} }
@Test
fun testViewWithNoCopy() = TorchTensorIntAlgebra{
val tensor = copyFromArray(intArrayOf(1,2,3,4,5,6), shape = intArrayOf(6))
val viewTensor = tensor.view(intArrayOf(2,3))
assertTrue(viewTensor.shape contentEquals intArrayOf(2,3))
viewTensor[intArrayOf(0,0)] = 10
assertEquals(tensor[intArrayOf(0)], 10)
}
} }