Add ojalgo module

CHANGELOG.md

@@ -4,15 +4,19 @@
 
 ### Added
 - Metropolis-Hastings sampler
+- Ojalgo `LinearSpace` implementation.
 
 ### Changed
+- attributes-kt moved to a separate project, and the version used is 0.3.0
+- Kotlin 2.1. Now use cross-compilation to deploy macOS targets.
+- Changed `origin` to `cmMatrix` in kmath-commons to avoid property name clash. Expose bidirectional conversion in `CMLinearSpace`
 
 ### Deprecated
 
 ### Removed
 
 ### Fixed
-- Fix EJML to properly treat vectors as columns
+- (BREAKING CHANGE) Fix EJML to properly treat vectors as columns
 
 ### Security
 

benchmarks/build.gradle.kts

@@ -53,7 +53,6 @@ kotlin {
                 implementation(project(":kmath-dimensions"))
                 implementation(project(":kmath-for-real"))
                 implementation(project(":kmath-tensors"))
-                implementation(project(":kmath-multik"))
                 implementation(libs.multik.default)
                 implementation(spclibs.kotlinx.benchmark.runtime)
             }
@@ -61,13 +60,14 @@ kotlin {
 
         val jvmMain by getting {
             dependencies {
-                implementation(project(":kmath-commons"))
+                implementation(projects.kmathCommons)
-                implementation(project(":kmath-ejml"))
+                implementation(projects.kmathEjml)
-                implementation(project(":kmath-nd4j"))
+                implementation(projects.kmathNd4j)
-                implementation(project(":kmath-kotlingrad"))
+                implementation(projects.kmathKotlingrad)
-                implementation(project(":kmath-viktor"))
+                implementation(projects.kmathViktor)
-//                implementation(project(":kmath-jafama"))
+                implementation(projects.kmathOjalgo)
                 implementation(projects.kmath.kmathTensorflow)
+                implementation(projects.kmathMultik)
                 implementation("org.tensorflow:tensorflow-core-platform:0.4.0")
                 implementation("org.nd4j:nd4j-native:1.0.0-M1")
                 //    uncomment if your system supports AVX2

benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/DotBenchmark.kt

@@ -10,10 +10,13 @@ import kotlinx.benchmark.Blackhole
 import kotlinx.benchmark.Scope
 import kotlinx.benchmark.State
 import space.kscience.kmath.commons.linear.CMLinearSpace
+import space.kscience.kmath.commons.linear.CMLinearSpace.dot
 import space.kscience.kmath.ejml.EjmlLinearSpaceDDRM
 import space.kscience.kmath.linear.Float64ParallelLinearSpace
 import space.kscience.kmath.linear.invoke
 import space.kscience.kmath.linear.linearSpace
+import space.kscience.kmath.ojalgo.Ojalgo
+import space.kscience.kmath.ojalgo.linearSpace
 import space.kscience.kmath.operations.Float64Field
 import space.kscience.kmath.tensorflow.produceWithTF
 import space.kscience.kmath.tensors.core.tensorAlgebra
@@ -70,6 +73,11 @@ internal class DotBenchmark {
         blackhole.consume(matrix1 dot matrix2)
     }
 
+    @Benchmark
+    fun ojalgoDot(blackhole: Blackhole) = Ojalgo.R064.linearSpace {
+        blackhole.consume(matrix1 dot matrix2)
+    }
+
     @Benchmark
     fun multikDot(blackhole: Blackhole) = with(multikAlgebra) {
         blackhole.consume(matrix1 dot matrix2)

benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/MatrixInverseBenchmark.kt

@@ -12,10 +12,9 @@ import kotlinx.benchmark.State
 import space.kscience.kmath.commons.linear.CMLinearSpace
 import space.kscience.kmath.commons.linear.lupSolver
 import space.kscience.kmath.ejml.EjmlLinearSpaceDDRM
-import space.kscience.kmath.linear.invoke
+import space.kscience.kmath.linear.*
-import space.kscience.kmath.linear.linearSpace
+import space.kscience.kmath.ojalgo.Ojalgo
-import space.kscience.kmath.linear.lupSolver
+import space.kscience.kmath.ojalgo.linearSpace
-import space.kscience.kmath.linear.parallel
 import space.kscience.kmath.operations.algebra
 import kotlin.random.Random
 
@@ -48,10 +47,14 @@ internal class MatrixInverseBenchmark {
         blackhole.consume(lupSolver().inverse(matrix))
     }
 
-
     @Benchmark
     fun ejmlInverse(blackhole: Blackhole) = EjmlLinearSpaceDDRM {
-        blackhole.consume(matrix.toEjml().inverted())
+        blackhole.consume(matrix.inverted())
+    }
+
+    @Benchmark
+    fun ojalgoInverse(blackhole: Blackhole) = Ojalgo.R064.linearSpace {
+        blackhole.consume(matrix.getOrComputeAttribute(Inverted))
     }
 
 }

build.gradle.kts

@@ -6,8 +6,6 @@ plugins {
     alias(spclibs.plugins.kotlinx.kover)
 }
 
-val attributesVersion by extra("0.2.0")
-
 allprojects {
     repositories {
         maven("https://repo.kotlin.link")
@@ -72,7 +70,7 @@ ksciencePublish {
         useSPCTeam()
     }
     repository("spc", "https://maven.sciprog.center/kscience")
-    sonatype("https://oss.sonatype.org")
+    central()
 }
 
 apiValidation.nonPublicMarkers.add("space.kscience.kmath.UnstableKMathAPI")

examples/src/main/kotlin/space/kscience/kmath/linear/eigenValueDecomposition.kt

@@ -5,6 +5,7 @@
 
 package space.kscience.kmath.linear
 
+import space.kscience.kmath.PerformancePitfall
 import space.kscience.kmath.commons.linear.CMLinearSpace
 import space.kscience.kmath.ejml.EjmlLinearSpaceDDRM
 import space.kscience.kmath.nd.StructureND
@@ -12,6 +13,7 @@ import space.kscience.kmath.operations.algebra
 import space.kscience.kmath.structures.Float64
 import kotlin.random.Random
 
+@OptIn(PerformancePitfall::class)
 fun main() {
     val dim = 46
 
@@ -21,7 +23,7 @@ fun main() {
 
     listOf(CMLinearSpace, EjmlLinearSpaceDDRM).forEach { algebra ->
         with(algebra) {
-            //create a simmetric matrix
+            //create a symmetric matrix
             val matrix = buildMatrix(dim, dim) { row, col ->
                 if (row >= col) u[row, col] else u[col, row]
             }

gradle.properties

@@ -10,7 +10,7 @@ org.gradle.workers.max=4
 kotlin.code.style=official
 kotlin.mpp.stability.nowarn=true
 kotlin.native.ignoreDisabledTargets=true
-org.jetbrains.dokka.experimental.gradle.pluginMode=V2EnabledWithHelpers
+org.jetbrains.dokka.experimental.gradle.pluginMode=V2Enabled
 kotlin.native.enableKlibsCrossCompilation=true
 
-toolsVersion=0.16.0-kotlin-2.1.0
+toolsVersion=0.16.1-kotlin-2.1.0

gradle/libs.versions.toml

@@ -4,6 +4,9 @@ commons-rng = "1.6"
 multik = "0.2.3"
 
 [libraries]
+attributes = "space.kscience:attributes-kt:0.3.0"
+
+commons-math = "org.apache.commons:commons-math3:3.6.1"
 
 commons-rng-simple = { module = "org.apache.commons:commons-rng-simple", version.ref = "commons-rng" }
 commons-rng-sampling = { module = "org.apache.commons:commons-rng-sampling", version.ref = "commons-rng" }
@@ -11,4 +14,6 @@ commons-rng-sampling = { module = "org.apache.commons:commons-rng-sampling", ver
 multik-core = { module = "org.jetbrains.kotlinx:multik-core", version.ref = "multik" }
 multik-default = { module = "org.jetbrains.kotlinx:multik-default", version.ref = "multik" }
 
+ojalgo = "org.ojalgo:ojalgo:55.1.0"
+
 [plugins]

gradle/wrapper/gradle-wrapper.properties

@@ -1,5 +1,5 @@
 distributionBase=GRADLE_USER_HOME
 distributionPath=wrapper/dists
-distributionUrl=https\://services.gradle.org/distributions/gradle-8.7-bin.zip
+distributionUrl=https\://services.gradle.org/distributions/gradle-8.11.1-bin.zip
 zipStoreBase=GRADLE_USER_HOME
 zipStorePath=wrapper/dists

kmath-ast/src/jsTest/kotlin/space/kscience/kmath/ast/utils.kt

@@ -22,7 +22,7 @@ import space.kscience.kmath.wasm.compile as wasmCompile
 import space.kscience.kmath.wasm.compileToExpression as wasmCompileToExpression
 
 @OptIn(UnstableKMathAPI::class)
-private object WasmCompilerTestContext : CompilerTestContext {
+internal object WasmCompilerTestContext : CompilerTestContext {
     override fun MST.compileToExpression(algebra: Int32Ring): Expression<Int> = wasmCompileToExpression(algebra)
     override fun MST.compile(algebra: Int32Ring, arguments: Map<Symbol, Int>): Int = wasmCompile(algebra, arguments)
     override fun MST.compileToExpression(algebra: Float64Field): Expression<Float64> = wasmCompileToExpression(algebra)
@@ -31,7 +31,7 @@ private object WasmCompilerTestContext : CompilerTestContext {
         wasmCompile(algebra, arguments)
 }
 
-private object ESTreeCompilerTestContext : CompilerTestContext {
+internal object ESTreeCompilerTestContext : CompilerTestContext {
     override fun MST.compileToExpression(algebra: Int32Ring): Expression<Int> = estreeCompileToExpression(algebra)
     override fun MST.compile(algebra: Int32Ring, arguments: Map<Symbol, Int>): Int = estreeCompile(algebra, arguments)
     override fun MST.compileToExpression(algebra: Float64Field): Expression<Float64> = estreeCompileToExpression(algebra)

kmath-commons/build.gradle.kts

@@ -13,7 +13,7 @@ kscience {
         api(projects.kmathOptimization)
         api(projects.kmathStat)
         api(projects.kmathFunctions)
-        api("org.apache.commons:commons-math3:3.6.1")
+        api(libs.commons.math)
     }
 }
 

kmath-commons/src/jvmMain/kotlin/space/kscience/kmath/commons/linear/CMLinearSpace.kt

@@ -6,8 +6,6 @@
 package space.kscience.kmath.commons.linear
 
 import org.apache.commons.math3.linear.*
-import org.apache.commons.math3.linear.LUDecomposition
-import org.apache.commons.math3.linear.SingularValueDecomposition
 import space.kscience.attributes.SafeType
 import space.kscience.kmath.UnstableKMathAPI
 import space.kscience.kmath.linear.*
@@ -23,27 +21,26 @@ import space.kscience.kmath.structures.Float64
 import space.kscience.kmath.structures.IntBuffer
 import space.kscience.kmath.structures.asBuffer
 
-public class CMMatrix(public val origin: RealMatrix) : Matrix<Float64> {
+@JvmInline
+public value class CMMatrix(public val cmMatrix: RealMatrix) : Matrix<Float64> {
 
-    override val rowNum: Int get() = origin.rowDimension
+    override val rowNum: Int get() = cmMatrix.rowDimension
-    override val colNum: Int get() = origin.columnDimension
+    override val colNum: Int get() = cmMatrix.columnDimension
 
-    override operator fun get(i: Int, j: Int): Double = origin.getEntry(i, j)
+    override operator fun get(i: Int, j: Int): Double = cmMatrix.getEntry(i, j)
 }
 
 @JvmInline
-public value class CMVector(public val origin: RealVector) : Point<Float64> {
+public value class CMVector(public val cmVector: RealVector) : Point<Float64> {
-    override val size: Int get() = origin.dimension
+    override val size: Int get() = cmVector.dimension
 
-    override operator fun get(index: Int): Double = origin.getEntry(index)
+    override operator fun get(index: Int): Double = cmVector.getEntry(index)
 
-    override operator fun iterator(): Iterator<Float64> = origin.toArray().iterator()
+    override operator fun iterator(): Iterator<Float64> = cmVector.toArray().iterator()
 
     override fun toString(): String = Buffer.toString(this)
 }
 
-public fun RealVector.toPoint(): CMVector = CMVector(this)
-
 public object CMLinearSpace : LinearSpace<Double, Float64Field> {
     override val elementAlgebra: Float64Field get() = Float64Field
 
@@ -59,52 +56,52 @@ public object CMLinearSpace : LinearSpace<Double, Float64Field> {
     }
 
     @OptIn(UnstableKMathAPI::class)
-    public fun Matrix<Float64>.toCM(): CMMatrix = when (val matrix = origin) {
+    public fun Matrix<Float64>.toCM(): RealMatrix = when (val matrix = origin) {
-        is CMMatrix -> matrix
+        is CMMatrix -> matrix.cmMatrix
         else -> {
             //TODO add feature analysis
             val array = Array(rowNum) { i -> DoubleArray(colNum) { j -> get(i, j) } }
-            Array2DRowRealMatrix(array).wrap()
+            Array2DRowRealMatrix(array)
         }
     }
 
-    public fun Point<Float64>.toCM(): CMVector = if (this is CMVector) this else {
+    public fun Point<Float64>.toCM(): RealVector = if (this is CMVector) cmVector else {
-        val array = DoubleArray(size) { this[it] }
+        val array = DoubleArray(size) { get(it) }
-        ArrayRealVector(array).wrap()
+        ArrayRealVector(array)
     }
 
-    internal fun RealMatrix.wrap(): CMMatrix = CMMatrix(this)
+    public fun RealMatrix.asMatrix(): CMMatrix = CMMatrix(this)
-    internal fun RealVector.wrap(): CMVector = CMVector(this)
+    public fun RealVector.asVector(): CMVector = CMVector(this)
 
     override fun buildVector(size: Int, initializer: Float64Field.(Int) -> Double): Point<Float64> =
-        ArrayRealVector(DoubleArray(size) { Float64Field.initializer(it) }).wrap()
+        ArrayRealVector(DoubleArray(size) { Float64Field.initializer(it) }).asVector()
 
     override fun Matrix<Float64>.plus(other: Matrix<Float64>): CMMatrix =
-        toCM().origin.add(other.toCM().origin).wrap()
+        toCM().add(other.toCM()).asMatrix()
 
     override fun Point<Float64>.plus(other: Point<Float64>): CMVector =
-        toCM().origin.add(other.toCM().origin).wrap()
+        toCM().add(other.toCM()).asVector()
 
     override fun Point<Float64>.minus(other: Point<Float64>): CMVector =
-        toCM().origin.subtract(other.toCM().origin).wrap()
+        toCM().subtract(other.toCM()).asVector()
 
     override fun Matrix<Float64>.dot(other: Matrix<Float64>): CMMatrix =
-        toCM().origin.multiply(other.toCM().origin).wrap()
+        toCM().multiply(other.toCM()).asMatrix()
 
     override fun Matrix<Float64>.dot(vector: Point<Float64>): CMVector =
-        toCM().origin.preMultiply(vector.toCM().origin).wrap()
+        toCM().preMultiply(vector.toCM()).asVector()
 
     override operator fun Matrix<Float64>.minus(other: Matrix<Float64>): CMMatrix =
-        toCM().origin.subtract(other.toCM().origin).wrap()
+        toCM().subtract(other.toCM()).asMatrix()
 
     override operator fun Matrix<Float64>.times(value: Double): CMMatrix =
-        toCM().origin.scalarMultiply(value).wrap()
+        toCM().scalarMultiply(value).asMatrix()
 
     override fun Double.times(m: Matrix<Float64>): CMMatrix =
         m * this
 
     override fun Point<Float64>.times(value: Double): CMVector =
-        toCM().origin.mapMultiply(value).wrap()
+        toCM().mapMultiply(value).asVector()
 
     override fun Double.times(v: Point<Float64>): CMVector =
         v * this
@@ -112,36 +109,38 @@ public object CMLinearSpace : LinearSpace<Double, Float64Field> {
     @OptIn(UnstableKMathAPI::class)
     override fun <V, A : StructureAttribute<V>> computeAttribute(structure: Structure2D<Float64>, attribute: A): V? {
 
-        val origin = structure.toCM().origin
+        val origin = structure.toCM()
 
         val raw: Any? = when (attribute) {
             IsDiagonal -> if (origin is DiagonalMatrix) Unit else null
-            Determinant -> LUDecomposition(origin).determinant
+            Determinant -> org.apache.commons.math3.linear.LUDecomposition(origin).determinant
+
+            Inverted -> org.apache.commons.math3.linear.LUDecomposition(origin).solver.inverse.asMatrix()
 
             LUP -> object : LupDecomposition<Float64> {
-                val lup by lazy { LUDecomposition(origin) }
+                val lup by lazy { org.apache.commons.math3.linear.LUDecomposition(origin) }
                 override val pivot: IntBuffer get() = lup.pivot.asBuffer()
-                override val l: Matrix<Float64> get() = lup.l.wrap()
+                override val l: Matrix<Float64> get() = lup.l.asMatrix().withAttribute(LowerTriangular)
-                override val u: Matrix<Float64> get() = lup.u.wrap()
+                override val u: Matrix<Float64> get() = lup.u.asMatrix().withAttribute(UpperTriangular)
             }
 
             Cholesky -> object : CholeskyDecomposition<Float64> {
                 val cmCholesky by lazy { org.apache.commons.math3.linear.CholeskyDecomposition(origin) }
-                override val l: Matrix<Float64> get() = cmCholesky.l.wrap()
+                override val l: Matrix<Float64> get() = cmCholesky.l.asMatrix()
             }
 
             QR -> object : QRDecomposition<Float64> {
                 val cmQr by lazy { org.apache.commons.math3.linear.QRDecomposition(origin) }
-                override val q: Matrix<Float64> get() = cmQr.q.wrap().withAttribute(OrthogonalAttribute)
+                override val q: Matrix<Float64> get() = cmQr.q.asMatrix().withAttribute(OrthogonalAttribute)
-                override val r: Matrix<Float64> get() = cmQr.r.wrap().withAttribute(UpperTriangular)
+                override val r: Matrix<Float64> get() = cmQr.r.asMatrix().withAttribute(UpperTriangular)
             }
 
             SVD -> object : space.kscience.kmath.linear.SingularValueDecomposition<Float64> {
-                val cmSvd by lazy { SingularValueDecomposition(origin) }
+                val cmSvd by lazy { org.apache.commons.math3.linear.SingularValueDecomposition(origin) }
 
-                override val u: Matrix<Float64> get() = cmSvd.u.wrap()
+                override val u: Matrix<Float64> get() = cmSvd.u.asMatrix()
-                override val s: Matrix<Float64> get() = cmSvd.s.wrap()
+                override val s: Matrix<Float64> get() = cmSvd.s.asMatrix()
-                override val v: Matrix<Float64> get() = cmSvd.v.wrap()
+                override val v: Matrix<Float64> get() = cmSvd.v.asMatrix()
                 override val singularValues: Point<Float64> get() = cmSvd.singularValues.asBuffer()
 
             }
@@ -149,8 +148,8 @@ public object CMLinearSpace : LinearSpace<Double, Float64Field> {
             EIG -> object : EigenDecomposition<Float64> {
                 val cmEigen by lazy { org.apache.commons.math3.linear.EigenDecomposition(origin) }
 
-                override val v: Matrix<Float64> get() = cmEigen.v.wrap()
+                override val v: Matrix<Float64> get() = cmEigen.v.asMatrix()
-                override val d: Matrix<Float64> get() = cmEigen.d.wrap()
+                override val d: Matrix<Float64> get() = cmEigen.d.asMatrix()
             }
 
             else -> null
@@ -161,8 +160,8 @@ public object CMLinearSpace : LinearSpace<Double, Float64Field> {
 
 }
 
-public operator fun CMMatrix.plus(other: CMMatrix): CMMatrix = CMMatrix(origin.add(other.origin))
+public operator fun CMMatrix.plus(other: CMMatrix): CMMatrix = CMMatrix(cmMatrix.add(other.cmMatrix))
 
-public operator fun CMMatrix.minus(other: CMMatrix): CMMatrix = CMMatrix(origin.subtract(other.origin))
+public operator fun CMMatrix.minus(other: CMMatrix): CMMatrix = CMMatrix(cmMatrix.subtract(other.cmMatrix))
 
-public infix fun CMMatrix.dot(other: CMMatrix): CMMatrix = CMMatrix(origin.multiply(other.origin))
+public infix fun CMMatrix.dot(other: CMMatrix): CMMatrix = CMMatrix(cmMatrix.multiply(other.cmMatrix))

kmath-commons/src/jvmMain/kotlin/space/kscience/kmath/commons/linear/CMSolver.kt

@@ -19,43 +19,44 @@ public enum class CMDecomposition {
     CHOLESKY
 }
 
-private fun CMLinearSpace.solver(
+private fun CMLinearSpace.cmSolver(
     a: Matrix<Float64>,
     decomposition: CMDecomposition = CMDecomposition.LUP,
 ): DecompositionSolver = when (decomposition) {
-    CMDecomposition.LUP -> LUDecomposition(a.toCM().origin).solver
+    CMDecomposition.LUP -> LUDecomposition(a.toCM()).solver
-    CMDecomposition.RRQR -> RRQRDecomposition(a.toCM().origin).solver
+    CMDecomposition.RRQR -> RRQRDecomposition(a.toCM()).solver
-    CMDecomposition.QR -> QRDecomposition(a.toCM().origin).solver
+    CMDecomposition.QR -> QRDecomposition(a.toCM()).solver
-    CMDecomposition.EIGEN -> EigenDecomposition(a.toCM().origin).solver
+    CMDecomposition.EIGEN -> EigenDecomposition(a.toCM()).solver
-    CMDecomposition.CHOLESKY -> CholeskyDecomposition(a.toCM().origin).solver
+    CMDecomposition.CHOLESKY -> CholeskyDecomposition(a.toCM()).solver
 }
 
 public fun CMLinearSpace.solve(
     a: Matrix<Float64>,
     b: Matrix<Float64>,
     decomposition: CMDecomposition = CMDecomposition.LUP,
-): CMMatrix = solver(a, decomposition).solve(b.toCM().origin).wrap()
+): CMMatrix = cmSolver(a, decomposition).solve(b.toCM()).asMatrix()
 
 public fun CMLinearSpace.solve(
     a: Matrix<Float64>,
     b: Point<Float64>,
     decomposition: CMDecomposition = CMDecomposition.LUP,
-): CMVector = solver(a, decomposition).solve(b.toCM().origin).toPoint()
+): CMVector = cmSolver(a, decomposition).solve(b.toCM()).asVector()
 
 public fun CMLinearSpace.inverse(
     a: Matrix<Float64>,
     decomposition: CMDecomposition = CMDecomposition.LUP,
-): CMMatrix = solver(a, decomposition).inverse.wrap()
+): CMMatrix = cmSolver(a, decomposition).inverse.asMatrix()
 
 
-public fun CMLinearSpace.solver(decomposition: CMDecomposition): LinearSolver<Float64> = object : LinearSolver<Float64> {
+public fun CMLinearSpace.solver(decomposition: CMDecomposition): LinearSolver<Float64> =
-    override fun solve(a: Matrix<Float64>, b: Matrix<Float64>): Matrix<Float64> =
+    object : LinearSolver<Float64> {
-        solver(a, decomposition).solve(b.toCM().origin).wrap()
+        override fun solve(a: Matrix<Float64>, b: Matrix<Float64>): Matrix<Float64> =
+            cmSolver(a, decomposition).solve(b.toCM()).asMatrix()
 
-    override fun solve(a: Matrix<Float64>, b: Point<Float64>): Point<Float64> =
+        override fun solve(a: Matrix<Float64>, b: Point<Float64>): Point<Float64> =
-        solver(a, decomposition).solve(b.toCM().origin).toPoint()
+            cmSolver(a, decomposition).solve(b.toCM()).asVector()
 
-    override fun inverse(matrix: Matrix<Float64>): Matrix<Float64> = solver(matrix, decomposition).inverse.wrap()
+        override fun inverse(matrix: Matrix<Float64>): Matrix<Float64> = cmSolver(matrix, decomposition).inverse.asMatrix()
-}
+    }
 
 public fun CMLinearSpace.lupSolver(): LinearSolver<Float64> = solver((CMDecomposition.LUP))

kmath-commons/src/jvmMain/kotlin/space/kscience/kmath/commons/optimization/CMOptimizer.kt

@@ -34,7 +34,7 @@ public object CMOptimizerEngine : OptimizationAttribute<() -> MultivariateOptimi
  * Specify a Commons-maths optimization engine
  */
 public fun AttributesBuilder<FunctionOptimization<Float64>>.cmEngine(optimizerBuilder: () -> MultivariateOptimizer) {
-    set(CMOptimizerEngine, optimizerBuilder)
+    CMOptimizerEngine(optimizerBuilder)
 }
 
 public object CMOptimizerData : SetAttribute<SymbolIndexer.() -> OptimizationData>

kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/LinearSolver.kt

@@ -6,7 +6,7 @@
 package space.kscience.kmath.linear
 
 /**
- * A group of methods to solve for *X* in equation *X = A<sup>&minus;1</sup> &middot; B*, where *A* and *B* are
+ * A group of methods to solve for $X$ in equation $X = A^{-1} \cdot B$, where $A$ and $B$ are
  * matrices or vectors.
  *
  * @param T the type of items.

kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/MatrixWrapper.kt

@@ -8,6 +8,7 @@ package space.kscience.kmath.linear
 import space.kscience.attributes.Attribute
 import space.kscience.attributes.Attributes
 import space.kscience.attributes.withAttribute
+import space.kscience.attributes.withFlag
 import space.kscience.kmath.UnstableKMathAPI
 import space.kscience.kmath.operations.Ring
 
@@ -47,7 +48,7 @@ public fun <T, A : Attribute<T>> Matrix<T>.withAttribute(
 public fun <T, A : Attribute<Unit>> Matrix<T>.withAttribute(
     attribute: A,
 ): MatrixWrapper<T> = if (this is MatrixWrapper) {
-    MatrixWrapper(origin, attributes.withAttribute(attribute))
+    MatrixWrapper(origin, attributes.withFlag(attribute))
 } else {
     MatrixWrapper(this, Attributes(attribute, Unit))
 }

kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/matrixAttributes.kt

@@ -15,7 +15,7 @@ import space.kscience.kmath.nd.StructureAttribute
  * A marker interface for algebras that operate on matrices
  * @param T type of matrix element
  */
-public interface MatrixScope<T> : AttributeScope<Matrix<T>>, WithType<T>
+public interface MatrixScope<T> : WithType<T>
 
 /**
  * A marker interface representing some properties of matrices or additional transformations of them. Features are used

kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/Structure1D.kt

@@ -117,12 +117,13 @@ internal class MutableBuffer1DWrapper<T>(val buffer: MutableBuffer<T>) : Mutable
 /**
  * Represent a [StructureND] as [Structure1D]. Throw error in case of dimension mismatch.
  */
-public fun <T> StructureND<T>.as1D(): Structure1D<T> = this as? Structure1D<T> ?: if (shape.size == 1) {
+public fun <T> StructureND<T>.as1D(): Structure1D<T> =
-    when (this) {
+    this as? Structure1D<T> ?: if (shape.size == 1) {
-        is BufferND -> Buffer1DWrapper(this.buffer)
+        when (this) {
-        else -> Structure1DWrapper(this)
+            is BufferND -> Buffer1DWrapper(this.buffer)
-    }
+            else -> Structure1DWrapper(this)
-} else error("Can't create 1d-structure from ${shape.size}d-structure")
+        }
+    } else error("Can't create 1d-structure from ${shape.size}d-structure")
 
 public fun <T> MutableStructureND<T>.as1D(): MutableStructure1D<T> =
     this as? MutableStructure1D<T> ?: if (shape.size == 1) {
@@ -133,13 +134,3 @@ public fun <T> MutableStructureND<T>.as1D(): MutableStructure1D<T> =
  * Represent this buffer as 1D structure
  */
 public fun <T> Buffer<T>.asND(): Structure1D<T> = Buffer1DWrapper(this)
-
-/**
- * Expose inner buffer of this [Structure1D] if possible
- */
-internal fun <T : Any> Structure1D<T>.asND(): Buffer<T> = when {
-    this is Buffer1DWrapper<T> -> buffer
-    this is Structure1DWrapper && structure is BufferND<T> -> structure.buffer
-    else -> this
-}
-

kmath-ejml/src/jvmMain/kotlin/space/kscience/kmath/ejml/implementations.kt

@@ -229,6 +229,7 @@ public object EjmlLinearSpaceDDRM : EjmlLinearSpace<Double, Float64Field, DMatri
             }
 
             Determinant -> CommonOps_DDRM.det(origin)
+
             SVD -> object : SingularValueDecomposition<Float64> {
                 val ejmlSvd by lazy {
                     DecompositionFactory_DDRM

kmath-jafama/build.gradle.kts

@@ -1,23 +0,0 @@
-plugins {
-    id("space.kscience.gradle.jvm")
-}
-
-description = "Jafama integration module"
-
-dependencies {
-    api(project(":kmath-core"))
-    api("net.jafama:jafama:2.3.2")
-}
-
-repositories {
-    mavenCentral()
-}
-
-readme {
-    maturity = space.kscience.gradle.Maturity.DEPRECATED
-    propertyByTemplate("artifact", rootProject.file("docs/templates/ARTIFACT-TEMPLATE.md"))
-
-    feature("jafama-double", "src/main/kotlin/space/kscience/kmath/jafama/") {
-        "Double ExtendedField implementations based on Jafama"
-    }
-}

kmath-ojalgo/build.gradle.kts

@@ -0,0 +1,22 @@
+plugins {
+    id("space.kscience.gradle.mpp")
+}
+
+description = "Ojalgo bindings for kmath"
+
+kscience {
+    jvm()
+    jvmMain {
+        api(projects.kmathCore)
+//        api(projects.kmathComplex)
+//        api(projects.kmathCoroutines)
+//        api(projects.kmathOptimization)
+//        api(projects.kmathStat)
+//        api(projects.kmathFunctions)
+        api(libs.ojalgo)
+    }
+}
+
+readme {
+    maturity = space.kscience.gradle.Maturity.PROTOTYPE
+}

kmath-ojalgo/src/jvmMain/kotlin/space/kscience/kmath/ojalgo/Ojalgo.kt

@@ -0,0 +1,35 @@
+/*
+ * Copyright 2018-2025 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.ojalgo
+
+import org.ojalgo.matrix.decomposition.*
+import org.ojalgo.matrix.store.PhysicalStore
+import org.ojalgo.matrix.store.R064Store
+import space.kscience.kmath.operations.Float64Field
+import space.kscience.kmath.operations.Ring
+import space.kscience.kmath.structures.Float64
+
+public class Ojalgo<T : Comparable<T>, A : Ring<T>>(
+    public val elementAlgebra: A,
+    public val storeFactory: PhysicalStore.Factory<T, *>,
+    public val lu: LU.Factory<T>,
+    public val cholesky: Cholesky.Factory<T>,
+    public val qr: QR.Factory<T>,
+    public val svd: SingularValue.Factory<T>,
+    public val eigen: Eigenvalue.Factory<T>
+) {
+    public companion object {
+        public val R064: Ojalgo<Float64, Float64Field> = Ojalgo(
+            elementAlgebra = Float64Field,
+            storeFactory = R064Store.FACTORY,
+            lu = LU.R064,
+            cholesky = Cholesky.R064,
+            qr = QR.R064,
+            svd = SingularValue.R064,
+            eigen = Eigenvalue.R064
+        )
+    }
+}

kmath-ojalgo/src/jvmMain/kotlin/space/kscience/kmath/ojalgo/OjalgoLinearSpace.kt

@@ -0,0 +1,180 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.ojalgo
+
+import org.ojalgo.matrix.store.MatrixStore
+import org.ojalgo.matrix.store.PhysicalStore
+import space.kscience.kmath.UnstableKMathAPI
+import space.kscience.kmath.linear.*
+import space.kscience.kmath.nd.StructureAttribute
+import space.kscience.kmath.operations.Ring
+import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.structures.IntBuffer
+import space.kscience.kmath.structures.asBuffer
+import space.kscience.kmath.structures.asList
+
+@JvmInline
+public value class OjalgoBuffer<T : Comparable<T>>(public val ojalgoMatrix: MatrixStore<T>) : Buffer<T> {
+    override val size: Int get() = ojalgoMatrix.size()
+
+    override fun get(index: Int): T = ojalgoMatrix.get(index.toLong())
+
+    override fun toString(): String = ojalgoMatrix.toString()
+}
+
+@JvmInline
+public value class OjalgoMatrix<T : Comparable<T>>(public val ojalgoVector: MatrixStore<T>) : Matrix<T> {
+    override val rowNum: Int get() = ojalgoVector.rowDim
+
+    override val colNum: Int get() = ojalgoVector.colDim
+
+    override fun get(i: Int, j: Int): T = ojalgoVector.get(i.toLong(), j.toLong())
+}
+
+
+public class OjalgoLinearSpace<T : Comparable<T>, A : Ring<T>>(
+    public val ojalgo: Ojalgo<T, A>
+) : LinearSpace<T, A> {
+
+    override val elementAlgebra: A
+        get() = ojalgo.elementAlgebra
+
+    public fun MatrixStore<T>.asMatrix(): OjalgoMatrix<T> = OjalgoMatrix(this)
+
+    public fun MatrixStore<T>.asVector(): OjalgoBuffer<T> = OjalgoBuffer(this)
+
+    /**
+     * If this matrix is [OjalgoMatrix] return it without conversion, otherwise create new [PhysicalStore]
+     */
+    @OptIn(UnstableKMathAPI::class)
+    public fun Matrix<T>.toOjalgo(): MatrixStore<T> = when (val matrix = origin) {
+        is OjalgoMatrix<T> -> matrix.ojalgoVector
+        else -> ojalgo.storeFactory.make(rowNum.toLong(), colNum.toLong()).apply {
+            for (row in 0 until rowNum) {
+                for (column in 0 until colNum) {
+                    set(row.toLong(), column.toLong(), get(row, column))
+                }
+            }
+        }
+
+    }
+
+    /**
+     * If this vector is [OjalgoBuffer] return it without conversion, otherwise create new [PhysicalStore]
+     */
+    public fun Point<T>.toOjalgo(): MatrixStore<T> =
+        (this as? OjalgoBuffer<T>)?.ojalgoMatrix ?: ojalgo.storeFactory.column(asList())
+
+    override fun buildMatrix(
+        rows: Int,
+        columns: Int,
+        initializer: A.(Int, Int) -> T
+    ): Matrix<T> {
+        val structure: MatrixStore<T> = ojalgo.storeFactory.make(rows.toLong(), columns.toLong()).apply {
+            for (row in 0 until rows) {
+                for (column in 0 until columns) {
+                    set(row.toLong(), column.toLong(), elementAlgebra.initializer(row, column))
+                }
+            }
+        }
+        return OjalgoMatrix(structure)
+    }
+
+    override fun buildVector(size: Int, initializer: A.(Int) -> T): Point<T> {
+        val structure: MatrixStore<T> = ojalgo.storeFactory.column(List(size) { elementAlgebra.initializer(it) })
+        return OjalgoBuffer(structure)
+    }
+
+    @OptIn(UnstableKMathAPI::class)
+    override fun <V, A : StructureAttribute<V>> computeAttribute(
+        structure: Matrix<T>,
+        attribute: A
+    ): V? {
+
+        val origin = structure.toOjalgo()
+
+        val raw: Any? = when (attribute) {
+            Determinant -> ojalgo.lu.make(origin).apply { decompose(origin) }.determinant
+
+            Inverted -> ojalgo.lu.make().apply { decompose(origin) }.inverse.asMatrix()
+
+            LUP -> object : LupDecomposition<T> {
+                val lup by lazy {
+                    ojalgo.lu.make(origin).apply { decompose(origin) }
+                }
+                override val pivot: IntBuffer get() = lup.pivotOrder.asBuffer()
+                override val l: Matrix<T> get() = lup.l.asMatrix().withAttribute(LowerTriangular)
+                override val u: Matrix<T> get() = lup.u.asMatrix().withAttribute(UpperTriangular)
+            }
+
+            Cholesky -> object : CholeskyDecomposition<T> {
+                val cholesky by lazy {
+                    ojalgo.cholesky.make(origin).apply { decompose(origin) }
+                }
+                override val l: Matrix<T> get() = cholesky.l.asMatrix()
+            }
+
+            QR -> object : QRDecomposition<T> {
+                val qr by lazy {
+                    ojalgo.qr.make(origin).apply { decompose(origin) }
+                }
+                override val q: Matrix<T> get() = qr.q.asMatrix().withAttribute(OrthogonalAttribute)
+                override val r: Matrix<T> get() = qr.r.asMatrix().withAttribute(UpperTriangular)
+            }
+
+            SVD -> object : SingularValueDecomposition<T> {
+                val svd by lazy {
+                    ojalgo.svd.make(origin).apply { decompose(origin) }
+                }
+
+                override val u: Matrix<T> get() = svd.u.asMatrix()
+                override val s: Matrix<T> get() = ojalgo.storeFactory.makeDiagonal(svd.singularValues).get().asMatrix()
+                override val v: Matrix<T> get() = svd.v.asMatrix()
+
+                override val singularValues: Point<T>
+                    get() = ojalgo.storeFactory.asFactory1D().make(svd.singularValues).asList().asBuffer()
+
+            }
+
+            EIG -> object : EigenDecomposition<T> {
+                val eigen by lazy {
+                    ojalgo.eigen.make(origin).apply { decompose(origin) }
+                }
+
+                override val v: Matrix<T> get() = eigen.v.asMatrix()
+                override val d: Matrix<T> get() = eigen.d.asMatrix()
+            }
+
+            else -> null
+        }
+        @Suppress("UNCHECKED_CAST")
+        return raw as V?
+    }
+
+    override fun Matrix<T>.times(value: T): OjalgoMatrix<T> = toOjalgo().multiply(value).asMatrix()
+
+    override fun Matrix<T>.dot(vector: Point<T>): OjalgoBuffer<T> = toOjalgo().multiply(vector.toOjalgo()).asVector()
+
+    override fun Matrix<T>.dot(other: Matrix<T>): OjalgoMatrix<T> = toOjalgo().multiply(other.toOjalgo()).asMatrix()
+
+    override fun Point<T>.times(value: T): OjalgoBuffer<T> = toOjalgo().multiply(value).asVector()
+
+    override fun Point<T>.minus(other: Point<T>): OjalgoBuffer<T> = toOjalgo().subtract(other.toOjalgo()).asVector()
+
+    override fun Matrix<T>.minus(other: Matrix<T>): OjalgoMatrix<T> = toOjalgo().subtract(other.toOjalgo()).asMatrix()
+
+    override fun Point<T>.plus(other: Point<T>): OjalgoBuffer<T> = toOjalgo().subtract(other.toOjalgo()).asVector()
+
+    override fun Matrix<T>.plus(other: Matrix<T>): OjalgoMatrix<T> = toOjalgo().add(other.toOjalgo()).asMatrix()
+
+    override fun Point<T>.unaryMinus(): OjalgoBuffer<T> = toOjalgo().negate().asVector()
+
+    override fun Matrix<T>.unaryMinus(): OjalgoMatrix<T> = toOjalgo().negate().asMatrix()
+
+}
+
+public val <T : Comparable<T>, A : Ring<T>> Ojalgo<T, A>.linearSpace: OjalgoLinearSpace<T, A>
+    get() = OjalgoLinearSpace(this)

kmath-ojalgo/src/jvmTest/kotlin/space/kscience/kmath/ojalgo/OjalgoMatrixTest.kt

@@ -0,0 +1,76 @@
+/*
+ * Copyright 2018-2025 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.ojalgo
+
+import org.junit.jupiter.api.Test
+import space.kscience.kmath.PerformancePitfall
+import space.kscience.kmath.UnstableKMathAPI
+import space.kscience.kmath.linear.*
+import space.kscience.kmath.nd.StructureND
+import space.kscience.kmath.structures.Float64
+import kotlin.test.assertEquals
+import kotlin.test.assertTrue
+
+
+@UnstableKMathAPI
+@OptIn(PerformancePitfall::class)
+@Suppress("UNUSED_VARIABLE")
+class OjalgoMatrixTest {
+
+    @Test
+    fun testTranspose() = with(Ojalgo.Companion.R064.linearSpace) {
+        val matrix = one(3, 3)
+        val transposed = matrix.transposed()
+        assertTrue { StructureND.Companion.contentEquals(matrix, transposed) }
+    }
+
+    @Test
+    fun testBuilder() = Ojalgo.Companion.R064.linearSpace {
+        val matrix = matrix(2, 3)(
+            1.0, 0.0, 0.0,
+            0.0, 1.0, 2.0
+        )
+
+        assertEquals(2.0, matrix[1, 2])
+    }
+
+    @Test
+    fun testMatrixExtension() = with(Ojalgo.Companion.R064.linearSpace) {
+        val transitionMatrix: Matrix<Float64> = VirtualMatrix(6, 6) { row, col ->
+            when {
+                col == 0 -> .50
+                row + 1 == col -> .50
+                row == 5 && col == 5 -> 1.0
+                else -> 0.0
+            }
+        }
+
+        infix fun Matrix<Float64>.pow(power: Int): Matrix<Float64> {
+            var res = this
+            repeat(power - 1) {
+                res = res dot this@pow
+            }
+            return res
+        }
+
+        val toTenthPower = transitionMatrix pow 10
+    }
+
+    @Test
+    fun test2DDot() = with(Ojalgo.Companion.R064.linearSpace) {
+        val firstMatrix = buildMatrix(2, 3) { i, j -> (i + j).toDouble() }
+        val secondMatrix = buildMatrix(3, 2) { i, j -> (i + j).toDouble() }
+
+//            val firstMatrix = produce(2, 3) { i, j -> (i + j).toDouble() }
+//            val secondMatrix = produce(3, 2) { i, j -> (i + j).toDouble() }
+        val result = firstMatrix dot secondMatrix
+        assertEquals(2, result.rowNum)
+        assertEquals(2, result.colNum)
+        assertEquals(8.0, result[0, 1])
+        assertEquals(8.0, result[1, 0])
+        assertEquals(14.0, result[1, 1])
+    }
+}

kmath-optimization/src/commonMain/kotlin/space/kscience/kmath/optimization/FunctionOptimization.kt

@@ -12,7 +12,7 @@ import space.kscience.kmath.expressions.Symbol
 public class OptimizationValue<V>(type: SafeType<V>) : PolymorphicAttribute<V>(type)
 
 public inline fun <reified T> AttributesBuilder<FunctionOptimization<T>>.value(value: T) {
-    set(OptimizationValue(safeTypeOf<T>()), value)
+    put(OptimizationValue(safeTypeOf<T>()), value)
 }
 
 public enum class OptimizationDirection {

kmath-optimization/src/commonMain/kotlin/space/kscience/kmath/optimization/OptimizationProblem.kt

@@ -24,7 +24,7 @@ public val <T> OptimizationProblem<T>.startPoint: Map<Symbol, T>
     get() = attributes[OptimizationStartPoint()] ?: error("Starting point not defined in $this")
 
 public fun <T> AttributesBuilder<OptimizationProblem<T>>.startAt(startingPoint: Map<Symbol, T>) {
-    set(OptimizationStartPoint(), startingPoint)
+    put(OptimizationStartPoint(), startingPoint)
 }
 
 
@@ -35,7 +35,7 @@ public class OptimizationCovariance<T> : OptimizationAttribute<NamedMatrix<T>>,
     PolymorphicAttribute<NamedMatrix<T>>(safeTypeOf())
 
 public fun <T> AttributesBuilder<OptimizationProblem<T>>.covariance(covariance: NamedMatrix<T>) {
-    set(OptimizationCovariance(), covariance)
+    put(OptimizationCovariance(), covariance)
 }
 
 
@@ -43,7 +43,7 @@ public class OptimizationResult<T>() : OptimizationAttribute<Map<Symbol, T>>,
     PolymorphicAttribute<Map<Symbol, T>>(safeTypeOf())
 
 public fun <T> AttributesBuilder<OptimizationProblem<T>>.result(result: Map<Symbol, T>) {
-    set(OptimizationResult(), result)
+    put(OptimizationResult(), result)
 }
 
 public val <T> OptimizationProblem<T>.resultOrNull: Map<Symbol, T>? get() = attributes[OptimizationResult()]

kmath-optimization/src/commonMain/kotlin/space/kscience/kmath/optimization/XYFit.kt

@@ -125,9 +125,9 @@ public suspend fun XYColumnarData<Double, Double, Double>.fitWith(
         this,
         modelExpression,
         attributes.modified<XYFit> {
-            set(OptimizationStartPoint(), startingPoint)
+            put(OptimizationStartPoint(), startingPoint)
-            if (!hasAny<OptimizationLog>()) {
+            if (!attributes.hasAny<OptimizationLog>()) {
-                set(OptimizationLog, Loggable.console)
+                put(OptimizationLog, Loggable.console)
             }
         },
         pointToCurveDistance,

settings.gradle.kts

@@ -38,7 +38,6 @@ dependencyResolutionManagement {
 
 include(
     ":test-utils",
-    ":attributes-kt",
     ":kmath-memory",
     ":kmath-complex",
     ":kmath-core",
@@ -61,6 +60,7 @@ include(
     ":kmath-tensors",
     ":kmath-jupyter",
     ":kmath-symja",
+    ":kmath-ojalgo",
     ":examples",
     ":benchmarks",
 )