Mercurial > repos > rv43 > tomo
comparison fit.py @ 49:26f99fdd8d61 draft
"planemo upload for repository https://github.com/rolfverberg/galaxytools commit 4f7738d02f4a3fd91373f43937ed311b6fe11a12"
author | rv43 |
---|---|
date | Thu, 28 Jul 2022 16:05:24 +0000 |
parents | |
children | 98a83f03d91b |
comparison
equal
deleted
inserted
replaced
48:059819ea1f0e | 49:26f99fdd8d61 |
---|---|
1 #!/usr/bin/env python3 | |
2 | |
3 # -*- coding: utf-8 -*- | |
4 """ | |
5 Created on Mon Dec 6 15:36:22 2021 | |
6 | |
7 @author: rv43 | |
8 """ | |
9 | |
10 import sys | |
11 import re | |
12 import logging | |
13 import numpy as np | |
14 | |
15 from asteval import Interpreter | |
16 from copy import deepcopy | |
17 #from lmfit import Minimizer | |
18 from lmfit import Model, Parameters | |
19 from lmfit.models import ConstantModel, LinearModel, QuadraticModel, PolynomialModel,\ | |
20 StepModel, RectangleModel, GaussianModel, LorentzianModel | |
21 | |
22 #from .general import * | |
23 from general import is_index, index_nearest, quickPlot | |
24 | |
25 # sigma = fwhm_factor*fwhm | |
26 fwhm_factor = { | |
27 'gaussian' : f'fwhm/(2*sqrt(2*log(2)))', | |
28 'lorentzian' : f'0.5*fwhm', | |
29 'splitlorentzian' : f'0.5*fwhm', # sigma = sigma_r | |
30 'voight' : f'0.2776*fwhm', # sigma = gamma | |
31 'pseudovoight' : f'0.5*fwhm'} # fraction = 0.5 | |
32 | |
33 # amplitude = height_factor*height*fwhm | |
34 height_factor = { | |
35 'gaussian' : f'height*fwhm*0.5*sqrt(pi/log(2))', | |
36 'lorentzian' : f'height*fwhm*0.5*pi', | |
37 'splitlorentzian' : f'height*fwhm*0.5*pi', # sigma = sigma_r | |
38 'voight' : f'3.334*height*fwhm', # sigma = gamma | |
39 'pseudovoight' : f'1.268*height*fwhm'} # fraction = 0.5 | |
40 | |
41 class Fit: | |
42 """Wrapper class for lmfit | |
43 """ | |
44 def __init__(self, x, y, models=None, **kwargs): | |
45 self._x = x | |
46 self._y = y | |
47 self._model = None | |
48 self._parameters = Parameters() | |
49 self._result = None | |
50 if models is not None: | |
51 if callable(models) or isinstance(models, str): | |
52 kwargs = self.add_model(models, **kwargs) | |
53 elif isinstance(models, (tuple, list)): | |
54 for model in models: | |
55 kwargs = self.add_model(model, **kwargs) | |
56 self.fit(**kwargs) | |
57 | |
58 @classmethod | |
59 def fit_data(cls, x, y, models, **kwargs): | |
60 return cls(x, y, models, **kwargs) | |
61 | |
62 @property | |
63 def best_errors(self): | |
64 if self._result is None: | |
65 return None | |
66 errors = {} | |
67 names = sorted(self._result.params) | |
68 for name in names: | |
69 par = self._result.params[name] | |
70 errors[name] = par.stderr | |
71 return errors | |
72 | |
73 @property | |
74 def best_fit(self): | |
75 if self._result is None: | |
76 return None | |
77 return self._result.best_fit | |
78 | |
79 @property | |
80 def best_parameters(self): | |
81 if self._result is None: | |
82 return None | |
83 parameters = [] | |
84 names = sorted(self._result.params) | |
85 for name in names: | |
86 par = self._result.params[name] | |
87 parameters.append({'name' : par.name, 'value' : par.value, 'error' : par.stderr, | |
88 'init_value' : par.init_value, 'min' : par.min, 'max' : par.max, | |
89 'vary' : par.vary, 'expr' : par.expr}) | |
90 return parameters | |
91 | |
92 @property | |
93 def best_values(self): | |
94 if self._result is None: | |
95 return None | |
96 return self._result.params.valuesdict() | |
97 | |
98 @property | |
99 def chisqr(self): | |
100 return self._result.chisqr | |
101 | |
102 @property | |
103 def covar(self): | |
104 return self._result.covar | |
105 | |
106 @property | |
107 def init_values(self): | |
108 if self._result is None: | |
109 return None | |
110 return self._result.init_values | |
111 | |
112 @property | |
113 def num_func_eval(self): | |
114 return self._result.nfev | |
115 | |
116 @property | |
117 def redchi(self): | |
118 return self._result.redchi | |
119 | |
120 @property | |
121 def residual(self): | |
122 return self._result.residual | |
123 | |
124 @property | |
125 def success(self): | |
126 if not self._result.success: | |
127 # print(f'ier = {self._result.ier}') | |
128 # print(f'lmdif_message = {self._result.lmdif_message}') | |
129 # print(f'message = {self._result.message}') | |
130 # print(f'nfev = {self._result.nfev}') | |
131 # print(f'redchi = {self._result.redchi}') | |
132 # print(f'success = {self._result.success}') | |
133 if self._result.ier == 0 or self._result.ier == 5: | |
134 logging.warning(f'ier = {self._result.ier}: {self._result.message}') | |
135 else: | |
136 logging.warning(f'ier = {self._result.ier}: {self._result.message}') | |
137 return True | |
138 # self.print_fit_report() | |
139 # self.plot() | |
140 return self._result.success | |
141 | |
142 @property | |
143 def var_names(self): | |
144 """Intended to be used with covar | |
145 """ | |
146 if self._result is None: | |
147 return None | |
148 return self._result.var_names | |
149 | |
150 def print_fit_report(self, show_correl=False): | |
151 if self._result is not None: | |
152 print(self._result.fit_report(show_correl=show_correl)) | |
153 | |
154 def add_parameter(self, **parameter): | |
155 if not isinstance(parameter, dict): | |
156 illegal_value(parameter, 'parameter', 'add_parameter') | |
157 return | |
158 self._parameters.add(**parameter) | |
159 | |
160 def add_model(self, model, prefix=None, parameters=None, **kwargs): | |
161 # Create the new model | |
162 # print('\nAt start adding model:') | |
163 # self._parameters.pretty_print() | |
164 if prefix is not None and not isinstance(prefix, str): | |
165 logging.warning('Ignoring illegal prefix: {model} {type(model)}') | |
166 prefix = None | |
167 if callable(model): | |
168 newmodel = Model(model, prefix=prefix) | |
169 elif isinstance(model, str): | |
170 if model == 'constant': | |
171 newmodel = ConstantModel(prefix=prefix) | |
172 elif model == 'linear': | |
173 newmodel = LinearModel(prefix=prefix) | |
174 elif model == 'quadratic': | |
175 newmodel = QuadraticModel(prefix=prefix) | |
176 elif model == 'gaussian': | |
177 newmodel = GaussianModel(prefix=prefix) | |
178 elif model == 'step': | |
179 form = kwargs.get('form') | |
180 if form is not None: | |
181 del kwargs['form'] | |
182 if form is None or form not in ('linear', 'atan', 'arctan', 'erf', 'logistic'): | |
183 logging.error(f'Illegal form parameter for build-in step model ({form})') | |
184 return kwargs | |
185 newmodel = StepModel(prefix=prefix, form=form) | |
186 elif model == 'rectangle': | |
187 form = kwargs.get('form') | |
188 if form is not None: | |
189 del kwargs['form'] | |
190 if form is None or form not in ('linear', 'atan', 'arctan', 'erf', 'logistic'): | |
191 logging.error(f'Illegal form parameter for build-in rectangle model ({form})') | |
192 return kwargs | |
193 newmodel = RectangleModel(prefix=prefix, form=form) | |
194 else: | |
195 logging.error('Unknown build-in fit model') | |
196 return kwargs | |
197 else: | |
198 illegal_value(model, 'model', 'add_model') | |
199 return kwargs | |
200 | |
201 # Add the new model to the current one | |
202 if self._model is None: | |
203 self._model = newmodel | |
204 else: | |
205 self._model += newmodel | |
206 if self._parameters is None: | |
207 self._parameters = newmodel.make_params() | |
208 else: | |
209 self._parameters += newmodel.make_params() | |
210 # print('\nAfter adding model:') | |
211 # self._parameters.pretty_print() | |
212 | |
213 # Initialize the model parameters | |
214 if prefix is None: | |
215 prefix = "" | |
216 if parameters is not None: | |
217 if not isinstance(parameters, (tuple, list)): | |
218 illegal_value(parameters, 'parameters', 'add_model') | |
219 return kwargs | |
220 for parameter in parameters: | |
221 if not isinstance(parameter, dict): | |
222 illegal_value(parameter, 'parameter in parameters', 'add_model') | |
223 return kwargs | |
224 parameter['name'] = prefix+parameter['name'] | |
225 self._parameters.add(**parameter) | |
226 for name, value in kwargs.items(): | |
227 if isinstance(value, (int, float)): | |
228 self._parameters.add(prefix+name, value=value) | |
229 # print('\nAt end add_model:') | |
230 # self._parameters.pretty_print() | |
231 | |
232 return kwargs | |
233 | |
234 def fit(self, interactive=False, guess=False, **kwargs): | |
235 if self._model is None: | |
236 logging.error('Undefined fit model') | |
237 return | |
238 # Current parameter values | |
239 pars = self._parameters.valuesdict() | |
240 # Apply parameter updates through keyword arguments | |
241 for par in set(pars) & set(kwargs): | |
242 pars[par] = kwargs.pop(par) | |
243 self._parameters[par].set(value=pars[par]) | |
244 # Check for uninitialized parameters | |
245 for par, value in pars.items(): | |
246 if value is None or np.isinf(value) or np.isnan(value): | |
247 if interactive: | |
248 self._parameters[par].set(value= | |
249 input_num(f'Enter an initial value for {par}: ')) | |
250 else: | |
251 self._parameters[par].set(value=1.0) | |
252 # print('\nAt start actual fit:') | |
253 # print(f'kwargs = {kwargs}') | |
254 # self._parameters.pretty_print() | |
255 # print(f'parameters:\n{self._parameters}') | |
256 # print(f'x = {self._x}') | |
257 # print(f'len(x) = {len(self._x)}') | |
258 # print(f'y = {self._y}') | |
259 # print(f'len(y) = {len(self._y)}') | |
260 if guess: | |
261 self._parameters = self._model.guess(self._y, x=self._x) | |
262 self._result = self._model.fit(self._y, self._parameters, x=self._x, **kwargs) | |
263 # print('\nAt end actual fit:') | |
264 # print(f'var_names:\n{self._result.var_names}') | |
265 # print(f'stderr:\n{np.sqrt(np.diagonal(self._result.covar))}') | |
266 # self._parameters.pretty_print() | |
267 # print(f'parameters:\n{self._parameters}') | |
268 | |
269 def plot(self): | |
270 if self._result is None: | |
271 return | |
272 components = self._result.eval_components() | |
273 plots = ((self._x, self._y, '.'), (self._x, self._result.best_fit, 'k-'), | |
274 (self._x, self._result.init_fit, 'g-')) | |
275 legend = ['data', 'best fit', 'init'] | |
276 if len(components) > 1: | |
277 for modelname, y in components.items(): | |
278 if isinstance(y, (int, float)): | |
279 y *= np.ones(len(self._x)) | |
280 plots += ((self._x, y, '--'),) | |
281 # if modelname[-1] == '_': | |
282 # legend.append(modelname[:-1]) | |
283 # else: | |
284 # legend.append(modelname) | |
285 quickPlot(plots, legend=legend, block=True) | |
286 | |
287 @staticmethod | |
288 def guess_init_peak(x, y, *args, center_guess=None, use_max_for_center=True): | |
289 """ Return a guess for the initial height, center and fwhm for a peak | |
290 """ | |
291 center_guesses = None | |
292 if len(x) != len(y): | |
293 logging.error(f'Illegal x and y lengths ({len(x)}, {len(y)}), skip initial guess') | |
294 return None, None, None | |
295 if isinstance(center_guess, (int, float)): | |
296 if len(args): | |
297 logging.warning('Ignoring additional arguments for single center_guess value') | |
298 elif isinstance(center_guess, (tuple, list, np.ndarray)): | |
299 if len(center_guess) == 1: | |
300 logging.warning('Ignoring additional arguments for single center_guess value') | |
301 if not isinstance(center_guess[0], (int, float)): | |
302 raise ValueError(f'Illegal center_guess type ({type(center_guess[0])})') | |
303 center_guess = center_guess[0] | |
304 else: | |
305 if len(args) != 1: | |
306 raise ValueError(f'Illegal number of arguments ({len(args)})') | |
307 n = args[0] | |
308 if not is_index(n, 0, len(center_guess)): | |
309 raise ValueError('Illegal argument') | |
310 center_guesses = center_guess | |
311 center_guess = center_guesses[n] | |
312 elif center_guess is not None: | |
313 raise ValueError(f'Illegal center_guess type ({type(center_guess)})') | |
314 | |
315 # Sort the inputs | |
316 index = np.argsort(x) | |
317 x = x[index] | |
318 y = y[index] | |
319 miny = y.min() | |
320 # print(f'miny = {miny}') | |
321 # print(f'x_range = {x[0]} {x[-1]} {len(x)}') | |
322 # print(f'y_range = {y[0]} {y[-1]} {len(y)}') | |
323 | |
324 # xx = x | |
325 # yy = y | |
326 # Set range for current peak | |
327 # print(f'center_guesses = {center_guesses}') | |
328 if center_guesses is not None: | |
329 if n == 0: | |
330 low = 0 | |
331 upp = index_nearest(x, (center_guesses[0]+center_guesses[1])/2) | |
332 elif n == len(center_guesses)-1: | |
333 low = index_nearest(x, (center_guesses[n-1]+center_guesses[n])/2) | |
334 upp = len(x) | |
335 else: | |
336 low = index_nearest(x, (center_guesses[n-1]+center_guesses[n])/2) | |
337 upp = index_nearest(x, (center_guesses[n]+center_guesses[n+1])/2) | |
338 # print(f'low = {low}') | |
339 # print(f'upp = {upp}') | |
340 x = x[low:upp] | |
341 y = y[low:upp] | |
342 # quickPlot(x, y, vlines=(x[0], center_guess, x[-1]), block=True) | |
343 | |
344 # Estimate FHHM | |
345 maxy = y.max() | |
346 # print(f'x_range = {x[0]} {x[-1]} {len(x)}') | |
347 # print(f'y_range = {y[0]} {y[-1]} {len(y)} {miny} {maxy}') | |
348 # print(f'center_guess = {center_guess}') | |
349 if center_guess is None: | |
350 center_index = np.argmax(y) | |
351 center = x[center_index] | |
352 height = maxy-miny | |
353 else: | |
354 if use_max_for_center: | |
355 center_index = np.argmax(y) | |
356 center = x[center_index] | |
357 if center_index < 0.1*len(x) or center_index > 0.9*len(x): | |
358 center_index = index_nearest(x, center_guess) | |
359 center = center_guess | |
360 else: | |
361 center_index = index_nearest(x, center_guess) | |
362 center = center_guess | |
363 height = y[center_index]-miny | |
364 # print(f'center_index = {center_index}') | |
365 # print(f'center = {center}') | |
366 # print(f'height = {height}') | |
367 half_height = miny+0.5*height | |
368 # print(f'half_height = {half_height}') | |
369 fwhm_index1 = 0 | |
370 for i in range(center_index, fwhm_index1, -1): | |
371 if y[i] < half_height: | |
372 fwhm_index1 = i | |
373 break | |
374 # print(f'fwhm_index1 = {fwhm_index1} {x[fwhm_index1]}') | |
375 fwhm_index2 = len(x)-1 | |
376 for i in range(center_index, fwhm_index2): | |
377 if y[i] < half_height: | |
378 fwhm_index2 = i | |
379 break | |
380 # print(f'fwhm_index2 = {fwhm_index2} {x[fwhm_index2]}') | |
381 # quickPlot((x,y,'o'), vlines=(x[fwhm_index1], center, x[fwhm_index2]), block=True) | |
382 if fwhm_index1 == 0 and fwhm_index2 < len(x)-1: | |
383 fwhm = 2*(x[fwhm_index2]-center) | |
384 elif fwhm_index1 > 0 and fwhm_index2 == len(x)-1: | |
385 fwhm = 2*(center-x[fwhm_index1]) | |
386 else: | |
387 fwhm = x[fwhm_index2]-x[fwhm_index1] | |
388 # print(f'fwhm_index1 = {fwhm_index1} {x[fwhm_index1]}') | |
389 # print(f'fwhm_index2 = {fwhm_index2} {x[fwhm_index2]}') | |
390 # print(f'fwhm = {fwhm}') | |
391 | |
392 # Return height, center and FWHM | |
393 # quickPlot((x,y,'o'), (xx,yy), vlines=(x[fwhm_index1], center, x[fwhm_index2]), block=True) | |
394 return height, center, fwhm | |
395 | |
396 | |
397 class FitMultipeak(Fit): | |
398 """Fit data with multiple peaks | |
399 """ | |
400 def __init__(self, x, y, normalize=True): | |
401 super().__init__(x, deepcopy(y)) | |
402 self._norm = None | |
403 self._fwhm_max = None | |
404 self._sigma_max = None | |
405 if normalize: | |
406 self._normalize() | |
407 #quickPlot((self._x,self._y), block=True) | |
408 | |
409 @classmethod | |
410 def fit_multipeak(cls, x, y, centers, peak_models='gaussian', center_exprs=None, fit_type=None, | |
411 background_order=None, fwhm_max=None, plot_components=None): | |
412 """Make sure that centers and fwhm_max are in the correct units and consistent with expr | |
413 for a uniform fit (fit_type == 'uniform') | |
414 """ | |
415 fit = cls(x, y) | |
416 success = fit.fit(centers, fit_type=fit_type, peak_models=peak_models, fwhm_max=fwhm_max, | |
417 center_exprs=center_exprs, background_order=background_order, | |
418 plot_components=plot_components) | |
419 if success: | |
420 return fit.best_fit, fit.residual, fit.best_values, fit.best_errors, fit.redchi, \ | |
421 fit.success | |
422 else: | |
423 return np.array([]), np.array([]), {}, {}, sys.float_info.max, False | |
424 | |
425 def fit(self, centers, fit_type=None, peak_models=None, center_exprs=None, fwhm_max=None, | |
426 background_order=None, plot_components=None, param_constraint=False): | |
427 self._fwhm_max = fwhm_max | |
428 # Create the multipeak model | |
429 self._create_model(centers, fit_type, peak_models, center_exprs, background_order, | |
430 param_constraint) | |
431 | |
432 # Perform the fit | |
433 try: | |
434 if param_constraint: | |
435 super().fit(fit_kws={'xtol' : 1.e-5, 'ftol' : 1.e-5, 'gtol' : 1.e-5}) | |
436 else: | |
437 super().fit() | |
438 except: | |
439 return False | |
440 | |
441 # Check for valid fit parameter results | |
442 fit_failure = self._check_validity() | |
443 success = True | |
444 if fit_failure: | |
445 if param_constraint: | |
446 logging.warning(' -> Should not happen with param_constraint set, fail the fit') | |
447 success = False | |
448 else: | |
449 logging.info(' -> Retry fitting with constraints') | |
450 self.fit(centers, fit_type, peak_models, center_exprs, fwhm_max=fwhm_max, | |
451 background_order=background_order, plot_components=plot_components, | |
452 param_constraint=True) | |
453 else: | |
454 # Renormalize the data and results | |
455 self._renormalize() | |
456 | |
457 # Print report and plot components if requested | |
458 if plot_components is not None: | |
459 self.print_fit_report() | |
460 self.plot() | |
461 | |
462 return success | |
463 | |
464 def _create_model(self, centers, fit_type=None, peak_models=None, center_exprs=None, | |
465 background_order=None, param_constraint=False): | |
466 """Create the multipeak model | |
467 """ | |
468 if isinstance(centers, (int, float)): | |
469 centers = [centers] | |
470 num_peaks = len(centers) | |
471 if peak_models is None: | |
472 peak_models = num_peaks*['gaussian'] | |
473 elif isinstance(peak_models, str): | |
474 peak_models = num_peaks*[peak_models] | |
475 if len(peak_models) != num_peaks: | |
476 raise ValueError(f'Inconsistent number of peaks in peak_models ({len(peak_models)} vs '+ | |
477 f'{num_peaks})') | |
478 if num_peaks == 1: | |
479 if fit_type is not None: | |
480 logging.warning('Ignoring fit_type input for fitting one peak') | |
481 fit_type = None | |
482 if center_exprs is not None: | |
483 logging.warning('Ignoring center_exprs input for fitting one peak') | |
484 center_exprs = None | |
485 else: | |
486 if fit_type == 'uniform': | |
487 if center_exprs is None: | |
488 center_exprs = [f'scale_factor*{cen}' for cen in centers] | |
489 if len(center_exprs) != num_peaks: | |
490 raise ValueError(f'Inconsistent number of peaks in center_exprs '+ | |
491 f'({len(center_exprs)} vs {num_peaks})') | |
492 elif fit_type == 'unconstrained' or fit_type is None: | |
493 if center_exprs is not None: | |
494 logging.warning('Ignoring center_exprs input for unconstrained fit') | |
495 center_exprs = None | |
496 else: | |
497 raise ValueError(f'Illegal fit_type in fit_multigaussian {fit_type}') | |
498 self._sigma_max = None | |
499 if param_constraint: | |
500 min_value = sys.float_info.min | |
501 if self._fwhm_max is not None: | |
502 self._sigma_max = np.zeros(num_peaks) | |
503 else: | |
504 min_value = None | |
505 | |
506 # Reset the fit | |
507 self._model = None | |
508 self._parameters = Parameters() | |
509 self._result = None | |
510 | |
511 # Add background model | |
512 if background_order is not None: | |
513 if background_order == 0: | |
514 self.add_model('constant', prefix='background', c=0.0) | |
515 elif background_order == 1: | |
516 self.add_model('linear', prefix='background', slope=0.0, intercept=0.0) | |
517 elif background_order == 2: | |
518 self.add_model('quadratic', prefix='background', a=0.0, b=0.0, c=0.0) | |
519 else: | |
520 raise ValueError(f'background_order = {background_order}') | |
521 | |
522 # Add peaks and guess initial fit parameters | |
523 ast = Interpreter() | |
524 if num_peaks == 1: | |
525 height_init, cen_init, fwhm_init = self.guess_init_peak(self._x, self._y) | |
526 if self._fwhm_max is not None and fwhm_init > self._fwhm_max: | |
527 fwhm_init = self._fwhm_max | |
528 ast(f'fwhm = {fwhm_init}') | |
529 ast(f'height = {height_init}') | |
530 sig_init = ast(fwhm_factor[peak_models[0]]) | |
531 amp_init = ast(height_factor[peak_models[0]]) | |
532 sig_max = None | |
533 if self._sigma_max is not None: | |
534 ast(f'fwhm = {self._fwhm_max}') | |
535 sig_max = ast(fwhm_factor[peak_models[0]]) | |
536 self._sigma_max[0] = sig_max | |
537 self.add_model(peak_models[0], parameters=( | |
538 {'name' : 'amplitude', 'value' : amp_init, 'min' : min_value}, | |
539 {'name' : 'center', 'value' : cen_init, 'min' : min_value}, | |
540 {'name' : 'sigma', 'value' : sig_init, 'min' : min_value, 'max' : sig_max})) | |
541 else: | |
542 if fit_type == 'uniform': | |
543 self.add_parameter(name='scale_factor', value=1.0) | |
544 for i in range(num_peaks): | |
545 height_init, cen_init, fwhm_init = self.guess_init_peak(self._x, self._y, i, | |
546 center_guess=centers) | |
547 if self._fwhm_max is not None and fwhm_init > self._fwhm_max: | |
548 fwhm_init = self._fwhm_max | |
549 ast(f'fwhm = {fwhm_init}') | |
550 ast(f'height = {height_init}') | |
551 sig_init = ast(fwhm_factor[peak_models[i]]) | |
552 amp_init = ast(height_factor[peak_models[i]]) | |
553 sig_max = None | |
554 if self._sigma_max is not None: | |
555 ast(f'fwhm = {self._fwhm_max}') | |
556 sig_max = ast(fwhm_factor[peak_models[i]]) | |
557 self._sigma_max[i] = sig_max | |
558 if fit_type == 'uniform': | |
559 self.add_model(peak_models[i], prefix=f'peak{i+1}_', parameters=( | |
560 {'name' : 'amplitude', 'value' : amp_init, 'min' : min_value}, | |
561 {'name' : 'center', 'expr' : center_exprs[i], 'min' : min_value}, | |
562 {'name' : 'sigma', 'value' : sig_init, 'min' : min_value, | |
563 'max' : sig_max})) | |
564 else: | |
565 self.add_model('gaussian', prefix=f'peak{i+1}_', parameters=( | |
566 {'name' : 'amplitude', 'value' : amp_init, 'min' : min_value}, | |
567 {'name' : 'center', 'value' : cen_init, 'min' : min_value}, | |
568 {'name' : 'sigma', 'value' : sig_init, 'min' : min_value, | |
569 'max' : sig_max})) | |
570 | |
571 def _check_validity(self): | |
572 """Check for valid fit parameter results | |
573 """ | |
574 fit_failure = False | |
575 index = re.compile(r'\d+') | |
576 for parameter in self.best_parameters: | |
577 name = parameter['name'] | |
578 if ((('amplitude' in name or 'height' in name) and parameter['value'] <= 0.0) or | |
579 (('sigma' in name or 'fwhm' in name) and parameter['value'] <= 0.0) or | |
580 ('center' in name and parameter['value'] <= 0.0) or | |
581 (name == 'scale_factor' and parameter['value'] <= 0.0)): | |
582 logging.info(f'Invalid fit result for {name} ({parameter["value"]})') | |
583 fit_failure = True | |
584 if 'sigma' in name and self._sigma_max is not None: | |
585 if name == 'sigma': | |
586 sigma_max = self._sigma_max[0] | |
587 else: | |
588 sigma_max = self._sigma_max[int(index.search(name).group())-1] | |
589 i = int(index.search(name).group())-1 | |
590 if parameter['value'] > sigma_max: | |
591 logging.info(f'Invalid fit result for {name} ({parameter["value"]})') | |
592 fit_failure = True | |
593 elif parameter['value'] == sigma_max: | |
594 logging.warning(f'Edge result on for {name} ({parameter["value"]})') | |
595 if 'fwhm' in name and self._fwhm_max is not None: | |
596 if parameter['value'] > self._fwhm_max: | |
597 logging.info(f'Invalid fit result for {name} ({parameter["value"]})') | |
598 fit_failure = True | |
599 elif parameter['value'] == self._fwhm_max: | |
600 logging.warning(f'Edge result on for {name} ({parameter["value"]})') | |
601 return fit_failure | |
602 | |
603 def _normalize(self): | |
604 """Normalize the data | |
605 """ | |
606 y_min = self._y.min() | |
607 self._norm = (y_min, self._y.max()-y_min) | |
608 if self._norm[1] == 0.0: | |
609 self._norm = None | |
610 else: | |
611 self._y = (self._y-self._norm[0])/self._norm[1] | |
612 | |
613 def _renormalize(self): | |
614 """Renormalize the data and results | |
615 """ | |
616 if self._norm is None: | |
617 return | |
618 self._y = self._norm[0]+self._norm[1]*self._y | |
619 self._result.best_fit = self._norm[0]+self._norm[1]*self._result.best_fit | |
620 for name in self._result.params: | |
621 par = self._result.params[name] | |
622 if 'amplitude' in name or 'height' in name or 'background' in name: | |
623 par.value *= self._norm[1] | |
624 if par.stderr is not None: | |
625 par.stderr *= self._norm[1] | |
626 if par.init_value is not None: | |
627 par.init_value *= self._norm[1] | |
628 if par.min is not None and not np.isinf(par.min): | |
629 par.min *= self._norm[1] | |
630 if par.max is not None and not np.isinf(par.max): | |
631 par.max *= self._norm[1] | |
632 if 'intercept' in name or 'backgroundc' in name: | |
633 par.value += self._norm[0] | |
634 if par.init_value is not None: | |
635 par.init_value += self._norm[0] | |
636 if par.min is not None and not np.isinf(par.min): | |
637 par.min += self._norm[0] | |
638 if par.max is not None and not np.isinf(par.max): | |
639 par.max += self._norm[0] | |
640 self._result.init_fit = self._norm[0]+self._norm[1]*self._result.init_fit | |
641 init_values = {} | |
642 for name in self._result.init_values: | |
643 init_values[name] = self._result.init_values[name] | |
644 if init_values[name] is None: | |
645 continue | |
646 if 'amplitude' in name or 'height' in name or 'background' in name: | |
647 init_values[name] *= self._norm[1] | |
648 if 'intercept' in name or 'backgroundc' in name: | |
649 init_values[name] += self._norm[0] | |
650 self._result.init_values = init_values | |
651 # Don't renormalized chisqr, it has no useful meaning in physical units | |
652 #self._result.chisqr *= self._norm[1]*self._norm[1] | |
653 if self._result.covar is not None: | |
654 for i, name in enumerate(self._result.var_names): | |
655 if 'amplitude' in name or 'height' in name or 'background' in name: | |
656 for j in range(len(self._result.var_names)): | |
657 if self._result.covar[i,j] is not None: | |
658 self._result.covar[i,j] *= self._norm[1] | |
659 if self._result.covar[j,i] is not None: | |
660 self._result.covar[j,i] *= self._norm[1] | |
661 # Don't renormalized redchi, it has no useful meaning in physical units | |
662 #self._result.redchi *= self._norm[1]*self._norm[1] | |
663 self._result.residual *= self._norm[1] |