modified: lmfit_test.py

	modified:   src/fitter/common.rs
	modified:   src/fitter/fit_handler.rs
	modified:   src/fitter/fit_settings.rs
	modified:   src/fitter/main_fitter.rs
	modified:   src/fitter/mod.rs
	new file:   src/fitter/models/exponential.rs
	modified:   src/fitter/models/linear.rs
	modified:   src/fitter/models/mod.rs
	new file:   src/fitter/models/powerlaw.rs
	new file:   src/fitter/models/quadratic.rs
	modified:   src/histoer/histo1d/histogram1d.rs

lmfit_test.py

@@ -107,118 +107,184 @@ def MultipleGaussianFit(x_data: list, y_data: list, peak_markers: list,
 
     return gaussian_params, background_params, x_data_line, y_data_line, fit_report
 
-def MultipleGaussianFitWithLinearBG(x_data: list, y_data: list, mean: list, 
-                        equal_sigma: bool = True, free_position: bool = True,
-                        bg_initial_guesses: list = (0, 0), bg_vary: list = (True, True, True)):
+# Multiple Gaussian fitting function
+def LinearFit(x_data: list, y_data: list, slope: list = ("slope", -np.inf, np.inf, 0.0, True), intercept = ("intercept", -np.inf, np.inf, 0.0, True)):
+    import lmfit
+    import numpy as np
 
-    model = lmfit.models.LinearModel(prefix='bg_')
-    params = model.make_params(slope=bg_initial_guesses[0], intercept=bg_initial_guesses[1])
-    params['bg_slope'].set(vary=bg_vary[0])
-    params['bg_intercept'].set(vary=bg_vary[1])
-
-    first_gaussian = lmfit.Model(gaussian, prefix=f'g0_')
+    # params = slope=[name, min, max, initial_guess, vary], intercept=[name, min, max, initial_guess, vary]
 
-    if model is None:
-        model = first_gaussian
+    model = lmfit.models.LinearModel()
+    params = model.make_params(slope=slope[3], intercept=intercept[3])
+    params['slope'].set(min=slope[1], max=slope[2], value=slope[3], vary=slope[4])
+    params['intercept'].set(min=intercept[1], max=intercept[2], value=intercept[3], vary=intercept[4])
+
+    result = model.fit(y_data, params, x=x_data)
+
+    print(result.fit_report())
+
+    # Extract Parameters
+    slope = float(result.params['slope'].value)
+    slope_err = result.params['slope'].stderr
+
+    if slope_err is None:
+        slope_err = float(0.0)
     else:
-        model += first_gaussian
-
-    params.update(first_gaussian.make_params(amplitude=1, mean=peak_markers[0], sigma=1))
-    params['g0_sigma'].set(min=0)  # Initial constraint for the first Gaussian's sigma
+        slope_err = float(slope_err)
+
+    intercept = float(result.params['intercept'].value)
 
-    # Fix the center of the first Gaussian if free_position=False
-    if not free_position:
-        params['g0_mean'].set(vary=False)
+    intercept_err = result.params['intercept'].stderr
+    if intercept_err is None:
+        intercept_err = float(0.0)
+    else:
+        intercept_err = float(intercept_err)
 
-    # Add additional Gaussians
-    for i, peak in enumerate(peak_markers[1:], start=1):
-        g = lmfit.Model(gaussian, prefix=f'g{i}_')
-        model += g
-        params.update(g.make_params(amplitude=1, mean=peak))
+    print(f"Slope: {slope} ± {slope_err}")
+    print(f"Intercept: {intercept} ± {intercept_err}")
 
-        # Set the sigma parameter depending on whether equal_sigma is True or False
-        if equal_sigma:
-            params[f'g{i}_sigma'].set(expr='g0_sigma')  # Constrain sigma to be the same as g1_sigma
-        else:
-            params[f'g{i}_sigma'].set(min=0)  # Allow different sigmas for each Gaussian
 
-        params[f'g{i}_amplitude'].set(min=0)
+    params = [
+        ('slope', slope, slope_err),
+        ('intercept', intercept, intercept_err)
+    ]
 
-        # Fix the center of the Gaussian if free_position=False
-        if not free_position:
-            params[f'g{i}_mean'].set(vary=False)
+    x = np.linspace(x_data[0], x_data[-1], 5 * len(x_data))
+    y = result.eval(x=x)
 
-    # Fit the model to the data
+    fit_report = str(result.fit_report())
+
+    return params, x, y, fit_report
+
+def QuadraticFit(x_data: list, y_data: list, a: list = ("a", -np.inf, np.inf, 0.0, True), b = ("b", -np.inf, np.inf, 0.0, True), c: list = ("a", -np.inf, np.inf, 0.0, True),):    
+    # params = [name, min, max, initial_guess, vary]
+
+    model = lmfit.models.QuadraticModel()
+    params = model.make_params(a=a[3], b=b[3], c=c[3])
+    params['a'].set(min=a[1], max=a[2], value=a[3], vary=a[4])
+    params['b'].set(min=b[1], max=b[2], value=b[3], vary=b[4])
+    params['c'].set(min=c[1], max=c[2], value=c[3], vary=c[4])
     result = model.fit(y_data, params, x=x_data)
 
-    print("\nFit Report:")
     print(result.fit_report())
 
-    # Create a list of native Python float tuples for Gaussian parameters
-    gaussian_params = []
-    for i in range(len(peak_markers)):
-        gaussian_params.append((
-            float(result.params[f'g{i}_amplitude'].value),
-            float(result.params[f'g{i}_amplitude'].stderr),
-            float(result.params[f'g{i}_mean'].value),
-            float(result.params[f'g{i}_mean'].stderr),
-            float(result.params[f'g{i}_sigma'].value),
-            float(result.params[f'g{i}_sigma'].stderr)
-        ))
+    # Extract Parameters
+    a = float(result.params['a'].value)
+    a_err = result.params['a'].stderr
+    if a_err is None:
+        a_err = float(0.0)
+    else:
+        a_err = float(a_err)
 
-    # Create a list of native Python float tuples for Background parameters
-    background_params = []
-    if bg_type != 'None':
-        for key in result.params:
-            if 'bg_' in key:
-                background_params.append((
-                    key,  # Keep the parameter name
-                    float(result.params[key].value),      # Convert the value to native Python float
-                    float(result.params[key].stderr)      # Convert the uncertainty to native Python float
-                ))
+    b = float(result.params['b'].value)
+    b_err = result.params['b'].stderr
+    if b_err is None:
+        b_err = float(0.0)
+    else:
+        b_err = float(b_err)
 
-    # Create smooth fit line with plenty of data points
-    x_data_line = np.linspace(x_data[0], x_data[-1], 5 * len(x_data))
-    y_data_line = result.eval(x=x_data_line)
+    c = float(result.params['c'].value)
+    c_err = result.params['c'].stderr
+    if c_err is None:
+        c_err = float(0.0)
+    else:
+        c_err = float(c_err)
+
+
+    params = [
+        ('a', a, a_err),
+        ('b', b, b_err),
+        ('c', c, c_err)
+    ]
+
+    x = np.linspace(x_data[0], x_data[-1], 5 * len(x_data))
+    y = result.eval(x=x)
 
     fit_report = str(result.fit_report())
 
-    return gaussian_params, background_params, x_data_line, y_data_line, fit_report
+    return params, x, y, fit_report
 
-# Multiple Gaussian fitting function
-def LinearFit(x_data: list, y_data: list, initial_guess: tuple = (0.0, 0.0), vary: tuple = (True, True)):
-    import lmfit
-    import numpy as np
+def PowerLawFit(x_data: list, y_data: list, amplitude: list = ("amplitude", -np.inf, np.inf, 0.0, True), exponent = ("exponent", -np.inf, np.inf, 0.0, True)):    
+    # params = [name, min, max, initial_guess, vary]
 
+    model = lmfit.models.PowerLawModel()
+    params = model.make_params(amplitude=amplitude[3], exponent=exponent[3])
+    params['amplitude'].set(min=amplitude[1], max=amplitude[2], value=amplitude[3], vary=amplitude[4])
+    params['exponent'].set(min=exponent[1], max=exponent[2], value=exponent[3], vary=exponent[4])
+    result = model.fit(y_data, params, x=x_data)
+
+    print(result.fit_report())
+
+    # Extract Parameters
+    amplitude = float(result.params['amplitude'].value)
+    amplitude_err = result.params['amplitude'].stderr
+    if amplitude_err is None:
+        amplitude_err = float(0.0)
+    else:
+        amplitude_err = float(amplitude_err)
 
-    model = lmfit.models.LinearModel()
-    params = model.make_params(slope=initial_guess[0], intercept=initial_guess[1])
-    params['slope'].set(vary=vary[0])
-    params['intercept'].set(vary=vary[1])
+    exponent = float(result.params['exponent'].value)
+    exponent_err = result.params['exponent'].stderr
+    if exponent_err is None:
+        exponent_err = float(0.0)
+    else:
+        exponent_err = float(exponent_err)
+
+    params = [
+        ('amplitude', amplitude, amplitude_err),
+        ('exponent', exponent, exponent_err)
+    ]
+
+    x = np.linspace(x_data[0], x_data[-1], 5 * len(x_data))
+    y = result.eval(x=x)
+
+    fit_report = str(result.fit_report())
+
+    return params, x, y, fit_report
 
+def ExponentialFit(x_data: list, y_data: list, amplitude: list = ("amplitude", -np.inf, np.inf, 0.0, True), decay = ("decay", -np.inf, np.inf, 0.0, True)):    
+    # params = [name, min, max, initial_guess, vary]
+
+    model = lmfit.models.ExponentialModel()
+    params = model.make_params(amplitude=amplitude[3], decay=decay[3])
+    params['amplitude'].set(min=amplitude[1], max=amplitude[2], value=amplitude[3], vary=amplitude[4])
+    params['decay'].set(min=decay[1], max=decay[2], value=decay[3], vary=decay[4])
     result = model.fit(y_data, params, x=x_data)
 
     print(result.fit_report())
 
     # Extract Parameters
-    slope = float(result.params['slope'].value)
-    slope_err = float(result.params['slope'].stderr)
-    intercept = float(result.params['intercept'].value)
-    intercept_err = float(result.params['intercept'].stderr)
+    amplitude = float(result.params['amplitude'].value)
+    amplitude_err = result.params['amplitude'].stderr
+    if amplitude_err is None:
+        amplitude_err = float(0.0)
+    else:
+        amplitude_err = float(amplitude_err)
+
+    decay = float(result.params['decay'].value)
+    decay_err = result.params['decay'].stderr
+    if decay_err is None:
+        decay_err = float(0.0)
+    else:
+        decay_err = float(decay_err)
 
     params = [
-        ('slope', slope, slope_err),
-        ('intercept', intercept, intercept_err)
+        ('amplitude', amplitude, amplitude_err),
+        ('decay', decay, decay_err)
     ]
 
     x = np.linspace(x_data[0], x_data[-1], 5 * len(x_data))
     y = result.eval(x=x)
 
-    return params, x, y
+    fit_report = str(result.fit_report())
+
+    return params, x, y, fit_report
+
 
 # Load the data using Polars
 df = pl.read_parquet("/Users/alconley/Projects/ICESPICE/207Bi/exp_data/207Bi_noICESPICE_f9mm_g0mm_run_13.parquet")
 
+
 df = df.with_columns([
     (pl.col("PIPS1000Energy") * 0.5395 + 2.5229).alias("PIPS1000EnergyCalibrated")
 ])
@@ -240,11 +306,15 @@ def LinearFit(x_data: list, y_data: list, initial_guess: tuple = (0.0, 0.0), var
 peak_markers = [554, 567]  # Replace with actual peak guesses
 
 # Fit the data
-gaussian_params, background_params, x_data_line, y_data_line = MultipleGaussianFit(bin_centers_filtered, counts_filtered, peak_markers, equal_sigma=True, free_position=False, bg_type='linear')
+# gaussian_params, background_params, x_data_line, y_data_line = MultipleGaussianFit(bin_centers_filtered, counts_filtered, peak_markers, equal_sigma=True, free_position=False, bg_type='linear')
+
+# slope = ("slope", -np.inf, np.inf, -2.0, False)
+# intercept = ("intercept", -np.inf, np.inf, 0.0, False)
+# LinearFit(bin_centers_filtered, counts_filtered, slope=slope, intercept=intercept)
 
-LinearFit(bin_centers_filtered, counts_filtered, initial_guess=(1000.0, 0.0), vary=(True, True))
+ExponentialFit(bin_centers_filtered, counts_filtered)
 
-# 
+# # 
 # # Plot the original data and the fit
 # plt.figure(figsize=(8, 6))
 

src/fitter/common.rs

@@ -1,19 +1,68 @@
-#[derive(Debug, Clone, serde::Deserialize, serde::Serialize)]
+#[derive(PartialEq, Debug, Clone, serde::Deserialize, serde::Serialize)]
 pub struct Data {
     pub x: Vec<f64>,
     pub y: Vec<f64>,
 }
 
-#[derive(Debug, Clone, serde::Deserialize, serde::Serialize)]
-pub struct Value {
+#[derive(PartialEq, Debug, Clone, serde::Deserialize, serde::Serialize)]
+pub struct Parameter {
     pub name: String,
-    pub value: f64,
-    pub uncertainity: f64,
+    pub min: f64,
+    pub max: f64,
+    pub initial_guess: f64,
+    pub vary: bool,
+    pub value: Option<f64>,
+    pub uncertainity: Option<f64>,
 }
 
-impl Value {
-    pub fn new(name: String, value: f64, uncertainity: f64) -> Self {
-        Value { name, value, uncertainity }
+impl Default for Parameter {
+    fn default() -> Self {
+        Parameter {
+            name: String::new(),
+            min: f64::NEG_INFINITY,
+            max: f64::INFINITY,
+            initial_guess: 0.0,
+            vary: true,
+            value: None,
+            uncertainity: None,
+        }
     }
-
-}
+}
+
+impl Parameter {
+    pub fn ui(&mut self, ui: &mut egui::Ui) {
+        ui.label(&self.name);
+        ui.add(
+            egui::DragValue::new(&mut self.initial_guess).speed(0.1), // .prefix("Initial Guess: ")
+                                                                      // .suffix(" a.u."),
+        )
+        .on_hover_text(format!("Initial guess for the {} parameter", self.name));
+
+        ui.add(
+            egui::DragValue::new(&mut self.min)
+                .speed(0.1)
+                // .prefix("Min: ")
+                .range(f64::NEG_INFINITY..=self.max), // .suffix(" a.u."),
+        )
+        .on_hover_text(format!("Minimum value for the {} parameter", self.name));
+
+        ui.add(
+            egui::DragValue::new(&mut self.max)
+                .speed(0.1)
+                // .prefix("Max: ")
+                .range(self.min..=f64::INFINITY), // .suffix(" a.u."),
+        )
+        .on_hover_text(format!("Maximum value for the {} parameter", self.name));
+
+        ui.checkbox(&mut self.vary, "").on_hover_text(format!(
+            "Allow the {} parameter to vary during the fitting process",
+            self.name
+        ));
+
+        if let Some(value) = self.value {
+            ui.separator();
+            ui.label(format!("{:.3}", value));
+            ui.label(format!("{:.3}", self.uncertainity.unwrap_or(0.0)));
+        }
+    }
+}