1
|
|
|
# Licensed under a 3-clause BSD style license - see LICENSE.rst |
2
|
|
|
import itertools |
3
|
|
|
import logging |
4
|
|
|
import numpy as np |
5
|
|
|
from gammapy.utils.pbar import progress_bar |
6
|
|
|
from gammapy.utils.table import table_from_row_data |
7
|
|
|
from .covariance import Covariance |
8
|
|
|
from .iminuit import ( |
9
|
|
|
confidence_iminuit, |
10
|
|
|
contour_iminuit, |
11
|
|
|
covariance_iminuit, |
12
|
|
|
optimize_iminuit, |
13
|
|
|
) |
14
|
|
|
from .scipy import confidence_scipy, optimize_scipy |
15
|
|
|
from .sherpa import optimize_sherpa |
16
|
|
|
|
17
|
|
|
__all__ = ["Fit"] |
18
|
|
|
|
19
|
|
|
log = logging.getLogger(__name__) |
20
|
|
|
|
21
|
|
|
|
22
|
|
|
class Registry: |
23
|
|
|
"""Registry of available backends for given tasks. |
24
|
|
|
|
25
|
|
|
Gives users the power to extend from their scripts. |
26
|
|
|
Used by `Fit` below. |
27
|
|
|
|
28
|
|
|
Not sure if we should call it "backend" or "method" or something else. |
29
|
|
|
Probably we will code up some methods, e.g. for profile analysis ourselves, |
30
|
|
|
using scipy or even just Python / Numpy? |
31
|
|
|
""" |
32
|
|
|
|
33
|
|
|
register = { |
34
|
|
|
"optimize": { |
35
|
|
|
"minuit": optimize_iminuit, |
36
|
|
|
"sherpa": optimize_sherpa, |
37
|
|
|
"scipy": optimize_scipy, |
38
|
|
|
}, |
39
|
|
|
"covariance": { |
40
|
|
|
"minuit": covariance_iminuit, |
41
|
|
|
# "sherpa": covariance_sherpa, |
42
|
|
|
# "scipy": covariance_scipy, |
43
|
|
|
}, |
44
|
|
|
"confidence": { |
45
|
|
|
"minuit": confidence_iminuit, |
46
|
|
|
# "sherpa": confidence_sherpa, |
47
|
|
|
"scipy": confidence_scipy, |
48
|
|
|
}, |
49
|
|
|
} |
50
|
|
|
|
51
|
|
|
@classmethod |
52
|
|
|
def get(cls, task, backend): |
53
|
|
|
if task not in cls.register: |
54
|
|
|
raise ValueError(f"Unknown task {task!r}") |
55
|
|
|
|
56
|
|
|
backend_options = cls.register[task] |
57
|
|
|
|
58
|
|
|
if backend not in backend_options: |
59
|
|
|
raise ValueError(f"Unknown backend {backend!r} for task {task!r}") |
60
|
|
|
|
61
|
|
|
return backend_options[backend] |
62
|
|
|
|
63
|
|
|
|
64
|
|
|
registry = Registry() |
65
|
|
|
|
66
|
|
|
|
67
|
|
|
class Fit: |
68
|
|
|
"""Fit class. |
69
|
|
|
|
70
|
|
|
The fit class provides a uniform interface to multiple fitting backends. |
71
|
|
|
Currently available: "minuit", "sherpa" and "scipy" |
72
|
|
|
|
73
|
|
|
Parameters |
74
|
|
|
---------- |
75
|
|
|
backend : {"minuit", "scipy" "sherpa"} |
76
|
|
|
Global backend used for fitting, default : minuit |
77
|
|
|
optimize_opts : dict |
78
|
|
|
Keyword arguments passed to the optimizer. For the `"minuit"` backend |
79
|
|
|
see https://iminuit.readthedocs.io/en/latest/api.html#iminuit.Minuit |
80
|
|
|
for a detailed description of the available options. If there is an entry |
81
|
|
|
'migrad_opts', those options will be passed to `iminuit.Minuit.migrad()`. |
82
|
|
|
|
83
|
|
|
For the `"sherpa"` backend you can from the options `method = {"simplex", "levmar", "moncar", "gridsearch"}` |
84
|
|
|
Those methods are described and compared in detail on |
85
|
|
|
http://cxc.cfa.harvard.edu/sherpa/methods/index.html. The available |
86
|
|
|
options of the optimization methods are described on the following |
87
|
|
|
pages in detail: |
88
|
|
|
|
89
|
|
|
* http://cxc.cfa.harvard.edu/sherpa/ahelp/neldermead.html |
90
|
|
|
* http://cxc.cfa.harvard.edu/sherpa/ahelp/montecarlo.html |
91
|
|
|
* http://cxc.cfa.harvard.edu/sherpa/ahelp/gridsearch.html |
92
|
|
|
* http://cxc.cfa.harvard.edu/sherpa/ahelp/levmar.html |
93
|
|
|
|
94
|
|
|
For the `"scipy"` backend the available options are described in detail here: |
95
|
|
|
https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html |
96
|
|
|
|
97
|
|
|
covariance_opts : dict |
98
|
|
|
Covariance options passed to the given backend. |
99
|
|
|
confidence_opts : dict |
100
|
|
|
Extra arguments passed to the backend. E.g. `iminuit.Minuit.minos` supports |
101
|
|
|
a ``maxcall`` option. For the scipy backend ``confidence_opts`` are forwarded |
102
|
|
|
to `~scipy.optimize.brentq`. If the confidence estimation fails, the bracketing |
103
|
|
|
interval can be adapted by modifying the the upper bound of the interval (``b``) value. |
104
|
|
|
store_trace : bool |
105
|
|
|
Whether to store the trace of the fit |
106
|
|
|
""" |
107
|
|
|
|
108
|
|
|
def __init__( |
109
|
|
|
self, |
110
|
|
|
backend="minuit", |
111
|
|
|
optimize_opts=None, |
112
|
|
|
covariance_opts=None, |
113
|
|
|
confidence_opts=None, |
114
|
|
|
store_trace=False, |
115
|
|
|
): |
116
|
|
|
self.store_trace = store_trace |
117
|
|
|
self.backend = backend |
118
|
|
|
|
119
|
|
|
if optimize_opts is None: |
120
|
|
|
optimize_opts = {"backend": backend} |
121
|
|
|
|
122
|
|
|
if covariance_opts is None: |
123
|
|
|
covariance_opts = {"backend": backend} |
124
|
|
|
|
125
|
|
|
if confidence_opts is None: |
126
|
|
|
confidence_opts = {"backend": backend} |
127
|
|
|
|
128
|
|
|
self.optimize_opts = optimize_opts |
129
|
|
|
self.covariance_opts = covariance_opts |
130
|
|
|
self.confidence_opts = confidence_opts |
131
|
|
|
self._minuit = None |
132
|
|
|
|
133
|
|
|
@property |
134
|
|
|
def minuit(self): |
135
|
|
|
"""Iminuit object""" |
136
|
|
|
return self._minuit |
137
|
|
|
|
138
|
|
|
@staticmethod |
139
|
|
|
def _parse_datasets(datasets): |
140
|
|
|
from gammapy.datasets import Datasets |
141
|
|
|
|
142
|
|
|
datasets = Datasets(datasets) |
143
|
|
|
return datasets, datasets.parameters |
144
|
|
|
|
145
|
|
|
def run(self, datasets): |
146
|
|
|
"""Run all fitting steps. |
147
|
|
|
|
148
|
|
|
Parameters |
149
|
|
|
---------- |
150
|
|
|
datasets : `Datasets` or list of `Dataset` |
151
|
|
|
Datasets to optimize. |
152
|
|
|
|
153
|
|
|
Returns |
154
|
|
|
------- |
155
|
|
|
fit_result : `FitResult` |
156
|
|
|
Fit result |
157
|
|
|
""" |
158
|
|
|
optimize_result = self.optimize(datasets=datasets) |
159
|
|
|
|
160
|
|
|
if self.backend not in registry.register["covariance"]: |
161
|
|
|
log.warning("No covariance estimate - not supported by this backend.") |
162
|
|
|
return optimize_result |
163
|
|
|
|
164
|
|
|
covariance_result = self.covariance(datasets=datasets) |
165
|
|
|
|
166
|
|
|
return FitResult( |
167
|
|
|
optimize_result=optimize_result, |
168
|
|
|
covariance_result=covariance_result, |
169
|
|
|
) |
170
|
|
|
|
171
|
|
|
def optimize(self, datasets): |
172
|
|
|
"""Run the optimization. |
173
|
|
|
|
174
|
|
|
Parameters |
175
|
|
|
---------- |
176
|
|
|
datasets : `Datasets` or list of `Dataset` |
177
|
|
|
Datasets to optimize. |
178
|
|
|
|
179
|
|
|
Returns |
180
|
|
|
------- |
181
|
|
|
optimize_result : `OptimizeResult` |
182
|
|
|
Optimization result |
183
|
|
|
""" |
184
|
|
|
datasets, parameters = self._parse_datasets(datasets=datasets) |
185
|
|
|
datasets.parameters.check_limits() |
186
|
|
|
|
187
|
|
|
parameters.autoscale() |
188
|
|
|
|
189
|
|
|
kwargs = self.optimize_opts.copy() |
190
|
|
|
backend = kwargs.pop("backend", self.backend) |
191
|
|
|
|
192
|
|
|
compute = registry.get("optimize", backend) |
193
|
|
|
# TODO: change this calling interface! |
194
|
|
|
# probably should pass a fit statistic, which has a model, which has parameters |
195
|
|
|
# and return something simpler, not a tuple of three things |
196
|
|
|
factors, info, optimizer = compute( |
197
|
|
|
parameters=parameters, |
198
|
|
|
function=datasets.stat_sum, |
199
|
|
|
store_trace=self.store_trace, |
200
|
|
|
**kwargs, |
201
|
|
|
) |
202
|
|
|
|
203
|
|
|
if backend == "minuit": |
204
|
|
|
self._minuit = optimizer |
205
|
|
|
kwargs["method"] = "migrad" |
206
|
|
|
|
207
|
|
|
trace = table_from_row_data(info.pop("trace")) |
208
|
|
|
|
209
|
|
|
if self.store_trace: |
210
|
|
|
idx = [ |
211
|
|
|
parameters.index(par) |
212
|
|
|
for par in parameters.unique_parameters.free_parameters |
213
|
|
|
] |
214
|
|
|
unique_names = np.array(datasets.models.parameters_unique_names)[idx] |
215
|
|
|
trace.rename_columns(trace.colnames[1:], list(unique_names)) |
216
|
|
|
|
217
|
|
|
# Copy final results into the parameters object |
218
|
|
|
parameters.set_parameter_factors(factors) |
219
|
|
|
parameters.check_limits() |
220
|
|
|
return OptimizeResult( |
221
|
|
|
parameters=parameters, |
222
|
|
|
total_stat=datasets.stat_sum(), |
223
|
|
|
backend=backend, |
224
|
|
|
method=kwargs.get("method", backend), |
225
|
|
|
trace=trace, |
226
|
|
|
**info, |
227
|
|
|
) |
228
|
|
|
|
229
|
|
|
def covariance(self, datasets): |
230
|
|
|
"""Estimate the covariance matrix. |
231
|
|
|
|
232
|
|
|
Assumes that the model parameters are already optimised. |
233
|
|
|
|
234
|
|
|
Parameters |
235
|
|
|
---------- |
236
|
|
|
datasets : `Datasets` or list of `Dataset` |
237
|
|
|
Datasets to optimize. |
238
|
|
|
|
239
|
|
|
Returns |
240
|
|
|
------- |
241
|
|
|
result : `CovarianceResult` |
242
|
|
|
Results |
243
|
|
|
""" |
244
|
|
|
datasets, parameters = self._parse_datasets(datasets=datasets) |
245
|
|
|
|
246
|
|
|
kwargs = self.covariance_opts.copy() |
247
|
|
|
kwargs["minuit"] = self.minuit |
248
|
|
|
backend = kwargs.pop("backend", self.backend) |
249
|
|
|
compute = registry.get("covariance", backend) |
250
|
|
|
|
251
|
|
|
with parameters.restore_status(): |
252
|
|
|
if self.backend == "minuit": |
253
|
|
|
method = "hesse" |
254
|
|
|
else: |
255
|
|
|
method = "" |
256
|
|
|
|
257
|
|
|
factor_matrix, info = compute( |
258
|
|
|
parameters=parameters, function=datasets.stat_sum, **kwargs |
259
|
|
|
) |
260
|
|
|
|
261
|
|
|
datasets.models.covariance = Covariance.from_factor_matrix( |
262
|
|
|
parameters=parameters, matrix=factor_matrix |
263
|
|
|
) |
264
|
|
|
|
265
|
|
|
# TODO: decide what to return, and fill the info correctly! |
266
|
|
|
return CovarianceResult( |
267
|
|
|
backend=backend, |
268
|
|
|
method=method, |
269
|
|
|
success=info["success"], |
270
|
|
|
message=info["message"], |
271
|
|
|
) |
272
|
|
|
|
273
|
|
|
def confidence(self, datasets, parameter, sigma=1, reoptimize=True): |
274
|
|
|
"""Estimate confidence interval. |
275
|
|
|
|
276
|
|
|
Extra ``kwargs`` are passed to the backend. |
277
|
|
|
E.g. `iminuit.Minuit.minos` supports a ``maxcall`` option. |
278
|
|
|
|
279
|
|
|
For the scipy backend ``kwargs`` are forwarded to `~scipy.optimize.brentq`. If the |
280
|
|
|
confidence estimation fails, the bracketing interval can be adapted by modifying the |
281
|
|
|
the upper bound of the interval (``b``) value. |
282
|
|
|
|
283
|
|
|
Parameters |
284
|
|
|
---------- |
285
|
|
|
datasets : `Datasets` or list of `Dataset` |
286
|
|
|
Datasets to optimize. |
287
|
|
|
parameter : `~gammapy.modeling.Parameter` |
288
|
|
|
Parameter of interest |
289
|
|
|
sigma : float |
290
|
|
|
Number of standard deviations for the confidence level |
291
|
|
|
reoptimize : bool |
292
|
|
|
Re-optimize other parameters, when computing the confidence region. |
293
|
|
|
|
294
|
|
|
Returns |
295
|
|
|
------- |
296
|
|
|
result : dict |
297
|
|
|
Dictionary with keys "errp", 'errn", "success" and "nfev". |
298
|
|
|
""" |
299
|
|
|
datasets, parameters = self._parse_datasets(datasets=datasets) |
300
|
|
|
|
301
|
|
|
kwargs = self.confidence_opts.copy() |
302
|
|
|
backend = kwargs.pop("backend", self.backend) |
303
|
|
|
|
304
|
|
|
compute = registry.get("confidence", backend) |
305
|
|
|
parameter = parameters[parameter] |
306
|
|
|
|
307
|
|
|
with parameters.restore_status(): |
308
|
|
|
result = compute( |
309
|
|
|
parameters=parameters, |
310
|
|
|
parameter=parameter, |
311
|
|
|
function=datasets.stat_sum, |
312
|
|
|
sigma=sigma, |
313
|
|
|
reoptimize=reoptimize, |
314
|
|
|
**kwargs, |
315
|
|
|
) |
316
|
|
|
|
317
|
|
|
result["errp"] *= parameter.scale |
318
|
|
|
result["errn"] *= parameter.scale |
319
|
|
|
return result |
320
|
|
|
|
321
|
|
|
def stat_profile(self, datasets, parameter, reoptimize=False): |
322
|
|
|
"""Compute fit statistic profile. |
323
|
|
|
|
324
|
|
|
The method used is to vary one parameter, keeping all others fixed. |
325
|
|
|
So this is taking a "slice" or "scan" of the fit statistic. |
326
|
|
|
|
327
|
|
|
Parameters |
328
|
|
|
---------- |
329
|
|
|
datasets : `Datasets` or list of `Dataset` |
330
|
|
|
Datasets to optimize. |
331
|
|
|
parameter : `~gammapy.modeling.Parameter` |
332
|
|
|
Parameter of interest. The specification for the scan, such as bounds |
333
|
|
|
and number of values is taken from the parameter object. |
334
|
|
|
reoptimize : bool |
335
|
|
|
Re-optimize other parameters, when computing the confidence region. |
336
|
|
|
|
337
|
|
|
Returns |
338
|
|
|
------- |
339
|
|
|
results : dict |
340
|
|
|
Dictionary with keys "values", "stat" and "fit_results". The latter contains an |
341
|
|
|
empty list, if `reoptimize` is set to False |
342
|
|
|
""" |
343
|
|
|
datasets, parameters = self._parse_datasets(datasets=datasets) |
344
|
|
|
parameter = parameters[parameter] |
345
|
|
|
values = parameter.scan_values |
346
|
|
|
|
347
|
|
|
stats = [] |
348
|
|
|
fit_results = [] |
349
|
|
|
with parameters.restore_status(): |
350
|
|
|
for value in progress_bar(values, desc="Scan values"): |
351
|
|
|
parameter.value = value |
352
|
|
|
if reoptimize: |
353
|
|
|
parameter.frozen = True |
354
|
|
|
result = self.optimize(datasets=datasets) |
355
|
|
|
stat = result.total_stat |
356
|
|
|
fit_results.append(result) |
357
|
|
|
else: |
358
|
|
|
stat = datasets.stat_sum() |
359
|
|
|
stats.append(stat) |
360
|
|
|
|
361
|
|
|
return { |
362
|
|
|
f"{parameter.name}_scan": values, |
363
|
|
|
"stat_scan": np.array(stats), |
364
|
|
|
"fit_results": fit_results, |
365
|
|
|
} |
366
|
|
|
|
367
|
|
|
def stat_surface(self, datasets, x, y, reoptimize=False): |
368
|
|
|
"""Compute fit statistic surface. |
369
|
|
|
|
370
|
|
|
The method used is to vary two parameters, keeping all others fixed. |
371
|
|
|
So this is taking a "slice" or "scan" of the fit statistic. |
372
|
|
|
|
373
|
|
|
Caveat: This method can be very computationally intensive and slow |
374
|
|
|
|
375
|
|
|
See also: `Fit.stat_contour` |
376
|
|
|
|
377
|
|
|
Parameters |
378
|
|
|
---------- |
379
|
|
|
datasets : `Datasets` or list of `Dataset` |
380
|
|
|
Datasets to optimize. |
381
|
|
|
x, y : `~gammapy.modeling.Parameter` |
382
|
|
|
Parameters of interest |
383
|
|
|
reoptimize : bool |
384
|
|
|
Re-optimize other parameters, when computing the confidence region. |
385
|
|
|
|
386
|
|
|
Returns |
387
|
|
|
------- |
388
|
|
|
results : dict |
389
|
|
|
Dictionary with keys "x_values", "y_values", "stat" and "fit_results". The latter contains an |
390
|
|
|
empty list, if `reoptimize` is set to False |
391
|
|
|
""" |
392
|
|
|
datasets, parameters = self._parse_datasets(datasets=datasets) |
393
|
|
|
|
394
|
|
|
x, y = parameters[x], parameters[y] |
395
|
|
|
|
396
|
|
|
stats = [] |
397
|
|
|
fit_results = [] |
398
|
|
|
|
399
|
|
|
with parameters.restore_status(): |
400
|
|
|
for x_value, y_value in progress_bar( |
401
|
|
|
itertools.product(x.scan_values, y.scan_values), desc="Trial values" |
402
|
|
|
): |
403
|
|
|
x.value, y.value = x_value, y_value |
404
|
|
|
|
405
|
|
|
if reoptimize: |
406
|
|
|
x.frozen, y.frozen = True, True |
407
|
|
|
result = self.optimize(datasets=datasets) |
408
|
|
|
stat = result.total_stat |
409
|
|
|
fit_results.append(result) |
410
|
|
|
else: |
411
|
|
|
stat = datasets.stat_sum() |
412
|
|
|
|
413
|
|
|
stats.append(stat) |
414
|
|
|
|
415
|
|
|
shape = (len(x.scan_values), len(y.scan_values)) |
416
|
|
|
stats = np.array(stats).reshape(shape) |
417
|
|
|
|
418
|
|
|
if reoptimize: |
419
|
|
|
fit_results = np.array(fit_results).reshape(shape) |
420
|
|
|
|
421
|
|
|
return { |
422
|
|
|
f"{x.name}_scan": x.scan_values, |
423
|
|
|
f"{y.name}_scan": y.scan_values, |
424
|
|
|
"stat_scan": stats, |
425
|
|
|
"fit_results": fit_results, |
426
|
|
|
} |
427
|
|
|
|
428
|
|
|
def stat_contour(self, datasets, x, y, numpoints=10, sigma=1): |
429
|
|
|
"""Compute stat contour. |
430
|
|
|
|
431
|
|
|
Calls ``iminuit.Minuit.mncontour``. |
432
|
|
|
|
433
|
|
|
This is a contouring algorithm for a 2D function |
434
|
|
|
which is not simply the fit statistic function. |
435
|
|
|
That 2D function is given at each point ``(par_1, par_2)`` |
436
|
|
|
by re-optimising all other free parameters, |
437
|
|
|
and taking the fit statistic at that point. |
438
|
|
|
|
439
|
|
|
Very compute-intensive and slow. |
440
|
|
|
|
441
|
|
|
Parameters |
442
|
|
|
---------- |
443
|
|
|
datasets : `Datasets` or list of `Dataset` |
444
|
|
|
Datasets to optimize. |
445
|
|
|
x, y : `~gammapy.modeling.Parameter` |
446
|
|
|
Parameters of interest |
447
|
|
|
numpoints : int |
448
|
|
|
Number of contour points |
449
|
|
|
sigma : float |
450
|
|
|
Number of standard deviations for the confidence level |
451
|
|
|
|
452
|
|
|
Returns |
453
|
|
|
------- |
454
|
|
|
result : dict |
455
|
|
|
Dictionary containing the parameter values defining the contour, with the |
456
|
|
|
boolean flag "success" and the info objects from ``mncontour``. |
457
|
|
|
""" |
458
|
|
|
datasets, parameters = self._parse_datasets(datasets=datasets) |
459
|
|
|
|
460
|
|
|
x = parameters[x] |
461
|
|
|
y = parameters[y] |
462
|
|
|
|
463
|
|
|
with parameters.restore_status(): |
464
|
|
|
result = contour_iminuit( |
465
|
|
|
parameters=parameters, |
466
|
|
|
function=datasets.stat_sum, |
467
|
|
|
x=x, |
468
|
|
|
y=y, |
469
|
|
|
numpoints=numpoints, |
470
|
|
|
sigma=sigma, |
471
|
|
|
) |
472
|
|
|
|
473
|
|
|
x_name = x.name |
474
|
|
|
y_name = y.name |
475
|
|
|
x = result["x"] * x.scale |
476
|
|
|
y = result["y"] * y.scale |
477
|
|
|
|
478
|
|
|
return { |
479
|
|
|
x_name: x, |
480
|
|
|
y_name: y, |
481
|
|
|
"success": result["success"], |
482
|
|
|
} |
483
|
|
|
|
484
|
|
|
|
485
|
|
|
class FitStepResult: |
486
|
|
|
"""Fit result base class""" |
487
|
|
|
|
488
|
|
|
def __init__(self, backend, method, success, message): |
489
|
|
|
self._success = success |
490
|
|
|
self._message = message |
491
|
|
|
self._backend = backend |
492
|
|
|
self._method = method |
493
|
|
|
|
494
|
|
|
@property |
495
|
|
|
def backend(self): |
496
|
|
|
"""Optimizer backend used for the fit.""" |
497
|
|
|
return self._backend |
498
|
|
|
|
499
|
|
|
@property |
500
|
|
|
def method(self): |
501
|
|
|
"""Optimizer method used for the fit.""" |
502
|
|
|
return self._method |
503
|
|
|
|
504
|
|
|
@property |
505
|
|
|
def success(self): |
506
|
|
|
"""Fit success status flag.""" |
507
|
|
|
return self._success |
508
|
|
|
|
509
|
|
|
@property |
510
|
|
|
def message(self): |
511
|
|
|
"""Optimizer status message.""" |
512
|
|
|
return self._message |
513
|
|
|
|
514
|
|
|
def __repr__(self): |
515
|
|
|
return ( |
516
|
|
|
f"{self.__class__.__name__}\n\n" |
517
|
|
|
f"\tbackend : {self.backend}\n" |
518
|
|
|
f"\tmethod : {self.method}\n" |
519
|
|
|
f"\tsuccess : {self.success}\n" |
520
|
|
|
f"\tmessage : {self.message}\n" |
521
|
|
|
) |
522
|
|
|
|
523
|
|
|
|
524
|
|
|
class CovarianceResult(FitStepResult): |
525
|
|
|
"""Covariance result object.""" |
526
|
|
|
|
527
|
|
|
pass |
528
|
|
|
|
529
|
|
|
|
530
|
|
|
class OptimizeResult(FitStepResult): |
531
|
|
|
"""Optimize result object.""" |
532
|
|
|
|
533
|
|
|
def __init__(self, parameters, nfev, total_stat, trace, **kwargs): |
534
|
|
|
self._parameters = parameters |
535
|
|
|
self._nfev = nfev |
536
|
|
|
self._total_stat = total_stat |
537
|
|
|
self._trace = trace |
538
|
|
|
super().__init__(**kwargs) |
539
|
|
|
|
540
|
|
|
@property |
541
|
|
|
def parameters(self): |
542
|
|
|
"""Best fit parameters""" |
543
|
|
|
return self._parameters |
544
|
|
|
|
545
|
|
|
@property |
546
|
|
|
def trace(self): |
547
|
|
|
"""Parameter trace from the optimisation""" |
548
|
|
|
return self._trace |
549
|
|
|
|
550
|
|
|
@property |
551
|
|
|
def nfev(self): |
552
|
|
|
"""Number of function evaluations.""" |
553
|
|
|
return self._nfev |
554
|
|
|
|
555
|
|
|
@property |
556
|
|
|
def total_stat(self): |
557
|
|
|
"""Value of the fit statistic at minimum.""" |
558
|
|
|
return self._total_stat |
559
|
|
|
|
560
|
|
|
def __repr__(self): |
561
|
|
|
str_ = super().__repr__() |
562
|
|
|
str_ += f"\tnfev : {self.nfev}\n" |
563
|
|
|
str_ += f"\ttotal stat : {self.total_stat:.2f}\n\n" |
564
|
|
|
return str_ |
565
|
|
|
|
566
|
|
|
|
567
|
|
|
class FitResult: |
568
|
|
|
"""Fit result class |
569
|
|
|
|
570
|
|
|
Parameters |
571
|
|
|
---------- |
572
|
|
|
optimize_result : `OptimizeResult` |
573
|
|
|
Result of the optimization step. |
574
|
|
|
covariance_result : `CovarianceResult` |
575
|
|
|
Result of the covariance step. |
576
|
|
|
""" |
577
|
|
|
|
578
|
|
|
def __init__(self, optimize_result=None, covariance_result=None): |
579
|
|
|
self._optimize_result = optimize_result |
580
|
|
|
self._covariance_result = covariance_result |
581
|
|
|
|
582
|
|
|
# TODO: is the convenience access needed? |
583
|
|
|
@property |
584
|
|
|
def parameters(self): |
585
|
|
|
"""Best fit parameters of the optimization step""" |
586
|
|
|
return self.optimize_result.parameters |
587
|
|
|
|
588
|
|
|
# TODO: is the convenience access needed? |
589
|
|
|
@property |
590
|
|
|
def total_stat(self): |
591
|
|
|
"""Total stat of the optimization step""" |
592
|
|
|
return self.optimize_result.total_stat |
593
|
|
|
|
594
|
|
|
# TODO: is the convenience access needed? |
595
|
|
|
@property |
596
|
|
|
def trace(self): |
597
|
|
|
"""Parameter trace of the optimisation step""" |
598
|
|
|
return self.optimize_result.trace |
599
|
|
|
|
600
|
|
|
# TODO: is the convenience access needed? |
601
|
|
|
@property |
602
|
|
|
def nfev(self): |
603
|
|
|
"""Number of function evaluations of the optimisation step""" |
604
|
|
|
return self.optimize_result.nfev |
605
|
|
|
|
606
|
|
|
# TODO: is the convenience access needed? |
607
|
|
|
@property |
608
|
|
|
def backend(self): |
609
|
|
|
"""Optimizer backend used for the fit.""" |
610
|
|
|
return self.optimize_result.backend |
611
|
|
|
|
612
|
|
|
# TODO: is the convenience access needed? |
613
|
|
|
@property |
614
|
|
|
def method(self): |
615
|
|
|
"""Optimizer method used for the fit.""" |
616
|
|
|
return self.optimize_result.method |
617
|
|
|
|
618
|
|
|
# TODO: is the convenience access needed? |
619
|
|
|
@property |
620
|
|
|
def message(self): |
621
|
|
|
"""Optimizer status message.""" |
622
|
|
|
return self.optimize_result.message |
623
|
|
|
|
624
|
|
|
@property |
625
|
|
|
def success(self): |
626
|
|
|
"""Total success flag""" |
627
|
|
|
success = self.optimize_result.success and self.covariance_result.success |
628
|
|
|
return success |
629
|
|
|
|
630
|
|
|
@property |
631
|
|
|
def optimize_result(self): |
632
|
|
|
"""Optimize result""" |
633
|
|
|
return self._optimize_result |
634
|
|
|
|
635
|
|
|
@property |
636
|
|
|
def covariance_result(self): |
637
|
|
|
"""Optimize result""" |
638
|
|
|
return self._optimize_result |
639
|
|
|
|
640
|
|
|
def __repr__(self): |
641
|
|
|
str_ = "" |
642
|
|
|
if self.optimize_result: |
643
|
|
|
str_ += str(self.optimize_result) |
644
|
|
|
|
645
|
|
|
if self.covariance_result: |
646
|
|
|
str_ += str(self.covariance_result) |
647
|
|
|
|
648
|
|
|
return str_ |
649
|
|
|
|