From c950303f79f9a76ff8fa524dd70d42762277cf05 Mon Sep 17 00:00:00 2001 From: eloiseyang <98988862+eloiseyang@users.noreply.github.com> Date: Fri, 2 Aug 2024 11:55:04 -0700 Subject: [PATCH] Fitting Base Class and Gaussian Model from PR 143 (#172) * Added MethodBase as an abstract class for all fit methods Added GaussianModel class using MethodBase Added ProjectionFit class Added tests for all the above including datasets Separated from 8c85a56 on phys-cgarnier/lcls-tools/dev * method_base.py flake8 * changed code base to be more pydantic this cleans up some things involving params, might have to remove priors still but its working * updated fit code to incorporate pydantic structures for method parameters * made changes as requested in the comment (minus removing _forward) * flake8 * committing GaussianFIt with property that returns beamsize * checked that GaussianFit works (roughly) and set processor and fit to have a default field, also set method in projection_fit.py to have a default file. This will make using the tool much easier. All that has to be down is call gaussian fit and give it an image to use * committing linted gaussian_fit.py * removed empty test for projection * more flake8 --------- Co-authored-by: Chris Garnier --- .../common/data_analysis/fit/gaussian_fit.py | 126 +++++++++++++++ .../common/data_analysis/fit/method_base.py | 148 ++++++++++++++++++ .../common/data_analysis/fit/methods.py | 103 ++++++++++++ .../common/data_analysis/fit/projection.py | 108 +++++++++++++ tests/datasets/fit/make_test_guassian.py | 15 ++ tests/datasets/fit/test_gaussian.npy | Bin 0 -> 8128 bytes .../common/data_analysis/fit/test_methods.py | 97 ++++++++++++ .../common/data_analysis/test_fit_gauss.py | 71 --------- 8 files changed, 597 insertions(+), 71 deletions(-) create mode 100644 lcls_tools/common/data_analysis/fit/gaussian_fit.py create mode 100644 lcls_tools/common/data_analysis/fit/method_base.py create mode 100644 lcls_tools/common/data_analysis/fit/methods.py create mode 100644 lcls_tools/common/data_analysis/fit/projection.py create mode 100644 tests/datasets/fit/make_test_guassian.py create mode 100644 tests/datasets/fit/test_gaussian.npy create mode 100644 tests/unit_tests/lcls_tools/common/data_analysis/fit/test_methods.py delete mode 100644 tests/unit_tests/lcls_tools/common/data_analysis/test_fit_gauss.py diff --git a/lcls_tools/common/data_analysis/fit/gaussian_fit.py b/lcls_tools/common/data_analysis/fit/gaussian_fit.py new file mode 100644 index 00000000..624780a1 --- /dev/null +++ b/lcls_tools/common/data_analysis/fit/gaussian_fit.py @@ -0,0 +1,126 @@ +import numpy as np +from pydantic import BaseModel, ConfigDict, PositiveFloat +from typing import Optional +from lcls_tools.common.data_analysis.fit.projection import ProjectionFit +from lcls_tools.common.image.processing import ImageProcessor + + +class GaussianFit(BaseModel): + # should rename to BeamsizeEvaluator or something ? + """ + Gaussian Fitting class that takes in a fitting tool (with method), + ImageProcessor, and Image and returns the beamsize and position. + """ + + model_config = ConfigDict(arbitrary_types_allowed=True) + processor: ImageProcessor = ImageProcessor() + fit: ProjectionFit = ProjectionFit() + min_log_intensity: float = 4.0 + bounding_box_half_width: PositiveFloat = 3.0 + apply_bounding_box_constraint: bool = True + image: Optional[np.ndarray] = None + + @property + def beamsize(self): + beamspot_chars = self.gaussian_fit(self.image) + beamsize = self.calculate_beamsize(beamspot_chars) + return beamsize + + def gaussian_fit(self, image): + self.processor.auto_process(image) + x_proj, y_proj = self.get_projections(image) + x_parameters = self.fit.fit_projection(x_proj) + y_parameters = self.fit.fit_projection(y_proj) + + return { + "centroid": np.array((x_parameters['mean'], + y_parameters['mean'])), + "rms_sizes": np.array((x_parameters['sigma'], + y_parameters['sigma'])), + "total_intensity": image.sum(), + } + + def get_projections(self, image): + x_projection = np.array(np.sum(image, axis=0)) + y_projection = np.array(np.sum(image, axis=1)) + return x_projection, y_projection + + def calculate_beamsize(self, beamspot_chars: dict): + ''' + Conditional beamsize calculation: + if condition1 (total_intensity is to small): then return NULL result. + else: do beamsize calculation and check that the beamspot is + within the ROI + ''' + if (np.log10(beamspot_chars['total_intensity']) + < self.min_log_intensity): + # TODO: character count is really not liking this one line + print(("log10 image intensity" + + f"{np.log10(beamspot_chars['total_intensity'])}" + + "below threshold")) + result = { + "Cx": np.NaN, + "Cy": np.NaN, + "Sx": np.NaN, + "Sy": np.NaN, + "bb_penalty": np.NaN, + "total_intensity": beamspot_chars['total_intensity'], + } + return result + else: + centroid = beamspot_chars['centroid'] + sizes = beamspot_chars['rms_sizes'] + + if np.all(~np.isnan(np.stack((centroid, sizes)))): + # get beam region bounding box + n_stds = self.bounding_box_half_width + + bounding_box_corner_pts = np.array( + ( + centroid - n_stds * sizes, + centroid + n_stds * sizes, + centroid - n_stds * sizes * np.array((-1, 1)), + centroid + n_stds * sizes * np.array((-1, 1)), + ) + ) + + if self.processor.roi is not None: + roi_radius = self.processor.roi.radius + else: + roi_radius = np.min(np.array(self.image.shape) / 1.5) + # TODO: maybe need to change this ^^^ + + # This is a check whether or not the beamspot is within + # the bounding box. + temp = bounding_box_corner_pts - np.array((roi_radius, + roi_radius)) + distances = np.linalg.norm(temp, axis=1) + bounding_box_penalty = np.max(distances) - roi_radius + + result = { + "Cx": centroid[0], + "Cy": centroid[1], + "Sx": sizes[0], + "Sy": sizes[1], + "bb_penalty": bounding_box_penalty, + "total_intensity": beamspot_chars["total_intensity"] + } + + # set results to none if the beam extends beyond the roi + # and the bounding box constraint is active + if (bounding_box_penalty > 0 + and self.apply_bounding_box_constraint): + for name in ["Cx", "Cy", "Sx", "Sy"]: + result[name] = np.NaN + + else: + result = { + "Cx": np.NaN, + "Cy": np.NaN, + "Sx": np.NaN, + "Sy": np.NaN, + "bb_penalty": np.NaN, + "total_intensity": beamspot_chars["total_intensity"] + } + + return result diff --git a/lcls_tools/common/data_analysis/fit/method_base.py b/lcls_tools/common/data_analysis/fit/method_base.py new file mode 100644 index 00000000..d686a082 --- /dev/null +++ b/lcls_tools/common/data_analysis/fit/method_base.py @@ -0,0 +1,148 @@ +from abc import ABC, abstractmethod +import numpy as np + +from pydantic import BaseModel, ConfigDict +from scipy.stats import rv_continuous + + +class Parameter(BaseModel): + model_config = ConfigDict(arbitrary_types_allowed=True) + bounds: list + initial_value: float = None + prior: rv_continuous = None + + +class ModelParameters(BaseModel): + model_config = ConfigDict(arbitrary_types_allowed=True) + name: str + parameters: dict[str, Parameter] + + @property + def bounds(self): + return np.vstack( + [ + np.array(parameter.bounds) + for parameter in self.parameters.values() + ] + ) + + @property + def initial_values(self): + return np.array( + [parameter.initial_value + for parameter in self.parameters.values()] + ) + + @initial_values.setter + def initial_values(self, initial_values: dict[str, float]): + for parameter, initial_value in initial_values.items(): + self.parameters[parameter].initial_value = initial_value + + @property + def priors(self): + return np.array( + [self.parameters[parameter].prior for parameter in self.parameters] + ) + + @priors.setter + def priors(self, priors: dict[str, float]): + for parameter, prior in priors.items(): + self.parameters[parameter].prior = prior + + +# TODO: define properties + +class MethodBase(ABC): + """ + Base abstract class for all fit methods, which serves as the bare minimum + skeleton code needed. Should be used only as a parent class to all method + models. + --------------------------- + Arguments: + param_names: list (list all of param names that the model will contain) + param_guesses: np.ndarray (array that contains a guess as to what + each param value is organized in the same order as param_names) + param_bounds: np.ndarray (array that contains the lower + and upper bound on for acceptable values of each parameter) + """ + + parameters: ModelParameters = None + + @abstractmethod + def find_init_values(self) -> list: + ... + + @abstractmethod + def find_priors(self, data: np.ndarray) -> dict: + ... + + def forward( + self, x: np.ndarray, method_parameter_dict: dict[str, float] + ) -> np.ndarray: + method_parameter_list = np.array( + [ + method_parameter_dict[parameter_name] + for parameter_name in self.parameters.parameters + ] + ) + return self._forward(x, method_parameter_list) + + @staticmethod + @abstractmethod + def _forward(x: np.ndarray, params: np.ndarray) -> np.ndarray: + ... + + def log_prior(self, method_parameter_dict: dict[str, rv_continuous]): + method_parameter_list = np.array( + [ + method_parameter_dict[parameter_name] + for parameter_name in self.parameters.parameters + ] + ) + return self._log_prior(method_parameter_list) + + @abstractmethod + def _log_prior(self, params: np.ndarray): + ... + + def log_likelihood(self, x: np.ndarray, y: np.ndarray, + method_parameter_dict: dict): + method_parameter_list = np.array( + [ + method_parameter_dict[parameter_name] + for parameter_name in self.parameters.parameters + ] + ) + return self._log_likelihood(x, y, method_parameter_list) + + def _log_likelihood( + self, x: np.ndarray, y: np.ndarray, method_parameter_list: np.ndarray + ): + return -np.sum((y - self._forward(x, method_parameter_list)) ** 2) + + def loss( + self, + method_parameter_list: np.ndarray, + x: np.ndarray, + y: np.ndarray, + use_priors: bool = False, + ): + loss_temp = -self._log_likelihood(x, y, method_parameter_list) + if use_priors: + loss_temp = loss_temp - self._log_prior(method_parameter_list) + return loss_temp + + @property + def profile_data(self): + """1D array typically projection data""" + return self._profile_data + + @profile_data.setter + def profile_data(self, profile_data): + if not isinstance(profile_data, np.ndarray): + raise TypeError("Input must be ndarray") + self._profile_data = profile_data + + self.find_init_values() + self.find_priors() + self.fitted_params_dict = {} diff --git a/lcls_tools/common/data_analysis/fit/methods.py b/lcls_tools/common/data_analysis/fit/methods.py new file mode 100644 index 00000000..9541e42b --- /dev/null +++ b/lcls_tools/common/data_analysis/fit/methods.py @@ -0,0 +1,103 @@ +import numpy as np +from scipy.stats import norm, gamma +from lcls_tools.common.data_analysis.fit.method_base import ( + MethodBase, + ModelParameters, + Parameter, +) + + +gaussian_parameters = ModelParameters( + name="Gaussian Parameters", + parameters={ + "mean": Parameter(bounds=[0.01, 1.0]), + "sigma": Parameter(bounds=[0.01, 5.0]), + "amplitude": Parameter(bounds=[0.01, 1.0]), + "offset": Parameter(bounds=[0.01, 1.0]), + }, +) + + +class GaussianModel(MethodBase): + """ + GaussianModel Class that finds initial parameter values for gaussian + distribution and builds probability density functions for the + likelyhood a parameter to be that value based on those initial + parameter values. Passing this class the variable profile_data + automatically updates the initial values and and probability + density functions to match that data. + """ + + parameters = gaussian_parameters + + def find_init_values(self) -> dict: + """Fit data without optimization, return values.""" + + data = self._profile_data + x = np.linspace(0, 1, len(data)) + offset = data.min() + 0.01 + amplitude = data.max() - offset + + weighted_mean = np.average(x, weights=data) + weighted_sigma = np.sqrt(np.cov(x, aweights=data)) + + init_values = { + "mean": weighted_mean, + "sigma": weighted_sigma, + "amplitude": amplitude, + "offset": offset, + } + self.parameters.initial_values = init_values + return init_values + + def find_priors(self) -> dict: + """ + Do initial guesses based on data and make distribution from that guess. + """ + # TODO: profile data setter method cals this and find values. + # should remove call to find values in that method. + # but since priors aren't supposed to be included (i think?) + # in this PR I will not update this. + init_values = self.find_init_values() + amplitude_mean = init_values["amplitude"] + amplitude_var = 0.05 + amplitude_alpha = (amplitude_mean**2) / amplitude_var + amplitude_beta = amplitude_mean / amplitude_var + amplitude_prior = gamma(amplitude_alpha, loc=0, + scale=1 / amplitude_beta) + + # Creating a normal distribution of points around the inital mean. + mean_prior = norm(init_values["mean"], 0.1) + sigma_alpha = 2.5 + sigma_beta = 5.0 + sigma_prior = gamma(sigma_alpha, loc=0, scale=1 / sigma_beta) + + # Creating a normal distribution of points around initial offset. + offset_prior = norm(init_values["offset"], 0.5) + priors = { + "mean": mean_prior, + "sigma": sigma_prior, + "amplitude": amplitude_prior, + "offset": offset_prior, + } + + self.parameters.priors = priors + return priors + + def _forward(self, x: np.array, method_parameter_list: np.array): + # Load distribution parameters + # needs to be array for scipy.minimize + mean = method_parameter_list[0] + sigma = method_parameter_list[1] + amplitude = method_parameter_list[2] + offset = method_parameter_list[3] + return (np.sqrt(2 * np.pi) * amplitude + * norm.pdf((x - mean) / sigma) + offset) + + def _log_prior(self, method_parameter_list: np.ndarray) -> float: + return np.sum( + [ + prior.logpdf(method_parameter_list[i]) + for i, prior in enumerate(self.parameters.priors) + ] + ) diff --git a/lcls_tools/common/data_analysis/fit/projection.py b/lcls_tools/common/data_analysis/fit/projection.py new file mode 100644 index 00000000..4b1d2907 --- /dev/null +++ b/lcls_tools/common/data_analysis/fit/projection.py @@ -0,0 +1,108 @@ +import numpy as np +import scipy.optimize +import scipy.signal +from pydantic import BaseModel, ConfigDict +from lcls_tools.common.data_analysis.fit.method_base import MethodBase +from lcls_tools.common.data_analysis.fit.methods import GaussianModel + + +class ProjectionFit(BaseModel): + """ + 1d fitting class that allows users to choose the model with which the fit + is performed, and if prior assumptions (bayesian regression) about + the data should be used when performing the fit. + Additionally there is an option to visualize the fitted data and priors. + -To perform a 1d fit, call fit_projection(projection_data={*data_to_fit*}) + ------------------------ + Arguments: + model: MethodBase (this argument is a child class object of method base + e.g GaussianModel & DoubleGaussianModel) + visualize_priors: bool (shows plots of the priors and init guess + distribution before fit) + use_priors: bool (incorporates prior distribution information into fit) + visualize_fit: bool (visualize the parameters as a function of the + forward function + from our model compared to distribution data) + """ + + # TODO: come up with better name + model_config = ConfigDict(arbitrary_types_allowed=True) + model: MethodBase = GaussianModel() + use_priors: bool = False + + def normalize(self, data: np.ndarray) -> np.ndarray: + """ + Normalize a 1d array by scaling and shifting data + s.t. data is between 0 and 1 + """ + data_copy = data.copy() + normalized_data = ((data_copy - np.min(data)) + / (np.max(data) - np.min(data))) + return normalized_data + + def unnormalize_model_params( + self, method_params_dict: dict, projection_data: np.ndarray + ) -> np.ndarray: + """ + Takes fitted and normalized params and returns them + to unnormalized values i.e the true fitted values of the distribution + """ + + projection_data_range = (np.max(projection_data) + - np.min(projection_data)) + length = len(projection_data) + for key, val in method_params_dict.items(): + if "sigma" in key or "mean" in key: + true_fitted_val = val * length + elif "offset" in key: + true_fitted_val = (val * projection_data_range + + np.min(projection_data)) + else: + true_fitted_val = val * projection_data_range + temp = {key: true_fitted_val} + method_params_dict.update(temp) + return method_params_dict + + def model_setup(self, projection_data=np.ndarray) -> None: + """sets up the model and init_values/priors""" + self.model.profile_data = projection_data + + def fit_model(self) -> scipy.optimize._optimize.OptimizeResult: + """ + Fits model params to distribution data and plots the fitted params + as a function of the model. + Returns optimizeResult object + """ + x = np.linspace(0, 1, len(self.model.profile_data)) + y = self.model.profile_data + + init_values = self.model.parameters.initial_values + bounds = self.model.parameters.bounds + res = scipy.optimize.minimize( + self.model.loss, + init_values, + args=(x, y, self.use_priors), + bounds=bounds, + method="Powell", + ) + return res + + def fit_projection(self, projection_data: np.ndarray) -> dict: + """ + type is dict[str, float] + Wrapper function that does all necessary steps to fit 1d array. + Returns a dictionary where the keys are the model params and their + values are the params fitted to the data + """ + assert len(projection_data.shape) == 1 + fitted_params_dict = {} + normalized_data = self.normalize(projection_data) + self.model_setup(projection_data=normalized_data) + res = self.fit_model() + + for i, param in enumerate(self.model.parameters.parameters): + fitted_params_dict[param] = (res.x)[i] + self.model.fitted_params_dict = fitted_params_dict.copy() + params_dict = self.unnormalize_model_params(fitted_params_dict, + projection_data) + return params_dict diff --git a/tests/datasets/fit/make_test_guassian.py b/tests/datasets/fit/make_test_guassian.py new file mode 100644 index 00000000..7ad4a63e --- /dev/null +++ b/tests/datasets/fit/make_test_guassian.py @@ -0,0 +1,15 @@ +import numpy as np +from scipy.stats import norm + +# import matplotlib.pyplot as plt + +# Generating 1000 point gaussian distribution +data = np.random.normal(size=1000) +fit = norm.pdf(data) + +# Plotting distribution +# plt.plot(data,fit, '.') +# plt.show() + +# Saving distribution +np.save("test_gaussian.npy", data) diff --git a/tests/datasets/fit/test_gaussian.npy b/tests/datasets/fit/test_gaussian.npy new file mode 100644 index 0000000000000000000000000000000000000000..a2026819e6d878dd1b7e407bbb7a9f7040a34d53 GIT binary patch literal 8128 zcmbVR_dgZh|2IOVl$KGXWMo9jTT#7kX&8}|5{i(L7q;6X{nA=OKib);2b3{r` zOzMuUy`%k28&g~R+ZO+Ke)gu7gT?B(!@Zk!7OQ*dgNllZa{J}P+{OOCYkl9pjm!ly z38IR`6Q+Wr!1hfMJ=lf;GtJ@Afgbbl;>?MbJkun|VJ5#w@m#=-@zS~aJ{7Qb&`_d2 zwijnioMkV^76XrVlZnv4Jan4ytO?y-i9PnTb`7g!V6&7RJe^+&x~}W@hkTvErG`x8 zGh-+24QeV(n#h9MQ*KcTh9l^kQ2f(TF$w%LTbsG4-%)w`Pg7+U4eviAKicijz!t$L zB9G3;psEm$>ambd7|!Wb^^GSHU8>~vJt0+rcIAeT&aeJL@EI8e-d%K@*dt(+Va-PT zVyh}7apE7`@BQGkC8q;wJk_3|(G;H3?p?nay##9wiq`6%ZwKnO5~q-i5tvi|#xN>O zMe(_|uV$Bq;qlNG^}X&>*jzg1l~1VyUEQ*m0tpjf*(WHO)G&<1do9oQ=u@C>I&t^D zmU667({u2$tiW>%yzv%0I0)V>&CWLYF&yJ=sFXTc0O3au3ss~KLJQ9m7p0mx;Fbz* zai;tM!kqtlKQawPt1kU!?wQ84s?E=q%$Q)$AJrsfTLy0@%3GF2NEozb+bgGu64X~a zupm6%1MGa}K8$@0z-g$@VNL8s^G%+1kzwVqQmD>BzB-8NYHtPhziELe_V9#{cc#!K z-k~9Fuom9a8veWciHPZfWb$$1Pk7YlMv&Y$3nMA-4>jk`;3$)k@Gg7^RKms$6CPxs z!gZ<1(34%b|Ce3v`rZokrVsU`g)hPDz1nr*hS|`uc8yQGe>NDn4J$UT;$C@oMEM&> z3ySBMdsRzEgO7veNkMJ~&S!YW&NTHv-SRt8*L@ts8GTxRu^=736u75O#Io>}S-=}E zTN*sp37$UJG>RkweQW#YCD2w92`=8jPBh(=ru$W?1xOM$e0F~+xG{;o^}&$__$2Pb zr%9Sa3!b&gQxy!5)r)yB6AtC#>L81IV?f9eO zYrS7?KB~Dzo!xY49vaR*BgosZ@Q54d(FCh3WXh4OzNC%7t}&S{Hy(~dekt3)*Y#{f z?Z1SrqoXWzNk2SV&zi#B3za@c^4SS{U;S=VqW(r5!5;_0kOGJ8Q-371k}){&_Qca$ z6ljX!($cbAMvt()zvwWE-h30G?fQSwA#oy~rb5RO%fQBZ&vD%B_jWsHc0HUL>~V0i z9mWJM$u&dKoP-LEt)VFmqu94c;{AsT8tOFYz23#u1oH+(-_Kt}OpZM^A996>m$wE4 zh|jF{hvRA1VHyin606q){YS?2n`-ZNcD5m-_?fEU3LSrZPL91YSpZP@O*0+DjWzz~+|)!H0s+2t#i<6xawa_V-D*vSuLV zc%|pot~v-QChyNYK?R+q$uxONBHW%lD|~)w5L+ql-6O;QAojU@lF?)%awmEh1QDq~ z`pT|bR>goXk)veaYtxvoJRZ&H9|ps|B@15{c7n6^n>1QbFS>S|&70?(!*MIy?B9~Z zaQmRa=C$hoaJZu2#9D)Kl#+WlTYP;G%^P()BX|eUCeD8E{-Pe7^GgYRerg7yt-shf z9E!m#rH8FEPuK{H)PHHe-~C1380kriF&^Si`^k!AuMCX2Aas2xVjjZ3|04b8#YGUi zQEtfOV1m`B9HAJ=A`H||HK{RXC)5VCQB4TtD06&4G!|&cug(8zpmYJpUKH297op>e zil$3fUNqy8E}6{>?`bgA*=U^>P>V-SG{lbTvCx`2@@v?NgJAY>N>=bb4V6-+twSQ` zu~Z_8@4xX`+`IddyU|=4+|c&oDSuDFLsj8nTRp}wU(@Fso7P|8^|S4r&z}R)#NfP@ z{cSK`x}-Sti-GQsM`G@@Wa3|T_p^snSHLHv#Wu*k2^38QwVwG+!upS&1wL8zpk0m$ zCAOXgmJvaJN?%N4hlj@+`SJ!>6Stihr>CvtZl_0yAE<`|aL= zhtYDn&zA%+=lW^tAHYs1UfA8W=5z+q!lUW}pR(|*l6=MqpGa)Vr7pkK8%4qMG;ghC zIv$zJy%K0Mg?Ey_ZlUy#;)|D!QmJ1T(X--AT%pngByc$dwXd_ssylI&-1CUd6pJHU ztr~G{%YWUWZcMO_)oXC(AIB3GAJ(O6_Ca=2=0bjK8`M5CZl1?~5U23ajPE)LYA&w3 zw$wWXHnsYr-)mUdzC=-|@~uDtB5~Va{S2@~CH=PrFHp5Ck8sx`1_q48)t;zDLY$1q z$c5G_Y;)b&9^n{=5nUIh!`DwkwX~hK>Lnh6wY*+Qqt_H3+2KPuFHc7^*`=Sudb7ZK zufL^WQwz$9eLPv1)eU32gLpIJ%0OON(p=)@B>d2}3b-tkkDDc#8&Bs9gWZg?V$pxA z^*wx>ayDrWqGvbscYkk&N6GSO$Jv`eSpC?*2H|qJ;_9kKWe&hq*Zl`Ytk{UmSdXN; z`+kF|v>dB89f5zRR8CSwDoPZ{#j-oF6GE51h+jD94J`t!z}K7TunOO=SaLU>3OSqf zd!+#nSY_DfJf4G%%kpu2Ogb#O+Yw^-)1Xe!JumV#4D|@a}FoFe7 z-44k#>c@cO$6r?7=6PTp@ouyF2O60CTs-qIr4v#r>ThL~<>03;%5vRTrjhn~?eC2} zO{g`?YZmI=1@EPkGJYq_!5#Ah&{H^u&#iT+Q=y{}mQtv|L(W8tkjqt48c`6Cc2rEU zWfr5j3#MasRAF(cX0z`g5w%OT2d`Q25DxZvH;QiWM>Tsx%cP%g@be|NLN}vDXow*# zm+>v*?HHlC43ibi+16KHU@`zldv=dL2xvf+TdrAQ*3(cXp|U^V_Y4l4_jfI_Xo45K zlh@65_QUr7q^mQj4LIg)#pZsy4Xm0>`yWp?Z8)G70=`+!|M(M3t#w^yFDF0SsbwD*(rf?4| zV!+a&Ed^tFsPpPozzz*2#!LP8`nCWIE2mZS=5mKnB~W%pUtlZJ<7Ag4(^zQgbHnNw z`xtg67(~D0nFQxpo21&XBzW0LbV$&m;{A0Am(7;OQS|fKJvL=opt4n9u&MVKC|!O? zHsR@oP@&L*pWIoPVs(Z^{MdxP@Ak{*{p*C6{js6L0%NGda^f-z8Gs;-+(X$bW6+pU zw#8wv1K5?rdLgI=HFodF?fo$fM*2)cQ)?Co??Xg1=%yl91YZx(*nhAf>F!N%X&lK9H zKePVV{1@YD6xVk*SHRQ1FTHDmTX2{8L(O2rMSK^=qg=YV2Nc-0GGDw-MUEj&4c}9x zc;aKr+v0sxSO~Ct>wmfvEr+?ag0pMUYR`G|5?K-qtz*71IuoT;2{|3m19XwF4q!r?Z|?|!7U*JTiA+<)!z^Pa|ITAr#-xeT%i&Yv$MUPO#I6AE^6fE6z=`dvDrj~hEf*wouFNX;dLG=DU-tx zUUj9lwXX`6`vwj3?({->ztFJfIvNz2PiFe74?);=$+&H|vw%aedm!In03-WXy0wi5 zQBuXqQndLWM%6vAm~fl~BY#Oxrc( zb_oeG51ssP4T8mv(|hB#HA18`S8Lo}!~x~VT>Z2WOpbCAT4OeX#{yZViK?sjn|BvV zYVU*Xa+hZP!WU75sVMg;bk+Y-OFZMthcQ@NQ@mxo30u^ZpKKo`!M%p!vk&$#@W_{v zUfYXAG~{DO6EquP{NwwW19FoX6VG?%F0&msi<}Y^%%y|f2UVj@BCTjxzxRxo1|8ET zJC*MKt_8o|?FK)a7IC~LtU9{B2VAGU|H1PvY_)uOI&*g=+$Xl4cC2a!54(NgOBRDr zXrxlmdVvl#CGI~W9RngtR0dx-FU~_eA!@20h%?hR~pklOJL$9h0Skzht zQOhA{Cshd?wH${S2^XfoYJS`xa_#L+nJS!oeaPV9wMv-KJyacji;lH}=KU&7^O#om zRq*6$-!Q1|vwCl5v15zP>%x(INPPFpFV}4h4)0i8uX}(5{`$I2n=+<=&!V*~hB=D| z{GRO6)oDeA)2XZrj~OrQ55mO9eQrF)oyeYh(N^GQ{i=`LI4}@b z50jj;Zv9=Yb6k5BHbt*?6)Dzph)Kk$-g>_Gsp%xaPC)+YUlq z#g9?=C^o)RabO7EsU;8E1tyJRJGCMOUE#4Z9L8i;fM#WcOUHYiSW?3+|n78yd!q z+fou=6w&~rcC%zIG^55kkKLhm?YL%}@#ug#9Y6JLN;BBe0mVm8DnGhL0!b?~TeS}i z*zib)bvC*ItVdhE^h;J@${B-$#-9gqRV4L?biac<+YMpa4QvE{vkPx>&aogUzI~`; zxD<{;sIT2B&*{8*{E7YW1nkrs2~B64gem!g1+_XN#+!VBu8}&dyl{feJ8uYl{k9tx zB~4=1tb75YIka*u+3cFW^d{E|IXqv_N;acr}if`Z!GTvlE5V`OS)Fil9zF?NUX332G?q zm*w|L2lWSjL0380iT(u)u@euQaeBUQ1JfuS^IGWUgIVS1JH~Y;(T{^Dw`_5fKwia_ zOZThOFAB_`ZoE>9G?c`G5_U7exmCzye|AP<+g*4Ooz0&X9@ARNpQ5R z@dycNYrY8bmbYVtaPlcB%QSenMa|(z%M7NAl&Jt8k(3;Lm(j3*5hU>=R{rI9i5y$dexHndt?V~tyD2d9 z=bsl%XBzGom8D!)qCi{8&K-Kc=HS=SPRaZStMiFAFNb{oJFs&@Ny0vvVc^NxOuuWI zg*M|9HdW3^2nlTYdGqcvj_wT4k4kKSO?$!YH2Y6X=l@fSHjC~D;|CdIhorkeqUk5;!)prOJh9~Ma+C%Srm3_G61|`vRLWoYjSiqZ3(0Kc(*~p+3Of@DY^Bij-Gi>rgiWUil)>OwS%5f5T z)5|Bajk4gb;YF2T8y586$P}nMK7i-PWBy5a4dJ<$CN740Byc!x`a*kvgoW2nz}S5r z;>P*lrP%OQ9Dk~R%>DWYM2&M2gTK#XzR0)j>3`=@#XqLTubBa29|}%&ohn1^`d#?jEf&P#>z zS~lQ_w)XMK2?pfVI!a9J8^W0V8J`veSK!RrfW4Q4%OPt>S@)(!6G)dA7#^&kW4-9? znGo_62)Vi+^iX6c8kdqb9j6V48%p!RjR@KUugtTjl$%;$En&6qUO9 zxDA~`3Onc<>d~om;44Q76?;wD#q#X?@k#Z^spUswaR1q1HZd{-AJ@oyi}y}Piq14Gse;!9Fw0l`s54{vaHlp#Abk4qGPbAG7lNE|Glsi8%Cy3 z@^p7u-|9T;4O6*p9@{IuiB-gQ;JD^a{a`c&io#D*v;ryk$E;)CXI1}h+n41>Tuuf{ zY4<>%KXZ74Z%3=H+!)-9(70A!TmmOVWHwH*M)AI4%f8$1hA=FC-#U(BWx$tluS$O= zV*Q4IjuNAOOys&L+*UP$A)XV`Je8y9^(G971Pba8N-Smmm;oxU@pk^eMbKm3i*A+c zhhr1u&D+nA;pO(IV_vI$5YoOf(A3BTnWA(?u;Lsf<(UeGCDa0Tl$_eqIs=#%m3;lo zFh2A!D9>^3#l^02^~xwR9>(80A4FFm-|g%gY5pqQS(lVV-dK-%BR4Ef6X(!S!fPwl zW&}R^#m+7Mo`c4A|4iz)S^QyetGYS280Lh_yNk(7n4T0jW%x7^x0fD07E&_|fvS;T zAHE-gI}LU|^h_$;$_|fF{4@jc5q&qbtcd7#w*IPT%n1Hb2*}juCu7uerlRrk5QZP> zIUm?Hj4Yi!@Ze?(e*3qp!ftvLaydZ!0ed<=J6TDynj>R!V^$LSE&`3LUt1-^PTUzS zzhz@;BQmR_e|sNZL7gu~NJwO%zR-b_6s%G)e<_oF78rer{l{?Tgq>HCXm)R^t&%L0a?z(k!^4ogEW55EBnuY+r{%*VIw(k z@j<^&<7XnS4D)tu6d2B}&oSWwUB&yYp+p{$Cw46@e4F)ZKa=1jcCUFVC|vQ z$|mg3{9acn47n9$}}!koQdg3EYU`MvYC zHp`|E2fnBu^uN`Dc3FJInI>eII(LWa{(1_`{GL}v-s#1CCFw>Ku~8uJ*=ZA{PQ{HO zvqgLCN>I_fq5q;I0}re*7MplJ4+^Gbry70!LaN_wO{!BTnpl#{E-5zS@8;sXgvB{H z7y9yyt`!X>a>gFN{xu4z1*K7PgdW^5?0Sm0xdqPGvC`U2=TLII&QkPC2UI474Vtf% z;HdiUnwo(>F@3uQ_;TWYrZk}& z51%_G{ybyVpPD~;s=B@wQaPQP?e*z!@cuEC=}r<1?|hP%Uoef)LYYm;{VgCTxl*h0 zr4c+^_@09K0%~`sa3`8gW0LaEseco*AivOl+C8%cH=W(*m-{Im9n~bIOADKzKt?T& zJ8>9R3^;gmd&^ON@!|Y+Ln2TVhYhyPa}gAkg!;wLOv6IL&V0W;1Gv#QaZWjG2&UbL zURtpu_(@pvY30lc^5!;u`c~6{N=nX-Tztte%-Q=e#}4s7f64@RRX^mq3B8gCBS9&V zaXX=X33Yy13MYIo!mM%?rRYOl5MlXA_Mqh&;24X)T$;uMFH19+Xw$GK!l!Cu**3@4(ntqf_!u zqbP4L#?>>PhMDb2Z7?%}r8Zv5wEJ^d*AQ_3A*BaieODYg;4%fmr@nr;xUUr&>9u`P zh71(Fzn8I|&kZ*_@f>z?n!q>q9Q6FaQdp~bh0(xG2iJ#EQCqJ}0GornQj@k9T92`9 za`$AQ^%sW}X{8Ft>AAIMmwOwEo|hS5I8DHoil7ynB*rR_)0H>0(m}H5Idg+k8-9K_ z_ifkb1#t5cQ+FqIfTKFE`-XM1Xv^QdIpiTHaX2w=Q=i9FsM)PK1 z+zU3n#*0%}8FOj@ZCS!Ub&h#PHCHD)~@;m%3z~3lapvP*VQ1znuS;0TD+e+ z8xh@nJlp!lfgCg8AQRk)zqT$p>HpV^i{^Ur!SBhCIk`gt4l^)Y=;gJT;|%18r+*vT zK7scGo}WlIq2bNrqJxjKTd|ktsOvk9C9R*3x4U76~;k%D;Ov6kQT=cworbyrqI3H>0{IG5uG|y#RrJiNrlEV+q z_9xS;da+gibl(^#y81uY-nlMKhDWc5HDV7I0Qbcsw~04~;m)7zYe(36 z!RRLU&Gq{Sko8}5JB6176SqWaH9pcYYI{a;abqWzK0nXI5A$H9U(#b3xSGfB&YZ05 z7=>eD_7Z2WwgYKLm-Jo7S-AI3T+Mle4u_o0pPtL?!zp96e}C)BU{6I#_|ut2*rOCF odZ(8TTiz) None: + self.data_location = "./tests/datasets/fit/" + self.data_filename = os.path.join(self.data_location, "test_gaussian.npy") + self.data = np.load(self.data_filename) + self.gaussian_model = GaussianModel() + self.gaussian_model.profile_data = self.data + # Decimal place to check calculation errors to + self.decimals = 2 + return super().setUp() + + @unittest.skip + def test_find_init_values(self): + init_dict = self.gaussian_model.find_init_values() + self.assertIsNotNone(init_dict) + # TODO: pick up test here, the init dict is not as expected + self.assertAlmostEqual(init_dict["mean"], -0.07874, places=self.decimals) + self.assertAlmostEqual(init_dict["sigma"], 1.0, places=self.decimals) + return + + @unittest.skip + def test_find_priors(self): + priors_dict = self.gaussian_model.find_priors() + self.assertIsNotNone(priors_dict) + return # TODO: flesh out this test + + @unittest.skip + def test_forward(self): + init_dict = self.gaussian_model.find_init_values() + fit = self.gaussian_model.forward(self.data, init_dict) + # TODO: better way to check this? w/o hardcode? + self.assertAlmostEqual(fit.max(), 0.4, places=self.decimals) + self.assertAlmostEqual(fit.min(), 0.0, places=self.decimals) + return + + # def super_gaussian(self, x, amp, mu, sig, P, offset): + # """Super Gaussian Function""" + # """Degree of P related to flatness of curve at peak""" + # return amp * np.exp((-abs(x - mu) ** (P)) / (2 * sig ** (P))) + offset + # def double_gaussian(self, x, amp, mu, sig, amp2, nu, rho, offset): + # return ( + # amp * np.exp(-np.power(x - mu, 2.0) / (2 * np.power(sig, 2.0))) + # + amp2 * np.exp(-np.power(x - nu, 2.0) / (2 * np.power(rho, 2.0))) + # ) + offset + # @unittest.skip("Assertion is not raised when it should be; fixing in issue #130.") + # def test_fit_tool_gaussian(self): + # # Test that the fitting tool can fit each type of gaussian distribution + # x_data = np.arange(500) + # # generated data for pure gaussian + # test_params = [3, 125, 45, 1.5] + # y_data = self.gaussian(x_data, *test_params) + # y_noise = np.random.normal(size=len(x_data), scale=0.04) + # y_test = y_data + y_noise + # fitting_tool = FittingTool(data=y_test) + # fits = fitting_tool.get_fit() + # self.assertIsInstance(fits, dict) + # for key, val in fits.items(): + # self.assertIsInstance(val, dict) + # self.assertIn("rmse", val) + # self.assertLessEqual(val["rmse"], 0.4) + # @unittest.skip("Assertion is not raised when it should be; fixing in issue #130.") + # def test_fit_tool_super_gaussian(self): + # x_data = np.arange(500) + # y_noise = np.random.normal(size=len(x_data), scale=0.04) + # test_params_super_gauss = [4, 215, 75, 4, 1] + # y_data_super_gauss = self.super_gaussian(x_data, *test_params_super_gauss) + # y_test_super_gauss = y_data_super_gauss + y_noise + # super_gauss_fitting_tool = FittingTool(data=y_test_super_gauss) + # s_fits = super_gauss_fitting_tool.get_fit() + # self.assertIsInstance(s_fits, dict) + # for key, val in s_fits.items(): + # self.assertIsInstance(val, dict) + # self.assertIn("rmse", val) + # if key == "super_gaussian": + # self.assertLessEqual(val["rmse"], 0.4) + # @unittest.skip("Assertion is not raised when it should be; fixing in issue #130.") + # def test_fit_tool_double_gaussian(self): + # # generated data for super gaussian + # x_data = np.arange(500) + # y_noise = np.random.normal(size=len(x_data), scale=0.04) + # test_params_double_gauss = [2, 100, 25, 10, 240, 25, 2] + # y_data_double_gauss = self.double_gaussian(x_data, *test_params_double_gauss) + # y_test_double_gauss = y_data_double_gauss + 3 * y_noise + # double_gauss_fitting_tool = FittingTool(data=y_test_double_gauss) + # d_fits = double_gauss_fitting_tool.get_fit() + # self.assertIsInstance(d_fits, dict) + # for key, val in d_fits.items(): + # self.assertIsInstance(val, dict) + # self.assertIn("rmse", val) + # if key == "double_gaussian": + # self.assertLessEqual(val["rmse"], 0.8) diff --git a/tests/unit_tests/lcls_tools/common/data_analysis/test_fit_gauss.py b/tests/unit_tests/lcls_tools/common/data_analysis/test_fit_gauss.py deleted file mode 100644 index 028f250b..00000000 --- a/tests/unit_tests/lcls_tools/common/data_analysis/test_fit_gauss.py +++ /dev/null @@ -1,71 +0,0 @@ -from lcls_tools.common.data_analysis.fitting_tool import FittingTool -import numpy as np -import unittest - - -class TestFitTool(unittest.TestCase): - def gaussian(self, x, amp, mu, sig, offset): - return amp * np.exp(-np.power(x - mu, 2.0) / (2 * np.power(sig, 2.0))) + offset - - def super_gaussian(self, x, amp, mu, sig, P, offset): - """Super Gaussian Function""" - """Degree of P related to flatness of curve at peak""" - return amp * np.exp((-abs(x - mu) ** (P)) / (2 * sig ** (P))) + offset - - def double_gaussian(self, x, amp, mu, sig, amp2, nu, rho, offset): - return ( - amp * np.exp(-np.power(x - mu, 2.0) / (2 * np.power(sig, 2.0))) - + amp2 * np.exp(-np.power(x - nu, 2.0) / (2 * np.power(rho, 2.0))) - ) + offset - - @unittest.skip("Assertion is not raised when it should be; fixing in issue #130.") - def test_fit_tool_gaussian(self): - # Test that the fitting tool can fit each type of gaussian distribution - x_data = np.arange(500) - - # generated data for pure gaussian - test_params = [3, 125, 45, 1.5] - y_data = self.gaussian(x_data, *test_params) - y_noise = np.random.normal(size=len(x_data), scale=0.04) - y_test = y_data + y_noise - - fitting_tool = FittingTool(data=y_test) - fits = fitting_tool.get_fit() - self.assertIsInstance(fits, dict) - for key, val in fits.items(): - self.assertIsInstance(val, dict) - self.assertIn("rmse", val) - self.assertLessEqual(val["rmse"], 0.4) - - @unittest.skip("Assertion is not raised when it should be; fixing in issue #130.") - def test_fit_tool_super_gaussian(self): - x_data = np.arange(500) - y_noise = np.random.normal(size=len(x_data), scale=0.04) - test_params_super_gauss = [4, 215, 75, 4, 1] - y_data_super_gauss = self.super_gaussian(x_data, *test_params_super_gauss) - y_test_super_gauss = y_data_super_gauss + y_noise - super_gauss_fitting_tool = FittingTool(data=y_test_super_gauss) - s_fits = super_gauss_fitting_tool.get_fit() - self.assertIsInstance(s_fits, dict) - for key, val in s_fits.items(): - self.assertIsInstance(val, dict) - self.assertIn("rmse", val) - if key == "super_gaussian": - self.assertLessEqual(val["rmse"], 0.4) - - @unittest.skip("Assertion is not raised when it should be; fixing in issue #130.") - def test_fit_tool_double_gaussian(self): - # generated data for super gaussian - x_data = np.arange(500) - y_noise = np.random.normal(size=len(x_data), scale=0.04) - test_params_double_gauss = [2, 100, 25, 10, 240, 25, 2] - y_data_double_gauss = self.double_gaussian(x_data, *test_params_double_gauss) - y_test_double_gauss = y_data_double_gauss + 3 * y_noise - double_gauss_fitting_tool = FittingTool(data=y_test_double_gauss) - d_fits = double_gauss_fitting_tool.get_fit() - self.assertIsInstance(d_fits, dict) - for key, val in d_fits.items(): - self.assertIsInstance(val, dict) - self.assertIn("rmse", val) - if key == "double_gaussian": - self.assertLessEqual(val["rmse"], 0.8)