CHANGELOG.md

@@ -13,6 +13,7 @@
 - Extended operations for ND4J fields
 - Jupyter Notebook integration module (kmath-jupyter)
 - `@PerformancePitfall` annotation to mark possibly slow API
+- Unified architecture for Integration and Optimization using features.
 - `BigInt` operation performance improvement and fixes by @zhelenskiy (#328)
 - Integration between `MST` and Symja `IExpr`
 
@@ -36,8 +37,11 @@
 - Remove Any restriction on polynomials
 - Add `out` variance to type parameters of `StructureND` and its implementations where possible
 - Rename `DifferentiableMstExpression` to `KotlingradExpression`
+- `FeatureSet` now accepts only `Feature`. It is possible to override keys and use interfaces.
+- Use `Symbol` factory function instead of `StringSymbol`
 
 ### Deprecated
+- Specialized `DoubleBufferAlgebra`
 
 ### Removed
 - Nearest in Domain. To be implemented in geometry package.

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

@@ -23,8 +23,8 @@ internal class DotBenchmark {
         const val dim = 1000
 
         //creating invertible matrix
-        val matrix1 = LinearSpace.real.buildMatrix(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 }
+        val matrix1 = LinearSpace.double.buildMatrix(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 }
-        val matrix2 = LinearSpace.real.buildMatrix(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 }
+        val matrix2 = LinearSpace.double.buildMatrix(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 }
 
         val cmMatrix1 = CMLinearSpace { matrix1.toCM() }
         val cmMatrix2 = CMLinearSpace { matrix2.toCM() }
@@ -63,7 +63,7 @@ internal class DotBenchmark {
 
     @Benchmark
     fun realDot(blackhole: Blackhole) {
-        LinearSpace.real {
+        LinearSpace.double {
             blackhole.consume(matrix1 dot matrix2)
         }
     }

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

@@ -14,8 +14,8 @@ import space.kscience.kmath.commons.linear.inverse
 import space.kscience.kmath.ejml.EjmlLinearSpaceDDRM
 import space.kscience.kmath.linear.InverseMatrixFeature
 import space.kscience.kmath.linear.LinearSpace
-import space.kscience.kmath.linear.inverseWithLup
 import space.kscience.kmath.linear.invoke
+import space.kscience.kmath.linear.lupSolver
 import space.kscience.kmath.nd.getFeature
 import kotlin.random.Random
 
@@ -25,7 +25,7 @@ internal class MatrixInverseBenchmark {
         private val random = Random(1224)
         private const val dim = 100
 
-        private val space = LinearSpace.real
+        private val space = LinearSpace.double
 
         //creating invertible matrix
         private val u = space.buildMatrix(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 }
@@ -35,7 +35,7 @@ internal class MatrixInverseBenchmark {
 
     @Benchmark
     fun kmathLupInversion(blackhole: Blackhole) {
-        blackhole.consume(LinearSpace.real.inverseWithLup(matrix))
+        blackhole.consume(LinearSpace.double.lupSolver().inverse(matrix))
     }
 
     @Benchmark

buildSrc/src/main/kotlin/space/kscience/kmath/ejml/codegen/ejmlCodegen.kt

@@ -203,7 +203,7 @@ public object EjmlLinearSpace${ops} : EjmlLinearSpace<${type}, ${kmathAlgebra},
     override fun ${type}.times(v: Point<${type}>): Ejml${type}Vector<${ejmlMatrixType}> = v * this
 
     @UnstableKMathAPI
-    override fun <F : StructureFeature> getFeature(structure: Matrix<${type}>, type: KClass<out F>): F? {
+    override fun <F : StructureFeature> computeFeature(structure: Matrix<${type}>, type: KClass<out F>): F? {
         structure.getFeature(type)?.let { return it }
         val origin = structure.toEjml().origin
 

examples/build.gradle.kts

@@ -20,6 +20,7 @@ dependencies {
     implementation(project(":kmath-coroutines"))
     implementation(project(":kmath-commons"))
     implementation(project(":kmath-complex"))
+    implementation(project(":kmath-optimization"))
     implementation(project(":kmath-stat"))
     implementation(project(":kmath-viktor"))
     implementation(project(":kmath-dimensions"))

examples/src/main/kotlin/space/kscience/kmath/fit/chiSquared.kt

@@ -3,16 +3,17 @@
  * Use of this source code is governed by the Apache 2.0 license that can be found in the LICENSE file.
  */
 
-package space.kscience.kmath.commons.fit
+package space.kscience.kmath.fit
 
 import kotlinx.html.br
 import kotlinx.html.h3
-import space.kscience.kmath.commons.optimization.chiSquared
+import space.kscience.kmath.commons.expressions.DSProcessor
-import space.kscience.kmath.commons.optimization.minimize
+import space.kscience.kmath.commons.optimization.CMOptimizer
 import space.kscience.kmath.distributions.NormalDistribution
+import space.kscience.kmath.expressions.binding
+import space.kscience.kmath.expressions.chiSquaredExpression
 import space.kscience.kmath.expressions.symbol
-import space.kscience.kmath.optimization.FunctionOptimization
+import space.kscience.kmath.optimization.*
-import space.kscience.kmath.optimization.OptimizationResult
 import space.kscience.kmath.real.DoubleVector
 import space.kscience.kmath.real.map
 import space.kscience.kmath.real.step
@@ -22,6 +23,7 @@ import space.kscience.kmath.structures.toList
 import space.kscience.plotly.*
 import space.kscience.plotly.models.ScatterMode
 import space.kscience.plotly.models.TraceValues
+import kotlin.math.abs
 import kotlin.math.pow
 import kotlin.math.sqrt
 
@@ -42,7 +44,7 @@ operator fun TraceValues.invoke(vector: DoubleVector) {
  */
 suspend fun main() {
     //A generator for a normally distributed values
-    val generator = NormalDistribution(2.0, 7.0)
+    val generator = NormalDistribution(0.0, 1.0)
 
     //A chain/flow of random values with the given seed
     val chain = generator.sample(RandomGenerator.default(112667))
@@ -53,7 +55,7 @@ suspend fun main() {
 
 
     //Perform an operation on each x value (much more effective, than numpy)
-    val y = x.map {
+    val y = x.map { it ->
         val value = it.pow(2) + it + 1
         value + chain.next() * sqrt(value)
     }
@@ -64,17 +66,21 @@ suspend fun main() {
     val yErr = y.map { sqrt(it) }//RealVector.same(x.size, sigma)
 
     // compute differentiable chi^2 sum for given model ax^2 + bx + c
-    val chi2 = FunctionOptimization.chiSquared(x, y, yErr) { x1 ->
+    val chi2 = DSProcessor.chiSquaredExpression(x, y, yErr) { arg ->
         //bind variables to autodiff context
         val a = bindSymbol(a)
         val b = bindSymbol(b)
         //Include default value for c if it is not provided as a parameter
         val c = bindSymbolOrNull(c) ?: one
-        a * x1.pow(2) + b * x1 + c
+        a * arg.pow(2) + b * arg + c
     }
 
     //minimize the chi^2 in given starting point. Derivatives are not required, they are already included.
-    val result: OptimizationResult<Double> = chi2.minimize(a to 1.5, b to 0.9, c to 1.0)
+    val result = chi2.optimizeWith(
+        CMOptimizer,
+        mapOf(a to 1.5, b to 0.9, c to 1.0),
+        FunctionOptimizationTarget.MINIMIZE
+    )
 
     //display a page with plot and numerical results
     val page = Plotly.page {
@@ -91,7 +97,7 @@ suspend fun main() {
             scatter {
                 mode = ScatterMode.lines
                 x(x)
-                y(x.map { result.point[a]!! * it.pow(2) + result.point[b]!! * it + 1 })
+                y(x.map { result.resultPoint[a]!! * it.pow(2) + result.resultPoint[b]!! * it + 1 })
                 name = "fit"
             }
         }
@@ -100,7 +106,7 @@ suspend fun main() {
             +"Fit result: $result"
         }
         h3 {
-            +"Chi2/dof = ${result.value / (x.size - 3)}"
+            +"Chi2/dof = ${result.resultValue / (x.size - 3)}"
         }
     }
 

examples/src/main/kotlin/space/kscience/kmath/fit/qowFit.kt

@@ -0,0 +1,106 @@
+/*
+ * Copyright 2018-2021 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.fit
+
+import kotlinx.html.br
+import kotlinx.html.h3
+import space.kscience.kmath.commons.expressions.DSProcessor
+import space.kscience.kmath.data.XYErrorColumnarData
+import space.kscience.kmath.distributions.NormalDistribution
+import space.kscience.kmath.expressions.Symbol
+import space.kscience.kmath.expressions.binding
+import space.kscience.kmath.expressions.symbol
+import space.kscience.kmath.optimization.QowOptimizer
+import space.kscience.kmath.optimization.chiSquaredOrNull
+import space.kscience.kmath.optimization.fitWith
+import space.kscience.kmath.optimization.resultPoint
+import space.kscience.kmath.real.map
+import space.kscience.kmath.real.step
+import space.kscience.kmath.stat.RandomGenerator
+import space.kscience.kmath.structures.asIterable
+import space.kscience.kmath.structures.toList
+import space.kscience.plotly.*
+import space.kscience.plotly.models.ScatterMode
+import kotlin.math.abs
+import kotlin.math.pow
+import kotlin.math.sqrt
+
+// Forward declaration of symbols that will be used in expressions.
+private val a by symbol
+private val b by symbol
+private val c by symbol
+
+
+/**
+ * Least squares fie with auto-differentiation. Uses `kmath-commons` and `kmath-for-real` modules.
+ */
+suspend fun main() {
+    //A generator for a normally distributed values
+    val generator = NormalDistribution(0.0, 1.0)
+
+    //A chain/flow of random values with the given seed
+    val chain = generator.sample(RandomGenerator.default(112667))
+
+
+    //Create a uniformly distributed x values like numpy.arrange
+    val x = 1.0..100.0 step 1.0
+
+
+    //Perform an operation on each x value (much more effective, than numpy)
+    val y = x.map { it ->
+        val value = it.pow(2) + it + 1
+        value + chain.next() * sqrt(value)
+    }
+    // this will also work, but less effective:
+    // val y = x.pow(2)+ x + 1 + chain.nextDouble()
+
+    // create same errors for all xs
+    val yErr = y.map { sqrt(abs(it)) }
+    require(yErr.asIterable().all { it > 0 }) { "All errors must be strictly positive" }
+
+    val result = XYErrorColumnarData.of(x, y, yErr).fitWith(
+        QowOptimizer,
+        DSProcessor,
+        mapOf(a to 0.9, b to 1.2, c to 2.0)
+    ) { arg ->
+        //bind variables to autodiff context
+        val a by binding
+        val b by binding
+        //Include default value for c if it is not provided as a parameter
+        val c = bindSymbolOrNull(c) ?: one
+        a * arg.pow(2) + b * arg + c
+    }
+
+    //display a page with plot and numerical results
+    val page = Plotly.page {
+        plot {
+            scatter {
+                mode = ScatterMode.markers
+                x(x)
+                y(y)
+                error_y {
+                    array = yErr.toList()
+                }
+                name = "data"
+            }
+            scatter {
+                mode = ScatterMode.lines
+                x(x)
+                y(x.map { result.model(result.resultPoint + (Symbol.x to it)) })
+                name = "fit"
+            }
+        }
+        br()
+        h3 {
+            +"Fit result: ${result.resultPoint}"
+        }
+        h3 {
+            +"Chi2/dof = ${result.chiSquaredOrNull!! / (x.size - 3)}"
+        }
+    }
+
+    page.makeFile()
+}

examples/src/main/kotlin/space/kscience/kmath/functions/matrixIntegration.kt

@@ -22,7 +22,7 @@ fun main(): Unit = DoubleField {
         }
 
         //Define a function in a nd space
-        val function: (Double) -> StructureND<Double> = { x: Double -> 3 * number(x).pow(2) + 2 * diagonal(x) + 1 }
+        val function: (Double) -> StructureND<Double> = { x: Double -> 3 * x.pow(2) + 2 * diagonal(x) + 1 }
 
         //get the result of the integration
         val result = gaussIntegrator.integrate(0.0..10.0, function = function)

examples/src/main/kotlin/space/kscience/kmath/structures/buffers.kt

@@ -0,0 +1,22 @@
+/*
+ * Copyright 2018-2021 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.structures
+
+import space.kscience.kmath.operations.DoubleField
+import space.kscience.kmath.operations.bufferAlgebra
+import space.kscience.kmath.operations.produce
+
+inline fun <reified R : Any> MutableBuffer.Companion.same(
+    n: Int,
+    value: R
+): MutableBuffer<R> = auto(n) { value }
+
+
+fun main() {
+    with(DoubleField.bufferAlgebra(5)) {
+        println(number(2.0) + produce(1, 2, 3, 4, 5))
+    }
+}

kmath-ast/src/commonMain/kotlin/space/kscience/kmath/ast/parser.kt

@@ -17,7 +17,7 @@ import com.github.h0tk3y.betterParse.lexer.regexToken
 import com.github.h0tk3y.betterParse.parser.ParseResult
 import com.github.h0tk3y.betterParse.parser.Parser
 import space.kscience.kmath.expressions.MST
-import space.kscience.kmath.expressions.StringSymbol
+import space.kscience.kmath.expressions.Symbol
 import space.kscience.kmath.operations.FieldOperations
 import space.kscience.kmath.operations.GroupOperations
 import space.kscience.kmath.operations.PowerOperations
@@ -43,7 +43,7 @@ public object ArithmeticsEvaluator : Grammar<MST>() {
     private val ws: Token by regexToken("\\s+".toRegex(), ignore = true)
 
     private val number: Parser<MST> by num use { MST.Numeric(text.toDouble()) }
-    private val singular: Parser<MST> by id use { StringSymbol(text) }
+    private val singular: Parser<MST> by id use { Symbol(text) }
 
     private val unaryFunction: Parser<MST> by (id and -lpar and parser(ArithmeticsEvaluator::subSumChain) and -rpar)
         .map { (id, term) -> MST.Unary(id.text, term) }

kmath-ast/src/jvmMain/kotlin/space/kscience/kmath/asm/internal/mapIntrinsics.kt

@@ -7,7 +7,6 @@
 
 package space.kscience.kmath.asm.internal
 
-import space.kscience.kmath.expressions.StringSymbol
 import space.kscience.kmath.expressions.Symbol
 
 /**
@@ -15,4 +14,4 @@ import space.kscience.kmath.expressions.Symbol
  *
  * @author Iaroslav Postovalov
  */
-internal fun <V> Map<Symbol, V>.getOrFail(key: String): V = getValue(StringSymbol(key))
+internal fun <V> Map<Symbol, V>.getOrFail(key: String): V = getValue(Symbol(key))

kmath-commons/build.gradle.kts

@@ -9,6 +9,7 @@ dependencies {
     api(project(":kmath-core"))
     api(project(":kmath-complex"))
     api(project(":kmath-coroutines"))
+    api(project(":kmath-optimization"))
     api(project(":kmath-stat"))
     api(project(":kmath-functions"))
     api("org.apache.commons:commons-math3:3.6.1")

kmath-commons/src/main/kotlin/space/kscience/kmath/commons/expressions/DerivativeStructureExpression.kt

@@ -103,12 +103,15 @@ public class DerivativeStructureField(
     override operator fun DerivativeStructure.minus(b: Number): DerivativeStructure = subtract(b.toDouble())
     override operator fun Number.plus(b: DerivativeStructure): DerivativeStructure = b + this
     override operator fun Number.minus(b: DerivativeStructure): DerivativeStructure = b - this
+}
 
-    public companion object :
+/**
-        AutoDiffProcessor<Double, DerivativeStructure, DerivativeStructureField, Expression<Double>> {
+ * Auto-diff processor based on Commons-math [DerivativeStructure]
-        override fun process(function: DerivativeStructureField.() -> DerivativeStructure): DifferentiableExpression<Double> =
+ */
-            DerivativeStructureExpression(function)
+public object DSProcessor : AutoDiffProcessor<Double, DerivativeStructure, DerivativeStructureField> {
-    }
+     override fun differentiate(
+        function: DerivativeStructureField.() -> DerivativeStructure,
+    ): DerivativeStructureExpression = DerivativeStructureExpression(function)
 }
 
 /**

kmath-commons/src/main/kotlin/space/kscience/kmath/commons/integration/CMGaussRuleIntegrator.kt

@@ -16,7 +16,7 @@ public class CMGaussRuleIntegrator(
     private var type: GaussRule = GaussRule.LEGANDRE,
 ) : UnivariateIntegrator<Double> {
 
-    override fun integrate(integrand: UnivariateIntegrand<Double>): UnivariateIntegrand<Double> {
+    override fun process(integrand: UnivariateIntegrand<Double>): UnivariateIntegrand<Double> {
         val range = integrand.getFeature<IntegrationRange>()?.range
             ?: error("Integration range is not provided")
         val integrator: GaussIntegrator = getIntegrator(range)
@@ -76,8 +76,8 @@ public class CMGaussRuleIntegrator(
             numPoints: Int = 100,
             type: GaussRule = GaussRule.LEGANDRE,
             function: (Double) -> Double,
-        ): Double = CMGaussRuleIntegrator(numPoints, type).integrate(
+        ): Double = CMGaussRuleIntegrator(numPoints, type).process(
             UnivariateIntegrand(function, IntegrationRange(range))
-        ).valueOrNull!!
+        ).value
     }
 }

kmath-commons/src/main/kotlin/space/kscience/kmath/commons/integration/CMIntegrator.kt

@@ -18,7 +18,7 @@ public class CMIntegrator(
     public val integratorBuilder: (Integrand) -> org.apache.commons.math3.analysis.integration.UnivariateIntegrator,
 ) : UnivariateIntegrator<Double> {
 
-    override fun integrate(integrand: UnivariateIntegrand<Double>): UnivariateIntegrand<Double> {
+    override fun process(integrand: UnivariateIntegrand<Double>): UnivariateIntegrand<Double> {
         val integrator = integratorBuilder(integrand)
         val maxCalls = integrand.getFeature<IntegrandMaxCalls>()?.maxCalls ?: defaultMaxCalls
         val remainingCalls = maxCalls - integrand.calls

kmath-commons/src/main/kotlin/space/kscience/kmath/commons/linear/CMMatrix.kt

@@ -95,7 +95,7 @@ public object CMLinearSpace : LinearSpace<Double, DoubleField> {
         v * this
 
     @UnstableKMathAPI
-    override fun <F : StructureFeature> getFeature(structure: Matrix<Double>, type: KClass<out F>): F? {
+    override fun <F : StructureFeature> computeFeature(structure: Matrix<Double>, type: KClass<out F>): F? {
         //Return the feature if it is intrinsic to the structure
         structure.getFeature(type)?.let { return it }
 
@@ -109,22 +109,22 @@ public object CMLinearSpace : LinearSpace<Double, DoubleField> {
                 LupDecompositionFeature<Double> {
                 private val lup by lazy { LUDecomposition(origin) }
                 override val determinant: Double by lazy { lup.determinant }
-                override val l: Matrix<Double> by lazy { CMMatrix(lup.l) + LFeature }
+                override val l: Matrix<Double> by lazy<Matrix<Double>> { CMMatrix(lup.l).withFeature(LFeature) }
-                override val u: Matrix<Double> by lazy { CMMatrix(lup.u) + UFeature }
+                override val u: Matrix<Double> by lazy<Matrix<Double>> { CMMatrix(lup.u).withFeature(UFeature) }
                 override val p: Matrix<Double> by lazy { CMMatrix(lup.p) }
             }
 
             CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature<Double> {
-                override val l: Matrix<Double> by lazy {
+                override val l: Matrix<Double> by lazy<Matrix<Double>> {
                     val cholesky = CholeskyDecomposition(origin)
-                    CMMatrix(cholesky.l) + LFeature
+                    CMMatrix(cholesky.l).withFeature(LFeature)
                 }
             }
 
             QRDecompositionFeature::class -> object : QRDecompositionFeature<Double> {
                 private val qr by lazy { QRDecomposition(origin) }
-                override val q: Matrix<Double> by lazy { CMMatrix(qr.q) + OrthogonalFeature }
+                override val q: Matrix<Double> by lazy<Matrix<Double>> { CMMatrix(qr.q).withFeature(OrthogonalFeature) }
-                override val r: Matrix<Double> by lazy { CMMatrix(qr.r) + UFeature }
+                override val r: Matrix<Double> by lazy<Matrix<Double>> { CMMatrix(qr.r).withFeature(UFeature) }
             }
 
             SingularValueDecompositionFeature::class -> object : SingularValueDecompositionFeature<Double> {

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

@@ -6,6 +6,7 @@
 package space.kscience.kmath.commons.linear
 
 import org.apache.commons.math3.linear.*
+import space.kscience.kmath.linear.LinearSolver
 import space.kscience.kmath.linear.Matrix
 import space.kscience.kmath.linear.Point
 
@@ -44,3 +45,12 @@ public fun CMLinearSpace.inverse(
     a: Matrix<Double>,
     decomposition: CMDecomposition = CMDecomposition.LUP,
 ): CMMatrix = solver(a, decomposition).inverse.wrap()
+
+
+public fun CMLinearSpace.solver(decomposition: CMDecomposition): LinearSolver<Double> = object : LinearSolver<Double> {
+    override fun solve(a: Matrix<Double>, b: Matrix<Double>): Matrix<Double> = solve(a, b, decomposition)
+
+    override fun solve(a: Matrix<Double>, b: Point<Double>): Point<Double> = solve(a, b, decomposition)
+
+    override fun inverse(matrix: Matrix<Double>): Matrix<Double> = inverse(matrix, decomposition)
+}

kmath-commons/src/main/kotlin/space/kscience/kmath/commons/optimization/CMOptimization.kt

@@ -1,126 +0,0 @@
-/*
- * Copyright 2018-2021 KMath contributors.
- * Use of this source code is governed by the Apache 2.0 license that can be found in the LICENSE file.
- */
-
-package space.kscience.kmath.commons.optimization
-
-import org.apache.commons.math3.optim.*
-import org.apache.commons.math3.optim.nonlinear.scalar.GoalType
-import org.apache.commons.math3.optim.nonlinear.scalar.MultivariateOptimizer
-import org.apache.commons.math3.optim.nonlinear.scalar.ObjectiveFunction
-import org.apache.commons.math3.optim.nonlinear.scalar.ObjectiveFunctionGradient
-import org.apache.commons.math3.optim.nonlinear.scalar.gradient.NonLinearConjugateGradientOptimizer
-import org.apache.commons.math3.optim.nonlinear.scalar.noderiv.AbstractSimplex
-import org.apache.commons.math3.optim.nonlinear.scalar.noderiv.NelderMeadSimplex
-import org.apache.commons.math3.optim.nonlinear.scalar.noderiv.SimplexOptimizer
-import space.kscience.kmath.expressions.*
-import space.kscience.kmath.misc.UnstableKMathAPI
-import space.kscience.kmath.optimization.*
-import kotlin.collections.set
-import kotlin.reflect.KClass
-
-public operator fun PointValuePair.component1(): DoubleArray = point
-public operator fun PointValuePair.component2(): Double = value
-
-@OptIn(UnstableKMathAPI::class)
-public class CMOptimization(
-    override val symbols: List<Symbol>,
-) : FunctionOptimization<Double>, NoDerivFunctionOptimization<Double>, SymbolIndexer, OptimizationFeature {
-
-    private val optimizationData: HashMap<KClass<out OptimizationData>, OptimizationData> = HashMap()
-    private var optimizerBuilder: (() -> MultivariateOptimizer)? = null
-    public var convergenceChecker: ConvergenceChecker<PointValuePair> = SimpleValueChecker(
-        DEFAULT_RELATIVE_TOLERANCE,
-        DEFAULT_ABSOLUTE_TOLERANCE,
-        DEFAULT_MAX_ITER
-    )
-
-    override var maximize: Boolean
-        get() = optimizationData[GoalType::class] == GoalType.MAXIMIZE
-        set(value) {
-            optimizationData[GoalType::class] = if (value) GoalType.MAXIMIZE else GoalType.MINIMIZE
-        }
-
-    public fun addOptimizationData(data: OptimizationData) {
-        optimizationData[data::class] = data
-    }
-
-    init {
-        addOptimizationData(MaxEval.unlimited())
-    }
-
-    public fun exportOptimizationData(): List<OptimizationData> = optimizationData.values.toList()
-
-    override fun initialGuess(map: Map<Symbol, Double>) {
-        addOptimizationData(InitialGuess(map.toDoubleArray()))
-    }
-
-    override fun function(expression: Expression<Double>) {
-        val objectiveFunction = ObjectiveFunction {
-            val args = it.toMap()
-            expression(args)
-        }
-        addOptimizationData(objectiveFunction)
-    }
-
-    override fun diffFunction(expression: DifferentiableExpression<Double>) {
-        function(expression)
-        val gradientFunction = ObjectiveFunctionGradient {
-            val args = it.toMap()
-            DoubleArray(symbols.size) { index ->
-                expression.derivative(symbols[index])(args)
-            }
-        }
-        addOptimizationData(gradientFunction)
-        if (optimizerBuilder == null) {
-            optimizerBuilder = {
-                NonLinearConjugateGradientOptimizer(
-                    NonLinearConjugateGradientOptimizer.Formula.FLETCHER_REEVES,
-                    convergenceChecker
-                )
-            }
-        }
-    }
-
-    public fun simplex(simplex: AbstractSimplex) {
-        addOptimizationData(simplex)
-        //Set optimization builder to simplex if it is not present
-        if (optimizerBuilder == null) {
-            optimizerBuilder = { SimplexOptimizer(convergenceChecker) }
-        }
-    }
-
-    public fun simplexSteps(steps: Map<Symbol, Double>) {
-        simplex(NelderMeadSimplex(steps.toDoubleArray()))
-    }
-
-    public fun goal(goalType: GoalType) {
-        addOptimizationData(goalType)
-    }
-
-    public fun optimizer(block: () -> MultivariateOptimizer) {
-        optimizerBuilder = block
-    }
-
-    override fun update(result: OptimizationResult<Double>) {
-        initialGuess(result.point)
-    }
-
-    override fun optimize(): OptimizationResult<Double> {
-        val optimizer = optimizerBuilder?.invoke() ?: error("Optimizer not defined")
-        val (point, value) = optimizer.optimize(*optimizationData.values.toTypedArray())
-        return OptimizationResult(point.toMap(), value, setOf(this))
-    }
-
-    public companion object : OptimizationProblemFactory<Double, CMOptimization> {
-        public const val DEFAULT_RELATIVE_TOLERANCE: Double = 1e-4
-        public const val DEFAULT_ABSOLUTE_TOLERANCE: Double = 1e-4
-        public const val DEFAULT_MAX_ITER: Int = 1000
-
-        override fun build(symbols: List<Symbol>): CMOptimization = CMOptimization(symbols)
-    }
-}
-
-public fun CMOptimization.initialGuess(vararg pairs: Pair<Symbol, Double>): Unit = initialGuess(pairs.toMap())
-public fun CMOptimization.simplexSteps(vararg pairs: Pair<Symbol, Double>): Unit = simplexSteps(pairs.toMap())

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

@@ -0,0 +1,145 @@
+/*
+ * Copyright 2018-2021 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.
+ */
+@file:OptIn(UnstableKMathAPI::class)
+package space.kscience.kmath.commons.optimization
+
+import org.apache.commons.math3.optim.*
+import org.apache.commons.math3.optim.nonlinear.scalar.GoalType
+import org.apache.commons.math3.optim.nonlinear.scalar.MultivariateOptimizer
+import org.apache.commons.math3.optim.nonlinear.scalar.ObjectiveFunction
+import org.apache.commons.math3.optim.nonlinear.scalar.ObjectiveFunctionGradient
+import org.apache.commons.math3.optim.nonlinear.scalar.gradient.NonLinearConjugateGradientOptimizer
+import org.apache.commons.math3.optim.nonlinear.scalar.noderiv.NelderMeadSimplex
+import org.apache.commons.math3.optim.nonlinear.scalar.noderiv.SimplexOptimizer
+import space.kscience.kmath.expressions.Symbol
+import space.kscience.kmath.expressions.SymbolIndexer
+import space.kscience.kmath.expressions.derivative
+import space.kscience.kmath.expressions.withSymbols
+import space.kscience.kmath.misc.UnstableKMathAPI
+import space.kscience.kmath.misc.log
+import space.kscience.kmath.optimization.*
+import kotlin.collections.set
+import kotlin.reflect.KClass
+
+public operator fun PointValuePair.component1(): DoubleArray = point
+public operator fun PointValuePair.component2(): Double = value
+
+public class CMOptimizerEngine(public val optimizerBuilder: () -> MultivariateOptimizer) : OptimizationFeature {
+    override fun toString(): String = "CMOptimizer($optimizerBuilder)"
+}
+
+/**
+ * Specify a Commons-maths optimization engine
+ */
+public fun FunctionOptimizationBuilder<Double>.cmEngine(optimizerBuilder: () -> MultivariateOptimizer) {
+    addFeature(CMOptimizerEngine(optimizerBuilder))
+}
+
+public class CMOptimizerData(public val data: List<SymbolIndexer.() -> OptimizationData>) : OptimizationFeature {
+    public constructor(vararg data: (SymbolIndexer.() -> OptimizationData)) : this(data.toList())
+
+    override fun toString(): String = "CMOptimizerData($data)"
+}
+
+/**
+ * Specify Commons-maths optimization data.
+ */
+public fun FunctionOptimizationBuilder<Double>.cmOptimizationData(data: SymbolIndexer.() -> OptimizationData) {
+    updateFeature<CMOptimizerData> {
+        val newData = (it?.data ?: emptyList()) + data
+        CMOptimizerData(newData)
+    }
+}
+
+public fun FunctionOptimizationBuilder<Double>.simplexSteps(vararg steps: Pair<Symbol, Double>) {
+    //TODO use convergence checker from features
+    cmEngine { SimplexOptimizer(CMOptimizer.defaultConvergenceChecker) }
+    cmOptimizationData { NelderMeadSimplex(mapOf(*steps).toDoubleArray()) }
+}
+
+@OptIn(UnstableKMathAPI::class)
+public object CMOptimizer : Optimizer<Double, FunctionOptimization<Double>> {
+
+    public const val DEFAULT_RELATIVE_TOLERANCE: Double = 1e-4
+    public const val DEFAULT_ABSOLUTE_TOLERANCE: Double = 1e-4
+    public const val DEFAULT_MAX_ITER: Int = 1000
+
+    public val defaultConvergenceChecker: SimpleValueChecker = SimpleValueChecker(
+        DEFAULT_RELATIVE_TOLERANCE,
+        DEFAULT_ABSOLUTE_TOLERANCE,
+        DEFAULT_MAX_ITER
+    )
+
+
+    override suspend fun optimize(
+        problem: FunctionOptimization<Double>,
+    ): FunctionOptimization<Double> {
+        val startPoint = problem.startPoint
+
+        val parameters = problem.getFeature<OptimizationParameters>()?.symbols
+            ?: problem.getFeature<OptimizationStartPoint<Double>>()?.point?.keys
+            ?: startPoint.keys
+
+
+        withSymbols(parameters) {
+            val convergenceChecker: ConvergenceChecker<PointValuePair> = SimpleValueChecker(
+                DEFAULT_RELATIVE_TOLERANCE,
+                DEFAULT_ABSOLUTE_TOLERANCE,
+                DEFAULT_MAX_ITER
+            )
+
+            val cmOptimizer: MultivariateOptimizer = problem.getFeature<CMOptimizerEngine>()?.optimizerBuilder?.invoke()
+                ?: NonLinearConjugateGradientOptimizer(
+                    NonLinearConjugateGradientOptimizer.Formula.FLETCHER_REEVES,
+                    convergenceChecker
+                )
+
+            val optimizationData: HashMap<KClass<out OptimizationData>, OptimizationData> = HashMap()
+
+            fun addOptimizationData(data: OptimizationData) {
+                optimizationData[data::class] = data
+            }
+
+            addOptimizationData(MaxEval.unlimited())
+            addOptimizationData(InitialGuess(startPoint.toDoubleArray()))
+
+            //fun exportOptimizationData(): List<OptimizationData> = optimizationData.values.toList()
+
+            val objectiveFunction = ObjectiveFunction {
+                val args = startPoint + it.toMap()
+                val res = problem.expression(args)
+                res
+            }
+            addOptimizationData(objectiveFunction)
+
+            val gradientFunction = ObjectiveFunctionGradient {
+                val args = startPoint + it.toMap()
+                val res = DoubleArray(symbols.size) { index ->
+                    problem.expression.derivative(symbols[index])(args)
+                }
+                res
+            }
+            addOptimizationData(gradientFunction)
+
+            val logger = problem.getFeature<OptimizationLog>()
+
+            for (feature in problem.features) {
+                when (feature) {
+                    is CMOptimizerData -> feature.data.forEach { dataBuilder ->
+                        addOptimizationData(dataBuilder())
+                    }
+                    is FunctionOptimizationTarget -> when (feature) {
+                        FunctionOptimizationTarget.MAXIMIZE -> addOptimizationData(GoalType.MAXIMIZE)
+                        FunctionOptimizationTarget.MINIMIZE -> addOptimizationData(GoalType.MINIMIZE)
+                    }
+                    else -> logger?.log { "The feature $feature is unused in optimization" }
+                }
+            }
+
+            val (point, value) = cmOptimizer.optimize(*optimizationData.values.toTypedArray())
+            return problem.withFeatures(OptimizationResult(point.toMap()), OptimizationValue(value))
+        }
+    }
+}

kmath-commons/src/main/kotlin/space/kscience/kmath/commons/optimization/cmFit.kt

@@ -1,73 +0,0 @@
-/*
- * Copyright 2018-2021 KMath contributors.
- * Use of this source code is governed by the Apache 2.0 license that can be found in the LICENSE file.
- */
-
-package space.kscience.kmath.commons.optimization
-
-import org.apache.commons.math3.analysis.differentiation.DerivativeStructure
-import space.kscience.kmath.commons.expressions.DerivativeStructureField
-import space.kscience.kmath.expressions.DifferentiableExpression
-import space.kscience.kmath.expressions.Expression
-import space.kscience.kmath.expressions.Symbol
-import space.kscience.kmath.optimization.FunctionOptimization
-import space.kscience.kmath.optimization.OptimizationResult
-import space.kscience.kmath.optimization.noDerivOptimizeWith
-import space.kscience.kmath.optimization.optimizeWith
-import space.kscience.kmath.structures.Buffer
-import space.kscience.kmath.structures.asBuffer
-
-/**
- * Generate a chi squared expression from given x-y-sigma data and inline model. Provides automatic differentiation
- */
-public fun FunctionOptimization.Companion.chiSquared(
-    x: Buffer<Double>,
-    y: Buffer<Double>,
-    yErr: Buffer<Double>,
-    model: DerivativeStructureField.(x: DerivativeStructure) -> DerivativeStructure,
-): DifferentiableExpression<Double> = chiSquared(DerivativeStructureField, x, y, yErr, model)
-
-/**
- * Generate a chi squared expression from given x-y-sigma data and inline model. Provides automatic differentiation
- */
-public fun FunctionOptimization.Companion.chiSquared(
-    x: Iterable<Double>,
-    y: Iterable<Double>,
-    yErr: Iterable<Double>,
-    model: DerivativeStructureField.(x: DerivativeStructure) -> DerivativeStructure,
-): DifferentiableExpression<Double> = chiSquared(
-    DerivativeStructureField,
-    x.toList().asBuffer(),
-    y.toList().asBuffer(),
-    yErr.toList().asBuffer(),
-    model
-)
-
-/**
- * Optimize expression without derivatives
- */
-public fun Expression<Double>.optimize(
-    vararg symbols: Symbol,
-    configuration: CMOptimization.() -> Unit,
-): OptimizationResult<Double> = noDerivOptimizeWith(CMOptimization, symbols = symbols, configuration)
-
-/**
- * Optimize differentiable expression
- */
-public fun DifferentiableExpression<Double>.optimize(
-    vararg symbols: Symbol,
-    configuration: CMOptimization.() -> Unit,
-): OptimizationResult<Double> = optimizeWith(CMOptimization, symbols = symbols, configuration)
-
-public fun DifferentiableExpression<Double>.minimize(
-    vararg startPoint: Pair<Symbol, Double>,
-    configuration: CMOptimization.() -> Unit = {},
-): OptimizationResult<Double> {
-    val symbols = startPoint.map { it.first }.toTypedArray()
-    return optimize(*symbols){
-        maximize = false
-        initialGuess(startPoint.toMap())
-        diffFunction(this@minimize)
-        configuration()
-    }
-}

kmath-commons/src/test/kotlin/space/kscience/kmath/commons/expressions/DerivativeStructureExpressionTest.kt

@@ -42,8 +42,8 @@ internal class AutoDiffTest {
     @Test
     fun autoDifTest() {
         val f = DerivativeStructureExpression {
-            val x by binding()
+            val x by binding
-            val y by binding()
+            val y by binding
             x.pow(2) + 2 * x * y + y.pow(2) + 1
         }
 

kmath-commons/src/test/kotlin/space/kscience/kmath/commons/optimization/OptimizeTest.kt

@@ -6,43 +6,42 @@
 package space.kscience.kmath.commons.optimization
 
 import kotlinx.coroutines.runBlocking
+import space.kscience.kmath.commons.expressions.DSProcessor
 import space.kscience.kmath.commons.expressions.DerivativeStructureExpression
 import space.kscience.kmath.distributions.NormalDistribution
+import space.kscience.kmath.expressions.Symbol.Companion.x
+import space.kscience.kmath.expressions.Symbol.Companion.y
+import space.kscience.kmath.expressions.chiSquaredExpression
 import space.kscience.kmath.expressions.symbol
-import space.kscience.kmath.optimization.FunctionOptimization
+import space.kscience.kmath.optimization.*
 import space.kscience.kmath.stat.RandomGenerator
+import space.kscience.kmath.structures.DoubleBuffer
+import space.kscience.kmath.structures.asBuffer
+import space.kscience.kmath.structures.map
 import kotlin.math.pow
 import kotlin.test.Test
 
 internal class OptimizeTest {
-    val x by symbol
-    val y by symbol
-
     val normal = DerivativeStructureExpression {
-        exp(-bindSymbol(x).pow(2) / 2) + exp(-bindSymbol(y)
+        exp(-bindSymbol(x).pow(2) / 2) + exp(-bindSymbol(y).pow(2) / 2)
-            .pow(2) / 2)
     }
 
     @Test
-    fun testGradientOptimization() {
+    fun testGradientOptimization() = runBlocking {
-        val result = normal.optimize(x, y) {
+        val result = normal.optimizeWith(CMOptimizer, x to 1.0, y to 1.0)
-            initialGuess(x to 1.0, y to 1.0)
+        println(result.resultPoint)
-            // no need to select optimizer. Gradient optimizer is used by default because gradients are provided by function.
+        println(result.resultValue)
-        }
-        println(result.point)
-        println(result.value)
     }
 
     @Test
-    fun testSimplexOptimization() {
+    fun testSimplexOptimization() = runBlocking {
-        val result = normal.optimize(x, y) {
+        val result = normal.optimizeWith(CMOptimizer, x to 1.0, y to 1.0) {
-            initialGuess(x to 1.0, y to 1.0)
             simplexSteps(x to 2.0, y to 0.5)
             //this sets simplex optimizer
         }
 
-        println(result.point)
+        println(result.resultPoint)
-        println(result.value)
+        println(result.resultValue)
     }
 
     @Test
@@ -54,21 +53,27 @@ internal class OptimizeTest {
         val sigma = 1.0
         val generator = NormalDistribution(0.0, sigma)
         val chain = generator.sample(RandomGenerator.default(112667))
-        val x = (1..100).map(Int::toDouble)
+        val x = (1..100).map(Int::toDouble).asBuffer()
 
         val y = x.map {
             it.pow(2) + it + 1 + chain.next()
         }
 
-        val yErr = List(x.size) { sigma }
+        val yErr = DoubleBuffer(x.size) { sigma }
 
-        val chi2 = FunctionOptimization.chiSquared(x, y, yErr) { x1 ->
+        val chi2 = DSProcessor.chiSquaredExpression(
+            x, y, yErr
+        ) { arg ->
             val cWithDefault = bindSymbolOrNull(c) ?: one
-            bindSymbol(a) * x1.pow(2) + bindSymbol(b) * x1 + cWithDefault
+            bindSymbol(a) * arg.pow(2) + bindSymbol(b) * arg + cWithDefault
         }
 
-        val result = chi2.minimize(a to 1.5, b to 0.9, c to 1.0)
+        val result: FunctionOptimization<Double> = chi2.optimizeWith(
+            CMOptimizer,
+            mapOf(a to 1.5, b to 0.9, c to 1.0),
+            FunctionOptimizationTarget.MINIMIZE
+        )
         println(result)
-        println("Chi2/dof = ${result.value / (x.size - 3)}")
+        println("Chi2/dof = ${result.resultValue / (x.size - 3)}")
     }
 }

kmath-core/src/commonMain/kotlin/space/kscience/kmath/data/ColumnarData.kt

@@ -25,6 +25,9 @@ public interface ColumnarData<out T> {
     public operator fun get(symbol: Symbol): Buffer<T>?
 }
 
+@UnstableKMathAPI
+public val ColumnarData<*>.indices: IntRange get() = 0 until size
+
 /**
  * A zero-copy method to represent a [Structure2D] as a two-column x-y data.
  * There could more than two columns in the structure.

kmath-core/src/commonMain/kotlin/space/kscience/kmath/data/XYColumnarData.kt

@@ -32,20 +32,43 @@ public interface XYColumnarData<out T, out X : T, out Y : T> : ColumnarData<T> {
         Symbol.y -> y
         else -> null
     }
-}
 
-@Suppress("FunctionName")
+    public companion object{
-@UnstableKMathAPI
+        @UnstableKMathAPI
-public fun <T, X : T, Y : T> XYColumnarData(x: Buffer<X>, y: Buffer<Y>): XYColumnarData<T, X, Y> {
+        public fun <T, X : T, Y : T> of(x: Buffer<X>, y: Buffer<Y>): XYColumnarData<T, X, Y> {
-    require(x.size == y.size) { "Buffer size mismatch. x buffer size is ${x.size}, y buffer size is ${y.size}" }
+            require(x.size == y.size) { "Buffer size mismatch. x buffer size is ${x.size}, y buffer size is ${y.size}" }
-    return object : XYColumnarData<T, X, Y> {
+            return object : XYColumnarData<T, X, Y> {
-        override val size: Int = x.size
+                override val size: Int = x.size
-        override val x: Buffer<X> = x
+                override val x: Buffer<X> = x
-        override val y: Buffer<Y> = y
+                override val y: Buffer<Y> = y
+            }
+        }
     }
 }
 
 
+/**
+ * Represent a [ColumnarData] as an [XYColumnarData]. The presence or respective columns is checked on creation.
+ */
+@UnstableKMathAPI
+public fun <T> ColumnarData<T>.asXYData(
+    xSymbol: Symbol,
+    ySymbol: Symbol,
+): XYColumnarData<T, T, T> = object : XYColumnarData<T, T, T> {
+    init {
+        requireNotNull(this@asXYData[xSymbol]){"The column with name $xSymbol is not present in $this"}
+        requireNotNull(this@asXYData[ySymbol]){"The column with name $ySymbol is not present in $this"}
+    }
+    override val size: Int get() = this@asXYData.size
+    override val x: Buffer<T> get() = this@asXYData[xSymbol]!!
+    override val y: Buffer<T> get() = this@asXYData[ySymbol]!!
+    override fun get(symbol: Symbol): Buffer<T>? = when (symbol) {
+        Symbol.x -> x
+        Symbol.y -> y
+        else -> this@asXYData.get(symbol)
+    }
+}
+
 /**
  * A zero-copy method to represent a [Structure2D] as a two-column x-y data.
  * There could more than two columns in the structure.

kmath-core/src/commonMain/kotlin/space/kscience/kmath/data/XYErrorColumnarData.kt

@@ -0,0 +1,44 @@
+/*
+ * Copyright 2018-2021 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.data
+
+import space.kscience.kmath.expressions.Symbol
+import space.kscience.kmath.misc.UnstableKMathAPI
+import space.kscience.kmath.structures.Buffer
+
+
+/**
+ * A [ColumnarData] with additional [Symbol.yError] column for an [Symbol.y] error
+ * Inherits [XYColumnarData].
+ */
+@UnstableKMathAPI
+public interface XYErrorColumnarData<T, out X : T, out Y : T> : XYColumnarData<T, X, Y> {
+    public val yErr: Buffer<Y>
+
+    override fun get(symbol: Symbol): Buffer<T> = when (symbol) {
+        Symbol.x -> x
+        Symbol.y -> y
+        Symbol.yError -> yErr
+        else -> error("A column for symbol $symbol not found")
+    }
+
+    public companion object {
+        public fun <T, X : T, Y : T> of(
+            x: Buffer<X>, y: Buffer<Y>, yErr: Buffer<Y>
+        ): XYErrorColumnarData<T, X, Y> {
+            require(x.size == y.size) { "Buffer size mismatch. x buffer size is ${x.size}, y buffer size is ${y.size}" }
+            require(y.size == yErr.size) { "Buffer size mismatch. y buffer size is ${x.size}, yErr buffer size is ${y.size}" }
+
+            return object : XYErrorColumnarData<T, X, Y> {
+                override val size: Int = x.size
+                override val x: Buffer<X> = x
+                override val y: Buffer<Y> = y
+                override val yErr: Buffer<Y> = yErr
+            }
+        }
+    }
+}
+

kmath-core/src/commonMain/kotlin/space/kscience/kmath/data/XYZColumnarData.kt

@@ -10,7 +10,7 @@ import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.structures.Buffer
 
 /**
- * A [XYColumnarData] with guaranteed [x], [y] and [z] columns designated by corresponding symbols.
+ * A [ColumnarData] with guaranteed [x], [y] and [z] columns designated by corresponding symbols.
  * Inherits [XYColumnarData].
  */
 @UnstableKMathAPI

kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/DifferentiableExpression.kt

@@ -5,6 +5,8 @@
 
 package space.kscience.kmath.expressions
 
+import space.kscience.kmath.operations.Algebra
+
 /**
  * Represents expression, which structure can be differentiated.
  *
@@ -64,7 +66,10 @@ public abstract class FirstDerivativeExpression<T> : DifferentiableExpression<T>
 
 /**
  * A factory that converts an expression in autodiff variables to a [DifferentiableExpression]
+ * @param T type of the constants for the expression
+ * @param I type of the actual expression state
+ * @param A type of expression algebra
  */
-public fun interface AutoDiffProcessor<T : Any, I : Any, out A : ExpressionAlgebra<T, I>, out R : Expression<T>> {
+public fun interface AutoDiffProcessor<T, I, out A : Algebra<I>> {
-    public fun process(function: A.() -> I): DifferentiableExpression<T>
+    public fun differentiate(function: A.() -> I): DifferentiableExpression<T>
 }

kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/Expression.kt

@@ -68,6 +68,7 @@ public interface ExpressionAlgebra<in T, E> : Algebra<E> {
 /**
  * Bind a symbol by name inside the [ExpressionAlgebra]
  */
-public fun <T, E> ExpressionAlgebra<T, E>.binding(): ReadOnlyProperty<Any?, E> = ReadOnlyProperty { _, property ->
+public val <T, E> ExpressionAlgebra<T, E>.binding: ReadOnlyProperty<Any?, E>
-    bindSymbol(property.name) ?: error("A variable with name ${property.name} does not exist")
+    get() = ReadOnlyProperty { _, property ->
-}
+        bindSymbol(property.name) ?: error("A variable with name ${property.name} does not exist")
+    }

kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/SimpleAutoDiff.kt

@@ -251,7 +251,9 @@ public class SimpleAutoDiffExpression<T : Any, F : Field<T>>(
 /**
  * Generate [AutoDiffProcessor] for [SimpleAutoDiffExpression]
  */
-public fun <T : Any, F : Field<T>> simpleAutoDiff(field: F): AutoDiffProcessor<T, AutoDiffValue<T>, SimpleAutoDiffField<T, F>, Expression<T>> =
+public fun <T : Any, F : Field<T>> simpleAutoDiff(
+    field: F
+): AutoDiffProcessor<T, AutoDiffValue<T>, SimpleAutoDiffField<T, F>> =
     AutoDiffProcessor { function ->
         SimpleAutoDiffExpression(field, function)
     }

kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/Symbol.kt

@@ -19,9 +19,12 @@ public interface Symbol : MST {
     public val identity: String
 
     public companion object {
-        public val x: StringSymbol = StringSymbol("x")
+        public val x: Symbol = Symbol("x")
-        public val y: StringSymbol = StringSymbol("y")
+        public val xError: Symbol = Symbol("x.error")
-        public val z: StringSymbol = StringSymbol("z")
+        public val y: Symbol = Symbol("y")
+        public val yError: Symbol = Symbol("y.error")
+        public val z: Symbol = Symbol("z")
+        public val zError: Symbol = Symbol("z.error")
     }
 }
 
@@ -29,10 +32,15 @@ public interface Symbol : MST {
  * A [Symbol] with a [String] identity
  */
 @JvmInline
-public value class StringSymbol(override val identity: String) : Symbol {
+internal value class StringSymbol(override val identity: String) : Symbol {
     override fun toString(): String = identity
 }
 
+/**
+ * Create s Symbols with a string identity
+ */
+public fun Symbol(identity: String): Symbol = StringSymbol(identity)
+
 /**
  * A delegate to create a symbol with a string identity in this scope
  */

kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/SymbolIndexer.kt

@@ -9,6 +9,7 @@ import space.kscience.kmath.linear.Point
 import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.nd.Structure2D
 import space.kscience.kmath.structures.BufferFactory
+import space.kscience.kmath.structures.DoubleBuffer
 import kotlin.contracts.InvocationKind
 import kotlin.contracts.contract
 import kotlin.jvm.JvmInline
@@ -47,6 +48,11 @@ public interface SymbolIndexer {
         return symbols.indices.associate { symbols[it] to get(it) }
     }
 
+    public fun <T> Point<T>.toMap(): Map<Symbol, T> {
+        require(size == symbols.size) { "The input array size for indexer should be ${symbols.size} but $size found" }
+        return symbols.indices.associate { symbols[it] to get(it) }
+    }
+
     public operator fun <T> Structure2D<T>.get(rowSymbol: Symbol, columnSymbol: Symbol): T =
         get(indexOf(rowSymbol), indexOf(columnSymbol))
 
@@ -56,6 +62,10 @@ public interface SymbolIndexer {
     public fun <T> Map<Symbol, T>.toPoint(bufferFactory: BufferFactory<T>): Point<T> =
         bufferFactory(symbols.size) { getValue(symbols[it]) }
 
+    public fun Map<Symbol, Double>.toPoint(): DoubleBuffer =
+        DoubleBuffer(symbols.size) { getValue(symbols[it]) }
+
+
     public fun Map<Symbol, Double>.toDoubleArray(): DoubleArray = DoubleArray(symbols.size) { getValue(symbols[it]) }
 }
 

kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/specialExpressions.kt

@@ -0,0 +1,53 @@
+/*
+ * Copyright 2018-2021 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.expressions
+
+import space.kscience.kmath.operations.ExtendedField
+import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.structures.asIterable
+import space.kscience.kmath.structures.indices
+import kotlin.jvm.JvmName
+
+/**
+ * Generate a chi squared expression from given x-y-sigma data and inline model. Provides automatic
+ * differentiation.
+ *
+ * **WARNING** All elements of [yErr] must be positive.
+ */
+@JvmName("genericChiSquaredExpression")
+public fun <T : Comparable<T>, I : Any, A> AutoDiffProcessor<T, I, A>.chiSquaredExpression(
+    x: Buffer<T>,
+    y: Buffer<T>,
+    yErr: Buffer<T>,
+    model: A.(I) -> I,
+): DifferentiableExpression<T> where A : ExtendedField<I>, A : ExpressionAlgebra<T, I> {
+    require(x.size == y.size) { "X and y buffers should be of the same size" }
+    require(y.size == yErr.size) { "Y and yErr buffer should of the same size" }
+
+    return differentiate {
+        var sum = zero
+
+        x.indices.forEach {
+            val xValue = const(x[it])
+            val yValue = const(y[it])
+            val yErrValue = const(yErr[it])
+            val modelValue = model(xValue)
+            sum += ((yValue - modelValue) / yErrValue).pow(2)
+        }
+
+        sum
+    }
+}
+
+public fun <I : Any, A> AutoDiffProcessor<Double, I, A>.chiSquaredExpression(
+    x: Buffer<Double>,
+    y: Buffer<Double>,
+    yErr: Buffer<Double>,
+    model: A.(I) -> I,
+): DifferentiableExpression<Double> where A : ExtendedField<I>, A : ExpressionAlgebra<Double, I> {
+    require(yErr.asIterable().all { it > 0.0 }) { "All errors must be strictly positive" }
+    return chiSquaredExpression<Double, I, A>(x, y, yErr, model)
+}

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

@@ -5,8 +5,6 @@
 
 package space.kscience.kmath.linear
 
-import space.kscience.kmath.nd.as1D
-
 /**
  * A group of methods to solve for *X* in equation *X = A<sup>&minus;1</sup> &middot; B*, where *A* and *B* are
  * matrices or vectors.
@@ -30,20 +28,3 @@ public interface LinearSolver<T : Any> {
     public fun inverse(matrix: Matrix<T>): Matrix<T>
 }
 
-/**
- * Convert matrix to vector if it is possible.
- */
-public fun <T : Any> Matrix<T>.asVector(): Point<T> =
-    if (this.colNum == 1)
-        as1D()
-    else
-        error("Can't convert matrix with more than one column to vector")
-
-/**
- * Creates an n &times; 1 [VirtualMatrix], where n is the size of the given buffer.
- *
- * @param T the type of elements contained in the buffer.
- * @receiver a buffer.
- * @return the new matrix.
- */
-public fun <T : Any> Point<T>.asMatrix(): VirtualMatrix<T> = VirtualMatrix(size, 1) { i, _ -> get(i) }

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

@@ -164,7 +164,7 @@ public interface LinearSpace<T : Any, out A : Ring<T>> {
     public operator fun T.times(v: Point<T>): Point<T> = v * this
 
     /**
-     * Get a feature of the structure in this scope. Structure features take precedence other context features.
+     * Compute a feature of the structure in this scope. Structure features take precedence other context features.
      *
      * @param F the type of feature.
      * @param structure the structure.
@@ -172,7 +172,7 @@ public interface LinearSpace<T : Any, out A : Ring<T>> {
      * @return a feature object or `null` if it isn't present.
      */
     @UnstableKMathAPI
-    public fun <F : StructureFeature> getFeature(structure: Matrix<T>, type: KClass<out F>): F? = structure.getFeature(type)
+    public fun <F : StructureFeature> computeFeature(structure: Matrix<T>, type: KClass<out F>): F? = structure.getFeature(type)
 
     public companion object {
 
@@ -184,7 +184,7 @@ public interface LinearSpace<T : Any, out A : Ring<T>> {
             bufferFactory: BufferFactory<T> = Buffer.Companion::boxing,
         ): LinearSpace<T, A> = BufferedLinearSpace(algebra, bufferFactory)
 
-        public val real: LinearSpace<Double, DoubleField> = buffered(DoubleField, ::DoubleBuffer)
+        public val double: LinearSpace<Double, DoubleField> = buffered(DoubleField, ::DoubleBuffer)
 
         /**
          * Automatic buffered matrix, unboxed if it is possible
@@ -202,9 +202,27 @@ public interface LinearSpace<T : Any, out A : Ring<T>> {
  * @return a feature object or `null` if it isn't present.
  */
 @UnstableKMathAPI
-public inline fun <T : Any, reified F : StructureFeature> LinearSpace<T, *>.getFeature(structure: Matrix<T>): F? =
+public inline fun <T : Any, reified F : StructureFeature> LinearSpace<T, *>.computeFeature(structure: Matrix<T>): F? =
-    getFeature(structure, F::class)
+    computeFeature(structure, F::class)
 
 
-public operator fun <LS : LinearSpace<*, *>, R> LS.invoke(block: LS.() -> R): R = run(block)
+public inline operator fun <LS : LinearSpace<*, *>, R> LS.invoke(block: LS.() -> R): R = run(block)
 
+
+/**
+ * Convert matrix to vector if it is possible.
+ */
+public fun <T : Any> Matrix<T>.asVector(): Point<T> =
+    if (this.colNum == 1)
+        as1D()
+    else
+        error("Can't convert matrix with more than one column to vector")
+
+/**
+ * Creates an n &times; 1 [VirtualMatrix], where n is the size of the given buffer.
+ *
+ * @param T the type of elements contained in the buffer.
+ * @receiver a buffer.
+ * @return the new matrix.
+ */
+public fun <T : Any> Point<T>.asMatrix(): VirtualMatrix<T> = VirtualMatrix(size, 1) { i, _ -> get(i) }

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

@@ -5,6 +5,7 @@
 
 package space.kscience.kmath.linear
 
+import space.kscience.kmath.misc.PerformancePitfall
 import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.nd.getFeature
 import space.kscience.kmath.operations.*
@@ -34,7 +35,7 @@ public class LupDecomposition<T : Any>(
             j == i -> elementContext.one
             else -> elementContext.zero
         }
-    } + LFeature
+    }.withFeature(LFeature)
 
 
     /**
@@ -44,7 +45,7 @@ public class LupDecomposition<T : Any>(
      */
     override val u: Matrix<T> = VirtualMatrix(lu.shape[0], lu.shape[1]) { i, j ->
         if (j >= i) lu[i, j] else elementContext.zero
-    } + UFeature
+    }.withFeature(UFeature)
 
     /**
      * Returns the P rows permutation matrix.
@@ -82,7 +83,7 @@ public fun <T : Comparable<T>> LinearSpace<T, Field<T>>.lup(
     val m = matrix.colNum
     val pivot = IntArray(matrix.rowNum)
 
-    //TODO just waits for KEEP-176
+    //TODO just waits for multi-receivers
     BufferAccessor2D(matrix.rowNum, matrix.colNum, factory).run {
         elementAlgebra {
             val lu = create(matrix)
@@ -156,10 +157,13 @@ public inline fun <reified T : Comparable<T>> LinearSpace<T, Field<T>>.lup(
     noinline checkSingular: (T) -> Boolean,
 ): LupDecomposition<T> = lup(MutableBuffer.Companion::auto, matrix, checkSingular)
 
-public fun LinearSpace<Double, DoubleField>.lup(matrix: Matrix<Double>): LupDecomposition<Double> =
+public fun LinearSpace<Double, DoubleField>.lup(
-    lup(::DoubleBuffer, matrix) { it < 1e-11 }
+    matrix: Matrix<Double>,
+    singularityThreshold: Double = 1e-11,
+): LupDecomposition<Double> =
+    lup(::DoubleBuffer, matrix) { it < singularityThreshold }
 
-public fun <T : Any> LupDecomposition<T>.solveWithLup(
+internal fun <T : Any> LupDecomposition<T>.solve(
     factory: MutableBufferFactory<T>,
     matrix: Matrix<T>,
 ): Matrix<T> {
@@ -207,41 +211,24 @@ public fun <T : Any> LupDecomposition<T>.solveWithLup(
     }
 }
 
-public inline fun <reified T : Any> LupDecomposition<T>.solveWithLup(matrix: Matrix<T>): Matrix<T> =
-    solveWithLup(MutableBuffer.Companion::auto, matrix)
-
 /**
- * Solves a system of linear equations *ax = b** using LUP decomposition.
+ * Produce a generic solver based on LUP decomposition
  */
+@PerformancePitfall()
 @OptIn(UnstableKMathAPI::class)
-public inline fun <reified T : Comparable<T>> LinearSpace<T, Field<T>>.solveWithLup(
+public fun <T : Comparable<T>, F : Field<T>> LinearSpace<T, F>.lupSolver(
-    a: Matrix<T>,
+    bufferFactory: MutableBufferFactory<T>,
-    b: Matrix<T>,
+    singularityCheck: (T) -> Boolean,
-    noinline bufferFactory: MutableBufferFactory<T> = MutableBuffer.Companion::auto,
+): LinearSolver<T> = object : LinearSolver<T> {
-    noinline checkSingular: (T) -> Boolean,
+    override fun solve(a: Matrix<T>, b: Matrix<T>): Matrix<T> {
-): Matrix<T> {
+        // Use existing decomposition if it is provided by matrix
-    // Use existing decomposition if it is provided by matrix
+        val decomposition = a.getFeature() ?: lup(bufferFactory, a, singularityCheck)
-    val decomposition = a.getFeature() ?: lup(bufferFactory, a, checkSingular)
+        return decomposition.solve(bufferFactory, b)
-    return decomposition.solveWithLup(bufferFactory, b)
+    }
+
+    override fun inverse(matrix: Matrix<T>): Matrix<T> = solve(matrix, one(matrix.rowNum, matrix.colNum))
 }
 
-public inline fun <reified T : Comparable<T>> LinearSpace<T, Field<T>>.inverseWithLup(
+@PerformancePitfall
-    matrix: Matrix<T>,
+public fun LinearSpace<Double, DoubleField>.lupSolver(singularityThreshold: Double = 1e-11): LinearSolver<Double> =
-    noinline bufferFactory: MutableBufferFactory<T> = MutableBuffer.Companion::auto,
+    lupSolver(::DoubleBuffer) { it < singularityThreshold }
-    noinline checkSingular: (T) -> Boolean,
-): Matrix<T> = solveWithLup(matrix, one(matrix.rowNum, matrix.colNum), bufferFactory, checkSingular)
-
-
-@OptIn(UnstableKMathAPI::class)
-public fun LinearSpace<Double, DoubleField>.solveWithLup(a: Matrix<Double>, b: Matrix<Double>): Matrix<Double> {
-    // Use existing decomposition if it is provided by matrix
-    val bufferFactory: MutableBufferFactory<Double> = ::DoubleBuffer
-    val decomposition: LupDecomposition<Double> = a.getFeature() ?: lup(bufferFactory, a) { it < 1e-11 }
-    return decomposition.solveWithLup(bufferFactory, b)
-}
-
-/**
- * Inverses a square matrix using LUP decomposition. Non square matrix will throw an error.
- */
-public fun LinearSpace<Double, DoubleField>.inverseWithLup(matrix: Matrix<Double>): Matrix<Double> =
-    solveWithLup(matrix, one(matrix.rowNum, matrix.colNum))

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

@@ -7,6 +7,8 @@ package space.kscience.kmath.linear
 
 import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.operations.Ring
+import space.kscience.kmath.structures.BufferAccessor2D
+import space.kscience.kmath.structures.MutableBuffer
 
 public class MatrixBuilder<T : Any, out A : Ring<T>>(
     public val linearSpace: LinearSpace<T, A>,
@@ -45,4 +47,31 @@ public inline fun <T : Any> LinearSpace<T, Ring<T>>.column(
     crossinline builder: (Int) -> T,
 ): Matrix<T> = buildMatrix(size, 1) { i, _ -> builder(i) }
 
-public fun <T : Any> LinearSpace<T, Ring<T>>.column(vararg values: T): Matrix<T> = column(values.size, values::get)
+public fun <T : Any> LinearSpace<T, Ring<T>>.column(vararg values: T): Matrix<T> = column(values.size, values::get)
+
+public object SymmetricMatrixFeature : MatrixFeature
+
+/**
+ * Naive implementation of a symmetric matrix builder, that adds a [SymmetricMatrixFeature] tag. The resulting matrix contains
+ * full `size^2` number of elements, but caches elements during calls to save [builder] calls. [builder] is always called in the
+ * upper triangle region meaning that `i <= j`
+ */
+public fun <T : Any, A : Ring<T>> MatrixBuilder<T, A>.symmetric(
+    builder: (i: Int, j: Int) -> T,
+): Matrix<T> {
+    require(columns == rows) { "In order to build symmetric matrix, number of rows $rows should be equal to number of columns $columns" }
+    return with(BufferAccessor2D<T?>(rows, rows, MutableBuffer.Companion::boxing)) {
+        val cache = factory(rows * rows) { null }
+        linearSpace.buildMatrix(rows, rows) { i, j ->
+            val cached = cache[i, j]
+            if (cached == null) {
+                val value = if (i <= j) builder(i, j) else builder(j, i)
+                cache[i, j] = value
+                cache[j, i] = value
+                value
+            } else {
+                cached
+            }
+        }.withFeature(SymmetricMatrixFeature)
+    }
+}

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

@@ -5,6 +5,7 @@
 
 package space.kscience.kmath.linear
 
+import space.kscience.kmath.misc.FeatureSet
 import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.nd.StructureFeature
 import space.kscience.kmath.nd.getFeature
@@ -18,7 +19,7 @@ import kotlin.reflect.KClass
  */
 public class MatrixWrapper<out T : Any> internal constructor(
     public val origin: Matrix<T>,
-    public val features: Set<MatrixFeature>,
+    public val features: FeatureSet<StructureFeature>,
 ) : Matrix<T> by origin {
 
     /**
@@ -28,8 +29,7 @@ public class MatrixWrapper<out T : Any> internal constructor(
     @UnstableKMathAPI
     @Suppress("UNCHECKED_CAST")
     override fun <F : StructureFeature> getFeature(type: KClass<out F>): F? =
-        features.singleOrNull(type::isInstance) as? F
+        features.getFeature(type) ?: origin.getFeature(type)
-            ?: origin.getFeature(type)
 
     override fun toString(): String = "MatrixWrapper(matrix=$origin, features=$features)"
 }
@@ -45,20 +45,23 @@ public val <T : Any> Matrix<T>.origin: Matrix<T>
 /**
  * Add a single feature to a [Matrix]
  */
-public operator fun <T : Any> Matrix<T>.plus(newFeature: MatrixFeature): MatrixWrapper<T> = if (this is MatrixWrapper) {
+public fun <T : Any> Matrix<T>.withFeature(newFeature: MatrixFeature): MatrixWrapper<T> = if (this is MatrixWrapper) {
-    MatrixWrapper(origin, features + newFeature)
+    MatrixWrapper(origin, features.with(newFeature))
 } else {
-    MatrixWrapper(this, setOf(newFeature))
+    MatrixWrapper(this, FeatureSet.of(newFeature))
 }
 
+@Deprecated("To be replaced by withFeature")
+public operator fun <T : Any> Matrix<T>.plus(newFeature: MatrixFeature): MatrixWrapper<T> = withFeature(newFeature)
+
 /**
  * Add a collection of features to a [Matrix]
  */
-public operator fun <T : Any> Matrix<T>.plus(newFeatures: Collection<MatrixFeature>): MatrixWrapper<T> =
+public fun <T : Any> Matrix<T>.withFeatures(newFeatures: Iterable<MatrixFeature>): MatrixWrapper<T> =
     if (this is MatrixWrapper) {
-        MatrixWrapper(origin, features + newFeatures)
+        MatrixWrapper(origin, features.with(newFeatures))
     } else {
-        MatrixWrapper(this, newFeatures.toSet())
+        MatrixWrapper(this, FeatureSet.of(newFeatures))
     }
 
 /**
@@ -69,7 +72,7 @@ public fun <T : Any> LinearSpace<T, Ring<T>>.one(
     columns: Int,
 ): Matrix<T> = VirtualMatrix(rows, columns) { i, j ->
     if (i == j) elementAlgebra.one else elementAlgebra.zero
-} + UnitFeature
+}.withFeature(UnitFeature)
 
 
 /**
@@ -80,7 +83,7 @@ public fun <T : Any> LinearSpace<T, Ring<T>>.zero(
     columns: Int,
 ): Matrix<T> = VirtualMatrix(rows, columns) { _, _ ->
     elementAlgebra.zero
-} + ZeroFeature
+}.withFeature(ZeroFeature)
 
 public class TransposedFeature<out T : Any>(public val original: Matrix<T>) : MatrixFeature
 
@@ -88,7 +91,7 @@ public class TransposedFeature<out T : Any>(public val original: Matrix<T>) : Ma
  * Create a virtual transposed matrix without copying anything. `A.transpose().transpose() === A`
  */
 @OptIn(UnstableKMathAPI::class)
-public fun <T : Any> Matrix<T>.transpose(): Matrix<T> = getFeature<TransposedFeature<T>>()?.original ?: (VirtualMatrix(
+public fun <T : Any> Matrix<T>.transpose(): Matrix<T> = getFeature<TransposedFeature<T>>()?.original ?: VirtualMatrix(
     colNum,
     rowNum,
-) { i, j -> get(j, i) } + TransposedFeature(this))
+) { i, j -> get(j, i) }.withFeature(TransposedFeature(this))

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

@@ -20,3 +20,6 @@ public class VirtualMatrix<out T : Any>(
 
     override operator fun get(i: Int, j: Int): T = generator(i, j)
 }
+
+public fun <T : Any> MatrixBuilder<T, *>.virtual(generator: (i: Int, j: Int) -> T): VirtualMatrix<T> =
+    VirtualMatrix(rows, columns, generator)

kmath-core/src/commonMain/kotlin/space/kscience/kmath/misc/Featured.kt

@@ -0,0 +1,61 @@
+/*
+ * Copyright 2018-2021 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.misc
+
+import kotlin.jvm.JvmInline
+import kotlin.reflect.KClass
+
+/**
+ * A entity that contains a set of features defined by their types
+ */
+public interface Featured<F : Any> {
+    public fun <T : F> getFeature(type: FeatureKey<T>): T?
+}
+
+public typealias FeatureKey<T> = KClass<out T>
+
+public interface Feature<F : Feature<F>> {
+
+    /**
+     * A key used for extraction
+     */
+    @Suppress("UNCHECKED_CAST")
+    public val key: FeatureKey<F>
+        get() = this::class as FeatureKey<F>
+}
+
+/**
+ * A container for a set of features
+ */
+@JvmInline
+public value class FeatureSet<F : Feature<F>> private constructor(public val features: Map<FeatureKey<F>, F>) : Featured<F> {
+    @Suppress("UNCHECKED_CAST")
+    override fun <T : F> getFeature(type: FeatureKey<T>): T? = features[type]?.let { it as T }
+
+    public inline fun <reified T : F> getFeature(): T? = getFeature(T::class)
+
+    public fun <T : F> with(feature: T, type: FeatureKey<F> = feature.key): FeatureSet<F> =
+        FeatureSet(features + (type to feature))
+
+    public fun with(other: FeatureSet<F>): FeatureSet<F> = FeatureSet(features + other.features)
+
+    public fun with(vararg otherFeatures: F): FeatureSet<F> =
+        FeatureSet(features + otherFeatures.associateBy { it.key })
+
+    public fun with(otherFeatures: Iterable<F>): FeatureSet<F> =
+        FeatureSet(features + otherFeatures.associateBy { it.key })
+
+    public operator fun iterator(): Iterator<F> = features.values.iterator()
+
+    override fun toString(): String = features.values.joinToString(prefix = "[ ", postfix = " ]")
+
+
+    public companion object {
+        public fun <F : Feature<F>> of(vararg features: F): FeatureSet<F> = FeatureSet(features.associateBy { it.key })
+        public fun <F : Feature<F>> of(features: Iterable<F>): FeatureSet<F> =
+            FeatureSet(features.associateBy { it.key })
+    }
+}

kmath-core/src/commonMain/kotlin/space/kscience/kmath/misc/annotations.kt

@@ -12,7 +12,7 @@ package space.kscience.kmath.misc
  * in some way that may break some code.
  */
 @MustBeDocumented
-@Retention(value = AnnotationRetention.BINARY)
+@Retention(value = AnnotationRetention.SOURCE)
 @RequiresOptIn("This API is unstable and could change in future", RequiresOptIn.Level.WARNING)
 public annotation class UnstableKMathAPI
 
@@ -21,7 +21,7 @@ public annotation class UnstableKMathAPI
  * slow-down in some cases. Refer to the documentation and benchmark it to be sure.
  */
 @MustBeDocumented
-@Retention(value = AnnotationRetention.BINARY)
+@Retention(value = AnnotationRetention.SOURCE)
 @RequiresOptIn(
     "Refer to the documentation to use this API in performance-critical code",
     RequiresOptIn.Level.WARNING

kmath-core/src/commonMain/kotlin/space/kscience/kmath/misc/logging.kt

@@ -0,0 +1,22 @@
+/*
+ * Copyright 2018-2021 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.misc
+
+import space.kscience.kmath.misc.Loggable.Companion.INFO
+
+public fun interface Loggable {
+    public fun log(tag: String, block: () -> String)
+
+    public companion object {
+        public const val INFO: String = "INFO"
+
+        public val console: Loggable = Loggable { tag, block ->
+            println("[$tag] ${block()}")
+        }
+    }
+}
+
+public fun Loggable.log(block: () -> String): Unit = log(INFO, block)

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

@@ -6,6 +6,8 @@
 package space.kscience.kmath.nd
 
 import space.kscience.kmath.linear.LinearSpace
+import space.kscience.kmath.misc.Feature
+import space.kscience.kmath.misc.Featured
 import space.kscience.kmath.misc.PerformancePitfall
 import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.operations.Ring
@@ -16,7 +18,7 @@ import kotlin.jvm.JvmName
 import kotlin.native.concurrent.ThreadLocal
 import kotlin.reflect.KClass
 
-public interface StructureFeature
+public interface StructureFeature : Feature<StructureFeature>
 
 /**
  * Represents n-dimensional structure i.e., multidimensional container of items of the same type and size. The number
@@ -27,7 +29,7 @@ public interface StructureFeature
  *
  * @param T the type of items.
  */
-public interface StructureND<out T> {
+public interface StructureND<out T> : Featured<StructureFeature> {
     /**
      * The shape of structure i.e., non-empty sequence of non-negative integers that specify sizes of dimensions of
      * this structure.
@@ -60,7 +62,7 @@ public interface StructureND<out T> {
      * If the feature is not present, `null` is returned.
      */
     @UnstableKMathAPI
-    public fun <F : StructureFeature> getFeature(type: KClass<out F>): F? = null
+    override fun <F : StructureFeature> getFeature(type: KClass<out F>): F? = null
 
     public companion object {
         /**

kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/BufferAlgebra.kt

@@ -0,0 +1,146 @@
+/*
+ * Copyright 2018-2021 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.operations
+
+import space.kscience.kmath.misc.UnstableKMathAPI
+import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.structures.BufferFactory
+import space.kscience.kmath.structures.DoubleBuffer
+
+/**
+ * An algebra over [Buffer]
+ */
+@UnstableKMathAPI
+public interface BufferAlgebra<T, A : Algebra<T>> : Algebra<Buffer<T>> {
+    public val bufferFactory: BufferFactory<T>
+    public val elementAlgebra: A
+
+    //TODO move to multi-receiver inline extension
+    public fun Buffer<T>.map(block: (T) -> T): Buffer<T> = bufferFactory(size) { block(get(it)) }
+
+    public fun Buffer<T>.zip(other: Buffer<T>, block: (left: T, right: T) -> T): Buffer<T> {
+        require(size == other.size) { "Incompatible buffer sizes. left: $size, right: ${other.size}" }
+        return bufferFactory(size) { block(this[it], other[it]) }
+    }
+
+    override fun unaryOperationFunction(operation: String): (arg: Buffer<T>) -> Buffer<T> {
+        val operationFunction = elementAlgebra.unaryOperationFunction(operation)
+        return { arg -> bufferFactory(arg.size) { operationFunction(arg[it]) } }
+    }
+
+    override fun binaryOperationFunction(operation: String): (left: Buffer<T>, right: Buffer<T>) -> Buffer<T> {
+        val operationFunction = elementAlgebra.binaryOperationFunction(operation)
+        return { left, right ->
+            bufferFactory(left.size) { operationFunction(left[it], right[it]) }
+        }
+    }
+}
+
+@UnstableKMathAPI
+public fun <T, A : TrigonometricOperations<T>> BufferAlgebra<T, A>.sin(arg: Buffer<T>): Buffer<T> =
+    arg.map(elementAlgebra::sin)
+
+@UnstableKMathAPI
+public fun <T, A : TrigonometricOperations<T>> BufferAlgebra<T, A>.cos(arg: Buffer<T>): Buffer<T> =
+    arg.map(elementAlgebra::cos)
+
+@UnstableKMathAPI
+public fun <T, A : TrigonometricOperations<T>> BufferAlgebra<T, A>.tan(arg: Buffer<T>): Buffer<T> =
+    arg.map(elementAlgebra::tan)
+
+@UnstableKMathAPI
+public fun <T, A : TrigonometricOperations<T>> BufferAlgebra<T, A>.asin(arg: Buffer<T>): Buffer<T> =
+    arg.map(elementAlgebra::asin)
+
+@UnstableKMathAPI
+public fun <T, A : TrigonometricOperations<T>> BufferAlgebra<T, A>.acos(arg: Buffer<T>): Buffer<T> =
+    arg.map(elementAlgebra::acos)
+
+@UnstableKMathAPI
+public fun <T, A : TrigonometricOperations<T>> BufferAlgebra<T, A>.atan(arg: Buffer<T>): Buffer<T> =
+    arg.map(elementAlgebra::atan)
+
+@UnstableKMathAPI
+public fun <T, A : ExponentialOperations<T>> BufferAlgebra<T, A>.exp(arg: Buffer<T>): Buffer<T> =
+    arg.map(elementAlgebra::exp)
+
+@UnstableKMathAPI
+public fun <T, A : ExponentialOperations<T>> BufferAlgebra<T, A>.ln(arg: Buffer<T>): Buffer<T> =
+    arg.map(elementAlgebra::ln)
+
+@UnstableKMathAPI
+public fun <T, A : ExponentialOperations<T>> BufferAlgebra<T, A>.sinh(arg: Buffer<T>): Buffer<T> =
+    arg.map(elementAlgebra::sinh)
+
+@UnstableKMathAPI
+public fun <T, A : ExponentialOperations<T>> BufferAlgebra<T, A>.cosh(arg: Buffer<T>): Buffer<T> =
+    arg.map(elementAlgebra::cosh)
+
+@UnstableKMathAPI
+public fun <T, A : ExponentialOperations<T>> BufferAlgebra<T, A>.tanh(arg: Buffer<T>): Buffer<T> =
+    arg.map(elementAlgebra::tanh)
+
+@UnstableKMathAPI
+public fun <T, A : ExponentialOperations<T>> BufferAlgebra<T, A>.asinh(arg: Buffer<T>): Buffer<T> =
+    arg.map(elementAlgebra::asinh)
+
+@UnstableKMathAPI
+public fun <T, A : ExponentialOperations<T>> BufferAlgebra<T, A>.acosh(arg: Buffer<T>): Buffer<T> =
+    arg.map(elementAlgebra::acosh)
+
+@UnstableKMathAPI
+public fun <T, A : ExponentialOperations<T>> BufferAlgebra<T, A>.atanh(arg: Buffer<T>): Buffer<T> =
+    arg.map(elementAlgebra::atanh)
+
+@UnstableKMathAPI
+public fun <T, A : PowerOperations<T>> BufferAlgebra<T, A>.pow(arg: Buffer<T>, pow: Number): Buffer<T> =
+    with(elementAlgebra) { arg.map { power(it, pow) } }
+
+
+@UnstableKMathAPI
+public class BufferField<T, A : Field<T>>(
+    override val bufferFactory: BufferFactory<T>,
+    override val elementAlgebra: A,
+    public val size: Int
+) : BufferAlgebra<T, A>, Field<Buffer<T>> {
+
+    public fun produce(vararg elements: T): Buffer<T> {
+        require(elements.size == size) { "Expected $size elements but found ${elements.size}" }
+        return bufferFactory(size) { elements[it] }
+    }
+
+    override val zero: Buffer<T> = bufferFactory(size) { elementAlgebra.zero }
+    override val one: Buffer<T> = bufferFactory(size) { elementAlgebra.one }
+
+
+    override fun add(a: Buffer<T>, b: Buffer<T>): Buffer<T> = a.zip(b, elementAlgebra::add)
+    override fun multiply(a: Buffer<T>, b: Buffer<T>): Buffer<T> = a.zip(b, elementAlgebra::multiply)
+    override fun divide(a: Buffer<T>, b: Buffer<T>): Buffer<T> = a.zip(b, elementAlgebra::divide)
+
+    override fun scale(a: Buffer<T>, value: Double): Buffer<T> = with(elementAlgebra) { a.map { scale(it, value) } }
+    override fun Buffer<T>.unaryMinus(): Buffer<T> = with(elementAlgebra) { map { -it } }
+
+    override fun unaryOperationFunction(operation: String): (arg: Buffer<T>) -> Buffer<T> {
+        return super<BufferAlgebra>.unaryOperationFunction(operation)
+    }
+
+    override fun binaryOperationFunction(operation: String): (left: Buffer<T>, right: Buffer<T>) -> Buffer<T> {
+        return super<BufferAlgebra>.binaryOperationFunction(operation)
+    }
+}
+
+//Double buffer specialization
+
+@UnstableKMathAPI
+public fun BufferField<Double, *>.produce(vararg elements: Number): Buffer<Double> {
+    require(elements.size == size) { "Expected $size elements but found ${elements.size}" }
+    return bufferFactory(size) { elements[it].toDouble() }
+}
+
+@UnstableKMathAPI
+public fun DoubleField.bufferAlgebra(size: Int): BufferField<Double, DoubleField> =
+    BufferField(::DoubleBuffer, DoubleField, size)
+

kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/BufferAccessor2D.kt

@@ -13,7 +13,7 @@ import space.kscience.kmath.nd.as2D
 /**
  * A context that allows to operate on a [MutableBuffer] as on 2d array
  */
-internal class BufferAccessor2D<T : Any>(
+internal class BufferAccessor2D<T>(
     val rowNum: Int,
     val colNum: Int,
     val factory: MutableBufferFactory<T>,

kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/DoubleBufferField.kt

@@ -5,13 +5,15 @@
 
 package space.kscience.kmath.structures
 
-import space.kscience.kmath.operations.ExtendedField
+import space.kscience.kmath.linear.Point
-import space.kscience.kmath.operations.ExtendedFieldOperations
+import space.kscience.kmath.misc.UnstableKMathAPI
+import space.kscience.kmath.operations.*
 import kotlin.math.*
 
 /**
  * [ExtendedFieldOperations] over [DoubleBuffer].
  */
+@Deprecated("To be replaced by generic BufferAlgebra")
 public object DoubleBufferFieldOperations : ExtendedFieldOperations<Buffer<Double>> {
     override fun Buffer<Double>.unaryMinus(): DoubleBuffer = if (this is DoubleBuffer) {
         DoubleBuffer(size) { -array[it] }
@@ -161,12 +163,17 @@ public object DoubleBufferFieldOperations : ExtendedFieldOperations<Buffer<Doubl
         DoubleBuffer(DoubleArray(arg.size) { ln(arg[it]) })
 }
 
+public object DoubleL2Norm : Norm<Point<Double>, Double> {
+    override fun norm(arg: Point<Double>): Double = sqrt(arg.fold(0.0) { acc: Double, d: Double -> acc + d.pow(2) })
+}
+
 /**
  * [ExtendedField] over [DoubleBuffer].
  *
  * @property size the size of buffers to operate on.
  */
-public class DoubleBufferField(public val size: Int) : ExtendedField<Buffer<Double>> {
+@Deprecated("To be replaced by generic BufferAlgebra")
+public class DoubleBufferField(public val size: Int) : ExtendedField<Buffer<Double>>, Norm<Buffer<Double>, Double> {
     override val zero: Buffer<Double> by lazy { DoubleBuffer(size) { 0.0 } }
     override val one: Buffer<Double> by lazy { DoubleBuffer(size) { 1.0 } }
 
@@ -274,4 +281,6 @@ public class DoubleBufferField(public val size: Int) : ExtendedField<Buffer<Doub
         require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
         return DoubleBufferFieldOperations.ln(arg)
     }
+
+    override fun norm(arg: Buffer<Double>): Double  = DoubleL2Norm.norm(arg)
 }

kmath-core/src/commonTest/kotlin/space/kscience/kmath/expressions/ExpressionFieldTest.kt

@@ -16,7 +16,7 @@ class ExpressionFieldTest {
     @Test
     fun testExpression() {
         val expression = with(FunctionalExpressionField(DoubleField)) {
-            val x by binding()
+            val x by binding
             x * x + 2 * x + one
         }
 
@@ -27,7 +27,7 @@ class ExpressionFieldTest {
     @Test
     fun separateContext() {
         fun <T> FunctionalExpressionField<T, *>.expression(): Expression<T> {
-            val x by binding()
+            val x by binding
             return x * x + 2 * x + one
         }
 
@@ -38,7 +38,7 @@ class ExpressionFieldTest {
     @Test
     fun valueExpression() {
         val expressionBuilder: FunctionalExpressionField<Double, *>.() -> Expression<Double> = {
-            val x by binding()
+            val x by binding
             x * x + 2 * x + one
         }
 

kmath-core/src/commonTest/kotlin/space/kscience/kmath/linear/DoubleLUSolverTest.kt

@@ -22,14 +22,14 @@ class DoubleLUSolverTest {
 
     @Test
     fun testInvertOne() {
-        val matrix = LinearSpace.real.one(2, 2)
+        val matrix = LinearSpace.double.one(2, 2)
-        val inverted = LinearSpace.real.inverseWithLup(matrix)
+        val inverted = LinearSpace.double.lupSolver().inverse(matrix)
         assertMatrixEquals(matrix, inverted)
     }
 
     @Test
     fun testDecomposition() {
-        LinearSpace.real.run {
+        LinearSpace.double.run {
             val matrix = matrix(2, 2)(
                 3.0, 1.0,
                 2.0, 3.0
@@ -46,14 +46,14 @@ class DoubleLUSolverTest {
 
     @Test
     fun testInvert() {
-        val matrix = LinearSpace.real.matrix(2, 2)(
+        val matrix = LinearSpace.double.matrix(2, 2)(
             3.0, 1.0,
             1.0, 3.0
         )
 
-        val inverted = LinearSpace.real.inverseWithLup(matrix)
+        val inverted = LinearSpace.double.lupSolver().inverse(matrix)
 
-        val expected = LinearSpace.real.matrix(2, 2)(
+        val expected = LinearSpace.double.matrix(2, 2)(
             0.375, -0.125,
             -0.125, 0.375
         )

kmath-core/src/commonTest/kotlin/space/kscience/kmath/linear/MatrixTest.kt

@@ -17,16 +17,17 @@ import kotlin.test.assertTrue
 @OptIn(PerformancePitfall::class)
 @Suppress("UNUSED_VARIABLE")
 class MatrixTest {
+
     @Test
     fun testTranspose() {
-        val matrix = LinearSpace.real.one(3, 3)
+        val matrix = LinearSpace.double.one(3, 3)
         val transposed = matrix.transpose()
         assertTrue { StructureND.contentEquals(matrix, transposed) }
     }
 
     @Test
     fun testBuilder() {
-        val matrix = LinearSpace.real.matrix(2, 3)(
+        val matrix = LinearSpace.double.matrix(2, 3)(
             1.0, 0.0, 0.0,
             0.0, 1.0, 2.0
         )
@@ -48,7 +49,7 @@ class MatrixTest {
         infix fun Matrix<Double>.pow(power: Int): Matrix<Double> {
             var res = this
             repeat(power - 1) {
-                res = LinearSpace.real.run { res dot this@pow }
+                res = LinearSpace.double.run { res dot this@pow }
             }
             return res
         }
@@ -61,7 +62,7 @@ class MatrixTest {
         val firstMatrix = StructureND.auto(2, 3) { (i, j) -> (i + j).toDouble() }.as2D()
         val secondMatrix = StructureND.auto(3, 2) { (i, j) -> (i + j).toDouble() }.as2D()
 
-        LinearSpace.real.run {
+        LinearSpace.double.run {
 //            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

kmath-core/src/commonTest/kotlin/space/kscience/kmath/structures/NumberNDFieldTest.kt

@@ -40,7 +40,7 @@ class NumberNDFieldTest {
     @Test
     fun testGeneration() {
 
-        val array = LinearSpace.real.buildMatrix(3, 3) { i, j ->
+        val array = LinearSpace.double.buildMatrix(3, 3) { i, j ->
             (i * 10 + j).toDouble()
         }
 

kmath-dimensions/src/commonMain/kotlin/space/kscience/kmath/dimensions/Wrappers.kt

@@ -151,7 +151,7 @@ public value class DMatrixContext<T : Any, out A : Ring<T>>(public val context:
         context.run { (this@transpose as Matrix<T>).transpose() }.coerce()
 
     public companion object {
-        public val real: DMatrixContext<Double, DoubleField> = DMatrixContext(LinearSpace.real)
+        public val real: DMatrixContext<Double, DoubleField> = DMatrixContext(LinearSpace.double)
     }
 }
 

kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/_generated.kt

@@ -1,995 +0,0 @@
-/*
- * Copyright 2018-2021 KMath contributors.
- * Use of this source code is governed by the Apache 2.0 license that can be found in the LICENSE file.
- */
-
-/* This file is generated with buildSrc/src/main/kotlin/space/kscience/kmath/ejml/codegen/ejmlCodegen.kt */
-
-package space.kscience.kmath.ejml
-
-import org.ejml.data.*
-import org.ejml.dense.row.CommonOps_DDRM
-import org.ejml.dense.row.CommonOps_FDRM
-import org.ejml.dense.row.factory.DecompositionFactory_DDRM
-import org.ejml.dense.row.factory.DecompositionFactory_FDRM
-import org.ejml.sparse.FillReducing
-import org.ejml.sparse.csc.CommonOps_DSCC
-import org.ejml.sparse.csc.CommonOps_FSCC
-import org.ejml.sparse.csc.factory.DecompositionFactory_DSCC
-import org.ejml.sparse.csc.factory.DecompositionFactory_FSCC
-import org.ejml.sparse.csc.factory.LinearSolverFactory_DSCC
-import org.ejml.sparse.csc.factory.LinearSolverFactory_FSCC
-import space.kscience.kmath.linear.*
-import space.kscience.kmath.linear.Matrix
-import space.kscience.kmath.misc.UnstableKMathAPI
-import space.kscience.kmath.nd.StructureFeature
-import space.kscience.kmath.operations.DoubleField
-import space.kscience.kmath.operations.FloatField
-import space.kscience.kmath.operations.invoke
-import space.kscience.kmath.structures.DoubleBuffer
-import space.kscience.kmath.structures.FloatBuffer
-import kotlin.reflect.KClass
-import kotlin.reflect.cast
-
-/**
- * [EjmlVector] specialization for [Double].
- */
-public class EjmlDoubleVector<out M : DMatrix>(override val origin: M) : EjmlVector<Double, M>(origin) {
-    init {
-        require(origin.numRows == 1) { "The origin matrix must have only one row to form a vector" }
-    }
-
-    override operator fun get(index: Int): Double = origin[0, index]
-}
-
-/**
- * [EjmlVector] specialization for [Float].
- */
-public class EjmlFloatVector<out M : FMatrix>(override val origin: M) : EjmlVector<Float, M>(origin) {
-    init {
-        require(origin.numRows == 1) { "The origin matrix must have only one row to form a vector" }
-    }
-
-    override operator fun get(index: Int): Float = origin[0, index]
-}
-
-/**
- * [EjmlMatrix] specialization for [Double].
- */
-public class EjmlDoubleMatrix<out M : DMatrix>(override val origin: M) : EjmlMatrix<Double, M>(origin) {
-    override operator fun get(i: Int, j: Int): Double = origin[i, j]
-}
-
-/**
- * [EjmlMatrix] specialization for [Float].
- */
-public class EjmlFloatMatrix<out M : FMatrix>(override val origin: M) : EjmlMatrix<Float, M>(origin) {
-    override operator fun get(i: Int, j: Int): Float = origin[i, j]
-}
-
-/**
- * [EjmlLinearSpace] implementation based on [CommonOps_DDRM], [DecompositionFactory_DDRM] operations and
- * [DMatrixRMaj] matrices.
- */
-public object EjmlLinearSpaceDDRM : EjmlLinearSpace<Double, DoubleField, DMatrixRMaj>() {
-    /**
-     * The [DoubleField] reference.
-     */
-    override val elementAlgebra: DoubleField get() = DoubleField
-
-    @Suppress("UNCHECKED_CAST")
-    override fun Matrix<Double>.toEjml(): EjmlDoubleMatrix<DMatrixRMaj> = when {
-        this is EjmlDoubleMatrix<*> && origin is DMatrixRMaj -> this as EjmlDoubleMatrix<DMatrixRMaj>
-        else -> buildMatrix(rowNum, colNum) { i, j -> get(i, j) }
-    }
-
-    @Suppress("UNCHECKED_CAST")
-    override fun Point<Double>.toEjml(): EjmlDoubleVector<DMatrixRMaj> = when {
-        this is EjmlDoubleVector<*> && origin is DMatrixRMaj -> this as EjmlDoubleVector<DMatrixRMaj>
-        else -> EjmlDoubleVector(DMatrixRMaj(size, 1).also {
-            (0 until it.numRows).forEach { row -> it[row, 0] = get(row) }
-        })
-    }
-
-    override fun buildMatrix(
-        rows: Int,
-        columns: Int,
-        initializer: DoubleField.(i: Int, j: Int) -> Double,
-    ): EjmlDoubleMatrix<DMatrixRMaj> = DMatrixRMaj(rows, columns).also {
-        (0 until rows).forEach { row ->
-            (0 until columns).forEach { col -> it[row, col] = elementAlgebra.initializer(row, col) }
-        }
-    }.wrapMatrix()
-
-    override fun buildVector(
-        size: Int,
-        initializer: DoubleField.(Int) -> Double,
-    ): EjmlDoubleVector<DMatrixRMaj> = EjmlDoubleVector(DMatrixRMaj(size, 1).also {
-        (0 until it.numRows).forEach { row -> it[row, 0] = elementAlgebra.initializer(row) }
-    })
-
-    private fun <T : DMatrix> T.wrapMatrix() = EjmlDoubleMatrix(this)
-    private fun <T : DMatrix> T.wrapVector() = EjmlDoubleVector(this)
-
-    override fun Matrix<Double>.unaryMinus(): Matrix<Double> = this * elementAlgebra { -one }
-
-    override fun Matrix<Double>.dot(other: Matrix<Double>): EjmlDoubleMatrix<DMatrixRMaj> {
-        val out = DMatrixRMaj(1, 1)
-        CommonOps_DDRM.mult(toEjml().origin, other.toEjml().origin, out)
-        return out.wrapMatrix()
-    }
-
-    override fun Matrix<Double>.dot(vector: Point<Double>): EjmlDoubleVector<DMatrixRMaj> {
-        val out = DMatrixRMaj(1, 1)
-        CommonOps_DDRM.mult(toEjml().origin, vector.toEjml().origin, out)
-        return out.wrapVector()
-    }
-
-    override operator fun Matrix<Double>.minus(other: Matrix<Double>): EjmlDoubleMatrix<DMatrixRMaj> {
-        val out = DMatrixRMaj(1, 1)
-
-        CommonOps_DDRM.add(
-            elementAlgebra.one, 
-            toEjml().origin, 
-            elementAlgebra { -one },
-            other.toEjml().origin, 
-            out,
-        )
-
-        return out.wrapMatrix()
-    }
-
-    override operator fun Matrix<Double>.times(value: Double): EjmlDoubleMatrix<DMatrixRMaj> {
-        val res = DMatrixRMaj(1, 1)
-        CommonOps_DDRM.scale(value, toEjml().origin, res)
-        return res.wrapMatrix()
-    }
-
-    override fun Point<Double>.unaryMinus(): EjmlDoubleVector<DMatrixRMaj> {
-        val res = DMatrixRMaj(1, 1)
-        CommonOps_DDRM.changeSign(toEjml().origin, res)
-        return res.wrapVector()
-    }
-
-    override fun Matrix<Double>.plus(other: Matrix<Double>): EjmlDoubleMatrix<DMatrixRMaj> {
-        val out = DMatrixRMaj(1, 1)
-        
-        CommonOps_DDRM.add(
-            elementAlgebra.one,
-            toEjml().origin,
-            elementAlgebra.one,
-            other.toEjml().origin, 
-            out,
-        )
-
-        return out.wrapMatrix()
-    }
-
-    override fun Point<Double>.plus(other: Point<Double>): EjmlDoubleVector<DMatrixRMaj> {
-        val out = DMatrixRMaj(1, 1)
-
-        CommonOps_DDRM.add(
-            elementAlgebra.one,
-            toEjml().origin,
-            elementAlgebra.one,
-            other.toEjml().origin,
-            out,
-        )
-
-        return out.wrapVector()
-    }
-
-    override fun Point<Double>.minus(other: Point<Double>): EjmlDoubleVector<DMatrixRMaj> {
-        val out = DMatrixRMaj(1, 1)
-
-        CommonOps_DDRM.add(
-            elementAlgebra.one,
-            toEjml().origin,
-            elementAlgebra { -one },
-            other.toEjml().origin,
-            out,
-        )
-
-        return out.wrapVector()
-    }
-
-    override fun Double.times(m: Matrix<Double>): EjmlDoubleMatrix<DMatrixRMaj> = m * this
-
-    override fun Point<Double>.times(value: Double): EjmlDoubleVector<DMatrixRMaj> {
-        val res = DMatrixRMaj(1, 1)
-        CommonOps_DDRM.scale(value, toEjml().origin, res)
-        return res.wrapVector()
-    }
-
-    override fun Double.times(v: Point<Double>): EjmlDoubleVector<DMatrixRMaj> = v * this
-
-    @UnstableKMathAPI
-    override fun <F : StructureFeature> getFeature(structure: Matrix<Double>, type: KClass<out F>): F? {
-        structure.getFeature(type)?.let { return it }
-        val origin = structure.toEjml().origin
-
-        return when (type) {
-            InverseMatrixFeature::class -> object : InverseMatrixFeature<Double> {
-                override val inverse: Matrix<Double> by lazy {
-                    val res = origin.copy()
-                    CommonOps_DDRM.invert(res)
-                    res.wrapMatrix()
-                }
-            }
-
-            DeterminantFeature::class -> object : DeterminantFeature<Double> {
-                override val determinant: Double by lazy { CommonOps_DDRM.det(origin) }
-            }
-
-            SingularValueDecompositionFeature::class -> object : SingularValueDecompositionFeature<Double> {
-                private val svd by lazy {
-                    DecompositionFactory_DDRM.svd(origin.numRows, origin.numCols, true, true, false)
-                        .apply { decompose(origin.copy()) }
-                }
-
-                override val u: Matrix<Double> by lazy { svd.getU(null, false).wrapMatrix() }
-                override val s: Matrix<Double> by lazy { svd.getW(null).wrapMatrix() }
-                override val v: Matrix<Double> by lazy { svd.getV(null, false).wrapMatrix() }
-                override val singularValues: Point<Double> by lazy { DoubleBuffer(svd.singularValues) }
-            }
-
-            QRDecompositionFeature::class -> object : QRDecompositionFeature<Double> {
-                private val qr by lazy {
-                    DecompositionFactory_DDRM.qr().apply { decompose(origin.copy()) }
-                }
-
-                override val q: Matrix<Double> by lazy {
-                    qr.getQ(null, false).wrapMatrix() + OrthogonalFeature
-                }
-
-                override val r: Matrix<Double> by lazy { qr.getR(null, false).wrapMatrix() + UFeature }
-            }
-
-            CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature<Double> {
-                override val l: Matrix<Double> by lazy {
-                    val cholesky =
-                        DecompositionFactory_DDRM.chol(structure.rowNum, true).apply { decompose(origin.copy()) }
-
-                    cholesky.getT(null).wrapMatrix() + LFeature
-                }
-            }
-
-            LupDecompositionFeature::class -> object : LupDecompositionFeature<Double> {
-                private val lup by lazy {
-                    DecompositionFactory_DDRM.lu(origin.numRows, origin.numCols).apply { decompose(origin.copy()) }
-                }
-
-                override val l: Matrix<Double> by lazy {
-                    lup.getLower(null).wrapMatrix() + LFeature
-                }
-
-                override val u: Matrix<Double> by lazy {
-                    lup.getUpper(null).wrapMatrix() + UFeature
-                }
-
-                override val p: Matrix<Double> by lazy { lup.getRowPivot(null).wrapMatrix() }
-            }
-
-            else -> null
-        }?.let(type::cast)
-    }
-
-    /**
-     * Solves for *x* in the following equation: *x = [a] <sup>-1</sup> &middot; [b]*.
-     *
-     * @param a the base matrix.
-     * @param b n by p matrix.
-     * @return the solution for *x* that is n by p.
-     */
-    public fun solve(a: Matrix<Double>, b: Matrix<Double>): EjmlDoubleMatrix<DMatrixRMaj> {
-        val res = DMatrixRMaj(1, 1)
-        CommonOps_DDRM.solve(DMatrixRMaj(a.toEjml().origin), DMatrixRMaj(b.toEjml().origin), res)
-        return res.wrapMatrix()
-    }
-
-    /**
-     * Solves for *x* in the following equation: *x = [a] <sup>-1</sup> &middot; [b]*.
-     *
-     * @param a the base matrix.
-     * @param b n by p vector.
-     * @return the solution for *x* that is n by p.
-     */
-    public fun solve(a: Matrix<Double>, b: Point<Double>): EjmlDoubleVector<DMatrixRMaj> {
-        val res = DMatrixRMaj(1, 1)
-        CommonOps_DDRM.solve(DMatrixRMaj(a.toEjml().origin), DMatrixRMaj(b.toEjml().origin), res)
-        return EjmlDoubleVector(res)
-    }
-}
-
-/**
- * [EjmlLinearSpace] implementation based on [CommonOps_FDRM], [DecompositionFactory_FDRM] operations and
- * [FMatrixRMaj] matrices.
- */
-public object EjmlLinearSpaceFDRM : EjmlLinearSpace<Float, FloatField, FMatrixRMaj>() {
-    /**
-     * The [FloatField] reference.
-     */
-    override val elementAlgebra: FloatField get() = FloatField
-
-    @Suppress("UNCHECKED_CAST")
-    override fun Matrix<Float>.toEjml(): EjmlFloatMatrix<FMatrixRMaj> = when {
-        this is EjmlFloatMatrix<*> && origin is FMatrixRMaj -> this as EjmlFloatMatrix<FMatrixRMaj>
-        else -> buildMatrix(rowNum, colNum) { i, j -> get(i, j) }
-    }
-
-    @Suppress("UNCHECKED_CAST")
-    override fun Point<Float>.toEjml(): EjmlFloatVector<FMatrixRMaj> = when {
-        this is EjmlFloatVector<*> && origin is FMatrixRMaj -> this as EjmlFloatVector<FMatrixRMaj>
-        else -> EjmlFloatVector(FMatrixRMaj(size, 1).also {
-            (0 until it.numRows).forEach { row -> it[row, 0] = get(row) }
-        })
-    }
-
-    override fun buildMatrix(
-        rows: Int,
-        columns: Int,
-        initializer: FloatField.(i: Int, j: Int) -> Float,
-    ): EjmlFloatMatrix<FMatrixRMaj> = FMatrixRMaj(rows, columns).also {
-        (0 until rows).forEach { row ->
-            (0 until columns).forEach { col -> it[row, col] = elementAlgebra.initializer(row, col) }
-        }
-    }.wrapMatrix()
-
-    override fun buildVector(
-        size: Int,
-        initializer: FloatField.(Int) -> Float,
-    ): EjmlFloatVector<FMatrixRMaj> = EjmlFloatVector(FMatrixRMaj(size, 1).also {
-        (0 until it.numRows).forEach { row -> it[row, 0] = elementAlgebra.initializer(row) }
-    })
-
-    private fun <T : FMatrix> T.wrapMatrix() = EjmlFloatMatrix(this)
-    private fun <T : FMatrix> T.wrapVector() = EjmlFloatVector(this)
-
-    override fun Matrix<Float>.unaryMinus(): Matrix<Float> = this * elementAlgebra { -one }
-
-    override fun Matrix<Float>.dot(other: Matrix<Float>): EjmlFloatMatrix<FMatrixRMaj> {
-        val out = FMatrixRMaj(1, 1)
-        CommonOps_FDRM.mult(toEjml().origin, other.toEjml().origin, out)
-        return out.wrapMatrix()
-    }
-
-    override fun Matrix<Float>.dot(vector: Point<Float>): EjmlFloatVector<FMatrixRMaj> {
-        val out = FMatrixRMaj(1, 1)
-        CommonOps_FDRM.mult(toEjml().origin, vector.toEjml().origin, out)
-        return out.wrapVector()
-    }
-
-    override operator fun Matrix<Float>.minus(other: Matrix<Float>): EjmlFloatMatrix<FMatrixRMaj> {
-        val out = FMatrixRMaj(1, 1)
-
-        CommonOps_FDRM.add(
-            elementAlgebra.one, 
-            toEjml().origin, 
-            elementAlgebra { -one },
-            other.toEjml().origin, 
-            out,
-        )
-
-        return out.wrapMatrix()
-    }
-
-    override operator fun Matrix<Float>.times(value: Float): EjmlFloatMatrix<FMatrixRMaj> {
-        val res = FMatrixRMaj(1, 1)
-        CommonOps_FDRM.scale(value, toEjml().origin, res)
-        return res.wrapMatrix()
-    }
-
-    override fun Point<Float>.unaryMinus(): EjmlFloatVector<FMatrixRMaj> {
-        val res = FMatrixRMaj(1, 1)
-        CommonOps_FDRM.changeSign(toEjml().origin, res)
-        return res.wrapVector()
-    }
-
-    override fun Matrix<Float>.plus(other: Matrix<Float>): EjmlFloatMatrix<FMatrixRMaj> {
-        val out = FMatrixRMaj(1, 1)
-        
-        CommonOps_FDRM.add(
-            elementAlgebra.one,
-            toEjml().origin,
-            elementAlgebra.one,
-            other.toEjml().origin, 
-            out,
-        )
-
-        return out.wrapMatrix()
-    }
-
-    override fun Point<Float>.plus(other: Point<Float>): EjmlFloatVector<FMatrixRMaj> {
-        val out = FMatrixRMaj(1, 1)
-
-        CommonOps_FDRM.add(
-            elementAlgebra.one,
-            toEjml().origin,
-            elementAlgebra.one,
-            other.toEjml().origin,
-            out,
-        )
-
-        return out.wrapVector()
-    }
-
-    override fun Point<Float>.minus(other: Point<Float>): EjmlFloatVector<FMatrixRMaj> {
-        val out = FMatrixRMaj(1, 1)
-
-        CommonOps_FDRM.add(
-            elementAlgebra.one,
-            toEjml().origin,
-            elementAlgebra { -one },
-            other.toEjml().origin,
-            out,
-        )
-
-        return out.wrapVector()
-    }
-
-    override fun Float.times(m: Matrix<Float>): EjmlFloatMatrix<FMatrixRMaj> = m * this
-
-    override fun Point<Float>.times(value: Float): EjmlFloatVector<FMatrixRMaj> {
-        val res = FMatrixRMaj(1, 1)
-        CommonOps_FDRM.scale(value, toEjml().origin, res)
-        return res.wrapVector()
-    }
-
-    override fun Float.times(v: Point<Float>): EjmlFloatVector<FMatrixRMaj> = v * this
-
-    @UnstableKMathAPI
-    override fun <F : StructureFeature> getFeature(structure: Matrix<Float>, type: KClass<out F>): F? {
-        structure.getFeature(type)?.let { return it }
-        val origin = structure.toEjml().origin
-
-        return when (type) {
-            InverseMatrixFeature::class -> object : InverseMatrixFeature<Float> {
-                override val inverse: Matrix<Float> by lazy {
-                    val res = origin.copy()
-                    CommonOps_FDRM.invert(res)
-                    res.wrapMatrix()
-                }
-            }
-
-            DeterminantFeature::class -> object : DeterminantFeature<Float> {
-                override val determinant: Float by lazy { CommonOps_FDRM.det(origin) }
-            }
-
-            SingularValueDecompositionFeature::class -> object : SingularValueDecompositionFeature<Float> {
-                private val svd by lazy {
-                    DecompositionFactory_FDRM.svd(origin.numRows, origin.numCols, true, true, false)
-                        .apply { decompose(origin.copy()) }
-                }
-
-                override val u: Matrix<Float> by lazy { svd.getU(null, false).wrapMatrix() }
-                override val s: Matrix<Float> by lazy { svd.getW(null).wrapMatrix() }
-                override val v: Matrix<Float> by lazy { svd.getV(null, false).wrapMatrix() }
-                override val singularValues: Point<Float> by lazy { FloatBuffer(svd.singularValues) }
-            }
-
-            QRDecompositionFeature::class -> object : QRDecompositionFeature<Float> {
-                private val qr by lazy {
-                    DecompositionFactory_FDRM.qr().apply { decompose(origin.copy()) }
-                }
-
-                override val q: Matrix<Float> by lazy {
-                    qr.getQ(null, false).wrapMatrix() + OrthogonalFeature
-                }
-
-                override val r: Matrix<Float> by lazy { qr.getR(null, false).wrapMatrix() + UFeature }
-            }
-
-            CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature<Float> {
-                override val l: Matrix<Float> by lazy {
-                    val cholesky =
-                        DecompositionFactory_FDRM.chol(structure.rowNum, true).apply { decompose(origin.copy()) }
-
-                    cholesky.getT(null).wrapMatrix() + LFeature
-                }
-            }
-
-            LupDecompositionFeature::class -> object : LupDecompositionFeature<Float> {
-                private val lup by lazy {
-                    DecompositionFactory_FDRM.lu(origin.numRows, origin.numCols).apply { decompose(origin.copy()) }
-                }
-
-                override val l: Matrix<Float> by lazy {
-                    lup.getLower(null).wrapMatrix() + LFeature
-                }
-
-                override val u: Matrix<Float> by lazy {
-                    lup.getUpper(null).wrapMatrix() + UFeature
-                }
-
-                override val p: Matrix<Float> by lazy { lup.getRowPivot(null).wrapMatrix() }
-            }
-
-            else -> null
-        }?.let(type::cast)
-    }
-
-    /**
-     * Solves for *x* in the following equation: *x = [a] <sup>-1</sup> &middot; [b]*.
-     *
-     * @param a the base matrix.
-     * @param b n by p matrix.
-     * @return the solution for *x* that is n by p.
-     */
-    public fun solve(a: Matrix<Float>, b: Matrix<Float>): EjmlFloatMatrix<FMatrixRMaj> {
-        val res = FMatrixRMaj(1, 1)
-        CommonOps_FDRM.solve(FMatrixRMaj(a.toEjml().origin), FMatrixRMaj(b.toEjml().origin), res)
-        return res.wrapMatrix()
-    }
-
-    /**
-     * Solves for *x* in the following equation: *x = [a] <sup>-1</sup> &middot; [b]*.
-     *
-     * @param a the base matrix.
-     * @param b n by p vector.
-     * @return the solution for *x* that is n by p.
-     */
-    public fun solve(a: Matrix<Float>, b: Point<Float>): EjmlFloatVector<FMatrixRMaj> {
-        val res = FMatrixRMaj(1, 1)
-        CommonOps_FDRM.solve(FMatrixRMaj(a.toEjml().origin), FMatrixRMaj(b.toEjml().origin), res)
-        return EjmlFloatVector(res)
-    }
-}
-
-/**
- * [EjmlLinearSpace] implementation based on [CommonOps_DSCC], [DecompositionFactory_DSCC] operations and
- * [DMatrixSparseCSC] matrices.
- */
-public object EjmlLinearSpaceDSCC : EjmlLinearSpace<Double, DoubleField, DMatrixSparseCSC>() {
-    /**
-     * The [DoubleField] reference.
-     */
-    override val elementAlgebra: DoubleField get() = DoubleField
-
-    @Suppress("UNCHECKED_CAST")
-    override fun Matrix<Double>.toEjml(): EjmlDoubleMatrix<DMatrixSparseCSC> = when {
-        this is EjmlDoubleMatrix<*> && origin is DMatrixSparseCSC -> this as EjmlDoubleMatrix<DMatrixSparseCSC>
-        else -> buildMatrix(rowNum, colNum) { i, j -> get(i, j) }
-    }
-
-    @Suppress("UNCHECKED_CAST")
-    override fun Point<Double>.toEjml(): EjmlDoubleVector<DMatrixSparseCSC> = when {
-        this is EjmlDoubleVector<*> && origin is DMatrixSparseCSC -> this as EjmlDoubleVector<DMatrixSparseCSC>
-        else -> EjmlDoubleVector(DMatrixSparseCSC(size, 1).also {
-            (0 until it.numRows).forEach { row -> it[row, 0] = get(row) }
-        })
-    }
-
-    override fun buildMatrix(
-        rows: Int,
-        columns: Int,
-        initializer: DoubleField.(i: Int, j: Int) -> Double,
-    ): EjmlDoubleMatrix<DMatrixSparseCSC> = DMatrixSparseCSC(rows, columns).also {
-        (0 until rows).forEach { row ->
-            (0 until columns).forEach { col -> it[row, col] = elementAlgebra.initializer(row, col) }
-        }
-    }.wrapMatrix()
-
-    override fun buildVector(
-        size: Int,
-        initializer: DoubleField.(Int) -> Double,
-    ): EjmlDoubleVector<DMatrixSparseCSC> = EjmlDoubleVector(DMatrixSparseCSC(size, 1).also {
-        (0 until it.numRows).forEach { row -> it[row, 0] = elementAlgebra.initializer(row) }
-    })
-
-    private fun <T : DMatrix> T.wrapMatrix() = EjmlDoubleMatrix(this)
-    private fun <T : DMatrix> T.wrapVector() = EjmlDoubleVector(this)
-
-    override fun Matrix<Double>.unaryMinus(): Matrix<Double> = this * elementAlgebra { -one }
-
-    override fun Matrix<Double>.dot(other: Matrix<Double>): EjmlDoubleMatrix<DMatrixSparseCSC> {
-        val out = DMatrixSparseCSC(1, 1)
-        CommonOps_DSCC.mult(toEjml().origin, other.toEjml().origin, out)
-        return out.wrapMatrix()
-    }
-
-    override fun Matrix<Double>.dot(vector: Point<Double>): EjmlDoubleVector<DMatrixSparseCSC> {
-        val out = DMatrixSparseCSC(1, 1)
-        CommonOps_DSCC.mult(toEjml().origin, vector.toEjml().origin, out)
-        return out.wrapVector()
-    }
-
-    override operator fun Matrix<Double>.minus(other: Matrix<Double>): EjmlDoubleMatrix<DMatrixSparseCSC> {
-        val out = DMatrixSparseCSC(1, 1)
-
-        CommonOps_DSCC.add(
-            elementAlgebra.one, 
-            toEjml().origin, 
-            elementAlgebra { -one },
-            other.toEjml().origin, 
-            out,
-            null, 
-            null,
-        )
-
-        return out.wrapMatrix()
-    }
-
-    override operator fun Matrix<Double>.times(value: Double): EjmlDoubleMatrix<DMatrixSparseCSC> {
-        val res = DMatrixSparseCSC(1, 1)
-        CommonOps_DSCC.scale(value, toEjml().origin, res)
-        return res.wrapMatrix()
-    }
-
-    override fun Point<Double>.unaryMinus(): EjmlDoubleVector<DMatrixSparseCSC> {
-        val res = DMatrixSparseCSC(1, 1)
-        CommonOps_DSCC.changeSign(toEjml().origin, res)
-        return res.wrapVector()
-    }
-
-    override fun Matrix<Double>.plus(other: Matrix<Double>): EjmlDoubleMatrix<DMatrixSparseCSC> {
-        val out = DMatrixSparseCSC(1, 1)
-        
-        CommonOps_DSCC.add(
-            elementAlgebra.one,
-            toEjml().origin,
-            elementAlgebra.one,
-            other.toEjml().origin, 
-            out,
-            null, 
-            null,
-        )
-
-        return out.wrapMatrix()
-    }
-
-    override fun Point<Double>.plus(other: Point<Double>): EjmlDoubleVector<DMatrixSparseCSC> {
-        val out = DMatrixSparseCSC(1, 1)
-
-        CommonOps_DSCC.add(
-            elementAlgebra.one,
-            toEjml().origin,
-            elementAlgebra.one,
-            other.toEjml().origin,
-            out,
-            null, 
-            null,
-        )
-
-        return out.wrapVector()
-    }
-
-    override fun Point<Double>.minus(other: Point<Double>): EjmlDoubleVector<DMatrixSparseCSC> {
-        val out = DMatrixSparseCSC(1, 1)
-
-        CommonOps_DSCC.add(
-            elementAlgebra.one,
-            toEjml().origin,
-            elementAlgebra { -one },
-            other.toEjml().origin,
-            out,
-            null, 
-            null,
-        )
-
-        return out.wrapVector()
-    }
-
-    override fun Double.times(m: Matrix<Double>): EjmlDoubleMatrix<DMatrixSparseCSC> = m * this
-
-    override fun Point<Double>.times(value: Double): EjmlDoubleVector<DMatrixSparseCSC> {
-        val res = DMatrixSparseCSC(1, 1)
-        CommonOps_DSCC.scale(value, toEjml().origin, res)
-        return res.wrapVector()
-    }
-
-    override fun Double.times(v: Point<Double>): EjmlDoubleVector<DMatrixSparseCSC> = v * this
-
-    @UnstableKMathAPI
-    override fun <F : StructureFeature> getFeature(structure: Matrix<Double>, type: KClass<out F>): F? {
-        structure.getFeature(type)?.let { return it }
-        val origin = structure.toEjml().origin
-
-        return when (type) {
-            QRDecompositionFeature::class -> object : QRDecompositionFeature<Double> {
-                private val qr by lazy {
-                    DecompositionFactory_DSCC.qr(FillReducing.NONE).apply { decompose(origin.copy()) }
-                }
-
-                override val q: Matrix<Double> by lazy {
-                    qr.getQ(null, false).wrapMatrix() + OrthogonalFeature
-                }
-
-                override val r: Matrix<Double> by lazy { qr.getR(null, false).wrapMatrix() + UFeature }
-            }
-
-            CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature<Double> {
-                override val l: Matrix<Double> by lazy {
-                    val cholesky =
-                        DecompositionFactory_DSCC.cholesky().apply { decompose(origin.copy()) }
-
-                    (cholesky.getT(null) as DMatrix).wrapMatrix() + LFeature
-                }
-            }
-
-            LUDecompositionFeature::class, DeterminantFeature::class, InverseMatrixFeature::class -> object :
-                LUDecompositionFeature<Double>, DeterminantFeature<Double>, InverseMatrixFeature<Double> {
-                private val lu by lazy {
-                    DecompositionFactory_DSCC.lu(FillReducing.NONE).apply { decompose(origin.copy()) }
-                }
-
-                override val l: Matrix<Double> by lazy {
-                    lu.getLower(null).wrapMatrix() + LFeature
-                }
-
-                override val u: Matrix<Double> by lazy {
-                    lu.getUpper(null).wrapMatrix() + UFeature
-                }
-
-                override val inverse: Matrix<Double> by lazy {
-                    var a = origin
-                    val inverse = DMatrixRMaj(1, 1)
-                    val solver = LinearSolverFactory_DSCC.lu(FillReducing.NONE)
-                    if (solver.modifiesA()) a = a.copy()
-                    val i = CommonOps_DDRM.identity(a.numRows)
-                    solver.solve(i, inverse)
-                    inverse.wrapMatrix()
-                }
-
-                override val determinant: Double by lazy { elementAlgebra.number(lu.computeDeterminant().real) }
-            }
-
-            else -> null
-        }?.let(type::cast)
-    }
-
-    /**
-     * Solves for *x* in the following equation: *x = [a] <sup>-1</sup> &middot; [b]*.
-     *
-     * @param a the base matrix.
-     * @param b n by p matrix.
-     * @return the solution for *x* that is n by p.
-     */
-    public fun solve(a: Matrix<Double>, b: Matrix<Double>): EjmlDoubleMatrix<DMatrixSparseCSC> {
-        val res = DMatrixSparseCSC(1, 1)
-        CommonOps_DSCC.solve(DMatrixSparseCSC(a.toEjml().origin), DMatrixSparseCSC(b.toEjml().origin), res)
-        return res.wrapMatrix()
-    }
-
-    /**
-     * Solves for *x* in the following equation: *x = [a] <sup>-1</sup> &middot; [b]*.
-     *
-     * @param a the base matrix.
-     * @param b n by p vector.
-     * @return the solution for *x* that is n by p.
-     */
-    public fun solve(a: Matrix<Double>, b: Point<Double>): EjmlDoubleVector<DMatrixSparseCSC> {
-        val res = DMatrixSparseCSC(1, 1)
-        CommonOps_DSCC.solve(DMatrixSparseCSC(a.toEjml().origin), DMatrixSparseCSC(b.toEjml().origin), res)
-        return EjmlDoubleVector(res)
-    }
-}
-
-/**
- * [EjmlLinearSpace] implementation based on [CommonOps_FSCC], [DecompositionFactory_FSCC] operations and
- * [FMatrixSparseCSC] matrices.
- */
-public object EjmlLinearSpaceFSCC : EjmlLinearSpace<Float, FloatField, FMatrixSparseCSC>() {
-    /**
-     * The [FloatField] reference.
-     */
-    override val elementAlgebra: FloatField get() = FloatField
-
-    @Suppress("UNCHECKED_CAST")
-    override fun Matrix<Float>.toEjml(): EjmlFloatMatrix<FMatrixSparseCSC> = when {
-        this is EjmlFloatMatrix<*> && origin is FMatrixSparseCSC -> this as EjmlFloatMatrix<FMatrixSparseCSC>
-        else -> buildMatrix(rowNum, colNum) { i, j -> get(i, j) }
-    }
-
-    @Suppress("UNCHECKED_CAST")
-    override fun Point<Float>.toEjml(): EjmlFloatVector<FMatrixSparseCSC> = when {
-        this is EjmlFloatVector<*> && origin is FMatrixSparseCSC -> this as EjmlFloatVector<FMatrixSparseCSC>
-        else -> EjmlFloatVector(FMatrixSparseCSC(size, 1).also {
-            (0 until it.numRows).forEach { row -> it[row, 0] = get(row) }
-        })
-    }
-
-    override fun buildMatrix(
-        rows: Int,
-        columns: Int,
-        initializer: FloatField.(i: Int, j: Int) -> Float,
-    ): EjmlFloatMatrix<FMatrixSparseCSC> = FMatrixSparseCSC(rows, columns).also {
-        (0 until rows).forEach { row ->
-            (0 until columns).forEach { col -> it[row, col] = elementAlgebra.initializer(row, col) }
-        }
-    }.wrapMatrix()
-
-    override fun buildVector(
-        size: Int,
-        initializer: FloatField.(Int) -> Float,
-    ): EjmlFloatVector<FMatrixSparseCSC> = EjmlFloatVector(FMatrixSparseCSC(size, 1).also {
-        (0 until it.numRows).forEach { row -> it[row, 0] = elementAlgebra.initializer(row) }
-    })
-
-    private fun <T : FMatrix> T.wrapMatrix() = EjmlFloatMatrix(this)
-    private fun <T : FMatrix> T.wrapVector() = EjmlFloatVector(this)
-
-    override fun Matrix<Float>.unaryMinus(): Matrix<Float> = this * elementAlgebra { -one }
-
-    override fun Matrix<Float>.dot(other: Matrix<Float>): EjmlFloatMatrix<FMatrixSparseCSC> {
-        val out = FMatrixSparseCSC(1, 1)
-        CommonOps_FSCC.mult(toEjml().origin, other.toEjml().origin, out)
-        return out.wrapMatrix()
-    }
-
-    override fun Matrix<Float>.dot(vector: Point<Float>): EjmlFloatVector<FMatrixSparseCSC> {
-        val out = FMatrixSparseCSC(1, 1)
-        CommonOps_FSCC.mult(toEjml().origin, vector.toEjml().origin, out)
-        return out.wrapVector()
-    }
-
-    override operator fun Matrix<Float>.minus(other: Matrix<Float>): EjmlFloatMatrix<FMatrixSparseCSC> {
-        val out = FMatrixSparseCSC(1, 1)
-
-        CommonOps_FSCC.add(
-            elementAlgebra.one, 
-            toEjml().origin, 
-            elementAlgebra { -one },
-            other.toEjml().origin, 
-            out,
-            null, 
-            null,
-        )
-
-        return out.wrapMatrix()
-    }
-
-    override operator fun Matrix<Float>.times(value: Float): EjmlFloatMatrix<FMatrixSparseCSC> {
-        val res = FMatrixSparseCSC(1, 1)
-        CommonOps_FSCC.scale(value, toEjml().origin, res)
-        return res.wrapMatrix()
-    }
-
-    override fun Point<Float>.unaryMinus(): EjmlFloatVector<FMatrixSparseCSC> {
-        val res = FMatrixSparseCSC(1, 1)
-        CommonOps_FSCC.changeSign(toEjml().origin, res)
-        return res.wrapVector()
-    }
-
-    override fun Matrix<Float>.plus(other: Matrix<Float>): EjmlFloatMatrix<FMatrixSparseCSC> {
-        val out = FMatrixSparseCSC(1, 1)
-        
-        CommonOps_FSCC.add(
-            elementAlgebra.one,
-            toEjml().origin,
-            elementAlgebra.one,
-            other.toEjml().origin, 
-            out,
-            null, 
-            null,
-        )
-
-        return out.wrapMatrix()
-    }
-
-    override fun Point<Float>.plus(other: Point<Float>): EjmlFloatVector<FMatrixSparseCSC> {
-        val out = FMatrixSparseCSC(1, 1)
-
-        CommonOps_FSCC.add(
-            elementAlgebra.one,
-            toEjml().origin,
-            elementAlgebra.one,
-            other.toEjml().origin,
-            out,
-            null, 
-            null,
-        )
-
-        return out.wrapVector()
-    }
-
-    override fun Point<Float>.minus(other: Point<Float>): EjmlFloatVector<FMatrixSparseCSC> {
-        val out = FMatrixSparseCSC(1, 1)
-
-        CommonOps_FSCC.add(
-            elementAlgebra.one,
-            toEjml().origin,
-            elementAlgebra { -one },
-            other.toEjml().origin,
-            out,
-            null, 
-            null,
-        )
-
-        return out.wrapVector()
-    }
-
-    override fun Float.times(m: Matrix<Float>): EjmlFloatMatrix<FMatrixSparseCSC> = m * this
-
-    override fun Point<Float>.times(value: Float): EjmlFloatVector<FMatrixSparseCSC> {
-        val res = FMatrixSparseCSC(1, 1)
-        CommonOps_FSCC.scale(value, toEjml().origin, res)
-        return res.wrapVector()
-    }
-
-    override fun Float.times(v: Point<Float>): EjmlFloatVector<FMatrixSparseCSC> = v * this
-
-    @UnstableKMathAPI
-    override fun <F : StructureFeature> getFeature(structure: Matrix<Float>, type: KClass<out F>): F? {
-        structure.getFeature(type)?.let { return it }
-        val origin = structure.toEjml().origin
-
-        return when (type) {
-            QRDecompositionFeature::class -> object : QRDecompositionFeature<Float> {
-                private val qr by lazy {
-                    DecompositionFactory_FSCC.qr(FillReducing.NONE).apply { decompose(origin.copy()) }
-                }
-
-                override val q: Matrix<Float> by lazy {
-                    qr.getQ(null, false).wrapMatrix() + OrthogonalFeature
-                }
-
-                override val r: Matrix<Float> by lazy { qr.getR(null, false).wrapMatrix() + UFeature }
-            }
-
-            CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature<Float> {
-                override val l: Matrix<Float> by lazy {
-                    val cholesky =
-                        DecompositionFactory_FSCC.cholesky().apply { decompose(origin.copy()) }
-
-                    (cholesky.getT(null) as FMatrix).wrapMatrix() + LFeature
-                }
-            }
-
-            LUDecompositionFeature::class, DeterminantFeature::class, InverseMatrixFeature::class -> object :
-                LUDecompositionFeature<Float>, DeterminantFeature<Float>, InverseMatrixFeature<Float> {
-                private val lu by lazy {
-                    DecompositionFactory_FSCC.lu(FillReducing.NONE).apply { decompose(origin.copy()) }
-                }
-
-                override val l: Matrix<Float> by lazy {
-                    lu.getLower(null).wrapMatrix() + LFeature
-                }
-
-                override val u: Matrix<Float> by lazy {
-                    lu.getUpper(null).wrapMatrix() + UFeature
-                }
-
-                override val inverse: Matrix<Float> by lazy {
-                    var a = origin
-                    val inverse = FMatrixRMaj(1, 1)
-                    val solver = LinearSolverFactory_FSCC.lu(FillReducing.NONE)
-                    if (solver.modifiesA()) a = a.copy()
-                    val i = CommonOps_FDRM.identity(a.numRows)
-                    solver.solve(i, inverse)
-                    inverse.wrapMatrix()
-                }
-
-                override val determinant: Float by lazy { elementAlgebra.number(lu.computeDeterminant().real) }
-            }
-
-            else -> null
-        }?.let(type::cast)
-    }
-
-    /**
-     * Solves for *x* in the following equation: *x = [a] <sup>-1</sup> &middot; [b]*.
-     *
-     * @param a the base matrix.
-     * @param b n by p matrix.
-     * @return the solution for *x* that is n by p.
-     */
-    public fun solve(a: Matrix<Float>, b: Matrix<Float>): EjmlFloatMatrix<FMatrixSparseCSC> {
-        val res = FMatrixSparseCSC(1, 1)
-        CommonOps_FSCC.solve(FMatrixSparseCSC(a.toEjml().origin), FMatrixSparseCSC(b.toEjml().origin), res)
-        return res.wrapMatrix()
-    }
-
-    /**
-     * Solves for *x* in the following equation: *x = [a] <sup>-1</sup> &middot; [b]*.
-     *
-     * @param a the base matrix.
-     * @param b n by p vector.
-     * @return the solution for *x* that is n by p.
-     */
-    public fun solve(a: Matrix<Float>, b: Point<Float>): EjmlFloatVector<FMatrixSparseCSC> {
-        val res = FMatrixSparseCSC(1, 1)
-        CommonOps_FSCC.solve(FMatrixSparseCSC(a.toEjml().origin), FMatrixSparseCSC(b.toEjml().origin), res)
-        return EjmlFloatVector(res)
-    }
-}
-

kmath-ejml/src/test/kotlin/space/kscience/kmath/ejml/EjmlMatrixTest.kt

@@ -11,7 +11,7 @@ import org.ejml.dense.row.RandomMatrices_DDRM
 import org.ejml.dense.row.factory.DecompositionFactory_DDRM
 import space.kscience.kmath.linear.DeterminantFeature
 import space.kscience.kmath.linear.LupDecompositionFeature
-import space.kscience.kmath.linear.getFeature
+import space.kscience.kmath.linear.computeFeature
 import space.kscience.kmath.misc.PerformancePitfall
 import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.nd.StructureND
@@ -59,9 +59,9 @@ internal class EjmlMatrixTest {
     fun features() {
         val m = randomMatrix
         val w = EjmlDoubleMatrix(m)
-        val det: DeterminantFeature<Double> = EjmlLinearSpaceDDRM.getFeature(w) ?: fail()
+        val det: DeterminantFeature<Double> = EjmlLinearSpaceDDRM.computeFeature(w) ?: fail()
         assertEquals(CommonOps_DDRM.det(m), det.determinant)
-        val lup: LupDecompositionFeature<Double> = EjmlLinearSpaceDDRM.getFeature(w) ?: fail()
+        val lup: LupDecompositionFeature<Double> = EjmlLinearSpaceDDRM.computeFeature(w) ?: fail()
 
         val ludecompositionF64 = DecompositionFactory_DDRM.lu(m.numRows, m.numCols)
             .also { it.decompose(m.copy()) }

kmath-for-real/src/commonMain/kotlin/space/kscience/kmath/real/RealMatrix.kt

@@ -32,18 +32,18 @@ import kotlin.math.pow
 public typealias RealMatrix = Matrix<Double>
 
 public fun realMatrix(rowNum: Int, colNum: Int, initializer: DoubleField.(i: Int, j: Int) -> Double): RealMatrix =
-    LinearSpace.real.buildMatrix(rowNum, colNum, initializer)
+    LinearSpace.double.buildMatrix(rowNum, colNum, initializer)
 
 @OptIn(UnstableKMathAPI::class)
 public fun realMatrix(rowNum: Int, colNum: Int): MatrixBuilder<Double, DoubleField> =
-    LinearSpace.real.matrix(rowNum, colNum)
+    LinearSpace.double.matrix(rowNum, colNum)
 
 public fun Array<DoubleArray>.toMatrix(): RealMatrix {
-    return LinearSpace.real.buildMatrix(size, this[0].size) { row, col -> this@toMatrix[row][col] }
+    return LinearSpace.double.buildMatrix(size, this[0].size) { row, col -> this@toMatrix[row][col] }
 }
 
 public fun Sequence<DoubleArray>.toMatrix(): RealMatrix = toList().let {
-    LinearSpace.real.buildMatrix(it.size, it[0].size) { row, col -> it[row][col] }
+    LinearSpace.double.buildMatrix(it.size, it[0].size) { row, col -> it[row][col] }
 }
 
 public fun RealMatrix.repeatStackVertical(n: Int): RealMatrix =
@@ -56,37 +56,37 @@ public fun RealMatrix.repeatStackVertical(n: Int): RealMatrix =
  */
 
 public operator fun RealMatrix.times(double: Double): RealMatrix =
-    LinearSpace.real.buildMatrix(rowNum, colNum) { row, col ->
+    LinearSpace.double.buildMatrix(rowNum, colNum) { row, col ->
         get(row, col) * double
     }
 
 public operator fun RealMatrix.plus(double: Double): RealMatrix =
-    LinearSpace.real.buildMatrix(rowNum, colNum) { row, col ->
+    LinearSpace.double.buildMatrix(rowNum, colNum) { row, col ->
         get(row, col) + double
     }
 
 public operator fun RealMatrix.minus(double: Double): RealMatrix =
-    LinearSpace.real.buildMatrix(rowNum, colNum) { row, col ->
+    LinearSpace.double.buildMatrix(rowNum, colNum) { row, col ->
         get(row, col) - double
     }
 
 public operator fun RealMatrix.div(double: Double): RealMatrix =
-    LinearSpace.real.buildMatrix(rowNum, colNum) { row, col ->
+    LinearSpace.double.buildMatrix(rowNum, colNum) { row, col ->
         get(row, col) / double
     }
 
 public operator fun Double.times(matrix: RealMatrix): RealMatrix =
-    LinearSpace.real.buildMatrix(matrix.rowNum, matrix.colNum) { row, col ->
+    LinearSpace.double.buildMatrix(matrix.rowNum, matrix.colNum) { row, col ->
         this@times * matrix[row, col]
     }
 
 public operator fun Double.plus(matrix: RealMatrix): RealMatrix =
-    LinearSpace.real.buildMatrix(matrix.rowNum, matrix.colNum) { row, col ->
+    LinearSpace.double.buildMatrix(matrix.rowNum, matrix.colNum) { row, col ->
         this@plus + matrix[row, col]
     }
 
 public operator fun Double.minus(matrix: RealMatrix): RealMatrix =
-    LinearSpace.real.buildMatrix(matrix.rowNum, matrix.colNum) { row, col ->
+    LinearSpace.double.buildMatrix(matrix.rowNum, matrix.colNum) { row, col ->
         this@minus - matrix[row, col]
     }
 
@@ -101,20 +101,20 @@ public operator fun Double.minus(matrix: RealMatrix): RealMatrix =
 
 @UnstableKMathAPI
 public operator fun RealMatrix.times(other: RealMatrix): RealMatrix =
-    LinearSpace.real.buildMatrix(rowNum, colNum) { row, col -> this@times[row, col] * other[row, col] }
+    LinearSpace.double.buildMatrix(rowNum, colNum) { row, col -> this@times[row, col] * other[row, col] }
 
 public operator fun RealMatrix.plus(other: RealMatrix): RealMatrix =
-    LinearSpace.real.run { this@plus + other }
+    LinearSpace.double.run { this@plus + other }
 
 public operator fun RealMatrix.minus(other: RealMatrix): RealMatrix =
-    LinearSpace.real.buildMatrix(rowNum, colNum) { row, col -> this@minus[row, col] - other[row, col] }
+    LinearSpace.double.buildMatrix(rowNum, colNum) { row, col -> this@minus[row, col] - other[row, col] }
 
 /*
  *  Operations on columns
  */
 
 public inline fun RealMatrix.appendColumn(crossinline mapper: (Buffer<Double>) -> Double): RealMatrix =
-    LinearSpace.real.buildMatrix(rowNum, colNum + 1) { row, col ->
+    LinearSpace.double.buildMatrix(rowNum, colNum + 1) { row, col ->
         if (col < colNum)
             get(row, col)
         else
@@ -122,7 +122,7 @@ public inline fun RealMatrix.appendColumn(crossinline mapper: (Buffer<Double>) -
     }
 
 public fun RealMatrix.extractColumns(columnRange: IntRange): RealMatrix =
-    LinearSpace.real.buildMatrix(rowNum, columnRange.count()) { row, col ->
+    LinearSpace.double.buildMatrix(rowNum, columnRange.count()) { row, col ->
         this@extractColumns[row, columnRange.first + col]
     }
 
@@ -155,14 +155,14 @@ public fun RealMatrix.max(): Double? = elements().map { (_, value) -> value }.ma
 public fun RealMatrix.average(): Double = elements().map { (_, value) -> value }.average()
 
 public inline fun RealMatrix.map(crossinline transform: (Double) -> Double): RealMatrix =
-    LinearSpace.real.buildMatrix(rowNum, colNum) { i, j ->
+    LinearSpace.double.buildMatrix(rowNum, colNum) { i, j ->
         transform(get(i, j))
     }
 
 /**
  * Inverse a square real matrix using LUP decomposition
  */
-public fun RealMatrix.inverseWithLup(): RealMatrix = LinearSpace.real.inverseWithLup(this)
+public fun RealMatrix.inverseWithLup(): RealMatrix = LinearSpace.double.lupSolver().inverse(this)
 
 //extended operations
 

kmath-for-real/src/commonMain/kotlin/space/kscience/kmath/real/RealVector.kt

@@ -8,11 +8,8 @@ package space.kscience.kmath.real
 import space.kscience.kmath.linear.Point
 import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.operations.Norm
-import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.structures.*
 import space.kscience.kmath.structures.MutableBuffer.Companion.double
-import space.kscience.kmath.structures.asBuffer
-import space.kscience.kmath.structures.fold
-import space.kscience.kmath.structures.indices
 import kotlin.math.pow
 import kotlin.math.sqrt
 
@@ -105,8 +102,4 @@ public fun DoubleVector.sum(): Double {
     return res
 }
 
-public object VectorL2Norm : Norm<DoubleVector, Double> {
+public val DoubleVector.norm: Double get() = DoubleL2Norm.norm(this)
-    override fun norm(arg: DoubleVector): Double = sqrt(arg.fold(0.0) { acc: Double, d: Double -> acc + d.pow(2) })
-}
-
-public val DoubleVector.norm: Double get() = VectorL2Norm.norm(this)

kmath-for-real/src/commonMain/kotlin/space/kscience/kmath/real/dot.kt

@@ -12,6 +12,6 @@ import space.kscience.kmath.linear.Matrix
 /**
  * Optimized dot product for real matrices
  */
-public infix fun Matrix<Double>.dot(other: Matrix<Double>): Matrix<Double> = LinearSpace.real.run {
+public infix fun Matrix<Double>.dot(other: Matrix<Double>): Matrix<Double> = LinearSpace.double.run {
     this@dot dot other
 }

kmath-for-real/src/commonTest/kotlin/space/kscience/kmath/real/DoubleMatrixTest.kt

@@ -65,7 +65,7 @@ internal class DoubleMatrixTest {
             4.0, 6.0, 2.0
         )
         val matrix2 = (matrix1 * 2.5 + 1.0 - 2.0) / 2.0
-        val expectedResult = LinearSpace.real.matrix(2, 3)(
+        val expectedResult = LinearSpace.double.matrix(2, 3)(
             0.75, -0.5, 3.25,
             4.5, 7.0, 2.0
         )
@@ -160,7 +160,7 @@ internal class DoubleMatrixTest {
 
     @Test
     fun testAllElementOperations() {
-        val matrix1 = LinearSpace.real.matrix(2, 4)(
+        val matrix1 = LinearSpace.double.matrix(2, 4)(
             -1.0, 0.0, 3.0, 15.0,
             4.0, -6.0, 7.0, -11.0
         )

kmath-for-real/src/commonTest/kotlin/space/kscience/kmath/real/DoubleVectorTest.kt

@@ -35,7 +35,7 @@ internal class DoubleVectorTest {
         val vector2 = DoubleBuffer(5) { 5 - it.toDouble() }
         val matrix1 = vector1.asMatrix()
         val matrix2 = vector2.asMatrix().transpose()
-        val product = LinearSpace.real.run { matrix1 dot matrix2 }
+        val product = LinearSpace.double.run { matrix1 dot matrix2 }
         assertEquals(5.0, product[1, 0])
         assertEquals(6.0, product[2, 2])
     }

kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/GaussIntegrator.kt

@@ -53,7 +53,7 @@ public class GaussIntegrator<T : Any>(
         }
     }
 
-    override fun integrate(integrand: UnivariateIntegrand<T>): UnivariateIntegrand<T> = with(algebra) {
+    override fun process(integrand: UnivariateIntegrand<T>): UnivariateIntegrand<T> = with(algebra) {
         val f = integrand.function
         val (points, weights) = buildRule(integrand)
         var res = zero
@@ -73,7 +73,8 @@ public class GaussIntegrator<T : Any>(
 }
 
 /**
- * Create a Gauss-Legendre integrator for this field.
+ * Create a Gauss integrator for this field. By default, uses Legendre rule to compute points and weights.
+ * Custom rules could be provided by [GaussIntegratorRuleFactory] feature.
  * @see [GaussIntegrator]
  */
 public val <T : Any> Field<T>.gaussIntegrator: GaussIntegrator<T> get() = GaussIntegrator(this)
@@ -97,7 +98,7 @@ public fun <T : Any> GaussIntegrator<T>.integrate(
     val ranges = UnivariateIntegrandRanges(
         (0 until intervals).map { i -> (range.start + rangeSize * i)..(range.start + rangeSize * (i + 1)) to order }
     )
-    return integrate(
+    return process(
         UnivariateIntegrand(
             function,
             IntegrationRange(range),

kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/Integrand.kt

@@ -5,15 +5,18 @@
 
 package space.kscience.kmath.integration
 
+import space.kscience.kmath.misc.Feature
+import space.kscience.kmath.misc.FeatureSet
+import space.kscience.kmath.misc.Featured
 import kotlin.reflect.KClass
 
-public interface IntegrandFeature {
+public interface IntegrandFeature : Feature<IntegrandFeature> {
     override fun toString(): String
 }
 
-public interface Integrand {
+public interface Integrand : Featured<IntegrandFeature> {
-    public val features: Set<IntegrandFeature>
+    public val features: FeatureSet<IntegrandFeature>
-    public fun <T : IntegrandFeature> getFeature(type: KClass<T>): T?
+    override fun <T : IntegrandFeature> getFeature(type: KClass<out T>): T? = features.getFeature(type)
 }
 
 public inline fun <reified T : IntegrandFeature> Integrand.getFeature(): T? = getFeature(T::class)

kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/Integrator.kt

@@ -12,5 +12,5 @@ public interface Integrator<I : Integrand> {
     /**
      * Runs one integration pass and return a new [Integrand] with a new set of features.
      */
-    public fun integrate(integrand: I): I
+    public fun process(integrand: I): I
 }

kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/MultivariateIntegrand.kt

@@ -6,29 +6,21 @@
 package space.kscience.kmath.integration
 
 import space.kscience.kmath.linear.Point
-import kotlin.reflect.KClass
+import space.kscience.kmath.misc.FeatureSet
 
 public class MultivariateIntegrand<T : Any> internal constructor(
-    private val featureMap: Map<KClass<*>, IntegrandFeature>,
+    override val features: FeatureSet<IntegrandFeature>,
     public val function: (Point<T>) -> T,
 ) : Integrand {
 
-    override val features: Set<IntegrandFeature> get() = featureMap.values.toSet()
-
-    @Suppress("UNCHECKED_CAST")
-    override fun <T : IntegrandFeature> getFeature(type: KClass<T>): T? = featureMap[type] as? T
-
-    public operator fun <F : IntegrandFeature> plus(pair: Pair<KClass<out F>, F>): MultivariateIntegrand<T> =
-        MultivariateIntegrand(featureMap + pair, function)
-
     public operator fun <F : IntegrandFeature> plus(feature: F): MultivariateIntegrand<T> =
-        plus(feature::class to feature)
+        MultivariateIntegrand(features.with(feature), function)
 }
 
 @Suppress("FunctionName")
 public fun <T : Any> MultivariateIntegrand(
     vararg features: IntegrandFeature,
     function: (Point<T>) -> T,
-): MultivariateIntegrand<T> = MultivariateIntegrand(features.associateBy { it::class }, function)
+): MultivariateIntegrand<T> = MultivariateIntegrand(FeatureSet.of(*features), function)
 
 public val <T : Any> MultivariateIntegrand<T>.value: T? get() = getFeature<IntegrandValue<T>>()?.value

kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/SimpsonIntegrator.kt

@@ -44,7 +44,7 @@ public class SimpsonIntegrator<T : Any>(
         return res
     }
 
-    override fun integrate(integrand: UnivariateIntegrand<T>): UnivariateIntegrand<T> {
+    override fun process(integrand: UnivariateIntegrand<T>): UnivariateIntegrand<T> {
         val ranges = integrand.getFeature<UnivariateIntegrandRanges>()
         return if (ranges != null) {
             val res = algebra.sum(ranges.ranges.map { integrateRange(integrand, it.first, it.second) })
@@ -90,7 +90,7 @@ public object DoubleSimpsonIntegrator : UnivariateIntegrator<Double> {
         return res
     }
 
-    override fun integrate(integrand: UnivariateIntegrand<Double>): UnivariateIntegrand<Double> {
+    override fun process(integrand: UnivariateIntegrand<Double>): UnivariateIntegrand<Double> {
         val ranges = integrand.getFeature<UnivariateIntegrandRanges>()
         return if (ranges != null) {
             val res = ranges.ranges.sumOf { integrateRange(integrand, it.first, it.second) }

kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/SplineIntegrator.kt

@@ -54,7 +54,7 @@ public class SplineIntegrator<T : Comparable<T>>(
     public val algebra: Field<T>,
     public val bufferFactory: MutableBufferFactory<T>,
 ) : UnivariateIntegrator<T> {
-    override fun integrate(integrand: UnivariateIntegrand<T>): UnivariateIntegrand<T> = algebra {
+    override fun process(integrand: UnivariateIntegrand<T>): UnivariateIntegrand<T> = algebra {
         val range = integrand.getFeature<IntegrationRange>()?.range ?: 0.0..1.0
 
         val interpolator: PolynomialInterpolator<T> = SplineInterpolator(algebra, bufferFactory)
@@ -83,7 +83,7 @@ public class SplineIntegrator<T : Comparable<T>>(
  */
 @UnstableKMathAPI
 public object DoubleSplineIntegrator : UnivariateIntegrator<Double> {
-    override fun integrate(integrand: UnivariateIntegrand<Double>): UnivariateIntegrand<Double> {
+    override fun process(integrand: UnivariateIntegrand<Double>): UnivariateIntegrand<Double> {
         val range = integrand.getFeature<IntegrationRange>()?.range ?: 0.0..1.0
         val interpolator: PolynomialInterpolator<Double> = SplineInterpolator(DoubleField, ::DoubleBuffer)
 

kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/UnivariateIntegrand.kt

@@ -5,33 +5,24 @@
 
 package space.kscience.kmath.integration
 
+import space.kscience.kmath.misc.FeatureSet
 import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.structures.Buffer
 import space.kscience.kmath.structures.DoubleBuffer
-import kotlin.reflect.KClass
 
 public class UnivariateIntegrand<T> internal constructor(
-    private val featureMap: Map<KClass<*>, IntegrandFeature>,
+    override val features: FeatureSet<IntegrandFeature>,
     public val function: (Double) -> T,
 ) : Integrand {
-
-    override val features: Set<IntegrandFeature> get() = featureMap.values.toSet()
-
-    @Suppress("UNCHECKED_CAST")
-    override fun <T : IntegrandFeature> getFeature(type: KClass<T>): T? = featureMap[type] as? T
-
-    public operator fun <F : IntegrandFeature> plus(pair: Pair<KClass<out F>, F>): UnivariateIntegrand<T> =
-        UnivariateIntegrand(featureMap + pair, function)
-
     public operator fun <F : IntegrandFeature> plus(feature: F): UnivariateIntegrand<T> =
-        plus(feature::class to feature)
+        UnivariateIntegrand(features.with(feature), function)
 }
 
 @Suppress("FunctionName")
 public fun <T : Any> UnivariateIntegrand(
     function: (Double) -> T,
     vararg features: IntegrandFeature,
-): UnivariateIntegrand<T> = UnivariateIntegrand(features.associateBy { it::class }, function)
+): UnivariateIntegrand<T> = UnivariateIntegrand(FeatureSet.of(*features), function)
 
 public typealias UnivariateIntegrator<T> = Integrator<UnivariateIntegrand<T>>
 
@@ -79,7 +70,7 @@ public val <T : Any> UnivariateIntegrand<T>.value: T get() = valueOrNull ?: erro
 public fun <T : Any> UnivariateIntegrator<T>.integrate(
     vararg features: IntegrandFeature,
     function: (Double) -> T,
-): UnivariateIntegrand<T> = integrate(UnivariateIntegrand(function, *features))
+): UnivariateIntegrand<T> = process(UnivariateIntegrand(function, *features))
 
 /**
  * A shortcut method to integrate a [function] in [range] with additional [features].
@@ -90,7 +81,7 @@ public fun <T : Any> UnivariateIntegrator<T>.integrate(
     range: ClosedRange<Double>,
     vararg features: IntegrandFeature,
     function: (Double) -> T,
-): UnivariateIntegrand<T> = integrate(UnivariateIntegrand(function, IntegrationRange(range), *features))
+): UnivariateIntegrand<T> = process(UnivariateIntegrand(function, IntegrationRange(range), *features))
 
 /**
  * A shortcut method to integrate a [function] in [range] with additional features.
@@ -107,5 +98,5 @@ public fun <T : Any> UnivariateIntegrator<T>.integrate(
         featureBuilder()
         add(IntegrationRange(range))
     }
-    return integrate(UnivariateIntegrand(function, *features.toTypedArray()))
+    return process(UnivariateIntegrand(function, *features.toTypedArray()))
 }

kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/Interpolator.kt

@@ -42,20 +42,20 @@ public fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolatePolynomials(
     x: Buffer<T>,
     y: Buffer<T>,
 ): PiecewisePolynomial<T> {
-    val pointSet = XYColumnarData(x, y)
+    val pointSet = XYColumnarData.of(x, y)
     return interpolatePolynomials(pointSet)
 }
 
 public fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolatePolynomials(
     data: Map<T, T>,
 ): PiecewisePolynomial<T> {
-    val pointSet = XYColumnarData(data.keys.toList().asBuffer(), data.values.toList().asBuffer())
+    val pointSet = XYColumnarData.of(data.keys.toList().asBuffer(), data.values.toList().asBuffer())
     return interpolatePolynomials(pointSet)
 }
 
 public fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolatePolynomials(
     data: List<Pair<T, T>>,
 ): PiecewisePolynomial<T> {
-    val pointSet = XYColumnarData(data.map { it.first }.asBuffer(), data.map { it.second }.asBuffer())
+    val pointSet = XYColumnarData.of(data.map { it.first }.asBuffer(), data.map { it.second }.asBuffer())
     return interpolatePolynomials(pointSet)
 }

kmath-kotlingrad/src/main/kotlin/space/kscience/kmath/kotlingrad/KotlingradExpression.kt

@@ -27,15 +27,26 @@ public class KotlingradExpression<T : Number, A : NumericAlgebra<T>>(
 ) : SpecialDifferentiableExpression<T, KotlingradExpression<T, A>> {
     override fun invoke(arguments: Map<Symbol, T>): T = mst.interpret(algebra, arguments)
 
-    override fun derivativeOrNull(symbols: List<Symbol>): KotlingradExpression<T, A> =
+    override fun derivativeOrNull(
-        KotlingradExpression(
+        symbols: List<Symbol>,
-            algebra,
+    ): KotlingradExpression<T, A> = KotlingradExpression(
-            symbols.map(Symbol::identity)
+        algebra,
-                .map(MstNumericAlgebra::bindSymbol)
+        symbols.map(Symbol::identity)
-                .map<Symbol, SVar<KMathNumber<T, A>>>(Symbol::toSVar)
+            .map(MstNumericAlgebra::bindSymbol)
-                .fold(mst.toSFun(), SFun<KMathNumber<T, A>>::d)
+            .map<Symbol, SVar<KMathNumber<T, A>>>(Symbol::toSVar)
-                .toMst(),
+            .fold(mst.toSFun(), SFun<KMathNumber<T, A>>::d)
-        )
+            .toMst(),
+    )
+}
+
+/**
+ * A diff processor using [MST] to Kotlingrad converter
+ */
+public class KotlingradProcessor<T : Number, A : NumericAlgebra<T>>(
+    public val algebra: A,
+) : AutoDiffProcessor<T, MST, MstExtendedField> {
+    override fun differentiate(function: MstExtendedField.() -> MST): DifferentiableExpression<T> =
+        MstExtendedField.function().toKotlingradExpression(algebra)
 }
 
 /**

kmath-optimization/build.gradle.kts

@@ -0,0 +1,20 @@
+plugins {
+    id("ru.mipt.npm.gradle.mpp")
+    id("ru.mipt.npm.gradle.native")
+}
+
+kscience {
+    useAtomic()
+}
+
+kotlin.sourceSets {
+    commonMain {
+        dependencies {
+            api(project(":kmath-coroutines"))
+        }
+    }
+}
+
+readme {
+    maturity = ru.mipt.npm.gradle.Maturity.EXPERIMENTAL
+}

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

@@ -0,0 +1,71 @@
+/*
+ * Copyright 2018-2021 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.optimization
+
+import space.kscience.kmath.expressions.DifferentiableExpression
+import space.kscience.kmath.expressions.Symbol
+import space.kscience.kmath.misc.FeatureSet
+
+public class OptimizationValue<T>(public val value: T) : OptimizationFeature {
+    override fun toString(): String = "Value($value)"
+}
+
+public enum class FunctionOptimizationTarget : OptimizationFeature {
+    MAXIMIZE,
+    MINIMIZE
+}
+
+public class FunctionOptimization<T>(
+    override val features: FeatureSet<OptimizationFeature>,
+    public val expression: DifferentiableExpression<T>,
+) : OptimizationProblem<T> {
+
+
+    override fun equals(other: Any?): Boolean {
+        if (this === other) return true
+        if (other == null || this::class != other::class) return false
+
+        other as FunctionOptimization<*>
+
+        if (features != other.features) return false
+        if (expression != other.expression) return false
+
+        return true
+    }
+
+    override fun hashCode(): Int {
+        var result = features.hashCode()
+        result = 31 * result + expression.hashCode()
+        return result
+    }
+
+    override fun toString(): String = "FunctionOptimization(features=$features)"
+}
+
+public fun <T> FunctionOptimization<T>.withFeatures(
+    vararg newFeature: OptimizationFeature,
+): FunctionOptimization<T> = FunctionOptimization(
+    features.with(*newFeature),
+    expression,
+)
+
+/**
+ * Optimizes differentiable expression using specific [optimizer] form given [startingPoint].
+ */
+public suspend fun <T : Any> DifferentiableExpression<T>.optimizeWith(
+    optimizer: Optimizer<T, FunctionOptimization<T>>,
+    startingPoint: Map<Symbol, T>,
+    vararg features: OptimizationFeature,
+): FunctionOptimization<T> {
+    val problem = FunctionOptimization<T>(FeatureSet.of(OptimizationStartPoint(startingPoint), *features), this)
+    return optimizer.optimize(problem)
+}
+
+public val <T> FunctionOptimization<T>.resultValueOrNull: T?
+    get() = getFeature<OptimizationResult<T>>()?.point?.let { expression(it) }
+
+public val <T> FunctionOptimization<T>.resultValue: T
+    get() = resultValueOrNull ?: error("Result is not present in $this")

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

@@ -0,0 +1,93 @@
+/*
+ * Copyright 2018-2021 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.optimization
+
+import space.kscience.kmath.data.XYColumnarData
+import space.kscience.kmath.expressions.DifferentiableExpression
+import space.kscience.kmath.expressions.Symbol
+import space.kscience.kmath.misc.FeatureSet
+
+public abstract class OptimizationBuilder<T, R : OptimizationProblem<T>> {
+    public val features: MutableList<OptimizationFeature> = ArrayList()
+
+    public fun addFeature(feature: OptimizationFeature) {
+        features.add(feature)
+    }
+
+    public inline fun <reified T : OptimizationFeature> updateFeature(update: (T?) -> T) {
+        val existing = features.find { it.key == T::class } as? T
+        val new = update(existing)
+        if (existing != null) {
+            features.remove(existing)
+        }
+        addFeature(new)
+    }
+
+    public abstract fun build(): R
+}
+
+public fun <T> OptimizationBuilder<T, *>.startAt(startingPoint: Map<Symbol, T>) {
+    addFeature(OptimizationStartPoint(startingPoint))
+}
+
+public class FunctionOptimizationBuilder<T>(
+    private val expression: DifferentiableExpression<T>,
+) : OptimizationBuilder<T, FunctionOptimization<T>>() {
+    override fun build(): FunctionOptimization<T> = FunctionOptimization(FeatureSet.of(features), expression)
+}
+
+public fun <T> FunctionOptimization(
+    expression: DifferentiableExpression<T>,
+    builder: FunctionOptimizationBuilder<T>.() -> Unit,
+): FunctionOptimization<T> = FunctionOptimizationBuilder(expression).apply(builder).build()
+
+public suspend fun <T> DifferentiableExpression<T>.optimizeWith(
+    optimizer: Optimizer<T, FunctionOptimization<T>>,
+    startingPoint: Map<Symbol, T>,
+    builder: FunctionOptimizationBuilder<T>.() -> Unit = {},
+): FunctionOptimization<T> {
+    val problem = FunctionOptimization<T>(this) {
+        startAt(startingPoint)
+        builder()
+    }
+    return optimizer.optimize(problem)
+}
+
+public suspend fun <T> DifferentiableExpression<T>.optimizeWith(
+    optimizer: Optimizer<T, FunctionOptimization<T>>,
+    vararg startingPoint: Pair<Symbol, T>,
+    builder: FunctionOptimizationBuilder<T>.() -> Unit = {},
+): FunctionOptimization<T> {
+    val problem = FunctionOptimization<T>(this) {
+        startAt(mapOf(*startingPoint))
+        builder()
+    }
+    return optimizer.optimize(problem)
+}
+
+
+public class XYOptimizationBuilder(
+    public val data: XYColumnarData<Double, Double, Double>,
+    public val model: DifferentiableExpression<Double>,
+) : OptimizationBuilder<Double, XYFit>() {
+
+    public var pointToCurveDistance: PointToCurveDistance = PointToCurveDistance.byY
+    public var pointWeight: PointWeight = PointWeight.byYSigma
+
+    override fun build(): XYFit = XYFit(
+        data,
+        model,
+        FeatureSet.of(features),
+        pointToCurveDistance,
+        pointWeight
+    )
+}
+
+public fun XYOptimization(
+    data: XYColumnarData<Double, Double, Double>,
+    model: DifferentiableExpression<Double>,
+    builder: XYOptimizationBuilder.() -> Unit,
+): XYFit = XYOptimizationBuilder(data, model).apply(builder).build()

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

@@ -0,0 +1,66 @@
+/*
+ * Copyright 2018-2021 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.optimization
+
+import space.kscience.kmath.expressions.DifferentiableExpression
+import space.kscience.kmath.expressions.Symbol
+import space.kscience.kmath.linear.Matrix
+import space.kscience.kmath.misc.*
+import kotlin.reflect.KClass
+
+public interface OptimizationFeature : Feature<OptimizationFeature> {
+    // enforce toString override
+    override fun toString(): String
+}
+
+public interface OptimizationProblem<T> : Featured<OptimizationFeature> {
+    public val features: FeatureSet<OptimizationFeature>
+    override fun <F : OptimizationFeature> getFeature(type: KClass<out F>): F? = features.getFeature(type)
+}
+
+public inline fun <reified F : OptimizationFeature> OptimizationProblem<*>.getFeature(): F? = getFeature(F::class)
+
+public open class OptimizationStartPoint<T>(public val point: Map<Symbol, T>) : OptimizationFeature {
+    override fun toString(): String = "StartPoint($point)"
+}
+
+
+public interface OptimizationPrior<T> : OptimizationFeature, DifferentiableExpression<T> {
+    override val key: FeatureKey<OptimizationFeature> get() = OptimizationPrior::class
+}
+
+public class OptimizationCovariance<T>(public val covariance: Matrix<T>) : OptimizationFeature {
+    override fun toString(): String = "Covariance($covariance)"
+}
+
+/**
+ * Get the starting point for optimization. Throws error if not defined.
+ */
+public val <T> OptimizationProblem<T>.startPoint: Map<Symbol, T>
+    get() = getFeature<OptimizationStartPoint<T>>()?.point
+        ?: error("Starting point not defined in $this")
+
+public open class OptimizationResult<T>(public val point: Map<Symbol, T>) : OptimizationFeature {
+    override fun toString(): String = "Result($point)"
+}
+
+public val <T> OptimizationProblem<T>.resultPointOrNull: Map<Symbol, T>?
+    get() = getFeature<OptimizationResult<T>>()?.point
+
+public val <T> OptimizationProblem<T>.resultPoint: Map<Symbol, T>
+    get() = resultPointOrNull ?: error("Result is not present in $this")
+
+public class OptimizationLog(private val loggable: Loggable) : Loggable by loggable, OptimizationFeature {
+    override fun toString(): String = "Log($loggable)"
+}
+
+public class OptimizationParameters(public val symbols: List<Symbol>) : OptimizationFeature {
+    public constructor(vararg symbols: Symbol) : this(listOf(*symbols))
+
+    override fun toString(): String = "Parameters($symbols)"
+}
+
+

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

@@ -0,0 +1,10 @@
+/*
+ * Copyright 2018-2021 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.optimization
+
+public interface Optimizer<T, P : OptimizationProblem<T>> {
+    public suspend fun optimize(problem: P): P
+}

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

@@ -0,0 +1,250 @@
+/*
+ * Copyright 2018-2021 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.optimization
+
+import space.kscience.kmath.expressions.DifferentiableExpression
+import space.kscience.kmath.expressions.Symbol
+import space.kscience.kmath.expressions.SymbolIndexer
+import space.kscience.kmath.expressions.derivative
+import space.kscience.kmath.linear.*
+import space.kscience.kmath.misc.UnstableKMathAPI
+import space.kscience.kmath.misc.log
+import space.kscience.kmath.operations.DoubleField
+import space.kscience.kmath.structures.DoubleBuffer
+import space.kscience.kmath.structures.DoubleL2Norm
+
+
+/**
+ * An optimizer based onf Fyodor Tkachev's quasi-optimal weights method.
+ * See [the article](http://arxiv.org/abs/physics/0604127).
+ */
+@UnstableKMathAPI
+public object QowOptimizer : Optimizer<Double, XYFit> {
+
+    private val linearSpace: LinearSpace<Double, DoubleField> = LinearSpace.double
+    private val solver: LinearSolver<Double> = linearSpace.lupSolver()
+
+    @OptIn(UnstableKMathAPI::class)
+    private class QoWeight(
+        val problem: XYFit,
+        val parameters: Map<Symbol, Double>,
+    ) : Map<Symbol, Double> by parameters, SymbolIndexer {
+        override val symbols: List<Symbol> = parameters.keys.toList()
+
+        val data get() = problem.data
+
+        /**
+         * Derivatives of the spectrum over parameters. First index in the point number, second one - index of parameter
+         */
+        val derivs: Matrix<Double> by lazy {
+            linearSpace.buildMatrix(problem.data.size, symbols.size) { d, s ->
+                problem.distance(d).derivative(symbols[s])(parameters)
+            }
+        }
+
+        /**
+         * Array of dispersions in each point
+         */
+        val dispersion: Point<Double> by lazy {
+            DoubleBuffer(problem.data.size) { d ->
+                1.0/problem.weight(d).invoke(parameters)
+            }
+        }
+
+        val prior: DifferentiableExpression<Double>? get() = problem.getFeature<OptimizationPrior<Double>>()
+
+        override fun toString(): String  = parameters.toString()
+    }
+
+    /**
+     * The signed distance from the model to the [d]-th point of data.
+     */
+    private fun QoWeight.distance(d: Int, parameters: Map<Symbol, Double>): Double = problem.distance(d)(parameters)
+
+
+    /**
+     * The derivative of [distance]
+     */
+    private fun QoWeight.distanceDerivative(symbol: Symbol, d: Int, parameters: Map<Symbol, Double>): Double =
+        problem.distance(d).derivative(symbol)(parameters)
+
+    /**
+     * Теоретическая ковариация весовых функций.
+     *
+     * D(\phi)=E(\phi_k(\theta_0) \phi_l(\theta_0))= disDeriv_k * disDeriv_l /sigma^2
+     */
+    private fun QoWeight.covarF(): Matrix<Double> =
+        linearSpace.matrix(size, size).symmetric { s1, s2 ->
+            (0 until data.size).sumOf { d -> derivs[d, s1] * derivs[d, s2] / dispersion[d] }
+        }
+
+    /**
+     * Экспериментальная ковариация весов. Формула (22) из
+     * http://arxiv.org/abs/physics/0604127
+     */
+    private fun QoWeight.covarFExp(theta: Map<Symbol, Double>): Matrix<Double> =
+        with(linearSpace) {
+            /*
+             * Важно! Если не делать предварителього вычисления этих производных, то
+             * количество вызывов функции будет dim^2 вместо dim Первый индекс -
+             * номер точки, второй - номер переменной, по которой берется производная
+             */
+            val eqvalues = linearSpace.buildMatrix(data.size, size) { d, s ->
+                distance(d, theta) * derivs[d, s] / dispersion[d]
+            }
+
+            buildMatrix(size, size) { s1, s2 ->
+                (0 until data.size).sumOf { d -> eqvalues[d, s2] * eqvalues[d, s1] }
+            }
+        }
+
+    /**
+     * Equation derivatives for Newton run
+     */
+    private fun QoWeight.getEqDerivValues(
+        theta: Map<Symbol, Double> = parameters,
+    ): Matrix<Double> = with(linearSpace) {
+        //Возвращает производную k-того Eq по l-тому параметру
+        //val res = Array(fitDim) { DoubleArray(fitDim) }
+        val sderiv = buildMatrix(data.size, size) { d, s ->
+            distanceDerivative(symbols[s], d, theta)
+        }
+
+        buildMatrix(size, size) { s1, s2 ->
+            val base = (0 until data.size).sumOf { d ->
+                require(dispersion[d] > 0)
+                sderiv[d, s2] * derivs[d, s1] / dispersion[d]
+            }
+            prior?.let { prior ->
+                //Check if this one is correct
+                val pi = prior(theta)
+                val deriv1 = prior.derivative(symbols[s1])(theta)
+                val deriv2 = prior.derivative(symbols[s2])(theta)
+                base + deriv1 * deriv2 / pi / pi
+            } ?: base
+        }
+    }
+
+
+    /**
+     * Значения уравнений метода квазиоптимальных весов
+     */
+    private fun QoWeight.getEqValues(theta: Map<Symbol, Double> = this): Point<Double> {
+        val distances = DoubleBuffer(data.size) { d -> distance(d, theta) }
+
+        return DoubleBuffer(size) { s ->
+            val base = (0 until data.size).sumOf { d -> distances[d] * derivs[d, s] / dispersion[d] }
+            //Поправка на априорную вероятность
+            prior?.let { prior ->
+                base - prior.derivative(symbols[s])(theta) / prior(theta)
+            } ?: base
+        }
+    }
+
+
+    private fun QoWeight.newtonianStep(
+        theta: Map<Symbol, Double>,
+        eqvalues: Point<Double>,
+    ): QoWeight = linearSpace {
+        with(this@newtonianStep) {
+            val start = theta.toPoint()
+            val invJacob = solver.inverse(this@newtonianStep.getEqDerivValues(theta))
+
+            val step = invJacob.dot(eqvalues)
+            return QoWeight(problem, theta + (start - step).toMap())
+        }
+    }
+
+    private fun QoWeight.newtonianRun(
+        maxSteps: Int = 100,
+        tolerance: Double = 0.0,
+        fast: Boolean = false,
+    ): QoWeight {
+
+        val logger = problem.getFeature<OptimizationLog>()
+
+        var dis: Double //discrepancy value
+        // Working with the full set of parameters
+        var par = problem.startPoint
+
+        logger?.log { "Starting newtonian iteration from: \n\t$par" }
+
+        var eqvalues = getEqValues(par) //Values of the weight functions
+
+        dis = DoubleL2Norm.norm(eqvalues) // discrepancy
+        logger?.log { "Starting discrepancy is $dis" }
+        var i = 0
+        var flag = false
+        while (!flag) {
+            i++
+            logger?.log { "Starting step number $i" }
+
+            val currentSolution = if (fast) {
+                //Берет значения матрицы в той точке, где считается вес
+                newtonianStep(this, eqvalues)
+            } else {
+                //Берет значения матрицы в точке par
+                newtonianStep(par, eqvalues)
+            }
+            // здесь должен стоять учет границ параметров
+            logger?.log { "Parameter values after step are: \n\t$currentSolution" }
+
+            eqvalues = getEqValues(currentSolution)
+            val currentDis = DoubleL2Norm.norm(eqvalues)// невязка после шага
+
+            logger?.log { "The discrepancy after step is: $currentDis." }
+
+            if (currentDis >= dis && i > 1) {
+                //дополнительно проверяем, чтобы был сделан хотя бы один шаг
+                flag = true
+                logger?.log { "The discrepancy does not decrease. Stopping iteration." }
+            } else {
+                par = currentSolution
+                dis = currentDis
+            }
+            if (i >= maxSteps) {
+                flag = true
+                logger?.log { "Maximum number of iterations reached. Stopping iteration." }
+            }
+            if (dis <= tolerance) {
+                flag = true
+                logger?.log { "Tolerance threshold is reached. Stopping iteration." }
+            }
+        }
+
+        return QoWeight(problem, par)
+    }
+
+    private fun QoWeight.covariance(): Matrix<Double> {
+        val logger = problem.getFeature<OptimizationLog>()
+
+        logger?.log {
+            """
+            Starting errors estimation using quasioptimal weights method. The starting weight is:
+                ${problem.startPoint}
+             """.trimIndent()
+        }
+
+        val covar = solver.inverse(getEqDerivValues())
+        //TODO fix eigenvalues check
+//        val decomposition = EigenDecomposition(covar.matrix)
+//        var valid = true
+//        for (lambda in decomposition.realEigenvalues) {
+//            if (lambda <= 0) {
+//                logger?.log { "The covariance matrix is not positive defined. Error estimation is not valid" }
+//                valid = false
+//            }
+//        }
+        return covar
+    }
+
+    override suspend fun optimize(problem: XYFit): XYFit {
+        val qowSteps = 2
+        val initialWeight = QoWeight(problem, problem.startPoint)
+        val res = initialWeight.newtonianRun()
+        return res.problem.withFeature(OptimizationResult(res.parameters))
+    }
+}

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

@@ -0,0 +1,146 @@
+/*
+ * Copyright 2018-2021 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.
+ */
+@file:OptIn(UnstableKMathAPI::class)
+
+package space.kscience.kmath.optimization
+
+import space.kscience.kmath.data.XYColumnarData
+import space.kscience.kmath.data.indices
+import space.kscience.kmath.expressions.*
+import space.kscience.kmath.misc.FeatureSet
+import space.kscience.kmath.misc.Loggable
+import space.kscience.kmath.misc.UnstableKMathAPI
+import space.kscience.kmath.operations.ExtendedField
+import space.kscience.kmath.operations.bindSymbol
+import kotlin.math.pow
+
+/**
+ * Specify the way to compute distance from point to the curve as DifferentiableExpression
+ */
+public interface PointToCurveDistance : OptimizationFeature {
+    public fun distance(problem: XYFit, index: Int): DifferentiableExpression<Double>
+
+    public companion object {
+        public val byY: PointToCurveDistance = object : PointToCurveDistance {
+            override fun distance(problem: XYFit, index: Int): DifferentiableExpression<Double> {
+                val x = problem.data.x[index]
+                val y = problem.data.y[index]
+
+                return object : DifferentiableExpression<Double> {
+                    override fun derivativeOrNull(
+                        symbols: List<Symbol>
+                    ): Expression<Double>? = problem.model.derivativeOrNull(symbols)?.let { derivExpression ->
+                        Expression { arguments ->
+                            derivExpression.invoke(arguments + (Symbol.x to x))
+                        }
+                    }
+
+                    override fun invoke(arguments: Map<Symbol, Double>): Double =
+                        problem.model(arguments + (Symbol.x to x)) - y
+                }
+            }
+
+            override fun toString(): String = "PointToCurveDistanceByY"
+        }
+    }
+}
+
+/**
+ * Compute a wight of the point. The more the weight, the more impact this point will have on the fit.
+ * By default, uses Dispersion^-1
+ */
+public interface PointWeight : OptimizationFeature {
+    public fun weight(problem: XYFit, index: Int): DifferentiableExpression<Double>
+
+    public companion object {
+        public fun bySigma(sigmaSymbol: Symbol): PointWeight = object : PointWeight {
+            override fun weight(problem: XYFit, index: Int): DifferentiableExpression<Double> =
+                object : DifferentiableExpression<Double> {
+                    override fun invoke(arguments: Map<Symbol, Double>): Double {
+                        return problem.data[sigmaSymbol]?.get(index)?.pow(-2) ?: 1.0
+                    }
+
+                    override fun derivativeOrNull(symbols: List<Symbol>): Expression<Double> = Expression { 0.0 }
+                }
+
+            override fun toString(): String = "PointWeightBySigma($sigmaSymbol)"
+
+        }
+
+        public val byYSigma: PointWeight = bySigma(Symbol.yError)
+    }
+}
+
+/**
+ * A fit problem for X-Y-Yerr data. Also known as "least-squares" problem.
+ */
+public class XYFit(
+    public val data: XYColumnarData<Double, Double, Double>,
+    public val model: DifferentiableExpression<Double>,
+    override val features: FeatureSet<OptimizationFeature>,
+    internal val pointToCurveDistance: PointToCurveDistance = PointToCurveDistance.byY,
+    internal val pointWeight: PointWeight = PointWeight.byYSigma,
+    public val xSymbol: Symbol = Symbol.x,
+) : OptimizationProblem<Double> {
+    public fun distance(index: Int): DifferentiableExpression<Double> = pointToCurveDistance.distance(this, index)
+
+    public fun weight(index: Int): DifferentiableExpression<Double> = pointWeight.weight(this, index)
+}
+
+public fun XYFit.withFeature(vararg features: OptimizationFeature): XYFit {
+    return XYFit(data, model, this.features.with(*features), pointToCurveDistance, pointWeight)
+}
+
+/**
+ * Fit given dta with
+ */
+public suspend fun <I : Any, A> XYColumnarData<Double, Double, Double>.fitWith(
+    optimizer: Optimizer<Double, XYFit>,
+    processor: AutoDiffProcessor<Double, I, A>,
+    startingPoint: Map<Symbol, Double>,
+    vararg features: OptimizationFeature = emptyArray(),
+    xSymbol: Symbol = Symbol.x,
+    pointToCurveDistance: PointToCurveDistance = PointToCurveDistance.byY,
+    pointWeight: PointWeight = PointWeight.byYSigma,
+    model: A.(I) -> I
+): XYFit where A : ExtendedField<I>, A : ExpressionAlgebra<Double, I> {
+    val modelExpression = processor.differentiate {
+        val x = bindSymbol(xSymbol)
+        model(x)
+    }
+
+    var actualFeatures = FeatureSet.of(*features, OptimizationStartPoint(startingPoint))
+
+    if (actualFeatures.getFeature<OptimizationLog>() == null) {
+        actualFeatures = actualFeatures.with(OptimizationLog(Loggable.console))
+    }
+    val problem = XYFit(
+        this,
+        modelExpression,
+        actualFeatures,
+        pointToCurveDistance,
+        pointWeight,
+        xSymbol
+    )
+    return optimizer.optimize(problem)
+}
+
+/**
+ * Compute chi squared value for completed fit. Return null for incomplete fit
+ */
+public val XYFit.chiSquaredOrNull: Double? get() {
+    val result = resultPointOrNull ?: return null
+
+    return data.indices.sumOf { index->
+
+        val x = data.x[index]
+        val y = data.y[index]
+        val yErr = data[Symbol.yError]?.get(index) ?: 1.0
+
+        val mu = model.invoke(result + (xSymbol to x) )
+
+        ((y - mu)/yErr).pow(2)
+    }
+}

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

@@ -0,0 +1,66 @@
+/*
+ * Copyright 2018-2021 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.optimization
+
+import space.kscience.kmath.data.XYColumnarData
+import space.kscience.kmath.data.indices
+import space.kscience.kmath.expressions.DifferentiableExpression
+import space.kscience.kmath.expressions.Expression
+import space.kscience.kmath.expressions.Symbol
+import space.kscience.kmath.expressions.derivative
+import space.kscience.kmath.misc.UnstableKMathAPI
+import kotlin.math.PI
+import kotlin.math.ln
+import kotlin.math.pow
+import kotlin.math.sqrt
+
+
+private val oneOver2Pi = 1.0 / sqrt(2 * PI)
+
+@UnstableKMathAPI
+internal fun XYFit.logLikelihood(): DifferentiableExpression<Double> = object : DifferentiableExpression<Double> {
+    override fun derivativeOrNull(symbols: List<Symbol>): Expression<Double> = Expression { arguments ->
+        data.indices.sumOf { index ->
+            val d = distance(index)(arguments)
+            val weight = weight(index)(arguments)
+            val weightDerivative = weight(index).derivative(symbols)(arguments)
+
+            //  -1 / (sqrt(2 PI) * sigma) + 2 (x-mu)/ 2 sigma^2 * d mu/ d theta - (x-mu)^2 / 2 * d w/ d theta
+            return@sumOf -oneOver2Pi * sqrt(weight) + //offset derivative
+                    d * model.derivative(symbols)(arguments) * weight - //model derivative
+                    d.pow(2) * weightDerivative / 2 //weight derivative
+        }
+    }
+
+    override fun invoke(arguments: Map<Symbol, Double>): Double {
+        return data.indices.sumOf { index ->
+            val d = distance(index)(arguments)
+            val weight = weight(index)(arguments)
+            //1/sqrt(2 PI sigma^2) - (x-mu)^2/ (2 * sigma^2)
+            oneOver2Pi * ln(weight) - d.pow(2) * weight
+        } / 2
+    }
+
+}
+
+/**
+ * Optimize given XY (least squares) [problem] using this function [Optimizer].
+ * The problem is treated as maximum likelihood problem and is done via maximizing logarithmic likelihood, respecting
+ * possible weight dependency on the model and parameters.
+ */
+@UnstableKMathAPI
+public suspend fun Optimizer<Double, FunctionOptimization<Double>>.maximumLogLikelihood(problem: XYFit): XYFit {
+    val functionOptimization = FunctionOptimization(problem.features, problem.logLikelihood())
+    val result = optimize(functionOptimization.withFeatures(FunctionOptimizationTarget.MAXIMIZE))
+    return XYFit(problem.data, problem.model, result.features)
+}
+
+@UnstableKMathAPI
+public suspend fun Optimizer<Double, FunctionOptimization<Double>>.maximumLogLikelihood(
+    data: XYColumnarData<Double, Double, Double>,
+    model: DifferentiableExpression<Double>,
+    builder: XYOptimizationBuilder.() -> Unit,
+): XYFit = maximumLogLikelihood(XYOptimization(data, model, builder))

kmath-stat/build.gradle.kts

@@ -1,6 +1,5 @@
 plugins {
-    kotlin("multiplatform")
+    id("ru.mipt.npm.gradle.mpp")
-    id("ru.mipt.npm.gradle.common")
     id("ru.mipt.npm.gradle.native")
 }
 

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

@@ -1,91 +0,0 @@
-/*
- * Copyright 2018-2021 KMath contributors.
- * Use of this source code is governed by the Apache 2.0 license that can be found in the LICENSE file.
- */
-
-package space.kscience.kmath.optimization
-
-import space.kscience.kmath.expressions.*
-import space.kscience.kmath.operations.ExtendedField
-import space.kscience.kmath.structures.Buffer
-import space.kscience.kmath.structures.indices
-
-/**
- * A likelihood function optimization problem with provided derivatives.
- */
-public interface FunctionOptimization<T : Any> : Optimization<T> {
-    /**
-     * The optimization direction. If true search for function maximum, if false, search for the minimum.
-     */
-    public var maximize: Boolean
-
-    /**
-     * Defines the initial guess for the optimization problem.
-     */
-    public fun initialGuess(map: Map<Symbol, T>)
-
-    /**
-     * Set a differentiable expression as objective function as function and gradient provider.
-     */
-    public fun diffFunction(expression: DifferentiableExpression<T>)
-
-    public companion object {
-        /**
-         * Generate a chi squared expression from given x-y-sigma data and inline model. Provides automatic
-         * differentiation.
-         */
-        public fun <T : Any, I : Any, A> chiSquared(
-            autoDiff: AutoDiffProcessor<T, I, A, Expression<T>>,
-            x: Buffer<T>,
-            y: Buffer<T>,
-            yErr: Buffer<T>,
-            model: A.(I) -> I,
-        ): DifferentiableExpression<T> where A : ExtendedField<I>, A : ExpressionAlgebra<T, I> {
-            require(x.size == y.size) { "X and y buffers should be of the same size" }
-            require(y.size == yErr.size) { "Y and yErr buffer should of the same size" }
-
-            return autoDiff.process {
-                var sum = zero
-
-                x.indices.forEach {
-                    val xValue = const(x[it])
-                    val yValue = const(y[it])
-                    val yErrValue = const(yErr[it])
-                    val modelValue = model(xValue)
-                    sum += ((yValue - modelValue) / yErrValue).pow(2)
-                }
-
-                sum
-            }
-        }
-    }
-}
-
-/**
- * Defines a chi-squared-based objective function.
- */
-public fun <T: Any, I : Any, A> FunctionOptimization<T>.chiSquared(
-    autoDiff: AutoDiffProcessor<T, I, A, Expression<T>>,
-    x: Buffer<T>,
-    y: Buffer<T>,
-    yErr: Buffer<T>,
-    model: A.(I) -> I,
-) where A : ExtendedField<I>, A : ExpressionAlgebra<T, I> {
-    val chiSquared = FunctionOptimization.chiSquared(autoDiff, x, y, yErr, model)
-    diffFunction(chiSquared)
-    maximize = false
-}
-
-/**
- * Optimizes differentiable expression using specific [OptimizationProblemFactory].
- */
-public fun <T : Any, F : FunctionOptimization<T>> DifferentiableExpression<T>.optimizeWith(
-    factory: OptimizationProblemFactory<T, F>,
-    vararg symbols: Symbol,
-    configuration: F.() -> Unit,
-): OptimizationResult<T> {
-    require(symbols.isNotEmpty()) { "Must provide a list of symbols for optimization" }
-    val problem = factory(symbols.toList(), configuration)
-    problem.diffFunction(this)
-    return problem.optimize()
-}

kmath-stat/src/commonMain/kotlin/space/kscience/kmath/optimization/NoDerivFunctionOptimization.kt

@@ -1,74 +0,0 @@
-/*
- * Copyright 2018-2021 KMath contributors.
- * Use of this source code is governed by the Apache 2.0 license that can be found in the LICENSE file.
- */
-
-package space.kscience.kmath.optimization
-
-import space.kscience.kmath.expressions.Expression
-import space.kscience.kmath.expressions.Symbol
-import space.kscience.kmath.structures.Buffer
-import space.kscience.kmath.structures.indices
-import kotlin.math.pow
-
-/**
- * A likelihood function optimization problem
- */
-public interface NoDerivFunctionOptimization<T : Any> : Optimization<T> {
-    /**
-     * The optimization direction. If `true` search for function maximum, search for the minimum otherwise.
-     */
-    public var maximize: Boolean
-
-    /**
-     * Define the initial guess for the optimization problem
-     */
-    public fun initialGuess(map: Map<Symbol, T>)
-
-    /**
-     * Set an objective function expression
-     */
-    public fun function(expression: Expression<T>)
-
-    public companion object {
-        /**
-         * Generate a chi squared expression from given x-y-sigma model represented by an expression. Does not provide derivatives
-         */
-        public fun chiSquared(
-            x: Buffer<Double>,
-            y: Buffer<Double>,
-            yErr: Buffer<Double>,
-            model: Expression<Double>,
-            xSymbol: Symbol = Symbol.x,
-        ): Expression<Double> {
-            require(x.size == y.size) { "X and y buffers should be of the same size" }
-            require(y.size == yErr.size) { "Y and yErr buffer should of the same size" }
-
-            return Expression { arguments ->
-                x.indices.sumOf {
-                    val xValue = x[it]
-                    val yValue = y[it]
-                    val yErrValue = yErr[it]
-                    val modifiedArgs = arguments + (xSymbol to xValue)
-                    val modelValue = model(modifiedArgs)
-                    ((yValue - modelValue) / yErrValue).pow(2)
-                }
-            }
-        }
-    }
-}
-
-
-/**
- * Optimize expression without derivatives using specific [OptimizationProblemFactory]
- */
-public fun <T : Any, F : NoDerivFunctionOptimization<T>> Expression<T>.noDerivOptimizeWith(
-    factory: OptimizationProblemFactory<T, F>,
-    vararg symbols: Symbol,
-    configuration: F.() -> Unit,
-): OptimizationResult<T> {
-    require(symbols.isNotEmpty()) { "Must provide a list of symbols for optimization" }
-    val problem = factory(symbols.toList(), configuration)
-    problem.function(this)
-    return problem.optimize()
-}

kmath-stat/src/commonMain/kotlin/space/kscience/kmath/optimization/Optimization.kt

@@ -1,49 +0,0 @@
-/*
- * Copyright 2018-2021 KMath contributors.
- * Use of this source code is governed by the Apache 2.0 license that can be found in the LICENSE file.
- */
-
-package space.kscience.kmath.optimization
-
-import space.kscience.kmath.expressions.Symbol
-
-public interface OptimizationFeature
-
-public class OptimizationResult<out T>(
-    public val point: Map<Symbol, T>,
-    public val value: T,
-    public val features: Set<OptimizationFeature> = emptySet(),
-) {
-    override fun toString(): String {
-        return "OptimizationResult(point=$point, value=$value)"
-    }
-}
-
-public operator fun <T> OptimizationResult<T>.plus(
-    feature: OptimizationFeature,
-): OptimizationResult<T> = OptimizationResult(point, value, features + feature)
-
-/**
- * A builder of optimization problems over [T] variables
- */
-public interface Optimization<T : Any> {
-
-    /**
-     * Update the problem from previous optimization run
-     */
-    public fun update(result: OptimizationResult<T>)
-
-    /**
-     * Make an optimization run
-     */
-    public fun optimize(): OptimizationResult<T>
-}
-
-public fun interface OptimizationProblemFactory<T : Any, out P : Optimization<T>> {
-    public fun build(symbols: List<Symbol>): P
-}
-
-public operator fun <T : Any, P : Optimization<T>> OptimizationProblemFactory<T, P>.invoke(
-    symbols: List<Symbol>,
-    block: P.() -> Unit,
-): P = build(symbols).apply(block)

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

@@ -1,41 +0,0 @@
-/*
- * Copyright 2018-2021 KMath contributors.
- * Use of this source code is governed by the Apache 2.0 license that can be found in the LICENSE file.
- */
-
-package space.kscience.kmath.optimization
-
-import space.kscience.kmath.data.ColumnarData
-import space.kscience.kmath.expressions.*
-import space.kscience.kmath.misc.UnstableKMathAPI
-import space.kscience.kmath.operations.ExtendedField
-import space.kscience.kmath.operations.Field
-
-@UnstableKMathAPI
-public interface XYFit<T : Any> : Optimization<T> {
-
-    public val algebra: Field<T>
-
-    /**
-     * Set X-Y data for this fit optionally including x and y errors
-     */
-    public fun data(
-        dataSet: ColumnarData<T>,
-        xSymbol: Symbol,
-        ySymbol: Symbol,
-        xErrSymbol: Symbol? = null,
-        yErrSymbol: Symbol? = null,
-    )
-
-    public fun model(model: (T) -> DifferentiableExpression<T>)
-
-    /**
-     * Set the differentiable model for this fit
-     */
-    public fun <I : Any, A> model(
-        autoDiff: AutoDiffProcessor<T, I, A, Expression<T>>,
-        modelFunction: A.(I) -> I,
-    ): Unit where A : ExtendedField<I>, A : ExpressionAlgebra<T, I> = model { arg ->
-        autoDiff.process { modelFunction(const(arg)) }
-    }
-}

kmath-stat/src/commonMain/tmp/QowFit.kt

@@ -0,0 +1,372 @@
+/*
+ * Copyright 2018-2021 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.optimization.qow
+
+import space.kscience.kmath.data.ColumnarData
+import space.kscience.kmath.data.XYErrorColumnarData
+import space.kscience.kmath.expressions.*
+import space.kscience.kmath.linear.*
+import space.kscience.kmath.misc.UnstableKMathAPI
+import space.kscience.kmath.operations.DoubleField
+import space.kscience.kmath.operations.Field
+import space.kscience.kmath.optimization.OptimizationFeature
+import space.kscience.kmath.optimization.OptimizationProblemFactory
+import space.kscience.kmath.optimization.OptimizationResult
+import space.kscience.kmath.optimization.XYOptimization
+import space.kscience.kmath.structures.DoubleBuffer
+import space.kscience.kmath.structures.DoubleL2Norm
+import kotlin.math.pow
+
+
+private typealias ParamSet = Map<Symbol, Double>
+
+@OptIn(UnstableKMathAPI::class)
+public class QowFit(
+    override val symbols: List<Symbol>,
+    private val space: LinearSpace<Double, DoubleField>,
+    private val solver: LinearSolver<Double>,
+) : XYOptimization<Double>, SymbolIndexer {
+
+    private var logger: FitLogger? = null
+
+    private var startingPoint: Map<Symbol, Double> = TODO()
+    private var covariance: Matrix<Double>? = TODO()
+    private val prior: DifferentiableExpression<Double, Expression<Double>>? = TODO()
+    private var data: XYErrorColumnarData<Double, Double, Double> = TODO()
+    private var model: DifferentiableExpression<Double, Expression<Double>> = TODO()
+
+    private val features = HashSet<OptimizationFeature>()
+
+    override fun update(result: OptimizationResult<Double>) {
+        TODO("Not yet implemented")
+    }
+
+    override val algebra: Field<Double>
+        get() = TODO("Not yet implemented")
+
+    override fun data(
+        dataSet: ColumnarData<Double>,
+        xSymbol: Symbol,
+        ySymbol: Symbol,
+        xErrSymbol: Symbol?,
+        yErrSymbol: Symbol?,
+    ) {
+        TODO("Not yet implemented")
+    }
+
+    override fun model(model: (Double) -> DifferentiableExpression<Double, *>) {
+        TODO("Not yet implemented")
+    }
+
+    private var x: Symbol = Symbol.x
+
+    /**
+     * The signed distance from the model to the [i]-th point of data.
+     */
+    private fun distance(i: Int, parameters: Map<Symbol, Double>): Double =
+        model(parameters + (x to data.x[i])) - data.y[i]
+
+
+    /**
+     * The derivative of [distance]
+     * TODO use expressions instead
+     */
+    private fun distanceDerivative(symbol: Symbol, i: Int, parameters: Map<Symbol, Double>): Double =
+        model.derivative(symbol)(parameters + (x to data.x[i]))
+
+    /**
+     * The dispersion of [i]-th data point
+     */
+    private fun getDispersion(i: Int, parameters: Map<Symbol, Double>): Double = data.yErr[i].pow(2)
+
+    private fun getCovariance(weight: QoWeight): Matrix<Double> = solver.inverse(getEqDerivValues(weight))
+
+    /**
+     * Теоретическая ковариация весовых функций.
+     *
+     * D(\phi)=E(\phi_k(\theta_0) \phi_l(\theta_0))= disDeriv_k * disDeriv_l /sigma^2
+     */
+    private fun covarF(weight: QoWeight): Matrix<Double> = space.buildSymmetricMatrix(symbols.size) { k, l ->
+        (0 until data.size).sumOf { i -> weight.derivs[k, i] * weight.derivs[l, i] / weight.dispersion[i] }
+    }
+
+    /**
+     * Экспериментальная ковариация весов. Формула (22) из
+     * http://arxiv.org/abs/physics/0604127
+     *
+     * @param source
+     * @param set
+     * @param fitPars
+     * @param weight
+     * @return
+     */
+    private fun covarFExp(weight: QoWeight, theta: Map<Symbol, Double>): Matrix<Double> = space.run {
+        /*
+         * Важно! Если не делать предварителього вычисления этих производных, то
+         * количество вызывов функции будет dim^2 вместо dim Первый индекс -
+         * номер точки, второй - номер переменной, по которой берется производная
+         */
+        val eqvalues = buildMatrix(data.size, symbols.size) { i, l ->
+            distance(i, theta) * weight.derivs[l, i] / weight.dispersion[i]
+        }
+
+        buildMatrix(symbols.size, symbols.size) { k, l ->
+            (0 until data.size).sumOf { i -> eqvalues[i, l] * eqvalues[i, k] }
+        }
+    }
+
+    /**
+     * производные уравнений для метода Ньютона
+     *
+     * @param source
+     * @param set
+     * @param fitPars
+     * @param weight
+     * @return
+     */
+    private fun getEqDerivValues(
+        weight: QoWeight, theta: Map<Symbol, Double> = weight.theta,
+    ): Matrix<Double> = space.run {
+        val fitDim = symbols.size
+        //Возвращает производную k-того Eq по l-тому параметру
+        val res = Array(fitDim) { DoubleArray(fitDim) }
+        val sderiv = buildMatrix(data.size, symbols.size) { i, l ->
+            distanceDerivative(symbols[l], i, theta)
+        }
+
+        buildMatrix(symbols.size, symbols.size) { k, l ->
+            val base = (0 until data.size).sumOf { i ->
+                require(weight.dispersion[i] > 0)
+                sderiv[i, l] * weight.derivs[k, i] / weight.dispersion[i]
+            }
+            prior?.let { prior ->
+                //Check if this one is correct
+                val pi = prior(theta)
+                val deriv1 = prior.derivative(symbols[k])(theta)
+                val deriv2 = prior.derivative(symbols[l])(theta)
+                base + deriv1 * deriv2 / pi / pi
+            } ?: base
+        }
+    }
+
+
+    /**
+     * Значения уравнений метода квазиоптимальных весов
+     *
+     * @param source
+     * @param set
+     * @param fitPars
+     * @param weight
+     * @return
+     */
+    private fun getEqValues(weight: QoWeight, theta: Map<Symbol, Double> = weight.theta): Point<Double> {
+        val distances = DoubleBuffer(data.size) { i -> distance(i, theta) }
+
+        return DoubleBuffer(symbols.size) { k ->
+            val base = (0 until data.size).sumOf { i -> distances[i] * weight.derivs[k, i] / weight.dispersion[i] }
+            //Поправка на априорную вероятность
+            prior?.let { prior ->
+                base - prior.derivative(symbols[k])(theta) / prior(theta)
+            } ?: base
+        }
+    }
+
+
+    /**
+     * The state of QOW fitter
+     * Created by Alexander Nozik on 17-Oct-16.
+     */
+    private inner class QoWeight(
+        val theta: Map<Symbol, Double>,
+    ) {
+
+        init {
+            require(data.size > 0) { "The state does not contain data" }
+        }
+
+        /**
+         * Derivatives of the spectrum over parameters. First index in the point number, second one - index of parameter
+         */
+        val derivs: Matrix<Double> by lazy {
+            space.buildMatrix(data.size, symbols.size) { i, k ->
+                distanceDerivative(symbols[k], i, theta)
+            }
+        }
+
+        /**
+         * Array of dispersions in each point
+         */
+        val dispersion: Point<Double> by lazy {
+            DoubleBuffer(data.size) { i -> getDispersion(i, theta) }
+        }
+
+    }
+
+    private fun newtonianStep(
+        weight: QoWeight,
+        par: Map<Symbol, Double>,
+        eqvalues: Point<Double>,
+    ): Map<Symbol, Double> = space.run {
+        val start = par.toPoint()
+        val invJacob = solver.inverse(getEqDerivValues(weight, par))
+
+        val step = invJacob.dot(eqvalues)
+        return par + (start - step).toMap()
+    }
+
+    private fun newtonianRun(
+        weight: QoWeight,
+        maxSteps: Int = 100,
+        tolerance: Double = 0.0,
+        fast: Boolean = false,
+    ): ParamSet {
+
+        var dis: Double//норма невязки
+        // Для удобства работаем всегда с полным набором параметров
+        var par = startingPoint
+
+        logger?.log { "Starting newtonian iteration from: \n\t$par" }
+
+        var eqvalues = getEqValues(weight, par)//значения функций
+
+        dis = DoubleL2Norm.norm(eqvalues)// невязка
+        logger?.log { "Starting discrepancy is $dis" }
+        var i = 0
+        var flag = false
+        while (!flag) {
+            i++
+            logger?.log { "Starting step number $i" }
+
+            val currentSolution = if (fast) {
+                //Берет значения матрицы в той точке, где считается вес
+                newtonianStep(weight, weight.theta, eqvalues)
+            } else {
+                //Берет значения матрицы в точке par
+                newtonianStep(weight, par, eqvalues)
+            }
+            // здесь должен стоять учет границ параметров
+            logger?.log { "Parameter values after step are: \n\t$currentSolution" }
+
+            eqvalues = getEqValues(weight, currentSolution)
+            val currentDis = DoubleL2Norm.norm(eqvalues)// невязка после шага
+
+            logger?.log { "The discrepancy after step is: $currentDis." }
+
+            if (currentDis >= dis && i > 1) {
+                //дополнительно проверяем, чтобы был сделан хотя бы один шаг
+                flag = true
+                logger?.log { "The discrepancy does not decrease. Stopping iteration." }
+            } else {
+                par = currentSolution
+                dis = currentDis
+            }
+            if (i >= maxSteps) {
+                flag = true
+                logger?.log { "Maximum number of iterations reached. Stopping iteration." }
+            }
+            if (dis <= tolerance) {
+                flag = true
+                logger?.log { "Tolerance threshold is reached. Stopping iteration." }
+            }
+        }
+
+        return par
+    }
+
+
+//
+//    override fun run(state: FitState, parentLog: History?, meta: Meta): FitResult {
+//        val log = Chronicle("QOW", parentLog)
+//        val action = meta.getString(FIT_STAGE_TYPE, TASK_RUN)
+//        log.report("QOW fit engine started task '{}'", action)
+//        return when (action) {
+//            TASK_SINGLE -> makeRun(state, log, meta)
+//            TASK_COVARIANCE -> generateErrors(state, log, meta)
+//            TASK_RUN -> {
+//                var res = makeRun(state, log, meta)
+//                res = makeRun(res.optState().get(), log, meta)
+//                generateErrors(res.optState().get(), log, meta)
+//            }
+//            else -> throw IllegalArgumentException("Unknown task")
+//        }
+//    }
+
+//    private fun makeRun(state: FitState, log: History, meta: Meta): FitResult {
+//        /*Инициализация объектов, задание исходных значений*/
+//        log.report("Starting fit using quasioptimal weights method.")
+//
+//        val fitPars = getFitPars(state, meta)
+//
+//        val curWeight = QoWeight(state, fitPars, state.parameters)
+//
+//        // вычисляем вес в allPar. Потом можно будет попробовать ручное задание веса
+//        log.report("The starting weight is: \n\t{}",
+//            MathUtils.toString(curWeight.theta))
+//
+//        //Стартовая точка такая же как и параметр веса
+//        /*Фитирование*/
+//        val res = newtonianRun(state, curWeight, log, meta)
+//
+//        /*Генерация результата*/
+//
+//        return FitResult.build(state.edit().setPars(res).build(), *fitPars)
+//    }
+
+    /**
+     * generateErrors.
+     */
+    private fun generateErrors(): Matrix<Double> {
+        logger?.log { """
+            Starting errors estimation using quasioptimal weights method. The starting weight is:
+                ${curWeight.theta}
+             """.trimIndent()}
+        val curWeight = QoWeight(startingPoint)
+
+        val covar = getCovariance(curWeight)
+
+        val decomposition = EigenDecomposition(covar.matrix)
+        var valid = true
+        for (lambda in decomposition.realEigenvalues) {
+            if (lambda <= 0) {
+                log.report("The covariance matrix is not positive defined. Error estimation is not valid")
+                valid = false
+            }
+        }
+    }
+
+
+    override suspend fun optimize(): OptimizationResult<Double> {
+        val curWeight = QoWeight(startingPoint)
+        logger?.log {
+            """
+            Starting fit using quasioptimal weights method. The starting weight is: 
+	            ${curWeight.theta}
+            """.trimIndent()
+        }
+        val res = newtonianRun(curWeight)
+    }
+
+
+    companion object : OptimizationProblemFactory<Double, QowFit> {
+        override fun build(symbols: List<Symbol>): QowFit {
+            TODO("Not yet implemented")
+        }
+
+
+        /**
+         * Constant `QOW_ENGINE_NAME="QOW"`
+         */
+        const val QOW_ENGINE_NAME = "QOW"
+
+        /**
+         * Constant `QOW_METHOD_FAST="fast"`
+         */
+        const val QOW_METHOD_FAST = "fast"
+
+
+    }
+}
+

kmath-stat/src/commonMain/tmp/minuit/AnalyticalGradientCalculator.kt

@@ -0,0 +1,61 @@
+/* 
+ * Copyright 2015 Alexander Nozik.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package ru.inr.mass.minuit
+
+import ru.inr.mass.maths.MultiFunction
+
+/**
+ *
+ * @version $Id$
+ */
+internal class AnalyticalGradientCalculator(fcn: MultiFunction?, state: MnUserTransformation, checkGradient: Boolean) :
+    GradientCalculator {
+    private val function: MultiFunction?
+    private val theCheckGradient: Boolean
+    private val theTransformation: MnUserTransformation
+    fun checkGradient(): Boolean {
+        return theCheckGradient
+    }
+
+    /** {@inheritDoc}  */
+    fun gradient(par: MinimumParameters): FunctionGradient {
+//      double[] grad = theGradCalc.gradientValue(theTransformation.andThen(par.vec()).data());
+        val point: DoubleArray = theTransformation.transform(par.vec()).toArray()
+        require(!(function.getDimension() !== theTransformation.parameters().size())) { "Invalid parameter size" }
+        val v: RealVector = ArrayRealVector(par.vec().getDimension())
+        for (i in 0 until par.vec().getDimension()) {
+            val ext: Int = theTransformation.extOfInt(i)
+            if (theTransformation.parameter(ext).hasLimits()) {
+                val dd: Double = theTransformation.dInt2Ext(i, par.vec().getEntry(i))
+                v.setEntry(i, dd * function.derivValue(ext, point))
+            } else {
+                v.setEntry(i, function.derivValue(ext, point))
+            }
+        }
+        return FunctionGradient(v)
+    }
+
+    /** {@inheritDoc}  */
+    fun gradient(par: MinimumParameters, grad: FunctionGradient?): FunctionGradient {
+        return gradient(par)
+    }
+
+    init {
+        function = fcn
+        theTransformation = state
+        theCheckGradient = checkGradient
+    }
+}

kmath-stat/src/commonMain/tmp/minuit/CombinedMinimizer.kt

@@ -0,0 +1,32 @@
+/* 
+ * Copyright 2015 Alexander Nozik.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package ru.inr.mass.minuit
+
+/**
+ *
+ * @version $Id$
+ */
+internal class CombinedMinimizer : ModularFunctionMinimizer() {
+    private val theMinBuilder: CombinedMinimumBuilder = CombinedMinimumBuilder()
+    private val theMinSeedGen: MnSeedGenerator = MnSeedGenerator()
+    override fun builder(): MinimumBuilder {
+        return theMinBuilder
+    }
+
+    override fun seedGenerator(): MinimumSeedGenerator {
+        return theMinSeedGen
+    }
+}

kmath-stat/src/commonMain/tmp/minuit/CombinedMinimumBuilder.kt

@@ -0,0 +1,58 @@
+/* 
+ * Copyright 2015 Alexander Nozik.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package ru.inr.mass.minuit
+
+import space.kscience.kmath.optimization.minuit.MINUITPlugin
+import space.kscience.kmath.optimization.minuit.MinimumSeed
+
+/**
+ *
+ * @version $Id$
+ */
+internal class CombinedMinimumBuilder : MinimumBuilder {
+    private val theSimplexMinimizer: SimplexMinimizer = SimplexMinimizer()
+    private val theVMMinimizer: VariableMetricMinimizer = VariableMetricMinimizer()
+
+    /** {@inheritDoc}  */
+    override fun minimum(
+        fcn: MnFcn?,
+        gc: GradientCalculator?,
+        seed: MinimumSeed?,
+        strategy: MnStrategy?,
+        maxfcn: Int,
+        toler: Double
+    ): FunctionMinimum {
+        val min: FunctionMinimum = theVMMinimizer.minimize(fcn!!, gc, seed, strategy, maxfcn, toler)
+        if (!min.isValid()) {
+            MINUITPlugin.logStatic("CombinedMinimumBuilder: migrad method fails, will try with simplex method first.")
+            val str = MnStrategy(2)
+            val min1: FunctionMinimum = theSimplexMinimizer.minimize(fcn, gc, seed, str, maxfcn, toler)
+            if (!min1.isValid()) {
+                MINUITPlugin.logStatic("CombinedMinimumBuilder: both migrad and simplex method fail.")
+                return min1
+            }
+            val seed1: MinimumSeed = theVMMinimizer.seedGenerator().generate(fcn, gc, min1.userState(), str)
+            val min2: FunctionMinimum = theVMMinimizer.minimize(fcn, gc, seed1, str, maxfcn, toler)
+            if (!min2.isValid()) {
+                MINUITPlugin.logStatic("CombinedMinimumBuilder: both migrad and method fails also at 2nd attempt.")
+                MINUITPlugin.logStatic("CombinedMinimumBuilder: return simplex minimum.")
+                return min1
+            }
+            return min2
+        }
+        return min
+    }
+}

kmath-stat/src/commonMain/tmp/minuit/ContoursError.kt

@@ -0,0 +1,150 @@
+/* 
+ * Copyright 2015 Alexander Nozik.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package ru.inr.mass.minuit
+
+/**
+ *
+ * ContoursError class.
+ *
+ * @author Darksnake
+ * @version $Id$
+ */
+class ContoursError internal constructor(
+    private val theParX: Int,
+    private val theParY: Int,
+    points: List<Range>,
+    xmnos: MinosError,
+    ymnos: MinosError,
+    nfcn: Int
+) {
+    private val theNFcn: Int
+    private val thePoints: List<Range> = points
+    private val theXMinos: MinosError
+    private val theYMinos: MinosError
+
+    /**
+     *
+     * nfcn.
+     *
+     * @return a int.
+     */
+    fun nfcn(): Int {
+        return theNFcn
+    }
+
+    /**
+     *
+     * points.
+     *
+     * @return a [List] object.
+     */
+    fun points(): List<Range> {
+        return thePoints
+    }
+
+    /**
+     * {@inheritDoc}
+     */
+    override fun toString(): String {
+        return MnPrint.toString(this)
+    }
+
+    /**
+     *
+     * xMinosError.
+     *
+     * @return a [hep.dataforge.MINUIT.MinosError] object.
+     */
+    fun xMinosError(): MinosError {
+        return theXMinos
+    }
+
+    /**
+     *
+     * xRange.
+     *
+     * @return
+     */
+    fun xRange(): Range {
+        return theXMinos.range()
+    }
+
+    /**
+     *
+     * xmin.
+     *
+     * @return a double.
+     */
+    fun xmin(): Double {
+        return theXMinos.min()
+    }
+
+    /**
+     *
+     * xpar.
+     *
+     * @return a int.
+     */
+    fun xpar(): Int {
+        return theParX
+    }
+
+    /**
+     *
+     * yMinosError.
+     *
+     * @return a [hep.dataforge.MINUIT.MinosError] object.
+     */
+    fun yMinosError(): MinosError {
+        return theYMinos
+    }
+
+    /**
+     *
+     * yRange.
+     *
+     * @return
+     */
+    fun yRange(): Range {
+        return theYMinos.range()
+    }
+
+    /**
+     *
+     * ymin.
+     *
+     * @return a double.
+     */
+    fun ymin(): Double {
+        return theYMinos.min()
+    }
+
+    /**
+     *
+     * ypar.
+     *
+     * @return a int.
+     */
+    fun ypar(): Int {
+        return theParY
+    }
+
+    init {
+        theXMinos = xmnos
+        theYMinos = ymnos
+        theNFcn = nfcn
+    }
+}

kmath-stat/src/commonMain/tmp/minuit/DavidonErrorUpdator.kt

@@ -0,0 +1,45 @@
+/* 
+ * Copyright 2015 Alexander Nozik.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package ru.inr.mass.minuit
+
+import org.apache.commons.math3.linear.RealVector
+import ru.inr.mass.minuit.*
+
+/**
+ *
+ * @version $Id$
+ */
+internal class DavidonErrorUpdator : MinimumErrorUpdator {
+    /** {@inheritDoc}  */
+    fun update(s0: MinimumState, p1: MinimumParameters, g1: FunctionGradient): MinimumError {
+        val V0: MnAlgebraicSymMatrix = s0.error().invHessian()
+        val dx: RealVector = MnUtils.sub(p1.vec(), s0.vec())
+        val dg: RealVector = MnUtils.sub(g1.getGradient(), s0.gradient().getGradient())
+        val delgam: Double = MnUtils.innerProduct(dx, dg)
+        val gvg: Double = MnUtils.similarity(dg, V0)
+        val vg: RealVector = MnUtils.mul(V0, dg)
+        var Vupd: MnAlgebraicSymMatrix =
+            MnUtils.sub(MnUtils.div(MnUtils.outerProduct(dx), delgam), MnUtils.div(MnUtils.outerProduct(vg), gvg))
+        if (delgam > gvg) {
+            Vupd = MnUtils.add(Vupd,
+                MnUtils.mul(MnUtils.outerProduct(MnUtils.sub(MnUtils.div(dx, delgam), MnUtils.div(vg, gvg))), gvg))
+        }
+        val sum_upd: Double = MnUtils.absoluteSumOfElements(Vupd)
+        Vupd = MnUtils.add(Vupd, V0)
+        val dcov: Double = 0.5 * (s0.error().dcovar() + sum_upd / MnUtils.absoluteSumOfElements(Vupd))
+        return MinimumError(Vupd, dcov)
+    }
+}

kmath-stat/src/commonMain/tmp/minuit/FunctionGradient.kt

@@ -0,0 +1,72 @@
+/* 
+ * Copyright 2015 Alexander Nozik.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package ru.inr.mass.minuit
+
+import org.apache.commons.math3.linear.ArrayRealVector
+
+/**
+ *
+ * @version $Id$
+ */
+class FunctionGradient {
+    private var theAnalytical = false
+    private var theG2ndDerivative: RealVector
+    private var theGStepSize: RealVector
+    private var theGradient: RealVector
+    private var theValid = false
+
+    constructor(n: Int) {
+        theGradient = ArrayRealVector(n)
+        theG2ndDerivative = ArrayRealVector(n)
+        theGStepSize = ArrayRealVector(n)
+    }
+
+    constructor(grd: RealVector) {
+        theGradient = grd
+        theG2ndDerivative = ArrayRealVector(grd.getDimension())
+        theGStepSize = ArrayRealVector(grd.getDimension())
+        theValid = true
+        theAnalytical = true
+    }
+
+    constructor(grd: RealVector, g2: RealVector, gstep: RealVector) {
+        theGradient = grd
+        theG2ndDerivative = g2
+        theGStepSize = gstep
+        theValid = true
+        theAnalytical = false
+    }
+
+    fun getGradient(): RealVector {
+        return theGradient
+    }
+
+    fun getGradientDerivative(): RealVector {
+        return theG2ndDerivative
+    }
+
+    fun getStep(): RealVector {
+        return theGStepSize
+    }
+
+    fun isAnalytical(): Boolean {
+        return theAnalytical
+    }
+
+    fun isValid(): Boolean {
+        return theValid
+    }
+}

kmath-stat/src/commonMain/tmp/minuit/FunctionMinimum.kt

@@ -0,0 +1,260 @@
+/* 
+ * Copyright 2015 Alexander Nozik.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package ru.inr.mass.minuit
+
+import ru.inr.mass.minuit.*
+import space.kscience.kmath.optimization.minuit.MinimumSeed
+
+/**
+ * Result of the minimization.
+ *
+ *
+ * The FunctionMinimum is the output of the minimizers and contains the
+ * minimization result. The methods
+ *
+ *  * userState(),
+ *  * userParameters() and
+ *  * userCovariance()
+ *
+ * are provided. These can be used as new input to a new minimization after some
+ * manipulation. The parameters and/or the FunctionMinimum can be printed using
+ * the toString() method or the MnPrint class.
+ *
+ * @author Darksnake
+ */
+class FunctionMinimum {
+    private var theAboveMaxEdm = false
+    private var theErrorDef: Double
+    private var theReachedCallLimit = false
+    private var theSeed: MinimumSeed
+    private var theStates: MutableList<MinimumState>
+    private var theUserState: MnUserParameterState
+
+    internal constructor(seed: MinimumSeed, up: Double) {
+        theSeed = seed
+        theStates = ArrayList()
+        theStates.add(MinimumState(seed.parameters(),
+            seed.error(),
+            seed.gradient(),
+            seed.parameters().fval(),
+            seed.nfcn()))
+        theErrorDef = up
+        theUserState = MnUserParameterState()
+    }
+
+    internal constructor(seed: MinimumSeed, states: MutableList<MinimumState>, up: Double) {
+        theSeed = seed
+        theStates = states
+        theErrorDef = up
+        theUserState = MnUserParameterState()
+    }
+
+    internal constructor(seed: MinimumSeed, states: MutableList<MinimumState>, up: Double, x: MnReachedCallLimit?) {
+        theSeed = seed
+        theStates = states
+        theErrorDef = up
+        theReachedCallLimit = true
+        theUserState = MnUserParameterState()
+    }
+
+    internal constructor(seed: MinimumSeed, states: MutableList<MinimumState>, up: Double, x: MnAboveMaxEdm?) {
+        theSeed = seed
+        theStates = states
+        theErrorDef = up
+        theAboveMaxEdm = true
+        theReachedCallLimit = false
+        theUserState = MnUserParameterState()
+    }
+
+    // why not
+    fun add(state: MinimumState) {
+        theStates.add(state)
+    }
+
+    /**
+     * returns the expected vertical distance to the minimum (EDM)
+     *
+     * @return a double.
+     */
+    fun edm(): Double {
+        return lastState().edm()
+    }
+
+    fun error(): MinimumError {
+        return lastState().error()
+    }
+
+    /**
+     *
+     *
+     * errorDef.
+     *
+     * @return a double.
+     */
+    fun errorDef(): Double {
+        return theErrorDef
+    }
+
+    /**
+     * Returns the function value at the minimum.
+     *
+     * @return a double.
+     */
+    fun fval(): Double {
+        return lastState().fval()
+    }
+
+    fun grad(): FunctionGradient {
+        return lastState().gradient()
+    }
+
+    fun hasAccurateCovar(): Boolean {
+        return state().error().isAccurate()
+    }
+
+    fun hasCovariance(): Boolean {
+        return state().error().isAvailable()
+    }
+
+    fun hasMadePosDefCovar(): Boolean {
+        return state().error().isMadePosDef()
+    }
+
+    fun hasPosDefCovar(): Boolean {
+        return state().error().isPosDef()
+    }
+
+    fun hasReachedCallLimit(): Boolean {
+        return theReachedCallLimit
+    }
+
+    fun hasValidCovariance(): Boolean {
+        return state().error().isValid()
+    }
+
+    fun hasValidParameters(): Boolean {
+        return state().parameters().isValid()
+    }
+
+    fun hesseFailed(): Boolean {
+        return state().error().hesseFailed()
+    }
+
+    fun isAboveMaxEdm(): Boolean {
+        return theAboveMaxEdm
+    }
+
+    /**
+     * In general, if this returns <CODE>true</CODE>, the minimizer did find a
+     * minimum without running into troubles. However, in some cases a minimum
+     * cannot be found, then the return value will be <CODE>false</CODE>.
+     * Reasons for the minimization to fail are
+     *
+     *  * the number of allowed function calls has been exhausted
+     *  * the minimizer could not improve the values of the parameters (and
+     * knowing that it has not converged yet)
+     *  * a problem with the calculation of the covariance matrix
+     *
+     * Additional methods for the analysis of the state at the minimum are
+     * provided.
+     *
+     * @return a boolean.
+     */
+    fun isValid(): Boolean {
+        return state().isValid() && !isAboveMaxEdm() && !hasReachedCallLimit()
+    }
+
+    private fun lastState(): MinimumState {
+        return theStates[theStates.size - 1]
+    }
+    // forward interface of last state
+    /**
+     * returns the total number of function calls during the minimization.
+     *
+     * @return a int.
+     */
+    fun nfcn(): Int {
+        return lastState().nfcn()
+    }
+
+    fun parameters(): MinimumParameters {
+        return lastState().parameters()
+    }
+
+    fun seed(): MinimumSeed {
+        return theSeed
+    }
+
+    fun state(): MinimumState {
+        return lastState()
+    }
+
+    fun states(): List<MinimumState> {
+        return theStates
+    }
+
+    /**
+     * {@inheritDoc}
+     *
+     * @return
+     */
+    override fun toString(): String {
+        return MnPrint.toString(this)
+    }
+
+    /**
+     *
+     *
+     * userCovariance.
+     *
+     * @return a [hep.dataforge.MINUIT.MnUserCovariance] object.
+     */
+    fun userCovariance(): MnUserCovariance {
+        if (!theUserState.isValid()) {
+            theUserState = MnUserParameterState(state(), errorDef(), seed().trafo())
+        }
+        return theUserState.covariance()
+    }
+
+    /**
+     *
+     *
+     * userParameters.
+     *
+     * @return a [hep.dataforge.MINUIT.MnUserParameters] object.
+     */
+    fun userParameters(): MnUserParameters {
+        if (!theUserState.isValid()) {
+            theUserState = MnUserParameterState(state(), errorDef(), seed().trafo())
+        }
+        return theUserState.parameters()
+    }
+
+    /**
+     * user representation of state at minimum
+     *
+     * @return a [hep.dataforge.MINUIT.MnUserParameterState] object.
+     */
+    fun userState(): MnUserParameterState {
+        if (!theUserState.isValid()) {
+            theUserState = MnUserParameterState(state(), errorDef(), seed().trafo())
+        }
+        return theUserState
+    }
+
+    internal class MnAboveMaxEdm
+    internal class MnReachedCallLimit
+}

kmath-stat/src/commonMain/tmp/minuit/GradientCalculator.kt

@@ -0,0 +1,41 @@
+/* 
+ * Copyright 2015 Alexander Nozik.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package ru.inr.mass.minuit
+
+/**
+ *
+ * @version $Id$
+ */
+interface GradientCalculator {
+    /**
+     *
+     * gradient.
+     *
+     * @param par a [hep.dataforge.MINUIT.MinimumParameters] object.
+     * @return a [hep.dataforge.MINUIT.FunctionGradient] object.
+     */
+    fun gradient(par: MinimumParameters?): FunctionGradient
+
+    /**
+     *
+     * gradient.
+     *
+     * @param par a [hep.dataforge.MINUIT.MinimumParameters] object.
+     * @param grad a [hep.dataforge.MINUIT.FunctionGradient] object.
+     * @return a [hep.dataforge.MINUIT.FunctionGradient] object.
+     */
+    fun gradient(par: MinimumParameters?, grad: FunctionGradient?): FunctionGradient
+}

kmath-stat/src/commonMain/tmp/minuit/HessianGradientCalculator.kt

@@ -0,0 +1,137 @@
+/* 
+ * Copyright 2015 Alexander Nozik.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package ru.inr.mass.minuit
+
+import org.apache.commons.math3.linear.ArrayRealVector
+import ru.inr.mass.minuit.*
+
+/**
+ *
+ * @version $Id$
+ */
+internal class HessianGradientCalculator(fcn: MnFcn, par: MnUserTransformation, stra: MnStrategy) : GradientCalculator {
+    private val theFcn: MnFcn = fcn
+    private val theStrategy: MnStrategy
+    private val theTransformation: MnUserTransformation
+    fun deltaGradient(par: MinimumParameters, gradient: FunctionGradient): Pair<FunctionGradient, RealVector> {
+        require(par.isValid()) { "parameters are invalid" }
+        val x: RealVector = par.vec().copy()
+        val grd: RealVector = gradient.getGradient().copy()
+        val g2: RealVector = gradient.getGradientDerivative()
+        val gstep: RealVector = gradient.getStep()
+        val fcnmin: Double = par.fval()
+        //   std::cout<<"fval: "<<fcnmin<<std::endl;
+        val dfmin: Double = 4.0 * precision().eps2() * (abs(fcnmin) + theFcn.errorDef())
+        val n: Int = x.getDimension()
+        val dgrd: RealVector = ArrayRealVector(n)
+
+        // initial starting values
+        for (i in 0 until n) {
+            val xtf: Double = x.getEntry(i)
+            val dmin: Double = 4.0 * precision().eps2() * (xtf + precision().eps2())
+            val epspri: Double = precision().eps2() + abs(grd.getEntry(i) * precision().eps2())
+            val optstp: Double = sqrt(dfmin / (abs(g2.getEntry(i)) + epspri))
+            var d: Double = 0.2 * abs(gstep.getEntry(i))
+            if (d > optstp) {
+                d = optstp
+            }
+            if (d < dmin) {
+                d = dmin
+            }
+            var chgold = 10000.0
+            var dgmin = 0.0
+            var grdold = 0.0
+            var grdnew = 0.0
+            for (j in 0 until ncycle()) {
+                x.setEntry(i, xtf + d)
+                val fs1: Double = theFcn.value(x)
+                x.setEntry(i, xtf - d)
+                val fs2: Double = theFcn.value(x)
+                x.setEntry(i, xtf)
+                //       double sag = 0.5*(fs1+fs2-2.*fcnmin);
+                grdold = grd.getEntry(i)
+                grdnew = (fs1 - fs2) / (2.0 * d)
+                dgmin = precision().eps() * (abs(fs1) + abs(fs2)) / d
+                if (abs(grdnew) < precision().eps()) {
+                    break
+                }
+                val change: Double = abs((grdold - grdnew) / grdnew)
+                if (change > chgold && j > 1) {
+                    break
+                }
+                chgold = change
+                grd.setEntry(i, grdnew)
+                if (change < 0.05) {
+                    break
+                }
+                if (abs(grdold - grdnew) < dgmin) {
+                    break
+                }
+                if (d < dmin) {
+                    break
+                }
+                d *= 0.2
+            }
+            dgrd.setEntry(i, max(dgmin, abs(grdold - grdnew)))
+        }
+        return Pair(FunctionGradient(grd, g2, gstep), dgrd)
+    }
+
+    fun fcn(): MnFcn {
+        return theFcn
+    }
+
+    fun gradTolerance(): Double {
+        return strategy().gradientTolerance()
+    }
+
+    /** {@inheritDoc}  */
+    fun gradient(par: MinimumParameters): FunctionGradient {
+        val gc = InitialGradientCalculator(theFcn, theTransformation, theStrategy)
+        val gra: FunctionGradient = gc.gradient(par)
+        return gradient(par, gra)
+    }
+
+    /** {@inheritDoc}  */
+    fun gradient(par: MinimumParameters, gradient: FunctionGradient): FunctionGradient {
+        return deltaGradient(par, gradient).getFirst()
+    }
+
+    fun ncycle(): Int {
+        return strategy().hessianGradientNCycles()
+    }
+
+    fun precision(): MnMachinePrecision {
+        return theTransformation.precision()
+    }
+
+    fun stepTolerance(): Double {
+        return strategy().gradientStepTolerance()
+    }
+
+    fun strategy(): MnStrategy {
+        return theStrategy
+    }
+
+    fun trafo(): MnUserTransformation {
+        return theTransformation
+    }
+
+    init {
+        theTransformation = par
+        theStrategy = stra
+    }
+}

kmath-stat/src/commonMain/tmp/minuit/InitialGradientCalculator.kt

@@ -0,0 +1,116 @@
+/* 
+ * Copyright 2015 Alexander Nozik.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package ru.inr.mass.minuit
+
+import org.apache.commons.math3.linear.ArrayRealVector
+import ru.inr.mass.minuit.*
+
+/**
+ * Calculating derivatives via finite differences
+ * @version $Id$
+ */
+internal class InitialGradientCalculator(fcn: MnFcn, par: MnUserTransformation, stra: MnStrategy) {
+    private val theFcn: MnFcn = fcn
+    private val theStrategy: MnStrategy
+    private val theTransformation: MnUserTransformation
+    fun fcn(): MnFcn {
+        return theFcn
+    }
+
+    fun gradTolerance(): Double {
+        return strategy().gradientTolerance()
+    }
+
+    fun gradient(par: MinimumParameters): FunctionGradient {
+        require(par.isValid()) { "Parameters are invalid" }
+        val n: Int = trafo().variableParameters()
+        require(n == par.vec().getDimension()) { "Parameters have invalid size" }
+        val gr: RealVector = ArrayRealVector(n)
+        val gr2: RealVector = ArrayRealVector(n)
+        val gst: RealVector = ArrayRealVector(n)
+
+        // initial starting values
+        for (i in 0 until n) {
+            val exOfIn: Int = trafo().extOfInt(i)
+            val `var`: Double = par.vec().getEntry(i) //parameter value
+            val werr: Double = trafo().parameter(exOfIn).error() //parameter error
+            val sav: Double = trafo().int2ext(i, `var`) //value after transformation
+            var sav2 = sav + werr //value after transfomation + error
+            if (trafo().parameter(exOfIn).hasLimits()) {
+                if (trafo().parameter(exOfIn).hasUpperLimit()
+                    && sav2 > trafo().parameter(exOfIn).upperLimit()
+                ) {
+                    sav2 = trafo().parameter(exOfIn).upperLimit()
+                }
+            }
+            var var2: Double = trafo().ext2int(exOfIn, sav2)
+            val vplu = var2 - `var`
+            sav2 = sav - werr
+            if (trafo().parameter(exOfIn).hasLimits()) {
+                if (trafo().parameter(exOfIn).hasLowerLimit()
+                    && sav2 < trafo().parameter(exOfIn).lowerLimit()
+                ) {
+                    sav2 = trafo().parameter(exOfIn).lowerLimit()
+                }
+            }
+            var2 = trafo().ext2int(exOfIn, sav2)
+            val vmin = var2 - `var`
+            val dirin: Double = 0.5 * (abs(vplu) + abs(vmin))
+            val g2: Double = 2.0 * theFcn.errorDef() / (dirin * dirin)
+            val gsmin: Double = 8.0 * precision().eps2() * (abs(`var`) + precision().eps2())
+            var gstep: Double = max(gsmin, 0.1 * dirin)
+            val grd = g2 * dirin
+            if (trafo().parameter(exOfIn).hasLimits()) {
+                if (gstep > 0.5) {
+                    gstep = 0.5
+                }
+            }
+            gr.setEntry(i, grd)
+            gr2.setEntry(i, g2)
+            gst.setEntry(i, gstep)
+        }
+        return FunctionGradient(gr, gr2, gst)
+    }
+
+    fun gradient(par: MinimumParameters, gra: FunctionGradient?): FunctionGradient {
+        return gradient(par)
+    }
+
+    fun ncycle(): Int {
+        return strategy().gradientNCycles()
+    }
+
+    fun precision(): MnMachinePrecision {
+        return theTransformation.precision()
+    }
+
+    fun stepTolerance(): Double {
+        return strategy().gradientStepTolerance()
+    }
+
+    fun strategy(): MnStrategy {
+        return theStrategy
+    }
+
+    fun trafo(): MnUserTransformation {
+        return theTransformation
+    }
+
+    init {
+        theTransformation = par
+        theStrategy = stra
+    }
+}

kmath-stat/src/commonMain/tmp/minuit/MINOSResult.kt

@@ -0,0 +1,70 @@
+/* 
+ * Copyright 2015 Alexander Nozik.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package space.kscience.kmath.optimization.minuit
+
+
+/**
+ * Контейнер для несимметричных оценок и доверительных интервалов
+ *
+ * @author Darksnake
+ * @version $Id: $Id
+ */
+class MINOSResult
+/**
+ *
+ * Constructor for MINOSResult.
+ *
+ * @param list an array of [String] objects.
+ */(private val names: Array<String>, private val errl: DoubleArray?, private val errp: DoubleArray?) :
+    IntervalEstimate {
+    fun getNames(): NameList {
+        return NameList(names)
+    }
+
+    fun getInterval(parName: String?): Pair<Value, Value> {
+        val index: Int = getNames().getNumberByName(parName)
+        return Pair(ValueFactory.of(errl!![index]), ValueFactory.of(errp!![index]))
+    }
+
+    val cL: Double
+        get() = 0.68
+
+    /** {@inheritDoc}  */
+    fun print(out: PrintWriter) {
+        if (errl != null || errp != null) {
+            out.println()
+            out.println("Assymetrical errors:")
+            out.println()
+            out.println("Name\tLower\tUpper")
+            for (i in 0 until getNames().size()) {
+                out.print(getNames().get(i))
+                out.print("\t")
+                if (errl != null) {
+                    out.print(errl[i])
+                } else {
+                    out.print("---")
+                }
+                out.print("\t")
+                if (errp != null) {
+                    out.print(errp[i])
+                } else {
+                    out.print("---")
+                }
+                out.println()
+            }
+        }
+    }
+}

kmath-stat/src/commonMain/tmp/minuit/MINUITFitter.kt

@@ -0,0 +1,205 @@
+/*
+ * Copyright 2018-2021 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.optimization.minuit
+
+import ru.inr.mass.minuit.*
+
+/**
+ *
+ *
+ * MINUITFitter class.
+ *
+ * @author Darksnake
+ * @version $Id: $Id
+ */
+class MINUITFitter : Fitter {
+    fun run(state: FitState, parentLog: History?, meta: Meta): FitResult {
+        val log = Chronicle("MINUIT", parentLog)
+        val action: String = meta.getString("action", TASK_RUN)
+        log.report("MINUIT fit engine started action '{}'", action)
+        return when (action) {
+            TASK_COVARIANCE -> runHesse(state, log, meta)
+            TASK_SINGLE, TASK_RUN -> runFit(state, log, meta)
+            else -> throw IllegalArgumentException("Unknown task")
+        }
+    }
+
+    @NotNull
+    fun getName(): String {
+        return MINUIT_ENGINE_NAME
+    }
+
+    /**
+     *
+     *
+     * runHesse.
+     *
+     * @param state a [hep.dataforge.stat.fit.FitState] object.
+     * @param log
+     * @return a [FitResult] object.
+     */
+    fun runHesse(state: FitState, log: History, meta: Meta?): FitResult {
+        val strategy: Int
+        strategy = Global.INSTANCE.getInt("MINUIT_STRATEGY", 2)
+        log.report("Generating errors using MnHesse 2-nd order gradient calculator.")
+        val fcn: MultiFunction
+        val fitPars: Array<String> = Fitter.Companion.getFitPars(state, meta)
+        val pars: ParamSet = state.getParameters()
+        fcn = MINUITUtils.getFcn(state, pars, fitPars)
+        val hesse = MnHesse(strategy)
+        val mnState: MnUserParameterState = hesse.calculate(fcn, MINUITUtils.getFitParameters(pars, fitPars))
+        val allPars: ParamSet = pars.copy()
+        for (fitPar in fitPars) {
+            allPars.setParValue(fitPar, mnState.value(fitPar))
+            allPars.setParError(fitPar, mnState.error(fitPar))
+        }
+        val newState: FitState.Builder = state.edit()
+        newState.setPars(allPars)
+        if (mnState.hasCovariance()) {
+            val mnCov: MnUserCovariance = mnState.covariance()
+            var j: Int
+            val cov = Array(mnState.variableParameters()) { DoubleArray(mnState.variableParameters()) }
+            for (i in 0 until mnState.variableParameters()) {
+                j = 0
+                while (j < mnState.variableParameters()) {
+                    cov[i][j] = mnCov.get(i, j)
+                    j++
+                }
+            }
+            newState.setCovariance(NamedMatrix(fitPars, cov), true)
+        }
+        return FitResult.build(newState.build(), fitPars)
+    }
+
+    fun runFit(state: FitState, log: History, meta: Meta): FitResult {
+        val minuit: MnApplication
+        log.report("Starting fit using Minuit.")
+        val strategy: Int
+        strategy = Global.INSTANCE.getInt("MINUIT_STRATEGY", 2)
+        var force: Boolean
+        force = Global.INSTANCE.getBoolean("FORCE_DERIVS", false)
+        val fitPars: Array<String> = Fitter.Companion.getFitPars(state, meta)
+        for (fitPar in fitPars) {
+            if (!state.modelProvidesDerivs(fitPar)) {
+                force = true
+                log.reportError("Model does not provide derivatives for parameter '{}'", fitPar)
+            }
+        }
+        if (force) {
+            log.report("Using MINUIT gradient calculator.")
+        }
+        val fcn: MultiFunction
+        val pars: ParamSet = state.getParameters().copy()
+        fcn = MINUITUtils.getFcn(state, pars, fitPars)
+        val method: String = meta.getString("method", MINUIT_MIGRAD)
+        when (method) {
+            MINUIT_MINOS, MINUIT_MINIMIZE -> minuit =
+                MnMinimize(fcn, MINUITUtils.getFitParameters(pars, fitPars), strategy)
+            MINUIT_SIMPLEX -> minuit = MnSimplex(fcn, MINUITUtils.getFitParameters(pars, fitPars), strategy)
+            else -> minuit = MnMigrad(fcn, MINUITUtils.getFitParameters(pars, fitPars), strategy)
+        }
+        if (force) {
+            minuit.setUseAnalyticalDerivatives(false)
+            log.report("Forced to use MINUIT internal derivative calculator!")
+        }
+
+//        minuit.setUseAnalyticalDerivatives(true);
+        val minimum: FunctionMinimum
+        val maxSteps: Int = meta.getInt("iterations", -1)
+        val tolerance: Double = meta.getDouble("tolerance", -1)
+        minimum = if (maxSteps > 0) {
+            if (tolerance > 0) {
+                minuit.minimize(maxSteps, tolerance)
+            } else {
+                minuit.minimize(maxSteps)
+            }
+        } else {
+            minuit.minimize()
+        }
+        if (!minimum.isValid()) {
+            log.report("Minimization failed!")
+        }
+        log.report("MINUIT run completed in {} function calls.", minimum.nfcn())
+
+        /*
+         * Генерация результата
+         */
+        val allPars: ParamSet = pars.copy()
+        for (fitPar in fitPars) {
+            allPars.setParValue(fitPar, minimum.userParameters().value(fitPar))
+            allPars.setParError(fitPar, minimum.userParameters().error(fitPar))
+        }
+        val newState: FitState.Builder = state.edit()
+        newState.setPars(allPars)
+        var valid: Boolean = minimum.isValid()
+        if (minimum.userCovariance().nrow() > 0) {
+            var j: Int
+            val cov = Array(minuit.variableParameters()) { DoubleArray(minuit.variableParameters()) }
+            if (cov[0].length == 1) {
+                cov[0][0] = minimum.userParameters().error(0) * minimum.userParameters().error(0)
+            } else {
+                for (i in 0 until minuit.variableParameters()) {
+                    j = 0
+                    while (j < minuit.variableParameters()) {
+                        cov[i][j] = minimum.userCovariance().get(i, j)
+                        j++
+                    }
+                }
+            }
+            newState.setCovariance(NamedMatrix(fitPars, cov), true)
+        }
+        if (method == MINUIT_MINOS) {
+            log.report("Starting MINOS procedure for precise error estimation.")
+            val minos = MnMinos(fcn, minimum, strategy)
+            var mnError: MinosError
+            val errl = DoubleArray(fitPars.size)
+            val errp = DoubleArray(fitPars.size)
+            for (i in fitPars.indices) {
+                mnError = minos.minos(i)
+                if (mnError.isValid()) {
+                    errl[i] = mnError.lower()
+                    errp[i] = mnError.upper()
+                } else {
+                    valid = false
+                }
+            }
+            val minosErrors = MINOSResult(fitPars, errl, errp)
+            newState.setInterval(minosErrors)
+        }
+        return FitResult.build(newState.build(), valid, fitPars)
+    }
+
+    companion object {
+        /**
+         * Constant `MINUIT_MIGRAD="MIGRAD"`
+         */
+        const val MINUIT_MIGRAD = "MIGRAD"
+
+        /**
+         * Constant `MINUIT_MINIMIZE="MINIMIZE"`
+         */
+        const val MINUIT_MINIMIZE = "MINIMIZE"
+
+        /**
+         * Constant `MINUIT_SIMPLEX="SIMPLEX"`
+         */
+        const val MINUIT_SIMPLEX = "SIMPLEX"
+
+        /**
+         * Constant `MINUIT_MINOS="MINOS"`
+         */
+        const val MINUIT_MINOS = "MINOS" //MINOS errors
+
+        /**
+         * Constant `MINUIT_HESSE="HESSE"`
+         */
+        const val MINUIT_HESSE = "HESSE" //HESSE errors
+
+        /**
+         * Constant `MINUIT_ENGINE_NAME="MINUIT"`
+         */
+        const val MINUIT_ENGINE_NAME = "MINUIT"
+    }
+}

kmath-stat/src/commonMain/tmp/minuit/MINUITPlugin.kt

@@ -0,0 +1,86 @@
+/*
+ * Copyright 2018-2021 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.optimization.minuit
+
+import hep.dataforge.context.*
+
+/**
+ * Мэнеджер для MINUITа. Пока не играет никакой активной роли кроме ведения
+ * внутреннего лога.
+ *
+ * @author Darksnake
+ * @version $Id: $Id
+ */
+@PluginDef(group = "hep.dataforge",
+    name = "MINUIT",
+    dependsOn = ["hep.dataforge:fitting"],
+    info = "The MINUIT fitter engine for DataForge fitting")
+class MINUITPlugin : BasicPlugin() {
+    fun attach(@NotNull context: Context?) {
+        super.attach(context)
+        clearStaticLog()
+    }
+
+    @Provides(Fitter.FITTER_TARGET)
+    fun getFitter(fitterName: String): Fitter? {
+        return if (fitterName == "MINUIT") {
+            MINUITFitter()
+        } else {
+            null
+        }
+    }
+
+    @ProvidesNames(Fitter.FITTER_TARGET)
+    fun listFitters(): List<String> {
+        return listOf("MINUIT")
+    }
+
+    fun detach() {
+        clearStaticLog()
+        super.detach()
+    }
+
+    class Factory : PluginFactory() {
+        fun build(meta: Meta?): Plugin {
+            return MINUITPlugin()
+        }
+
+        fun getType(): java.lang.Class<out Plugin?> {
+            return MINUITPlugin::class.java
+        }
+    }
+
+    companion object {
+        /**
+         * Constant `staticLog`
+         */
+        private val staticLog: Chronicle? = Chronicle("MINUIT-STATIC", Global.INSTANCE.getHistory())
+
+        /**
+         *
+         *
+         * clearStaticLog.
+         */
+        fun clearStaticLog() {
+            staticLog.clear()
+        }
+
+        /**
+         *
+         *
+         * logStatic.
+         *
+         * @param str  a [String] object.
+         * @param pars a [Object] object.
+         */
+        fun logStatic(str: String?, vararg pars: Any?) {
+            checkNotNull(staticLog) { "MINUIT log is not initialized." }
+            staticLog.report(str, pars)
+            LoggerFactory.getLogger("MINUIT").info(String.format(str, *pars))
+            //        Out.out.printf(str,pars);
+//        Out.out.println();
+        }
+    }
+}

kmath-stat/src/commonMain/tmp/minuit/MINUITUtils.kt

@@ -0,0 +1,121 @@
+/*
+ * Copyright 2018-2021 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.optimization.minuit
+
+import hep.dataforge.MINUIT.FunctionMinimum
+
+internal object MINUITUtils {
+    fun getFcn(source: FitState, allPar: ParamSet, fitPars: Array<String>): MultiFunction {
+        return MnFunc(source, allPar, fitPars)
+    }
+
+    fun getFitParameters(set: ParamSet, fitPars: Array<String>): MnUserParameters {
+        val pars = MnUserParameters()
+        var i: Int
+        var par: Param
+        i = 0
+        while (i < fitPars.size) {
+            par = set.getByName(fitPars[i])
+            pars.add(fitPars[i], par.getValue(), par.getErr())
+            if (par.getLowerBound() > Double.NEGATIVE_INFINITY && par.getUpperBound() < Double.POSITIVE_INFINITY) {
+                pars.setLimits(i, par.getLowerBound(), par.getUpperBound())
+            } else if (par.getLowerBound() > Double.NEGATIVE_INFINITY) {
+                pars.setLowerLimit(i, par.getLowerBound())
+            } else if (par.getUpperBound() < Double.POSITIVE_INFINITY) {
+                pars.setUpperLimit(i, par.getUpperBound())
+            }
+            i++
+        }
+        return pars
+    }
+
+    fun getValueSet(allPar: ParamSet, names: Array<String>, values: DoubleArray): ParamSet {
+        assert(values.size == names.size)
+        assert(allPar.getNames().contains(names))
+        val vector: ParamSet = allPar.copy()
+        for (i in values.indices) {
+            vector.setParValue(names[i], values[i])
+        }
+        return vector
+    }
+
+    fun isValidArray(ar: DoubleArray): Boolean {
+        for (i in ar.indices) {
+            if (java.lang.Double.isNaN(ar[i])) {
+                return false
+            }
+        }
+        return true
+    }
+
+    /**
+     *
+     *
+     * printMINUITResult.
+     *
+     * @param out a [PrintWriter] object.
+     * @param minimum a [hep.dataforge.MINUIT.FunctionMinimum] object.
+     */
+    fun printMINUITResult(out: PrintWriter, minimum: FunctionMinimum?) {
+        out.println()
+        out.println("***MINUIT INTERNAL FIT INFORMATION***")
+        out.println()
+        MnPrint.print(out, minimum)
+        out.println()
+        out.println("***END OF MINUIT INTERNAL FIT INFORMATION***")
+        out.println()
+    }
+
+    internal class MnFunc(source: FitState, allPar: ParamSet, fitPars: Array<String>) : MultiFunction {
+        var source: FitState
+        var allPar: ParamSet
+        var fitPars: Array<String>
+        fun value(doubles: DoubleArray): Double {
+            assert(isValidArray(doubles))
+            assert(doubles.size == fitPars.size)
+            return -2 * source.getLogProb(getValueSet(allPar, fitPars, doubles))
+            //                    source.getChi2(getValueSet(allPar, fitPars, doubles));
+        }
+
+        @Throws(NotDefinedException::class)
+        fun derivValue(n: Int, doubles: DoubleArray): Double {
+            assert(isValidArray(doubles))
+            assert(doubles.size == getDimension())
+            val set: ParamSet = getValueSet(allPar, fitPars, doubles)
+
+//            double res;
+//            double d, s, deriv;
+//
+//            res = 0;
+//            for (int i = 0; i < source.getDataNum(); i++) {
+//                d = source.getDis(i, set);
+//                s = source.getDispersion(i, set);
+//                if (source.modelProvidesDerivs(fitPars[n])) {
+//                    deriv = source.getDisDeriv(fitPars[n], i, set);
+//                } else {
+//                    throw new NotDefinedException();
+//                    // Такого не должно быть, поскольку мы где-то наверху должы были проверить, что производные все есть.
+//                }
+//                res += 2 * d * deriv / s;
+//            }
+            return -2 * source.getLogProbDeriv(fitPars[n], set)
+        }
+
+        fun getDimension(): Int {
+            return fitPars.size
+        }
+
+        fun providesDeriv(n: Int): Boolean {
+            return source.modelProvidesDerivs(fitPars[n])
+        }
+
+        init {
+            this.source = source
+            this.allPar = allPar
+            this.fitPars = fitPars
+            assert(source.getModel().getNames().contains(fitPars))
+        }
+    }
+}

kmath-stat/src/commonMain/tmp/minuit/MinimumBuilder.kt

@@ -0,0 +1,45 @@
+/* 
+ * Copyright 2015 Alexander Nozik.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package ru.inr.mass.minuit
+
+import space.kscience.kmath.optimization.minuit.MinimumSeed
+
+/**
+ *
+ * @version $Id$
+ */
+interface MinimumBuilder {
+    /**
+     *
+     * minimum.
+     *
+     * @param fcn a [hep.dataforge.MINUIT.MnFcn] object.
+     * @param gc a [hep.dataforge.MINUIT.GradientCalculator] object.
+     * @param seed a [hep.dataforge.MINUIT.MinimumSeed] object.
+     * @param strategy a [hep.dataforge.MINUIT.MnStrategy] object.
+     * @param maxfcn a int.
+     * @param toler a double.
+     * @return a [hep.dataforge.MINUIT.FunctionMinimum] object.
+     */
+    fun minimum(
+        fcn: MnFcn?,
+        gc: GradientCalculator?,
+        seed: MinimumSeed?,
+        strategy: MnStrategy?,
+        maxfcn: Int,
+        toler: Double
+    ): FunctionMinimum
+}

kmath-stat/src/commonMain/tmp/minuit/MinimumError.kt

@@ -0,0 +1,155 @@
+/* 
+ * Copyright 2015 Alexander Nozik.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package ru.inr.mass.minuit
+
+import space.kscience.kmath.optimization.minuit.MINUITPlugin
+
+/**
+ * MinimumError keeps the inverse 2nd derivative (inverse Hessian) used for
+ * calculating the parameter step size (-V*g) and for the covariance update
+ * (ErrorUpdator). The covariance matrix is equal to twice the inverse Hessian.
+ *
+ * @version $Id$
+ */
+class MinimumError {
+    private var theAvailable = false
+    private var theDCovar: Double
+    private var theHesseFailed = false
+    private var theInvertFailed = false
+    private var theMadePosDef = false
+    private var theMatrix: MnAlgebraicSymMatrix
+    private var thePosDef = false
+    private var theValid = false
+
+    constructor(n: Int) {
+        theMatrix = MnAlgebraicSymMatrix(n)
+        theDCovar = 1.0
+    }
+
+    constructor(mat: MnAlgebraicSymMatrix, dcov: Double) {
+        theMatrix = mat
+        theDCovar = dcov
+        theValid = true
+        thePosDef = true
+        theAvailable = true
+    }
+
+    constructor(mat: MnAlgebraicSymMatrix, x: MnHesseFailed?) {
+        theMatrix = mat
+        theDCovar = 1.0
+        theValid = false
+        thePosDef = false
+        theMadePosDef = false
+        theHesseFailed = true
+        theInvertFailed = false
+        theAvailable = true
+    }
+
+    constructor(mat: MnAlgebraicSymMatrix, x: MnMadePosDef?) {
+        theMatrix = mat
+        theDCovar = 1.0
+        theValid = false
+        thePosDef = false
+        theMadePosDef = true
+        theHesseFailed = false
+        theInvertFailed = false
+        theAvailable = true
+    }
+
+    constructor(mat: MnAlgebraicSymMatrix, x: MnInvertFailed?) {
+        theMatrix = mat
+        theDCovar = 1.0
+        theValid = false
+        thePosDef = true
+        theMadePosDef = false
+        theHesseFailed = false
+        theInvertFailed = true
+        theAvailable = true
+    }
+
+    constructor(mat: MnAlgebraicSymMatrix, x: MnNotPosDef?) {
+        theMatrix = mat
+        theDCovar = 1.0
+        theValid = false
+        thePosDef = false
+        theMadePosDef = false
+        theHesseFailed = false
+        theInvertFailed = false
+        theAvailable = true
+    }
+
+    fun dcovar(): Double {
+        return theDCovar
+    }
+
+    fun hesseFailed(): Boolean {
+        return theHesseFailed
+    }
+
+    fun hessian(): MnAlgebraicSymMatrix {
+        return try {
+            val tmp: MnAlgebraicSymMatrix = theMatrix.copy()
+            tmp.invert()
+            tmp
+        } catch (x: SingularMatrixException) {
+            MINUITPlugin.logStatic("BasicMinimumError inversion fails; return diagonal matrix.")
+            val tmp = MnAlgebraicSymMatrix(theMatrix.nrow())
+            var i = 0
+            while (i < theMatrix.nrow()) {
+                tmp[i, i] = 1.0 / theMatrix[i, i]
+                i++
+            }
+            tmp
+        }
+    }
+
+    fun invHessian(): MnAlgebraicSymMatrix {
+        return theMatrix
+    }
+
+    fun invertFailed(): Boolean {
+        return theInvertFailed
+    }
+
+    fun isAccurate(): Boolean {
+        return theDCovar < 0.1
+    }
+
+    fun isAvailable(): Boolean {
+        return theAvailable
+    }
+
+    fun isMadePosDef(): Boolean {
+        return theMadePosDef
+    }
+
+    fun isPosDef(): Boolean {
+        return thePosDef
+    }
+
+    fun isValid(): Boolean {
+        return theValid
+    }
+
+    fun matrix(): MnAlgebraicSymMatrix {
+        return MnUtils.mul(theMatrix, 2)
+    }
+
+    internal class MnHesseFailed
+    internal class MnInvertFailed
+    internal class MnMadePosDef
+    internal class MnNotPosDef
+}

kmath-stat/src/commonMain/tmp/minuit/MinimumErrorUpdator.kt

@@ -0,0 +1,33 @@
+/* 
+ * Copyright 2015 Alexander Nozik.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package ru.inr.mass.minuit
+
+/**
+ *
+ * @version $Id$
+ */
+internal interface MinimumErrorUpdator {
+    /**
+     *
+     * update.
+     *
+     * @param state a [hep.dataforge.MINUIT.MinimumState] object.
+     * @param par a [hep.dataforge.MINUIT.MinimumParameters] object.
+     * @param grad a [hep.dataforge.MINUIT.FunctionGradient] object.
+     * @return a [hep.dataforge.MINUIT.MinimumError] object.
+     */
+    fun update(state: MinimumState?, par: MinimumParameters?, grad: FunctionGradient?): MinimumError?
+}

kmath-stat/src/commonMain/tmp/minuit/MinimumParameters.kt

@@ -0,0 +1,70 @@
+/* 
+ * Copyright 2015 Alexander Nozik.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package ru.inr.mass.minuit
+
+import org.apache.commons.math3.linear.ArrayRealVector
+
+/**
+ *
+ * @version $Id$
+ */
+class MinimumParameters {
+    private var theFVal = 0.0
+    private var theHasStep = false
+    private var theParameters: RealVector
+    private var theStepSize: RealVector
+    private var theValid = false
+
+    constructor(n: Int) {
+        theParameters = ArrayRealVector(n)
+        theStepSize = ArrayRealVector(n)
+    }
+
+    constructor(avec: RealVector, fval: Double) {
+        theParameters = avec
+        theStepSize = ArrayRealVector(avec.getDimension())
+        theFVal = fval
+        theValid = true
+    }
+
+    constructor(avec: RealVector, dirin: RealVector, fval: Double) {
+        theParameters = avec
+        theStepSize = dirin
+        theFVal = fval
+        theValid = true
+        theHasStep = true
+    }
+
+    fun dirin(): RealVector {
+        return theStepSize
+    }
+
+    fun fval(): Double {
+        return theFVal
+    }
+
+    fun hasStepSize(): Boolean {
+        return theHasStep
+    }
+
+    fun isValid(): Boolean {
+        return theValid
+    }
+
+    fun vec(): RealVector {
+        return theParameters
+    }
+}

kmath-stat/src/commonMain/tmp/minuit/MinimumSeed.kt

@@ -0,0 +1,66 @@
+/* 
+ * Copyright 2015 Alexander Nozik.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package space.kscience.kmath.optimization.minuit
+
+import ru.inr.mass.minuit.*
+
+/**
+ *
+ * @version $Id$
+ */
+class MinimumSeed(state: MinimumState, trafo: MnUserTransformation) {
+    private val theState: MinimumState = state
+    private val theTrafo: MnUserTransformation = trafo
+    private val theValid: Boolean = true
+    val edm: Double get() = state().edm()
+
+    fun error(): MinimumError {
+        return state().error()
+    }
+
+    fun fval(): Double {
+        return state().fval()
+    }
+
+    fun gradient(): FunctionGradient {
+        return state().gradient()
+    }
+
+    fun isValid(): Boolean {
+        return theValid
+    }
+
+    fun nfcn(): Int {
+        return state().nfcn()
+    }
+
+    fun parameters(): MinimumParameters {
+        return state().parameters()
+    }
+
+    fun precision(): MnMachinePrecision {
+        return theTrafo.precision()
+    }
+
+    fun state(): MinimumState {
+        return theState
+    }
+
+    fun trafo(): MnUserTransformation {
+        return theTrafo
+    }
+
+}

kmath-stat/src/commonMain/tmp/minuit/MinimumSeedGenerator.kt

@@ -0,0 +1,39 @@
+/* 
+ * Copyright 2015 Alexander Nozik.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package ru.inr.mass.minuit
+
+import space.kscience.kmath.optimization.minuit.MinimumSeed
+
+/**
+ * base class for seed generators (starting values); the seed generator prepares
+ * initial starting values from the input (MnUserParameterState) for the
+ * minimization;
+ *
+ * @version $Id$
+ */
+interface MinimumSeedGenerator {
+    /**
+     *
+     * generate.
+     *
+     * @param fcn a [hep.dataforge.MINUIT.MnFcn] object.
+     * @param calc a [hep.dataforge.MINUIT.GradientCalculator] object.
+     * @param user a [hep.dataforge.MINUIT.MnUserParameterState] object.
+     * @param stra a [hep.dataforge.MINUIT.MnStrategy] object.
+     * @return a [hep.dataforge.MINUIT.MinimumSeed] object.
+     */
+    fun generate(fcn: MnFcn?, calc: GradientCalculator?, user: MnUserParameterState?, stra: MnStrategy?): MinimumSeed
+}