#pylint: disable = C0111
from scipy.optimize import curve_fit
import numpy as np
import warnings
import progressbar
import matplotlib.pyplot as plt
from scipy.stats import linregress
import pickle
import rot_en
from scipy.constants import e,k
from scipy.ndimage.filters import median_filter, percentile_filter, gaussian_filter
#Aby to slo picklovat, nesmi se definovat az ve zdedene funkci
def default_func_Alpha(pwr, alpha, beta):
return alpha * pwr / (1 + beta*pwr)
#Aby to slo picklovat, nesmi se definovat az ve zdedene funkci
def default_func_FluoVsPwr(E, alpha, B, t):
if alpha < 0 or B < 0 or t < 0:
ret = np.empty(E.shape)
ret[:] = np.inf
return ret
ret = alpha*(E - B*E**2/(B*E + t) * ((np.exp(-(B*E+t))-1)/(B*E+t)+1))
return ret
def func_FluoVsPwr_dalpha(E, alpha, B, t):
return -B*E**2*(1 + (np.exp(-B*E - t) - 1)/(B*E + t))/(B*E + t) + E
def func_FluoVsPwr_dB(E, alpha, B, t):
return alpha*(B*E**3*(1 + (np.exp(-B*E - t) - 1)/(B*E + t))/(B*E + t)**2 - B*E**2*(-E*np.exp(-B*E - t)/(B*E + t) - E*(np.exp(-B*E - t) - 1)/(B*E + t)**2)/(B*E + t) - E**2*(1 + (np.exp(-B*E - t) - 1)/(B*E + t))/(B*E + t))
def func_FluoVsPwr_dt(E, alpha, B, t):
return B*E**2*alpha*(B*E + t + (B*E + t)*np.exp(B*E + t) - 2*np.exp(B*E + t) + 2)*np.exp(-B*E - t)/(B*E + t)**3
def linfunc(x, slope, intercept):
return slope*x + intercept
def default_func_tau(time, amplitude, tau, offset):
if tau<0 or amplitude<0:
ret = np.empty(time.shape)
ret[:] = np.inf
return ret
ret = amplitude * np.exp(-time/tau) + offset
return ret
def func_damplitude(time, amplitude, tau):
return np.exp(-time/tau)
def func_dtau(time, amplitude, tau):
return amplitude * time / tau**2 * func_damplitude(time, amplitude, tau)
def func_doffset(time, amplitude, tau):
return np.ones(time.shape)
[docs]def read_lifbase_table(filename):
"""
Read database file exported from LIFBASE.
Args:
**filename**: (*str*) path to the file created by LIFBASE to be read
Returns:
dictionary of coefficients with keys ``N P1 P2 Q1 Q2 R1 R2 O12 Q12 P12 R21 Q21 S21``
Example of use::
coefs = read_lifbase_table('abs_00.txt')
print coefs['R1'][coefs['N']==4]
"""
keys = 'N P1 P2 Q1 Q2 R1 R2 O12 Q12 P12 R21 Q21 S21'.split()
arr = np.genfromtxt(filename, delimiter = ',', skip_header = 4)
coefs = {}
for i,head in enumerate(keys):
coefs[head] = arr[:,i]
return coefs
[docs]class Crossfit(object):
"""Parent class for fitting planar LIF data. Expects series of 2D images
(triDData) and some x-axis (x) that enables fitting of function (func)
to the measured data for each pixel.
Attributes:
**data**: (3D numpy array) a series of 2D images (triDData)
independent_variable: (numpy array) the independent variable for fitting.
May be 1D with the length equal to the number of frames of data, 2D with
the first dimension equal to the first dimension of self.data and the second
dimension equal to the second or third dimension of self.data or 3D with the
same shape as self.data
"""
def __init__(self, triDData, independet_variable, func=None, info={}):
self.data = triDData
self.independent_variable = independet_variable
self.func = func
self.result = {}
self.params = []
self.info=info
[docs] def default_func(self):
"""
To be provided in inherited classes...
"""
raise NotImplemented()
[docs] def func_from_keys(self):
"""
To be provided in inherited classes...
"""
raise NotImplemented()
[docs] def dfunc_from_keys(self):
"""
To be provided in inherited classes...
"""
raise NotImplemented()
[docs] def get_params(self, y, x):
"""
To be provided in inherited classes...
"""
raise NotImplemented()
[docs] def fit_at_yx(self, y, x, keys='all'):
"""Find the best estimate of parameters of self.default_func for data at (y, x) with scipy.optimize.curve_fit().
Args:
**y, x**: (*int*) pixel coordinates
**keys**: (*list*) of keywords describing the argument to be
optimized or simply 'all'. The keywords should be provided by
inherited classes in their definition.
"""
p_vals = []
for key in keys:
p_vals.append(self.result[key][y, x])
func = lambda time, *params: self.func_from_keys(keys, params, y, x)
dfunc = lambda params, time, yvals, f: self.dfunc_from_keys(keys, params, y, x)
#pylint: disable=W0632
indep = self.format_x(y, x)
try:
popt, pcov = curve_fit(func, indep, self.data[:,y,x], p0=p_vals,maxfev=2000, Dfun=dfunc, col_deriv=1)
except RuntimeError:
n=len(keys)
popt = np.empty(n)
popt[:] = np.nan
pcov = np.empty((n,n))
pcov[:,:] = np.nan
#pylint: enable=W0632
for i,key in enumerate(keys):
self.result[key][y, x] = popt[i]
#if key == 'tau':
#self.result[key][y, x] *= maxtime
try:
self.result['d'+key][y,x] = pcov[i,i]**0.5
except TypeError:
self.result['d'+key][y,x] = np.inf
self.result['cov', (y,x)] = pcov
self.result['cov_order'] = keys
[docs] def fit(self, keys='all'):
"""
call self.fit_at_yx() for each pixel in a loop.
"""
if keys == 'all':
keys = self.params
#keys = ['amplitude', 'tau', 'offset']
zmax, ymax, xmax = self.data.shape
pbar = progressbar.ProgressBar(ymax*xmax)
for y in xrange(ymax):
for x in xrange(xmax):
self.fit_at_yx(y, x, keys)
pbar.animate(y*xmax+x + 1)
[docs] def show_fit(self, y, x, save_to_file=None, xlabel='', ylabel='', draw=False, show=True):
"""Method for plotting the the cross-section of triDData at position
y, x vs. the independent_variable for that particular position,
together with fitted funcion.
Args:
**y, x**: pixel coordinates
**save_to_file**: If a string is given, it represents the name of
the file the figure will be saved to. If None (dafault), no file will be written.
**xlabel, ylabel**: (*strings*) passed directly to matplotlib
**draw**: defaults to False. if True, matplotlib.draw() will be called to show the
figure. May be used in loop to create an animation, as ax.clear() is also called.
**show**: defaults to True - the figure will be shown by matplotlib.show()
Returns: matplotlib.figure object.
"""
xfit = self.independent_var_to_plot(y, x)
p = self.get_params(y, x)
yfit = self.func(xfit, *p)
fig, ax = plt.subplots(1, 1)
xmeas = self.format_x(y, x)
ax.plot(xmeas, self.data[:, y, x], '+')
ax.plot(xfit, yfit)
ax.set_xlabel(xlabel)
ax.set_ylabel(ylabel)
ax.set_title(str(p))
if save_to_file:
fig.savefig(save_to_file)
elif draw:
plt.draw()
plt.pause(0.0001)
plt.gca().clear()
elif show==True:
plt.show()
return fig
[docs] def independent_var_to_plot(self, y, x, numpoints=100):
"""Returns linspace(xmin, xmax, numpoints), where xmin, xmax are the
minimum and maximum of the independent_variable at position y,
x. Intended to yield x-axis values for theoretical data form
the fits.
Args:
**y, x**: pixel coordinates
**numpoints**: (*integer*) defaults to 100. Number of equally spaced points to be returned
Returns:
1D numpy array from numpy.linspace
"""
xmeas = self.format_x(y, x)
m1 = np.min(xmeas)
m2 = np.max(xmeas)
xfit = np.linspace(m1, m2, numpoints)
return xfit
[docs] def get_covariance_matrix(self, y, x):
"""calculates and returns the real variance-covariance matrix from the
scipu.optimize.curve_fit's output by multiplying it with
(len(data)-len(parameters)) / sum_of_squares
"""
indep = self.format_x(y, x)
p = self.get_params(y, x)
dep = self.func(indep, *p)
sumsq = np.sum((self.data[:, y, x] - dep)**2)
try:
cov = np.zeros(self.result['cov', (y, x)].shape)
cov = self.result['cov', (y, x)]*(len(indep) - cov.shape[0]) / sumsq
except (TypeError,AttributeError) as e:
msg = 'get_covariance_matrix(): ' + str(e)
warnings.warn(msg, Warning)
cov = None
return cov
[docs] def get_corr_matrix(self, y, x):
"""return matrix of correlation coefficients. See result['cov_order']
for assignign to the right variables.
"""
mat = self.get_covariance_matrix(y, x)
corr = np.zeros((mat.shape))
for param1 in xrange(mat.shape[0]):
for param2 in xrange(mat.shape[1]):
corr[param1, param2] = mat[param1, param2]/np.sqrt(mat[param1, param1] * mat[param2, param2])
return corr
def save(self, filename):
with open(filename, 'wb') as output:
pickle.dump(self, output, pickle.HIGHEST_PROTOCOL)
@staticmethod
def load(filename):
with open(filename, 'rb') as inp:
return_val = pickle.load(inp)
return return_val
[docs]class Alpha(Crossfit):
"""Class for working with planar LIF measurements of fluorescence
vs. laser power.
Attributes:
**func**: default is :math:`\\alpha E_L/(1 + \\beta E_L)`. The fit()
method assumes only this form and works with linearized version for
inverted values of laser power. Changing this attribute is currently
pointless, but might become useful in future versions.
**M**: (*float*) the smaller this number, the shorter range of values is
permitted for refit(). For decent data M=2 can lead to cca 20%
overestimation of alpha and M = 1 to cca 10%
overestimation. Underestimation is unlikely.
"""
def __init__(self, triDData, laser_pwr, params=None, func=None, info={}):
Crossfit.__init__(self, triDData, laser_pwr ,
func=func, info=info)
self.raw = {}
self.M=None
self.betamedx=True
if self.func is None:
self.func = default_func_Alpha
def get_params(self, y, x):
return [self.result['alpha'][y, x], self.result['beta'][y, x]]
# @staticmethod
# def default_dfunc(params, pwr, y, f):
# alpha = params[0]
# beta = params[1]
# return [pwr/(1+beta*pwr), -alpha *pwr**2/ (1+beta*pwr)**2]
[docs] def invert(self, data):
"""Calculate inverse of given array and mark values that are invalid
(i.e. == np.nan, == np.inf or < 0).
Args:
**data**: numpy array
Returns:
(*tuple*) inverted, mask
inverted: inverted VALID data. The invalid values are omitted
mask: a boolean array of the shape of the original data that
can be used to filter arrays of the same size according to the invalid values in the original data.
"""
inverted = 1./data
mask = np.ones(data.shape) == 1
mask[np.isinf(inverted)] = False
mask[inverted<0] = False
mask[np.isnan(inverted)] = False
return inverted[mask], mask
[docs] def fit(self):
"""Fitting the measured fluorescence vs. laser power data
pixelwise. Assumes the dependence fluorescence = alpha * pwr / (1 +
beta*pwr). No additive offset is allowed. The results of the fit are
stored in the dictionary self.result with keys 'alpha', 'beta',
'delta_alpha'. The values of delta_beta (estimate of standard
deviation of beta) can be obtained by the method delta_beta(y, x)
(single value) or delta_beta_map() (a 2D array for all pixels)
Returns:
**alpha**: 2D numpy array of the alpha parameters
"""
zmax, ymax, xmax = self.data.shape
pbar = progressbar.ProgressBar(ymax*xmax)
self.result['alpha'] = np.zeros((ymax, xmax))
self.result['beta'] = np.zeros((ymax, xmax))
self.result['delta_alpha'] = np.zeros((ymax, xmax))
for y in xrange(ymax):
for x in xrange(xmax):
pwr = self.format_x(y, x)
maxpwr = max(pwr)
pwr /= maxpwr ##rescale for better fit
data_to_fit, mask = self.invert(self.data[:,y,x])
pwr_to_fit = pwr[mask]**(-1)
slope, intercept, r, p, stderr = linregress(pwr_to_fit, data_to_fit)
pwr *= maxpwr ## undo rescaling
self.result['alpha'][y, x] = 1./(slope*maxpwr)
self.result['beta'][y, x] = intercept * self.result['alpha'][y, x]
self.result['delta_alpha'][y,x] = stderr/(slope**2*maxpwr)
self.raw[(y, x)] = {'slope':slope, 'intercept': intercept, 'r':r, 'p':p, 'stderr':stderr}
pbar.animate(y*xmax+x + 1)
return self.result['alpha']
[docs] def refit(self,M=1, betamedx=True):
"""Advised to run after Alpha.fit(). Re-runs the fit with the
limitation that the laser_power * median(self.result['beta']) < M,
(i.e., restricts the used values of laser power and ignores the large
ones) to force obeying the condition of *weak* saturation.
Args:
**M**: (*float*) the upper limit of fitting range, i.e. only values with beta*pwr < M are permitted
**betamedx**: (*bool*) determines the method to extract beta. If True (default), median of all valid betas is used. If False, the beta on the particular pixel is used.
"""
self.M = M
self.betamedx=betamedx
zmax, ymax, xmax = self.data.shape
pbar = progressbar.ProgressBar(ymax*xmax)
self.result['delta_beta'] = self.delta_beta_map()
if self.betamedx:
valid_beta = (self.result['delta_beta'] > 0) * (self.result['delta_beta'] < 100)
betamed = np.median(self.result['beta'][valid_beta])
for y in xrange(ymax):
for x in xrange(xmax):
pwr = self.format_x(y, x)
pwrarr = np.array(pwr)
#mask_pwr = pwrarr*self.result['beta'][y,x] < 0.4
if not self.betamedx:
betamed = self.result['beta'][y,x]
mask_pwr = pwrarr*betamed < M
# do not recalculate if there are too few points left
if sum(mask_pwr) < 5:
self.result['alpha'][y, x] = np.nan
self.result['beta'][y,x] = np.nan
self.result['delta_alpha'][y,x] = np.nan
self.raw[(y, x)] = {'slope':np.nan, 'intercept': np.nan,
'r':np.nan, 'p':np.nan, 'stderr':np.nan}
else:
pwrarr = pwrarr[mask_pwr]
maxpwr = max(pwrarr)
pwrarr /= maxpwr ##rescale for better fit
data_to_fit, mask = self.invert(self.data[mask_pwr,y,x])
pwr_to_fit = pwrarr[mask]**(-1)
slope, intercept, r, p, stderr = linregress(pwr_to_fit, data_to_fit)
pwrarr *= maxpwr ## undo rescaling
self.result['alpha'][y, x] = 1./(slope*maxpwr)
self.result['beta'][y, x] = intercept * self.result['alpha'][y, x]
self.result['delta_alpha'][y,x] = stderr/(slope**2*maxpwr)
self.raw[(y, x)] = {'slope':slope, 'intercept': intercept, 'r':r, 'p':p, 'stderr':stderr}
pbar.animate(y*xmax+x + 1)
return self.result['alpha']
[docs] def show_refit(self, y, x, save_to_file=None,draw = False, betax=True, show=True):
"""
Similar to show fit in the parent Crossfit class, but highlights
the values really used in refit() method.
Args:
**y, x**: pixel coordinates
**betax**: (*bool*) if False, the x-axis of the plot is in the units of laser power. If True (default), the x-axis is multiplied by beta.
**save_to_file**: If a string is
given, it represents the name of the file the figure will be saved to. If None (default), no file is written.
**xlabel, ylabel**: (*strings*) passed directly to matplotlib
**draw**: if True, matplotlib.draw() and clear() are called. May be used in a loop to create an animation.
**show**: if True (default), matplotlib.show() is called.
Returns:
matplotlib.figure object.
"""
pwr = self.format_x(y, x)
valid_beta = (self.result['beta'] > 0) * (self.result['delta_beta'] < 100)
if self.betamedx:
betamed = np.median(self.result['beta'][valid_beta])
else:
betamed = self.result['beta'][y,x]
#mask_pwr = pwr*self.result['beta'][y,x] < 0.4
mask_pwr = pwr * betamed < self.M
if betax:
pwr *= betamed
pwr_used = pwr[mask_pwr]
pwrtpl = self.independent_var_to_plot(y, x)
fig, ax = plt.subplots(1, 1)
ax.set_title('used: ' + str(sum(mask_pwr)))
ax.plot(pwr, self.data[:,y,x], '+')
#print pwr_used
ax.plot(pwr_used, self.data[mask_pwr, y, x], 'o', linewidth=2)
p = self.get_params(y, x)
plt.plot(pwrtpl, self.func(pwrtpl/betamed, *p))
if betax:
pwr /= betamed
if save_to_file:
fig.savefig(save_to_file)
elif draw:
plt.draw()
plt.pause(0.001)
plt.gca().clear()
elif show:
plt.show()
return fig
[docs] def delta_beta(self, y, x):
"""Calculate and retun the estimate of standard deviation of the beta parameter
from the fit at pixel given by y,x coordinate.
"""
pwr = self.format_x(y, x)
var_pwr = np.var(pwr)
var_fl = np.var(self.data[:,y, x])
n = len(pwr)
return var_fl * np.sqrt(n) / (n-2) * np.sqrt(1 + np.mean(pwr)**2/var_pwr)
[docs] def beta_Imax(self):
"""Calculate and return the 2D map of values beta(y,x) *
max(laser_power(y,x)). Served for checking the condition of
weak saturation.
"""
zmax, ymax, xmax = self.data.shape
beta_Imax = np.zeros((ymax, xmax))
try:
beta = self.result['beta']
except KeyError:
msg = 'Alphas.beta_Imax: no result for \'beta\' found!'
warnings.warn()
return None
for y in xrange(ymax):
for x in xrange(xmax):
pwr = self.format_x(y, x)
maxpwr = max(pwr)
beta_Imax[y, x] = maxpwr * beta[y,x]
return beta_Imax
[docs] def print_result(self, y, x):
"""Returns a formatted string with the values of alpha +/-
delta_alpha and beta +/- delta_beta at pixel given by
coordinates y,x
"""
ret = 'alpha'
ret += '\t{:.4e}'.format(self.result['alpha'][y, x])
ret += '\t+/-\t'
ret += '\t{:.4e}'.format(self.result['delta_alpha'][y,x])
#ret += '\t{:.4e}'.format(self.delta_alpha(y, x))
ret += '\n'
ret += 'beta'
ret += '\t{:.4e}'.format(self.result['beta'][y,x])
ret += '\t+/-\t'
ret += '\t{:.4e}'.format(self.delta_beta(y,x))
return ret
[docs] def delta_beta_map(self, relative=True):
"""Runs delta_beta for each pixel and returns a 2D map of standard deviations.
Args:
**relative**: (*bool*) defaults to True, in which case returns the
relative standard deviation in percent. Otherwise the absolute values
in units of inverted laser power are returned.
"""
ymax, xmax = self.result['beta'].shape
ret = np.zeros((ymax, xmax))
for y in xrange(ymax):
for x in xrange(xmax):
ret[y, x] = self.delta_beta(y, x)
if relative:
ret[y, x] /= self.result['beta'][y, x] / 100
return ret
[docs]class Tau(Crossfit):
"""Class for working with planar LIF measurements of lfuorescence vs. delay after the laser pulse.
Attributes:
**func**: a function used to fit the measured dependent. Defaults to
:math:`A {\\rm e}^{(-t/ \\tau)} + C`.
**dfunc_dict**: dictionary of derivatives of func with respect to its
parameters. The keys should be 'amplitude', 'tau' and 'offset' and the
values should be the respective functions taking the arguments
(time_axis, amplitude, tau) and returning an array with the same shape
as the time_axis.
**params**: initial guess of the respective values. If left to None,
the default values of amplitude = 1e3, tau = 10 and offset = 1e3
are filled from start to the self.result arrays
**result**: dictionary with keys 'amplitude', 'tau', 'offset' (fit
results), 'damplitude', 'dtau', 'doffset' (standard error
estimates), ('cov', (y, x)) storing "covariance matrices" returned
by scipy.optimize.curve_fit and 'cov_order' storing the order of
the parameters used to calculate the "covariance matrices".
"""
def __init__(self, triDData, time, params=None, func=None, dfunc_dict=None, info={}):
Crossfit.__init__(self, triDData, time, func=func, info=info)
if self.func is None:
self.func = default_func_tau
self.dfunc_dict=dfunc_dict
if self.dfunc_dict is None:
self.dfunc_dict = {'amplitude':func_damplitude,
'tau':func_dtau,
'offset':func_doffset}
zmax, ymax, xmax = self.data.shape
self.result['amplitude'] = np.empty((ymax, xmax))
self.result['tau'] = np.empty((ymax, xmax))
self.result['offset'] = np.empty((ymax, xmax))
self.params = ['amplitude', 'tau', 'offset']
if params is None:
ampl = 1e3
tau = 10
offset = 1e3
else:
ampl = params[0]
tau = params[1]
offset = params[2]
self.result['amplitude'][:,:] = ampl
self.result['tau'][:,:] = tau
self.result['offset'][:,:] = offset
self.result['damplitude'] = np.zeros((ymax, xmax))
self.result['dtau'] = np.zeros((ymax, xmax))
self.result['doffset'] = np.zeros((ymax, xmax))
def func_from_keys(self, keys, values, y, x):
#print y, x
for key, p in zip(keys,values):
self.result[key][y, x] = p
return self.func(self.independent_variable,
self.result['amplitude'][y, x],
self.result['tau'][y, x],
self.result['offset'][y, x])
def dfunc_from_keys(self, keys, values, y, x):
for key, p in zip(keys,values):
self.result[key][y, x] = p
ret = []
params = (self.independent_variable, self.result['amplitude'][y, x],
self.result['tau'][y, x])
for key in keys:
ret.append(self.dfunc_dict[key](*params))
return ret
def get_params(self, y, x):
return [self.result['amplitude'][y, x], self.result['tau'][y, x], self.result['offset'][y, x]]
[docs]class FluoVsPwr(Crossfit):
"""
Replacement of ALpha with more rigorous function. Is not very practical for real-world fitting.
"""
def __init__(self, triDData, laser_pwr, params=None, info={}):
Crossfit.__init__(self, triDData, laser_pwr, info=info)
zmax, ymax, xmax = self.data.shape
if self.func is None:
self.func = default_func_FluoVsPwr
self.result['alpha'] = np.empty((ymax, xmax))
self.result['B'] = np.empty((ymax, xmax))
self.result['t'] = np.empty((ymax, xmax))
self.params = ['alpha', 'B', 't']
self.dfunc_dict = {'alpha':func_FluoVsPwr_dalpha,
'B':func_FluoVsPwr_dB,
't':func_FluoVsPwr_dt}
if params is None:
alpha = 1e9
B = 1e5
t = 1
else:
alpha = params[0]
B = params[1]
t = params[2]
self.result['alpha'][:,:] = alpha
self.result['B'][:,:] = B
self.result['t'][:,:] = t
self.result['dalpha'] = np.zeros((ymax, xmax))
self.result['dB'] = np.zeros((ymax, xmax))
self.result['dt'] = np.zeros((ymax, xmax))
def func_from_keys(self, keys, values, y, x):
for key, p in zip(keys,values):
self.result[key][y, x] = p
return self.func(self.format_x(y,x),
self.result['alpha'][y, x],
self.result['B'][y, x],
self.result['t'][y, x])
def dfunc_from_keys(self, keys, values, y, x):
for key, p in zip(keys,values):
self.result[key][y, x] = p
ret = []
params = (self.format_x(y, x), self.result['alpha'][y, x],
self.result['B'][y, x], self.result['t'][y, x])
for key in keys:
ret.append(self.dfunc_dict[key](*params))
return ret
def get_params(self, y, x):
return [self.result['alpha'][y, x], self.result['B'][y, x], self.result['t'][y, x]]
class Temperature(Alpha):
def __init__(self, dict_of_Alphas, B_coefs, smooth = False, **kwargs):
"""Args:
dict_of_Alphas: (*dictionary*) of Alpha objects as values and a
tuple specifying the rotational line used for excitation, e.g. ('R1',2)
B_coefs: a dictionary of B coeficients, as returned by
rot_en.read_lifbase_table(). I.e., dictionary of columns with keys
'N P1 P2 Q1 Q2 R1 R2 O12 Q12 P12 R21 Q21 S21', where N is rotational
number.
**kwrgs: passed directly to smooth_alphas and clear_outliers
"""
self.Bs=B_coefs
self.dict_of_Alphas = dict_of_Alphas
self.Es = {}
self.ls = None
self.process_alphas(smooth = smooth, **kwargs)
evals = np.array(self.Es.values())
evals.sort()
self.strip_alphas()
self.smooth_options = {}
Alpha.__init__(self, self.ls, np.array(evals), func=linfunc)
def smooth_alphas(self, **kwargs):
"""
Call a smoothing algoritm from scipy.ndimage.
**kwargs:
**perc**: (*float*) if given, percentile_filter(alpha, perc) will be called
**gaussian**: if True, clear_outliers and gaussian_filter will be called.
"""
self.smooth_options = kwargs
perc = kwargs.pop("perc", None)
gauss = kwargs.pop('gaussian', False)
for key in self.dict_of_Alphas:
alpha = self.dict_of_Alphas[key]
if perc is None and not gauss:
alpha.result['alpha_smooth'] = median_filter(alpha.result['alpha'], **kwargs)
elif perc is None and gauss:
self.clear_outliers(alpha, **kwargs)
alpha.result['alpha_smooth'] = gaussian_filter(alpha.result['alpha'], **kwargs)
else:
alpha.result['alpha_smooth'] = percentile_filter(alpha.result['alpha'], perc, **kwargs)
@staticmethod
def clear_outliers(alpha, **kwargs):
"""
Replace outlying values by minimal/maximal allowed value prior to calculating gaussian filter.
**kwargs:
**perc**: used to calculate the allowed maximum as a percentile of given data. Defaults to 99.9
**min**: allowed minimum. Defaults to 0.
"""
perc = kwargs.pop('perc', 99.9)
low = kwargs.pop('min', 0)
alpha.result['alpha'][alpha.result['alpha'] < low] = low
lim = np.percentile(alpha.result['alpha'], perc)
alpha.result['alpha'][alpha.result['alpha'] > lim ]= lim
def strip_alphas(self):
"""
Free some memory by discarding data no longer needed.
"""
##smazat z alf to, co neni treba drzet
for key in self.dict_of_Alphas:
self.dict_of_Alphas[key].data = None
self.dict_of_Alphas[key].raw = None
def get_params(self, y, x):
return self.raw[(y, x)]['slope'], self.raw[(y,x)]['intercept']
def extract_Ns(self):
Ns = []
for key in self.dict_of_Alphas.keys():
Ns.append(key[1])
return np.array(Ns, dtype=int)
def get_rot_en(self):
for key in self.dict_of_Alphas.keys():
N = key[1]
manifold = int(key[0][-1])
#self.Es[(N, manifold)] = rot_en.E_N_Pi(N, manifold)
self.Es[key] = rot_en.E_N_Pi(N, manifold) / e #energy in eV
return self.Es
def process_alphas(self, smooth=False, **kwargs):
if smooth:
self.smooth_alphas(**kwargs)
alpha_key = 'alpha_smooth'
else:
alpha_key = 'alpha'
y, x = self.dict_of_Alphas.values()[0].result['alpha'].shape
ls = np.zeros((len(self.dict_of_Alphas), y, x) )
for i,key in enumerate(self.order_keys_by_EN()):
g = 2*(key[1] + 0.5*(-1)**(1+int(key[0][-1]))) + 1
ls[i] = np.log(self.dict_of_Alphas[key].result[alpha_key] /
(g * self.Bs[key[0]][self.Bs['N']==key[1]]))
self.ls = ls
return self.ls
def order_keys_by_EN(self):
Es = self.get_rot_en()
Evals = np.array(Es.values())
idcs = Evals.argsort()
keys = Es.keys()
keys_sorted=[]
for i in idcs:
key = keys[i]
keys_sorted.append(key)
return keys_sorted
def fit(self):
zmax, ymax, xmax = self.data.shape
pbar = progressbar.ProgressBar(ymax*xmax)
self.result['T'] = np.zeros((ymax, xmax))
self.result['dT'] = np.zeros((ymax, xmax))
evals = self.Es.values()
evals.sort()
evals = np.array(evals)
for y in xrange(ymax):
for x in xrange(xmax):
mask_to_fit = np.isnan(self.ls[:,y,x])
slope, intercept, r, p, stderr = linregress(evals[~mask_to_fit], self.ls[~mask_to_fit,y,x])
self.result['T'][y,x] = -1./(slope*k)*e
self.result['dT'][y,x] = stderr / (k*slope**2) * e
self.raw[(y, x)] = {'slope':slope, 'intercept': intercept, 'r':r, 'p':p, 'stderr':stderr}
pbar.animate(y*xmax+x + 1)
return self.result['T']
[docs]class WaitingCorr(object):
"""Calculate correction for the fact that the detection begins after
the laser pulse. The fluorescence temporal shape is calculated as
a convolution of single-exponential and a*t**b * exp(-c*t).
"""
def __init__(self, tau_map, delta_tau_map=None, **kwargs):
self.taus = tau_map
self.laser_b = kwargs.pop('laser_b', 12.6)
self.laser_c = kwargs.pop('laser_c', 2.44)
self.t_max = kwargs.pop('t_max',500)
self.T_start_ns = kwargs.pop('T_start_ns',19)
self.t_points = kwargs.pop('t_points',1000)
self.dtaus = delta_tau_map
if self.dtaus is None:
self.dtaus = np.zeros(self.taus.shape)
self.time = np.linspace(0,self.t_max,self.t_points)
self.result = np.ones(self.taus.shape)
def laser_pulse(self):
return (self.time)**self.laser_b * np.exp(-self.laser_c*self.time)
def fluorescence(self, tau):
ex = np.exp(-self.time/tau)
ret = np.convolve(ex, self.laser_pulse(), mode='full')
return ret[:len(self.time)]
def find_tstart_idx(self):
return np.where(self.time>self.T_start_ns)[0][0]
def corr(self, tau):
total = np.trapz(self.fluorescence(tau), self.time)
tstart_idx = self.find_tstart_idx()
after = np.trapz(self.fluorescence(tau)[tstart_idx:],
self.time[tstart_idx:])
return total/after
def process(self):
ymax, xmax = self.taus.shape
pbar = progressbar.ProgressBar(ymax*xmax)
for y in xrange(ymax):
for x in xrange(xmax):
if (self.dtaus[y,x] < 50) and (self.taus[y,x] >
self.T_start_ns/6.) and (self.taus[y,x] <
self.T_start_ns*1e3):
self.result[y, x] = self.corr(self.taus[y,x])
pbar.animate(y*xmax+x + 1)
return self.result
def show(self, y, x, save_to_file=None, xlabel='time [ns]',
ylabel = 'modelled fluorescence [arb. units]'):
fluo = self.fluorescence(self.taus[y, x])
pulse = self.laser_pulse()
fig, ax = plt.subplots(1,1)
ax.plot(self.time, fluo/np.max(fluo))
ax.plot(self.time, pulse/np.max(pulse))
ax.vlines(self.T_start_ns, 0, 1)
ax.set_xlabel(xlabel)
ax.set_ylabel(ylabel)
if save_to_file:
fig.savefig(save_to_file)
else:
plt.show()
return fig
def save(self, filename):
with open(filename, 'wb') as output:
pickle.dump(self, output, pickle.HIGHEST_PROTOCOL)
@staticmethod
def load(filename):
with open(filename, 'rb') as inp:
return_val = pickle.load(inp)
return return_val