[WIP] porting model

build.gradle.kts

@@ -12,8 +12,8 @@ allprojects {
     version = "0.1.0-SNAPSHOT"
 }
 
-val dataforgeVersion by extra("0.4.0-dev-2")
+val dataforgeVersion by extra("0.4.0")
-val kmathVersion by extra("0.3.0-dev-3")
+val kmathVersion by extra("0.3.0-dev-8")
 
 apiValidation{
     validationDisabled = true

docs/diagrams/online.puml

@@ -0,0 +1,27 @@
+@startuml
+participant Control
+participant Detector
+participant HV
+
+Control -> Detector: check connection
+
+Control -> HV: set voltage
+activate HV
+HV -> Control: voltage set
+deactivate HV
+
+Control -> Detector: start measurement
+activate Detector
+
+Control -->o Detector: cancel measurement
+
+Detector --> Control: asynchronous info
+
+Detector -> Control: measurement result (binary)
+deactivate Detector
+
+HV -> Control: dump HV log
+
+Control -> Control: combine and store data
+
+@enduml

numass-model/build.gradle.kts

@@ -13,6 +13,7 @@ kotlin.sourceSets {
         dependencies {
             api("space.kscience:dataforge-meta:$dataforgeVersion")
             api("space.kscience:kmath-for-real:$kmathVersion")
+            api("space.kscience:kmath-functions:$kmathVersion")
         }
     }
     jvmMain{

numass-model/src/commonMain/kotlin/ru/inr/mass/models/AbstractParametric.kt

@@ -1,48 +0,0 @@
-/* 
- * 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 hep.dataforge.stat.parametric
-
-import hep.dataforge.exceptions.NameNotFoundException
-
-abstract class AbstractParametric : AbstractNamedSet {
-    constructor(names: NameList?) : super(names) {}
-    constructor(list: Array<String?>?) : super(list) {}
-    constructor(set: NameSetContainer?) : super(set) {}
-
-    /**
-     * Provide default value for parameter with name `name`. Throws
-     * NameNotFound if no default found for given parameter.
-     *
-     * @param name
-     * @return
-     */
-    protected fun getDefaultParameter(name: String?): Double {
-        throw NameNotFoundException(name)
-    }
-
-    /**
-     * Extract value from input vector using default value if required parameter
-     * not found
-     *
-     * @param name
-     * @param set
-     * @return
-     */
-    protected fun getParameter(name: String?, set: Values): Double {
-        //FIXME add default value
-        return set.getDouble(name)
-    }
-}

numass-model/src/commonMain/kotlin/ru/inr/mass/models/AbstractParametricBiFunction.kt

@@ -1,22 +0,0 @@
-/*
- * To change this license header, choose License Headers in Project Properties.
- * To change this template file, choose Tools | Templates
- * and open the template in the editor.
- */
-package hep.dataforge.stat.parametric
-
-import hep.dataforge.exceptions.NotDefinedException
-
-abstract class AbstractParametricBiFunction : AbstractParametric, ParametricBiFunction {
-    constructor(names: NameList?) : super(names) {}
-    constructor(list: Array<String?>?) : super(list) {}
-    constructor(set: NameSetContainer?) : super(set) {}
-
-    fun derivValue(parName: String?, x: Double, y: Double, set: Values?): Double {
-        return if (!getNames().contains(parName)) {
-            0
-        } else {
-            throw NotDefinedException()
-        }
-    }
-}

numass-model/src/commonMain/kotlin/ru/inr/mass/models/AbstractParametricFunction.kt

@@ -1,34 +0,0 @@
-/* 
- * 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 hep.dataforge.stat.parametric
-
-import hep.dataforge.names.NameList
-
-typealias NameList = List<String>
-
-/**
- *
- *
- * Abstract AbstractNamedSpectrum class.
- *
- * @author Alexander Nozik
- * @version $Id: $Id
- */
-abstract class AbstractParametricFunction : AbstractParametric, ParametricFunction {
-    constructor(names: NameList?) : super(names) {}
-    constructor(vararg list: String?) : super(list) {}
-    constructor(set: NameSetContainer?) : super(set) {}
-}

numass-model/src/commonMain/kotlin/ru/inr/mass/models/AbstractParametricValue.kt

@@ -1,14 +0,0 @@
-/*
- * To change this license header, choose License Headers in Project Properties.
- * To change this template file, choose Tools | Templates
- * and open the template in the editor.
- */
-package hep.dataforge.stat.parametric
-
-import hep.dataforge.names.NameList
-
-abstract class AbstractParametricValue : AbstractParametric, ParametricValue {
-    constructor(names: NameList?) : super(names) {}
-    constructor(list: Array<String?>?) : super(list) {}
-    constructor(set: NameSetContainer?) : super(set) {}
-}

numass-model/src/commonMain/kotlin/ru/inr/mass/models/NBkgSpectrum.kt

@@ -15,61 +15,32 @@
  */
 package ru.inr.mass.models
 
-import hep.dataforge.names.NamesUtils.combineNamesWithEquals
+import space.kscience.kmath.misc.Symbol
-import hep.dataforge.stat.parametric.AbstractParametricFunction
+import space.kscience.kmath.misc.symbol
-import hep.dataforge.stat.parametric.ParametricFunction
-import hep.dataforge.utils.MultiCounter
-import hep.dataforge.values.ValueProvider
-import hep.dataforge.values.Values
-import space.kscience.kmath.expressions.Symbol
-
-typealias Values = Map<Symbol, Double>
 
 /**
  *
- * @author Darksnake
  */
-open class NBkgSpectrum(private val source: ParametricFunction) : AbstractParametricFunction(*combineNamesWithEquals(source.namesAsArray(), *list)) {
+public class NBkgSpectrum(public val source: Spectrum) : DifferentiableSpectrum {
-
+    override fun invoke(x: Double, arguments: Map<Symbol, Double>): Double {
-    override fun derivValue(parName: String, x: Double, set: Values): Double {
+        val normValue = arguments[norm] ?: 1.0
-        return when (parName) {
+        val bkgValue = arguments[bkg] ?: 0.0
-            "N" -> source.value(x, set)
+        return normValue * source(x, arguments) + bkgValue
-            "bkg" -> 1.0
-            else -> getN(set) * source.derivValue(parName, x, set)
-        }
     }
 
-    private fun getBkg(set: ValueProvider): Double {
+    override fun derivativeOrNull(symbols: List<Symbol>): Spectrum? = when {
-        return set.getDouble("bkg")
+        symbols.isEmpty() -> this
+        symbols.size == 1 -> when (symbols.first()) {
+            norm -> Spectrum { x, arguments -> source(x, arguments) + (arguments[bkg] ?: 0.0) }
+            bkg -> Spectrum { x, arguments -> (arguments[norm] ?: 1.0) * source(x, arguments) }
+            else -> (source as? DifferentiableSpectrum)?.derivativeOrNull(symbols)?.let { NBkgSpectrum(it) }
+        }
+        else -> null
     }
 
-    private fun getN(set: ValueProvider): Double {
+    public companion object {
-        return set.getDouble("N")
+        public val norm: Symbol by symbol
-    }
+        public val bkg: Symbol by symbol
-
-    override fun providesDeriv(name: String): Boolean {
-        return when (name) {
-            "N","bkg" -> true
-            else -> this.source.providesDeriv(name)
-        }
-    }
-
-    override fun value(x: Double, set: Values): Double {
-        this.counter.increase("value")
-        return getN(set) * source.value(x, set) + getBkg(set)
-    }
-
-    override fun getDefaultParameter(name: String): Double {
-        return when (name) {
-            "bkg" -> 0.0
-            "N" -> 1.0
-            else -> super.getDefaultParameter(name)
-        }
-    }
-
-    companion object {
-
-        private val list = arrayOf("N", "bkg")
     }
 
 }

numass-model/src/commonMain/kotlin/ru/inr/mass/models/NumassBeta.kt

@@ -5,25 +5,20 @@
  */
 package inr.numass.models.sterile
 
-import hep.dataforge.exceptions.NotDefinedException
+import ru.inr.mass.models.DifferentiableKernel
-import hep.dataforge.stat.parametric.AbstractParametricBiFunction
+import ru.inr.mass.models.Kernel
-import hep.dataforge.stat.parametric.AbstractParametricFunction
+import space.kscience.kmath.misc.StringSymbol
-import hep.dataforge.stat.parametric.ParametricFunction
+import space.kscience.kmath.misc.Symbol
-import hep.dataforge.values.Values
+import space.kscience.kmath.misc.symbol
-
+import kotlin.math.*
-import java.lang.Math.*
 
 /**
  * A bi-function for beta-spectrum calculation taking energy and final state as
  * input.
  *
- *
  * dissertation p.33
- *
- *
- * @author [Alexander Nozik](mailto:altavir@gmail.com)
  */
-class NumassBeta : AbstractParametricBiFunction(list) {
+public class NumassBeta : DifferentiableKernel {
 
     /**
      * Beta spectrum derivative
@@ -35,7 +30,6 @@ class NumassBeta : AbstractParametricBiFunction(list) {
      * @return
      * @throws NotDefinedException
      */
-    @Throws(NotDefinedException::class)
     private fun derivRoot(n: Int, E0: Double, mnu2: Double, E: Double): Double {
         val D = E0 - E//E0-E
         if (D == 0.0) {
@@ -58,7 +52,7 @@ class NumassBeta : AbstractParametricBiFunction(list) {
             if (E >= E0 + mu) {
                 0.0
             } else {
-                val root = sqrt(Math.max(D * D - mnu2, 0.0))
+                val root = sqrt(max(D * D - mnu2, 0.0))
                 val exp = exp(-1 - D / mu)
                 when (n) {
                     0 -> factor(E) * (D * (D + mu * exp) / root + root * (1 - exp))
@@ -79,29 +73,28 @@ class NumassBeta : AbstractParametricBiFunction(list) {
      * @return
      * @throws NotDefinedException
      */
-    @Throws(NotDefinedException::class)
+    private fun derivRootsterile(symbol: Symbol, E: Double, E0: Double, pars: Map<Symbol, Double>): Double {
-    private fun derivRootsterile(name: String, E: Double, E0: Double, pars: Values): Double {
+        val mnu2Value = pars.getValue(mnu2)
-        val mnu2 = getParameter("mnu2", pars)
+        val msterile2Value = pars.getValue(msterile2)
-        val mst2 = getParameter("msterile2", pars)
+        val u2Value = pars.getValue(u2)
-        val u2 = getParameter("U2", pars)
 
-        return when (name) {
+        return when (symbol) {
-            "E0" -> {
+            e0 -> {
-                if (u2 == 0.0) {
+                if (u2Value == 0.0) {
-                    derivRoot(0, E0, mnu2, E)
+                    derivRoot(0, E0, mnu2Value, E)
                 } else {
-                    u2 * derivRoot(0, E0, mst2, E) + (1 - u2) * derivRoot(0, E0, mnu2, E)
+                    u2Value * derivRoot(0, E0, msterile2Value, E) + (1 - u2Value) * derivRoot(0, E0, mnu2Value, E)
                 }
             }
-            "mnu2" -> (1 - u2) * derivRoot(1, E0, mnu2, E)
+            mnu2 -> (1 - u2Value) * derivRoot(1, E0, mnu2Value, E)
-            "msterile2" -> {
+            msterile2 -> {
-                if (u2 == 0.0) {
+                if (u2Value == 0.0) {
                     0.0
                 } else {
-                    u2 * derivRoot(1, E0, mst2, E)
+                    u2Value * derivRoot(1, E0, msterile2Value, E)
                 }
             }
-            "U2" -> root(E0, mst2, E) - root(E0, mnu2, E)
+            u2 -> root(E0, msterile2Value, E) - root(E0, mnu2Value, E)
             else -> 0.0
         }
 
@@ -119,18 +112,14 @@ class NumassBeta : AbstractParametricBiFunction(list) {
         val eTot = E + me
         val pe = sqrt(E * (E + 2.0 * me))
         val ve = pe / eTot
-        val yfactor = 2.0 * 2.0 * 1.0 / 137.039 * Math.PI
+        val yfactor = 2.0 * 2.0 * 1.0 / 137.039 * PI
         val y = yfactor / ve
         val fn = y / abs(1.0 - exp(-y))
         val fermi = fn * (1.002037 - 0.001427 * ve)
-        val res = fermi * pe * eTot
+        val res: Double = fermi * pe * eTot
         return K * res
     }
 
-    override fun providesDeriv(name: String): Boolean {
-        return true
-    }
-
     /**
      * Bare beta spectrum with Mainz negative mass correction
      *
@@ -142,14 +131,14 @@ class NumassBeta : AbstractParametricBiFunction(list) {
     private fun root(E0: Double, mnu2: Double, E: Double): Double {
         //bare beta-spectrum
         val delta = E0 - E
-        val bare = factor(E) * delta * sqrt(Math.max(delta * delta - mnu2, 0.0))
+        val bare = factor(E) * delta * sqrt(max(delta * delta - mnu2, 0.0))
         return when {
-            mnu2 >= 0 -> Math.max(bare, 0.0)
+            mnu2 >= 0 -> max(bare, 0.0)
             delta == 0.0 -> 0.0
             delta + 0.812 * sqrt(-mnu2) <= 0 -> 0.0              //sqrt(0.66)
             else -> {
                 val aux = sqrt(-mnu2 * 0.66) / delta
-                Math.max(bare * (1 + aux * exp(-1 - 1 / aux)), 0.0)
+                max(bare * (1 + aux * exp(-1 - 1 / aux)), 0.0)
             }
         }
     }
@@ -162,10 +151,10 @@ class NumassBeta : AbstractParametricBiFunction(list) {
      * @param pars
      * @return
      */
-    private fun rootsterile(E: Double, E0: Double, pars: Values): Double {
+    private fun rootsterile(E: Double, E0: Double, pars: Map<Symbol, Double>): Double {
-        val mnu2 = getParameter("mnu2", pars)
+        val mnu2 = pars.getValue(mnu2)
-        val mst2 = getParameter("msterile2", pars)
+        val mst2 = pars.getValue(msterile2)
-        val u2 = getParameter("U2", pars)
+        val u2 = pars.getValue(u2)
 
         return if (u2 == 0.0) {
             root(E0, mnu2, E)
@@ -175,51 +164,65 @@ class NumassBeta : AbstractParametricBiFunction(list) {
 // P(rootsterile)+ (1-P)root
     }
 
-    override fun getDefaultParameter(name: String): Double {
+    override val x: Symbol = StringSymbol("fs")
-        return when (name) {
+    override val y: Symbol = StringSymbol("eIn")
-            "mnu2", "U2", "msterile2" -> 0.0
+
-            else -> super.getDefaultParameter(name)
+    override fun invoke(x: Double, y: Double, arguments: Map<Symbol, Double>): Double {
-        }
+        val e0 = arguments.getValue(e0)
+        return rootsterile(y, e0 - x, arguments)
     }
 
-    override fun derivValue(parName: String, fs: Double, eIn: Double, pars: Values): Double {
+    override fun derivativeOrNull(symbols: List<Symbol>): Kernel? = when (symbols.size) {
-        val e0 = getParameter("E0", pars)
+        0 -> this
-        return derivRootsterile(parName, eIn, e0 - fs, pars)
+        1 -> Kernel { fs, eIn, arguments ->
+            val e0 = arguments.getValue(e0)
+            derivRootsterile(symbols.first(), eIn, e0 - fs, arguments)
         }
-
+        else -> null
-    override fun value(fs: Double, eIn: Double, pars: Values): Double {
-        val e0 = getParameter("E0", pars)
-        return rootsterile(eIn, e0 - fs, pars)
-    }
-
-    /**
-     * Get univariate spectrum with given final state
-     */
-    fun getSpectrum(fs: Double = 0.0): ParametricFunction {
-        return BetaSpectrum(fs);
-    }
-
-    inner class BetaSpectrum(val fs: Double) : AbstractParametricFunction(*list) {
-
-        override fun providesDeriv(name: String): Boolean {
-            return this@NumassBeta.providesDeriv(name)
-        }
-
-        override fun derivValue(parName: String, x: Double, set: Values): Double {
-            return this@NumassBeta.derivValue(parName, fs, x, set)
-        }
-
-        override fun value(x: Double, set: Values): Double {
-            return this@NumassBeta.value(fs, x, set)
-        }
-
     }
+//
+//    override fun getDefaultParameter(name: String): Double {
+//        return when (name) {
+//            "mnu2", "U2", "msterile2" -> 0.0
+//            else -> super.getDefaultParameter(name)
+//        }
+//    }
+//
+//    override fun derivValue(parName: String, fs: Double, eIn: Double, pars: Values): Double {
+//        val e0 = getParameter("E0", pars)
+//        return derivRootsterile(parName, eIn, e0 - fs, pars)
+//    }
+//
+//    /**
+//     * Get univariate spectrum with given final state
+//     */
+//    fun getSpectrum(fs: Double = 0.0): ParametricFunction {
+//        return BetaSpectrum(fs);
+//    }
+//
+//    inner class BetaSpectrum(val fs: Double) : AbstractParametricFunction(*list) {
+//
+//        override fun providesDeriv(name: String): Boolean {
+//            return this@NumassBeta.providesDeriv(name)
+//        }
+//
+//        override fun derivValue(parName: String, x: Double, set: Values): Double {
+//            return this@NumassBeta.derivValue(parName, fs, x, set)
+//        }
+//
+//        override fun value(x: Double, set: Values): Double {
+//            return this@NumassBeta.value(fs, x, set)
+//        }
+//
+//    }
 
 
-    companion object {
+    public companion object {
-
+        private const val K: Double = 1E-23
-        private const val K = 1E-23
+        public val e0: Symbol by symbol
-        private val list = arrayOf("E0", "mnu2", "msterile2", "U2")
+        public val mnu2: Symbol by symbol
+        public val msterile2: Symbol by symbol
+        public val u2: Symbol by symbol
     }
 
 }

numass-model/src/commonMain/kotlin/ru/inr/mass/models/NumassResolution.kt

@@ -5,83 +5,78 @@
  */
 package inr.numass.models.sterile
 
-import hep.dataforge.context.Context
+import ru.inr.mass.models.BivariateFunction
-import hep.dataforge.maths.functions.FunctionLibrary
+import ru.inr.mass.models.DifferentiableKernel
-import hep.dataforge.meta.Meta
+import ru.inr.mass.models.Kernel
-import hep.dataforge.stat.parametric.AbstractParametricBiFunction
+import space.kscience.kmath.misc.Symbol
-import hep.dataforge.values.Values
+import space.kscience.kmath.misc.symbol
-import inr.numass.models.ResolutionFunction
+import kotlin.math.sqrt
-import inr.numass.utils.ExpressionUtils
-import org.apache.commons.math3.analysis.BivariateFunction
-import java.lang.Math.sqrt
-import java.util.*
 
 /**
  * @author [Alexander Nozik](mailto:altavir@gmail.com)
  */
-class NumassResolution(context: Context, meta: Meta) : AbstractParametricBiFunction(list) {
+public class NumassResolution(
+    public val resA: Double = 8.3e-5,
+    public val resB: Double = 0.0,
+    public val tailFunction: BivariateFunction =  { _, _ -> 1.0 },
+) : DifferentiableKernel {
 
-    private val resA: Double = meta.getDouble("A", 8.3e-5)
+//    private val tailFunction: Kernel = when {
-    private val resB = meta.getDouble("B", 0.0)
+//        meta["tail"] != null -> {
-    private val tailFunction: BivariateFunction = when {
+//            val tailFunctionStr = meta["tail"].string
-        meta.hasValue("tail") -> {
+//            if (tailFunctionStr.startsWith("function::")) {
-            val tailFunctionStr = meta.getString("tail")
+//                FunctionLibrary.buildFrom(context).buildBivariateFunction(tailFunctionStr.substring(10))
-            if (tailFunctionStr.startsWith("function::")) {
+//            } else {
-                FunctionLibrary.buildFrom(context).buildBivariateFunction(tailFunctionStr.substring(10))
+//                BivariateFunction { E, U ->
-            } else {
+//                    val binding = HashMap<String, Any>()
-                BivariateFunction { E, U ->
+//                    binding["E"] = E
-                    val binding = HashMap<String, Any>()
+//                    binding["U"] = U
-                    binding["E"] = E
+//                    binding["D"] = E - U
-                    binding["U"] = U
+//                    ExpressionUtils.function(tailFunctionStr, binding)
-                    binding["D"] = E - U
+//                }
-                    ExpressionUtils.function(tailFunctionStr, binding)
+//            }
-                }
+//        }
-            }
+//        meta.hasValue("tailAlpha") -> {
-        }
+//            //add polynomial function here
-        meta.hasValue("tailAlpha") -> {
+//            val alpha = meta.getDouble("tailAlpha")
-            //add polynomial function here
+//            val beta = meta.getDouble("tailBeta", 0.0)
-            val alpha = meta.getDouble("tailAlpha")
+//            BivariateFunction { E: Double, U: Double -> 1 - (E - U) * (alpha + E / 1000.0 * beta) / 1000.0 }
-            val beta = meta.getDouble("tailBeta", 0.0)
+//
-            BivariateFunction { E: Double, U: Double -> 1 - (E - U) * (alpha + E / 1000.0 * beta) / 1000.0 }
+//        }
-
+//        else -> ResolutionFunction.getConstantTail()
-        }
+//    }
-        else -> ResolutionFunction.getConstantTail()
-    }
-
-    override fun derivValue(parName: String, x: Double, y: Double, set: Values): Double {
-        return 0.0
-    }
 
     private fun getValueFast(E: Double, U: Double): Double {
         val delta = resA * E
         return when {
             E - U < 0 -> 0.0
-            E - U > delta -> tailFunction.value(E, U)
+            E - U > delta -> tailFunction(E, U)
             else -> (E - U) / delta
         }
     }
 
-    override fun providesDeriv(name: String): Boolean {
+    override fun derivativeOrNull(symbols: List<Symbol>): Kernel = Kernel { _, _, _ -> 0.0 }
-        return true
-    }
 
-    override fun value(E: Double, U: Double, set: Values): Double {
+    override val x: Symbol get() = e
-        assert(resA > 0)
+    override val y: Symbol get() = u
+
+    override fun invoke(E: Double, U: Double, arguments: Map<Symbol, Double>): Double {
         if (resB <= 0) {
             return this.getValueFast(E, U)
         }
-        assert(resB > 0)
         val delta = resA * E
         return when {
             E - U < 0 -> 0.0
-            E - U > delta -> tailFunction.value(E, U)
+            E - U > delta -> tailFunction(E, U)
             else -> (1 - sqrt(1 - (E - U) / E * resB)) / (1 - sqrt(1 - resA * resB))
         }
     }
 
-    companion object {
 
-        private val list = arrayOf<String>() //leaving
+    public companion object {
+        public val e: Symbol by symbol
+        public val u: Symbol by symbol
+
     }
 
 }

numass-model/src/commonMain/kotlin/ru/inr/mass/models/NumassTransmission.kt

@@ -5,62 +5,63 @@
  */
 package inr.numass.models.sterile
 
-import hep.dataforge.context.Context
+import ru.inr.mass.models.DifferentiableKernel
-import hep.dataforge.maths.functions.FunctionLibrary
+import ru.inr.mass.models.Kernel
-import hep.dataforge.meta.Meta
+import ru.inr.mass.models.UnivariateFunction
-import hep.dataforge.stat.parametric.AbstractParametricBiFunction
+import space.kscience.kmath.integration.GaussIntegrator
-import hep.dataforge.values.Values
+import space.kscience.kmath.integration.integrate
-import inr.numass.models.misc.LossCalculator
+import space.kscience.kmath.misc.Symbol
-import inr.numass.utils.ExpressionUtils
+import space.kscience.kmath.misc.symbol
-import org.apache.commons.math3.analysis.BivariateFunction
+import space.kscience.kmath.operations.DoubleField
-import org.slf4j.LoggerFactory
+import kotlin.jvm.Synchronized
-import java.util.*
+import kotlin.math.*
-
 
 
 /**
  * @author [Alexander Nozik](mailto:altavir@gmail.com)
  */
-class NumassTransmission(context: Context, meta: Meta) : AbstractParametricBiFunction(list) {
+public class NumassTransmission(public val trapFunc: Kernel, private val adjustX: Boolean = false) :
-    private val trapFunc: BivariateFunction
+    DifferentiableKernel {
+    //    private val trapFunc: Kernel =         if (meta.hasValue("trapping")) {
+//        val trapFuncStr = meta.getString("trapping")
+//        trapFunc = if (trapFuncStr.startsWith("function::")) {
+//            FunctionLibrary.buildFrom(context).buildBivariateFunction(trapFuncStr.substring(10))
+//        } else {
+//            BivariateFunction { Ei: Double, Ef: Double ->
+//                val binding = HashMap<String, Any>()
+//                binding["Ei"] = Ei
+//                binding["Ef"] = Ef
+//                ExpressionUtils.function(trapFuncStr, binding)
+//            }
+//        }
+//    } else {
+//        LoggerFactory.getLogger(javaClass).warn("Trapping function not defined. Using default")
+//        trapFunc = FunctionLibrary.buildFrom(context).buildBivariateFunction("numass.trap.nominal")
+//    }
     //private val lossCache = HashMap<Double, UnivariateFunction>()
+    override fun derivativeOrNull(symbols: List<Symbol>): Kernel? = when (symbols.size) {
+        0 -> this
+        1 -> when (symbols.first()) {
+            trap -> trapFunc
+            thickness -> Kernel { eIn, eOut, arguments ->
+                val thickness = arguments[thickness] ?: 0.0
+                val probs = getLossProbDerivs(thickness)
 
-    private val adjustX: Boolean = meta.getBoolean("adjustX", false)
+                var sum = 0.0
-
+                for (i in 1 until probs.size) {
-    init {
+                    sum += probs[i] * getLossValue(i, eIn, eOut)
-        if (meta.hasValue("trapping")) {
-            val trapFuncStr = meta.getString("trapping")
-            trapFunc = if (trapFuncStr.startsWith("function::")) {
-                FunctionLibrary.buildFrom(context).buildBivariateFunction(trapFuncStr.substring(10))
-            } else {
-                BivariateFunction { Ei: Double, Ef: Double ->
-                    val binding = HashMap<String, Any>()
-                    binding["Ei"] = Ei
-                    binding["Ef"] = Ef
-                    ExpressionUtils.function(trapFuncStr, binding)
                 }
+                sum
             }
-        } else {
+            else -> null
-            LoggerFactory.getLogger(javaClass).warn("Trapping function not defined. Using default")
-            trapFunc = FunctionLibrary.buildFrom(context).buildBivariateFunction("numass.trap.nominal")
         }
+        else -> null
     }
 
-    override fun derivValue(parName: String, eIn: Double, eOut: Double, set: Values): Double {
+    override fun invoke(x: Double, y: Double, arguments: Map<Symbol, Double>): Double {
-        return when (parName) {
-            "trap" -> trapFunc.value(eIn, eOut)
-            "X" -> LossCalculator.getTotalLossDeriv(set, eIn, eOut)
-            else -> super.derivValue(parName, eIn, eOut, set)
-        }
-    }
-
-    override fun providesDeriv(name: String): Boolean {
-        return true
-    }
-
-    override fun value(eIn: Double, eOut: Double, set: Values): Double {
         // loss part
-        val loss = LossCalculator.getTotalLossValue(set, eIn, eOut)
+        val thickness = arguments[thickness] ?: 0.0
+        val loss = getTotalLossValue(thickness, x, y)
         //        double loss;
         //
         //        if(eIn-eOut >= 300){
@@ -74,13 +75,387 @@ class NumassTransmission(context: Context, meta: Meta) : AbstractParametricBiFun
         //        }
 
         //trapping part
-        val trap = getParameter("trap", set) * trapFunc.value(eIn, eOut)
+        val trap = arguments.getOrElse(trap) { 1.0 } * trapFunc(x, y, arguments)
         return loss + trap
     }
 
-    companion object {
+    public companion object {
+        public val trap: Symbol by symbol
+        public val thickness: Symbol by symbol
 
-        private val list = arrayOf("trap", "X")
+        private val cache = HashMap<Int, UnivariateFunction>()
+
+        private const val ION_POTENTIAL = 15.4//eV
+
+
+        private fun getX(arguments: Map<Symbol, Double>, eIn: Double, adjustX: Boolean = false): Double {
+            return if (adjustX) {
+                //From our article
+                arguments.getValue(thickness) * ln(eIn / ION_POTENTIAL) * eIn * ION_POTENTIAL / 1.9580741410115568e6
+            } else {
+                arguments.getValue(thickness)
+            }
+        }
+
+        private fun p0(eIn: Double, set: Map<Symbol, Double>): Double = getLossProbability(0, getX(set, eIn))
+
+        private fun getGunLossProbabilities(X: Double): List<Double> {
+            val res = ArrayList<Double>()
+            var prob: Double
+            if (X > 0) {
+                prob = exp(-X)
+            } else {
+                // если x ==0, то выживает только нулевой член, первый равен 1
+                res.add(1.0)
+                return res
+            }
+            res.add(prob)
+
+            var n = 0
+            while (prob > SCATTERING_PROBABILITY_THRESHOLD) {
+                /*
+                * prob(n) = prob(n-1)*X/n;
+                 */
+                n++
+                prob *= X / n
+                res.add(prob)
+            }
+
+            return res
+        }
+
+        fun getGunZeroLossProb(x: Double): Double {
+            return exp(-x)
+        }
+
+        private fun getCachedSpectrum(order: Int): UnivariateFunction {
+            return when {
+                order <= 0 -> error("Non-positive loss cache order")
+                order == 1 -> singleScatterFunction
+                else -> cache.getOrPut(order) {
+                    //LoggerFactory.getLogger(javaClass).debug("Scatter cache of order {} not found. Updating", order)
+                    getNextLoss(getMargin(order), getCachedSpectrum(order - 1))
+                }
+            }
+        }
+
+        /**
+         * Ленивое рекурсивное вычисление функции потерь через предыдущие
+         *
+         * @param order
+         * @return
+         */
+        private fun getLoss(order: Int): UnivariateFunction {
+            return getCachedSpectrum(order)
+        }
+
+        private fun getLossProbDerivs(x: Double): List<Double> {
+            val res = ArrayList<Double>()
+            val probs = getLossProbabilities(x)
+
+            var delta = exp(-x)
+            res.add((delta - probs[0]) / x)
+            for (i in 1 until probs.size) {
+                delta *= x / i
+                res.add((delta - probs[i]) / x)
+            }
+
+            return res
+        }
+
+        /**
+         * рекурсивно вычисляем все вероятности, котрорые выше порога
+         *
+         *
+         * дисер, стр.48
+         *
+         * @param X
+         * @return
+         */
+        private fun calculateLossProbabilities(x: Double): List<Double> {
+            val res = ArrayList<Double>()
+            var prob: Double
+            if (x > 0) {
+                prob = 1 / x * (1 - exp(-x))
+            } else {
+                // если x ==0, то выживает только нулевой член, первый равен нулю
+                res.add(1.0)
+                return res
+            }
+            res.add(prob)
+
+            while (prob > SCATTERING_PROBABILITY_THRESHOLD) {
+                /*
+            * prob(n) = prob(n-1)-1/n! * X^n * exp(-X);
+             */
+                var delta = exp(-x)
+                for (i in 1 until res.size + 1) {
+                    delta *= x / i
+                }
+                prob -= delta / x
+                res.add(prob)
+            }
+
+            return res
+        }
+
+        fun getLossProbabilities(x: Double): List<Double> = lossProbCache.getOrPut(x) { calculateLossProbabilities(x) }
+
+        fun getLossProbability(order: Int, X: Double): Double {
+            if (order == 0) {
+                return if (X > 0) {
+                    1 / X * (1 - exp(-X))
+                } else {
+                    1.0
+                }
+            }
+            val probs = getLossProbabilities(X)
+            return if (order >= probs.size) {
+                0.0
+            } else {
+                probs[order]
+            }
+        }
+
+        fun getLossValue(order: Int, Ei: Double, Ef: Double): Double {
+            return when {
+                Ei - Ef < 5.0 -> 0.0
+                Ei - Ef >= getMargin(order) -> 0.0
+                else -> getLoss(order).invoke(Ei - Ef)
+            }
+        }
+
+        /**
+         * функция потерь с произвольными вероятностями рассеяния
+         *
+         * @param probs
+         * @param Ei
+         * @param Ef
+         * @return
+         */
+        fun getLossValue(probs: List<Double>, Ei: Double, Ef: Double): Double {
+            var sum = 0.0
+            for (i in 1 until probs.size) {
+                sum += probs[i] * getLossValue(i, Ei, Ef)
+            }
+            return sum
+        }
+
+        /**
+         * граница интегрирования
+         *
+         * @param order
+         * @return
+         */
+        private fun getMargin(order: Int): Double {
+            return 80 + order * 50.0
+        }
+
+        /**
+         * генерирует кэшированную функцию свертки loss со спектром однократных
+         * потерь
+         *
+         * @param loss
+         * @return
+         */
+        @Synchronized
+        private fun getNextLoss(margin: Double, loss: UnivariateFunction): UnivariateFunction {
+            val res = { x: Double ->
+                integrator.integrate(5.0..margin) { y ->
+                    loss(x - y) * singleScatterFunction(y)
+                }
+            }
+
+            return FunctionCaching.cacheUnivariateFunction(0.0, margin, 200, res)
+
+        }
+
+        /**
+         * Значение полной производной функции потерь с учетом всех неисчезающих
+         * порядков
+         *
+         * @param X
+         * @param eIn
+         * @param eOut
+         * @return
+         */
+        private fun getTotalLossDeriv(X: Double, eIn: Double, eOut: Double): Double {
+            val probs = getLossProbDerivs(X)
+
+            var sum = 0.0
+            for (i in 1 until probs.size) {
+                sum += probs[i] * getLossValue(i, eIn, eOut)
+            }
+            return sum
+        }
+
+        /**
+         * Значение полной функции потерь с учетом всех неисчезающих порядков
+         *
+         * @param x
+         * @param Ei
+         * @param Ef
+         * @return
+         */
+        fun getTotalLossValue(x: Double, Ei: Double, Ef: Double): Double {
+            return if (x == 0.0) {
+                0.0
+            } else {
+                val probs = getLossProbabilities(x)
+                (1 until probs.size).sumByDouble { i ->
+                    probs[i] * getLossValue(i, Ei, Ef)
+                }
+            }
+        }
+
+
+        /**
+         * порог по вероятности, до которого вычисляются компоненты функции потерь
+         */
+        private const val SCATTERING_PROBABILITY_THRESHOLD = 1e-3
+        private val integrator = GaussIntegrator(DoubleField)
+        private val lossProbCache = Misc.getLRUCache<Double, List<Double>>(100)
+
+
+        private val A1 = 0.204
+        private val A2 = 0.0556
+        private val b = 14.0
+        private val pos1 = 12.6
+        private val pos2 = 14.3
+        private val w1 = 1.85
+        private val w2 = 12.5
+
+        public val singleScatterFunction: UnivariateFunction = { eps: Double ->
+            when {
+                eps <= 0 -> 0.0
+                eps <= b -> {
+                    val z = eps - pos1
+                    A1 * exp(-2.0 * z * z / w1 / w1)
+                }
+                else -> {
+                    val z = 4.0 * (eps - pos2) * (eps - pos2)
+                    A2 / (1 + z / w2 / w2)
+                }
+            }
+        }
+
+
+        /**
+         * A generic loss function for numass experiment in "Lobashev"
+         * parameterization
+         *
+         * @param exPos
+         * @param ionPos
+         * @param exW
+         * @param ionW
+         * @param exIonRatio
+         * @return
+         */
+        public fun getSingleScatterFunction(
+            exPos: Double,
+            ionPos: Double,
+            exW: Double,
+            ionW: Double,
+            exIonRatio: Double,
+        ): UnivariateFunction {
+            val func: UnivariateFunction = { eps: Double ->
+                if (eps <= 0) {
+                    0.0
+                } else {
+                    val z1 = eps - exPos
+                    val ex = exIonRatio * exp(-2.0 * z1 * z1 / exW / exW)
+
+                    val z = 4.0 * (eps - ionPos) * (eps - ionPos)
+                    val ion = 1 / (1 + z / ionW / ionW)
+
+                    if (eps < exPos) {
+                        ex
+                    } else {
+                        max(ex, ion)
+                    }
+                }
+            }
+
+            val cutoff = 25.0
+            //caclulating lorentz integral analythically
+            val tailNorm = (atan((ionPos - cutoff) * 2.0 / ionW) + 0.5 * PI) * ionW / 2.0
+            val norm: Double = integrator.integrate(range = 0.0..cutoff, function = func) + tailNorm
+            return { e -> func(e) / norm }
+        }
+
+
+        public val exPos: Symbol by symbol
+        public val ionPos: Symbol by symbol
+        public val exW: Symbol by symbol
+        public val ionW: Symbol by symbol
+        public val exIonRatio: Symbol by symbol
+
+        public fun getSingleScatterFunction(set: Map<Symbol, Double>): UnivariateFunction {
+            val exPos = set.getValue(exPos)
+            val ionPos = set.getValue(ionPos)
+            val exW = set.getValue(exW)
+            val ionW = set.getValue(ionW)
+            val exIonRatio = set.getValue(exIonRatio)
+
+            return getSingleScatterFunction(exPos, ionPos, exW, ionW, exIonRatio)
+        }
+
+        public val trapFunction: (Double, Double) -> Double = { Ei: Double, Ef: Double ->
+            val eps = Ei - Ef
+            if (eps > 10) {
+                1.86e-04 * exp(-eps / 25.0) + 5.5e-05
+            } else {
+                0.0
+            }
+        }
+
+//        fun plotScatter(frame: PlotFrame, set: Values) {
+//            //"X", "shift", "exPos", "ionPos", "exW", "ionW", "exIonRatio"
+//
+//            //        JFreeChartFrame frame = JFreeChartFrame.drawFrame("Differential scattering crosssection", null);
+//            val X = set.getDouble("X")
+//
+//            val exPos = set.getDouble("exPos")
+//
+//            val ionPos = set.getDouble("ionPos")
+//
+//            val exW = set.getDouble("exW")
+//
+//            val ionW = set.getDouble("ionW")
+//
+//            val exIonRatio = set.getDouble("exIonRatio")
+//
+//            val scatterFunction = getSingleScatterFunction(exPos, ionPos, exW, ionW, exIonRatio)
+//
+//            if (set.names.contains("X")) {
+//                val probs = LossCalculator.getGunLossProbabilities(set.getDouble("X"))
+//                val single = { e: Double -> probs[1] * scatterFunction.value(e) }
+//                frame.add(XYFunctionPlot.plot("Single scattering", 0.0, 100.0, 1000) { x: Double -> single(x) })
+//
+//                for (i in 2 until probs.size) {
+//                    val scatter = { e: Double -> probs[i] * LossCalculator.getLossValue(i, e, 0.0) }
+//                    frame.add(XYFunctionPlot.plot(i.toString() + " scattering", 0.0, 100.0, 1000) { x: Double -> scatter(x) })
+//                }
+//
+//                val total = UnivariateFunction { eps ->
+//                    if (probs.size == 1) {
+//                        return@UnivariateFunction 0.0
+//                    }
+//                    var sum = probs[1] * scatterFunction.value(eps)
+//                    for (i in 2 until probs.size) {
+//                        sum += probs[i] * LossCalculator.getLossValue(i, eps, 0.0)
+//                    }
+//                    return@UnivariateFunction sum
+//                }
+//
+//                frame.add(XYFunctionPlot.plot("Total loss", 0.0, 100.0, 1000) { x: Double -> total.value(x) })
+//
+//            } else {
+//
+//                frame.add(XYFunctionPlot.plot("Differential cross-section", 0.0, 100.0, 2000) { x: Double -> scatterFunction.value(x) })
+//            }
+//
+//        }
     }
 
 }

numass-model/src/commonMain/kotlin/ru/inr/mass/models/ParametricBiFunction.kt

@@ -1,40 +0,0 @@
-/*
- * To change this license header, choose License Headers in Project Properties.
- * To change this template file, choose Tools | Templates
- * and open the template in the editor.
- */
-package hep.dataforge.stat.parametric
-
-import hep.dataforge.exceptions.NotDefinedException
-
-/**
- *
- * @author Alexander Nozik
- */
-interface ParametricBiFunction : NameSetContainer {
-    fun derivValue(parName: String?, x: Double, y: Double, set: Values?): Double
-    fun value(x: Double, y: Double, set: Values?): Double
-    fun providesDeriv(name: String?): Boolean
-    fun derivative(parName: String?): ParametricBiFunction? {
-        return if (providesDeriv(parName)) {
-            object : ParametricBiFunction {
-                override fun derivValue(parName: String?, x: Double, y: Double, set: Values?): Double {
-                    throw NotDefinedException()
-                }
-
-                override fun value(x: Double, y: Double, set: Values?): Double {
-                    return this@ParametricBiFunction.derivValue(parName, x, y, set)
-                }
-
-                override fun providesDeriv(name: String?): Boolean {
-                    return !names.contains(name)
-                }
-
-                val names: NameList
-                    get() = this@ParametricBiFunction.getNames()
-            }
-        } else {
-            throw NotDefinedException()
-        }
-    }
-}

numass-model/src/commonMain/kotlin/ru/inr/mass/models/ParametricBiFunctionCache.kt

@@ -1,39 +0,0 @@
-/*
- * Copyright  2018 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 inr.numass.models.sterile
-
-import hep.dataforge.names.NameList
-import hep.dataforge.stat.parametric.ParametricBiFunction
-import hep.dataforge.values.Values
-
-class ParametricBiFunctionCache(val function: ParametricBiFunction): ParametricBiFunction {
-    override fun derivValue(parName: String?, x: Double, y: Double, set: Values?): Double {
-        TODO("not implemented") //To change body of created functions use File | Settings | File Templates.
-    }
-
-    override fun getNames(): NameList {
-        TODO("not implemented") //To change body of created functions use File | Settings | File Templates.
-    }
-
-    override fun value(x: Double, y: Double, set: Values?): Double {
-        TODO("not implemented") //To change body of created functions use File | Settings | File Templates.
-    }
-
-    override fun providesDeriv(name: String?): Boolean {
-        TODO("not implemented") //To change body of created functions use File | Settings | File Templates.
-    }
-}

numass-model/src/commonMain/kotlin/ru/inr/mass/models/ParametricFunction.kt

@@ -1,58 +0,0 @@
-/* 
- * 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 hep.dataforge.stat.parametric
-
-import hep.dataforge.exceptions.NotDefinedException
-
-/**
- *
- *
- * NamedSpectrum interface.
- *
- * @author Alexander Nozik
- * @version $Id: $Id
- */
-interface ParametricFunction : NameSetContainer {
-    fun derivValue(parName: String?, x: Double, set: Values?): Double
-    fun value(x: Double, set: Values?): Double
-    fun providesDeriv(name: String?): Boolean
-    fun derivative(parName: String?): ParametricFunction? {
-        return if (providesDeriv(parName)) {
-            object : ParametricFunction {
-                override fun derivValue(parName: String?, x: Double, set: Values?): Double {
-                    return if (names.contains(parName)) {
-                        throw NotDefinedException()
-                    } else {
-                        0
-                    }
-                }
-
-                override fun value(x: Double, set: Values?): Double {
-                    return this@ParametricFunction.derivValue(parName, x, set)
-                }
-
-                override fun providesDeriv(name: String?): Boolean {
-                    return !names.contains(name)
-                }
-
-                val names: NameList
-                    get() = this@ParametricFunction.getNames()
-            }
-        } else {
-            throw NotDefinedException()
-        }
-    }
-}

numass-model/src/commonMain/kotlin/ru/inr/mass/models/ParametricValue.kt

@@ -1,52 +0,0 @@
-/* 
- * 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 hep.dataforge.stat.parametric
-
-import hep.dataforge.exceptions.NotDefinedException
-
-/**
- * A function mapping parameter set to real value
- *
- * @author Alexander Nozik
- */
-interface ParametricValue : NameSetContainer {
-    /**
-     * Value
-     * @param pars
-     * @return
-     */
-    fun value(pars: Values?): Double
-
-    /**
-     * Partial derivative value for given parameter
-     * @param derivParName
-     * @param pars
-     * @return
-     */
-    fun derivValue(derivParName: String?, pars: Values?): Double {
-        throw NotDefinedException()
-    }
-
-    /**
-     * Returns true if this object provides explicit analytical value derivative for given parameter
-     *
-     * @param name a [String] object.
-     * @return a boolean.
-     */
-    fun providesDeriv(name: String?): Boolean {
-        return false
-    }
-}

numass-model/src/commonMain/kotlin/ru/inr/mass/models/Spectrum.kt

@@ -0,0 +1,38 @@
+package ru.inr.mass.models
+
+import space.kscience.kmath.expressions.DifferentiableExpression
+import space.kscience.kmath.expressions.Expression
+import space.kscience.kmath.misc.Symbol
+
+public fun interface Spectrum : Expression<Double> {
+    public val abscissa: Symbol get() = Symbol.x
+
+    public operator fun invoke(x: Double, arguments: Map<Symbol, Double>): Double
+
+    override fun invoke(arguments: Map<Symbol, Double>): Double =
+        invoke(arguments[abscissa] ?: error("Argument $abscissa not found in arguments"), arguments)
+}
+
+public interface DifferentiableSpectrum : DifferentiableExpression<Double, Spectrum>, Spectrum
+
+public fun interface Kernel : Expression<Double> {
+    public val x: Symbol get() = Symbol.x
+    public val y: Symbol get() = Symbol.y
+
+    public operator fun invoke(x: Double, y: Double, arguments: Map<Symbol, Double>): Double
+
+    override fun invoke(arguments: Map<Symbol, Double>): Double {
+        val xValue = arguments[x] ?: error("$x value not found in arguments")
+        val yValue = arguments[y] ?: error("$y value not found in arguments")
+        return invoke(xValue, yValue, arguments)
+    }
+}
+
+public interface DifferentiableKernel : DifferentiableExpression<Double, Kernel>, Kernel
+
+public fun <T> Expression<T>.withDefault(default: Map<Symbol, T>): Expression<T> = Expression { args ->
+    invoke(default + args)
+}
+
+public typealias UnivariateFunction = (Double) -> Double
+public typealias BivariateFunction = (Double, Double) -> Double

numass-model/src/jvmMain/kotlin/ru/inr/mass/maths/DerivativeCalculator.kt

@@ -13,9 +13,11 @@
  * See the License for the specific language governing permissions and
  * limitations under the License.
  */
-package hep.dataforge.stat.parametric
+package ru.inr.mass.maths
 
 import hep.dataforge.stat.fit.Param
+import hep.dataforge.stat.parametric.ParametricUtils
+import hep.dataforge.stat.parametric.ParametricValue
 
 /**
  *