From d9ebadd22ab0d92f8e98ed1bbec77f90ea9b4ce9 Mon Sep 17 00:00:00 2001 From: Iaroslav Postovalov Date: Sun, 24 Jan 2021 01:42:20 +0700 Subject: [PATCH] Implement LUP decomposition in GSL module --- .../kmath/gsl/codegen/matricesCodegen.kt | 41 ++- .../kmath/gsl/codegen/vectorsCodegen.kt | 20 +- .../kscience/kmath/linear/MatrixContext.kt | 48 ++- .../kscience/kmath/structures/Buffers.kt | 11 + .../src/nativeInterop/cinterop/libgsl.def | 2 +- .../kotlin/kscience/kmath/gsl/GslComplex.kt | 48 ++- .../kotlin/kscience/kmath/gsl/GslMatrix.kt | 25 +- .../kscience/kmath/gsl/GslMatrixContext.kt | 163 +++++++++-- .../kscience/kmath/gsl/GslMemoryHolder.kt | 20 +- .../kotlin/kscience/kmath/gsl/GslVector.kt | 4 +- .../kotlin/kscience/kmath/gsl/_Matrices.kt | 274 +++++++++++++----- .../kotlin/kscience/kmath/gsl/_Vectors.kt | 112 +++---- .../nativeTest/kotlin/GslMatrixRealTest.kt | 20 ++ 13 files changed, 568 insertions(+), 220 deletions(-) diff --git a/buildSrc/src/main/kotlin/kscience/kmath/gsl/codegen/matricesCodegen.kt b/buildSrc/src/main/kotlin/kscience/kmath/gsl/codegen/matricesCodegen.kt index 0dd880f7b..fd00a0255 100644 --- a/buildSrc/src/main/kotlin/kscience/kmath/gsl/codegen/matricesCodegen.kt +++ b/buildSrc/src/main/kotlin/kscience/kmath/gsl/codegen/matricesCodegen.kt @@ -14,38 +14,53 @@ private fun KtPsiFactory.createMatrixClass( kotlinTypeName: String, kotlinTypeAlias: String = kotlinTypeName ) { + fun fn(pattern: String) = fn(pattern, cTypeName) val className = "Gsl${kotlinTypeAlias}Matrix" val structName = sn("gsl_matrixR", cTypeName) @Language("kotlin") val text = """internal class $className( - override val nativeHandle: CPointer<$structName>, - scope: DeferScope + override val rawNativeHandle: CPointer<$structName>, + scope: AutofreeScope, ) : GslMatrix<$kotlinTypeName, $structName>(scope) { override val rowNum: Int - get() = nativeHandleChecked().pointed.size1.toInt() + get() = nativeHandle.pointed.size1.toInt() override val colNum: Int - get() = nativeHandleChecked().pointed.size2.toInt() + get() = nativeHandle.pointed.size2.toInt() + + override val rows: Buffer> + get() = VirtualBuffer(rowNum) { r -> + Gsl${kotlinTypeAlias}Vector( + ${fn("gsl_matrixRrow")}(nativeHandle, r.toULong()).placeTo(scope).pointed.vector.ptr, + scope, + ).apply { shouldBeFreed = false } + } + + override val columns: Buffer> + get() = VirtualBuffer(rowNum) { c -> + Gsl${kotlinTypeAlias}Vector( + ${fn("gsl_matrixRcolumn")}(nativeHandle, c.toULong()).placeTo(scope).pointed.vector.ptr, + scope, + ).apply { shouldBeFreed = false } + } override operator fun get(i: Int, j: Int): $kotlinTypeName = - ${fn("gsl_matrixRget", cTypeName)}(nativeHandleChecked(), i.toULong(), j.toULong()) + ${fn("gsl_matrixRget")}(nativeHandle, i.toULong(), j.toULong()) override operator fun set(i: Int, j: Int, value: ${kotlinTypeName}): Unit = - ${fn("gsl_matrixRset", cTypeName)}(nativeHandleChecked(), i.toULong(), j.toULong(), value) + ${fn("gsl_matrixRset")}(nativeHandle, i.toULong(), j.toULong(), value) override fun copy(): $className { - val new = requireNotNull(${fn("gsl_matrixRalloc", cTypeName)}(rowNum.toULong(), colNum.toULong())) - ${fn("gsl_matrixRmemcpy", cTypeName)}(new, nativeHandleChecked()) + val new = requireNotNull(${fn("gsl_matrixRalloc")}(rowNum.toULong(), colNum.toULong())) + ${fn("gsl_matrixRmemcpy")}(new, nativeHandle) return $className(new, scope) } - override fun close(): Unit = ${fn("gsl_matrixRfree", cTypeName)}(nativeHandleChecked()) + override fun close(): Unit = ${fn("gsl_matrixRfree")}(nativeHandle) override fun equals(other: Any?): Boolean { if (other is $className) - return ${ - fn("gsl_matrixRequal", cTypeName) - }(nativeHandleChecked(), other.nativeHandleChecked()) == 1 + return ${fn("gsl_matrixRequal")}(nativeHandle, other.nativeHandle) == 1 return super.equals(other) } @@ -64,7 +79,7 @@ fun matricesCodegen(outputFile: String, project: Project = createProject()) { f += createNewLine(2) f += createImportDirective(ImportPath.fromString("kotlinx.cinterop.*")) f += createNewLine(1) - f += createImportDirective(ImportPath.fromString("kscience.kmath.linear.*")) + f += createImportDirective(ImportPath.fromString("kscience.kmath.structures.*")) f += createNewLine(1) f += createImportDirective(ImportPath.fromString("org.gnu.gsl.*")) f += createNewLine(2) diff --git a/buildSrc/src/main/kotlin/kscience/kmath/gsl/codegen/vectorsCodegen.kt b/buildSrc/src/main/kotlin/kscience/kmath/gsl/codegen/vectorsCodegen.kt index 87acd258c..513d05727 100644 --- a/buildSrc/src/main/kotlin/kscience/kmath/gsl/codegen/vectorsCodegen.kt +++ b/buildSrc/src/main/kotlin/kscience/kmath/gsl/codegen/vectorsCodegen.kt @@ -14,38 +14,36 @@ private fun KtPsiFactory.createVectorClass( kotlinTypeName: String, kotlinTypeAlias: String = kotlinTypeName ) { + fun fn(pattern: String) = fn(pattern, cTypeName) val className = "Gsl${kotlinTypeAlias}Vector" val structName = sn("gsl_vectorR", cTypeName) @Language("kotlin") val text = - """internal class $className(override val nativeHandle: CPointer<$structName>, scope: DeferScope) : + """internal class $className(override val rawNativeHandle: CPointer<$structName>, scope: AutofreeScope) : GslVector<$kotlinTypeName, $structName>(scope) { override val size: Int - get() = nativeHandleChecked().pointed.size.toInt() + get() = nativeHandle.pointed.size.toInt() override operator fun get(index: Int): $kotlinTypeName = - ${fn("gsl_vectorRget", cTypeName)}(nativeHandleChecked(), index.toULong()) + ${fn("gsl_vectorRget")}(nativeHandle, index.toULong()) override operator fun set(index: Int, value: $kotlinTypeName): Unit = - ${fn("gsl_vectorRset", cTypeName)}(nativeHandleChecked(), index.toULong(), value) + ${fn("gsl_vectorRset")}(nativeHandle, index.toULong(), value) override fun copy(): $className { - val new = requireNotNull(${fn("gsl_vectorRalloc", cTypeName)}(size.toULong())) - ${fn("gsl_vectorRmemcpy", cTypeName)}(new, nativeHandleChecked()) + val new = requireNotNull(${fn("gsl_vectorRalloc")}(size.toULong())) + ${fn("gsl_vectorRmemcpy")}(new, nativeHandle) return ${className}(new, scope) } override fun equals(other: Any?): Boolean { if (other is $className) - return ${ - fn("gsl_vectorRequal", - cTypeName) - }(nativeHandleChecked(), other.nativeHandleChecked()) == 1 + return ${fn("gsl_vectorRequal")}(nativeHandle, other.nativeHandle) == 1 return super.equals(other) } - override fun close(): Unit = ${fn("gsl_vectorRfree", cTypeName)}(nativeHandleChecked()) + override fun close(): Unit = ${fn("gsl_vectorRfree")}(nativeHandle) }""" f += createClass(text) diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixContext.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixContext.kt index 40c7cd5a4..ad1674c43 100644 --- a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixContext.kt +++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixContext.kt @@ -1,5 +1,6 @@ package kscience.kmath.linear +import kscience.kmath.misc.UnstableKMathAPI import kscience.kmath.operations.Ring import kscience.kmath.operations.SpaceOperations import kscience.kmath.operations.invoke @@ -8,18 +9,22 @@ import kscience.kmath.structures.Buffer import kscience.kmath.structures.BufferFactory import kscience.kmath.structures.Matrix import kscience.kmath.structures.asSequence +import kotlin.reflect.KClass /** - * Basic operations on matrices. Operates on [Matrix] + * Basic operations on matrices. Operates on [Matrix]. + * + * @param T the type of items in the matrices. + * @param M the type of operated matrices. */ public interface MatrixContext> : SpaceOperations> { /** - * Produce a matrix with this context and given dimensions + * Produces a matrix with this context and given dimensions. */ public fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> T): M /** - * Produce a point compatible with matrix space (and possibly optimized for it) + * Produces a point compatible with matrix space (and possibly optimized for it). */ public fun point(size: Int, initializer: (Int) -> T): Point = Buffer.boxing(size, initializer) @@ -66,8 +71,19 @@ public interface MatrixContext> : SpaceOperations): M = m * this - public companion object { + /** + * Gets a feature from the matrix. This function may return some additional features than + * [kscience.kmath.structures.NDStructure.getFeature]. + * + * @param F the type of feature. + * @param m the matrix. + * @param type the [KClass] instance of [F]. + * @return a feature object or `null` if it isn't present. + */ + @UnstableKMathAPI + public fun getFeature(m: Matrix, type: KClass): F? = m.getFeature(type) + public companion object { /** * A structured matrix with custom buffer */ @@ -84,9 +100,31 @@ public interface MatrixContext> : SpaceOperations MatrixContext.getFeature(m: Matrix): F? = + getFeature(m, F::class) + +/** + * Partial implementation of [MatrixContext] for matrices of [Ring]. + * + * @param T the type of items in the matrices. + * @param R the type of ring of matrix elements. + * @param M the type of operated matrices. + */ public interface GenericMatrixContext, out M : Matrix> : MatrixContext { /** - * The ring context for matrix elements + * The ring instance of [T]. */ public val elementContext: R diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/Buffers.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/Buffers.kt index bfec6f871..ff19a4103 100644 --- a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/Buffers.kt +++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/Buffers.kt @@ -241,11 +241,22 @@ public class ArrayBuffer(private val array: Array) : MutableBuffer { override fun copy(): MutableBuffer = ArrayBuffer(array.copyOf()) } + /** * Returns an [ArrayBuffer] that wraps the original array. */ public fun Array.asBuffer(): ArrayBuffer = ArrayBuffer(this) +/** + * Creates a new [ArrayBuffer] with the specified [size], where each element is calculated by calling the specified + * [init] function. + * + * The function [init] is called for each array element sequentially starting from the first one. + * It should return the value for an array element given its index. + */ +public inline fun ArrayBuffer(size: Int, init: (Int) -> T): ArrayBuffer = + Array(size) { i -> init(i) }.asBuffer() + /** * Immutable wrapper for [MutableBuffer]. * diff --git a/kmath-gsl/src/nativeInterop/cinterop/libgsl.def b/kmath-gsl/src/nativeInterop/cinterop/libgsl.def index bd3a861f0..c50dd41cd 100644 --- a/kmath-gsl/src/nativeInterop/cinterop/libgsl.def +++ b/kmath-gsl/src/nativeInterop/cinterop/libgsl.def @@ -1,5 +1,5 @@ package=org.gnu.gsl -headers=gsl/gsl_blas.h gsl/gsl_linalg.h +headers=gsl/gsl_blas.h gsl/gsl_linalg.h gsl/gsl_permute_matrix.h compilerOpts.linux=-I/usr/include compilerOpts.osx=-I/usr/local -I/usr/local/include linkerOpts.linux=-L/usr/lib64 -L/usr/lib/x86_64-linux-gnu -lgsl diff --git a/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslComplex.kt b/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslComplex.kt index 0233ed8bc..6881bb0e4 100644 --- a/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslComplex.kt +++ b/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslComplex.kt @@ -2,6 +2,8 @@ package kscience.kmath.gsl import kotlinx.cinterop.* import kscience.kmath.operations.Complex +import kscience.kmath.structures.Buffer +import kscience.kmath.structures.VirtualBuffer import org.gnu.gsl.* internal fun CValue.toKMath(): Complex = useContents { Complex(dat[0], dat[1]) } @@ -11,58 +13,74 @@ internal fun Complex.toGsl(): CValue = cValue { dat[1] = im } -internal class GslComplexMatrix(override val nativeHandle: CPointer, scope: DeferScope) : +internal class GslComplexMatrix(override val rawNativeHandle: CPointer, scope: AutofreeScope) : GslMatrix(scope) { override val rowNum: Int - get() = nativeHandleChecked().pointed.size1.toInt() + get() = nativeHandle.pointed.size1.toInt() override val colNum: Int - get() = nativeHandleChecked().pointed.size2.toInt() + get() = nativeHandle.pointed.size2.toInt() + + override val rows: Buffer> + get() = VirtualBuffer(rowNum) { r -> + GslComplexVector( + gsl_matrix_complex_row(nativeHandle, r.toULong()).placeTo(scope).pointed.vector.ptr, + scope, + ).apply { shouldBeFreed = false } + } + + override val columns: Buffer> + get() = VirtualBuffer(rowNum) { c -> + GslComplexVector( + gsl_matrix_complex_column(nativeHandle, c.toULong()).placeTo(scope).pointed.vector.ptr, + scope, + ).apply { shouldBeFreed = false } + } override operator fun get(i: Int, j: Int): Complex = - gsl_matrix_complex_get(nativeHandleChecked(), i.toULong(), j.toULong()).toKMath() + gsl_matrix_complex_get(nativeHandle, i.toULong(), j.toULong()).toKMath() override operator fun set(i: Int, j: Int, value: Complex): Unit = - gsl_matrix_complex_set(nativeHandleChecked(), i.toULong(), j.toULong(), value.toGsl()) + gsl_matrix_complex_set(nativeHandle, i.toULong(), j.toULong(), value.toGsl()) override fun copy(): GslComplexMatrix { val new = requireNotNull(gsl_matrix_complex_alloc(rowNum.toULong(), colNum.toULong())) - gsl_matrix_complex_memcpy(new, nativeHandleChecked()) + gsl_matrix_complex_memcpy(new, nativeHandle) return GslComplexMatrix(new, scope) } - override fun close(): Unit = gsl_matrix_complex_free(nativeHandleChecked()) + override fun close(): Unit = gsl_matrix_complex_free(nativeHandle) override fun equals(other: Any?): Boolean { if (other is GslComplexMatrix) - return gsl_matrix_complex_equal(nativeHandleChecked(), other.nativeHandleChecked()) == 1 + return gsl_matrix_complex_equal(nativeHandle, other.nativeHandle) == 1 return super.equals(other) } } -internal class GslComplexVector(override val nativeHandle: CPointer, scope: DeferScope) : +internal class GslComplexVector(override val rawNativeHandle: CPointer, scope: AutofreeScope) : GslVector(scope) { override val size: Int - get() = nativeHandleChecked().pointed.size.toInt() + get() = nativeHandle.pointed.size.toInt() - override fun get(index: Int): Complex = gsl_vector_complex_get(nativeHandleChecked(), index.toULong()).toKMath() + override fun get(index: Int): Complex = gsl_vector_complex_get(nativeHandle, index.toULong()).toKMath() override fun set(index: Int, value: Complex): Unit = - gsl_vector_complex_set(nativeHandleChecked(), index.toULong(), value.toGsl()) + gsl_vector_complex_set(nativeHandle, index.toULong(), value.toGsl()) override fun copy(): GslComplexVector { val new = requireNotNull(gsl_vector_complex_alloc(size.toULong())) - gsl_vector_complex_memcpy(new, nativeHandleChecked()) + gsl_vector_complex_memcpy(new, nativeHandle) return GslComplexVector(new, scope) } override fun equals(other: Any?): Boolean { if (other is GslComplexVector) - return gsl_vector_complex_equal(nativeHandleChecked(), other.nativeHandleChecked()) == 1 + return gsl_vector_complex_equal(nativeHandle, other.nativeHandle) == 1 return super.equals(other) } - override fun close(): Unit = gsl_vector_complex_free(nativeHandleChecked()) + override fun close(): Unit = gsl_vector_complex_free(nativeHandle) } diff --git a/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslMatrix.kt b/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslMatrix.kt index 13283946e..c28236a67 100644 --- a/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslMatrix.kt +++ b/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslMatrix.kt @@ -1,33 +1,40 @@ package kscience.kmath.gsl +import kotlinx.cinterop.AutofreeScope import kotlinx.cinterop.CStructVar -import kotlinx.cinterop.DeferScope import kscience.kmath.structures.Matrix import kscience.kmath.structures.NDStructure +import kscience.kmath.structures.asSequence /** * Wraps gsl_matrix_* objects from GSL. */ -public abstract class GslMatrix internal constructor(scope: DeferScope) : +public abstract class GslMatrix internal constructor(scope: AutofreeScope) : GslMemoryHolder(scope), Matrix { internal abstract operator fun set(i: Int, j: Int, value: T) internal abstract fun copy(): GslMatrix public override fun equals(other: Any?): Boolean { - return NDStructure.equals(this, other as? NDStructure<*> ?: return false) + return NDStructure.contentEquals(this, other as? NDStructure<*> ?: return false) } public override fun hashCode(): Int { var ret = 7 - val nRows = rowNum - val nCols = colNum - ret = ret * 31 + nRows - ret = ret * 31 + nCols + ret = ret * 31 + rowNum + ret = ret * 31 + colNum - for (row in 0 until nRows) - for (col in 0 until nCols) + for (row in 0 until rowNum) + for (col in 0 until colNum) ret = ret * 31 + (11 * (row + 1) + 17 * (col + 1)) * this[row, col].hashCode() return ret } + + public override fun toString(): String = if (rowNum <= 5 && colNum <= 5) + "Matrix(rowsNum = $rowNum, colNum = $colNum)\n" + + rows.asSequence().joinToString(prefix = "(", postfix = ")", separator = "\n ") { buffer -> + buffer.asSequence().joinToString(separator = "\t") { it.toString() } + } + else + "Matrix(rowsNum = $rowNum, colNum = $colNum)" } diff --git a/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslMatrixContext.kt b/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslMatrixContext.kt index 95437c8c7..c32d852ba 100644 --- a/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslMatrixContext.kt +++ b/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslMatrixContext.kt @@ -1,31 +1,33 @@ package kscience.kmath.gsl -import kotlinx.cinterop.CStructVar -import kotlinx.cinterop.DeferScope -import kotlinx.cinterop.memScoped -import kotlinx.cinterop.pointed -import kscience.kmath.linear.MatrixContext -import kscience.kmath.linear.Point +import kotlinx.cinterop.* +import kscience.kmath.linear.* +import kscience.kmath.misc.UnstableKMathAPI import kscience.kmath.operations.Complex import kscience.kmath.operations.ComplexField import kscience.kmath.operations.toComplex import kscience.kmath.structures.Matrix +import kscience.kmath.structures.toList import org.gnu.gsl.* +import kotlin.reflect.KClass +import kotlin.reflect.cast internal inline fun GslMatrix.fill(initializer: (Int, Int) -> T): GslMatrix = apply { - (0 until rowNum).forEach { row -> (0 until colNum).forEach { col -> this[row, col] = initializer(row, col) } } + for (row in 0 until rowNum) { + for (col in 0 until colNum) this[row, col] = initializer(row, col) + } } internal inline fun GslVector.fill(initializer: (Int) -> T): GslVector = - apply { (0 until size).forEach { index -> this[index] = initializer(index) } } + apply { + for (index in 0 until size) this[index] = initializer(index) + } /** * Represents matrix context implementing where all the operations are delegated to GSL. */ -public abstract class GslMatrixContext internal constructor( - internal val scope: DeferScope, -) : MatrixContext> { +public abstract class GslMatrixContext : MatrixContext> { init { ensureHasGslErrorHandler() } @@ -52,19 +54,22 @@ public abstract class GslMatrixContext T): GslMatrix = produceDirtyMatrix(rows, columns).fill(initializer) - public override fun point(size: Int, initializer: (Int) -> T): Point = + public override fun point(size: Int, initializer: (Int) -> T): GslVector = produceDirtyVector(size).fill(initializer) } /** * Represents [Double] matrix context implementing where all the operations are delegated to GSL. */ -public class GslRealMatrixContext(scope: DeferScope) : GslMatrixContext(scope) { - override fun produceDirtyMatrix(rows: Int, columns: Int): GslMatrix = - GslRealMatrix(nativeHandle = requireNotNull(gsl_matrix_alloc(rows.toULong(), columns.toULong())), scope = scope) +public class GslRealMatrixContext(internal val scope: AutofreeScope) : + GslMatrixContext() { + override fun produceDirtyMatrix(rows: Int, columns: Int): GslMatrix = GslRealMatrix( + rawNativeHandle = requireNotNull(gsl_matrix_alloc(rows.toULong(), columns.toULong())), + scope = scope, + ) override fun produceDirtyVector(size: Int): GslVector = - GslRealVector(nativeHandle = requireNotNull(gsl_vector_alloc(size.toULong())), scope = scope) + GslRealVector(rawNativeHandle = requireNotNull(gsl_vector_alloc(size.toULong())), scope = scope) public override fun Matrix.dot(other: Matrix): GslMatrix { val x = toGsl().nativeHandle @@ -105,6 +110,58 @@ public class GslRealMatrixContext(scope: DeferScope) : GslMatrixContext getFeature(m: Matrix, type: KClass): F? = when (type) { + LupDecompositionFeature::class -> object : LupDecompositionFeature { + private val lup by lazy { + val lu = m.toGsl().copy() + val n = lu.rowNum + val perm = gsl_permutation_alloc(n.toULong()) + + memScoped { + val i = alloc() + gsl_linalg_LU_decomp(lu.nativeHandle, perm, i.ptr) + } + + val one = produce(n, n) { i, j -> if (i == j) 1.0 else 0.0 } + val p = produce(n, n) { _, _ -> 0.0 } + + for (j in 0 until perm!!.pointed.size.toInt()) { + println("7 $j") + val k = gsl_permutation_get(perm, j.toULong()).toInt() + val col = one.columns[k] + println(p.toGsl()) + println(col.toList()) + gsl_matrix_set_col(p.nativeHandle, j.toULong(), col.toGsl().nativeHandle) + } + + gsl_permutation_free(perm) + lu to p + } + + override val l by lazy { + VirtualMatrix(lup.first.shape[0], lup.first.shape[1]) { i, j -> + when { + j < i -> lup.first[i, j] + j == i -> 1.0 + else -> 0.0 + } + } + LFeature + } + + override val u by lazy { + VirtualMatrix( + lup.first.shape[0], + lup.first.shape[1], + ) { i, j -> if (j >= i) lup.first[i, j] else 0.0 } + UFeature + } + + override val p by lazy(lup::second) + } + else -> super.getFeature(m, type) + }?.let(type::cast) } /** @@ -116,22 +173,22 @@ public fun GslRealMatrixContext(block: GslRealMatrixContext.() -> R): R = /** * Represents [Float] matrix context implementing where all the operations are delegated to GSL. */ -public class GslFloatMatrixContext(scope: DeferScope) : - GslMatrixContext(scope) { +public class GslFloatMatrixContext(internal val scope: AutofreeScope) : + GslMatrixContext() { override fun produceDirtyMatrix(rows: Int, columns: Int): GslMatrix = GslFloatMatrix( - nativeHandle = requireNotNull(gsl_matrix_float_alloc(rows.toULong(), columns.toULong())), + rawNativeHandle = requireNotNull(gsl_matrix_float_alloc(rows.toULong(), columns.toULong())), scope = scope, ) override fun produceDirtyVector(size: Int): GslVector = - GslFloatVector(nativeHandle = requireNotNull(value = gsl_vector_float_alloc(size.toULong())), scope = scope) + GslFloatVector(rawNativeHandle = requireNotNull(value = gsl_vector_float_alloc(size.toULong())), scope = scope) public override fun Matrix.dot(other: Matrix): GslMatrix { val x = toGsl().nativeHandle val a = other.toGsl().nativeHandle val result = requireNotNull(gsl_matrix_float_calloc(a.pointed.size1, a.pointed.size2)) gsl_blas_sgemm(CblasNoTrans, CblasNoTrans, 1f, x, a, 1f, result) - return GslFloatMatrix(nativeHandle = result, scope = scope) + return GslFloatMatrix(rawNativeHandle = result, scope = scope) } public override fun Matrix.dot(vector: Point): GslVector { @@ -139,7 +196,7 @@ public class GslFloatMatrixContext(scope: DeferScope) : val a = vector.toGsl().nativeHandle val result = requireNotNull(gsl_vector_float_calloc(a.pointed.size)) gsl_blas_sgemv(CblasNoTrans, 1f, x, a, 1f, result) - return GslFloatVector(nativeHandle = result, scope = scope) + return GslFloatVector(rawNativeHandle = result, scope = scope) } public override fun Matrix.times(value: Float): GslMatrix { @@ -176,22 +233,22 @@ public fun GslFloatMatrixContext(block: GslFloatMatrixContext.() -> R): R = /** * Represents [Complex] matrix context implementing where all the operations are delegated to GSL. */ -public class GslComplexMatrixContext(scope: DeferScope) : - GslMatrixContext(scope) { +public class GslComplexMatrixContext(internal val scope: AutofreeScope) : + GslMatrixContext() { override fun produceDirtyMatrix(rows: Int, columns: Int): GslMatrix = GslComplexMatrix( - nativeHandle = requireNotNull(gsl_matrix_complex_alloc(rows.toULong(), columns.toULong())), + rawNativeHandle = requireNotNull(gsl_matrix_complex_alloc(rows.toULong(), columns.toULong())), scope = scope, ) override fun produceDirtyVector(size: Int): GslVector = - GslComplexVector(nativeHandle = requireNotNull(gsl_vector_complex_alloc(size.toULong())), scope = scope) + GslComplexVector(rawNativeHandle = requireNotNull(gsl_vector_complex_alloc(size.toULong())), scope = scope) public override fun Matrix.dot(other: Matrix): GslMatrix { val x = toGsl().nativeHandle val a = other.toGsl().nativeHandle val result = requireNotNull(gsl_matrix_complex_calloc(a.pointed.size1, a.pointed.size2)) gsl_blas_zgemm(CblasNoTrans, CblasNoTrans, ComplexField.one.toGsl(), x, a, ComplexField.one.toGsl(), result) - return GslComplexMatrix(nativeHandle = result, scope = scope) + return GslComplexMatrix(rawNativeHandle = result, scope = scope) } public override fun Matrix.dot(vector: Point): GslVector { @@ -225,6 +282,58 @@ public class GslComplexMatrixContext(scope: DeferScope) : gsl_matrix_complex_sub(g1.nativeHandle, b.toGsl().nativeHandle) return g1 } + + @Suppress("IMPLICIT_CAST_TO_ANY") + @UnstableKMathAPI + public override fun getFeature(m: Matrix, type: KClass): F? = when (type) { + LupDecompositionFeature::class -> object : LupDecompositionFeature { + private val lup by lazy { + val lu = m.toGsl().copy() + val n = lu.rowNum + val perm = gsl_permutation_alloc(n.toULong()) + + memScoped { + val i = alloc() + gsl_linalg_complex_LU_decomp(lu.nativeHandle, perm, i.ptr) + } + + val one = produce(n, n) { i, j -> if (i == j) 1.0.toComplex() else 0.0.toComplex() } + val p = produce(n, n) { _, _ -> 0.0.toComplex() } + + for (j in 0 until perm!!.pointed.size.toInt()) { + println("7 $j") + val k = gsl_permutation_get(perm, j.toULong()).toInt() + val col = one.columns[k] + println(p.toGsl()) + println(col.toList()) + gsl_matrix_complex_set_col(p.nativeHandle, j.toULong(), col.toGsl().nativeHandle) + } + + gsl_permutation_free(perm) + lu to p + } + + override val l by lazy { + VirtualMatrix(lup.first.shape[0], lup.first.shape[1]) { i, j -> + when { + j < i -> lup.first[i, j] + j == i -> 1.0.toComplex() + else -> 0.0.toComplex() + } + } + LFeature + } + + override val u by lazy { + VirtualMatrix( + lup.first.shape[0], + lup.first.shape[1], + ) { i, j -> if (j >= i) lup.first[i, j] else 0.0.toComplex() } + UFeature + } + + override val p by lazy(lup::second) + } + else -> super.getFeature(m, type) + }?.let(type::cast) } /** diff --git a/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslMemoryHolder.kt b/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslMemoryHolder.kt index d68a52a62..46d1ba60b 100644 --- a/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslMemoryHolder.kt +++ b/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslMemoryHolder.kt @@ -1,5 +1,6 @@ package kscience.kmath.gsl +import kotlinx.cinterop.AutofreeScope import kotlinx.cinterop.CPointer import kotlinx.cinterop.CStructVar import kotlinx.cinterop.DeferScope @@ -12,23 +13,26 @@ import kotlinx.cinterop.DeferScope * * @param scope the scope where this object is declared. */ -public abstract class GslMemoryHolder internal constructor(internal val scope: DeferScope) { - internal abstract val nativeHandle: CPointer +public abstract class GslMemoryHolder internal constructor(internal val scope: AutofreeScope) { + internal abstract val rawNativeHandle: CPointer private var isClosed: Boolean = false + internal val nativeHandle: CPointer + get() { + check(!isClosed) { "The use of GSL object that is closed." } + return rawNativeHandle + } + + internal var shouldBeFreed = true + init { ensureHasGslErrorHandler() scope.defer { - close() + if (shouldBeFreed) close() isClosed = true } } - internal fun nativeHandleChecked(): CPointer { - check(!isClosed) { "The use of GSL object that is closed." } - return nativeHandle - } - internal abstract fun close() } diff --git a/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslVector.kt b/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslVector.kt index 0033c47a1..8bc1240b4 100644 --- a/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslVector.kt +++ b/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslVector.kt @@ -1,13 +1,13 @@ package kscience.kmath.gsl +import kotlinx.cinterop.AutofreeScope import kotlinx.cinterop.CStructVar -import kotlinx.cinterop.DeferScope import kscience.kmath.linear.Point /** * Wraps gsl_vector_* objects from GSL. */ -public abstract class GslVector internal constructor(scope: DeferScope) : +public abstract class GslVector internal constructor(scope: AutofreeScope) : GslMemoryHolder(scope), Point { internal abstract operator fun set(index: Int, value: T) internal abstract fun copy(): GslVector diff --git a/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/_Matrices.kt b/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/_Matrices.kt index 07bc8af79..1a826d493 100644 --- a/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/_Matrices.kt +++ b/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/_Matrices.kt @@ -1,260 +1,388 @@ package kscience.kmath.gsl import kotlinx.cinterop.* -import kscience.kmath.linear.* +import kscience.kmath.structures.* import org.gnu.gsl.* internal class GslRealMatrix( - override val nativeHandle: CPointer, - scope: DeferScope + override val rawNativeHandle: CPointer, + scope: AutofreeScope, ) : GslMatrix(scope) { override val rowNum: Int - get() = nativeHandleChecked().pointed.size1.toInt() + get() = nativeHandle.pointed.size1.toInt() override val colNum: Int - get() = nativeHandleChecked().pointed.size2.toInt() + get() = nativeHandle.pointed.size2.toInt() + + override val rows: Buffer> + get() = VirtualBuffer(rowNum) { r -> + GslRealVector( + gsl_matrix_row(nativeHandle, r.toULong()).placeTo(scope).pointed.vector.ptr, + scope, + ).apply { shouldBeFreed = false } + } + + override val columns: Buffer> + get() = VirtualBuffer(rowNum) { c -> + GslRealVector( + gsl_matrix_column(nativeHandle, c.toULong()).placeTo(scope).pointed.vector.ptr, + scope, + ).apply { shouldBeFreed = false } + } override operator fun get(i: Int, j: Int): Double = - gsl_matrix_get(nativeHandleChecked(), i.toULong(), j.toULong()) + gsl_matrix_get(nativeHandle, i.toULong(), j.toULong()) override operator fun set(i: Int, j: Int, value: Double): Unit = - gsl_matrix_set(nativeHandleChecked(), i.toULong(), j.toULong(), value) + gsl_matrix_set(nativeHandle, i.toULong(), j.toULong(), value) override fun copy(): GslRealMatrix { val new = requireNotNull(gsl_matrix_alloc(rowNum.toULong(), colNum.toULong())) - gsl_matrix_memcpy(new, nativeHandleChecked()) + gsl_matrix_memcpy(new, nativeHandle) return GslRealMatrix(new, scope) } - override fun close(): Unit = gsl_matrix_free(nativeHandleChecked()) + override fun close(): Unit = gsl_matrix_free(nativeHandle) override fun equals(other: Any?): Boolean { if (other is GslRealMatrix) - return gsl_matrix_equal(nativeHandleChecked(), other.nativeHandleChecked()) == 1 + return gsl_matrix_equal(nativeHandle, other.nativeHandle) == 1 return super.equals(other) } } internal class GslFloatMatrix( - override val nativeHandle: CPointer, - scope: DeferScope + override val rawNativeHandle: CPointer, + scope: AutofreeScope, ) : GslMatrix(scope) { override val rowNum: Int - get() = nativeHandleChecked().pointed.size1.toInt() + get() = nativeHandle.pointed.size1.toInt() override val colNum: Int - get() = nativeHandleChecked().pointed.size2.toInt() + get() = nativeHandle.pointed.size2.toInt() + + override val rows: Buffer> + get() = VirtualBuffer(rowNum) { r -> + GslFloatVector( + gsl_matrix_float_row(nativeHandle, r.toULong()).placeTo(scope).pointed.vector.ptr, + scope, + ).apply { shouldBeFreed = false } + } + + override val columns: Buffer> + get() = VirtualBuffer(rowNum) { c -> + GslFloatVector( + gsl_matrix_float_column(nativeHandle, c.toULong()).placeTo(scope).pointed.vector.ptr, + scope, + ).apply { shouldBeFreed = false } + } override operator fun get(i: Int, j: Int): Float = - gsl_matrix_float_get(nativeHandleChecked(), i.toULong(), j.toULong()) + gsl_matrix_float_get(nativeHandle, i.toULong(), j.toULong()) override operator fun set(i: Int, j: Int, value: Float): Unit = - gsl_matrix_float_set(nativeHandleChecked(), i.toULong(), j.toULong(), value) + gsl_matrix_float_set(nativeHandle, i.toULong(), j.toULong(), value) override fun copy(): GslFloatMatrix { val new = requireNotNull(gsl_matrix_float_alloc(rowNum.toULong(), colNum.toULong())) - gsl_matrix_float_memcpy(new, nativeHandleChecked()) + gsl_matrix_float_memcpy(new, nativeHandle) return GslFloatMatrix(new, scope) } - override fun close(): Unit = gsl_matrix_float_free(nativeHandleChecked()) + override fun close(): Unit = gsl_matrix_float_free(nativeHandle) override fun equals(other: Any?): Boolean { if (other is GslFloatMatrix) - return gsl_matrix_float_equal(nativeHandleChecked(), other.nativeHandleChecked()) == 1 + return gsl_matrix_float_equal(nativeHandle, other.nativeHandle) == 1 return super.equals(other) } } internal class GslShortMatrix( - override val nativeHandle: CPointer, - scope: DeferScope + override val rawNativeHandle: CPointer, + scope: AutofreeScope, ) : GslMatrix(scope) { override val rowNum: Int - get() = nativeHandleChecked().pointed.size1.toInt() + get() = nativeHandle.pointed.size1.toInt() override val colNum: Int - get() = nativeHandleChecked().pointed.size2.toInt() + get() = nativeHandle.pointed.size2.toInt() + + override val rows: Buffer> + get() = VirtualBuffer(rowNum) { r -> + GslShortVector( + gsl_matrix_short_row(nativeHandle, r.toULong()).placeTo(scope).pointed.vector.ptr, + scope, + ).apply { shouldBeFreed = false } + } + + override val columns: Buffer> + get() = VirtualBuffer(rowNum) { c -> + GslShortVector( + gsl_matrix_short_column(nativeHandle, c.toULong()).placeTo(scope).pointed.vector.ptr, + scope, + ).apply { shouldBeFreed = false } + } override operator fun get(i: Int, j: Int): Short = - gsl_matrix_short_get(nativeHandleChecked(), i.toULong(), j.toULong()) + gsl_matrix_short_get(nativeHandle, i.toULong(), j.toULong()) override operator fun set(i: Int, j: Int, value: Short): Unit = - gsl_matrix_short_set(nativeHandleChecked(), i.toULong(), j.toULong(), value) + gsl_matrix_short_set(nativeHandle, i.toULong(), j.toULong(), value) override fun copy(): GslShortMatrix { val new = requireNotNull(gsl_matrix_short_alloc(rowNum.toULong(), colNum.toULong())) - gsl_matrix_short_memcpy(new, nativeHandleChecked()) + gsl_matrix_short_memcpy(new, nativeHandle) return GslShortMatrix(new, scope) } - override fun close(): Unit = gsl_matrix_short_free(nativeHandleChecked()) + override fun close(): Unit = gsl_matrix_short_free(nativeHandle) override fun equals(other: Any?): Boolean { if (other is GslShortMatrix) - return gsl_matrix_short_equal(nativeHandleChecked(), other.nativeHandleChecked()) == 1 + return gsl_matrix_short_equal(nativeHandle, other.nativeHandle) == 1 return super.equals(other) } } internal class GslUShortMatrix( - override val nativeHandle: CPointer, - scope: DeferScope + override val rawNativeHandle: CPointer, + scope: AutofreeScope, ) : GslMatrix(scope) { override val rowNum: Int - get() = nativeHandleChecked().pointed.size1.toInt() + get() = nativeHandle.pointed.size1.toInt() override val colNum: Int - get() = nativeHandleChecked().pointed.size2.toInt() + get() = nativeHandle.pointed.size2.toInt() + + override val rows: Buffer> + get() = VirtualBuffer(rowNum) { r -> + GslUShortVector( + gsl_matrix_ushort_row(nativeHandle, r.toULong()).placeTo(scope).pointed.vector.ptr, + scope, + ).apply { shouldBeFreed = false } + } + + override val columns: Buffer> + get() = VirtualBuffer(rowNum) { c -> + GslUShortVector( + gsl_matrix_ushort_column(nativeHandle, c.toULong()).placeTo(scope).pointed.vector.ptr, + scope, + ).apply { shouldBeFreed = false } + } override operator fun get(i: Int, j: Int): UShort = - gsl_matrix_ushort_get(nativeHandleChecked(), i.toULong(), j.toULong()) + gsl_matrix_ushort_get(nativeHandle, i.toULong(), j.toULong()) override operator fun set(i: Int, j: Int, value: UShort): Unit = - gsl_matrix_ushort_set(nativeHandleChecked(), i.toULong(), j.toULong(), value) + gsl_matrix_ushort_set(nativeHandle, i.toULong(), j.toULong(), value) override fun copy(): GslUShortMatrix { val new = requireNotNull(gsl_matrix_ushort_alloc(rowNum.toULong(), colNum.toULong())) - gsl_matrix_ushort_memcpy(new, nativeHandleChecked()) + gsl_matrix_ushort_memcpy(new, nativeHandle) return GslUShortMatrix(new, scope) } - override fun close(): Unit = gsl_matrix_ushort_free(nativeHandleChecked()) + override fun close(): Unit = gsl_matrix_ushort_free(nativeHandle) override fun equals(other: Any?): Boolean { if (other is GslUShortMatrix) - return gsl_matrix_ushort_equal(nativeHandleChecked(), other.nativeHandleChecked()) == 1 + return gsl_matrix_ushort_equal(nativeHandle, other.nativeHandle) == 1 return super.equals(other) } } internal class GslLongMatrix( - override val nativeHandle: CPointer, - scope: DeferScope + override val rawNativeHandle: CPointer, + scope: AutofreeScope, ) : GslMatrix(scope) { override val rowNum: Int - get() = nativeHandleChecked().pointed.size1.toInt() + get() = nativeHandle.pointed.size1.toInt() override val colNum: Int - get() = nativeHandleChecked().pointed.size2.toInt() + get() = nativeHandle.pointed.size2.toInt() + + override val rows: Buffer> + get() = VirtualBuffer(rowNum) { r -> + GslLongVector( + gsl_matrix_long_row(nativeHandle, r.toULong()).placeTo(scope).pointed.vector.ptr, + scope, + ).apply { shouldBeFreed = false } + } + + override val columns: Buffer> + get() = VirtualBuffer(rowNum) { c -> + GslLongVector( + gsl_matrix_long_column(nativeHandle, c.toULong()).placeTo(scope).pointed.vector.ptr, + scope, + ).apply { shouldBeFreed = false } + } override operator fun get(i: Int, j: Int): Long = - gsl_matrix_long_get(nativeHandleChecked(), i.toULong(), j.toULong()) + gsl_matrix_long_get(nativeHandle, i.toULong(), j.toULong()) override operator fun set(i: Int, j: Int, value: Long): Unit = - gsl_matrix_long_set(nativeHandleChecked(), i.toULong(), j.toULong(), value) + gsl_matrix_long_set(nativeHandle, i.toULong(), j.toULong(), value) override fun copy(): GslLongMatrix { val new = requireNotNull(gsl_matrix_long_alloc(rowNum.toULong(), colNum.toULong())) - gsl_matrix_long_memcpy(new, nativeHandleChecked()) + gsl_matrix_long_memcpy(new, nativeHandle) return GslLongMatrix(new, scope) } - override fun close(): Unit = gsl_matrix_long_free(nativeHandleChecked()) + override fun close(): Unit = gsl_matrix_long_free(nativeHandle) override fun equals(other: Any?): Boolean { if (other is GslLongMatrix) - return gsl_matrix_long_equal(nativeHandleChecked(), other.nativeHandleChecked()) == 1 + return gsl_matrix_long_equal(nativeHandle, other.nativeHandle) == 1 return super.equals(other) } } internal class GslULongMatrix( - override val nativeHandle: CPointer, - scope: DeferScope + override val rawNativeHandle: CPointer, + scope: AutofreeScope, ) : GslMatrix(scope) { override val rowNum: Int - get() = nativeHandleChecked().pointed.size1.toInt() + get() = nativeHandle.pointed.size1.toInt() override val colNum: Int - get() = nativeHandleChecked().pointed.size2.toInt() + get() = nativeHandle.pointed.size2.toInt() + + override val rows: Buffer> + get() = VirtualBuffer(rowNum) { r -> + GslULongVector( + gsl_matrix_ulong_row(nativeHandle, r.toULong()).placeTo(scope).pointed.vector.ptr, + scope, + ).apply { shouldBeFreed = false } + } + + override val columns: Buffer> + get() = VirtualBuffer(rowNum) { c -> + GslULongVector( + gsl_matrix_ulong_column(nativeHandle, c.toULong()).placeTo(scope).pointed.vector.ptr, + scope, + ).apply { shouldBeFreed = false } + } override operator fun get(i: Int, j: Int): ULong = - gsl_matrix_ulong_get(nativeHandleChecked(), i.toULong(), j.toULong()) + gsl_matrix_ulong_get(nativeHandle, i.toULong(), j.toULong()) override operator fun set(i: Int, j: Int, value: ULong): Unit = - gsl_matrix_ulong_set(nativeHandleChecked(), i.toULong(), j.toULong(), value) + gsl_matrix_ulong_set(nativeHandle, i.toULong(), j.toULong(), value) override fun copy(): GslULongMatrix { val new = requireNotNull(gsl_matrix_ulong_alloc(rowNum.toULong(), colNum.toULong())) - gsl_matrix_ulong_memcpy(new, nativeHandleChecked()) + gsl_matrix_ulong_memcpy(new, nativeHandle) return GslULongMatrix(new, scope) } - override fun close(): Unit = gsl_matrix_ulong_free(nativeHandleChecked()) + override fun close(): Unit = gsl_matrix_ulong_free(nativeHandle) override fun equals(other: Any?): Boolean { if (other is GslULongMatrix) - return gsl_matrix_ulong_equal(nativeHandleChecked(), other.nativeHandleChecked()) == 1 + return gsl_matrix_ulong_equal(nativeHandle, other.nativeHandle) == 1 return super.equals(other) } } internal class GslIntMatrix( - override val nativeHandle: CPointer, - scope: DeferScope + override val rawNativeHandle: CPointer, + scope: AutofreeScope, ) : GslMatrix(scope) { override val rowNum: Int - get() = nativeHandleChecked().pointed.size1.toInt() + get() = nativeHandle.pointed.size1.toInt() override val colNum: Int - get() = nativeHandleChecked().pointed.size2.toInt() + get() = nativeHandle.pointed.size2.toInt() + + override val rows: Buffer> + get() = VirtualBuffer(rowNum) { r -> + GslIntVector( + gsl_matrix_int_row(nativeHandle, r.toULong()).placeTo(scope).pointed.vector.ptr, + scope, + ).apply { shouldBeFreed = false } + } + + override val columns: Buffer> + get() = VirtualBuffer(rowNum) { c -> + GslIntVector( + gsl_matrix_int_column(nativeHandle, c.toULong()).placeTo(scope).pointed.vector.ptr, + scope, + ).apply { shouldBeFreed = false } + } override operator fun get(i: Int, j: Int): Int = - gsl_matrix_int_get(nativeHandleChecked(), i.toULong(), j.toULong()) + gsl_matrix_int_get(nativeHandle, i.toULong(), j.toULong()) override operator fun set(i: Int, j: Int, value: Int): Unit = - gsl_matrix_int_set(nativeHandleChecked(), i.toULong(), j.toULong(), value) + gsl_matrix_int_set(nativeHandle, i.toULong(), j.toULong(), value) override fun copy(): GslIntMatrix { val new = requireNotNull(gsl_matrix_int_alloc(rowNum.toULong(), colNum.toULong())) - gsl_matrix_int_memcpy(new, nativeHandleChecked()) + gsl_matrix_int_memcpy(new, nativeHandle) return GslIntMatrix(new, scope) } - override fun close(): Unit = gsl_matrix_int_free(nativeHandleChecked()) + override fun close(): Unit = gsl_matrix_int_free(nativeHandle) override fun equals(other: Any?): Boolean { if (other is GslIntMatrix) - return gsl_matrix_int_equal(nativeHandleChecked(), other.nativeHandleChecked()) == 1 + return gsl_matrix_int_equal(nativeHandle, other.nativeHandle) == 1 return super.equals(other) } } internal class GslUIntMatrix( - override val nativeHandle: CPointer, - scope: DeferScope + override val rawNativeHandle: CPointer, + scope: AutofreeScope, ) : GslMatrix(scope) { override val rowNum: Int - get() = nativeHandleChecked().pointed.size1.toInt() + get() = nativeHandle.pointed.size1.toInt() override val colNum: Int - get() = nativeHandleChecked().pointed.size2.toInt() + get() = nativeHandle.pointed.size2.toInt() + + override val rows: Buffer> + get() = VirtualBuffer(rowNum) { r -> + GslUIntVector( + gsl_matrix_uint_row(nativeHandle, r.toULong()).placeTo(scope).pointed.vector.ptr, + scope, + ).apply { shouldBeFreed = false } + } + + override val columns: Buffer> + get() = VirtualBuffer(rowNum) { c -> + GslUIntVector( + gsl_matrix_uint_column(nativeHandle, c.toULong()).placeTo(scope).pointed.vector.ptr, + scope, + ).apply { shouldBeFreed = false } + } override operator fun get(i: Int, j: Int): UInt = - gsl_matrix_uint_get(nativeHandleChecked(), i.toULong(), j.toULong()) + gsl_matrix_uint_get(nativeHandle, i.toULong(), j.toULong()) override operator fun set(i: Int, j: Int, value: UInt): Unit = - gsl_matrix_uint_set(nativeHandleChecked(), i.toULong(), j.toULong(), value) + gsl_matrix_uint_set(nativeHandle, i.toULong(), j.toULong(), value) override fun copy(): GslUIntMatrix { val new = requireNotNull(gsl_matrix_uint_alloc(rowNum.toULong(), colNum.toULong())) - gsl_matrix_uint_memcpy(new, nativeHandleChecked()) + gsl_matrix_uint_memcpy(new, nativeHandle) return GslUIntMatrix(new, scope) } - override fun close(): Unit = gsl_matrix_uint_free(nativeHandleChecked()) + override fun close(): Unit = gsl_matrix_uint_free(nativeHandle) override fun equals(other: Any?): Boolean { if (other is GslUIntMatrix) - return gsl_matrix_uint_equal(nativeHandleChecked(), other.nativeHandleChecked()) == 1 + return gsl_matrix_uint_equal(nativeHandle, other.nativeHandle) == 1 return super.equals(other) } diff --git a/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/_Vectors.kt b/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/_Vectors.kt index d7ed4718b..5b5236242 100644 --- a/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/_Vectors.kt +++ b/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/_Vectors.kt @@ -3,219 +3,219 @@ package kscience.kmath.gsl import kotlinx.cinterop.* import org.gnu.gsl.* -internal class GslRealVector(override val nativeHandle: CPointer, scope: DeferScope) : +internal class GslRealVector(override val rawNativeHandle: CPointer, scope: AutofreeScope) : GslVector(scope) { override val size: Int - get() = nativeHandleChecked().pointed.size.toInt() + get() = nativeHandle.pointed.size.toInt() override operator fun get(index: Int): Double = - gsl_vector_get(nativeHandleChecked(), index.toULong()) + gsl_vector_get(nativeHandle, index.toULong()) override operator fun set(index: Int, value: Double): Unit = - gsl_vector_set(nativeHandleChecked(), index.toULong(), value) + gsl_vector_set(nativeHandle, index.toULong(), value) override fun copy(): GslRealVector { val new = requireNotNull(gsl_vector_alloc(size.toULong())) - gsl_vector_memcpy(new, nativeHandleChecked()) + gsl_vector_memcpy(new, nativeHandle) return GslRealVector(new, scope) } override fun equals(other: Any?): Boolean { if (other is GslRealVector) - return gsl_vector_equal(nativeHandleChecked(), other.nativeHandleChecked()) == 1 + return gsl_vector_equal(nativeHandle, other.nativeHandle) == 1 return super.equals(other) } - override fun close(): Unit = gsl_vector_free(nativeHandleChecked()) + override fun close(): Unit = gsl_vector_free(nativeHandle) } -internal class GslFloatVector(override val nativeHandle: CPointer, scope: DeferScope) : +internal class GslFloatVector(override val rawNativeHandle: CPointer, scope: AutofreeScope) : GslVector(scope) { override val size: Int - get() = nativeHandleChecked().pointed.size.toInt() + get() = nativeHandle.pointed.size.toInt() override operator fun get(index: Int): Float = - gsl_vector_float_get(nativeHandleChecked(), index.toULong()) + gsl_vector_float_get(nativeHandle, index.toULong()) override operator fun set(index: Int, value: Float): Unit = - gsl_vector_float_set(nativeHandleChecked(), index.toULong(), value) + gsl_vector_float_set(nativeHandle, index.toULong(), value) override fun copy(): GslFloatVector { val new = requireNotNull(gsl_vector_float_alloc(size.toULong())) - gsl_vector_float_memcpy(new, nativeHandleChecked()) + gsl_vector_float_memcpy(new, nativeHandle) return GslFloatVector(new, scope) } override fun equals(other: Any?): Boolean { if (other is GslFloatVector) - return gsl_vector_float_equal(nativeHandleChecked(), other.nativeHandleChecked()) == 1 + return gsl_vector_float_equal(nativeHandle, other.nativeHandle) == 1 return super.equals(other) } - override fun close(): Unit = gsl_vector_float_free(nativeHandleChecked()) + override fun close(): Unit = gsl_vector_float_free(nativeHandle) } -internal class GslShortVector(override val nativeHandle: CPointer, scope: DeferScope) : +internal class GslShortVector(override val rawNativeHandle: CPointer, scope: AutofreeScope) : GslVector(scope) { override val size: Int - get() = nativeHandleChecked().pointed.size.toInt() + get() = nativeHandle.pointed.size.toInt() override operator fun get(index: Int): Short = - gsl_vector_short_get(nativeHandleChecked(), index.toULong()) + gsl_vector_short_get(nativeHandle, index.toULong()) override operator fun set(index: Int, value: Short): Unit = - gsl_vector_short_set(nativeHandleChecked(), index.toULong(), value) + gsl_vector_short_set(nativeHandle, index.toULong(), value) override fun copy(): GslShortVector { val new = requireNotNull(gsl_vector_short_alloc(size.toULong())) - gsl_vector_short_memcpy(new, nativeHandleChecked()) + gsl_vector_short_memcpy(new, nativeHandle) return GslShortVector(new, scope) } override fun equals(other: Any?): Boolean { if (other is GslShortVector) - return gsl_vector_short_equal(nativeHandleChecked(), other.nativeHandleChecked()) == 1 + return gsl_vector_short_equal(nativeHandle, other.nativeHandle) == 1 return super.equals(other) } - override fun close(): Unit = gsl_vector_short_free(nativeHandleChecked()) + override fun close(): Unit = gsl_vector_short_free(nativeHandle) } -internal class GslUShortVector(override val nativeHandle: CPointer, scope: DeferScope) : +internal class GslUShortVector(override val rawNativeHandle: CPointer, scope: AutofreeScope) : GslVector(scope) { override val size: Int - get() = nativeHandleChecked().pointed.size.toInt() + get() = nativeHandle.pointed.size.toInt() override operator fun get(index: Int): UShort = - gsl_vector_ushort_get(nativeHandleChecked(), index.toULong()) + gsl_vector_ushort_get(nativeHandle, index.toULong()) override operator fun set(index: Int, value: UShort): Unit = - gsl_vector_ushort_set(nativeHandleChecked(), index.toULong(), value) + gsl_vector_ushort_set(nativeHandle, index.toULong(), value) override fun copy(): GslUShortVector { val new = requireNotNull(gsl_vector_ushort_alloc(size.toULong())) - gsl_vector_ushort_memcpy(new, nativeHandleChecked()) + gsl_vector_ushort_memcpy(new, nativeHandle) return GslUShortVector(new, scope) } override fun equals(other: Any?): Boolean { if (other is GslUShortVector) - return gsl_vector_ushort_equal(nativeHandleChecked(), other.nativeHandleChecked()) == 1 + return gsl_vector_ushort_equal(nativeHandle, other.nativeHandle) == 1 return super.equals(other) } - override fun close(): Unit = gsl_vector_ushort_free(nativeHandleChecked()) + override fun close(): Unit = gsl_vector_ushort_free(nativeHandle) } -internal class GslLongVector(override val nativeHandle: CPointer, scope: DeferScope) : +internal class GslLongVector(override val rawNativeHandle: CPointer, scope: AutofreeScope) : GslVector(scope) { override val size: Int - get() = nativeHandleChecked().pointed.size.toInt() + get() = nativeHandle.pointed.size.toInt() override operator fun get(index: Int): Long = - gsl_vector_long_get(nativeHandleChecked(), index.toULong()) + gsl_vector_long_get(nativeHandle, index.toULong()) override operator fun set(index: Int, value: Long): Unit = - gsl_vector_long_set(nativeHandleChecked(), index.toULong(), value) + gsl_vector_long_set(nativeHandle, index.toULong(), value) override fun copy(): GslLongVector { val new = requireNotNull(gsl_vector_long_alloc(size.toULong())) - gsl_vector_long_memcpy(new, nativeHandleChecked()) + gsl_vector_long_memcpy(new, nativeHandle) return GslLongVector(new, scope) } override fun equals(other: Any?): Boolean { if (other is GslLongVector) - return gsl_vector_long_equal(nativeHandleChecked(), other.nativeHandleChecked()) == 1 + return gsl_vector_long_equal(nativeHandle, other.nativeHandle) == 1 return super.equals(other) } - override fun close(): Unit = gsl_vector_long_free(nativeHandleChecked()) + override fun close(): Unit = gsl_vector_long_free(nativeHandle) } -internal class GslULongVector(override val nativeHandle: CPointer, scope: DeferScope) : +internal class GslULongVector(override val rawNativeHandle: CPointer, scope: AutofreeScope) : GslVector(scope) { override val size: Int - get() = nativeHandleChecked().pointed.size.toInt() + get() = nativeHandle.pointed.size.toInt() override operator fun get(index: Int): ULong = - gsl_vector_ulong_get(nativeHandleChecked(), index.toULong()) + gsl_vector_ulong_get(nativeHandle, index.toULong()) override operator fun set(index: Int, value: ULong): Unit = - gsl_vector_ulong_set(nativeHandleChecked(), index.toULong(), value) + gsl_vector_ulong_set(nativeHandle, index.toULong(), value) override fun copy(): GslULongVector { val new = requireNotNull(gsl_vector_ulong_alloc(size.toULong())) - gsl_vector_ulong_memcpy(new, nativeHandleChecked()) + gsl_vector_ulong_memcpy(new, nativeHandle) return GslULongVector(new, scope) } override fun equals(other: Any?): Boolean { if (other is GslULongVector) - return gsl_vector_ulong_equal(nativeHandleChecked(), other.nativeHandleChecked()) == 1 + return gsl_vector_ulong_equal(nativeHandle, other.nativeHandle) == 1 return super.equals(other) } - override fun close(): Unit = gsl_vector_ulong_free(nativeHandleChecked()) + override fun close(): Unit = gsl_vector_ulong_free(nativeHandle) } -internal class GslIntVector(override val nativeHandle: CPointer, scope: DeferScope) : +internal class GslIntVector(override val rawNativeHandle: CPointer, scope: AutofreeScope) : GslVector(scope) { override val size: Int - get() = nativeHandleChecked().pointed.size.toInt() + get() = nativeHandle.pointed.size.toInt() override operator fun get(index: Int): Int = - gsl_vector_int_get(nativeHandleChecked(), index.toULong()) + gsl_vector_int_get(nativeHandle, index.toULong()) override operator fun set(index: Int, value: Int): Unit = - gsl_vector_int_set(nativeHandleChecked(), index.toULong(), value) + gsl_vector_int_set(nativeHandle, index.toULong(), value) override fun copy(): GslIntVector { val new = requireNotNull(gsl_vector_int_alloc(size.toULong())) - gsl_vector_int_memcpy(new, nativeHandleChecked()) + gsl_vector_int_memcpy(new, nativeHandle) return GslIntVector(new, scope) } override fun equals(other: Any?): Boolean { if (other is GslIntVector) - return gsl_vector_int_equal(nativeHandleChecked(), other.nativeHandleChecked()) == 1 + return gsl_vector_int_equal(nativeHandle, other.nativeHandle) == 1 return super.equals(other) } - override fun close(): Unit = gsl_vector_int_free(nativeHandleChecked()) + override fun close(): Unit = gsl_vector_int_free(nativeHandle) } -internal class GslUIntVector(override val nativeHandle: CPointer, scope: DeferScope) : +internal class GslUIntVector(override val rawNativeHandle: CPointer, scope: AutofreeScope) : GslVector(scope) { override val size: Int - get() = nativeHandleChecked().pointed.size.toInt() + get() = nativeHandle.pointed.size.toInt() override operator fun get(index: Int): UInt = - gsl_vector_uint_get(nativeHandleChecked(), index.toULong()) + gsl_vector_uint_get(nativeHandle, index.toULong()) override operator fun set(index: Int, value: UInt): Unit = - gsl_vector_uint_set(nativeHandleChecked(), index.toULong(), value) + gsl_vector_uint_set(nativeHandle, index.toULong(), value) override fun copy(): GslUIntVector { val new = requireNotNull(gsl_vector_uint_alloc(size.toULong())) - gsl_vector_uint_memcpy(new, nativeHandleChecked()) + gsl_vector_uint_memcpy(new, nativeHandle) return GslUIntVector(new, scope) } override fun equals(other: Any?): Boolean { if (other is GslUIntVector) - return gsl_vector_uint_equal(nativeHandleChecked(), other.nativeHandleChecked()) == 1 + return gsl_vector_uint_equal(nativeHandle, other.nativeHandle) == 1 return super.equals(other) } - override fun close(): Unit = gsl_vector_uint_free(nativeHandleChecked()) + override fun close(): Unit = gsl_vector_uint_free(nativeHandle) } diff --git a/kmath-gsl/src/nativeTest/kotlin/GslMatrixRealTest.kt b/kmath-gsl/src/nativeTest/kotlin/GslMatrixRealTest.kt index b4b446462..c171f6d0a 100644 --- a/kmath-gsl/src/nativeTest/kotlin/GslMatrixRealTest.kt +++ b/kmath-gsl/src/nativeTest/kotlin/GslMatrixRealTest.kt @@ -3,6 +3,8 @@ package kscience.kmath.gsl import kscience.kmath.linear.RealMatrixContext import kscience.kmath.operations.invoke import kscience.kmath.structures.Matrix +import kscience.kmath.structures.asIterable +import kscience.kmath.structures.toList import kotlin.random.Random import kotlin.test.Test import kotlin.test.assertEquals @@ -39,4 +41,22 @@ internal class GslMatrixRealTest { assertEquals(mat, mat3) assertEquals(mat2, mat3) } + + @Test + fun rows() = GslRealMatrixContext { + val mat = produce(2, 2) { i, j -> i.toDouble() + j } + + mat.rows.asIterable().zip(listOf(listOf(0.0, 1.0), listOf(1.0, 2.0))).forEach { (a, b) -> + assertEquals(a.toList(), b) + } + } + + @Test + fun columns() = GslRealMatrixContext { + val mat = produce(2, 2) { i, j -> i.toDouble() + j } + + mat.columns.asIterable().zip(listOf(listOf(0.0, 1.0), listOf(1.0, 2.0))).forEach { (a, b) -> + assertEquals(a.toList(), b) + } + } }