Fix derivative for the spectrum and introduce fit example

This commit is contained in:
Alexander Nozik 2023-01-24 21:07:04 +03:00
parent 26f5e7461c
commit 43d79db6bb
11 changed files with 42 additions and 33 deletions

View File

@ -1,4 +1,5 @@
import space.kscience.gradle.*
import space.kscience.gradle.useApache2Licence
import space.kscience.gradle.useSPCTeam
plugins {
id("space.kscience.gradle.project")
@ -17,7 +18,7 @@ allprojects {
val dataforgeVersion by extra("0.6.0-dev-15")
val tablesVersion: String by extra("0.2.0-dev-3")
val kmathVersion by extra("0.3.1-dev-6")
val kmathVersion by extra("0.3.1-dev-9")
val visionForgeVersion: String by rootProject.extra("0.3.0-dev-6")

View File

@ -13,4 +13,4 @@ org.gradle.parallel=true
org.gradle.jvmargs=-XX:MaxMetaspaceSize=1G
toolsVersion=0.13.3-kotlin-1.7.20
compose.version=1.2.1
compose.version=1.2.2

View File

@ -24,6 +24,7 @@ kotlin {
implementation(project(":numass-data-model"))
implementation("space.kscience:visionforge-core:$visionForgeVersion")
implementation("space.kscience:visionforge-plotly:$visionForgeVersion")
implementation(compose.runtime)
}
}
jvmMain {

View File

@ -2,11 +2,8 @@ package ru.inr.mass.data.server
import space.kscience.dataforge.context.Context
import space.kscience.dataforge.misc.DFExperimental
import space.kscience.visionforge.html.HtmlVisionFragment
import space.kscience.visionforge.html.ResourceLocation
import space.kscience.visionforge.html.VisionPage
import space.kscience.visionforge.html.importScriptHeader
import space.kscience.visionforge.makeFile
import space.kscience.visionforge.html.*
import space.kscience.visionforge.visionManager
import java.awt.Desktop
import java.nio.file.Path
@ -19,7 +16,7 @@ public fun Context.makeNumassWebFile(
show: Boolean = true,
content: HtmlVisionFragment,
): Unit {
val actualPath = VisionPage(this, content = content).makeFile(path) { actualPath: Path ->
val actualPath = VisionPage(visionManager, content = content).makeFile(path) { actualPath: Path ->
mapOf(
"title" to VisionPage.title(title),
"numassWeb" to VisionPage.importScriptHeader("js/numass-web.js", resourceLocation, actualPath)

View File

@ -32,8 +32,8 @@ public class NBkgSpectrum(public val source: Spectrum) : DifferentiableSpectrum
override fun derivativeOrNull(symbols: List<Symbol>): Spectrum? = when {
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) }
norm -> Spectrum { x, arguments -> source(x, arguments) }
bkg -> Spectrum { x, _ -> 1.0 }
else -> (source as? DifferentiableSpectrum)?.derivativeOrNull(symbols)?.let { NBkgSpectrum(it) }
}
else -> null

View File

@ -3,6 +3,8 @@
* To change this template file, choose Tools | Templates
* and open the template in the editor.
*/
@file:Suppress("PARAMETER_NAME_CHANGED_ON_OVERRIDE")
package ru.inr.mass.models
import space.kscience.kmath.expressions.Symbol

View File

@ -12,6 +12,8 @@ import kotlin.math.sqrt
/**
* @author [Alexander Nozik](mailto:altavir@gmail.com)
*/
@Suppress("PARAMETER_NAME_CHANGED_ON_OVERRIDE")
public class NumassResolution(
public val resA: Double = 8.3e-5,
public val resB: Double = 0.0,
@ -43,12 +45,12 @@ public class NumassResolution(
// else -> ResolutionFunction.getConstantTail()
// }
private fun getValueFast(E: Double, U: Double): Double {
val delta = resA * E
private fun getValueFast(e: Double, u: Double): Double {
val delta = resA * e
return when {
E - U < 0 -> 0.0
E - U > delta -> tailFunction(E, U)
else -> (E - U) / delta
e - u < 0 -> 0.0
e - u > delta -> tailFunction(e, u)
else -> (e - u) / delta
}
}
@ -57,15 +59,15 @@ public class NumassResolution(
override val x: Symbol get() = e
override val y: Symbol get() = u
override fun invoke(E: Double, U: Double, arguments: Map<Symbol, Double>): Double {
override fun invoke(e: Double, u: Double, arguments: Map<Symbol, Double>): Double {
if (resB <= 0) {
return this.getValueFast(E, U)
return this.getValueFast(e, u)
}
val delta = resA * E
val delta = resA * e
return when {
E - U < 0 -> 0.0
E - U > delta -> tailFunction(E, U)
else -> (1 - sqrt(1 - (E - U) / E * resB)) / (1 - sqrt(1 - resA * resB))
e - u < 0 -> 0.0
e - u > delta -> tailFunction(e, u)
else -> (1 - sqrt(1 - (e - u) / e * resB)) / (1 - sqrt(1 - resA * resB))
}
}

View File

@ -274,13 +274,13 @@ public class NumassTransmission(
* Значение полной производной функции потерь с учетом всех неисчезающих
* порядков
*
* @param X
* @param thickness
* @param eIn
* @param eOut
* @return
*/
private fun getTotalLossDeriv(X: Double, eIn: Double, eOut: Double): Double {
val probs = getLossProbDerivs(X)
private fun getTotalLossDeriv(thickness: Double, eIn: Double, eOut: Double): Double {
val probs = getLossProbDerivs(thickness)
var sum = 0.0
for (i in 1 until probs.size) {
@ -399,8 +399,8 @@ public class NumassTransmission(
// return getSingleScatterFunction(exPos, ionPos, exW, ionW, exIonRatio)
// }
public val trapFunction: (Double, Double) -> Double = { Ei: Double, Ef: Double ->
val eps = Ei - Ef
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 {

View File

@ -5,6 +5,7 @@ import ru.inr.mass.workspace.buffer
import ru.inr.mass.workspace.fitWith
import ru.inr.mass.workspace.generate
import space.kscience.kmath.expressions.Symbol
import space.kscience.kmath.expressions.derivative
import space.kscience.kmath.misc.UnstableKMathAPI
import space.kscience.kmath.operations.asSequence
import space.kscience.kmath.optimization.*
@ -28,10 +29,14 @@ suspend fun main() {
NumassBeta.e0 to 18575.0,
NumassBeta.msterile2 to 1000.0.pow(2),
NumassBeta.u2 to 1e-2,
NumassTransmission.thickness to 0.1,
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 }
@ -41,11 +46,12 @@ suspend fun main() {
val fit: XYFit = generatedData.fitWith(
optimizer = QowOptimizer,
modelExpression = spectrum,
startingPoint = args,
OptimizationParameters(NBkgSpectrum.norm, NBkgSpectrum.bkg, NumassBeta.e0)
startingPoint = args + mapOf(NBkgSpectrum.norm to 8.1e5),
OptimizationParameters(NBkgSpectrum.norm, NBkgSpectrum.bkg, NumassBeta.e0),
OptimizationIterations(20)
)
val fitResult = fit.resultPoint
println("Chi squared/dof: ${fit.chiSquaredOrNull}/${fit.dof}")
Plotly.plot {
scatter {
@ -66,7 +72,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, fitResult) }
y.numbers = x.doubles.map { spectrum(it, args + fit.resultPoint) }
}
}.makeFile()
}

View File

@ -11,8 +11,8 @@ import space.kscience.kmath.misc.FeatureSet
import space.kscience.kmath.misc.Loggable
import space.kscience.kmath.misc.UnstableKMathAPI
import space.kscience.kmath.optimization.*
import space.kscience.kmath.random.RandomGenerator
import space.kscience.kmath.samplers.PoissonSampler
import space.kscience.kmath.stat.RandomGenerator
import space.kscience.kmath.stat.next
import space.kscience.kmath.structures.asBuffer
import kotlin.math.sqrt