Commit d2461113 authored by Christian Schneider's avatar Christian Schneider
Browse files
parents 17b458fe c2cd0583
......@@ -29,10 +29,17 @@ def lorentz_function(x, x0, gamma, offset, amplitude):
\mathrm{offset}
"""
return amplitude*(gamma/2)**2/((x-x0)**2 + (gamma/2)**2) + offset
return amplitude * (gamma / 2) ** 2 / (
(x - x0) ** 2 + (gamma / 2) ** 2) + offset
def exp_func(x, tau, amp, vert_offset, hor_offset):
return amp*np.exp(-(x+hor_offset)/tau) + vert_offset
return amp * np.exp(-(x + hor_offset) / tau) + vert_offset
def T2_func(x, tau, f, phase, amp, vert_offset, hor_offset):
return (amp * np.cos(2 * np.pi * f * (x + hor_offset) + phase) *
np.exp(-(x + hor_offset) / tau) + vert_offset)
class fit_plugin(object):
......@@ -45,31 +52,32 @@ class fit_plugin(object):
data : datamodule
Datamodule with stored data
"""
self.data = data
self.data = data
def fit_routine(self, gmodel, params, plot=True, print_results=True,
plot_init=False):
"""Fit and plot/print results"""
try:
result = gmodel.fit(self.data.y, params, x=self.data.x)
except:
raise('Fit did not converge. Try to tune guess parameters')
# Plot result ############################################################
raise ('Fit did not converge. Try to tune guess parameters')
# Plot result ##########################################################
if plot:
p = bp.figure(plot_width=800, plot_height=400)
# add both a line and circles on the same plot
p.circle(self.data.x, self.data.y, size=4)
p.line(self.data.x, result.best_fit, color='firebrick', line_width=2)
p.line(self.data.x, result.best_fit, color='firebrick',
line_width=2)
if plot_init:
p.line(self.data.x, result.init_fit, color='green', line_width=2)
p.line(self.data.x, result.init_fit, color='green',
line_width=2)
bp.show(p)
# Parameters to give back ################################################
# Parameters to give back ##############################################
r_sq = 1 - result.residual.var() / np.var(self.data.y)
# Print results
......@@ -84,25 +92,42 @@ class fit_plugin(object):
self.r_squared = r_sq
# Save fitresult to datamodule
f_dict = {'Value': [result.params[k].value for k in result.params.keys()],
'Error': [result.params[k].stderr for k in result.params.keys()]}
f_dict = {
'Value': [result.params[k].value for k in result.params.keys()],
'Error': [result.params[k].stderr for k in result.params.keys()]}
self.data.fitresults = pd.DataFrame(data=f_dict,
index=list(result.params.keys()),
columns=['Value', 'Error'])
# Qubit fit functions ######################################################
def T1(self, tau0=None, amp0=None, vert_offset0=None, hor_offset0=None,
plot=True, print_results=True, plot_init=False):
"""Fit exponential decay to data"""
"""Fit exponential decay to data
Parameters
-----------
tau0 : int, float
Decay time in units of x axis
amp0 : float
Amplitude in units of y axis
vert_offset0 : float
Y Offset in y units
hor_offset : float
X Offset in x units
plot : bool
Plot fit after successful fitting
print_results : bool
Print results after successful fitting
plot_init : bool
Plot initial guess for parameters
"""
# Guess values #########################################################
# Vertical offset
if vert_offset0 is None:
vert_offset_guess = np.mean(self.data.y)
else:
vert_offset_guess = vert_offset0
# Horizontal offset
if hor_offset0 is None:
hor_offset_guess = self.data.x[0]
......@@ -114,25 +139,107 @@ class fit_plugin(object):
amp_guess = self.data.y[0] - self.data.y[-1]
else:
amp_guess = amp0
# Decay time
if tau0 is None:
tau_guess = (np.abs(self.data.x[-1] - self.data.x[0]))/10
tau_guess = (np.abs(self.data.x[-1] - self.data.x[0])) / 10
else:
tau_guess = tau0
# Create LM Fit model ##################################################
gmodel = Model(exp_func)
# Init pars to guess values ############################################
params = gmodel.make_params(tau=tau_guess,
amp=amp_guess,
vert_offset=vert_offset_guess,
hor_offset=hor_offset_guess)
# Fitting ##############################################################
self.fit_routine(gmodel, params, plot, print_results, plot_init)
def T2(self, tau0=None, f0=None, phase0=None, amp0=None, vert_offset0=None,
hor_offset0=None, plot=True, print_results=True, plot_init=False):
"""Fit T2 decay to data
Function used for fitting
math:`amp \cos{2\pi f (x-hor_offset) + \phase} \exp{-(
x-hor_offset)/\tau) + vert_offset`
"""
# Guess values #########################################################
# Vertical offset
# Guess: mean value
if vert_offset0 is None:
vert_offset_guess = np.mean(self.data.y)
else:
vert_offset_guess = vert_offset0
# Horizontal offset
# Guess: first element in x
if hor_offset0 is None:
hor_offset_guess = self.data.x[0]
else:
hor_offset_guess = hor_offset0
# Decay time
# Guess: A tenth of x range
if tau0 is None:
tau_guess = (np.abs(self.data.x[-1] - self.data.x[0])) / 10
else:
tau_guess = tau0
# -- Do fft for guessing frequency, amplitude and phase ----------------
freqs = np.fft.fftfreq(self.data.x.size,
self.data.x[1] - self.data.x[0])
fft = np.fft.fft(self.data.y)
# Remove DC
freqs = freqs[1:]
fft = fft[1:]
# Find frequencies (biggest peaks in spectrum)
srtd = np.argsort(np.abs(fft))
# ----------------------------------------------------------------------
# Amplitude
# Guess from fft
if amp0 is None:
amp_guess = np.abs(fft)[srtd[-1]] / fft.size + np.abs(fft)[
srtd[-2]] / fft.size
amp_guess *= 5 # Guessing that cos decaying over a fifth of xrange
else:
amp_guess = amp0
# Phase
# Guess from fft
if phase0 is None:
phase_guess = np.abs(
np.angle(fft[srtd[-1]]) - np.angle(fft[srtd[-2]]))
else:
phase_guess = phase0
# Frequency
# Guess from fft
if f0 is None:
f_guess = freqs[srtd[-1]]
else:
f_guess = f0
# Create LM Fit model ##################################################
gmodel = Model(T2_func)
# Init pars to guess values ############################################
params = gmodel.make_params(tau=tau_guess,
f=f_guess,
phase=phase_guess,
amp=amp_guess,
vert_offset=vert_offset_guess,
hor_offset=hor_offset_guess)
# Fitting ##############################################################
self.fit_routine(gmodel, params, plot, print_results, plot_init)
def lorentzian(self, x0=None, gamma=None, amplitude=None, offset=None,
plot=True, print_results=True, plot_init=False):
......@@ -149,45 +256,47 @@ class fit_plugin(object):
"""
# Guess values ###########################################################
# x0
if x0 is None:
# Guess using find_peaks
peaks, _ = find_peaks(np.abs(self.data.y), distance=len(self.data.y))
peaks, _ = find_peaks(np.abs(self.data.y),
distance=len(self.data.y))
x0_guess = self.data.x[peaks[0]]
else:
x0_guess = x0
# Ampltiude
if amplitude is None:
# Guess using find_peaks
peaks, _ = find_peaks(np.abs(self.data.y), distance=len(self.data.y))
peaks, _ = find_peaks(np.abs(self.data.y),
distance=len(self.data.y))
amp_guess = self.data.y[peaks[0]]
else:
amp_guess = amplitude
# Width
if gamma is None:
# Guess a fraction of range
gamma_guess = np.abs(self.data.x[-1] - self.data.x[0])/100
gamma_guess = np.abs(self.data.x[-1] - self.data.x[0]) / 100
else:
gamma_guess = gamma
if offset is None:
# Guess mean value
offset_guess = np.mean(self.data.y)
else:
offset_guess = offset
# Create LM Fit model ####################################################
gmodel = Model(lorentz_function)
# Init pars to guess values ##############################################
params = gmodel.make_params(x0=x0_guess,
amplitude=amp_guess,
gamma=gamma_guess,
offset=offset_guess)
# Fitting ################################################################
self.fit_routine(gmodel, params, plot, print_results, plot_init)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment