From b9015d8d4173e1bf6fbda9e2ce5a23fb0ba094c1 Mon Sep 17 00:00:00 2001 From: alconley Date: Mon, 21 Oct 2024 21:53:48 -0400 Subject: [PATCH] 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 | 224 +++++++++++++++++--------- src/fitter/common.rs | 69 ++++++-- src/fitter/fit_handler.rs | 26 +-- src/fitter/fit_settings.rs | 66 +++++++- src/fitter/main_fitter.rs | 212 +++++++++++++----------- src/fitter/mod.rs | 2 +- src/fitter/models/exponential.rs | 217 +++++++++++++++++++++++++ src/fitter/models/linear.rs | 185 ++++++++++++++------- src/fitter/models/mod.rs | 7 +- src/fitter/models/powerlaw.rs | 217 +++++++++++++++++++++++++ src/fitter/models/quadratic.rs | 248 +++++++++++++++++++++++++++++ src/histoer/histo1d/histogram1d.rs | 113 +++++++------ 12 files changed, 1264 insertions(+), 322 deletions(-) create mode 100644 src/fitter/models/exponential.rs create mode 100644 src/fitter/models/powerlaw.rs create mode 100644 src/fitter/models/quadratic.rs diff --git a/lmfit_test.py b/lmfit_test.py index aa803c7..a314417 100644 --- a/lmfit_test.py +++ b/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)) diff --git a/src/fitter/common.rs b/src/fitter/common.rs index 220e84b..cc31525 100644 --- a/src/fitter/common.rs +++ b/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, pub y: Vec, } -#[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, + pub uncertainity: Option, } -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, + } } - -} \ No newline at end of file +} + +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))); + } + } +} diff --git a/src/fitter/fit_handler.rs b/src/fitter/fit_handler.rs index c551cb1..718584a 100644 --- a/src/fitter/fit_handler.rs +++ b/src/fitter/fit_handler.rs @@ -39,7 +39,6 @@ impl Fits { self.stored_fits.push(temp_fit.clone()); } - } pub fn set_log(&mut self, log_y: bool, log_x: bool) { @@ -55,7 +54,6 @@ impl Fits { pub fn set_stored_fits_background_color(&mut self, color: egui::Color32) { for fit in &mut self.stored_fits { fit.background_line.color = color; - } } @@ -218,22 +216,6 @@ impl Fits { } } - pub fn fit_lines_ui(&mut self, ui: &mut egui::Ui) { - egui::ScrollArea::vertical().show(ui, |ui| { - ui.vertical(|ui| { - ui.label("Fit Lines"); - - if let Some(temp_fit) = &mut self.temp_fit { - temp_fit.lines_ui(ui); - } - - for fit in &mut self.stored_fits { - fit.lines_ui(ui); - } - }); - }); - } - pub fn fit_context_menu_ui(&mut self, ui: &mut egui::Ui) { ui.menu_button("Fits", |ui| { self.save_and_load_ui(ui); @@ -251,11 +233,9 @@ impl Fits { ui.separator(); - egui::ScrollArea::vertical() - .max_height(300.0) - .show(ui, |ui| { - self.fit_lines_ui(ui); - }); + if let Some(temp_fit) = &mut self.temp_fit { + temp_fit.fit_result_ui(ui); + } }); } } diff --git a/src/fitter/fit_settings.rs b/src/fitter/fit_settings.rs index 4f18b78..f0e029d 100644 --- a/src/fitter/fit_settings.rs +++ b/src/fitter/fit_settings.rs @@ -1,4 +1,8 @@ -use super::main_fitter::Model; +use crate::fitter::main_fitter::BackgroundModel; +use crate::fitter::models::exponential::ExponentialParameters; +use crate::fitter::models::linear::LinearParameters; +use crate::fitter::models::powerlaw::PowerLawParameters; +use crate::fitter::models::quadratic::QuadraticParameters; #[derive(Debug, Clone, serde::Deserialize, serde::Serialize)] pub struct FitSettings { @@ -9,7 +13,11 @@ pub struct FitSettings { pub fit_stats_height: f32, pub equal_stddev: bool, pub free_position: bool, - pub background_model: Model, + pub background_model: BackgroundModel, + pub linear_params: LinearParameters, + pub quadratic_params: QuadraticParameters, + pub power_law_params: PowerLawParameters, + pub exponential_params: ExponentialParameters, } impl Default for FitSettings { @@ -22,7 +30,11 @@ impl Default for FitSettings { fit_stats_height: 0.0, equal_stddev: true, free_position: true, - background_model: Model::Linear, + background_model: BackgroundModel::Linear(LinearParameters::default()), + linear_params: LinearParameters::default(), + quadratic_params: QuadraticParameters::default(), + power_law_params: PowerLawParameters::default(), + exponential_params: ExponentialParameters::default(), } } } @@ -68,11 +80,51 @@ impl FitSettings { ui.separator(); - ui.heading("Background Fit Models"); + ui.heading("Background Models"); ui.horizontal(|ui| { - ui.radio_value(&mut self.background_model, Model::Linear, "Linear"); - ui.radio_value(&mut self.background_model, Model::None, "None"); + ui.radio_value( + &mut self.background_model, + BackgroundModel::Linear(self.linear_params.clone()), + "Linear", + ); + ui.radio_value( + &mut self.background_model, + BackgroundModel::Quadratic(self.quadratic_params.clone()), + "Quadratic", + ); + ui.radio_value( + &mut self.background_model, + BackgroundModel::PowerLaw(self.power_law_params.clone()), + "Power Law", + ); + ui.radio_value( + &mut self.background_model, + BackgroundModel::Exponential(self.exponential_params.clone()), + "Exponential", + ); + ui.radio_value(&mut self.background_model, BackgroundModel::None, "None"); }); - + + if let BackgroundModel::Linear(params) = &mut self.background_model { + params.ui(ui); + self.linear_params = params.clone(); + } + + if let BackgroundModel::Quadratic(params) = &mut self.background_model { + params.ui(ui); + self.quadratic_params = params.clone(); + } + + if let BackgroundModel::PowerLaw(params) = &mut self.background_model { + params.ui(ui); + self.power_law_params = params.clone(); + } + + if let BackgroundModel::Exponential(params) = &mut self.background_model { + params.ui(ui); + self.exponential_params = params.clone(); + } + + ui.separator(); } } diff --git a/src/fitter/main_fitter.rs b/src/fitter/main_fitter.rs index b378ca8..8866837 100644 --- a/src/fitter/main_fitter.rs +++ b/src/fitter/main_fitter.rs @@ -1,29 +1,51 @@ use crate::egui_plot_stuff::egui_line::EguiLine; -use super::models::gaussian::{Background, GaussianFitter}; -use super::common::{Data, Value}; -use super::models::linear::LinearFitter; +// use super::models::gaussian::{Background, GaussianFitter}; +use super::common::Data; +use super::models::exponential::{ExponentialFitter, ExponentialParameters}; +use super::models::linear::{LinearFitter, LinearParameters}; +use super::models::powerlaw::{PowerLawFitter, PowerLawParameters}; +use super::models::quadratic::{QuadraticFitter, QuadraticParameters}; #[derive(Debug, Clone, serde::Deserialize, serde::Serialize, PartialEq)] -pub enum Model { +pub enum FitModel { Gaussian(Vec, bool, bool, f64), // initial peak locations, free sigma, free position, bin width - Linear, // slope, intercept, vary parameters None, } #[derive(Debug, Clone, serde::Deserialize, serde::Serialize)] pub enum FitResult { - Gaussian(GaussianFitter), + // Gaussian(GaussianFitter), +} + +#[derive(PartialEq, Debug, Clone, serde::Deserialize, serde::Serialize)] +pub enum BackgroundModel { + Linear(LinearParameters), + Quadratic(QuadraticParameters), + PowerLaw(PowerLawParameters), + Exponential(ExponentialParameters), + None, +} + +#[derive(Debug, Clone, serde::Deserialize, serde::Serialize)] +pub enum BackgroundResult { Linear(LinearFitter), + Quadratic(QuadraticFitter), + PowerLaw(PowerLawFitter), + Exponential(ExponentialFitter), } #[derive(Debug, Clone, serde::Deserialize, serde::Serialize)] pub struct Fitter { pub name: String, + pub data: Data, - pub model: Model, - pub background: Model, - pub background_result: Option, - pub result: Option, + + pub background_model: BackgroundModel, + pub background_result: Option, + + pub fit_model: FitModel, + pub fit_result: Option, + pub background_line: EguiLine, pub composition_line: EguiLine, pub decomposition_lines: Vec, @@ -31,117 +53,99 @@ pub struct Fitter { impl Fitter { // Constructor to create a new Fitter with empty data and specified model - pub fn new(model: Model, data: Data) -> Self { + pub fn new(data: Data) -> Self { Fitter { name: "Fit".to_string(), + data, - model, - background: Model::Linear, + + background_model: BackgroundModel::None, background_result: None, - result: None, + + fit_model: FitModel::None, + fit_result: None, + background_line: EguiLine::new(egui::Color32::GREEN), composition_line: EguiLine::new(egui::Color32::DARK_BLUE), decomposition_lines: Vec::new(), } } - pub fn fit(&mut self) { - // Perform the fit based on the model - match &self.model { - Model::Gaussian(peak_markers, equal_sigma, free_position, bin_width) => { - // Extract background parameters if the background model is fitted - let background_parameters = match &self.background_result { - Some(FitResult::Linear(fit)) => { - // Extract background parameters and fit points from LinearFitter - Some(Background { - model: Model::Linear, - parameters: Some(vec![ - fit.slope.clone().unwrap(), - fit.intercept.clone().unwrap(), - ]), - varying: Some(vec![fit.vary_parameters.0, fit.vary_parameters.1]), - fit_points: fit.fit_points.clone(), - }) - } - Some(FitResult::Gaussian(_)) => { - None // Background is Gaussian, not relevant here + pub fn fit_background(&mut self) { + log::info!("Fitting background"); + match &self.background_model { + BackgroundModel::Linear(params) => { + let mut fit = LinearFitter::new(self.data.clone()); + + // Copy the parameters from the background model to the LinearFitter + fit.paramaters = params.clone(); // Ensure LinearParameters are used + + // Perform the fit + match fit.lmfit() { + Ok(_) => { + self.background_line.points = fit.fit_points.clone(); + self.background_result = Some(BackgroundResult::Linear(fit)); } - None => { - None // No background result present + Err(e) => { + eprintln!("Error: {}", e); } - }; - - // Create the GaussianFitter instance with the pre-existing background parameters - let mut fit = GaussianFitter::new( - self.data.clone(), - peak_markers.clone(), - self.background.clone(), - background_parameters, // Pass background parameters if they exist - *equal_sigma, - *free_position, - *bin_width, - ); - - // Perform the fit using lmfit + } + } + BackgroundModel::Quadratic(params) => { + let mut fit = QuadraticFitter::new(self.data.clone()); + + // Copy the parameters from the background model to the QuadraticFitter + fit.paramaters = params.clone(); // Ensure QuadraticParameters are used + + // Perform the fit match fit.lmfit() { Ok(_) => { - // Process Gaussian fit result and update lines - self.composition_line.points = fit.composition_points.clone(); - self.result = Some(FitResult::Gaussian(fit)); + self.background_line.points = fit.fit_points.clone(); + self.background_result = Some(BackgroundResult::Quadratic(fit)); } Err(e) => { eprintln!("Error: {}", e); } } } - Model::Linear => { - // Perform Linear fit - let mut fit = LinearFitter::new(self.data.clone()); + BackgroundModel::PowerLaw(params) => { + let mut fit = PowerLawFitter::new(self.data.clone()); + + // Copy the parameters from the background model to the PowerLawFitter + fit.paramaters = params.clone(); // Ensure PowerLawParameters are used + + // Perform the fit match fit.lmfit() { Ok(_) => { - self.composition_line.points = fit.fit_points.clone(); - self.result = Some(FitResult::Linear(fit)); + self.background_line.points = fit.fit_points.clone(); + self.background_result = Some(BackgroundResult::PowerLaw(fit)); } Err(e) => { eprintln!("Error: {}", e); } } } - Model::None => { - // No fitting required for 'None' - } - } - } + BackgroundModel::Exponential(params) => { + let mut fit = ExponentialFitter::new(self.data.clone()); + // Copy the parameters from the background model to the ExponentialFitter + fit.paramaters = params.clone(); // Ensure ExponentialParameters are used - pub fn fit_background(&mut self) { - log::info!("Fitting background"); - // Perform the fit based on the background model - match &self.background { - Model::Gaussian(_, _, _, _) => { - // Add Gaussian background fit logic here if needed - } - Model::Linear => { - log::info!("Fitting background with a linear line using `lmfit`."); - // Perform Linear fit for background - let mut fit = LinearFitter::new(self.data.clone()); - log::info!("{:?}", fit); + // Perform the fit match fit.lmfit() { Ok(_) => { self.background_line.points = fit.fit_points.clone(); - self.background_result = Some(FitResult::Linear(fit)); + self.background_result = Some(BackgroundResult::Exponential(fit)); } Err(e) => { eprintln!("Error: {}", e); } } } - Model::None => { - // No background fitting required for 'None' + BackgroundModel::None => { log::info!("No background fitting required for 'None'"); } } - log::info!("Finished fitting background"); } @@ -175,24 +179,52 @@ impl Fitter { pub fn set_name(&mut self, name: String) { self.composition_line.name = format!("{}-Composition", name); - + for (i, line) in self.decomposition_lines.iter_mut().enumerate() { line.name = format!("{}-Peak {}", name, i); } - self.background_line.name = format!("{}-Background", name); + self.background_line.name = format!("{}-Background", name); } - pub fn lines_ui(&mut self, ui: &mut egui::Ui) { - self.background_line.menu_button(ui); - - self.composition_line.menu_button(ui); + pub fn fit_result_ui(&mut self, ui: &mut egui::Ui) { + ui.collapsing(self.name.clone(), |ui| { + egui::ScrollArea::vertical() + .min_scrolled_height(300.0) + .show(ui, |ui| { + if let Some(fit_result) = &self.fit_result { + self.composition_line.menu_button(ui); + } - for line in &mut self.decomposition_lines { - line.menu_button(ui); - } + for line in &mut self.decomposition_lines { + line.menu_button(ui); + } - ui.separator(); + ui.separator(); + + if let Some(background_result) = &self.background_result { + ui.label("Background"); + match background_result { + BackgroundResult::Linear(fit) => { + fit.ui(ui); + } + BackgroundResult::Quadratic(fit) => { + fit.ui(ui); + } + BackgroundResult::PowerLaw(fit) => { + fit.ui(ui); + } + BackgroundResult::Exponential(fit) => { + fit.ui(ui); + } + } + ui.horizontal(|ui| { + ui.label("Line"); + self.background_line.menu_button(ui); + }); + } + }); + }); } // Draw the background, decomposition, and composition lines @@ -219,4 +251,4 @@ impl Fitter { self.background_line.log_y = log_y; self.background_line.log_x = log_x; } -} \ No newline at end of file +} diff --git a/src/fitter/mod.rs b/src/fitter/mod.rs index ea42c5a..64ad653 100644 --- a/src/fitter/mod.rs +++ b/src/fitter/mod.rs @@ -1,5 +1,5 @@ +pub mod common; pub mod fit_handler; pub mod fit_settings; pub mod main_fitter; pub mod models; -pub mod common; \ No newline at end of file diff --git a/src/fitter/models/exponential.rs b/src/fitter/models/exponential.rs new file mode 100644 index 0000000..c528748 --- /dev/null +++ b/src/fitter/models/exponential.rs @@ -0,0 +1,217 @@ +use crate::fitter::common::{Data, Parameter}; +use pyo3::{prelude::*, types::PyModule}; + +#[derive(PartialEq, Debug, Clone, serde::Deserialize, serde::Serialize)] +pub struct ExponentialParameters { + pub amplitude: Parameter, + pub decay: Parameter, +} + +impl Default for ExponentialParameters { + fn default() -> Self { + ExponentialParameters { + amplitude: Parameter { + name: "amplitude".to_string(), + ..Default::default() + }, + decay: Parameter { + name: "decay".to_string(), + ..Default::default() + }, + } + } +} + +impl ExponentialParameters { + pub fn ui(&mut self, ui: &mut egui::Ui) { + ui.horizontal(|ui| { + ui.label("Fit Parameters"); + if ui.small_button("Reset").clicked() { + *self = ExponentialParameters::default(); + } + }); + // create a grid for the param + egui::Grid::new("Exponential_params_grid") + .striped(true) + .num_columns(5) + .show(ui, |ui| { + ui.label("Parameter"); + ui.label("Initial Guess"); + ui.label("Min"); + ui.label("Max"); + ui.label("Vary"); + ui.end_row(); + self.amplitude.ui(ui); + ui.end_row(); + self.decay.ui(ui); + }); + } +} + +#[derive(Debug, Clone, serde::Deserialize, serde::Serialize)] +pub struct ExponentialFitter { + pub data: Data, + pub paramaters: ExponentialParameters, + pub fit_points: Vec<[f64; 2]>, + pub fit_report: String, +} + +impl ExponentialFitter { + pub fn new(data: Data) -> Self { + ExponentialFitter { + data, + paramaters: ExponentialParameters::default(), + fit_points: Vec::new(), + fit_report: String::new(), + } + } + + pub fn lmfit(&mut self) -> PyResult<()> { + log::info!("Fitting data with a Exponential line using `lmfit`."); + Python::with_gil(|py| { + // let sys = py.import_bound("sys")?; + // let version: String = sys.getattr("version")?.extract()?; + // let executable: String = sys.getattr("executable")?.extract()?; + // println!("Using Python version: {}", version); + // println!("Python executable: {}", executable); + + // Check if the `uproot` module can be imported + match py.import_bound("lmfit") { + Ok(_) => { + // println!("Successfully imported `lmfit` module."); + } + Err(_) => { + eprintln!("Error: `lmfit` module could not be found. Make sure you are using the correct Python environment with `lmfit` installed."); + return Err(PyErr::new::( + "`lmfit` module not available", + )); + } + } + + match py.import_bound("numpy") { + Ok(_) => { + // println!("Successfully imported `lmfit` module."); + } + Err(_) => { + eprintln!("Error: `numpy` module could not be found. Make sure you are using the correct Python environment with `numpy` installed."); + return Err(PyErr::new::( + "`numpy` module not available", + )); + } + } + + // Define the Python code as a module + let code = r#" +import lmfit +import numpy as np + +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 + 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 = [ + ('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) + + fit_report = str(result.fit_report()) + + return params, x, y, fit_report +"#; + + // Compile the Python code into a module + let module = PyModule::from_code_bound(py, code, "Exponential.py", "Exponential")?; + + let x_data = self.data.x.clone(); + let y_data = self.data.y.clone(); + let amplitude_para = ( + self.paramaters.amplitude.name.clone(), + self.paramaters.amplitude.min, + self.paramaters.amplitude.max, + self.paramaters.amplitude.initial_guess, + self.paramaters.amplitude.vary, + ); + let decay_para = ( + self.paramaters.decay.name.clone(), + self.paramaters.decay.min, + self.paramaters.decay.max, + self.paramaters.decay.initial_guess, + self.paramaters.decay.vary, + ); + + let result = module.getattr("ExponentialFit")?.call1(( + x_data, + y_data, + amplitude_para, + decay_para, + ))?; + + let params = result.get_item(0)?.extract::>()?; + let x = result.get_item(1)?.extract::>()?; + let y = result.get_item(2)?.extract::>()?; + let fit_report = result.get_item(3)?.extract::()?; + + self.paramaters.amplitude.value = Some(params[0].1); + self.paramaters.amplitude.uncertainity = Some(params[0].2); + self.paramaters.decay.value = Some(params[1].1); + self.paramaters.decay.uncertainity = Some(params[1].2); + + self.fit_points = x.iter().zip(y.iter()).map(|(&x, &y)| [x, y]).collect(); + self.fit_report = fit_report; + + Ok(()) + }) + } + + pub fn ui(&self, ui: &mut egui::Ui) { + // add menu button for the fit report + ui.horizontal(|ui| { + if let Some(amplitude) = &self.paramaters.amplitude.value { + ui.label(format!( + "amplitude: {:.3} ± {:.3}", + amplitude, + self.paramaters.amplitude.uncertainity.unwrap_or(0.0) + )); + } + ui.separator(); + if let Some(decay) = &self.paramaters.decay.value { + ui.label(format!( + "decay: {:.3} ± {:.3}", + decay, + self.paramaters.decay.uncertainity.unwrap_or(0.0) + )); + } + ui.separator(); + ui.menu_button("Fit Report", |ui| { + ui.horizontal_wrapped(|ui| { + ui.label(self.fit_report.clone()); + }); + }); + }); + } +} diff --git a/src/fitter/models/linear.rs b/src/fitter/models/linear.rs index 93a8f8f..b75d12f 100644 --- a/src/fitter/models/linear.rs +++ b/src/fitter/models/linear.rs @@ -1,14 +1,57 @@ +use crate::fitter::common::{Data, Parameter}; use pyo3::{prelude::*, types::PyModule}; -use crate::fitter::common::{Data, Value}; +#[derive(PartialEq, Debug, Clone, serde::Deserialize, serde::Serialize)] +pub struct LinearParameters { + pub slope: Parameter, + pub intercept: Parameter, +} + +impl Default for LinearParameters { + fn default() -> Self { + LinearParameters { + slope: Parameter { + name: "slope".to_string(), + ..Default::default() + }, + intercept: Parameter { + name: "intercept".to_string(), + ..Default::default() + }, + } + } +} + +impl LinearParameters { + pub fn ui(&mut self, ui: &mut egui::Ui) { + ui.horizontal(|ui| { + ui.label("Fit Parameters"); + if ui.small_button("Reset").clicked() { + *self = LinearParameters::default(); + } + }); + // create a grid for the param + egui::Grid::new("linear_params_grid") + .striped(true) + .num_columns(5) + .show(ui, |ui| { + ui.label("Parameter"); + ui.label("Initial Guess"); + ui.label("Min"); + ui.label("Max"); + ui.label("Vary"); + ui.end_row(); + self.slope.ui(ui); + ui.end_row(); + self.intercept.ui(ui); + }); + } +} #[derive(Debug, Clone, serde::Deserialize, serde::Serialize)] pub struct LinearFitter { pub data: Data, - pub initial_guess: (f64, f64), - pub vary_parameters: (bool, bool), - pub slope: Option, - pub intercept: Option, + pub paramaters: LinearParameters, pub fit_points: Vec<[f64; 2]>, pub fit_report: String, } @@ -17,10 +60,7 @@ impl LinearFitter { pub fn new(data: Data) -> Self { LinearFitter { data, - initial_guess: (0.0, 0.0), - vary_parameters: (true, true), - slope: None, - intercept: None, + paramaters: LinearParameters::default(), fit_points: Vec::new(), fit_report: String::new(), } @@ -48,17 +88,32 @@ impl LinearFitter { } } + match py.import_bound("numpy") { + Ok(_) => { + // println!("Successfully imported `lmfit` module."); + } + Err(_) => { + eprintln!("Error: `numpy` module could not be found. Make sure you are using the correct Python environment with `numpy` installed."); + return Err(PyErr::new::( + "`numpy` module not available", + )); + } + } + // Define the Python code as a module let code = r#" +import lmfit +import numpy as np -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 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)): + # params + # slope=[name, min, max, initial_guess, vary] + # intercept = [name, min, max, initial_guess, vary] 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]) + 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) @@ -66,9 +121,20 @@ def LinearFit(x_data: list, y_data: list, initial_guess: tuple = (0.0, 0.0), var # Extract Parameters slope = float(result.params['slope'].value) - slope_err = float(result.params['slope'].stderr) + slope_err = result.params['slope'].stderr + + if slope_err is None: + slope_err = float(0.0) + else: + slope_err = float(slope_err) + intercept = float(result.params['intercept'].value) - intercept_err = float(result.params['intercept'].stderr) + + intercept_err = result.params['intercept'].stderr + if intercept_err is None: + intercept_err = float(0.0) + else: + intercept_err = float(intercept_err) params = [ ('slope', slope, slope_err), @@ -80,29 +146,44 @@ def LinearFit(x_data: list, y_data: list, initial_guess: tuple = (0.0, 0.0), var fit_report = str(result.fit_report()) - print(fit_report) - return params, x, y, fit_report "#; // Compile the Python code into a module - let module = - PyModule::from_code_bound(py, code, "linear.py", "linear")?; + let module = PyModule::from_code_bound(py, code, "linear.py", "linear")?; let x_data = self.data.x.clone(); let y_data = self.data.y.clone(); - let initial_guess = self.initial_guess; - let vary = self.vary_parameters; - - let result = module.getattr("LinearFit")?.call1((x_data, y_data, initial_guess, vary))?; + let slope_para = ( + self.paramaters.slope.name.clone(), + self.paramaters.slope.min, + self.paramaters.slope.max, + self.paramaters.slope.initial_guess, + self.paramaters.slope.vary, + ); + let intercept_para = ( + self.paramaters.intercept.name.clone(), + self.paramaters.intercept.min, + self.paramaters.intercept.max, + self.paramaters.intercept.initial_guess, + self.paramaters.intercept.vary, + ); + + let result = + module + .getattr("LinearFit")? + .call1((x_data, y_data, slope_para, intercept_para))?; let params = result.get_item(0)?.extract::>()?; let x = result.get_item(1)?.extract::>()?; let y = result.get_item(2)?.extract::>()?; let fit_report = result.get_item(3)?.extract::()?; - self.slope = Some(Value::new(params[0].0.clone(), params[0].1, params[0].2)); - self.intercept = Some(Value::new(params[1].0.clone(), params[1].1, params[1].2)); + self.paramaters.slope.value = Some(params[0].1); + self.paramaters.slope.uncertainity = Some(params[0].2); + self.paramaters.intercept.value = Some(params[1].1); + self.paramaters.intercept.uncertainity = Some(params[1].2); + self.fit_points = x.iter().zip(y.iter()).map(|(&x, &y)| [x, y]).collect(); self.fit_report = fit_report; @@ -110,36 +191,30 @@ def LinearFit(x_data: list, y_data: list, initial_guess: tuple = (0.0, 0.0), var }) } - pub fn ui(&mut self, ui: &mut egui::Ui) { + pub fn ui(&self, ui: &mut egui::Ui) { + // add menu button for the fit report ui.horizontal(|ui| { - ui.label("Initial Guess:"); - ui.add(egui::DragValue::new(&mut self.initial_guess.0).speed(0.1).prefix("Slope: ")); - ui.checkbox(&mut self.vary_parameters.0, "Vary"); - + if let Some(slope) = &self.paramaters.slope.value { + ui.label(format!( + "Slope: {:.3} ± {:.3}", + slope, + self.paramaters.slope.uncertainity.unwrap_or(0.0) + )); + } ui.separator(); - - ui.add(egui::DragValue::new(&mut self.initial_guess.1).speed(0.1).prefix("Intercept: ")); - ui.checkbox(&mut self.vary_parameters.1, "Vary"); - }); - - ui.horizontal(|ui| { - if ui.button("Fit").clicked() { - match self.lmfit() { - Ok(_) => { - log::info!("Successfully fit data with a linear line."); - } - Err(e) => { - eprintln!("Error: {}", e); - } - } + if let Some(intercept) = &self.paramaters.intercept.value { + ui.label(format!( + "Intercept: {:.3} ± {:.3}", + intercept, + self.paramaters.intercept.uncertainity.unwrap_or(0.0) + )); } - ui.label("Best Fit:"); - ui.monospace(format!("y = {:.3}x + {:.3}", self.slope.as_ref().map_or(0.0, |v| v.value), self.intercept.as_ref().map_or(0.0, |v| v.value))); + ui.separator(); + ui.menu_button("Fit Report", |ui| { + ui.horizontal_wrapped(|ui| { + ui.label(self.fit_report.clone()); + }); + }); }); - - ui.separator(); - - ui.label("Fit Report:"); - ui.monospace(self.fit_report.clone()); } -} \ No newline at end of file +} diff --git a/src/fitter/models/mod.rs b/src/fitter/models/mod.rs index f6692ca..5da31e4 100644 --- a/src/fitter/models/mod.rs +++ b/src/fitter/models/mod.rs @@ -1,2 +1,5 @@ -pub mod gaussian; -pub mod linear; \ No newline at end of file +// pub mod gaussian; +pub mod exponential; +pub mod linear; +pub mod powerlaw; +pub mod quadratic; diff --git a/src/fitter/models/powerlaw.rs b/src/fitter/models/powerlaw.rs new file mode 100644 index 0000000..e70696e --- /dev/null +++ b/src/fitter/models/powerlaw.rs @@ -0,0 +1,217 @@ +use crate::fitter::common::{Data, Parameter}; +use pyo3::{prelude::*, types::PyModule}; + +#[derive(PartialEq, Debug, Clone, serde::Deserialize, serde::Serialize)] +pub struct PowerLawParameters { + pub amplitude: Parameter, + pub exponent: Parameter, +} + +impl Default for PowerLawParameters { + fn default() -> Self { + PowerLawParameters { + amplitude: Parameter { + name: "amplitude".to_string(), + ..Default::default() + }, + exponent: Parameter { + name: "exponent".to_string(), + ..Default::default() + }, + } + } +} + +impl PowerLawParameters { + pub fn ui(&mut self, ui: &mut egui::Ui) { + ui.horizontal(|ui| { + ui.label("Fit Parameters"); + if ui.small_button("Reset").clicked() { + *self = PowerLawParameters::default(); + } + }); + // create a grid for the param + egui::Grid::new("PowerLaw_params_grid") + .striped(true) + .num_columns(5) + .show(ui, |ui| { + ui.label("Parameter"); + ui.label("Initial Guess"); + ui.label("Min"); + ui.label("Max"); + ui.label("Vary"); + ui.end_row(); + self.amplitude.ui(ui); + ui.end_row(); + self.exponent.ui(ui); + }); + } +} + +#[derive(Debug, Clone, serde::Deserialize, serde::Serialize)] +pub struct PowerLawFitter { + pub data: Data, + pub paramaters: PowerLawParameters, + pub fit_points: Vec<[f64; 2]>, + pub fit_report: String, +} + +impl PowerLawFitter { + pub fn new(data: Data) -> Self { + PowerLawFitter { + data, + paramaters: PowerLawParameters::default(), + fit_points: Vec::new(), + fit_report: String::new(), + } + } + + pub fn lmfit(&mut self) -> PyResult<()> { + log::info!("Fitting data with a PowerLaw line using `lmfit`."); + Python::with_gil(|py| { + // let sys = py.import_bound("sys")?; + // let version: String = sys.getattr("version")?.extract()?; + // let executable: String = sys.getattr("executable")?.extract()?; + // println!("Using Python version: {}", version); + // println!("Python executable: {}", executable); + + // Check if the `uproot` module can be imported + match py.import_bound("lmfit") { + Ok(_) => { + // println!("Successfully imported `lmfit` module."); + } + Err(_) => { + eprintln!("Error: `lmfit` module could not be found. Make sure you are using the correct Python environment with `lmfit` installed."); + return Err(PyErr::new::( + "`lmfit` module not available", + )); + } + } + + match py.import_bound("numpy") { + Ok(_) => { + // println!("Successfully imported `lmfit` module."); + } + Err(_) => { + eprintln!("Error: `numpy` module could not be found. Make sure you are using the correct Python environment with `numpy` installed."); + return Err(PyErr::new::( + "`numpy` module not available", + )); + } + } + + // Define the Python code as a module + let code = r#" +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) + + 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 +"#; + + // Compile the Python code into a module + let module = PyModule::from_code_bound(py, code, "powerlaw.py", "powerlaw")?; + + let x_data = self.data.x.clone(); + let y_data = self.data.y.clone(); + let amplitude_para = ( + self.paramaters.amplitude.name.clone(), + self.paramaters.amplitude.min, + self.paramaters.amplitude.max, + self.paramaters.amplitude.initial_guess, + self.paramaters.amplitude.vary, + ); + let exponent_para = ( + self.paramaters.exponent.name.clone(), + self.paramaters.exponent.min, + self.paramaters.exponent.max, + self.paramaters.exponent.initial_guess, + self.paramaters.exponent.vary, + ); + + let result = module.getattr("PowerLawFit")?.call1(( + x_data, + y_data, + amplitude_para, + exponent_para, + ))?; + + let params = result.get_item(0)?.extract::>()?; + let x = result.get_item(1)?.extract::>()?; + let y = result.get_item(2)?.extract::>()?; + let fit_report = result.get_item(3)?.extract::()?; + + self.paramaters.amplitude.value = Some(params[0].1); + self.paramaters.amplitude.uncertainity = Some(params[0].2); + self.paramaters.exponent.value = Some(params[1].1); + self.paramaters.exponent.uncertainity = Some(params[1].2); + + self.fit_points = x.iter().zip(y.iter()).map(|(&x, &y)| [x, y]).collect(); + self.fit_report = fit_report; + + Ok(()) + }) + } + + pub fn ui(&self, ui: &mut egui::Ui) { + // add menu button for the fit report + ui.horizontal(|ui| { + if let Some(amplitude) = &self.paramaters.amplitude.value { + ui.label(format!( + "amplitude: {:.3} ± {:.3}", + amplitude, + self.paramaters.amplitude.uncertainity.unwrap_or(0.0) + )); + } + ui.separator(); + if let Some(exponent) = &self.paramaters.exponent.value { + ui.label(format!( + "exponent: {:.3} ± {:.3}", + exponent, + self.paramaters.exponent.uncertainity.unwrap_or(0.0) + )); + } + ui.separator(); + ui.menu_button("Fit Report", |ui| { + ui.horizontal_wrapped(|ui| { + ui.label(self.fit_report.clone()); + }); + }); + }); + } +} diff --git a/src/fitter/models/quadratic.rs b/src/fitter/models/quadratic.rs new file mode 100644 index 0000000..824e5cd --- /dev/null +++ b/src/fitter/models/quadratic.rs @@ -0,0 +1,248 @@ +use crate::fitter::common::{Data, Parameter}; +use pyo3::{prelude::*, types::PyModule}; + +#[derive(PartialEq, Debug, Clone, serde::Deserialize, serde::Serialize)] +pub struct QuadraticParameters { + pub a: Parameter, + pub b: Parameter, + pub c: Parameter, +} + +impl Default for QuadraticParameters { + fn default() -> Self { + QuadraticParameters { + a: Parameter { + name: "a".to_string(), + ..Default::default() + }, + b: Parameter { + name: "b".to_string(), + ..Default::default() + }, + c: Parameter { + name: "c".to_string(), + ..Default::default() + }, + } + } +} + +impl QuadraticParameters { + pub fn ui(&mut self, ui: &mut egui::Ui) { + ui.horizontal(|ui| { + ui.label("Fit Parameters"); + if ui.small_button("Reset").clicked() { + *self = QuadraticParameters::default(); + } + }); + // create a grid for the param + egui::Grid::new("linear_params_grid") + .striped(true) + .num_columns(5) + .show(ui, |ui| { + ui.label("Parameter"); + ui.label("Initial Guess"); + ui.label("Min"); + ui.label("Max"); + ui.label("Vary"); + ui.end_row(); + self.a.ui(ui); + ui.end_row(); + self.b.ui(ui); + ui.end_row(); + self.c.ui(ui); + }); + } +} + +#[derive(Debug, Clone, serde::Deserialize, serde::Serialize)] +pub struct QuadraticFitter { + pub data: Data, + pub paramaters: QuadraticParameters, + pub fit_points: Vec<[f64; 2]>, + pub fit_report: String, +} + +impl QuadraticFitter { + pub fn new(data: Data) -> Self { + QuadraticFitter { + data, + paramaters: QuadraticParameters::default(), + fit_points: Vec::new(), + fit_report: String::new(), + } + } + + pub fn lmfit(&mut self) -> PyResult<()> { + log::info!("Fitting data with a linear line using `lmfit`."); + Python::with_gil(|py| { + // let sys = py.import_bound("sys")?; + // let version: String = sys.getattr("version")?.extract()?; + // let executable: String = sys.getattr("executable")?.extract()?; + // println!("Using Python version: {}", version); + // println!("Python executable: {}", executable); + + // Check if the `uproot` module can be imported + match py.import_bound("lmfit") { + Ok(_) => { + // println!("Successfully imported `lmfit` module."); + } + Err(_) => { + eprintln!("Error: `lmfit` module could not be found. Make sure you are using the correct Python environment with `lmfit` installed."); + return Err(PyErr::new::( + "`lmfit` module not available", + )); + } + } + + match py.import_bound("numpy") { + Ok(_) => { + // println!("Successfully imported `lmfit` module."); + } + Err(_) => { + eprintln!("Error: `numpy` module could not be found. Make sure you are using the correct Python environment with `numpy` installed."); + return Err(PyErr::new::( + "`numpy` module not available", + )); + } + } + + // Define the Python code as a module + let code = r#" +import lmfit +import numpy as np + +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(result.fit_report()) + + # 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) + + 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) + + 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 params, x, y, fit_report +"#; + + // Compile the Python code into a module + let module = PyModule::from_code_bound(py, code, "quadratic.py", "quadratic")?; + + let x_data = self.data.x.clone(); + let y_data = self.data.y.clone(); + let a_para = ( + self.paramaters.a.name.clone(), + self.paramaters.a.min, + self.paramaters.a.max, + self.paramaters.a.initial_guess, + self.paramaters.a.vary, + ); + let b_para = ( + self.paramaters.b.name.clone(), + self.paramaters.b.min, + self.paramaters.b.max, + self.paramaters.b.initial_guess, + self.paramaters.b.vary, + ); + let c_para = ( + self.paramaters.c.name.clone(), + self.paramaters.c.min, + self.paramaters.c.max, + self.paramaters.c.initial_guess, + self.paramaters.c.vary, + ); + + let result = module + .getattr("QuadraticFit")? + .call1((x_data, y_data, a_para, b_para, c_para))?; + + let params = result.get_item(0)?.extract::>()?; + let x = result.get_item(1)?.extract::>()?; + let y = result.get_item(2)?.extract::>()?; + let fit_report = result.get_item(3)?.extract::()?; + + self.paramaters.a.value = Some(params[0].1); + self.paramaters.a.uncertainity = Some(params[0].2); + self.paramaters.b.value = Some(params[1].1); + self.paramaters.b.uncertainity = Some(params[1].2); + self.paramaters.c.value = Some(params[2].1); + self.paramaters.c.uncertainity = Some(params[2].2); + + self.fit_points = x.iter().zip(y.iter()).map(|(&x, &y)| [x, y]).collect(); + self.fit_report = fit_report; + + Ok(()) + }) + } + + pub fn ui(&self, ui: &mut egui::Ui) { + // add menu button for the fit report + ui.horizontal(|ui| { + if let Some(a) = &self.paramaters.a.value { + ui.label(format!( + "a: {:.3} ± {:.3}", + a, + self.paramaters.a.uncertainity.unwrap_or(0.0) + )); + } + ui.separator(); + if let Some(b) = &self.paramaters.b.value { + ui.label(format!( + "b: {:.3} ± {:.3}", + b, + self.paramaters.b.uncertainity.unwrap_or(0.0) + )); + } + ui.separator(); + if let Some(c) = &self.paramaters.c.value { + ui.label(format!( + "c: {:.3} ± {:.3}", + c, + self.paramaters.c.uncertainity.unwrap_or(0.0) + )); + } + ui.separator(); + ui.menu_button("Fit Report", |ui| { + ui.horizontal_wrapped(|ui| { + ui.label(self.fit_report.clone()); + }); + }); + }); + } +} diff --git a/src/histoer/histo1d/histogram1d.rs b/src/histoer/histo1d/histogram1d.rs index f684233..d960768 100644 --- a/src/histoer/histo1d/histogram1d.rs +++ b/src/histoer/histo1d/histogram1d.rs @@ -1,9 +1,9 @@ -use egui::Vec2b; use super::plot_settings::PlotSettings; use crate::egui_plot_stuff::egui_line::EguiLine; -use crate::fitter::fit_handler::Fits; -use crate::fitter::main_fitter::{Model, Fitter}; use crate::fitter::common::Data; +use crate::fitter::fit_handler::Fits; +use crate::fitter::main_fitter::{BackgroundModel, FitModel, Fitter}; +use egui::Vec2b; #[derive(Debug, Clone, serde::Deserialize, serde::Serialize)] pub struct Histogram { @@ -148,67 +148,66 @@ impl Histogram { .filter_map(|&pos| self.get_bin_count_and_center(pos)) .unzip(); - let mut fitter = Fitter::new( - self.fits.settings.background_model.clone(), - Data { - x: x_data, - y: y_data, - }, - ); + let mut fitter = Fitter::new(Data { + x: x_data, + y: y_data, + }); + + fitter.background_model = self.fits.settings.background_model.clone(); fitter.fit_background(); - fitter.background_line.name = format!("{} Temp Background", self.name); - - self.fits.temp_fit = Some(fitter); + fitter.name = format!("{} Temp Fit", self.name); + fitter.set_name(self.name.clone()); + self.fits.temp_fit = Some(fitter); } pub fn fit_gaussians(&mut self) { - let region_marker_positions = self.plot_settings.markers.get_region_marker_positions(); - if region_marker_positions.len() != 2 { - log::error!("Need to set two region markers to fit the histogram"); - return; - } - - // check to see if there are at least 2 background markers - if self.plot_settings.markers.get_background_marker_positions().len() >= 2 { - self.fit_background(); - } - - self.plot_settings - .markers - .remove_peak_markers_outside_region(); - let peak_positions = self.plot_settings.markers.get_peak_marker_positions(); - - let (start_x, end_x) = (region_marker_positions[0], region_marker_positions[1]); - - let data = Data { - x: self.get_bin_centers_between(start_x, end_x), - y: self.get_bin_counts_between(start_x, end_x), - }; - - let equal_stdev = self.fits.settings.equal_stddev; - let free_position = self.fits.settings.free_position; - let bin_width = self.bin_width; - - let mut fitter = Fitter::new( - Model::Gaussian(peak_positions.clone(), equal_stdev, free_position, bin_width), - data, - ); - - // Check to see if there is a temp background fit and force the background parameters if it exists - if let Some(temp_fit) = &self.fits.temp_fit { - fitter.background_line = temp_fit.background_line.clone(); - fitter.background_result = temp_fit.background_result.clone(); - fitter.background = temp_fit.model.clone(); - } - - self.fits.temp_fit = None; - - fitter.fit(); - fitter.set_name(self.name.clone()); - self.fits.temp_fit = Some(fitter); + // let region_marker_positions = self.plot_settings.markers.get_region_marker_positions(); + // if region_marker_positions.len() != 2 { + // log::error!("Need to set two region markers to fit the histogram"); + // return; + // } + + // // check to see if there are at least 2 background markers + // if self.plot_settings.markers.get_background_marker_positions().len() >= 2 { + // self.fit_background(); + // } + + // self.plot_settings + // .markers + // .remove_peak_markers_outside_region(); + // let peak_positions = self.plot_settings.markers.get_peak_marker_positions(); + + // let (start_x, end_x) = (region_marker_positions[0], region_marker_positions[1]); + + // let data = Data { + // x: self.get_bin_centers_between(start_x, end_x), + // y: self.get_bin_counts_between(start_x, end_x), + // }; + + // let equal_stdev = self.fits.settings.equal_stddev; + // let free_position = self.fits.settings.free_position; + // let bin_width = self.bin_width; + + // let mut fitter = Fitter::new( + // Model::Gaussian(peak_positions.clone(), equal_stdev, free_position, bin_width), + // data, + // ); + + // // Check to see if there is a temp background fit and force the background parameters if it exists + // if let Some(temp_fit) = &self.fits.temp_fit { + // fitter.background_line = temp_fit.background_line.clone(); + // fitter.background_result = temp_fit.background_result.clone(); + // fitter.background = temp_fit.model.clone(); + // } + + // self.fits.temp_fit = None; + + // fitter.fit(); + // fitter.set_name(self.name.clone()); + // self.fits.temp_fit = Some(fitter); } // Draw the histogram, fit lines, markers, and stats