diff --git a/lmfit_test.py b/lmfit_test.py new file mode 100644 index 0000000..fb76955 --- /dev/null +++ b/lmfit_test.py @@ -0,0 +1,125 @@ +import lmfit +import numpy as np +import polars as pl +import matplotlib.pyplot as plt + +def gaussian(x, amplitude, mean, sigma): + return amplitude * np.exp(-(x - mean)**2 / (2 * sigma**2)) + +# Multiple Gaussian fitting function +def MultipleGaussianFit(x_data: list, y_data: list, peak_markers: list, equal_sigma:bool=True, free_position:bool=True, background_type:str='linear'): + + # Initialize the model with or without a background based on the flag + if background_type == 'linear': + model = lmfit.models.LinearModel(prefix='bg_') + params = model.make_params(slope=0, intercept=0) # Initial guesses for linear background + elif background_type == 'quadratic': + model = lmfit.models.QuadraticModel(prefix='bg_') + params = model.make_params(a=0, b=1, c=0) # Initial guesses for quadratic background + elif background_type == 'exponential': + model = lmfit.models.ExponentialModel(prefix='bg_') + params = model.make_params(amplitude=1, decay=100) # Initial guesses for exponential background + elif background_type == 'powerlaw': + model = lmfit.models.PowerLawModel(prefix='bg_') + params = model.make_params(amplitude=1, exponent=-0.5) # Initial guesses for power-law background + elif background_type is None: + model = None + params = lmfit.Parameters() + else: + raise ValueError("Unsupported background model") + + first_gaussian = lmfit.Model(gaussian, prefix=f'g0_') + + if model is None: + model = first_gaussian + 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 + + # Fix the center of the first Gaussian if free_position=False + if not free_position: + params['g0_mean'].set(vary=False) + + # Add additional Gaussians + for i, peak in enumerate(peak_markers[1:], start=1): + print(f'Adding Gaussian {i} at {peak}') + g = lmfit.Model(gaussian, prefix=f'g{i}_') + model += g + params.update(g.make_params(amplitude=1, mean=peak)) + + # 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) + + # Fix the center of the Gaussian if free_position=False + if not free_position: + params[f'g{i}_mean'].set(vary=False) + + # Fit the model to the data + result = model.fit(y_data, params, x=x_data) + + print(result.fit_report()) + + # Create an array of (name, value, uncertainty) tuples + param_array = [(param, result.params[param].value, result.params[param].stderr) for param in result.params] + + + + x_result = x_data + y_result = result.best_fit + + composition_points = [x_data, result.best_fit] + + + # Return the fit result + return param_array, + +# Load the data using Polars +df = pl.read_parquet("/home/alconley/git_workspace/ICESPICE/207Bi/exp_data/207Bi_noICESPICE_f9mm_g0mm_run_13.parquet") + +df = df.with_columns([ + (pl.col("PIPS1000Energy") * 0.5395 + 2.5229).alias("PIPS1000EnergyCalibrated") +]) + + +# Create a histogram from the x_data using np.histogram +counts, bin_edges = np.histogram(df["PIPS1000EnergyCalibrated"], bins=1200, range=(0, 1200)) + +# Get the bin centers +bin_centers = 0.5 * (bin_edges[:-1] + bin_edges[1:]) + +# Filter the bin centers and counts between 1018 and 1044 +fit_mask = (bin_centers >= 520) & (bin_centers <= 600) +bin_centers_filtered = bin_centers[fit_mask] +counts_filtered = counts[fit_mask] + + +# Define the peak markers (this is where you mark the initial guesses for the peak positions) +peak_markers = [554, 567] # Replace with actual peak guesses + +# Fit the data +param = MultipleGaussianFit(bin_centers_filtered, counts_filtered, peak_markers, equal_sigma=True, free_position=False, background_type='powerlaw') + +# # Plot the original data and the fit +# plt.figure(figsize=(8, 6)) + +# # Plot the original data +# plt.step(bin_centers, counts, where="mid", label="Data") + +# # Plot the Gaussian fit only in the selected range +# plt.plot(bin_centers_filtered, result.best_fit, color="red") + +# # Add labels and legend +# plt.xlabel("PIPS1000Energy") +# plt.ylabel("Counts") +# plt.title("Gaussian Fit to Data") +# plt.legend() + +# # Show the plot +# plt.show() \ No newline at end of file diff --git a/src/fitter/background_fitter.rs b/src/fitter/background_fitter.rs deleted file mode 100644 index 13c2187..0000000 --- a/src/fitter/background_fitter.rs +++ /dev/null @@ -1,132 +0,0 @@ -use crate::egui_plot_stuff::egui_line::EguiLine; - -use super::main_fitter::{FitModel, FitResult}; -use super::models::double_exponential::DoubleExponentialFitter; -use super::models::exponential::ExponentialFitter; -use super::models::polynomial::PolynomialFitter; - -#[derive(Debug, Clone, serde::Deserialize, serde::Serialize)] -pub struct BackgroundFitter { - pub x_data: Vec, - pub y_data: Vec, - pub model: FitModel, - pub result: Option, - pub fit_line: EguiLine, -} - -impl BackgroundFitter { - pub fn new(x_data: Vec, y_data: Vec, model: FitModel) -> Self { - BackgroundFitter { - x_data, - y_data, - model, - result: None, - fit_line: EguiLine::new(egui::Color32::GREEN), - } - } - - pub fn fit(&mut self) { - match self.model { - FitModel::Gaussian(_, _, _, _) => { - log::error!("Gaussian background fitting not implemented"); - } - - FitModel::Polynomial(degree) => { - log::info!("Fitting polynomial of degree {}", degree); - let mut polynomial_fitter = PolynomialFitter::new(degree); - polynomial_fitter.x_data.clone_from(&self.x_data); - polynomial_fitter.y_data.clone_from(&self.y_data); - polynomial_fitter.fit(); - - // Update the fit line - if polynomial_fitter.coefficients.is_some() { - self.fit_line - .points - .clone_from(&polynomial_fitter.fit_line.points); - } - - self.fit_line.name = "Background".to_string(); - - self.result = Some(FitResult::Polynomial(polynomial_fitter)); - } - - FitModel::Exponential(initial_b_guess) => { - log::info!( - "Fitting exponential with initial b guess {}", - initial_b_guess - ); - let mut exponential_fitter = ExponentialFitter::new(initial_b_guess); - exponential_fitter.x_data.clone_from(&self.x_data); - exponential_fitter.y_data.clone_from(&self.y_data); - exponential_fitter.fit(); - - // Update the fit line - if exponential_fitter.coefficients.is_some() { - self.fit_line - .points - .clone_from(&exponential_fitter.fit_line.points); - } - - self.fit_line.name = "Background".to_string(); - - self.result = Some(FitResult::Exponential(exponential_fitter)); - } - - FitModel::DoubleExponential(initial_b_guess, initial_d_guess) => { - log::info!( - "Fitting double exponential with initial b guess {} and initial d guess {}", - initial_b_guess, - initial_d_guess - ); - let mut double_exponential_fitter = - DoubleExponentialFitter::new(initial_b_guess, initial_d_guess); - double_exponential_fitter.x_data.clone_from(&self.x_data); - double_exponential_fitter.y_data.clone_from(&self.y_data); - double_exponential_fitter.fit(); - - // Update the fit line - if double_exponential_fitter.coefficients.is_some() { - self.fit_line - .points - .clone_from(&double_exponential_fitter.fit_line.points); - - self.fit_line.name = "Background".to_string(); - - self.result = Some(FitResult::DoubleExponential(double_exponential_fitter)); - } - } - } - } - - pub fn fitter_stats(&self, ui: &mut egui::Ui) { - if let Some(fit) = &self.result { - match fit { - FitResult::Gaussian(fit) => fit.fit_params_ui(ui), - FitResult::Polynomial(fit) => fit.fit_params_ui(ui), - FitResult::Exponential(fit) => fit.fit_params_ui(ui), - FitResult::DoubleExponential(fit) => fit.fit_params_ui(ui), - } - } - } - - pub fn draw(&self, plot_ui: &mut egui_plot::PlotUi) { - self.fit_line.draw(plot_ui); - } - - pub fn subtract_background(&self, x_data: Vec, y_data: Vec) -> Vec { - if let Some(fit) = &self.result { - match fit { - FitResult::Polynomial(fitter) => fitter.subtract_background(x_data, y_data), - FitResult::Exponential(fitter) => fitter.subtract_background(x_data, y_data), - FitResult::DoubleExponential(fitter) => fitter.subtract_background(x_data, y_data), - _ => { - log::error!("Gaussian background fitting not implemented"); - vec![0.0; x_data.len()] - } - } - } else { - log::error!("No fit result available for background subtraction"); - vec![0.0; x_data.len()] - } - } -} diff --git a/src/fitter/fit_handler.rs b/src/fitter/fit_handler.rs index 9a99414..d2d251b 100644 --- a/src/fitter/fit_handler.rs +++ b/src/fitter/fit_handler.rs @@ -3,14 +3,14 @@ use rfd::FileDialog; use std::fs::File; use std::io::{Read, Write}; -use super::background_fitter::BackgroundFitter; +// use super::background_fitter::BackgroundFitter; use super::fit_settings::FitSettings; use super::main_fitter::Fitter; #[derive(Debug, Clone, serde::Deserialize, serde::Serialize)] pub struct Fits { pub temp_fit: Option, - pub temp_background_fit: Option, + // pub temp_background_fit: Option, pub stored_fits: Vec, pub settings: FitSettings, } @@ -25,7 +25,7 @@ impl Fits { pub fn new() -> Self { Fits { temp_fit: None, - temp_background_fit: None, + // temp_background_fit: None, stored_fits: Vec::new(), settings: FitSettings::default(), } @@ -33,7 +33,7 @@ impl Fits { pub fn store_temp_fit(&mut self) { if let Some(temp_fit) = &mut self.temp_fit.take() { - temp_fit.set_background_color(egui::Color32::DARK_GREEN); + // temp_fit.set_background_color(egui::Color32::DARK_GREEN); temp_fit.set_composition_color(egui::Color32::DARK_BLUE); temp_fit.set_decomposition_color(egui::Color32::from_rgb(150, 0, 255)); @@ -42,7 +42,7 @@ impl Fits { self.stored_fits.push(temp_fit.clone()); } - self.temp_background_fit = None; + // self.temp_background_fit = None; } pub fn set_log(&mut self, log_y: bool, log_x: bool) { @@ -50,23 +50,23 @@ impl Fits { temp_fit.set_log(log_y, log_x); } - if let Some(temp_background_fit) = &mut self.temp_background_fit { - temp_background_fit.fit_line.log_y = log_y; - temp_background_fit.fit_line.log_x = log_x; - } + // if let Some(temp_background_fit) = &mut self.temp_background_fit { + // temp_background_fit.fit_line.log_y = log_y; + // temp_background_fit.fit_line.log_x = log_x; + // } for fit in &mut self.stored_fits { fit.set_log(log_y, log_x); } } - pub fn set_stored_fits_background_color(&mut self, color: egui::Color32) { - for fit in &mut self.stored_fits { - if let Some(background) = &mut fit.background { - background.fit_line.color = color; - } - } - } + // pub fn set_stored_fits_background_color(&mut self, color: egui::Color32) { + // for fit in &mut self.stored_fits { + // if let Some(background) = &mut fit.background { + // background.fit_line.color = color; + // } + // } + // } pub fn set_stored_fits_composition_color(&mut self, color: egui::Color32) { for fit in &mut self.stored_fits { @@ -86,13 +86,13 @@ impl Fits { if let Some(temp_fit) = &mut self.temp_fit { temp_fit.show_decomposition(self.settings.show_decomposition); temp_fit.show_composition(self.settings.show_composition); - temp_fit.show_background(self.settings.show_background); + // temp_fit.show_background(self.settings.show_background); } for fit in &mut self.stored_fits { fit.show_decomposition(self.settings.show_decomposition); fit.show_composition(self.settings.show_composition); - fit.show_background(self.settings.show_background); + // fit.show_background(self.settings.show_background); } } @@ -128,7 +128,7 @@ impl Fits { serde_json::from_str(&contents).expect("Failed to deserialize fits"); self.stored_fits.extend(loaded_fits.stored_fits); // Append loaded fits to current stored fits self.temp_fit = loaded_fits.temp_fit; // override temp_fit - self.temp_background_fit = loaded_fits.temp_background_fit; // override temp_background_fit + // self.temp_background_fit = loaded_fits.temp_background_fit; // override temp_background_fit } Err(e) => { log::error!("Error opening file: {:?}", e); @@ -153,7 +153,7 @@ impl Fits { pub fn remove_temp_fits(&mut self) { self.temp_fit = None; - self.temp_background_fit = None; + // self.temp_background_fit = None; } pub fn draw(&mut self, plot_ui: &mut egui_plot::PlotUi) { @@ -163,9 +163,9 @@ impl Fits { temp_fit.draw(plot_ui); } - if let Some(temp_background_fit) = &self.temp_background_fit { - temp_background_fit.draw(plot_ui); - } + // if let Some(temp_background_fit) = &self.temp_background_fit { + // temp_background_fit.draw(plot_ui); + // } for fit in &mut self.stored_fits.iter() { fit.draw(plot_ui); @@ -193,9 +193,9 @@ impl Fits { if self.temp_fit.is_some() { ui.label("Current"); - if let Some(temp_fit) = &self.temp_fit { - temp_fit.fitter_stats(ui); - } + // if let Some(temp_fit) = &self.temp_fit { + // temp_fit.fitter_stats(ui); + // } } if !self.stored_fits.is_empty() { @@ -211,7 +211,7 @@ impl Fits { ui.separator(); }); - fit.fitter_stats(ui); + // fit.fitter_stats(ui); } } }); diff --git a/src/fitter/fit_settings.rs b/src/fitter/fit_settings.rs index 0ceab9e..ce6cdd6 100644 --- a/src/fitter/fit_settings.rs +++ b/src/fitter/fit_settings.rs @@ -9,7 +9,7 @@ pub struct FitSettings { pub fit_stats_height: f32, pub free_stddev: bool, pub free_position: bool, - pub background_model: FitModel, + // pub background_model: FitModel, pub background_poly_degree: usize, pub background_single_exp_initial_guess: f64, pub background_double_exp_initial_guess: (f64, f64), @@ -25,7 +25,7 @@ impl Default for FitSettings { fit_stats_height: 0.0, free_stddev: false, free_position: true, - background_model: FitModel::Polynomial(1), + // background_model: FitModel::Polynomial(1), background_poly_degree: 1, background_single_exp_initial_guess: 200.0, background_double_exp_initial_guess: (200.0, 800.0), @@ -75,70 +75,70 @@ impl FitSettings { ui.separator(); ui.heading("Background Fit Models"); - ui.label("Polynomial"); - ui.horizontal(|ui| { - ui.radio_value( - &mut self.background_model, - FitModel::Polynomial(1), - "Linear", - ); - ui.radio_value( - &mut self.background_model, - FitModel::Polynomial(2), - "Quadratic", - ); - ui.radio_value(&mut self.background_model, FitModel::Polynomial(3), "Cubic"); - ui.radio_value( - &mut self.background_model, - FitModel::Polynomial(self.background_poly_degree), - "n", - ); - ui.add( - egui::DragValue::new(&mut self.background_poly_degree) - .speed(1) - .prefix("Degree: ") - .range(1.0..=f32::INFINITY), - ); - }); - - ui.label("Exponential"); - ui.horizontal(|ui| { - ui.radio_value( - &mut self.background_model, - FitModel::Exponential(self.background_single_exp_initial_guess), - "Single", - ); - - ui.add( - egui::DragValue::new(&mut self.background_single_exp_initial_guess) - .speed(10) - .prefix("b: "), - ); - - ui.radio_value( - &mut self.background_model, - FitModel::DoubleExponential( - self.background_double_exp_initial_guess.0, - self.background_double_exp_initial_guess.1, - ), - "Double", - ); - - ui.add( - egui::DragValue::new(&mut self.background_double_exp_initial_guess.0) - .speed(10) - .prefix("b: ") - .range(0.0..=f64::INFINITY), - ); - - ui.add( - egui::DragValue::new(&mut self.background_double_exp_initial_guess.1) - .speed(10) - .prefix("d: ") - .range(0.0..=f64::INFINITY), - ); - }); - - ui.separator(); + // ui.label("Polynomial"); + // ui.horizontal(|ui| { + // ui.radio_value( + // &mut self.background_model, + // FitModel::Polynomial(1), + // "Linear", + // ); + // ui.radio_value( + // &mut self.background_model, + // FitModel::Polynomial(2), + // "Quadratic", + // ); + // ui.radio_value(&mut self.background_model, FitModel::Polynomial(3), "Cubic"); + // ui.radio_value( + // &mut self.background_model, + // FitModel::Polynomial(self.background_poly_degree), + // "n", + // ); + // ui.add( + // egui::DragValue::new(&mut self.background_poly_degree) + // .speed(1) + // .prefix("Degree: ") + // .range(1.0..=f32::INFINITY), + // ); + // }); + + // ui.label("Exponential"); + // ui.horizontal(|ui| { + // ui.radio_value( + // &mut self.background_model, + // FitModel::Exponential(self.background_single_exp_initial_guess), + // "Single", + // ); + + // ui.add( + // egui::DragValue::new(&mut self.background_single_exp_initial_guess) + // .speed(10) + // .prefix("b: "), + // ); + + // ui.radio_value( + // &mut self.background_model, + // FitModel::DoubleExponential( + // self.background_double_exp_initial_guess.0, + // self.background_double_exp_initial_guess.1, + // ), + // "Double", + // ); + + // ui.add( + // egui::DragValue::new(&mut self.background_double_exp_initial_guess.0) + // .speed(10) + // .prefix("b: ") + // .range(0.0..=f64::INFINITY), + // ); + + // ui.add( + // egui::DragValue::new(&mut self.background_double_exp_initial_guess.1) + // .speed(10) + // .prefix("d: ") + // .range(0.0..=f64::INFINITY), + // ); + // }); + + // ui.separator(); } } diff --git a/src/fitter/main_fitter.rs b/src/fitter/main_fitter.rs index d902b24..f8514e8 100644 --- a/src/fitter/main_fitter.rs +++ b/src/fitter/main_fitter.rs @@ -1,26 +1,24 @@ -use super::models::double_exponential::DoubleExponentialFitter; -use super::models::exponential::ExponentialFitter; -use super::models::gaussian::GaussianFitter; -use super::models::polynomial::PolynomialFitter; +// use super::models::double_exponential::DoubleExponentialFitter; +// use super::models::exponential::ExponentialFitter; +// use super::models::gaussian::GaussianFitter; +// use super::models::polynomial::PolynomialFitter; use crate::egui_plot_stuff::egui_line::EguiLine; +use super::models::gaussian::{GaussianFitter, BackgroundModels}; -use crate::fitter::background_fitter::BackgroundFitter; +// use crate::fitter::background_fitter::BackgroundFitter; #[derive(Debug, Clone, serde::Deserialize, serde::Serialize, PartialEq)] pub enum FitModel { Gaussian(Vec, bool, bool, f64), // put the initial peak locations in here, free sigma, free position - Polynomial(usize), // the degree of the polynomial: 1 for linear, 2 for quadratic, etc. - Exponential(f64), // the initial guess for the exponential decay constant - DoubleExponential(f64, f64), // the initial guess for the exponential decay constants + // Polynomial(usize), // the degree of the polynomial: 1 for linear, 2 for quadratic, etc. + // Exponential(f64), // the initial guess for the exponential decay constant + // DoubleExponential(f64, f64), // the initial guess for the exponential decay constants } #[derive(Debug, Clone, serde::Deserialize, serde::Serialize)] pub enum FitResult { Gaussian(GaussianFitter), - Polynomial(PolynomialFitter), - Exponential(ExponentialFitter), - DoubleExponential(DoubleExponentialFitter), } #[derive(Debug, Clone, serde::Deserialize, serde::Serialize)] pub struct Fitter { @@ -28,8 +26,8 @@ pub struct Fitter { pub x_data: Vec, pub y_data: Vec, pub y_err: Option>, - pub background: Option, pub model: FitModel, + pub background_model: BackgroundModels, pub result: Option, pub decomposition_lines: Vec, pub composition_line: EguiLine, @@ -37,39 +35,20 @@ pub struct Fitter { impl Fitter { // Constructor to create a new Fitter with empty data and specified model - pub fn new(model: FitModel, background: Option) -> Self { + pub fn new(model: FitModel) -> Self { Fitter { name: "Fit".to_string(), x_data: Vec::new(), y_data: Vec::new(), y_err: None, - background, model, + background_model: BackgroundModels::Linear, result: None, decomposition_lines: Vec::new(), composition_line: EguiLine::default(), } } - fn subtract_background(&self) -> Vec { - if let Some(bg_fitter) = &self.background { - match &bg_fitter.result { - Some(FitResult::Polynomial(fitter)) => { - fitter.subtract_background(self.x_data.clone(), self.y_data.clone()) - } - Some(FitResult::Exponential(fitter)) => { - fitter.subtract_background(self.x_data.clone(), self.y_data.clone()) - } - Some(FitResult::DoubleExponential(fitter)) => { - fitter.subtract_background(self.x_data.clone(), self.y_data.clone()) - } - _ => self.y_data.clone(), - } - } else { - self.y_data.clone() - } - } - pub fn get_peak_markers(&self) -> Vec { if let Some(FitResult::Gaussian(fit)) = &self.result { fit.peak_markers.clone() @@ -81,151 +60,39 @@ impl Fitter { } pub fn fit(&mut self) { - // Fit the background if it's defined and there is no background result - if let Some(bg_fitter) = &mut self.background { - if bg_fitter.result.is_none() { - bg_fitter.fit(); - } - } - - // Perform the background subtraction if necessary - let y_data_corrected = self.subtract_background(); // Perform the fit based on the model match &self.model { - FitModel::Gaussian(peak_markers, free_stddev, free_position, bin_width) => { + FitModel::Gaussian(peak_markers, equal_sigma, free_position, bin_width ) => { // Perform Gaussian fit let mut fit = GaussianFitter::new( self.x_data.clone(), - y_data_corrected, + self.y_data.clone(), peak_markers.clone(), - *free_stddev, + self.background_model.clone(), + *equal_sigma, *free_position, *bin_width, ); - fit.multi_gauss_fit(); - - // get the fit_lines and store them in the decomposition_lines - let decomposition_default_color = egui::Color32::from_rgb(255, 0, 255); - if let Some(fit_lines) = &fit.fit_lines { - for (i, line) in fit_lines.iter().enumerate() { - let mut fit_line = EguiLine::new(decomposition_default_color); - fit_line.name = format!("Peak {}", i); - - fit_line.points.clone_from(line); - fit_line.name_in_legend = false; - fit_line.width = 1.0; - self.decomposition_lines.push(fit_line); - } - } - - // calculate the composition line - if let Some(background) = &self.background { - match &background.result { - Some(FitResult::Polynomial(fitter)) => { - if let Some(coef) = &fitter.coefficients { - let composition_points = - fit.composition_fit_points_polynomial(coef.clone()); - let mut line = EguiLine::new(egui::Color32::BLUE); - line.name = "Composition".to_string(); - line.points = composition_points; - line.width = 1.0; - self.composition_line = line; - } - } - Some(FitResult::Exponential(fitter)) => { - if let Some(coef) = &fitter.coefficients { - let a = coef.a.value; - let b = coef.b.value; - let composition_points = - fit.composition_fit_points_exponential(a, b); - let mut line = EguiLine::new(egui::Color32::BLUE); - line.name = "Composition".to_string(); - line.points = composition_points; - line.width = 1.0; - self.composition_line = line; - } - } - Some(FitResult::DoubleExponential(fitter)) => { - if let Some(coef) = &fitter.coefficients { - let a = coef.a.value; - let b = coef.b.value; - let c = coef.c.value; - let d = coef.d.value; - let composition_points = - fit.composition_fit_points_double_exponential(a, b, c, d); - let mut line = EguiLine::new(egui::Color32::BLUE); - line.name = "Composition".to_string(); - line.points = composition_points; - line.width = 1.0; - self.composition_line = line; - } - } - _ => {} - } - // if let Some((slope, intercept)) = background.get_slope_intercept() { - // let composition_points = - // fit.composition_fit_points_linear_bg(slope, intercept); - - // let mut line = EguiLine::new(egui::Color32::BLUE); - // line.name = "Composition".to_string(); - // line.points = composition_points; - // self.composition_line = line; - // } - } - - self.result = Some(FitResult::Gaussian(fit)); - } - - FitModel::Polynomial(degree) => { - // Perform Polynomial fit - let mut fit = PolynomialFitter::new(*degree); - fit.x_data.clone_from(&self.x_data); - fit.y_data.clone_from(&y_data_corrected); - fit.fit(); - - self.result = Some(FitResult::Polynomial(fit)); - } - - FitModel::Exponential(initial_b_guess) => { - // Perform Exponential fit - let mut fit = ExponentialFitter::new(*initial_b_guess); - fit.x_data.clone_from(&self.x_data); - fit.y_data.clone_from(&y_data_corrected); - fit.fit(); - - self.result = Some(FitResult::Exponential(fit)); - } - - FitModel::DoubleExponential(initial_b_guess, initial_d_guess) => { - // Perform Double Exponential fit - let mut fit = DoubleExponentialFitter::new(*initial_b_guess, *initial_d_guess); - fit.x_data.clone_from(&self.x_data); - fit.y_data.clone_from(&y_data_corrected); - fit.fit(); - - self.result = Some(FitResult::DoubleExponential(fit)); + fit.fit_with_lmfit(); } } } - pub fn fitter_stats(&self, ui: &mut egui::Ui) { - if let Some(fit) = &self.result { - match fit { - FitResult::Gaussian(fit) => fit.fit_params_ui(ui), - FitResult::Polynomial(fit) => fit.fit_params_ui(ui), - FitResult::Exponential(fit) => fit.fit_params_ui(ui), - FitResult::DoubleExponential(fit) => fit.fit_params_ui(ui), - } - } - } + // pub fn fitter_stats(&self, ui: &mut egui::Ui) { + // if let Some(fit) = &self.result { + // match fit { + // FitResult::Gaussian(fit) => fit.fit_params_ui(ui), + // } + // } + // } - pub fn set_background_color(&mut self, color: egui::Color32) { - if let Some(background) = &mut self.background { - background.fit_line.color = color; - } - } + // pub fn set_background_color(&mut self, color: egui::Color32) { + // if let Some(background) = &mut self.background { + // background.fit_line.color = color; + // } + // } pub fn set_composition_color(&mut self, color: egui::Color32) { self.composition_line.color = color; @@ -247,11 +114,11 @@ impl Fitter { self.composition_line.draw = show; } - pub fn show_background(&mut self, show: bool) { - if let Some(background) = &mut self.background { - background.fit_line.draw = show; - } - } + // pub fn show_background(&mut self, show: bool) { + // if let Some(background) = &mut self.background { + // background.fit_line.draw = show; + // } + // } pub fn set_name(&mut self, name: String) { self.composition_line.name = format!("{}-Composition", name); @@ -260,15 +127,15 @@ impl Fitter { line.name = format!("{}-Peak {}", name, i); } - if let Some(background) = &mut self.background { - background.fit_line.name = format!("{}-Background", name); - } + // if let Some(background) = &mut self.background { + // background.fit_line.name = format!("{}-Background", name); + // } } pub fn lines_ui(&mut self, ui: &mut egui::Ui) { - if let Some(background) = &mut self.background { - background.fit_line.menu_button(ui); - } + // if let Some(background) = &mut self.background { + // background.fit_line.menu_button(ui); + // } self.composition_line.menu_button(ui); @@ -286,10 +153,10 @@ impl Fitter { line.draw(plot_ui); } - // Draw the background if it exists - if let Some(background) = &self.background { - background.draw(plot_ui); - } + // // Draw the background if it exists + // if let Some(background) = &self.background { + // background.draw(plot_ui); + // } // Draw the composition line self.composition_line.draw(plot_ui); @@ -302,10 +169,10 @@ impl Fitter { line.log_x = log_x; } - if let Some(background) = &mut self.background { - background.fit_line.log_y = log_y; - background.fit_line.log_x = log_x; - } + // if let Some(background) = &mut self.background { + // background.fit_line.log_y = log_y; + // background.fit_line.log_x = log_x; + // } self.composition_line.log_y = log_y; self.composition_line.log_x = log_x; diff --git a/src/fitter/mod.rs b/src/fitter/mod.rs index 4e987af..8ff4085 100644 --- a/src/fitter/mod.rs +++ b/src/fitter/mod.rs @@ -1,5 +1,4 @@ -pub mod background_fitter; pub mod fit_handler; pub mod fit_settings; pub mod main_fitter; -pub mod models; +pub mod models; \ No newline at end of file diff --git a/src/fitter/models/double_exponential.rs b/src/fitter/models/double_exponential.rs deleted file mode 100644 index 53d5f9f..0000000 --- a/src/fitter/models/double_exponential.rs +++ /dev/null @@ -1,250 +0,0 @@ -use crate::egui_plot_stuff::egui_line::EguiLine; - -use nalgebra::DVector; -use varpro::model::builder::SeparableModelBuilder; -use varpro::solvers::levmar::{LevMarProblemBuilder, LevMarSolver}; - -#[derive(Debug, Clone, serde::Deserialize, serde::Serialize)] -pub struct Coefficient { - pub value: f64, - pub uncertainty: f64, -} - -#[derive(Debug, Clone, serde::Deserialize, serde::Serialize)] -pub struct Coefficients { - pub a: Coefficient, - pub b: Coefficient, - pub c: Coefficient, - pub d: Coefficient, -} - -#[derive(Debug, Clone, serde::Deserialize, serde::Serialize)] -pub struct DoubleExponentialFitter { - pub x_data: Vec, - pub y_data: Vec, - pub weights: Vec, - pub initial_b_guess: f64, - pub initial_d_guess: f64, - pub coefficients: Option, - pub fit_line: EguiLine, -} - -impl DoubleExponentialFitter { - /// Creates a new ExponentialFitter with the given data. - pub fn new(initial_b_guess: f64, initial_d_guess: f64) -> Self { - let mut fit_line = EguiLine::new(egui::Color32::GREEN); - fit_line.name = "Double Exponential Fit".to_string(); - fit_line.width = 1.0; - - DoubleExponentialFitter { - x_data: Vec::new(), - y_data: Vec::new(), - weights: Vec::new(), - initial_b_guess, - initial_d_guess, - coefficients: None, - fit_line, - } - } - - fn exponential(x: &DVector, b: f64) -> DVector { - x.map(|x_val| (-x_val / b).exp()) - } - - fn exponential_pd_b(x: &DVector, b: f64) -> DVector { - x.map(|x_val| (x_val / b.powi(2)) * (-x_val / b).exp()) - } - - fn exponential_pd_d(x: &DVector, d: f64) -> DVector { - x.map(|x_val| (x_val / d.powi(2)) * (-x_val / d).exp()) - } - - pub fn fit(&mut self) { - let x_data = DVector::from_vec(self.x_data.clone()); - let y_data = DVector::from_vec(self.y_data.clone()); - // let weights = DVector::from_vec(self.weights.clone()); - - if x_data.len() < 4 { - log::error!("Not enough data points to fit exponential"); - return; - } - - let parameter_names: Vec = vec!["b".to_string(), "d".to_string()]; - - let intitial_parameters = vec![self.initial_b_guess, self.initial_d_guess]; - - let builder_proxy = SeparableModelBuilder::::new(parameter_names) - .initial_parameters(intitial_parameters) - .independent_variable(x_data) - .function(&["b"], Self::exponential) - .partial_deriv("b", Self::exponential_pd_b) - .function(&["d"], Self::exponential) - .partial_deriv("d", Self::exponential_pd_d); - - let model = match builder_proxy.build() { - Ok(model) => model, - Err(err) => { - log::error!("Error building model: {}", err); - return; - } - }; - - let problem = match LevMarProblemBuilder::new(model) - .observations(y_data) - // .weights(weights) - .build() - { - Ok(problem) => problem, - Err(err) => { - log::error!("Error building problem: {}", err); - return; - } - }; - - if let Ok((fit_result, fit_statistics)) = - LevMarSolver::default().fit_with_statistics(problem) - { - log::info!("fit_result: {:?}", fit_result); - log::info!("fit_statistics: {:?}", fit_statistics); - log::info!( - "Weighted residuals: {:?}", - fit_statistics.weighted_residuals() - ); - log::info!( - "Regression standard error: {:?}", - fit_statistics.regression_standard_error() - ); - log::info!( - "Covariance matrix: {:?}\n", - fit_statistics.covariance_matrix() - ); - - let nonlinear_parameters = fit_result.nonlinear_parameters(); - log::info!("nonlinear_parameters: {:?}", nonlinear_parameters); - - let nonlinear_variances = fit_statistics.nonlinear_parameters_variance(); - - let linear_coefficients = fit_result.linear_coefficients(); - - let linear_coefficients = match linear_coefficients { - Some(coefficients) => coefficients, - None => { - log::error!("No linear coefficients found"); - return; - } - }; - - log::info!("linear_coefficients: {:?}", linear_coefficients); - - let linear_variances = fit_statistics.linear_coefficients_variance(); - - let parameter_a = linear_coefficients[0]; - let parameter_a_variance = linear_variances[0]; - let parameter_a_uncertainity = parameter_a_variance.sqrt(); - - let parameter_b = nonlinear_parameters[0]; - let parameter_b_variance = nonlinear_variances[0]; - let parameter_b_uncertainity = parameter_b_variance.sqrt(); - - let parameter_c = linear_coefficients[1]; - let parameter_c_variance = linear_variances[1]; - let parameter_c_uncertainity = parameter_c_variance.sqrt(); - - let parameter_d = nonlinear_parameters[1]; - let parameter_d_variance = nonlinear_variances[1]; - let parameter_d_uncertainity = parameter_d_variance.sqrt(); - - self.coefficients = Some(Coefficients { - a: Coefficient { - value: parameter_a, - uncertainty: parameter_a_uncertainity, - }, - b: Coefficient { - value: parameter_b, - uncertainty: parameter_b_uncertainity, - }, - c: Coefficient { - value: parameter_c, - uncertainty: parameter_c_uncertainity, - }, - d: Coefficient { - value: parameter_d, - uncertainty: parameter_d_uncertainity, - }, - }); - - self.compute_fit_points(); - } - } - - fn compute_fit_points(&mut self) { - if let Some(coefficients) = &self.coefficients { - let a = coefficients.a.value; - let b = coefficients.b.value; - let c = coefficients.c.value; - let d = coefficients.d.value; - - let x_min = self.x_data.iter().cloned().fold(f64::INFINITY, f64::min); - let x_max = self - .x_data - .iter() - .cloned() - .fold(f64::NEG_INFINITY, f64::max); - - let number_points = 1000; - for i in 0..number_points { - let x = x_min + (x_max - x_min) / (number_points as f64) * (i as f64); - let y = a * (-x / b).exp() + c * (-x / d).exp(); - self.fit_line.add_point(x, y); - } - } - } - - pub fn subtract_background(&self, x_data: Vec, y_data: Vec) -> Vec { - if let Some(coefficients) = &self.coefficients { - let a = coefficients.a.value; - let b = coefficients.b.value; - let c = coefficients.c.value; - let d = coefficients.d.value; - - let mut y_data = y_data.clone(); - - for (i, x) in x_data.iter().enumerate() { - let y = a * (-x / b).exp() + c * (-x / d).exp(); - y_data[i] -= y; - } - - y_data - } else { - y_data - } - } - - fn _draw(&self, plot_ui: &mut egui_plot::PlotUi) { - self.fit_line.draw(plot_ui); - } - - pub fn fit_params_ui(&self, ui: &mut egui::Ui) { - ui.label("Coefficients:"); - if let Some(coef) = &self.coefficients { - ui.label(format!( - "a: {:.3} ± {:.3}", - coef.a.value, coef.a.uncertainty - )); - ui.label(format!( - "b: {:.3} ± {:.3}", - coef.b.value, coef.b.uncertainty - )); - ui.label(format!( - "c: {:.3} ± {:.3}", - coef.c.value, coef.c.uncertainty - )); - ui.label(format!( - "d: {:.3} ± {:.3}", - coef.d.value, coef.d.uncertainty - )); - } else { - ui.label("No coefficients found"); - } - } -} diff --git a/src/fitter/models/exponential.rs b/src/fitter/models/exponential.rs deleted file mode 100644 index beb73e1..0000000 --- a/src/fitter/models/exponential.rs +++ /dev/null @@ -1,209 +0,0 @@ -use crate::egui_plot_stuff::egui_line::EguiLine; - -use nalgebra::DVector; -use varpro::model::builder::SeparableModelBuilder; -use varpro::solvers::levmar::{LevMarProblemBuilder, LevMarSolver}; - -#[derive(Debug, Clone, serde::Deserialize, serde::Serialize)] -pub struct Coefficient { - pub value: f64, - pub uncertainty: f64, -} - -#[derive(Debug, Clone, serde::Deserialize, serde::Serialize)] -pub struct Coefficients { - pub a: Coefficient, - pub b: Coefficient, -} - -#[derive(Debug, Clone, serde::Deserialize, serde::Serialize)] -pub struct ExponentialFitter { - pub x_data: Vec, - pub y_data: Vec, - pub weights: Vec, - pub initial_b_guess: f64, - pub coefficients: Option, - pub fit_line: EguiLine, -} - -impl ExponentialFitter { - /// Creates a new ExponentialFitter with the given data. - pub fn new(initial_b_guess: f64) -> Self { - let mut fit_line = EguiLine::new(egui::Color32::GREEN); - fit_line.name = "Exponential Fit".to_string(); - fit_line.width = 1.0; - - ExponentialFitter { - x_data: Vec::new(), - y_data: Vec::new(), - weights: Vec::new(), - initial_b_guess, - coefficients: None, - fit_line, - } - } - - fn exponential(x: &DVector, b: f64) -> DVector { - x.map(|x_val| (-x_val / b).exp()) - } - - fn exponential_pd_b(x: &DVector, b: f64) -> DVector { - x.map(|x_val| (x_val / b.powi(2)) * (-x_val / b).exp()) - } - - pub fn fit(&mut self) { - let x_data = DVector::from_vec(self.x_data.clone()); - let y_data = DVector::from_vec(self.y_data.clone()); - // let weights = DVector::from_vec(self.weights.clone()); - - if x_data.len() < 3 { - log::error!("Not enough data points to fit exponential"); - return; - } - - let parameter_names: Vec = vec!["b".to_string()]; - - let intitial_parameters = vec![self.initial_b_guess]; - - let builder_proxy = SeparableModelBuilder::::new(parameter_names) - .initial_parameters(intitial_parameters) - .independent_variable(x_data) - .function(&["b"], Self::exponential) - .partial_deriv("b", Self::exponential_pd_b); - - let model = match builder_proxy.build() { - Ok(model) => model, - Err(err) => { - log::error!("Error building model: {}", err); - return; - } - }; - - let problem = match LevMarProblemBuilder::new(model) - .observations(y_data) - // .weights(weights) - .build() - { - Ok(problem) => problem, - Err(err) => { - log::error!("Error building problem: {}", err); - return; - } - }; - - if let Ok((fit_result, fit_statistics)) = - LevMarSolver::default().fit_with_statistics(problem) - { - log::info!("fit_result: {:?}", fit_result); - log::info!("fit_statistics: {:?}", fit_statistics); - log::info!( - "Weighted residuals: {:?}", - fit_statistics.weighted_residuals() - ); - log::info!( - "Regression standard error: {:?}", - fit_statistics.regression_standard_error() - ); - log::info!( - "Covariance matrix: {:?}\n", - fit_statistics.covariance_matrix() - ); - - let nonlinear_parameters = fit_result.nonlinear_parameters(); - let nonlinear_variances = fit_statistics.nonlinear_parameters_variance(); - - let linear_coefficients = fit_result.linear_coefficients(); - - let linear_coefficients = match linear_coefficients { - Some(coefficients) => coefficients, - None => { - log::error!("No linear coefficients found"); - return; - } - }; - let linear_variances = fit_statistics.linear_coefficients_variance(); - - let parameter_a = linear_coefficients[0]; - let parameter_a_variance = linear_variances[0]; - let parameter_a_uncertainity = parameter_a_variance.sqrt(); - - let parameter_b = nonlinear_parameters[0]; - let parameter_b_variance = nonlinear_variances[0]; - let parameter_b_uncertainity = parameter_b_variance.sqrt(); - - self.coefficients = Some(Coefficients { - a: Coefficient { - value: parameter_a, - uncertainty: parameter_a_uncertainity, - }, - b: Coefficient { - value: parameter_b, - uncertainty: parameter_b_uncertainity, - }, - }); - - self.compute_fit_points(); - } - } - - fn compute_fit_points(&mut self) { - if let Some(coefficients) = &self.coefficients { - let a = coefficients.a.value; - let b = coefficients.b.value; - - let x_min = self.x_data.iter().cloned().fold(f64::INFINITY, f64::min); - let x_max = self - .x_data - .iter() - .cloned() - .fold(f64::NEG_INFINITY, f64::max); - - let number_points = 1000; - for i in 0..number_points { - let x = x_min + (x_max - x_min) / (number_points as f64) * (i as f64); - let y = a * (-x / b).exp(); - self.fit_line.add_point(x, y); - } - } - } - - pub fn subtract_background(&self, x_data: Vec, y_data: Vec) -> Vec { - if let Some(coef) = &self.coefficients { - if coef.a.value == 0.0 { - log::error!("No coefficients found for exponential fit"); - return y_data; - } - - let mut y_data = y_data.clone(); - - for (i, x) in x_data.iter().enumerate() { - let y = coef.a.value * (-x / coef.b.value).exp(); - y_data[i] -= y; - } - - y_data - } else { - y_data - } - } - - fn _draw(&self, plot_ui: &mut egui_plot::PlotUi) { - self.fit_line.draw(plot_ui); - } - - pub fn fit_params_ui(&self, ui: &mut egui::Ui) { - ui.label("Coefficients:"); - if let Some(coef) = &self.coefficients { - ui.label(format!( - "a: {:.3} ± {:.3}", - coef.a.value, coef.a.uncertainty - )); - ui.label(format!( - "b: {:.3} ± {:.3}", - coef.b.value, coef.b.uncertainty - )); - } else { - ui.label("No coefficients found"); - } - } -} diff --git a/src/fitter/models/gaussian.rs b/src/fitter/models/gaussian.rs index 47c8824..9278bac 100644 --- a/src/fitter/models/gaussian.rs +++ b/src/fitter/models/gaussian.rs @@ -1,661 +1,176 @@ -use nalgebra::DVector; -use prettytable::{row, Table}; -use varpro::model::builder::SeparableModelBuilder; -use varpro::model::SeparableModel; -use varpro::solvers::levmar::{FitResult, LevMarProblemBuilder, LevMarSolver}; -use varpro::statistics::FitStatistics; // Import prettytable - -#[derive(Default, Clone, Debug, serde::Serialize, serde::Deserialize)] -pub struct Value { - pub value: f64, - pub uncertainty: f64, -} - -#[derive(Default, Clone, Debug, serde::Serialize, serde::Deserialize)] -pub struct GaussianParams { - pub amplitude: Value, - pub mean: Value, - pub sigma: Value, - pub fwhm: Value, - pub area: Value, +use pyo3::{prelude::*, types::PyModule}; + +#[derive(Debug, Clone, serde::Deserialize, serde::Serialize)] +pub enum BackgroundModels { + None, + Linear, + Quadratic, + Exponential, + Power, } -impl GaussianParams { - pub fn new(amplitude: Value, mean: Value, sigma: Value, bin_width: f64) -> Option { - if sigma.value < 0.0 { - log::error!("Sigma value is negative"); - return None; - } - - let fwhm = Self::calculate_fwhm(sigma.value); - let fwhm_uncertainty = Self::fwhm_uncertainty(sigma.uncertainty); - - let area = Self::calculate_area(amplitude.value, sigma.value, bin_width); - if area < 0.0 { - log::error!("Area is negative"); - return None; - } - let area_uncertainty = Self::area_uncertainty(amplitude.clone(), sigma.clone()); - - Some(GaussianParams { - amplitude, - mean, - sigma, - fwhm: Value { - value: fwhm, - uncertainty: fwhm_uncertainty, - }, - area: Value { - value: area, - uncertainty: area_uncertainty, - }, - }) - } - - // Method to calculate FWHM - fn calculate_fwhm(sigma: f64) -> f64 { - 2.0 * (2.0 * f64::ln(2.0)).sqrt() * sigma - } - - // Method to calculate FWHM uncertainty - fn fwhm_uncertainty(sigma_uncertainty: f64) -> f64 { - 2.0 * (2.0 * f64::ln(2.0)).sqrt() * sigma_uncertainty - } - - // Method to calculate area - fn calculate_area(amplitude: f64, sigma: f64, bin_width: f64) -> f64 { - amplitude * sigma * (2.0 * std::f64::consts::PI).sqrt() / bin_width - } - - // Method to calculate area uncertainty - fn area_uncertainty(amplitude: Value, sigma: Value) -> f64 { - let two_pi_sqrt = (2.0 * std::f64::consts::PI).sqrt(); - ((sigma.value * two_pi_sqrt * amplitude.uncertainty).powi(2) - + (amplitude.value * two_pi_sqrt * sigma.uncertainty).powi(2)) - .sqrt() - } - - pub fn params_ui(&self, ui: &mut egui::Ui) { - ui.label(format!( - "{:.2} ± {:.2}", - self.mean.value, self.mean.uncertainty - )); - ui.label(format!( - "{:.2} ± {:.2}", - self.fwhm.value, self.fwhm.uncertainty - )); - ui.label(format!( - "{:.2} ± {:.2}", - self.area.value, self.area.uncertainty - )); - } - - pub fn fit_line_points(&self) -> Vec<[f64; 2]> { - let num_points = 1000; - let start = self.mean.value - 5.0 * self.sigma.value; // Adjust start and end to be +/- 5 sigma from the mean - let end = self.mean.value + 5.0 * self.sigma.value; - let step = (end - start) / num_points as f64; - - (0..num_points) - .map(|i| { - let x = start + step * i as f64; - let y = self.amplitude.value - * (-((x - self.mean.value).powi(2)) / (2.0 * self.sigma.value.powi(2))).exp(); - [x, y] - }) - .collect() - } -} - -#[derive(Default, Clone, Debug, serde::Serialize, serde::Deserialize)] +#[derive(Debug, Clone, serde::Deserialize, serde::Serialize)] pub struct GaussianFitter { - x: Vec, - y: Vec, + pub x_data: Vec, + pub y_data: Vec, pub peak_markers: Vec, - pub fit_params: Option>, - pub fit_lines: Option>>, - pub free_stddev: bool, // false = fit all the gaussians with the same sigma - pub free_position: bool, // false = fix the position of the gaussians to the peak_markers + pub background_model: BackgroundModels, + pub equal_stdev: bool, + pub free_position: bool, pub bin_width: f64, } impl GaussianFitter { pub fn new( - x: Vec, - y: Vec, + x_data: Vec, + y_data: Vec, peak_markers: Vec, - free_stddev: bool, + background_model: BackgroundModels, + equal_stdev: bool, free_position: bool, bin_width: f64, ) -> Self { Self { - x, - y, + x_data, + y_data, peak_markers, - fit_params: None, - fit_lines: None, - free_stddev, + background_model, + equal_stdev, free_position, bin_width, } } - fn gaussian(x: &DVector, mean: f64, sigma: f64) -> DVector { - x.map(|x_val| (-((x_val - mean).powi(2)) / (2.0 * sigma.powi(2))).exp()) - } - - fn gaussian_pd_mean(x: &DVector, mean: f64, sigma: f64) -> DVector { - x.map(|x_val| { - (x_val - mean) / sigma.powi(2) - * (-((x_val - mean).powi(2)) / (2.0 * sigma.powi(2))).exp() - }) - } - - fn gaussian_pd_std_dev(x: &DVector, mean: f64, sigma: f64) -> DVector { - x.map(|x_val| { - let exponent = -((x_val - mean).powi(2)) / (2.0 * sigma.powi(2)); - (x_val - mean).powi(2) / sigma.powi(3) * exponent.exp() - }) - } - - fn average_sigma(&self) -> f64 { - let min_x = self.x.iter().cloned().fold(f64::INFINITY, f64::min); - let max_x = self.x.iter().cloned().fold(f64::NEG_INFINITY, f64::max); - let range = max_x - min_x; - - range / (5.0 * self.peak_markers.len() as f64) - } - - pub fn multi_gauss_fit(&mut self) { - self.fit_params = None; - self.fit_lines = None; - - if self.x.len() != self.y.len() { - log::error!("x_data and y_data must have the same length"); - return; - } - - if self.peak_markers.is_empty() { - // set peak marker at the x value when y is the maximum value of the data - let max_y = self.y.iter().cloned().fold(f64::NEG_INFINITY, f64::max); - let max_x = self - .x - .iter() - .cloned() - .zip(self.y.iter().cloned()) - .find(|&(_, y)| y == max_y) - .map(|(x, _)| x) - .unwrap_or_default(); - self.peak_markers.push(max_x); - } - - let x_data = DVector::from_vec(self.x.clone()); - let y_data = DVector::from_vec(self.y.clone()); - - let (parameter_names, initial_guesses) = match (self.free_stddev, self.free_position) { - (true, true) => self.generate_free_stddev_free_position(), - (false, true) => self.generate_fixed_stddev_free_position(), - (false, false) => self.generate_fixed_stddev_fixed_position(), - (true, false) => self.generate_free_stddev_fixed_position(), - }; - - let mut builder_proxy = SeparableModelBuilder::::new(parameter_names) - .initial_parameters(initial_guesses) - .independent_variable(x_data); - - for i in 0..self.peak_markers.len() { - let sigma_param = if self.free_stddev { - format!("sigma{}", i) - } else { - "sigma".to_string() - }; - - let mean_param = format!("mean{}", i); - - // If `mean` is free, include its function and derivatives - if self.free_position { - builder_proxy = builder_proxy - .function(&[mean_param.clone(), sigma_param.clone()], Self::gaussian) - .partial_deriv(mean_param, Self::gaussian_pd_mean) - .partial_deriv(sigma_param.clone(), Self::gaussian_pd_std_dev); - } else { - // If `mean` is fixed, only include the function and sigma's derivative - let fixed_mean = self.peak_markers[i]; - builder_proxy = builder_proxy - .function( - [sigma_param.clone()], - move |x: &DVector, sigma: f64| { - x.map(|x_val| { - (-((x_val - fixed_mean).powi(2)) / (2.0 * sigma.powi(2))).exp() - }) - }, - ) - .partial_deriv(sigma_param, move |x: &DVector, sigma: f64| { - x.map(|x_val| { - let exponent = -((x_val - fixed_mean).powi(2)) / (2.0 * sigma.powi(2)); - (x_val - fixed_mean).powi(2) / sigma.powi(3) * exponent.exp() - }) - }); - } - } - - let model = match builder_proxy.build() { - Ok(model) => model, - Err(e) => { - log::error!("Failed to build model: {:?}", e); - return; - } - }; - - let problem = match LevMarProblemBuilder::new(model) - .observations(y_data) - .build() - { - Ok(problem) => problem, - Err(e) => { - log::error!("Failed to build problem: {:?}", e); - return; - } - }; - - match LevMarSolver::default().fit_with_statistics(problem) { - Ok((fit_result, fit_statistics)) => { - self.process_fit_result(fit_result, fit_statistics); - } - Err(e) => { - log::error!("Failed to fit model: {:?}", e); - } - } - } - - fn generate_free_stddev_free_position(&self) -> (Vec, Vec) { - let mut parameter_names = Vec::new(); - let mut initial_guesses = Vec::new(); - - for (i, &mean) in self.peak_markers.iter().enumerate() { - parameter_names.push(format!("mean{}", i)); - parameter_names.push(format!("sigma{}", i)); - initial_guesses.push(mean); - initial_guesses.push(self.average_sigma()); - } - - (parameter_names, initial_guesses) - } - - fn generate_fixed_stddev_free_position(&self) -> (Vec, Vec) { - let mut parameter_names = Vec::new(); - let mut initial_guesses = Vec::new(); - - for (i, &mean) in self.peak_markers.iter().enumerate() { - parameter_names.push(format!("mean{}", i)); - initial_guesses.push(mean); - } - - parameter_names.push("sigma".to_string()); - initial_guesses.push(self.average_sigma()); - - (parameter_names, initial_guesses) - } - - fn generate_fixed_stddev_fixed_position(&self) -> (Vec, Vec) { - let parameter_names = vec!["sigma".to_string()]; - let initial_guesses = vec![self.average_sigma()]; - - (parameter_names, initial_guesses) - } - - fn generate_free_stddev_fixed_position(&self) -> (Vec, Vec) { - let mut parameter_names = Vec::new(); - let mut initial_guesses = Vec::new(); - - for i in 0..self.peak_markers.len() { - parameter_names.push(format!("sigma{}", i)); - initial_guesses.push(self.average_sigma()); - } - - (parameter_names, initial_guesses) - } - - fn process_fit_result( - &mut self, - fit_result: FitResult, false>, - fit_statistics: FitStatistics>, - ) { - let nonlinear_parameters = fit_result.nonlinear_parameters(); // DVector - let nonlinear_variances = fit_statistics.nonlinear_parameters_variance(); // DVector - let linear_coefficients = match fit_result.linear_coefficients() { - Some(coefficients) => coefficients, - None => { - log::error!("Failed to get linear coefficients."); - return; - } - }; - let linear_variances = fit_statistics.linear_coefficients_variance(); // DVector - - // Check if sizes match - if nonlinear_parameters.len() != nonlinear_variances.len() { - log::error!( - "Mismatch in sizes: nonlinear_parameters ({}), nonlinear_variances ({}).", - nonlinear_parameters.len(), - nonlinear_variances.len() - ); - return; - } - - if nonlinear_parameters.is_empty() || linear_coefficients.is_empty() { - log::error!("Empty fit result: no nonlinear or linear parameters."); - return; - } + pub fn fit_with_lmfit(&mut self) -> PyResult<()> { + 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); - let mut params: Vec = Vec::new(); - let sigma_value = if self.free_stddev { - None // Free sigma, so extract it for each Gaussian - } else { - match self.extract_sigma(&nonlinear_parameters, &nonlinear_variances) { - Some(sigma) => Some(sigma), - None => { - log::error!("Failed to extract sigma value."); - return; + // Check if the `uproot` module can be imported + match py.import_bound("lmfit") { + Ok(_) => { + // println!("Successfully imported `lmfit` module."); } - } - }; - - for (i, &litude) in linear_coefficients.iter().enumerate() { - let mean_value = if self.free_position { - match self.extract_mean(i, &nonlinear_parameters, &nonlinear_variances) { - Some(mean) => mean, - None => { - log::error!("Failed to extract mean value for Gaussian {}.", i); - return; - } + 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", + )); } - } else { - Value { - value: self.peak_markers.get(i).copied().unwrap_or_default(), - uncertainty: 0.0, // Fixed position, so uncertainty is 0 - } - }; - - let sigma_value = sigma_value.clone().unwrap_or_else(|| { - self.extract_sigma_for_gaussian(i, &nonlinear_parameters, &nonlinear_variances) - .unwrap_or_else(|| { - log::error!("Failed to extract sigma for Gaussian {}.", i); - Value { - value: 0.0, - uncertainty: 0.0, - } - }) - }); - - let amplitude_variance = linear_variances.get(i).copied().unwrap_or(0.0); - - if let Some(gaussian_params) = GaussianParams::new( - Value { - value: amplitude, - uncertainty: amplitude_variance.sqrt(), - }, - mean_value, - sigma_value, - self.bin_width, - ) { - params.push(gaussian_params); - } else { - log::error!( - "Invalid Gaussian parameters for Gaussian {}. Removing peak and trying again.", - i - ); - self.peak_markers.remove(i); - self.multi_gauss_fit(); // Retry the fit after removing this Gaussian - return; } - } - - self.peak_markers.clear(); - for mean in ¶ms { - self.peak_markers.push(mean.mean.value); - } - - self.fit_params = Some(params); - self.get_fit_lines(); - - self.print_fit_statistics(&fit_statistics); - self.print_peak_info(); - } - - fn extract_mean( - &self, - index: usize, - nonlinear_parameters: &DVector, - nonlinear_variances: &DVector, - ) -> Option { - if index >= nonlinear_parameters.len() || index >= nonlinear_variances.len() { - log::error!("Index out of bounds when extracting mean."); - return None; - } - Some(Value { - value: nonlinear_parameters[index], - uncertainty: nonlinear_variances[index].sqrt(), - }) - } - - fn extract_sigma( - &self, - nonlinear_parameters: &DVector, - nonlinear_variances: &DVector, - ) -> Option { - let sigma_index = nonlinear_parameters.len() - 1; // Sigma is the last parameter when fixed - if sigma_index >= nonlinear_variances.len() { - log::error!("Sigma index out of bounds when extracting sigma."); - return None; - } - Some(Value { - value: nonlinear_parameters[sigma_index], - uncertainty: nonlinear_variances[sigma_index].sqrt(), - }) - } - fn extract_sigma_for_gaussian( - &self, - index: usize, - nonlinear_parameters: &DVector, - nonlinear_variances: &DVector, - ) -> Option { - let sigma_index = index * 2 + 1; // Free sigma, after each mean - if sigma_index >= nonlinear_parameters.len() || sigma_index >= nonlinear_variances.len() { - log::error!( - "Index out of bounds when extracting sigma for Gaussian {}.", - index - ); - return None; - } - Some(Value { - value: nonlinear_parameters[sigma_index], - uncertainty: nonlinear_variances[sigma_index].sqrt(), - }) - } - - pub fn print_peak_info(&self) { - if let Some(fit_params) = &self.fit_params { - // Create a new table - let mut table = Table::new(); - - // Add the header row - table.add_row(row!["Index", "Amplitude", "Mean", "Sigma", "FWHM", "Area"]); - - // Add each peak's parameters as rows - for (i, params) in fit_params.iter().enumerate() { - table.add_row(row![ - i, - format!( - "{:.3} ± {:.3}", - params.amplitude.value, params.amplitude.uncertainty - ), - format!("{:.3} ± {:.3}", params.mean.value, params.mean.uncertainty), - format!( - "{:.3} ± {:.3}", - params.sigma.value, params.sigma.uncertainty - ), - format!("{:.3} ± {:.3}", params.fwhm.value, params.fwhm.uncertainty), - format!("{:.3} ± {:.3}", params.area.value, params.area.uncertainty), - ]); - } - - // Print the table to the terminal - table.printstd(); - } else { - println!("No fit parameters available to display."); - } - } - - pub fn print_fit_statistics(&self, fit_statistics: &FitStatistics>) { - // Print covariance matrix - // let covariance_matrix = fit_statistics.covariance_matrix(); - // println!("Covariance Matrix:"); - // Self::print_matrix(covariance_matrix); - - // // Print correlation matrix - // let correlation_matrix = fit_statistics.calculate_correlation_matrix(); - // println!("Correlation Matrix:"); - // Self::print_matrix(&correlation_matrix); - - // Print regression standard error - let regression_standard_error = fit_statistics.regression_standard_error(); - println!( - "Regression Standard Error: {:.6}", - regression_standard_error - ); - - // Print reduced chi-squared - let reduced_chi2 = fit_statistics.reduced_chi2(); - println!("Reduced Chi-Squared: {:.6}", reduced_chi2); - } - - // // Helper function to print a matrix using prettytable - // fn print_matrix(matrix: &DMatrix) { - // let mut table = Table::new(); + // Define the Python code as a module + let code = r#" +import lmfit +import numpy as np + +def gaussian(x, amplitude, mean, sigma): + return amplitude * np.exp(-(x - mean)**2 / (2 * sigma**2)) + +# Multiple Gaussian fitting function +def MultipleGaussianFit(x_data: list, y_data: list, peak_markers: list, equal_sigma:bool=True, free_position:bool=True, background_type:str='linear'): + + # Initialize the model with or without a background based on the flag + if background_type == 'linear': + model = lmfit.models.LinearModel(prefix='bg_') + params = model.make_params(slope=0, intercept=0) # Initial guesses for linear background + elif background_type == 'quadratic': + model = lmfit.models.QuadraticModel(prefix='bg_') + params = model.make_params(a=0, b=1, c=0) # Initial guesses for quadratic background + elif background_type == 'exponential': + model = lmfit.models.ExponentialModel(prefix='bg_') + params = model.make_params(amplitude=1, decay=100) # Initial guesses for exponential background + elif background_type == 'powerlaw': + model = lmfit.models.PowerLawModel(prefix='bg_') + params = model.make_params(amplitude=1, exponent=-0.5) # Initial guesses for power-law background + elif background_type is None: + model = None + params = lmfit.Parameters() + else: + raise ValueError("Unsupported background model") + + first_gaussian = lmfit.Model(gaussian, prefix=f'g0_') + + if model is None: + model = first_gaussian + 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 + + # Fix the center of the first Gaussian if free_position=False + if not free_position: + params['g0_mean'].set(vary=False) + + # 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)) + + # 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) + + # Fix the center of the Gaussian if free_position=False + if not free_position: + params[f'g{i}_mean'].set(vary=False) + + # Fit the model to the data + result = model.fit(y_data, params, x=x_data) + + print(result.fit_report()) + + # Create an array of (name, value, uncertainty) tuples + param_array = [(param, result.params[param].value, result.params[param].stderr) for param in result.params] + + # Return the fit result + return param_array +"#; + + // Compile the Python code into a module + let module = + PyModule::from_code_bound(py, code, "gaussian.py", "gaussian")?; + + let x_data = self.x_data.clone(); + let y_data = self.y_data.clone(); + let peak_markers = self.peak_markers.clone(); + let equal_sigma = self.equal_stdev; + let free_position = self.free_position; + let background_model = match self.background_model { + BackgroundModels::None => "None", + BackgroundModels::Linear => "linear", + BackgroundModels::Quadratic => "quadratic", + BackgroundModels::Exponential => "exponential", + BackgroundModels::Power => "powerlaw", + }; - // // Iterate over the rows of the matrix - // for row in matrix.row_iter() { - // let mut table_row = Vec::new(); - // for val in row.iter() { - // table_row.push(cell!(format!("{:.6}", val))); - // } - // // Add the row using add_row, but without row! macro - // table.add_row(prettytable::Row::new(table_row)); - // } + let result = module.getattr("MultipleGaussianFit")?.call1((x_data, y_data, peak_markers, equal_sigma, free_position, background_model))?; + println!("result: {:?}", result); - // // Print the table - // table.printstd(); - // } + // [(String, f64, f64), ...] - pub fn get_fit_lines(&mut self) { - if let Some(fit_params) = &self.fit_params { - let mut fit_lines = Vec::new(); + let length: usize = result.len()?; + for i in 0..length { + let item = result.get_item(i)?; + let name: String = item.get_item(0)?.extract()?; + let value: f64 = item.get_item(1)?.extract()?; + let uncertainty: f64 = item.get_item(2)?.extract()?; - for params in fit_params.iter() { - let line = params.fit_line_points(); - fit_lines.push(line); + println!("name: {}, value: {}, uncertainty: {}", name, value, uncertainty); } - self.fit_lines = Some(fit_lines); - } else { - self.fit_lines = None; - } - } - - pub fn composition_fit_points_polynomial(&self, coef: Vec) -> Vec<[f64; 2]> { - // coef = [c0, c1, c2, ...] c0 + c1*x + c2*x^2 + ... - let num_points = 3000; - let min_x = self.x.iter().cloned().fold(f64::INFINITY, f64::min); - let max_x = self.x.iter().cloned().fold(f64::NEG_INFINITY, f64::max); - let step = (max_x - min_x) / num_points as f64; - - (0..=num_points) - .map(|i| { - let x = min_x + step * i as f64; - let y_gauss = self.fit_params.as_ref().map_or(0.0, |params| { - params.iter().fold(0.0, |sum, param| { - sum + param.amplitude.value - * (-((x - param.mean.value).powi(2)) - / (2.0 * param.sigma.value.powi(2))) - .exp() - }) - }); - let y_background = coef - .iter() - .enumerate() - .fold(0.0, |sum, (j, c)| sum + c * x.powi(j as i32)); - let y_total = y_gauss + y_background; - [x, y_total] - }) - .collect() - } - - pub fn composition_fit_points_exponential(&self, a: f64, b: f64) -> Vec<[f64; 2]> { - let num_points = 3000; - let min_x = self.x.iter().cloned().fold(f64::INFINITY, f64::min); - let max_x = self.x.iter().cloned().fold(f64::NEG_INFINITY, f64::max); - let step = (max_x - min_x) / num_points as f64; - - (0..=num_points) - .map(|i| { - let x = min_x + step * i as f64; - let y_gauss = self.fit_params.as_ref().map_or(0.0, |params| { - params.iter().fold(0.0, |sum, param| { - sum + param.amplitude.value - * (-((x - param.mean.value).powi(2)) - / (2.0 * param.sigma.value.powi(2))) - .exp() - }) - }); - let y_background = a * (-x / b).exp(); - let y_total = y_gauss + y_background; - [x, y_total] - }) - .collect() - } - - pub fn composition_fit_points_double_exponential( - &self, - a: f64, - b: f64, - c: f64, - d: f64, - ) -> Vec<[f64; 2]> { - let num_points = 3000; - let min_x = self.x.iter().cloned().fold(f64::INFINITY, f64::min); - let max_x = self.x.iter().cloned().fold(f64::NEG_INFINITY, f64::max); - let step = (max_x - min_x) / num_points as f64; - (0..=num_points) - .map(|i| { - let x = min_x + step * i as f64; - let y_gauss = self.fit_params.as_ref().map_or(0.0, |params| { - params.iter().fold(0.0, |sum, param| { - sum + param.amplitude.value - * (-((x - param.mean.value).powi(2)) - / (2.0 * param.sigma.value.powi(2))) - .exp() - }) - }); - let y_background = a * (-x / b).exp() + c * (-x / d).exp(); - let y_total = y_gauss + y_background; - [x, y_total] - }) - .collect() - } - - pub fn fit_params_ui(&self, ui: &mut egui::Ui) { - if let Some(fit_params) = &self.fit_params { - for (i, params) in fit_params.iter().enumerate() { - if i != 0 { - ui.label(""); - } - - ui.label(format!("{}", i)); - params.params_ui(ui); - ui.end_row(); - } - } + Ok(()) + }) } -} + +} \ No newline at end of file diff --git a/src/fitter/models/mod.rs b/src/fitter/models/mod.rs index 818d0a0..1e195d9 100644 --- a/src/fitter/models/mod.rs +++ b/src/fitter/models/mod.rs @@ -1,4 +1 @@ -pub mod double_exponential; -pub mod exponential; pub mod gaussian; -pub mod polynomial; diff --git a/src/fitter/models/polynomial.rs b/src/fitter/models/polynomial.rs deleted file mode 100644 index b3d5d21..0000000 --- a/src/fitter/models/polynomial.rs +++ /dev/null @@ -1,125 +0,0 @@ -use crate::egui_plot_stuff::egui_line::EguiLine; -use compute::predict::PolynomialRegressor; - -#[derive(Debug, Clone, serde::Deserialize, serde::Serialize)] -pub struct PolynomialFitter { - pub x_data: Vec, - pub y_data: Vec, - pub degree: usize, - pub coefficients: Option>, - pub fit_line: EguiLine, -} - -impl PolynomialFitter { - /// Creates a new PolynomialFitter with the given data. - pub fn new(degree: usize) -> Self { - let mut fit_line = EguiLine::new(egui::Color32::GREEN); - fit_line.name = "Polynomial Fit".to_string(); - - PolynomialFitter { - x_data: Vec::new(), - y_data: Vec::new(), - degree, - coefficients: None, - fit_line, - } - } - - pub fn fit(&mut self) { - let mut regressor = PolynomialRegressor::new(self.degree); - - if self.x_data.len() < self.degree + 1 { - log::error!( - "Not enough data points to fit polynomial of degree {}", - self.degree - ); - return; - } - regressor.fit(&self.x_data, &self.y_data); - - self.coefficients = Some(regressor.coef.clone()); - self.compute_fit_points(); - - log::info!("Polynomial fit coefficients: {:?}", regressor.coef); - } - - pub fn subtract_background(&self, x_data: Vec, y_data: Vec) -> Vec { - if let Some(coef) = &self.coefficients { - if coef.is_empty() { - log::error!("No coefficients found for polynomial fit"); - return y_data; - } - - let mut y_data = y_data.clone(); - - for (i, x) in x_data.iter().enumerate() { - let y = coef - .iter() - .enumerate() - .fold(0.0, |acc, (j, c)| acc + c * x.powi(j as i32)); - y_data[i] -= y; - } - - y_data - } else { - y_data - } - } - - fn compute_fit_points(&mut self) { - if let Some(coef) = &self.coefficients { - if coef.is_empty() { - log::error!("No coefficients found for polynomial fit"); - return; - } - - self.fit_line.clear_points(); - - // get the min and max x values - let (x_min, x_max) = self - .x_data - .iter() - .fold((f64::INFINITY, f64::NEG_INFINITY), |(min, max), &x| { - (min.min(x), max.max(x)) - }); - - let number_points = 1000; - for i in 0..number_points { - let x = x_min + (x_max - x_min) / (number_points as f64) * (i as f64); - let y = coef - .iter() - .enumerate() - .fold(0.0, |acc, (j, c)| acc + c * x.powi(j as i32)); - self.fit_line.add_point(x, y); - } - } - } - - pub fn _draw(&self, plot_ui: &mut egui_plot::PlotUi) { - self.fit_line.draw(plot_ui); - } - - pub fn fit_params_ui(&self, ui: &mut egui::Ui) { - // ui.horizontal(|ui| { - // ui.label("Polynomial degree:"); - // ui.add(egui::DragValue::new(&mut self.degree).speed(1.0)); - // }); - - // if ui.button("Fit").clicked() { - // self.fit(); - // } - - ui.label("Coefficients:"); - if let Some(coef) = &self.coefficients { - if coef.is_empty() { - ui.label("No coefficients found"); - } else { - for (i, coef) in coef.iter().enumerate() { - ui.label(format!("c{}: {}", i, coef)); - } - } - } else { - ui.label("No coefficients found"); - } - } -} diff --git a/src/histoer/histo1d/context_menu.rs b/src/histoer/histo1d/context_menu.rs index 784a270..a60891e 100644 --- a/src/histoer/histo1d/context_menu.rs +++ b/src/histoer/histo1d/context_menu.rs @@ -12,11 +12,11 @@ impl Histogram { // Add find peaks button ui.separator(); ui.heading("Peak Finder"); - if ui.button("Detect Peaks") - .on_hover_text("Takes the settings (adjust below) and finds peaks in the spectrum\nIf there are background markers, it will fit a background before it finds the peaks in between the min and max values. Likewise for region markers.\nKeybind: o").clicked() { - self.find_peaks(); - } - self.plot_settings.find_peaks_settings.menu_button(ui); + // if ui.button("Detect Peaks") + // .on_hover_text("Takes the settings (adjust below) and finds peaks in the spectrum\nIf there are background markers, it will fit a background before it finds the peaks in between the min and max values. Likewise for region markers.\nKeybind: o").clicked() { + // self.find_peaks(); + // } + // self.plot_settings.find_peaks_settings.menu_button(ui); ui.separator(); ui.heading("Rebin"); diff --git a/src/histoer/histo1d/histogram1d.rs b/src/histoer/histo1d/histogram1d.rs index a4442f3..d05210e 100644 --- a/src/histoer/histo1d/histogram1d.rs +++ b/src/histoer/histo1d/histogram1d.rs @@ -2,10 +2,12 @@ use egui::Vec2b; use super::plot_settings::PlotSettings; use crate::egui_plot_stuff::egui_line::EguiLine; -use crate::fitter::background_fitter::BackgroundFitter; +// use crate::fitter::background_fitter::BackgroundFitter; use crate::fitter::fit_handler::Fits; use crate::fitter::main_fitter::{FitModel, Fitter}; +use crate::fitter::models::gaussian::{GaussianFitter, BackgroundModels}; + #[derive(Debug, Clone, serde::Deserialize, serde::Serialize)] pub struct Histogram { pub name: String, @@ -134,28 +136,28 @@ impl Histogram { }) } - pub fn fit_background(&mut self) { - self.fits.remove_temp_fits(); + // pub fn fit_background(&mut self) { + // self.fits.remove_temp_fits(); - let marker_positions = self.plot_settings.markers.get_background_marker_positions(); - if marker_positions.len() < 2 { - log::error!("Need to set at least two background markers to fit the histogram"); - return; - } + // let marker_positions = self.plot_settings.markers.get_background_marker_positions(); + // if marker_positions.len() < 2 { + // log::error!("Need to set at least two background markers to fit the histogram"); + // return; + // } - let (x_data, y_data): (Vec, Vec) = marker_positions - .iter() - .filter_map(|&pos| self.get_bin_count_and_center(pos)) - .unzip(); + // let (x_data, y_data): (Vec, Vec) = marker_positions + // .iter() + // .filter_map(|&pos| self.get_bin_count_and_center(pos)) + // .unzip(); - // let mut background_fitter = BackgroundFitter::new(x_data, y_data, FitModel::Linear); - let mut background_fitter = - BackgroundFitter::new(x_data, y_data, self.fits.settings.background_model.clone()); - background_fitter.fit(); + // // let mut background_fitter = BackgroundFitter::new(x_data, y_data, FitModel::Linear); + // let mut background_fitter = + // BackgroundFitter::new(x_data, y_data, self.fits.settings.background_model.clone()); + // background_fitter.fit(); - background_fitter.fit_line.name = format!("{} Temp Background", self.name); - self.fits.temp_background_fit = Some(background_fitter); - } + // background_fitter.fit_line.name = format!("{} Temp Background", self.name); + // self.fits.temp_background_fit = Some(background_fitter); + // } pub fn fit_gaussians(&mut self) { let region_marker_positions = self.plot_settings.markers.get_region_marker_positions(); @@ -169,43 +171,56 @@ impl Histogram { .remove_peak_markers_outside_region(); let peak_positions = self.plot_settings.markers.get_peak_marker_positions(); - if self.fits.temp_background_fit.is_none() { - if self.plot_settings.markers.background_markers.len() <= 1 { - for position in region_marker_positions.iter() { - self.plot_settings.markers.add_background_marker(*position); - } - } - self.fit_background(); - } - - let mut fitter = Fitter::new( - FitModel::Gaussian( - peak_positions, - self.fits.settings.free_stddev, - self.fits.settings.free_position, - self.bin_width, - ), - self.fits.temp_background_fit.clone(), - ); + // if self.fits.temp_background_fit.is_none() { + // if self.plot_settings.markers.background_markers.len() <= 1 { + // for position in region_marker_positions.iter() { + // self.plot_settings.markers.add_background_marker(*position); + // } + // } + // self.fit_background(); + // } + + // let mut fitter = Fitter::new( + // FitModel::Gaussian( + // peak_positions, + // self.fits.settings.free_stddev, + // self.fits.settings.free_position, + // self.bin_width, + // ), + // // self.fits.temp_background_fit.clone(), + // ); let (start_x, end_x) = (region_marker_positions[0], region_marker_positions[1]); - fitter.x_data = self.get_bin_centers_between(start_x, end_x); - fitter.y_data = self.get_bin_counts_between(start_x, end_x); + // fitter.x_data = self.get_bin_centers_between(start_x, end_x); + // fitter.y_data = self.get_bin_counts_between(start_x, end_x); - fitter.fit(); - fitter.set_name(self.name.clone()); + let x_data = self.get_bin_centers_between(start_x, end_x); + let y_data = self.get_bin_counts_between(start_x, end_x); + let background_model = BackgroundModels::Linear; + let equal_stdev = self.fits.settings.free_stddev; + let free_position = self.fits.settings.free_position; + let bin_width = self.bin_width; - // clear peak markers and add the new peak markers - self.plot_settings.markers.clear_peak_markers(); - let peak_values = fitter.get_peak_markers(); - for peak in peak_values { - self.plot_settings.markers.add_peak_marker(peak); - } + let mut fitter = GaussianFitter::new(x_data, y_data, peak_positions.clone(), background_model, equal_stdev, free_position, bin_width); + + fitter.fit_with_lmfit(); + + // fitter.fit(); + + // fitter.set_name(self.name.clone()); + + // // clear peak markers and add the new peak markers + // self.plot_settings.markers.clear_peak_markers(); + + // let peak_values = fitter.get_peak_markers(); + // for peak in peak_values { + // self.plot_settings.markers.add_peak_marker(peak); + // } - self.fits.temp_fit = Some(fitter); + // self.fits.temp_fit = Some(fitter); } // Draw the histogram, fit lines, markers, and stats diff --git a/src/histoer/histo1d/keybinds.rs b/src/histoer/histo1d/keybinds.rs index 35832be..d68e560 100644 --- a/src/histoer/histo1d/keybinds.rs +++ b/src/histoer/histo1d/keybinds.rs @@ -41,9 +41,9 @@ impl Histogram { self.fits.remove_temp_fits(); } - if ui.input(|i| i.key_pressed(egui::Key::G)) { - self.fit_background(); - } + // if ui.input(|i| i.key_pressed(egui::Key::G)) { + // self.fit_background(); + // } if ui.input(|i| i.key_pressed(egui::Key::F)) { self.fit_gaussians(); @@ -61,9 +61,9 @@ impl Histogram { self.plot_settings.egui_settings.log_y = !self.plot_settings.egui_settings.log_y; } - if ui.input(|i| i.key_pressed(egui::Key::O)) { - self.find_peaks(); - } + // if ui.input(|i| i.key_pressed(egui::Key::O)) { + // self.find_peaks(); + // } } } diff --git a/src/histoer/histo1d/peak_finder.rs b/src/histoer/histo1d/peak_finder.rs index 2eb6a54..4a8f59e 100644 --- a/src/histoer/histo1d/peak_finder.rs +++ b/src/histoer/histo1d/peak_finder.rs @@ -5,81 +5,81 @@ use super::histogram1d::Histogram; impl Histogram { // Add a function to find peaks - pub fn find_peaks(&mut self) { - // Clear the peak markers - self.plot_settings.markers.clear_peak_markers(); - - let region_marker_positions = self.plot_settings.markers.get_region_marker_positions(); - let mut background_marker_positions = - self.plot_settings.markers.get_background_marker_positions(); - - // Sort background markers - background_marker_positions.sort_by(|a, b| a.partial_cmp(b).unwrap()); - - let mut peaks_found_with_background = false; - let mut peaks_found_with_region = false; - - let y_data: Vec = if background_marker_positions.len() >= 2 { - self.fit_background(); - - // Extract the y data from the temp background fit - if let Some(temp_background) = &self.fits.temp_background_fit { - // if there are region markers, use the data between them - let (start_x, end_x) = if region_marker_positions.len() == 2 { - peaks_found_with_region = true; - (region_marker_positions[0], region_marker_positions[1]) - } else { - peaks_found_with_background = true; - ( - background_marker_positions[0], - background_marker_positions[background_marker_positions.len() - 1], - ) - }; - - let x_data = self.get_bin_centers_between(start_x, end_x); - let y_data = self.get_bin_counts_between(start_x, end_x); - - // Put the data in the background fitter to subtract the background - temp_background.subtract_background(x_data, y_data) - } else { - log::error!("Failed to fit background"); - return; - } - } else if region_marker_positions.len() == 2 { - let (start_x, end_x) = (region_marker_positions[0], region_marker_positions[1]); - peaks_found_with_region = true; - self.get_bin_counts_between(start_x, end_x) - } else { - self.bins.iter().map(|&count| count as f64).collect() - }; - - let peaks = self.plot_settings.find_peaks_settings.find_peaks(y_data); - - // Add peak markers at detected peaks - for peak in &peaks { - let peak_position = peak.middle_position(); - log::info!("Peak at position: {}", peak_position); - // Adjust peak position relative to the first background marker - if peaks_found_with_background { - let adjusted_peak_position = - self.bin_width * peak_position as f64 + background_marker_positions[0]; - self.plot_settings - .markers - .add_peak_marker(adjusted_peak_position); - } else if peaks_found_with_region { - let adjusted_peak_position = - self.bin_width * peak_position as f64 + region_marker_positions[0]; - self.plot_settings - .markers - .add_peak_marker(adjusted_peak_position); - } else { - let adjusted_peak_position = self.bin_width * peak_position as f64 + self.range.0; - self.plot_settings - .markers - .add_peak_marker(adjusted_peak_position); - } - } - } + // pub fn find_peaks(&mut self) { + // // Clear the peak markers + // self.plot_settings.markers.clear_peak_markers(); + + // let region_marker_positions = self.plot_settings.markers.get_region_marker_positions(); + // let mut background_marker_positions = + // self.plot_settings.markers.get_background_marker_positions(); + + // // Sort background markers + // background_marker_positions.sort_by(|a, b| a.partial_cmp(b).unwrap()); + + // let mut peaks_found_with_background = false; + // let mut peaks_found_with_region = false; + + // let y_data: Vec = if background_marker_positions.len() >= 2 { + // // self.fit_background(); + + // // Extract the y data from the temp background fit + // if let Some(temp_background) = &self.fits.temp_background_fit { + // // if there are region markers, use the data between them + // let (start_x, end_x) = if region_marker_positions.len() == 2 { + // peaks_found_with_region = true; + // (region_marker_positions[0], region_marker_positions[1]) + // } else { + // peaks_found_with_background = true; + // ( + // background_marker_positions[0], + // background_marker_positions[background_marker_positions.len() - 1], + // ) + // }; + + // let x_data = self.get_bin_centers_between(start_x, end_x); + // let y_data = self.get_bin_counts_between(start_x, end_x); + + // // Put the data in the background fitter to subtract the background + // temp_background.subtract_background(x_data, y_data) + // } else { + // log::error!("Failed to fit background"); + // return; + // } + // } else if region_marker_positions.len() == 2 { + // let (start_x, end_x) = (region_marker_positions[0], region_marker_positions[1]); + // peaks_found_with_region = true; + // self.get_bin_counts_between(start_x, end_x) + // } else { + // self.bins.iter().map(|&count| count as f64).collect() + // }; + + // let peaks = self.plot_settings.find_peaks_settings.find_peaks(y_data); + + // // Add peak markers at detected peaks + // for peak in &peaks { + // let peak_position = peak.middle_position(); + // log::info!("Peak at position: {}", peak_position); + // // Adjust peak position relative to the first background marker + // if peaks_found_with_background { + // let adjusted_peak_position = + // self.bin_width * peak_position as f64 + background_marker_positions[0]; + // self.plot_settings + // .markers + // .add_peak_marker(adjusted_peak_position); + // } else if peaks_found_with_region { + // let adjusted_peak_position = + // self.bin_width * peak_position as f64 + region_marker_positions[0]; + // self.plot_settings + // .markers + // .add_peak_marker(adjusted_peak_position); + // } else { + // let adjusted_peak_position = self.bin_width * peak_position as f64 + self.range.0; + // self.plot_settings + // .markers + // .add_peak_marker(adjusted_peak_position); + // } + // } + // } } #[derive(Debug, Clone, serde::Serialize, serde::Deserialize)] diff --git a/src/histogram_scripter/manual_histogram_script.rs b/src/histogram_scripter/manual_histogram_script.rs index e1297ef..fdc88d1 100644 --- a/src/histogram_scripter/manual_histogram_script.rs +++ b/src/histogram_scripter/manual_histogram_script.rs @@ -7,14 +7,14 @@ use std::f64::consts::PI; #[allow(clippy::all)] pub fn manual_add_histograms(h: &mut Histogrammer, lf: LazyFrame) { - // let lf = lf.with_column( - // (lit(-1.77049e-06)*col("PIPS1000Energy")*col("PIPS1000Energy") + lit(0.544755003513083)*col("PIPS1000Energy") + lit(-1.36822594543883)).alias("PIPS1000EnergyCalibrated") - // ); + let lf = lf.with_column( + (lit(-1.77049e-06)*col("PIPS1000Energy")*col("PIPS1000Energy") + lit(0.544755003513083)*col("PIPS1000Energy") + lit(-1.36822594543883)).alias("PIPS1000EnergyCalibrated") + ); - // h.add_fill_hist1d("PIPS1000Energy", &lf, "PIPS1000Energy", 16384, (0.0, 16384.0), None); - // h.add_fill_hist1d("PIPS1000EnergyCalibrated", &lf, "PIPS1000EnergyCalibrated", 1200, (0.0, 1200.0), None); + h.add_fill_hist1d("PIPS1000Energy", &lf, "PIPS1000Energy", 16384, (0.0, 16384.0), None); + h.add_fill_hist1d("PIPS1000EnergyCalibrated", &lf, "PIPS1000EnergyCalibrated", 1200, (0.0, 1200.0), None); - // /* + /* // Declare all the lazyframes that will be used let lf = lf.with_columns(vec![ @@ -311,5 +311,5 @@ pub fn manual_add_histograms(h: &mut Histogrammer, lf: LazyFrame) { */ - // */ + */ }