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] |
