|
1
|
|
|
# -*- coding: utf-8 -*- |
|
2
|
|
|
|
|
3
|
|
|
""" |
|
4
|
|
|
This file contains methods for gaussian-like fitting, these methods |
|
5
|
|
|
are imported by class FitLogic. |
|
6
|
|
|
|
|
7
|
|
|
Qudi is free software: you can redistribute it and/or modify |
|
8
|
|
|
it under the terms of the GNU General Public License as published by |
|
9
|
|
|
the Free Software Foundation, either version 3 of the License, or |
|
10
|
|
|
(at your option) any later version. |
|
11
|
|
|
|
|
12
|
|
|
Qudi is distributed in the hope that it will be useful, |
|
13
|
|
|
but WITHOUT ANY WARRANTY; without even the implied warranty of |
|
14
|
|
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
|
15
|
|
|
GNU General Public License for more details. |
|
16
|
|
|
|
|
17
|
|
|
You should have received a copy of the GNU General Public License |
|
18
|
|
|
along with Qudi. If not, see <http://www.gnu.org/licenses/>. |
|
19
|
|
|
|
|
20
|
|
|
Copyright (c) the Qudi Developers. See the COPYRIGHT.txt file at the |
|
21
|
|
|
top-level directory of this distribution and at <https://github.com/Ulm-IQO/qudi/> |
|
22
|
|
|
""" |
|
23
|
|
|
|
|
24
|
|
|
|
|
25
|
|
|
import logging |
|
26
|
|
|
logger = logging.getLogger(__name__) |
|
27
|
|
|
import numpy as np |
|
28
|
|
|
from lmfit.models import Model, GaussianModel, ConstantModel |
|
29
|
|
|
from lmfit import Parameters |
|
30
|
|
|
|
|
31
|
|
|
from scipy.interpolate import InterpolatedUnivariateSpline |
|
32
|
|
|
from scipy.signal import gaussian |
|
33
|
|
|
from scipy.ndimage import filters |
|
34
|
|
|
|
|
35
|
|
|
############################################################################ |
|
36
|
|
|
# # |
|
37
|
|
|
# 1D gaussian model # |
|
38
|
|
|
# # |
|
39
|
|
|
############################################################################ |
|
40
|
|
|
|
|
41
|
|
|
def make_gaussian_model(self): |
|
42
|
|
|
""" This method creates a model of a gaussian with an offset. |
|
43
|
|
|
|
|
44
|
|
|
@return tuple: (object model, object params) |
|
45
|
|
|
|
|
46
|
|
|
Explanation of the objects: |
|
47
|
|
|
object lmfit.model.CompositeModel model: |
|
48
|
|
|
A model the lmfit module will use for that fit. Here a |
|
49
|
|
|
gaussian model. Returns an object of the class |
|
50
|
|
|
lmfit.model.CompositeModel. |
|
51
|
|
|
|
|
52
|
|
|
object lmfit.parameter.Parameters params: |
|
53
|
|
|
It is basically an OrderedDict, so a dictionary, with keys |
|
54
|
|
|
denoting the parameters as string names and values which are |
|
55
|
|
|
lmfit.parameter.Parameter (without s) objects, keeping the |
|
56
|
|
|
information about the current value. |
|
57
|
|
|
The used model has the Parameter with the meaning: |
|
58
|
|
|
'amplitude' : amplitude |
|
59
|
|
|
'center' : center |
|
60
|
|
|
'sigm' : sigma |
|
61
|
|
|
'fwhm' : full width half maximum |
|
62
|
|
|
'c' : offset |
|
63
|
|
|
|
|
64
|
|
|
For further information have a look in: |
|
65
|
|
|
http://cars9.uchicago.edu/software/python/lmfit/builtin_models.html#models.GaussianModel |
|
66
|
|
|
""" |
|
67
|
|
|
|
|
68
|
|
|
model = GaussianModel() + ConstantModel() |
|
69
|
|
|
params = model.make_params() |
|
70
|
|
|
|
|
71
|
|
|
return model, params |
|
72
|
|
|
|
|
73
|
|
|
def make_gaussian_fit(self, axis=None, data=None, add_parameters=None): |
|
74
|
|
|
""" This method performes a 1D gaussian fit on the provided data. |
|
75
|
|
|
|
|
76
|
|
|
@param array[] axis: axis values |
|
77
|
|
|
@param array[] x_data: data |
|
78
|
|
|
@param dict add_parameters: Additional parameters which will substitute the |
|
79
|
|
|
estimated parameters/bounds |
|
80
|
|
|
|
|
81
|
|
|
@return object result: lmfit.model.ModelFit object, all parameters |
|
82
|
|
|
provided about the fitting, like: success, |
|
83
|
|
|
initial fitting values, best fitting values, data |
|
84
|
|
|
with best fit with given axis,... |
|
85
|
|
|
""" |
|
86
|
|
|
|
|
87
|
|
|
mod_final, params = self.make_gaussian_model() |
|
88
|
|
|
|
|
89
|
|
|
error, params = self.estimate_gaussian(axis, data, params) |
|
90
|
|
|
|
|
91
|
|
|
# auxiliary variables |
|
92
|
|
|
stepsize = abs(axis[1] - axis[0]) |
|
93
|
|
|
n_steps = len(axis) |
|
94
|
|
|
|
|
95
|
|
|
# Define constraints |
|
96
|
|
|
params['center'].min = (axis[0]) - n_steps * stepsize |
|
97
|
|
|
params['center'].max = (axis[-1]) + n_steps * stepsize |
|
98
|
|
|
params['amplitude'].min = 100 # that is already noise from APD |
|
99
|
|
|
params['amplitude'].max = data.max() * params['sigma'].value * np.sqrt(2 * np.pi) |
|
100
|
|
|
params['sigma'].min = stepsize |
|
101
|
|
|
params['sigma'].max = 3 * (axis[-1] - axis[0]) |
|
102
|
|
|
params['c'].min = 100 # that is already noise from APD |
|
103
|
|
|
params['c'].max = data.max() * params['sigma'].value * np.sqrt(2 * np.pi) |
|
104
|
|
|
|
|
105
|
|
|
# overwrite values of additional parameters |
|
106
|
|
|
if add_parameters is not None: |
|
107
|
|
|
params = self._substitute_parameter(parameters=params, |
|
108
|
|
|
update_dict=add_parameters) |
|
109
|
|
|
try: |
|
110
|
|
|
result = mod_final.fit(data, x=axis, params=params) |
|
111
|
|
|
except: |
|
112
|
|
|
logger.warning('The 1D gaussian fit did not work.') |
|
113
|
|
|
result = mod_final.fit(data, x=axis, params=params) |
|
114
|
|
|
print(result.message) |
|
115
|
|
|
|
|
116
|
|
|
return result |
|
117
|
|
|
|
|
118
|
|
|
def estimate_gaussian(self, x_axis=None, data=None, params=None): |
|
119
|
|
|
""" This method provides a one dimensional gaussian function. |
|
120
|
|
|
|
|
121
|
|
|
@param array x_axis: x values |
|
122
|
|
|
@param array data: value of each data point corresponding to x values |
|
123
|
|
|
@param Parameters object params: object includes parameter dictionary which can be set |
|
124
|
|
|
|
|
125
|
|
|
@return tuple (error, params): |
|
126
|
|
|
|
|
127
|
|
|
Explanation of the return parameter: |
|
128
|
|
|
int error: error code (0:OK, -1:error) |
|
129
|
|
|
Parameters object params: set parameters of initial values |
|
130
|
|
|
""" |
|
131
|
|
|
|
|
132
|
|
|
error = 0 |
|
133
|
|
|
# check if parameters make sense |
|
134
|
|
|
parameters = [x_axis, data] |
|
135
|
|
|
for var in parameters: |
|
136
|
|
|
if not isinstance(var, (frozenset, list, set, tuple, np.ndarray)): |
|
137
|
|
|
logger.error('Given parameter is no array.') |
|
138
|
|
|
error = -1 |
|
139
|
|
|
elif len(np.shape(var)) != 1: |
|
140
|
|
|
logger.error('Given parameter is no one dimensional array.') |
|
141
|
|
|
error = -1 |
|
142
|
|
|
if not isinstance(params, Parameters): |
|
143
|
|
|
logger.error('Parameters object is not valid in estimate_gaussian.') |
|
144
|
|
|
error = -1 |
|
145
|
|
|
|
|
146
|
|
|
# If the estimator is not good enough one can start improvement with |
|
147
|
|
|
# a convolution |
|
148
|
|
|
|
|
149
|
|
|
# set parameters |
|
150
|
|
|
params['center'].value = x_axis[np.argmax(data)] |
|
151
|
|
|
params['sigma'].value = (x_axis.max() - x_axis.min()) / 3. |
|
152
|
|
|
params['amplitude'].value = (data.max() - data.min()) * (params['sigma'].value * np.sqrt(2 * np.pi)) |
|
153
|
|
|
params['c'].value = data.min() |
|
154
|
|
|
|
|
155
|
|
|
return error, params |
|
156
|
|
|
|
|
157
|
|
|
############################################################################ |
|
158
|
|
|
# # |
|
159
|
|
|
# 2D gaussian model # |
|
160
|
|
|
# # |
|
161
|
|
|
############################################################################ |
|
162
|
|
|
1 |
|
163
|
|
|
|
|
164
|
|
|
def make_twoDgaussian_fit(self, axis=None, data=None, |
|
165
|
|
|
add_parameters=None, estimator="estimate_twoDgaussian_MLE"): |
|
166
|
|
|
""" This method performes a 2D gaussian fit on the provided data. |
|
167
|
|
|
|
|
168
|
|
|
@param array[] axis: axis values |
|
169
|
|
|
@param array[] x_data: data |
|
170
|
|
|
@param dict add_parameters: Additional parameters |
|
171
|
|
|
@param string estimator: the string should contain the name of the function you want to use to estimate |
|
172
|
|
|
the parameters. The default estimator is estimate_twoDgaussian_MLE. |
|
173
|
|
|
|
|
174
|
|
|
@return object result: lmfit.model.ModelFit object, all parameters |
|
175
|
|
|
provided about the fitting, like: success, |
|
176
|
|
|
initial fitting values, best fitting values, data |
|
177
|
|
|
with best fit with given axis,... |
|
178
|
|
|
""" |
|
179
|
|
|
|
|
180
|
|
|
x_axis, y_axis = axis |
|
181
|
|
|
|
|
182
|
|
|
if estimator is "estimate_twoDgaussian_MLE": |
|
183
|
|
|
error, \ |
|
184
|
|
|
amplitude, \ |
|
185
|
|
|
x_zero, \ |
|
186
|
|
|
y_zero, \ |
|
187
|
|
|
sigma_x, \ |
|
188
|
|
|
sigma_y, \ |
|
189
|
|
|
theta, \ |
|
190
|
|
|
offset = self.estimate_twoDgaussian_MLE(x_axis=x_axis, |
|
191
|
|
|
y_axis=y_axis, data=data) |
|
192
|
|
|
else: |
|
193
|
|
|
error, \ |
|
194
|
|
|
amplitude, \ |
|
195
|
|
|
x_zero, \ |
|
196
|
|
|
y_zero, \ |
|
197
|
|
|
sigma_x, \ |
|
198
|
|
|
sigma_y, \ |
|
199
|
|
|
theta, \ |
|
200
|
|
|
offset = globals()[estimator](0, x_axis=x_axis, |
|
201
|
|
|
y_axis=y_axis, data=data) |
|
202
|
|
|
|
|
203
|
|
|
mod, params = self.make_twoDgaussian_model() |
|
204
|
|
|
|
|
205
|
|
|
#auxiliary variables |
|
206
|
|
|
stepsize_x=x_axis[1]-x_axis[0] |
|
207
|
|
|
stepsize_y=y_axis[1]-y_axis[0] |
|
208
|
|
|
n_steps_x=len(x_axis) |
|
209
|
|
|
n_steps_y=len(y_axis) |
|
210
|
|
|
|
|
211
|
|
|
#When I was sitting in the train coding and my girlfiend was sitting next to me she said: "Look it looks like an animal!" - is it a fox or a rabbit??? |
|
212
|
|
|
|
|
213
|
|
|
#Defining standard parameters |
|
214
|
|
|
# (Name, Value, Vary, Min, Max, Expr) |
|
215
|
|
|
params.add_many(('amplitude', amplitude, True, 100, 1e7, None), |
|
216
|
|
|
( 'sigma_x', sigma_x, True, 1*(stepsize_x) , 3*(x_axis[-1]-x_axis[0]), None), |
|
217
|
|
|
( 'sigma_y', sigma_y, True, 1*(stepsize_y) , 3*(y_axis[-1]-y_axis[0]) , None), |
|
218
|
|
|
( 'x_zero', x_zero, True, (x_axis[0])-n_steps_x*stepsize_x , x_axis[-1]+n_steps_x*stepsize_x, None), |
|
219
|
|
|
( 'y_zero', y_zero, True, (y_axis[0])-n_steps_y*stepsize_y , (y_axis[-1])+n_steps_y*stepsize_y, None), |
|
220
|
|
|
( 'theta', 0., True, 0. , np.pi, None), |
|
221
|
|
|
( 'offset', offset, True, 0, 1e7, None)) |
|
222
|
|
|
|
|
223
|
|
|
|
|
224
|
|
|
# redefine values of additional parameters |
|
225
|
|
|
if add_parameters is not None: |
|
226
|
|
|
params=self._substitute_parameter(parameters=params, |
|
227
|
|
|
update_dict=add_parameters) |
|
228
|
|
|
|
|
229
|
|
|
try: |
|
230
|
|
|
result=mod.fit(data, x=axis,params=params) |
|
231
|
|
|
except: |
|
232
|
|
|
result=mod.fit(data, x=axis,params=params) |
|
233
|
|
|
logger.warning('The 2D gaussian fit did not work: {0}'.format( |
|
234
|
|
|
result.message)) |
|
235
|
|
|
|
|
236
|
|
|
return result |
|
237
|
|
|
|
|
238
|
|
|
|
|
239
|
|
|
def make_twoDgaussian_model(self): |
|
240
|
|
|
""" This method creates a model of the 2D gaussian function. |
|
241
|
|
|
|
|
242
|
|
|
The parameters are: 'amplitude', 'center', 'sigm, 'fwhm' and offset |
|
243
|
|
|
'c'. For function see: |
|
244
|
|
|
|
|
245
|
|
|
@return lmfit.model.CompositeModel model: Returns an object of the |
|
246
|
|
|
class CompositeModel |
|
247
|
|
|
@return lmfit.parameter.Parameters params: Returns an object of the |
|
248
|
|
|
class Parameters with all |
|
249
|
|
|
parameters for the |
|
250
|
|
|
gaussian model. |
|
251
|
|
|
|
|
252
|
|
|
""" |
|
253
|
|
|
|
|
254
|
|
|
def twoDgaussian_function(x, amplitude, x_zero, y_zero, sigma_x, sigma_y, |
|
255
|
|
|
theta, offset): |
|
256
|
|
|
# FIXME: x_data_tuple: dimension of arrays |
|
257
|
|
|
|
|
258
|
|
|
""" This method provides a two dimensional gaussian function. |
|
259
|
|
|
|
|
260
|
|
|
Function taken from: |
|
261
|
|
|
http://stackoverflow.com/questions/21566379/fitting-a-2d-gaussian-function-using-scipy-optimize-curve-fit-valueerror-and-m/21566831 |
|
262
|
|
|
|
|
263
|
|
|
Question from: http://stackoverflow.com/users/2097737/bland & http://stackoverflow.com/users/3273102/kokomoking |
|
264
|
|
|
& http://stackoverflow.com/users/2767207/jojodmo |
|
265
|
|
|
Answer: http://stackoverflow.com/users/1461210/ali-m & http://stackoverflow.com/users/5234/mrjrdnthms |
|
266
|
|
|
|
|
267
|
|
|
@param array[k][M] x_data_tuple: array which is (k,M)-shaped, x and y |
|
268
|
|
|
values |
|
269
|
|
|
@param float or int amplitude: Amplitude of gaussian |
|
270
|
|
|
@param float or int x_zero: x value of maximum |
|
271
|
|
|
@param float or int y_zero: y value of maximum |
|
272
|
|
|
@param float or int sigma_x: standard deviation in x direction |
|
273
|
|
|
@param float or int sigma_y: standard deviation in y direction |
|
274
|
|
|
@param float or int theta: angle for eliptical gaussians |
|
275
|
|
|
@param float or int offset: offset |
|
276
|
|
|
|
|
277
|
|
|
@return callable function: returns the function |
|
278
|
|
|
""" |
|
279
|
|
|
|
|
280
|
|
|
(u, v) = x |
|
281
|
|
|
x_zero = float(x_zero) |
|
282
|
|
|
y_zero = float(y_zero) |
|
283
|
|
|
|
|
284
|
|
|
a = (np.cos(theta) ** 2) / (2 * sigma_x ** 2) \ |
|
285
|
|
|
+ (np.sin(theta) ** 2) / (2 * sigma_y ** 2) |
|
286
|
|
|
b = -(np.sin(2 * theta)) / (4 * sigma_x ** 2) \ |
|
287
|
|
|
+ (np.sin(2 * theta)) / (4 * sigma_y ** 2) |
|
288
|
|
|
c = (np.sin(theta) ** 2) / (2 * sigma_x ** 2) \ |
|
289
|
|
|
+ (np.cos(theta) ** 2) / (2 * sigma_y ** 2) |
|
290
|
|
|
g = offset + amplitude * np.exp(- (a * ((u - x_zero) ** 2) \ |
|
291
|
|
|
+ 2 * b * (u - x_zero) * (v - y_zero) \ |
|
292
|
|
|
+ c * ((v - y_zero) ** 2))) |
|
293
|
|
|
return g.ravel() |
|
294
|
|
|
|
|
295
|
|
|
model = Model(twoDgaussian_function) |
|
296
|
|
|
params = model.make_params() |
|
297
|
|
|
|
|
298
|
|
|
return model, params |
|
299
|
|
|
|
|
300
|
|
|
|
|
301
|
|
|
def estimate_twoDgaussian(self, x_axis=None, y_axis=None, data=None): |
|
302
|
|
|
# TODO:Make clever estimator |
|
303
|
|
|
# FIXME: 1D array x_axis, y_axis, 2D data??? |
|
304
|
|
|
""" This method provides a two dimensional gaussian function. |
|
305
|
|
|
|
|
306
|
|
|
@param array x_axis: x values |
|
307
|
|
|
@param array y_axis: y values |
|
308
|
|
|
@param array data: value of each data point corresponding to |
|
309
|
|
|
x and y values |
|
310
|
|
|
|
|
311
|
|
|
@return float amplitude: estimated amplitude |
|
312
|
|
|
@return float x_zero: estimated x value of maximum |
|
313
|
|
|
@return float y_zero: estimated y value of maximum |
|
314
|
|
|
@return float sigma_x: estimated standard deviation in x direction |
|
315
|
|
|
@return float sigma_y: estimated standard deviation in y direction |
|
316
|
|
|
@return float theta: estimated angle for eliptical gaussians |
|
317
|
|
|
@return float offset: estimated offset |
|
318
|
|
|
@return int error: error code (0:OK, -1:error) |
|
319
|
|
|
""" |
|
320
|
|
|
|
|
321
|
|
|
# #needed me 1 hour to think about, but not needed in the end...maybe needed at a later point |
|
322
|
|
|
# len_x=np.where(x_axis[0]==x_axis)[0][1] |
|
323
|
|
|
# len_y=len(data)/len_x |
|
324
|
|
|
|
|
325
|
|
|
|
|
326
|
|
|
amplitude = float(data.max() - data.min()) |
|
327
|
|
|
|
|
328
|
|
|
x_zero = x_axis[data.argmax()] |
|
329
|
|
|
y_zero = y_axis[data.argmax()] |
|
330
|
|
|
|
|
331
|
|
|
sigma_x = (x_axis.max() - x_axis.min()) / 3. |
|
332
|
|
|
sigma_y = (y_axis.max() - y_axis.min()) / 3. |
|
333
|
|
|
theta = 0.0 |
|
334
|
|
|
offset = float(data.min()) |
|
335
|
|
|
error = 0 |
|
336
|
|
|
# check for sensible values |
|
337
|
|
|
parameters = [x_axis, y_axis, data] |
|
338
|
|
|
for var in parameters: |
|
339
|
|
|
# FIXME: Why don't you check earlier? |
|
340
|
|
|
# FIXME: Check for 1D array, 2D |
|
341
|
|
|
if not isinstance(var, (frozenset, list, set, tuple, np.ndarray)): |
|
342
|
|
|
logger.error('Given parameter is not an array.') |
|
343
|
|
|
amplitude = 0. |
|
344
|
|
|
x_zero = 0. |
|
345
|
|
|
y_zero = 0. |
|
346
|
|
|
sigma_x = 0. |
|
347
|
|
|
sigma_y = 0. |
|
348
|
|
|
theta = 0.0 |
|
349
|
|
|
offset = 0. |
|
350
|
|
|
error = -1 |
|
351
|
|
|
|
|
352
|
|
|
return error, amplitude, x_zero, y_zero, sigma_x, sigma_y, theta, offset |
|
353
|
|
|
|
|
354
|
|
View Code Duplication |
def estimate_twoDgaussian_MLE(self, x_axis=None, y_axis=None, data=None): |
|
|
|
|
|
|
355
|
|
|
# TODO: Make good estimates for sigma_x, sigma_y and theta |
|
356
|
|
|
""" This method provides an estimate for the parameters characterizing a |
|
357
|
|
|
two dimensional gaussian. It is based on the maximum likelihood estimation |
|
358
|
|
|
(at the moment only for the x_zero and y_zero values). |
|
359
|
|
|
|
|
360
|
|
|
@param array x_axis: x values |
|
361
|
|
|
@param array y_axis: y values |
|
362
|
|
|
@param array data: value of each data point corresponding to |
|
363
|
|
|
x and y values |
|
364
|
|
|
@return tuple parameters: estimated value of parameters based on the data |
|
365
|
|
|
""" |
|
366
|
|
|
|
|
367
|
|
|
amplitude = float(data.max() - data.min()) |
|
368
|
|
|
|
|
369
|
|
|
x_zero = np.sum(x_axis * data) / np.sum(data) |
|
370
|
|
|
y_zero = np.sum(y_axis * data) / np.sum(data) |
|
371
|
|
|
|
|
372
|
|
|
sigma_x = (x_axis.max() - x_axis.min()) / 3. |
|
373
|
|
|
sigma_y = (y_axis.max() - y_axis.min()) / 3. |
|
374
|
|
|
theta = 0.0 |
|
375
|
|
|
offset = float(data.min()) |
|
376
|
|
|
error = 0 |
|
377
|
|
|
# check for sensible values |
|
378
|
|
|
parameters = [x_axis, y_axis, data] |
|
379
|
|
|
for var in parameters: |
|
380
|
|
|
# FIXME: Why don't you check earlier? |
|
381
|
|
|
# FIXME: Check for 1D array, 2D |
|
382
|
|
|
if not isinstance(var, (frozenset, list, set, tuple, np.ndarray)): |
|
383
|
|
|
logger.error('Given parameter is not an array.') |
|
384
|
|
|
amplitude = 0. |
|
385
|
|
|
x_zero = 0. |
|
386
|
|
|
y_zero = 0. |
|
387
|
|
|
sigma_x = 0. |
|
388
|
|
|
sigma_y = 0. |
|
389
|
|
|
theta = 0.0 |
|
390
|
|
|
offset = 0. |
|
391
|
|
|
error = -1 |
|
392
|
|
|
|
|
393
|
|
|
return error, amplitude, x_zero, y_zero, sigma_x, sigma_y, theta, offset |
|
394
|
|
|
|
|
395
|
|
View Code Duplication |
def estimate_twoDgaussian_MLE(self, x_axis=None, y_axis=None, data=None): |
|
|
|
|
|
|
396
|
|
|
# TODO: Make good estimates for sigma_x, sigma_y and theta |
|
397
|
|
|
""" This method provides an estimate for the parameters characterizing a |
|
398
|
|
|
two dimensional gaussian. It is based on the maximum likelihood estimation |
|
399
|
|
|
(at the moment only for the x_zero and y_zero values). |
|
400
|
|
|
|
|
401
|
|
|
@param array x_axis: x values |
|
402
|
|
|
@param array y_axis: y values |
|
403
|
|
|
@param array data: value of each data point corresponding to |
|
404
|
|
|
x and y values |
|
405
|
|
|
@return tuple parameters: estimated value of parameters based on the data |
|
406
|
|
|
""" |
|
407
|
|
|
|
|
408
|
|
|
amplitude = float(data.max() - data.min()) |
|
409
|
|
|
|
|
410
|
|
|
x_zero = np.sum(x_axis * data) / np.sum(data) |
|
411
|
|
|
y_zero = np.sum(y_axis * data) / np.sum(data) |
|
412
|
|
|
|
|
413
|
|
|
sigma_x = (x_axis.max() - x_axis.min()) / 3. |
|
414
|
|
|
sigma_y = (y_axis.max() - y_axis.min()) / 3. |
|
415
|
|
|
theta = 0.0 |
|
416
|
|
|
offset = float(data.min()) |
|
417
|
|
|
error = 0 |
|
418
|
|
|
# check for sensible values |
|
419
|
|
|
parameters = [x_axis, y_axis, data] |
|
420
|
|
|
for var in parameters: |
|
421
|
|
|
# FIXME: Why don't you check earlier? |
|
422
|
|
|
# FIXME: Check for 1D array, 2D |
|
423
|
|
|
if not isinstance(var, (frozenset, list, set, tuple, np.ndarray)): |
|
424
|
|
|
logger.error('Given parameter is not an array.') |
|
425
|
|
|
amplitude = 0. |
|
426
|
|
|
x_zero = 0. |
|
427
|
|
|
y_zero = 0. |
|
428
|
|
|
sigma_x = 0. |
|
429
|
|
|
sigma_y = 0. |
|
430
|
|
|
theta = 0.0 |
|
431
|
|
|
offset = 0. |
|
432
|
|
|
error = -1 |
|
433
|
|
|
|
|
434
|
|
|
return error, amplitude, x_zero, y_zero, sigma_x, sigma_y, theta, offset |
|
435
|
|
|
|
|
436
|
|
|
|
|
437
|
|
|
############################################################################ |
|
438
|
|
|
# # |
|
439
|
|
|
# Double Gaussian Model # |
|
440
|
|
|
# # |
|
441
|
|
|
############################################################################ |
|
442
|
|
|
|
|
443
|
|
|
def make_multiplegaussian_model(self, no_of_gauss=None): |
|
444
|
|
|
""" This method creates a model of multiple gaussians with an offset. The |
|
445
|
|
|
parameters are: 'amplitude', 'center', 'sigma', 'fwhm' and offset |
|
446
|
|
|
'c'. For function see: |
|
447
|
|
|
http://cars9.uchicago.edu/software/python/lmfit/builtin_models.html#models.LorentzianModel |
|
448
|
|
|
|
|
449
|
|
|
@return lmfit.model.CompositeModel model: Returns an object of the |
|
450
|
|
|
class CompositeModel |
|
451
|
|
|
@return lmfit.parameter.Parameters params: Returns an object of the |
|
452
|
|
|
class Parameters with all |
|
453
|
|
|
parameters for the |
|
454
|
|
|
lorentzian model. |
|
455
|
|
|
""" |
|
456
|
|
|
|
|
457
|
|
|
model = ConstantModel() |
|
458
|
|
|
for ii in range(no_of_gauss): |
|
459
|
|
|
model += GaussianModel(prefix='gaussian{0}_'.format(ii)) |
|
460
|
|
|
|
|
461
|
|
|
params = model.make_params() |
|
462
|
|
|
|
|
463
|
|
|
return model, params |
|
464
|
|
|
|
|
465
|
|
|
|
|
466
|
|
|
def estimate_doublegaussian_gatedcounter(self, x_axis=None, data=None, params=None, |
|
467
|
|
|
threshold_fraction=0.4, |
|
468
|
|
|
minimal_threshold=0.1, |
|
469
|
|
|
sigma_threshold_fraction=0.2): |
|
470
|
|
|
""" This method provides a an estimator for a double gaussian fit with the parameters |
|
471
|
|
|
coming from the physical properties of an experiment done in gated counter: |
|
472
|
|
|
- positive peak |
|
473
|
|
|
- no values below 0 |
|
474
|
|
|
- rather broad overlapping funcitons |
|
475
|
|
|
|
|
476
|
|
|
@param array x_axis: x values |
|
477
|
|
|
@param array data: value of each data point corresponding to |
|
478
|
|
|
x values |
|
479
|
|
|
@param Parameters object params: Needed parameters |
|
480
|
|
|
@param float threshold : Threshold to find second gaussian |
|
481
|
|
|
@param float minimal_threshold: Threshold is lowerd to minimal this |
|
482
|
|
|
value as a fraction |
|
483
|
|
|
@param float sigma_threshold_fraction: Threshold for detecting |
|
484
|
|
|
the end of the peak |
|
485
|
|
|
|
|
486
|
|
|
@return int error: error code (0:OK, -1:error) |
|
487
|
|
|
@return Parameters object params: estimated values |
|
488
|
|
|
""" |
|
489
|
|
|
|
|
490
|
|
|
error = 0 |
|
491
|
|
|
|
|
492
|
|
|
data_smooth = self.gaussian_smoothing(data=data) |
|
493
|
|
|
|
|
494
|
|
|
# search for double gaussian |
|
495
|
|
|
|
|
496
|
|
|
error, \ |
|
497
|
|
|
sigma0_argleft, dip0_arg, sigma0_argright, \ |
|
498
|
|
|
sigma1_argleft, dip1_arg, sigma1_argright = \ |
|
499
|
|
|
self._search_double_dip(x_axis, |
|
500
|
|
|
data_smooth * (-1), |
|
501
|
|
|
threshold_fraction, |
|
502
|
|
|
minimal_threshold, |
|
503
|
|
|
sigma_threshold_fraction, |
|
504
|
|
|
make_prints=False |
|
505
|
|
|
) |
|
506
|
|
|
|
|
507
|
|
|
# set offset to zero |
|
508
|
|
|
params['c'].value = 0.0 |
|
509
|
|
|
|
|
510
|
|
|
params['gaussian0_center'].value = x_axis[dip0_arg] |
|
511
|
|
|
|
|
512
|
|
|
# integral of data corresponds to sqrt(2) * Amplitude * Sigma |
|
513
|
|
|
function = InterpolatedUnivariateSpline(x_axis, data_smooth, k=1) |
|
514
|
|
|
Integral = function.integral(x_axis[0], x_axis[-1]) |
|
515
|
|
|
|
|
516
|
|
|
amp_0 = data_smooth[dip0_arg] - params['c'].value |
|
517
|
|
|
amp_1 = data_smooth[dip1_arg] - params['c'].value |
|
518
|
|
|
|
|
519
|
|
|
params['gaussian0_sigma'].value = Integral / (amp_0 + amp_1) / np.sqrt(2 * np.pi) |
|
520
|
|
|
params['gaussian0_amplitude'].value = amp_0 * params['gaussian0_sigma'].value * np.sqrt(2 * np.pi) |
|
521
|
|
|
|
|
522
|
|
|
params['gaussian1_center'].value = x_axis[dip1_arg] |
|
523
|
|
|
params['gaussian1_sigma'].value = Integral / (amp_0 + amp_1) / np.sqrt(2 * np.pi) |
|
524
|
|
|
params['gaussian1_amplitude'].value = amp_1 * params['gaussian1_sigma'].value * np.sqrt(2 * np.pi) |
|
525
|
|
|
|
|
526
|
|
|
return error, params |
|
527
|
|
|
|
|
528
|
|
|
|
|
529
|
|
|
def estimate_doublegaussian_odmr(self, x_axis=None, data=None, params=None, |
|
530
|
|
|
threshold_fraction=0.4, |
|
531
|
|
|
minimal_threshold=0.1, |
|
532
|
|
|
sigma_threshold_fraction=0.2): |
|
533
|
|
|
""" This method provides a an estimator for a double gaussian fit with the parameters |
|
534
|
|
|
coming from the physical properties of an experiment done in gated counter: |
|
535
|
|
|
- positive peak |
|
536
|
|
|
- no values below 0 |
|
537
|
|
|
- rather broad overlapping funcitons |
|
538
|
|
|
|
|
539
|
|
|
@param array x_axis: x values |
|
540
|
|
|
@param array data: value of each data point corresponding to |
|
541
|
|
|
x values |
|
542
|
|
|
@param Parameters object params: Needed parameters |
|
543
|
|
|
@param float threshold : Threshold to find second gaussian |
|
544
|
|
|
@param float minimal_threshold: Threshold is lowerd to minimal this |
|
545
|
|
|
value as a fraction |
|
546
|
|
|
@param float sigma_threshold_fraction: Threshold for detecting |
|
547
|
|
|
the end of the peak |
|
548
|
|
|
|
|
549
|
|
|
@return int error: error code (0:OK, -1:error) |
|
550
|
|
|
@return Parameters object params: estimated values |
|
551
|
|
|
""" |
|
552
|
|
|
|
|
553
|
|
|
error, \ |
|
554
|
|
|
params['gaussian0_amplitude'].value, \ |
|
555
|
|
|
params['gaussian1_amplitude'].value, \ |
|
556
|
|
|
params['gaussian0_center'].value, \ |
|
557
|
|
|
params['gaussian1_center'].value, \ |
|
558
|
|
|
params['gaussian0_sigma'].value, \ |
|
559
|
|
|
params['gaussian1_sigma'].value, \ |
|
560
|
|
|
params['c'].value = self.estimate_doublelorentz(x_axis, data) |
|
561
|
|
|
|
|
562
|
|
|
return error, params |
|
563
|
|
|
|
|
564
|
|
|
|
|
565
|
|
|
def make_doublegaussian_fit(self, axis=None, data=None, |
|
566
|
|
|
add_parameters=None, |
|
567
|
|
|
estimator='gated_counter', |
|
568
|
|
|
threshold_fraction=0.4, |
|
569
|
|
|
minimal_threshold=0.2, |
|
570
|
|
|
sigma_threshold_fraction=0.3): |
|
571
|
|
|
""" This method performes a 1D double gaussian fit on the provided data. |
|
572
|
|
|
|
|
573
|
|
|
@param array [] axis: axis values |
|
574
|
|
|
@param array[] data: data |
|
575
|
|
|
@param dictionary add_parameters: Additional parameters |
|
576
|
|
|
@param float threshold_fraction : Threshold to find second gaussian |
|
577
|
|
|
@param float minimal_threshold: Threshold is lowerd to minimal this |
|
578
|
|
|
value as a fraction |
|
579
|
|
|
@param float sigma_threshold_fraction: Threshold for detecting |
|
580
|
|
|
the end of the peak |
|
581
|
|
|
|
|
582
|
|
|
@return lmfit.model.ModelFit result: All parameters provided about |
|
583
|
|
|
the fitting, like: success, |
|
584
|
|
|
initial fitting values, best |
|
585
|
|
|
fitting values, data with best |
|
586
|
|
|
fit with given axis,... |
|
587
|
|
|
|
|
588
|
|
|
""" |
|
589
|
|
|
|
|
590
|
|
|
model, params = self.make_multiplegaussian_model(no_of_gauss=2) |
|
591
|
|
|
|
|
592
|
|
|
if estimator == 'gated_counter': |
|
593
|
|
|
error, params = self.estimate_doublegaussian_gatedcounter(axis, data, params, |
|
594
|
|
|
threshold_fraction, |
|
595
|
|
|
minimal_threshold, |
|
596
|
|
|
sigma_threshold_fraction) |
|
597
|
|
|
# Defining constraints |
|
598
|
|
|
params['c'].min = 0.0 |
|
599
|
|
|
|
|
600
|
|
|
params['gaussian0_amplitude'].min = 0.0 |
|
601
|
|
|
params['gaussian1_amplitude'].min = 0.0 |
|
602
|
|
|
|
|
603
|
|
|
elif estimator == 'odmr_dip': |
|
604
|
|
|
error, params = self.estimate_doublegaussian_odmr(axis, data, params, |
|
605
|
|
|
threshold_fraction, |
|
606
|
|
|
minimal_threshold, |
|
607
|
|
|
sigma_threshold_fraction) |
|
608
|
|
|
|
|
609
|
|
|
# redefine values of additional parameters |
|
610
|
|
|
if add_parameters is not None: |
|
611
|
|
|
params = self._substitute_parameter(parameters=params, |
|
612
|
|
|
update_dict=add_parameters) |
|
613
|
|
|
try: |
|
614
|
|
|
result = model.fit(data, x=axis, params=params) |
|
615
|
|
|
except: |
|
616
|
|
|
result = model.fit(data, x=axis, params=params) |
|
617
|
|
|
logger.warning('The double gaussian fit did not work: {0}'.format( |
|
618
|
|
|
result.message)) |
|
619
|
|
|
|
|
620
|
|
|
return result |
|
621
|
|
|
|
Ulm-IQO /
qudi
Loading history...