From 6650353683d1bd2fc6d90ff546ccdd88b298c877 Mon Sep 17 00:00:00 2001 From: chernov Date: Wed, 29 Nov 2023 12:57:15 +0300 Subject: [PATCH] add cli parser to fit.kt --- numass-workspace/build.gradle.kts | 4 +- .../main/kotlin/ru/inr/mass/scripts/fit.kt | 124 ++++++++++++------ .../ru/inr/mass/scripts/get-spectrum.kt | 75 ++++++++--- 3 files changed, 148 insertions(+), 55 deletions(-) diff --git a/numass-workspace/build.gradle.kts b/numass-workspace/build.gradle.kts index 4ed4d30..dff89f5 100644 --- a/numass-workspace/build.gradle.kts +++ b/numass-workspace/build.gradle.kts @@ -2,7 +2,7 @@ plugins { id("space.kscience.gradle.jvm") id("com.github.johnrengelman.shadow") version "7.1.2" `maven-publish` -// application + application } kotlin { @@ -32,5 +32,5 @@ kscience{ } //application { -// mainClass.set("ru.inr.mass.scripts.Spectrum_serverKt") +// mainClass.set("ru.inr.mass.scripts.Get_spectrumKt") //} \ No newline at end of file diff --git a/numass-workspace/src/main/kotlin/ru/inr/mass/scripts/fit.kt b/numass-workspace/src/main/kotlin/ru/inr/mass/scripts/fit.kt index 0fdf198..fd71e94 100644 --- a/numass-workspace/src/main/kotlin/ru/inr/mass/scripts/fit.kt +++ b/numass-workspace/src/main/kotlin/ru/inr/mass/scripts/fit.kt @@ -1,5 +1,17 @@ package ru.inr.mass.scripts +import com.github.ajalt.clikt.core.CliktCommand +import com.github.ajalt.clikt.core.context +import com.github.ajalt.clikt.output.MordantHelpFormatter +import com.github.ajalt.clikt.parameters.arguments.argument +import com.github.ajalt.clikt.parameters.arguments.help +import com.github.ajalt.clikt.parameters.options.default +import com.github.ajalt.clikt.parameters.options.help +import com.github.ajalt.clikt.parameters.options.option +import com.github.ajalt.clikt.parameters.options.prompt +import com.github.ajalt.clikt.parameters.types.double +import com.github.ajalt.clikt.parameters.types.int +import kotlinx.coroutines.runBlocking import ru.inr.mass.models.* import ru.inr.mass.workspace.buffer import ru.inr.mass.workspace.fitWith @@ -37,47 +49,79 @@ fun parse(filename: String): XYErrorColumnarData { errors.map { it }.asBuffer() ) } +class Args : CliktCommand() { -//@OptIn(UnstableKMathAPI::class) -@UnstableKMathAPI -suspend fun main() { + init { + context { + helpFormatter = { MordantHelpFormatter(it, showDefaultValues = true) } + } + } -// val fitParams = FitParams().main(args); + val spectrum: String by argument().help("path to spectrum TSV file") + val norm: Double by option().double().default(8.311751921319484E8).help("NBkgSpectrum.norm") + val bkg: Double by option().double().default(1612.626961212948).help("NBkgSpectrum.bkg") + + val mnu2: Double by option().double().default(0.0).help("NumassBeta.mnu2") + val e0: Double by option().double().default(18572.0).help("NumassBeta.e0") + val msterile2: Double by option().double().default(0.0.pow(2)).help("NumassBeta.msterile2") + val u2: Double by option().double().default(1e-2).help("NumassBeta.u2") + + val thickness: Double by option().double().default(0.3).help("NumassTransmission.thickness") + val trap: Double by option().double().default(1.0).help("NumassTransmission.trap") + + @UnstableKMathAPI + override fun run() { + val fitParams: Map = mapOf( + NBkgSpectrum.norm to norm, + NBkgSpectrum.bkg to bkg, + NumassBeta.mnu2 to mnu2, + NumassBeta.e0 to e0, + NumassBeta.msterile2 to msterile2, + NumassBeta.u2 to u2, + NumassTransmission.thickness to thickness, + NumassTransmission.trap to trap + ) + + runBlocking { + process(spectrum, fitParams) + } + } +} + +fun main(args: Array) = Args().main(args) + +suspend fun fitNumassSpectrum(data: XYErrorColumnarData, initial: Map): XYFit { val spectrum: NBkgSpectrum = SterileNeutrinoSpectrum( // fss = FSS.default, - transmission = NumassTransmission(NumassTransmission.trapFunction), + transmission = NumassTransmission(NumassTransmission.trapFunction), // resolution = NumassResolution(1.7e-4) ).withNBkg() - val args: Map = mapOf( - NBkgSpectrum.norm to 8.84e8, - NBkgSpectrum.bkg to -1.702, - NumassBeta.mnu2 to 0.0, - NumassBeta.e0 to 18572.0, - NumassBeta.msterile2 to 0.0.pow(2), - NumassBeta.u2 to 1e-2, - NumassTransmission.thickness to 0.3, - NumassTransmission.trap to 1.0 - ) - -// listOf(NBkgSpectrum.norm, NBkgSpectrum.bkg, NumassBeta.e0).forEach { -// println("$it: ${spectrum.derivative(it).invoke(14000.0, args)}") -// } - -// val timePerPoint = 30.0 -// -// val strategy = (12000.0..19000.0 step 100.0).asSequence().associateWith { timePerPoint } - - val data = parse("/home/chernov/projects/notebooks/manual-fit/06.10.2023") - - val fit: XYFit = data.fitWith( + return data.fitWith( optimizer = QowOptimizer, modelExpression = spectrum, - startingPoint = args, - OptimizationParameters(NumassBeta.e0, NBkgSpectrum.norm, NBkgSpectrum.bkg), + startingPoint = initial, + OptimizationParameters(NumassBeta.e0, NBkgSpectrum.norm, NBkgSpectrum.bkg, +// NumassTransmission.trap + ), OptimizationIterations(20) ) +} + + +@UnstableKMathAPI +suspend fun process(spectrumFile: String, fitParams: Map) { + + val data = parse(spectrumFile) + + val spectrum: NBkgSpectrum = SterileNeutrinoSpectrum( +// fss = FSS.default, + transmission = NumassTransmission(NumassTransmission.trapFunction), +// resolution = NumassResolution(1.7e-4) + ).withNBkg() + + val fit = fitNumassSpectrum(data, fitParams) println(fit.resultPoint) @@ -95,11 +139,15 @@ suspend fun main() { scatter { name = "Initial" mode = ScatterMode.lines - x.buffer = 12000.0..18600.0 step 50.0 - y.numbers = x.doubles.map { spectrum(it, args + fit.resultPoint) } - File("/home/chernov/projects/notebooks/manual-fit/04.10.2023-fit").printWriter().use { out -> - y.numbers.forEach { - out.println(it) + x.buffer = data.x + y.numbers = x.doubles.map { spectrum(it, fitParams + fit.resultPoint) } + File(spectrumFile + "-fit").printWriter().use { out -> + out.println("U_sp\tcounts\tresidials") + + data.indices.map{ + val value = spectrum(data.x[it], fitParams + fit.resultPoint) + val dif = data.y[it] - value + out.println("${data.x[it]}\t${data.y[it]}\t${dif/data.yErr[it]}") } } } @@ -107,7 +155,7 @@ suspend fun main() { name = "Fit" mode = ScatterMode.lines x.buffer = 12000.0..18600.0 step 10.0 - y.numbers = x.doubles.map { spectrum(it, args + fit.resultPoint) } + y.numbers = x.doubles.map { spectrum(it, fitParams + fit.resultPoint) } } } plot{ @@ -116,7 +164,7 @@ suspend fun main() { mode = ScatterMode.markers x.buffer = data.x y.numbers = data.indices.map{ - val value = spectrum(data.x[it], args + fit.resultPoint) + val value = spectrum(data.x[it], fitParams + fit.resultPoint) val dif = data.y[it] - value dif/data.yErr[it] } @@ -125,7 +173,7 @@ suspend fun main() { plot{ histogram { x.numbers = data.indices.map{ - val value = spectrum(data.x[it], args + fit.resultPoint) + val value = spectrum(data.x[it], fitParams + fit.resultPoint) val dif = data.y[it] - value dif/data.yErr[it] } @@ -138,7 +186,7 @@ suspend fun main() { mode = ScatterMode.lines x.buffer = data.x y.numbers = data.indices.map { - val value = spectrum(data.x[it], args + fit.resultPoint) + val value = spectrum(data.x[it], fitParams + fit.resultPoint) val dif = data.y[it] - value dif / data.yErr[it] } diff --git a/numass-workspace/src/main/kotlin/ru/inr/mass/scripts/get-spectrum.kt b/numass-workspace/src/main/kotlin/ru/inr/mass/scripts/get-spectrum.kt index a556d0b..552fc24 100644 --- a/numass-workspace/src/main/kotlin/ru/inr/mass/scripts/get-spectrum.kt +++ b/numass-workspace/src/main/kotlin/ru/inr/mass/scripts/get-spectrum.kt @@ -1,31 +1,76 @@ package ru.inr.mass.scripts +import com.github.ajalt.clikt.core.CliktCommand +import com.github.ajalt.clikt.core.context +import com.github.ajalt.clikt.output.MordantHelpFormatter +import com.github.ajalt.clikt.parameters.arguments.argument +import com.github.ajalt.clikt.parameters.arguments.help +import com.github.ajalt.clikt.parameters.options.default +import com.github.ajalt.clikt.parameters.options.help +import com.github.ajalt.clikt.parameters.options.option +import com.github.ajalt.clikt.parameters.types.double import kotlinx.coroutines.Dispatchers import kotlinx.coroutines.async import kotlinx.coroutines.runBlocking import ru.inr.mass.models.* +import space.kscience.kmath.data.XYErrorColumnarData +import space.kscience.kmath.UnstableKMathAPI import space.kscience.kmath.expressions.Symbol import space.kscience.kmath.operations.toTypedArray import space.kscience.kmath.real.plus import space.kscience.kmath.real.step +import space.kscience.plotly.models.HoverMode +import kotlin.math.pow -fun main(args: Array) { - val cli = args.map { it.toDouble() } +private class CalcSpectrum : CliktCommand() { - val args: Map = mapOf( - NBkgSpectrum.norm to cli[0], - NBkgSpectrum.bkg to cli[1], - NumassBeta.mnu2 to cli[2], - NumassBeta.e0 to cli[3], - NumassBeta.msterile2 to cli[4], - NumassBeta.u2 to cli[5], - NumassTransmission.thickness to cli[6], - NumassTransmission.trap to cli[7] - ) + init { + context { + helpFormatter = { MordantHelpFormatter(it, showDefaultValues = true) } + } + } + + val spectrum: String by argument().help("path to spectrum TSV file") + + val norm: Double by option().double().default(8.84e8).help("NBkgSpectrum.norm") + val bkg: Double by option().double().default(-1.702).help("NBkgSpectrum.bkg") + + val mnu2: Double by option().double().default(0.0).help("NumassBeta.mnu2") + val e0: Double by option().double().default(18572.0).help("NumassBeta.e0") + val msterile2: Double by option().double().default(0.0.pow(2)).help("NumassBeta.msterile2") + val u2: Double by option().double().default(1e-2).help("NumassBeta.u2") + + val thickness: Double by option().double().default(0.3).help("NumassTransmission.thickness") + val trap: Double by option().double().default(1.0).help("NumassTransmission.trap") + + @UnstableKMathAPI + override fun run() { + val fitParams: Map = mapOf( + NBkgSpectrum.norm to norm, + NBkgSpectrum.bkg to bkg, + NumassBeta.mnu2 to mnu2, + NumassBeta.e0 to e0, + NumassBeta.msterile2 to msterile2, + NumassBeta.u2 to u2, + NumassTransmission.thickness to thickness, + NumassTransmission.trap to trap + ) + + runBlocking { + calcSpectrum(spectrum, fitParams) + } + } +} + +fun main(args: Array) = CalcSpectrum().main(args) + +suspend fun calcSpectrum(spectrumFile: String, fitParams: Map) { + + val data = parse(spectrumFile) + + val range = data.x.toTypedArray(); runBlocking(Dispatchers.Default) { - val range = (12000.0..16000.0 step 50.0).toTypedArray() - .plus((16100.0..18600.0 step 100.0).toTypedArray()); val handles = range.map { @@ -34,7 +79,7 @@ fun main(args: Array) { transmission = NumassTransmission(NumassTransmission.trapFunction), ).withNBkg() - val res = spectrum.invoke(it, args) + val res = spectrum.invoke(it, fitParams) return@async res } }