Update spectrum computation

2021-05-23 20:01:36 +03:00 · 2021-05-23 20:01:36 +03:00 · 5e1a4c9477
commit 5e1a4c9477
parent 0a5812f0fa
14 changed files with 187 additions and 113 deletions
--- a/build.gradle.kts
+++ b/build.gradle.kts
@ -9,17 +9,13 @@ allprojects {
    }

    group = "ru.inr.mass"
-    version = "0.1.0-SNAPSHOT"
+    version = "0.1.0-dev-1"
 }

 val dataforgeVersion by extra("0.4.0")
-val kmathVersion by extra("0.3.0-dev-10")
-
-apiValidation{
-    validationDisabled = true
-}
+val kmathVersion by extra("0.3.0-dev-13")

 ksciencePublish{
-    configurePublications("https://mipt-npm.jetbrains.space/p/numass/code/numass/")
+    vcs("https://mipt-npm.jetbrains.space/p/numass/code/numass/")
    space("https://maven.pkg.jetbrains.space/mipt-npm/p/numass/maven")
 }
--- a/numass-model/src/commonMain/kotlin/ru/inr/mass/models/FSS.kt
+++ b/numass-model/src/commonMain/kotlin/ru/inr/mass/models/FSS.kt
@ -16,13 +16,10 @@
 package ru.inr.mass.models

 import space.kscience.kmath.data.ColumnarData
-import space.kscience.kmath.misc.Symbol
+import space.kscience.kmath.expressions.Symbol
+import space.kscience.kmath.expressions.symbol
 import space.kscience.kmath.misc.UnstableKMathAPI
-import space.kscience.kmath.misc.symbol
-import space.kscience.kmath.real.div
-import space.kscience.kmath.real.sum
 import space.kscience.kmath.structures.Buffer
-import space.kscience.kmath.structures.DoubleBuffer


@OptIn(UnstableKMathAPI::class)
@ -39,18 +36,5 @@ public class FSS(public val es: Buffer<Double>, public val ps: Buffer<Double>) :
    public companion object {
        public val p: Symbol by symbol
        public val e: Symbol by symbol
-
-        public val default: FSS by lazy {
-            val stream = FSS::class.java.getResourceAsStream("/data/FS.txt") ?: error("Default FS resource not found")
-            stream.use { inputStream ->
-                val data = inputStream.bufferedReader().lineSequence().map {
-                    val (e, p) = it.split("\t")
-                    e.toDouble() to p.toDouble()
-                }.toList()
-                val es = DoubleBuffer(data.size) { data[it].first }
-                val ps = DoubleBuffer(data.size) { data[it].second }
-                FSS(es, ps / ps.sum())
-            }
-        }
    }
 }
--- a/numass-model/src/commonMain/kotlin/ru/inr/mass/models/NBkgSpectrum.kt
+++ b/numass-model/src/commonMain/kotlin/ru/inr/mass/models/NBkgSpectrum.kt
@ -15,8 +15,9 @@
 */
 package ru.inr.mass.models

-import space.kscience.kmath.misc.Symbol
-import space.kscience.kmath.misc.symbol
+import space.kscience.kmath.expressions.Symbol
+import space.kscience.kmath.expressions.symbol
+

 /**
 *
--- a/numass-model/src/commonMain/kotlin/ru/inr/mass/models/NumassBeta.kt
+++ b/numass-model/src/commonMain/kotlin/ru/inr/mass/models/NumassBeta.kt
@ -3,13 +3,11 @@
 * To change this template file, choose Tools | Templates
 * and open the template in the editor.
 */
-package inr.numass.models.sterile
+package ru.inr.mass.models

-import ru.inr.mass.models.DifferentiableKernel
-import ru.inr.mass.models.Kernel
-import space.kscience.kmath.misc.StringSymbol
-import space.kscience.kmath.misc.Symbol
-import space.kscience.kmath.misc.symbol
+import space.kscience.kmath.expressions.StringSymbol
+import space.kscience.kmath.expressions.Symbol
+import space.kscience.kmath.expressions.symbol
 import kotlin.math.*

 /**
--- a/numass-model/src/commonMain/kotlin/ru/inr/mass/models/NumassResolution.kt
+++ b/numass-model/src/commonMain/kotlin/ru/inr/mass/models/NumassResolution.kt
@ -3,12 +3,10 @@
 * To change this template file, choose Tools | Templates
 * and open the template in the editor.
 */
-package inr.numass.models.sterile
+package ru.inr.mass.models

-import ru.inr.mass.models.DifferentiableKernel
-import ru.inr.mass.models.Kernel
-import space.kscience.kmath.misc.Symbol
-import space.kscience.kmath.misc.symbol
+import space.kscience.kmath.expressions.Symbol
+import space.kscience.kmath.expressions.symbol
 import kotlin.math.sqrt

 /**
--- a/numass-model/src/commonMain/kotlin/ru/inr/mass/models/NumassTransmission.kt
+++ b/numass-model/src/commonMain/kotlin/ru/inr/mass/models/NumassTransmission.kt
@ -3,27 +3,21 @@
 * To change this template file, choose Tools | Templates
 * and open the template in the editor.
 */
-package inr.numass.models.sterile
+package ru.inr.mass.models

-import ru.inr.mass.models.DifferentiableKernel
-import ru.inr.mass.models.Kernel
-import ru.inr.mass.models.cache
+import space.kscience.kmath.expressions.Symbol
+import space.kscience.kmath.expressions.symbol
 import space.kscience.kmath.functions.PiecewisePolynomial
 import space.kscience.kmath.functions.UnivariateFunction
-import space.kscience.kmath.functions.value
+import space.kscience.kmath.functions.asFunction
 import space.kscience.kmath.integration.GaussIntegrator
 import space.kscience.kmath.integration.integrate
 import space.kscience.kmath.integration.value
-import space.kscience.kmath.misc.Symbol
-import space.kscience.kmath.misc.symbol
 import space.kscience.kmath.operations.DoubleField
+import kotlin.jvm.Synchronized
 import kotlin.math.*


-public fun PiecewisePolynomial<Double>.asFunction(defaultValue: Double = 0.0): UnivariateFunction<Double> = {
-    value(DoubleField, it) ?: defaultValue
-}
-
 /**
 * @author [Alexander Nozik](mailto:altavir@gmail.com)
 */
@ -144,7 +138,7 @@ public class NumassTransmission(
                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)).asFunction()
+                    getNextLoss(getMargin(order), getCachedSpectrum(order - 1)).asFunction(DoubleField, 0.0)
                }
            }
        }
--- a/numass-model/src/commonMain/kotlin/ru/inr/mass/models/Spectrum.kt
+++ b/numass-model/src/commonMain/kotlin/ru/inr/mass/models/Spectrum.kt
@ -2,7 +2,7 @@ package ru.inr.mass.models

 import space.kscience.kmath.expressions.DifferentiableExpression
 import space.kscience.kmath.expressions.Expression
-import space.kscience.kmath.misc.Symbol
+import space.kscience.kmath.expressions.Symbol

 public fun interface Spectrum : Expression<Double> {
    public val abscissa: Symbol get() = Symbol.x
--- a/numass-model/src/commonMain/kotlin/ru/inr/mass/models/SterileNeutrinoSpectrum.kt
+++ b/numass-model/src/commonMain/kotlin/ru/inr/mass/models/SterileNeutrinoSpectrum.kt
@ -3,23 +3,25 @@
 * To change this template file, choose Tools | Templates
 * and open the template in the editor.
 */
-package inr.numass.models.sterile
+package ru.inr.mass.models


-import inr.numass.models.sterile.NumassBeta.e0
-import inr.numass.models.sterile.NumassBeta.mnu2
-import inr.numass.models.sterile.NumassBeta.msterile2
-import inr.numass.models.sterile.NumassBeta.u2
-import inr.numass.models.sterile.NumassTransmission.Companion.thickness
-import inr.numass.models.sterile.NumassTransmission.Companion.trap
-import ru.inr.mass.models.*
+import ru.inr.mass.models.NumassBeta.e0
+import ru.inr.mass.models.NumassBeta.mnu2
+import ru.inr.mass.models.NumassBeta.msterile2
+import ru.inr.mass.models.NumassBeta.u2
+import ru.inr.mass.models.NumassTransmission.Companion.thickness
+import ru.inr.mass.models.NumassTransmission.Companion.trap
+import space.kscience.kmath.expressions.Symbol
 import space.kscience.kmath.expressions.derivative
+import space.kscience.kmath.integration.UnivariateIntegrandRanges
+import space.kscience.kmath.integration.gaussIntegrator
 import space.kscience.kmath.integration.integrate
-import space.kscience.kmath.integration.integrator
 import space.kscience.kmath.integration.value
-import space.kscience.kmath.misc.Symbol
 import space.kscience.kmath.operations.DoubleField
-import kotlin.math.min
+import space.kscience.kmath.real.step
+import space.kscience.kmath.structures.toDoubleArray
+

 /**
 * @param source variables:Eo offset,Ein; parameters: "mnu2", "msterile2", "U2"
@ -30,7 +32,7 @@ public class SterileNeutrinoSpectrum(
    public val source: DifferentiableKernel = NumassBeta,
    public val transmission: DifferentiableKernel = NumassTransmission(),
    public val resolution: DifferentiableKernel = NumassResolution(),
-    public val fss: FSS? = FSS.default,
+    public val fss: FSS? = null,
 ) : DifferentiableSpectrum {

    /**
@ -91,9 +93,16 @@ public class SterileNeutrinoSpectrum(
 //            getHighDensityIntegrator()
 //        }

-        return DoubleField.integrator.integrate(u..eMax) { eIn ->
+        return DoubleField.gaussIntegrator.integrate(u..eMax, generateRanges(
+            u..eMax,
+            u + 2.0,
+            u + 7.0,
+            u + 15.0,
+            u + 30.0,
+            *((u + 50)..(u + 6000) step 25.0).toDoubleArray()
+        )) { eIn ->
            sumByFSS(eIn, sourceFunction, arguments) * transResFunction(eIn, u, arguments)
-        }.value ?: error("Integration failed")
+        }.value
    }

    private fun sumByFSS(eIn: Double, sourceFunction: Kernel, arguments: Map<Symbol, Double>): Double {
@ -121,20 +130,35 @@ public class SterileNeutrinoSpectrum(
            }
        }

-        private fun lossRes(transFunc: Kernel, eIn: Double, u: Double, arguments: Map<Symbol, Double>): Double {
-            val integrand = { eOut: Double -> transFunc(eIn, eOut, arguments) * resolution(eOut, u, arguments) }
-
-            val border = u + 30
-            val firstPart = DoubleField.integrator.integrate(u..min(eIn, border), function = integrand).value
-            val secondPart: Double = if (eIn > border) {
-                DoubleField.integrator.integrate(border..eIn, function = integrand).value
-            } else {
-                0.0
-            }
-            return firstPart + secondPart
-        }
-
+        private fun lossRes(
+            transFunc: Kernel,
+            eIn: Double,
+            u: Double,
+            arguments: Map<Symbol, Double>,
+        ): Double = DoubleField.gaussIntegrator.integrate(u..eIn, generateRanges(
+            u..eIn,
+            *((u + 25)..(u + 6000) step 25.0).toDoubleArray()
+        )) { eOut: Double ->
+            transFunc(eIn, eOut, arguments) * resolution(eOut, u, arguments)
+        }.value
    }

    public companion object
 }
+
+
+internal fun generateRanges(
+    range: ClosedFloatingPointRange<Double>,
+    vararg borders: Double,
+    points: Int = 5,
+): UnivariateIntegrandRanges {
+    if (borders.isEmpty() || borders.first() > range.endInclusive) return UnivariateIntegrandRanges(range to points)
+    val ranges = listOf(
+        range.start,
+        *borders.filter { it in range }.sorted().toTypedArray(),
+        range.endInclusive
+    ).zipWithNext { l, r ->
+        l..r to points
+    }
+    return UnivariateIntegrandRanges(ranges)
+}
--- a/numass-model/src/commonMain/kotlin/ru/inr/mass/models/functionCaching.kt
+++ b/numass-model/src/commonMain/kotlin/ru/inr/mass/models/functionCaching.kt
--- a/numass-model/src/commonTest/kotlin/ru/inr/mass/models/TestGenerateRanges.kt
+++ b/numass-model/src/commonTest/kotlin/ru/inr/mass/models/TestGenerateRanges.kt
@ -0,0 +1,14 @@
+package ru.inr.mass.models
+
+import kotlin.test.Test
+import kotlin.test.assertEquals
+
+internal class TestGenerateRanges {
+    @Test
+    fun simpleRanges() {
+        val ranges = generateRanges(0.0..100.0, 10.0, 55.0, 120.0)
+        assertEquals(3, ranges.ranges.size)
+        assertEquals(55.0..100.0, ranges.ranges.last().first)
+        assertEquals(10.0..55.0, ranges.ranges[1].first)
+    }
+}
--- a/numass-model/src/jvmMain/kotlin/ru/inr/mass/models/defaultFss.kt
+++ b/numass-model/src/jvmMain/kotlin/ru/inr/mass/models/defaultFss.kt
@ -0,0 +1,20 @@
+package ru.inr.mass.models
+
+import space.kscience.kmath.real.div
+import space.kscience.kmath.real.sum
+import space.kscience.kmath.structures.DoubleBuffer
+
+private val defaultFss: FSS by lazy {
+    val stream = FSS::class.java.getResourceAsStream("/data/FS.txt") ?: error("Default FS resource not found")
+    stream.use { inputStream ->
+        val data = inputStream.bufferedReader().lineSequence().map {
+            val (e, p) = it.split("\t")
+            e.toDouble() to p.toDouble()
+        }.toList()
+        val es = DoubleBuffer(data.size) { data[it].first }
+        val ps = DoubleBuffer(data.size) { data[it].second }
+        FSS(es, ps / ps.sum())
+    }
+}
+
+public val FSS.Companion.default: FSS get() = defaultFss
--- a/numass-workspace/src/main/kotlin/ru/inr/mass/scripts/sterileSpectrum.kt
+++ b/numass-workspace/src/main/kotlin/ru/inr/mass/scripts/sterileSpectrum.kt
@ -1,28 +1,27 @@
 package ru.inr.mass.scripts

-import inr.numass.models.sterile.NumassBeta.e0
-import inr.numass.models.sterile.NumassBeta.mnu2
-import inr.numass.models.sterile.NumassBeta.msterile2
-import inr.numass.models.sterile.NumassBeta.u2
-import inr.numass.models.sterile.NumassResolution
-import inr.numass.models.sterile.NumassTransmission
-import inr.numass.models.sterile.NumassTransmission.Companion.thickness
-import inr.numass.models.sterile.NumassTransmission.Companion.trap
-import inr.numass.models.sterile.SterileNeutrinoSpectrum
 import kotlinx.html.code
+import ru.inr.mass.models.*
 import ru.inr.mass.models.NBkgSpectrum.Companion.bkg
 import ru.inr.mass.models.NBkgSpectrum.Companion.norm
-import ru.inr.mass.models.withNBkg
+import ru.inr.mass.models.NumassBeta.e0
+import ru.inr.mass.models.NumassBeta.mnu2
+import ru.inr.mass.models.NumassBeta.msterile2
+import ru.inr.mass.models.NumassBeta.u2
+import ru.inr.mass.models.NumassTransmission.Companion.thickness
+import ru.inr.mass.models.NumassTransmission.Companion.trap
 import ru.inr.mass.workspace.buffer
-import space.kscience.kmath.misc.Symbol
+import space.kscience.kmath.expressions.Symbol
 import space.kscience.kmath.real.step
 import space.kscience.plotly.*
+import space.kscience.plotly.models.AxisType
 import space.kscience.plotly.models.ScatterMode
 import space.kscience.plotly.models.appendXY
 import kotlin.math.pow
+import kotlin.system.measureTimeMillis

 fun main() {
-    val spectrum = SterileNeutrinoSpectrum().withNBkg()
+    val spectrum = SterileNeutrinoSpectrum(fss = FSS.default).withNBkg()

    val args: Map<Symbol, Double> = mapOf(
        norm to 8e5,
@ -36,23 +35,32 @@ fun main() {
    )

    Plotly.page {
+        val spectrumTime = measureTimeMillis {
+            plot {
+                scatter {
+                    name = "Computed spectrum"
+                    mode = ScatterMode.lines
+                    x.buffer = 14000.0..18600.0 step 10.0
+                    y.numbers = x.doubles.map { spectrum(it, args) }
+                }
+                layout {
+                    title = "Sterile neutrino spectrum"
+                    yaxis.type = AxisType.log
+                }
+            }
+        }
+        println("Spectrum with 460 points computed in $spectrumTime millis")
+
        plot {
            scatter {
-                name = "Computed spectrum"
-                mode = ScatterMode.lines
-                x.buffer = 14000.0..18600.0 step 10.0
-                y.numbers = x.doubles.map { spectrum(it, args) }
-            }
-            scatter {
-                name = "Old spectrum"
                mode = ScatterMode.markers
                javaClass.getResource("/old-spectrum.dat").readText().lines().map {
-                    val (u, w) = it.split("\t")
-                    appendXY(u.toDouble(), w.toDouble())
+                    val (u, w) = it.split("\t").map { it.toDouble() }
+                    appendXY(u, w / spectrum(u, args) - 1.0)
                }
            }
            layout {
-                title = "Sterile neutrino spectrum"
+                title = "Sterile neutrino old/new ratio"
            }
        }
        plot {
@ -79,19 +87,8 @@ fun main() {
        }

        code {
-            +"""
-            norm to 8e5,
-            bkg to 2.0,
-            mnu2 to 0.0,
-            e0 to 18575.0,
-            msterile2 to 1000.0.pow(2),
-            u2 to 1e-2,
-            thickness to 0.1,
-            trap to 1.0 
-            """.trimIndent()
+            +args.toString()
        }
    }.makeFile()

-    println(spectrum(14000.0, args))
-
 }
--- a/numass-workspace/src/test/kotlin/ru/inr/mass/TestSterileSpectrumShape.kt
+++ b/numass-workspace/src/test/kotlin/ru/inr/mass/TestSterileSpectrumShape.kt
@ -0,0 +1,47 @@
+package ru.inr.mass
+
+import ru.inr.mass.models.*
+import space.kscience.kmath.misc.Symbol
+import kotlin.math.abs
+import kotlin.math.pow
+import kotlin.math.sqrt
+import kotlin.test.Ignore
+import kotlin.test.Test
+import kotlin.test.assertTrue
+
+internal class TestSterileSpectrumShape {
+    val args: Map<Symbol, Double> = mapOf(
+        NBkgSpectrum.norm to 8e5,
+        NBkgSpectrum.bkg to 2.0,
+        NumassBeta.mnu2 to 0.0,
+        NumassBeta.e0 to 18575.0,
+        NumassBeta.msterile2 to 1000.0.pow(2),
+        NumassBeta.u2 to 1e-2,
+        NumassTransmission.thickness to 0.1,
+        NumassTransmission.trap to 1.0
+    )
+
+    val spectrum = SterileNeutrinoSpectrum(fss = FSS.default).withNBkg()
+
+    val data by lazy {
+        javaClass.getResource("/old-spectrum.dat").readText().lines().map {
+            val (u, oldValue) = it.split("\t")
+            u.toDouble() to oldValue.toDouble()
+        }
+    }
+
+    @Test
+    @Ignore
+    fun checkAbsoluteDif() {
+        val t = 1e6
+
+        val res = data.sumOf { (u, oldValue) ->
+            val newValue = spectrum(u, args)
+            abs(newValue - oldValue) * sqrt(t) / sqrt(oldValue)
+        }
+
+        println(res / data.size)
+
+        assertTrue { res / data.size < 0.1 }
+    }
+}
--- a/settings.gradle.kts
+++ b/settings.gradle.kts
@ -4,9 +4,10 @@ pluginManagement {
        mavenCentral()
        jcenter()
        gradlePluginPortal()
+        mavenLocal()
    }

-    val toolsVersion = "0.9.7"
+    val toolsVersion = "0.9.9"
    val kotlinVersion = "1.5.0"

    plugins {