refactor: optimise SterileNeutrinoSpectrumRWall

numass-workspace/src/main/kotlin/ru/inr/mass/scripts/SterileNeutrinoSpectrumRWall.kt

@@ -0,0 +1,80 @@
+package ru.inr.mass.scripts
+
+import ru.inr.mass.models.DifferentiableKernel
+import ru.inr.mass.models.DifferentiableSpectrum
+import ru.inr.mass.models.NBkgSpectrum
+import ru.inr.mass.models.NumassBeta
+import ru.inr.mass.models.NumassResolution
+import ru.inr.mass.models.NumassTransmission
+import ru.inr.mass.models.Spectrum
+import ru.inr.mass.models.SterileNeutrinoSpectrum
+import ru.inr.mass.models.SterileNeutrinoSpectrumFast
+import ru.inr.mass.workspace.fitWith
+import ru.inr.mass.workspace.freeParameters
+import ru.inr.mass.workspace.iterations
+import space.kscience.kmath.UnstableKMathAPI
+import space.kscience.kmath.data.XYErrorColumnarData
+import space.kscience.kmath.expressions.Symbol
+import space.kscience.kmath.expressions.symbol
+import space.kscience.kmath.misc.Loggable
+import space.kscience.kmath.optimization.OptimizationLog
+import space.kscience.kmath.optimization.QowOptimizer
+import space.kscience.kmath.optimization.XYFit
+
+class SterileNeutrinoSpectrumRWall(
+    private val wallModel: WallModel,
+    private val trapModel: ITrap,
+    private val adiabaticityModel: AdiabaticityModel,
+    source: DifferentiableKernel = NumassBeta,
+    parallel: Boolean = true,
+) : DifferentiableSpectrum {
+
+    val spectrum = when (parallel) {
+        true -> SterileNeutrinoSpectrumFast(
+            source = source, transmission = NumassTransmission(
+                trapFunc = { ei, ef, _ -> trapModel.value(ei, ei - ef) },
+                adjustX = false,
+            ), resolution = NumassResolution(1.7e-4, tailFunction = { e, u -> adiabaticityModel(e, u) }), fss = null
+        )
+
+        false -> SterileNeutrinoSpectrum(
+            source = source, transmission = NumassTransmission(
+                trapFunc = { ei, ef, _ -> trapModel.value(ei, ei - ef) },
+                adjustX = false,
+            ), resolution = NumassResolution(1.7e-4, tailFunction = { e, u -> adiabaticityModel(e, u) }), fss = null
+        )
+    }
+
+    override fun invoke(u: Double, arguments: Map<Symbol, Double>): Double {
+        val rw = arguments[rearWall]!!
+        return spectrum.invoke(u, arguments) * (1.0 + rw * wallModel(u))
+    }
+
+    override fun derivativeOrNull(symbols: List<Symbol>): Spectrum? {
+        if (symbols.isEmpty()) return this
+        return when (symbols.singleOrNull() ?: TODO("First derivatives only")) {
+            rearWall -> Spectrum { u, args -> wallModel(u) * spectrum.invoke(u, args) }
+            else -> spectrum.derivativeOrNull(symbols)
+        }
+    }
+
+    companion object {
+        val rearWall: Symbol by symbol
+
+        @OptIn(UnstableKMathAPI::class)
+        suspend fun fit(
+            spectrum: NBkgSpectrum,
+            data: XYErrorColumnarData<Double, Double, Double>,
+            fitParams: List<Symbol>,
+            initial: Map<Symbol, Double>,
+            logger: Loggable
+        ): XYFit {
+            return data.fitWith(
+                optimizer = QowOptimizer, modelExpression = spectrum, startingPoint = initial, attributesBuilder = {
+                    freeParameters(*fitParams.toTypedArray())
+                    iterations(20)
+                    OptimizationLog(logger)
+                })
+        }
+    }
+}

numass-workspace/src/main/kotlin/ru/inr/mass/scripts/fit-custom.kt

@@ -12,110 +12,68 @@ import com.github.ajalt.clikt.parameters.options.option
 import com.github.ajalt.clikt.parameters.types.double
 import com.github.ajalt.clikt.parameters.types.enum
 import kotlinx.coroutines.runBlocking
-import kotlinx.html.*
+import kotlinx.html.h3
+import kotlinx.html.h4
+import kotlinx.html.p
+import kotlinx.html.pre
 import ru.inr.mass.models.*
+import ru.inr.mass.scripts.SterileNeutrinoSpectrumRWall.Companion.rearWall
 import ru.inr.mass.workspace.buffer
-import ru.inr.mass.workspace.fitWith
-import ru.inr.mass.workspace.freeParameters
-import ru.inr.mass.workspace.iterations
 import space.kscience.kmath.UnstableKMathAPI
 import space.kscience.kmath.data.XYErrorColumnarData
 import space.kscience.kmath.data.indices
 import space.kscience.kmath.expressions.Symbol
-import space.kscience.kmath.expressions.symbol
 import space.kscience.kmath.misc.Loggable
+import space.kscience.kmath.operations.toTypedArray
 import space.kscience.kmath.optimization.*
-import space.kscience.kmath.structures.toDoubleArray
 import space.kscience.plotly.*
 import space.kscience.plotly.models.ScatterMode
-import java.io.File
-import java.nio.file.Path
 import kotlin.io.path.Path
 import kotlin.io.path.nameWithoutExtension
-import kotlin.math.pow
 
-val rearWall: Symbol by symbol
-
-
-class SterileNeutrinoSpectrumRWallMP(
-    source: DifferentiableKernel = NumassBeta,
-    transmission: DifferentiableKernel = NumassTransmission(),
-    resolution: DifferentiableKernel = NumassResolution(),
-    val wallFunc: WallModel,
-    fss: FSS? = null
-) : SterileNeutrinoSpectrumFast(source, transmission, resolution, fss) {
-
-    override fun invoke(u: Double, arguments: Map<Symbol, Double>): Double {
-        val rearWall = arguments[rearWall]!!
-        return super.invoke(u, arguments) * (1.0 + rearWall * wallFunc(u))
-    }
-
-    override fun derivativeOrNull(symbols: List<Symbol>): Spectrum? {
-        if (symbols.isEmpty()) return this
-        return when (symbols.singleOrNull() ?: TODO("First derivatives only")) {
-            rearWall -> Spectrum { u, arguments ->
-                wallFunc(u) * super.invoke(u, arguments)
-            }
-
-            else -> super.derivativeOrNull(symbols)
-        }
-    }
-}
-
-private class CustomArgs : CliktCommand() {
+private class FitCustom : CliktCommand(name = "fit-custom") {
 
     init {
-        context {
-            helpFormatter = { MordantHelpFormatter(it, showDefaultValues = true) }
-        }
+        context { helpFormatter = { MordantHelpFormatter(it, showDefaultValues = true) } }
     }
 
-    val spectrum: String by argument().help("path to spectrum TSV file")
-    val full: String? by option().help("path to full spectrum (for test)")
-    val postfix: String? by option().help("save file postfix")
+    val spectrum by argument().help("path to spectrum TSV file")
+    val full by option().help("path to full spectrum (for residuals/test)")
+    val postfix by option().help("save file postfix")
 
-    val norm: Double by option().double().default(8.311751921319484E8).help("NBkgSpectrum.norm")
+    val norm by option().double().default(8.311751921319484E8).help("NBkgSpectrum.norm")
+    val bkg by option().double().default(1612.626961212948).help("NBkgSpectrum.bkg")
+    val fixBkg by option().flag().help("Don't fit bkg")
 
-    val bkg: Double by option().double().default(1612.626961212948).help("NBkgSpectrum.bkg")
-    val fixBkg: Boolean by option().flag().help("Don't fit bkg variable")
+    val mnu2 by option().double().default(0.0).help("NumassBeta.mnu2")
 
-    val mnu2: Double by option().double().default(0.0).help("NumassBeta.mnu2")
+    val e0 by option().double().default(18572.0).help("NumassBeta.e0")
+    val fixE0 by option().flag().help("Don't fit e0")
 
-    val e0: Double by option().double().default(18572.0).help("NumassBeta.e0")
-    val fixE0: Boolean by option().flag().help("Don't fit e0 variable")
+    val msterile2 by option().double().default(0.0).help("NumassBeta.msterile2")
 
-    val msterile2: Double by option().double().default(0.0.pow(2)).help("NumassBeta.msterile2")
+    val u2 by option().double().default(1e-2).help("NumassBeta.u2")
+    val fixU2 by option().flag().help("Don't fit u2")
 
-    val u2: Double by option().double().default(1e-2).help("NumassBeta.u2")
-    val fixU2: Boolean by option().flag().help("Don't fit u2 variable")
+    val thickness by option().double().default(0.3).help("NumassTransmission.thickness")
 
-    val thickness: Double by option().double().default(0.3).help("NumassTransmission.thickness")
+    val trap by option().double().default(1.0).help("NumassTransmission.trap")
+    val trapFunc by option().enum<TrapModel>().default(TrapModel.FINE_5DAY).help("trap model")
+    val fixTrap by option().flag().help("Don't fit trap")
 
-    val trap: Double by option().double().default(1.0).help("NumassTransmission.trap")
-    val trapFunc: TrapModel by option().enum<TrapModel>().default(TrapModel.FINE_5DAY)
-        .help(
-            "trap model"
-        )
-    val fixTrap: Boolean by option().flag().help("Don't fit trap variable")
+    val rwall by option().double().default(0.17).help("rear wall factor")
+    val wallFunc by option().enum<WallModel>().default(WallModel.TSV_2025_12_01).help("wall model")
+    val fixRwall by option().flag().help("Don't fit rwall")
 
-
-    val rwall: Double by option().double().default(0.17).help("rear wall")
-    val wallFunc: WallModel by option().enum<WallModel>().default(WallModel.TSV_2025_12_01)
-        .help(
-            "wall model"
-        )
-    val fixRwall: Boolean by option().flag().help("Don't fit rwall variable")
-
-    // TODO: convert to enum
-    val adiabacityFunc: AdiabaticityModel by option().enum<AdiabaticityModel>().default(AdiabaticityModel.`2024_11`)
-        .help(
-            "adiabacity model"
-        )
+    val adiabacityFunc by option().enum<AdiabaticityModel>().default(AdiabaticityModel.`2024_11`)
+        .help("adiabacity tail model")
 
     @UnstableKMathAPI
-    override fun run() {
+    override fun run() = runBlocking {
+        val data = parse(spectrum)
+        val testData = full?.let { parse(it) } ?: data
 
-        val fitParams: Map<Symbol, Double> = mapOf(
+        val initialParams = mapOf(
             NBkgSpectrum.norm to norm,
             NBkgSpectrum.bkg to bkg,
             NumassBeta.mnu2 to mnu2,
@@ -127,276 +85,110 @@ private class CustomArgs : CliktCommand() {
             rearWall to rwall
         )
 
-        runBlocking {
-            processCustom(
-                spectrum,
-                full,
-                postfix,
-                fitParams,
-                fixBkg,
-                fixU2,
-                fixTrap,
-                fixE0,
-                trapFunc,
-                fixRwall,
-                wallFunc,
-                cliArgs = this@CustomArgs,
-                adiabacityFunc = adiabacityFunc
-            )
+        val spectrumModel = SterileNeutrinoSpectrumRWall(
+            wallModel = wallFunc, trapModel = trapFunc, adiabaticityModel = adiabacityFunc, parallel = true
+        ).withNBkg()
+
+        val fitVars = mutableListOf(NBkgSpectrum.norm)
+        if (!fixE0) fitVars += NumassBeta.e0
+        if (!fixBkg) fitVars += NBkgSpectrum.bkg
+        if (!fixU2) fitVars += NumassBeta.u2
+        if (!fixTrap) fitVars += NumassTransmission.trap
+        if (!fixRwall) fitVars += rearWall
+
+        val logMessages = mutableListOf<String>()
+        val logger = Loggable { tag, msg ->
+            val line = "[$tag] ${msg()}"
+            logMessages += line
+            println(line)
         }
+
+        val fit = SterileNeutrinoSpectrumRWall.fit(
+            spectrumModel,
+            data,
+            fitVars,
+            initialParams,
+            logger
+        )
+
+        val dataPath = Path(spectrum)
+        val outName = dataPath.nameWithoutExtension + (postfix?.let { "-$it" } ?: "") + ".html"
+        val reportPath = dataPath.parent.resolve(outName)
+
+        buildReportPage(
+            spectrumFile = spectrum,
+            postfix = postfix,
+            data = data,
+            testData = testData,
+            model = spectrumModel,
+            initial = initialParams,
+            fit = fit,
+            logMessages = logMessages,
+            cliArgs = this@FitCustom
+        ).makeFile(reportPath, show = false)
+
+        println("Report saved → $reportPath")
+        openInBrowser("file:///$reportPath")
     }
 }
 
-suspend fun fitNumassSpectrumRWall(
-    spectrum: NBkgSpectrum,
-    data: XYErrorColumnarData<Double, Double, Double>,
-    fitParams: List<Symbol>,
-    initial: Map<Symbol, Double>,
-    logger: Loggable
-): XYFit {
-    return data.fitWith(
-        optimizer = QowOptimizer, modelExpression = spectrum, startingPoint = initial, attributesBuilder = {
-            freeParameters(*fitParams.toTypedArray())
-            iterations(20)
-            OptimizationLog(logger)
-        })
-}
-
 @UnstableKMathAPI
-suspend fun processCustom(
+fun buildReportPage(
     spectrumFile: String,
-    full: String?,
     postfix: String?,
-    fitParams: Map<Symbol, Double>,
-    fixBkg: Boolean,
-    fixU2: Boolean,
-    fixTrap: Boolean,
-    fixE0: Boolean,
-    trapFunc: ITrap,
-    fixRwall: Boolean,
-    wallFunc: WallModel,
-    cliArgs: Any? = null,
-    adiabacityFunc: AdiabaticityModel
-) {
-    println(cliArgs)
-
-    val postfix = if (postfix != null) "-$postfix" else ""
-
-    val data = parse(spectrumFile)
-    val testData = if (full != null) parse(full) else data
-
-    val spectrum: NBkgSpectrum = SterileNeutrinoSpectrumRWallMP(
-//        fss = FSS.default,
-        wallFunc = wallFunc,
-        transmission = NumassTransmission(
-//            NumassTransmission.trapFunction,
-            trapFunc = { ei, ef, _ ->
-                val delta = ei - ef
-                trapFunc.value(ei, delta)
-            },
-            adjustX = false,
-        ), resolution = NumassResolution(1.7e-4, tailFunction = { e, u ->
-            adiabacityFunc(e, u)
-        })
-    ).withNBkg()
-
-    val logMessages = emptyList<String>().toMutableList()
-    val dumpLogger = Loggable { tag, block ->
-        logMessages += "[$tag] ${block()}"
-        println("[$tag] ${block()}")
-    }
-
-    val fitVars = mutableListOf(
-        NBkgSpectrum.norm,
-//        NBkgSpectrum.bkg,
-//        rearWall,
-//        NumassBeta.u2,
-//        NumassBeta.mnu2,
-//        NumassBeta.msterile2,
-//        NumassTransmission.trap,
-    )
-
-    if (!fixE0) {
-        fitVars += NumassBeta.e0
-    }
-
-    if (!fixBkg) {
-        fitVars += NBkgSpectrum.bkg
-    }
-
-    if (!fixU2) {
-        fitVars += NumassBeta.u2
-    }
-
-    if (!fixTrap) {
-        fitVars += NumassTransmission.trap
-    }
-
-    if (!fixRwall) {
-        fitVars += rearWall
-    }
-
-    val fit = fitNumassSpectrumRWall(spectrum, data, fitVars, fitParams, dumpLogger)
-
-    val dataPath = Path(spectrumFile)
-    val reportPath = Path(
-        dataPath.parent.toString(), dataPath.nameWithoutExtension + postfix + ".html"
-    )
-
-    val report =
-        reportPage(spectrumFile, postfix, data, spectrum, fitParams, fit, testData, dataPath, cliArgs, logMessages)
-
-    report.makeFile(path = reportPath, show = false)
-    openInBrowser("file:///$reportPath")
-}
-
-@UnstableKMathAPI
-fun reportPage(
-    spectrumFile: String,
-    postfix: String,
     data: XYErrorColumnarData<Double, Double, Double>,
-    spectrum: NBkgSpectrum,
-    fitParams: Map<Symbol, Double>,
-    fit: XYFit,
     testData: XYErrorColumnarData<Double, Double, Double>,
-    dataPath: Path,
-    cliArgs: Any?,
-    logMessages: MutableList<String>
-): PlotlyPage {
-    val report = Plotly.page {
-        h3 {
-            +"Fit for $spectrumFile${postfix.drop(1).takeIf { it.isNotEmpty() }?.let { " ($it)" } ?: ""}"
-        }
-        plot {
-            scatter {
-                name = "Data"
-                mode = ScatterMode.markers
-                x.buffer = data.x
-                y.buffer = data.y
-            }
-            scatter {
-                name = "Initial"
-                mode = ScatterMode.lines
-                x.buffer = data.x
-                y.numbers = x.doubles.map { spectrum(it, fitParams + fit.result) }
-            }
-            scatter {
-                name = "Fit"
-                mode = ScatterMode.lines
-                x.buffer = data.x
-                y.numbers = data.x.toDoubleArray().map { spectrum(it, fitParams + fit.result) }
-            }
-        }
-        plot {
-            layout {
-                title = "Residuals"
-            }
-            scatter {
-                name = "Residuals"
-                mode = ScatterMode.markers
-                x.buffer = testData.x
-                y.numbers = testData.indices.map {
-                    val value = spectrum(testData.x[it], fitParams + fit.result)
-                    val dif = testData.y[it] - value
-                    dif / testData.yErr[it]
-                }
+    model: NBkgSpectrum,
+    initial: Map<Symbol, Double>,
+    fit: XYFit,
+    logMessages: List<String>,
+    cliArgs: Any?
+): PlotlyPage = Plotly.page {
+    h3 { +"Fit for $spectrumFile${postfix?.let { " ($it)" } ?: ""}" }
 
-                File(dataPath.parent.toString(), dataPath.nameWithoutExtension + postfix + ".fit.tsv").printWriter()
-                    .use { out ->
-                        out.println("U_sp\tcounts\tresidials")
-                        testData.indices.forEach {
-                            val value = spectrum(testData.x[it], fitParams + fit.result)
-                            val dif = testData.y[it] - value
-                            val residial = dif / testData.yErr[it]
-                            out.println("${testData.x[it]}\t${value}\t$residial")
-                        }
-                    }
-            }
+    plot {
+        scatter { name = "Data"; mode = ScatterMode.markers; x.buffer = data.x; y.buffer = data.y }
+        scatter {
+            name = "Initial"; mode = ScatterMode.lines; x.buffer = data.x; y.numbers =
+            data.x.toTypedArray().map { model(it, initial) }
         }
-        plot {
-            layout {
-                title = "Residuals distribution"
-            }
-            histogram {
-                x.numbers = data.indices.map {
-                    val value = spectrum(data.x[it], fitParams + fit.result)
-                    val dif = data.y[it] - value
-                    dif / data.yErr[it]
-                }
-                name = "Res histo"
-            }
+        scatter {
+            name = "Fit"; mode = ScatterMode.lines; x.buffer = data.x; y.numbers =
+            data.x.toTypedArray().map { model(it, initial + fit.result) }
         }
-        p {
-            +"command args = $cliArgs"
-        }
-        br()
-        p {
-            +"initial params = $fitParams"
-        }
-        br()
-        p {
-            +"fit result = ${fit.result}"
-        }
-        br()
-        p {
-            +"chi^2 = ${fit.chiSquaredOrNull}"
-        }
-        br()
-        p {
-            +"dof = ${fit.dof}"
-        }
-        h4 {
-            +"Fit log:"
-        }
-        pre {
-            logMessages.forEach {
-                +it
-                +"\n"
-            }
-        }
-        h4 {
-            +"Residuals:"
-        }
-        table {
-            tr {
-                th { +"U_sp" }
-                th { +"Counts" }
-                th { +"residials" }
-            }
-            testData.indices.forEach {
-                tr {
-                    td { +testData.x[it].toString() }
-                    td {
-                        val value = spectrum(testData.x[it], fitParams + fit.result)
-                        +value.toString()
-                    }
-                    td {
-                        val value = spectrum(testData.x[it], fitParams + fit.result)
-                        val dif = testData.y[it] - value
-                        +(dif / testData.yErr[it]).toString()
-                    }
-                }
-            }
-        }
-        h4 {
-            +"Data:"
-        }
-        table {
-            tr {
-                th { +"U_sp" }
-                th { +"counts" }
-                th { +"residials" }
-            }
-            data.indices.forEach {
-                tr {
-                    td { +data.x[it].toString() }
-                    td { +data.y[it].toString() }
-                    td { +data.yErr[it].toString() }
-                }
+    }
+
+    plot {
+        layout.title = "Residuals (σ)"
+        scatter {
+            name = "Residuals"
+            mode = ScatterMode.markers
+            x.buffer = testData.x
+            y.numbers = testData.indices.map {
+                val pred = model(testData.x[it], initial + fit.result)
+                (testData.y[it] - pred) / testData.yErr[it]
             }
         }
     }
-    return report
+
+    plot {
+        layout.title = "Residuals distribution"
+        histogram {
+            name = "Res histo"
+            x.numbers = data.indices.map {
+                val pred = model(data.x[it], initial + fit.result)
+                (data.y[it] - pred) / data.yErr[it]
+            }
+        }
+    }
+
+    p { +"CLI: $cliArgs" }
+    p { +"Initial: $initial" }
+    p { +"Result: ${fit.result}" }
+    p { +"χ² = ${fit.chiSquaredOrNull ?: "—"},  dof = ${fit.dof}" }
+
+    h4 { +"Log" }
+    pre { logMessages.forEach { +"$it\n" } }
 }
 
-fun main(args: Array<String>) = CustomArgs().main(args)
+fun main(args: Array<String>) = FitCustom().main(args)

numass-workspace/src/main/kotlin/ru/inr/mass/scripts/get-spectrum.kt

@@ -12,6 +12,7 @@ import com.github.ajalt.clikt.parameters.types.enum
 import kotlinx.coroutines.Dispatchers
 import kotlinx.coroutines.runBlocking
 import ru.inr.mass.models.*
+import ru.inr.mass.scripts.SterileNeutrinoSpectrumRWall.Companion.rearWall
 import space.kscience.kmath.UnstableKMathAPI
 import space.kscience.kmath.expressions.Symbol
 import space.kscience.kmath.operations.toTypedArray
@@ -55,12 +56,6 @@ private class CalcSpectrum : CliktCommand() {
             "wall model"
         )
 
-    val fixRwall: Boolean by option().flag().help("Don't fit rwall variable")
-
-    val noAdiabacity: Boolean by option().flag().help("Use const 1 instead of adiabacity function")
-
-    val fixE0: Boolean by option().flag().help("Don't fit e0 variable")
-
     // TODO: convert to enum
     val adiabacityFunc: AdiabaticityModel by option().enum<AdiabaticityModel>().default(AdiabaticityModel.`2024_11`)
         .help(
@@ -94,22 +89,15 @@ suspend fun calcSpectrum(fitParams: Map<Symbol, Double>, trapFunc: ITrap, wallFu
 
     runBlocking(Dispatchers.Default) {
 
-        val spectrum: NBkgSpectrum = SterileNeutrinoSpectrumRWallMP(
-            wallFunc = wallFunc,
-            transmission = NumassTransmission(
-                trapFunc = { ei, ef, _ ->
-                    val delta = ei - ef
-                    trapFunc.value(ei, delta)
-                },
-                adjustX = false,
-            ),
-            resolution = NumassResolution(1.7e-4, tailFunction = { e, u -> AdiabaticityModel.`2024_11`(e, u) })
+        val spectrum: NBkgSpectrum = SterileNeutrinoSpectrumRWall(
+            wallModel = wallFunc,
+            trapFunc,
+            AdiabaticityModel.`2024_11`,
+            parallel = true,
         ).withNBkg()
 
         val out = range.map {
 //            async(Dispatchers.Default) {
-
-
             val res = spectrum.invoke(it, fitParams)
 //                return@async res
 //            }

numass-workspace/src/main/kotlin/ru/inr/mass/scripts/grid-fit.kt

@@ -12,10 +12,9 @@ import com.github.ajalt.clikt.parameters.options.option
 import com.github.ajalt.clikt.parameters.options.multiple
 import com.github.ajalt.clikt.parameters.types.double
 import com.github.ajalt.clikt.parameters.types.enum
-import kotlinx.coroutines.coroutineScope
-import kotlinx.coroutines.launch
 import kotlinx.coroutines.runBlocking
 import ru.inr.mass.models.*
+import ru.inr.mass.scripts.SterileNeutrinoSpectrumRWall.Companion.rearWall
 import space.kscience.kmath.UnstableKMathAPI
 import space.kscience.kmath.expressions.Symbol
 import space.kscience.kmath.misc.Loggable
@@ -27,33 +26,6 @@ import kotlin.io.path.createDirectories
 import kotlin.io.path.nameWithoutExtension
 import kotlin.math.log10
 import kotlin.math.pow
-import kotlinx.coroutines.Dispatchers
-import kotlinx.coroutines.withContext
-
-class SterileNeutrinoSpectrumRWall(
-    source: DifferentiableKernel = NumassBeta,
-    transmission: DifferentiableKernel = NumassTransmission(),
-    resolution: DifferentiableKernel = NumassResolution(),
-    val wallFunc: WallModel,
-    fss: FSS? = null
-) : SterileNeutrinoSpectrum(source, transmission, resolution, fss) {
-
-    override fun invoke(u: Double, arguments: Map<Symbol, Double>): Double {
-        val rearWall = arguments[rearWall]!!
-        return super.invoke(u, arguments) * (1.0 + rearWall * wallFunc(u))
-    }
-
-    override fun derivativeOrNull(symbols: List<Symbol>): Spectrum? {
-        if (symbols.isEmpty()) return this
-        return when (symbols.singleOrNull() ?: TODO("First derivatives only")) {
-            rearWall -> Spectrum { u, arguments ->
-                wallFunc(u) * super.invoke(u, arguments)
-            }
-
-            else -> super.derivativeOrNull(symbols)
-        }
-    }
-}
 
 private class GridArgs : CliktCommand() {
 
@@ -122,16 +94,8 @@ private class GridArgs : CliktCommand() {
             rearWall to rwall
         )
 
-        val spectrumModel: NBkgSpectrum = SterileNeutrinoSpectrumRWallMP(
-            wallFunc = wallFunc, transmission = NumassTransmission(
-                trapFunc = { ei, ef, _ ->
-                    val delta = ei - ef
-                    trapFunc.value(ei, delta)
-                },
-                adjustX = false,
-            ), resolution = NumassResolution(1.7e-4, tailFunction = { e, u ->
-                adiabacityFunc(e, u)
-            })
+        val spectrumModel: NBkgSpectrum = SterileNeutrinoSpectrumRWall(
+            wallModel = wallFunc, trapModel = trapFunc, adiabaticityModel = adiabacityFunc, parallel = true
         ).withNBkg()
 
         val fitVars = mutableListOf(NBkgSpectrum.norm)
@@ -231,7 +195,7 @@ private class GridArgs : CliktCommand() {
                         println("$current/$total [$taskId] $msg")
                     }
 
-                    val fit = fitNumassSpectrumRWall(
+                    val fit = SterileNeutrinoSpectrumRWall.fit(
                         spectrumModel, data, fitVars, fitParams, dumpLogger
                     )
 
@@ -243,17 +207,16 @@ private class GridArgs : CliktCommand() {
                         val reportFileName = "${mStr}-${u2Str}.html"
                         val reportPath = reportsDir.resolve(reportFileName)
 
-                        reportPage(
+                        buildReportPage(
                             spectrum,
                             postfix,
                             data,
+                            testData = dataForView,
                             spectrumModel,
                             fitParams,
                             fit,
-                            dataForView,
-                            dataPath,
+                            logMessages,
                             this@GridArgs,
-                            logMessages
                         ).makeFile(path = reportPath, show = false)
                     }
 

numass-workspace/src/main/kotlin/ru/inr/mass/scripts/trap-compare-rearwall.kt

@@ -11,6 +11,7 @@ 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.scripts.SterileNeutrinoSpectrumRWall.Companion.rearWall
 import ru.inr.mass.workspace.buffer
 import space.kscience.kmath.expressions.Symbol
 import space.kscience.kmath.real.step
@@ -19,8 +20,9 @@ import space.kscience.plotly.models.ScatterMode
 
 fun main() {
 
-    val trapInterpolator = TrapModel.FINE_5DAY
-    val wallInterpolator = WallModel.TSV_2025_12_01
+    val trapModel = TrapModel.FINE_5DAY
+    val wallModel = WallModel.TSV_2025_12_01
+    val adiabacityModel = AdiabaticityModel.vasya
 
     val range = 2000.0..18600.0 step 50.0
 
@@ -45,19 +47,8 @@ fun main() {
     val rwall3 = args.toMutableMap();
     rwall3[rearWall] = 3.0
 
-    val spectrumC: NBkgSpectrum = SterileNeutrinoSpectrumRWallMP(
-        fss = null,
-        transmission = NumassTransmission(
-            trapFunc = { ei, ef, _ ->
-                val delta = ei - ef
-                trapInterpolator.value(ei, delta)
-            },
-            adjustX = false,
-        ),
-        wallFunc = wallInterpolator,
-        resolution = NumassResolution(1.7e-4, tailFunction = {
-                e,u -> AdiabaticityModel.vasya(e,u)
-        })
+    val spectrumC: NBkgSpectrum = SterileNeutrinoSpectrumRWall(
+        wallModel, trapModel, adiabacityModel
     ).withNBkg()
 
     Plotly.page {
@@ -86,11 +77,11 @@ fun main() {
         }
         br()
         p {
-            +"trapping model = $trapInterpolator"
+            +"trapping model = $trapModel"
         }
         br()
         p {
-            +"wall model = $wallInterpolator"
+            +"wall model = $wallModel"
         }
         br()
     }.makeFile()
@@ -103,6 +94,12 @@ fun main() {
         val r2 = spectrumC(it, rwall2)
         val r3 = spectrumC(it, rwall3)
 
-        println("${it}\t${spectrumC(it, args)}\t${spectrumC(it, rwall1)}\t${spectrumC(it, rwall2)}\t${spectrumC(it, rwall3)}")
+        println(
+            "${it}\t${spectrumC(it, args)}\t${spectrumC(it, rwall1)}\t${spectrumC(it, rwall2)}\t${
+                spectrumC(
+                    it, rwall3
+                )
+            }"
+        )
     }
 }

numass-workspace/src/main/kotlin/ru/inr/mass/scripts/trap-compare-spectrums.kt

@@ -9,6 +9,7 @@ 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.scripts.SterileNeutrinoSpectrumRWall.Companion.rearWall
 import ru.inr.mass.workspace.buffer
 import space.kscience.kmath.expressions.Symbol
 import space.kscience.kmath.real.step
@@ -16,8 +17,6 @@ import space.kscience.plotly.*
 import space.kscience.plotly.models.ScatterMode
 
 fun main() {
-
-    val trapInterpolator = TrapModel.FINE_5DAY
     val range = 2000.0..18600.0 step 50.0
 
     val args: Map<Symbol, Double> = mapOf(
@@ -41,19 +40,8 @@ fun main() {
     val argsTrap3 = args.toMutableMap();
     argsTrap3[trap] = 3.0
 
-    val spectrumC: NBkgSpectrum = SterileNeutrinoSpectrumRWallMP(
-        fss = null,
-        transmission = NumassTransmission(
-            trapFunc = { ei, ef, _ ->
-                val delta = ei - ef
-                trapInterpolator.value(ei, delta)
-            },
-            adjustX = false,
-        ),
-        wallFunc = WallModel.TSV_2025_12_01,
-        resolution = NumassResolution(1.7e-4, tailFunction = {
-                e,u -> AdiabaticityModel.vasya(e,u)
-        })
+    val spectrumC: NBkgSpectrum = SterileNeutrinoSpectrumRWall(
+        WallModel.TSV_2025_12_01, TrapModel.FINE_5DAY, AdiabaticityModel.vasya
     ).withNBkg()
 
     Plotly.page {
@@ -81,6 +69,12 @@ fun main() {
 
     println("HV\ttrap=${args[trap]}\ttrap=${argsTrap1[trap]}\ttrap=${argsTrap2[trap]}\ttrap=${argsTrap3[trap]}")
     range.iterator().forEach {
-        println("${it}\t${spectrumC(it, args)}\t${spectrumC(it, argsTrap1)}\t${spectrumC(it, argsTrap2)}\t${spectrumC(it, argsTrap3)}")
+        println(
+            "${it}\t${spectrumC(it, args)}\t${spectrumC(it, argsTrap1)}\t${
+                spectrumC(
+                    it, argsTrap2
+                )
+            }\t${spectrumC(it, argsTrap3)}"
+        )
     }
 }

numass-workspace/src/main/kotlin/ru/inr/mass/scripts/wall-compare-spectrums.kt

@@ -9,6 +9,7 @@ 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.scripts.SterileNeutrinoSpectrumRWall.Companion.rearWall
 import ru.inr.mass.workspace.buffer
 import space.kscience.kmath.expressions.Symbol
 import space.kscience.kmath.real.step
@@ -35,19 +36,8 @@ fun main() {
     val argsRWall = args.toMutableMap();
     argsRWall[rearWall] = 0.04
 
-    val spectrumC: NBkgSpectrum = SterileNeutrinoSpectrumRWallMP(
-        fss = null,
-        transmission = NumassTransmission(
-            trapFunc = { ei, ef, _ ->
-                val delta = ei - ef
-                trapInterpolator.value(ei, delta)
-            },
-            adjustX = false,
-        ),
-        wallFunc = WallModel.GEANT_2025_03_10,
-        resolution = NumassResolution(1.7e-4, tailFunction = {
-                e,u -> AdiabaticityModel.`2024_11`(e,u)
-        })
+    val spectrumC: NBkgSpectrum = SterileNeutrinoSpectrumRWall(
+        WallModel.GEANT_2025_03_10, trapInterpolator, AdiabaticityModel.`2024_11`
     ).withNBkg()
 
     Plotly.page {