add cli parser to fit.kt

This commit is contained in:
chernov 2023-11-29 12:57:15 +03:00
parent bc77af263b
commit 6650353683
3 changed files with 148 additions and 55 deletions

View File

@ -2,7 +2,7 @@ plugins {
id("space.kscience.gradle.jvm") id("space.kscience.gradle.jvm")
id("com.github.johnrengelman.shadow") version "7.1.2" id("com.github.johnrengelman.shadow") version "7.1.2"
`maven-publish` `maven-publish`
// application application
} }
kotlin { kotlin {
@ -32,5 +32,5 @@ kscience{
} }
//application { //application {
// mainClass.set("ru.inr.mass.scripts.Spectrum_serverKt") // mainClass.set("ru.inr.mass.scripts.Get_spectrumKt")
//} //}

View File

@ -1,5 +1,17 @@
package ru.inr.mass.scripts 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.models.*
import ru.inr.mass.workspace.buffer import ru.inr.mass.workspace.buffer
import ru.inr.mass.workspace.fitWith import ru.inr.mass.workspace.fitWith
@ -37,47 +49,79 @@ fun parse(filename: String): XYErrorColumnarData<Double, Double, Double> {
errors.map { it }.asBuffer() errors.map { it }.asBuffer()
) )
} }
class Args : CliktCommand() {
//@OptIn(UnstableKMathAPI::class) init {
@UnstableKMathAPI context {
suspend fun main() { 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<Symbol, Double> = 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<String>) = Args().main(args)
suspend fun fitNumassSpectrum(data: XYErrorColumnarData<Double, Double, Double>, initial: Map<Symbol, Double>): XYFit {
val spectrum: NBkgSpectrum = SterileNeutrinoSpectrum( val spectrum: NBkgSpectrum = SterileNeutrinoSpectrum(
// fss = FSS.default, // fss = FSS.default,
transmission = NumassTransmission(NumassTransmission.trapFunction), transmission = NumassTransmission(NumassTransmission.trapFunction),
// resolution = NumassResolution(1.7e-4) // resolution = NumassResolution(1.7e-4)
).withNBkg() ).withNBkg()
val args: Map<Symbol, Double> = mapOf( return data.fitWith(
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(
optimizer = QowOptimizer, optimizer = QowOptimizer,
modelExpression = spectrum, modelExpression = spectrum,
startingPoint = args, startingPoint = initial,
OptimizationParameters(NumassBeta.e0, NBkgSpectrum.norm, NBkgSpectrum.bkg), OptimizationParameters(NumassBeta.e0, NBkgSpectrum.norm, NBkgSpectrum.bkg,
// NumassTransmission.trap
),
OptimizationIterations(20) OptimizationIterations(20)
) )
}
@UnstableKMathAPI
suspend fun process(spectrumFile: String, fitParams: Map<Symbol, Double>) {
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) println(fit.resultPoint)
@ -95,11 +139,15 @@ suspend fun main() {
scatter { scatter {
name = "Initial" name = "Initial"
mode = ScatterMode.lines mode = ScatterMode.lines
x.buffer = 12000.0..18600.0 step 50.0 x.buffer = data.x
y.numbers = x.doubles.map { spectrum(it, args + fit.resultPoint) } y.numbers = x.doubles.map { spectrum(it, fitParams + fit.resultPoint) }
File("/home/chernov/projects/notebooks/manual-fit/04.10.2023-fit").printWriter().use { out -> File(spectrumFile + "-fit").printWriter().use { out ->
y.numbers.forEach { out.println("U_sp\tcounts\tresidials")
out.println(it)
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" name = "Fit"
mode = ScatterMode.lines mode = ScatterMode.lines
x.buffer = 12000.0..18600.0 step 10.0 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{ plot{
@ -116,7 +164,7 @@ suspend fun main() {
mode = ScatterMode.markers mode = ScatterMode.markers
x.buffer = data.x x.buffer = data.x
y.numbers = data.indices.map{ 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 val dif = data.y[it] - value
dif/data.yErr[it] dif/data.yErr[it]
} }
@ -125,7 +173,7 @@ suspend fun main() {
plot{ plot{
histogram { histogram {
x.numbers = data.indices.map{ 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 val dif = data.y[it] - value
dif/data.yErr[it] dif/data.yErr[it]
} }
@ -138,7 +186,7 @@ suspend fun main() {
mode = ScatterMode.lines mode = ScatterMode.lines
x.buffer = data.x x.buffer = data.x
y.numbers = data.indices.map { 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 val dif = data.y[it] - value
dif / data.yErr[it] dif / data.yErr[it]
} }

View File

@ -1,31 +1,76 @@
package ru.inr.mass.scripts 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.Dispatchers
import kotlinx.coroutines.async import kotlinx.coroutines.async
import kotlinx.coroutines.runBlocking import kotlinx.coroutines.runBlocking
import ru.inr.mass.models.* 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.expressions.Symbol
import space.kscience.kmath.operations.toTypedArray import space.kscience.kmath.operations.toTypedArray
import space.kscience.kmath.real.plus import space.kscience.kmath.real.plus
import space.kscience.kmath.real.step import space.kscience.kmath.real.step
import space.kscience.plotly.models.HoverMode
import kotlin.math.pow
fun main(args: Array<String>) { private class CalcSpectrum : CliktCommand() {
val cli = args.map { it.toDouble() }
val args: Map<Symbol, Double> = mapOf( init {
NBkgSpectrum.norm to cli[0], context {
NBkgSpectrum.bkg to cli[1], helpFormatter = { MordantHelpFormatter(it, showDefaultValues = true) }
NumassBeta.mnu2 to cli[2], }
NumassBeta.e0 to cli[3], }
NumassBeta.msterile2 to cli[4],
NumassBeta.u2 to cli[5], val spectrum: String by argument().help("path to spectrum TSV file")
NumassTransmission.thickness to cli[6],
NumassTransmission.trap to cli[7] 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<Symbol, Double> = 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<String>) = CalcSpectrum().main(args)
suspend fun calcSpectrum(spectrumFile: String, fitParams: Map<Symbol, Double>) {
val data = parse(spectrumFile)
val range = data.x.toTypedArray();
runBlocking(Dispatchers.Default) { 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 { val handles = range.map {
@ -34,7 +79,7 @@ fun main(args: Array<String>) {
transmission = NumassTransmission(NumassTransmission.trapFunction), transmission = NumassTransmission(NumassTransmission.trapFunction),
).withNBkg() ).withNBkg()
val res = spectrum.invoke(it, args) val res = spectrum.invoke(it, fitParams)
return@async res return@async res
} }
} }