QOW is working more or less

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,19 +3,17 @@
  * 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.commons.fit
+package space.kscience.kmath.fit
 
 import kotlinx.html.br
 import kotlinx.html.h3
 import space.kscience.kmath.commons.expressions.DSProcessor
 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.FunctionOptimizationTarget
+import space.kscience.kmath.optimization.*
-import space.kscience.kmath.optimization.optimizeWith
-import space.kscience.kmath.optimization.resultPoint
-import space.kscience.kmath.optimization.resultValue
 import space.kscience.kmath.real.DoubleVector
 import space.kscience.kmath.real.map
 import space.kscience.kmath.real.step
@@ -25,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
 
@@ -45,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))
@@ -56,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)
     }

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

@@ -0,0 +1,105 @@
+/*
+ * 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.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 + 100
+        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 1.0, b to 1.2, c to 99.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.resultValue / (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)

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/optimization/CMOptimizer.kt

@@ -18,6 +18,7 @@ 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
@@ -108,15 +109,17 @@ public object CMOptimizer : Optimizer<Double, FunctionOptimization<Double>> {
 
             val objectiveFunction = ObjectiveFunction {
                 val args = startPoint + it.toMap()
-                problem.expression(args)
+                val res = problem.expression(args)
+                res
             }
             addOptimizationData(objectiveFunction)
 
             val gradientFunction = ObjectiveFunctionGradient {
                 val args = startPoint + it.toMap()
-                DoubleArray(symbols.size) { index ->
+                val res = DoubleArray(symbols.size) { index ->
                     problem.expression.derivative(symbols[index])(args)
                 }
+                res
             }
             addOptimizationData(gradientFunction)
 

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

@@ -15,6 +15,7 @@ import space.kscience.kmath.expressions.chiSquaredExpression
 import space.kscience.kmath.expressions.symbol
 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
@@ -58,7 +59,7 @@ internal class OptimizeTest {
             it.pow(2) + it + 1 + chain.next()
         }
 
-        val yErr = List(x.size) { sigma }.asBuffer()
+        val yErr = DoubleBuffer(x.size) { sigma }
 
         val chi2 = DSProcessor.chiSquaredExpression(
             x, y, yErr

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

@@ -32,11 +32,10 @@ public interface XYColumnarData<out T, out X : T, out Y : T> : ColumnarData<T> {
         Symbol.y -> y
         else -> null
     }
-}
 
-@Suppress("FunctionName")
+    public companion object{
         @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}" }
             return object : XYColumnarData<T, X, Y> {
                 override val size: Int = x.size
@@ -44,6 +43,9 @@ public fun <T, X : T, Y : T> XYColumnarData(x: Buffer<X>, y: Buffer<Y>): XYColum
                 override val y: Buffer<Y> = y
             }
         }
+    }
+}
+
 
 /**
  * Represent a [ColumnarData] as an [XYColumnarData]. The presence or respective columns is checked on creation.

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

@@ -5,15 +5,13 @@
 
 package space.kscience.kmath.data
 
-import space.kscience.kmath.data.XYErrorColumnarData.Companion
 import space.kscience.kmath.expressions.Symbol
-import space.kscience.kmath.expressions.symbol
 import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.structures.Buffer
 
 
 /**
- * A [ColumnarData] with additional [Companion.yErr] column for an [Symbol.y] error
+ * A [ColumnarData] with additional [Symbol.yError] column for an [Symbol.y] error
  * Inherits [XYColumnarData].
  */
 @UnstableKMathAPI
@@ -23,11 +21,24 @@ public interface XYErrorColumnarData<T, out X : T, out Y : T> : XYColumnarData<T
     override fun get(symbol: Symbol): Buffer<T> = when (symbol) {
         Symbol.x -> x
         Symbol.y -> y
-        Companion.yErr -> yErr
+        Symbol.yError -> yErr
         else -> error("A column for symbol $symbol not found")
     }
 
     public companion object {
-        public val yErr: Symbol by symbol
+        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/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>
+    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/specialExpressions.kt

@@ -7,13 +7,18 @@ 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.
  */
-public fun <T : Any, I : Any, A> AutoDiffProcessor<T, I, A>.chiSquaredExpression(
+@JvmName("genericChiSquaredExpression")
+public fun <T : Comparable<T>, I : Any, A> AutoDiffProcessor<T, I, A>.chiSquaredExpression(
     x: Buffer<T>,
     y: Buffer<T>,
     yErr: Buffer<T>,
@@ -36,3 +41,13 @@ public fun <T : Any, I : Any, A> AutoDiffProcessor<T, I, A>.chiSquaredExpression
         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/misc/Featured.kt

@@ -5,6 +5,7 @@
 
 package space.kscience.kmath.misc
 
+import kotlin.jvm.JvmInline
 import kotlin.reflect.KClass
 
 /**
@@ -29,7 +30,8 @@ public interface Feature<F : Feature<F>> {
 /**
  * A container for a set of features
  */
-public class FeatureSet<F : Feature<F>> private constructor(public val features: Map<FeatureKey<F>, F>) : Featured<F> {
+@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 }
 
@@ -48,6 +50,9 @@ public class FeatureSet<F : Feature<F>> private constructor(public val features:
 
     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> =

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

@@ -5,10 +5,18 @@
 
 package space.kscience.kmath.misc
 
-public interface Loggable {
+import space.kscience.kmath.misc.Loggable.Companion.INFO
-    public fun log(tag: String = INFO, block: () -> String)
+
+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/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-functions/src/commonMain/kotlin/space/kscience/kmath/integration/GaussIntegrator.kt

@@ -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)

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-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

@@ -5,11 +5,9 @@
 
 package space.kscience.kmath.optimization
 
-import space.kscience.kmath.expressions.*
+import space.kscience.kmath.expressions.DifferentiableExpression
+import space.kscience.kmath.expressions.Symbol
 import space.kscience.kmath.misc.FeatureSet
-import space.kscience.kmath.operations.ExtendedField
-import space.kscience.kmath.structures.Buffer
-import space.kscience.kmath.structures.indices
 
 public class OptimizationValue<T>(public val value: T) : OptimizationFeature {
     override fun toString(): String = "Value($value)"
@@ -25,7 +23,26 @@ public class FunctionOptimization<T>(
     public val expression: DifferentiableExpression<T>,
 ) : OptimizationProblem<T> {
 
-    public companion object
+
+    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(

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

@@ -72,15 +72,15 @@ public suspend fun <T> DifferentiableExpression<T>.optimizeWith(
 public class XYOptimizationBuilder(
     public val data: XYColumnarData<Double, Double, Double>,
     public val model: DifferentiableExpression<Double>,
-) : OptimizationBuilder<Double, XYOptimization>() {
+) : OptimizationBuilder<Double, XYFit>() {
 
     public var pointToCurveDistance: PointToCurveDistance = PointToCurveDistance.byY
     public var pointWeight: PointWeight = PointWeight.byYSigma
 
-    override fun build(): XYOptimization = XYOptimization(
+    override fun build(): XYFit = XYFit(
-        FeatureSet.of(features),
         data,
         model,
+        FeatureSet.of(features),
         pointToCurveDistance,
         pointWeight
     )
@@ -90,4 +90,4 @@ public fun XYOptimization(
     data: XYColumnarData<Double, Double, Double>,
     model: DifferentiableExpression<Double>,
     builder: XYOptimizationBuilder.() -> Unit,
-): XYOptimization = XYOptimizationBuilder(data, model).apply(builder).build()
+): XYFit = XYOptimizationBuilder(data, model).apply(builder).build()

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

@@ -11,6 +11,7 @@ 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
@@ -21,14 +22,14 @@ import space.kscience.kmath.structures.DoubleL2Norm
  * See [the article](http://arxiv.org/abs/physics/0604127).
  */
 @UnstableKMathAPI
-public class QowOptimizer : Optimizer<Double, XYOptimization> {
+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 inner class QoWeight(
+    private class QoWeight(
-        val problem: XYOptimization,
+        val problem: XYFit,
         val parameters: Map<Symbol, Double>,
     ) : Map<Symbol, Double> by parameters, SymbolIndexer {
         override val symbols: List<Symbol> = parameters.keys.toList()
@@ -39,8 +40,8 @@ public class QowOptimizer : Optimizer<Double, XYOptimization> {
          * 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) { i, k ->
+            linearSpace.buildMatrix(problem.data.size, symbols.size) { d, s ->
-                problem.distance(i).derivative(symbols[k])(parameters)
+                problem.distance(d).derivative(symbols[s])(parameters)
             }
         }
 
@@ -48,25 +49,27 @@ public class QowOptimizer : Optimizer<Double, XYOptimization> {
          * Array of dispersions in each point
          */
         val dispersion: Point<Double> by lazy {
-            DoubleBuffer(problem.data.size) { i ->
+            DoubleBuffer(problem.data.size) { d ->
-                problem.weight(i).invoke(parameters)
+                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 [i]-th point of data.
+     * The signed distance from the model to the [d]-th point of data.
      */
-    private fun QoWeight.distance(i: Int, parameters: Map<Symbol, Double>): Double = problem.distance(i)(parameters)
+    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, i: Int, parameters: Map<Symbol, Double>): Double =
+    private fun QoWeight.distanceDerivative(symbol: Symbol, d: Int, parameters: Map<Symbol, Double>): Double =
-        problem.distance(i).derivative(symbol)(parameters)
+        problem.distance(d).derivative(symbol)(parameters)
 
     /**
      * Теоретическая ковариация весовых функций.
@@ -74,8 +77,8 @@ public class QowOptimizer : Optimizer<Double, XYOptimization> {
      * 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 { k, l ->
+        linearSpace.matrix(size, size).symmetric { s1, s2 ->
-            (0 until data.size).sumOf { i -> derivs[k, i] * derivs[l, i] / dispersion[i] }
+            (0 until data.size).sumOf { d -> derivs[d, s1] * derivs[d, s2] / dispersion[d] }
         }
 
     /**
@@ -89,12 +92,12 @@ public class QowOptimizer : Optimizer<Double, XYOptimization> {
              * количество вызывов функции будет dim^2 вместо dim Первый индекс -
              * номер точки, второй - номер переменной, по которой берется производная
              */
-            val eqvalues = linearSpace.buildMatrix(data.size, size) { i, l ->
+            val eqvalues = linearSpace.buildMatrix(data.size, size) { d, s ->
-                distance(i, theta) * derivs[l, i] / dispersion[i]
+                distance(d, theta) * derivs[d, s] / dispersion[d]
             }
 
-            buildMatrix(size, size) { k, l ->
+            buildMatrix(size, size) { s1, s2 ->
-                (0 until data.size).sumOf { i -> eqvalues[i, l] * eqvalues[i, k] }
+                (0 until data.size).sumOf { d -> eqvalues[d, s2] * eqvalues[d, s1] }
             }
         }
 
@@ -106,20 +109,20 @@ public class QowOptimizer : Optimizer<Double, XYOptimization> {
     ): Matrix<Double> = with(linearSpace) {
         //Возвращает производную k-того Eq по l-тому параметру
         //val res = Array(fitDim) { DoubleArray(fitDim) }
-        val sderiv = buildMatrix(data.size, size) { i, l ->
+        val sderiv = buildMatrix(data.size, size) { d, s ->
-            distanceDerivative(symbols[l], i, theta)
+            distanceDerivative(symbols[s], d, theta)
         }
 
-        buildMatrix(size, size) { k, l ->
+        buildMatrix(size, size) { s1, s2 ->
-            val base = (0 until data.size).sumOf { i ->
+            val base = (0 until data.size).sumOf { d ->
-                require(dispersion[i] > 0)
+                require(dispersion[d] > 0)
-                sderiv[i, l] * derivs[k, i] / dispersion[i]
+                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[k])(theta)
+                val deriv1 = prior.derivative(symbols[s1])(theta)
-                val deriv2 = prior.derivative(symbols[l])(theta)
+                val deriv2 = prior.derivative(symbols[s2])(theta)
                 base + deriv1 * deriv2 / pi / pi
             } ?: base
         }
@@ -130,13 +133,13 @@ public class QowOptimizer : Optimizer<Double, XYOptimization> {
      * Значения уравнений метода квазиоптимальных весов
      */
     private fun QoWeight.getEqValues(theta: Map<Symbol, Double> = this): Point<Double> {
-        val distances = DoubleBuffer(data.size) { i -> distance(i, theta) }
+        val distances = DoubleBuffer(data.size) { d -> distance(d, theta) }
 
-        return DoubleBuffer(size) { k ->
+        return DoubleBuffer(size) { s ->
-            val base = (0 until data.size).sumOf { i -> distances[i] * derivs[k, i] / dispersion[i] }
+            val base = (0 until data.size).sumOf { d -> distances[d] * derivs[d, s] / dispersion[d] }
             //Поправка на априорную вероятность
             prior?.let { prior ->
-                base - prior.derivative(symbols[k])(theta) / prior(theta)
+                base - prior.derivative(symbols[s])(theta) / prior(theta)
             } ?: base
         }
     }
@@ -163,15 +166,15 @@ public class QowOptimizer : Optimizer<Double, XYOptimization> {
 
         val logger = problem.getFeature<OptimizationLog>()
 
-        var dis: Double//норма невязки
+        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)//значения функций
+        var eqvalues = getEqValues(par) //Values of the weight functions
 
-        dis = DoubleL2Norm.norm(eqvalues)// невязка
+        dis = DoubleL2Norm.norm(eqvalues) // discrepancy
         logger?.log { "Starting discrepancy is $dis" }
         var i = 0
         var flag = false
@@ -238,7 +241,8 @@ public class QowOptimizer : Optimizer<Double, XYOptimization> {
         return covar
     }
 
-    override suspend fun optimize(problem: XYOptimization): XYOptimization {
+    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,125 @@
+/*
+ * 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.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,
+) : 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
+    )
+    return optimizer.optimize(problem)
+}

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)(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/XYOptimization.kt

@@ -1,189 +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/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.DifferentiableExpression
-import space.kscience.kmath.expressions.Expression
-import space.kscience.kmath.expressions.Symbol
-import space.kscience.kmath.expressions.derivative
-import space.kscience.kmath.misc.FeatureSet
-import space.kscience.kmath.misc.UnstableKMathAPI
-import kotlin.math.PI
-import kotlin.math.ln
-import kotlin.math.pow
-import kotlin.math.sqrt
-
-/**
- * Specify the way to compute distance from point to the curve as DifferentiableExpression
- */
-public interface PointToCurveDistance : OptimizationFeature {
-    public fun distance(problem: XYOptimization, index: Int): DifferentiableExpression<Double>
-
-    public companion object {
-        public val byY: PointToCurveDistance = object : PointToCurveDistance {
-            override fun distance(problem: XYOptimization, 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)
-
-                    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: XYOptimization, index: Int): DifferentiableExpression<Double>
-
-    public companion object {
-        public fun bySigma(sigmaSymbol: Symbol): PointWeight = object : PointWeight {
-            override fun weight(problem: XYOptimization, 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)
-    }
-}
-
-/**
- * An optimization for XY data.
- */
-public class XYOptimization(
-    override val features: FeatureSet<OptimizationFeature>,
-    public val data: XYColumnarData<Double, Double, Double>,
-    public val model: DifferentiableExpression<Double>,
-    internal val pointToCurveDistance: PointToCurveDistance = PointToCurveDistance.byY,
-    internal val pointWeight: PointWeight = PointWeight.byYSigma,
-) : 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 XYOptimization.withFeature(vararg features: OptimizationFeature): XYOptimization {
-    return XYOptimization(this.features.with(*features), data, model, pointToCurveDistance, pointWeight)
-}
-
-private val oneOver2Pi = 1.0 / sqrt(2 * PI)
-
-internal fun XYOptimization.likelihood(): 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)(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.
- */
-public suspend fun Optimizer<Double, FunctionOptimization<Double>>.maximumLogLikelihood(problem: XYOptimization): XYOptimization {
-    val functionOptimization = FunctionOptimization(problem.features, problem.likelihood())
-    val result = optimize(functionOptimization.withFeatures(FunctionOptimizationTarget.MAXIMIZE))
-    return XYOptimization(result.features, problem.data, problem.model)
-}
-
-public suspend fun Optimizer<Double, FunctionOptimization<Double>>.maximumLogLikelihood(
-    data: XYColumnarData<Double, Double, Double>,
-    model: DifferentiableExpression<Double>,
-    builder: XYOptimizationBuilder.() -> Unit,
-): XYOptimization = maximumLogLikelihood(XYOptimization(data, model, builder))
-
-//public suspend fun XYColumnarData<Double, Double, Double>.fitWith(
-//    optimizer: XYOptimization,
-//    problemBuilder: XYOptimizationBuilder.() -> Unit = {},
-//
-//)
-
-
-//
-//@UnstableKMathAPI
-//public interface XYFit<T> : OptimizationProblem {
-//
-//    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)) }
-//    }
-//}
-
-//
-///**
-// * Define 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)
-//    function(chiSquared)
-//    maximize = false
-//}

settings.gradle.kts

@@ -26,6 +26,7 @@ include(
     ":kmath-histograms",
     ":kmath-commons",
     ":kmath-viktor",
+    ":kmath-optimization",
     ":kmath-stat",
     ":kmath-nd4j",
     ":kmath-dimensions",