Skip to content

GitLab

Commit 831cbad2 authored by Unknown's avatar Unknown
Browse files

More PEP8 cleanup

parent 66229074
Expand all Hide whitespace changes
Inline Side-by-side
pycroscopy/analysis/model.py
View file @ 831cbad2
......@@ -37,6 +37,7 @@ class Model(object):
None
"""
def __init__(self, h5_main, variables=['Frequency'], parallel=True):
"""
For now, we assume that the guess dataset has not been generated for this dataset but we will relax this
......@@ -90,7 +91,7 @@ class Model(object):
if self._maxCpus == 1:
self._parallel = False
self._maxMemoryMB = getAvailableMem() / 1024**2 # in Mb
self._maxMemoryMB = getAvailableMem() / 1024 ** 2 # in Mb
self._maxDataChunk = int(self._maxMemoryMB / self._maxCpus)
......@@ -161,7 +162,6 @@ class Model(object):
print('Finished reading all data!')
self.data = None
def _get_guess_chunk(self):
"""
Returns a chunk of guess dataset corresponding to the main dataset.
......@@ -252,11 +252,11 @@ class Model(object):
"""
warn('Please override the _create_fit_datasets specific to your model')
self.fit = None # replace with actual h5 dataset
self.fit = None # replace with actual h5 dataset
pass
def do_guess(self, processors=None, strategy='wavelet_peaks',
options={"peak_widths": np.array([10, 200]), "peak_step":20}):
options={"peak_widths": np.array([10, 200]), "peak_step": 20}):
"""
Parameters
......
pycroscopy/analysis/utils/atom_finding.py
View file @ 831cbad2
......@@ -34,6 +34,7 @@ atom_dtype = np.dtype({'names': ['x', 'y', 'type'],
atom_coeff_dtype = np.dtype({'names': ['Amplitude', 'x', 'y', 'Sigma'],
'formats': [np.float32, np.float32, np.float32, np.float32]})
def multi_gauss_surface_fit(coef_mat, s_mat):
"""
Evaluates the provided coefficients for N gaussian peaks to generate a 2D matrix
......@@ -233,7 +234,7 @@ def fit_atom_positions_parallel(parm_dict, fitting_parms, num_cores=None):
all_atom_guesses = parm_dict['atom_pos_guess']
t_start = tm.time()
num_cores = recommendCores(all_atom_guesses.shape[0], requested_cores=num_cores, lengthy_computation=False)
if num_cores>1:
if num_cores > 1:
pool = mp.Pool(processes=num_cores)
parm_list = itt.izip(range(all_atom_guesses.shape[0]), itt.repeat(parm_dict), itt.repeat(fitting_parms))
chunk = int(all_atom_guesses.shape[0] / num_cores)
......@@ -316,7 +317,7 @@ def fit_atom_positions_dset(h5_grp, fitting_parms=None, num_cores=None):
fitting_results = fit_atom_positions_parallel(parm_dict, fitting_parms, num_cores=num_cores)
# Make datasets to write back to file:
guess_parms = np.zeros(shape=(num_atoms, num_nearest_neighbors+1), dtype=atom_coeff_dtype)
guess_parms = np.zeros(shape=(num_atoms, num_nearest_neighbors + 1), dtype=atom_coeff_dtype)
fit_parms = np.zeros(shape=guess_parms.shape, dtype=guess_parms.dtype)
for atom_ind, single_atom_results in enumerate(fitting_results):
......
pycroscopy/analysis/utils/atom_finding_general_gaussian.py
View file @ 831cbad2
......@@ -44,9 +44,9 @@ def do_fit(single_parm):
coef_fit_mat = np.reshape(plsq.x, (-1, 7))
#if verbose:
# if verbose:
# return coef_guess_mat, lb_mat, ub_mat, coef_fit_mat, fit_region, s_mat, plsq
#else:
# else:
return coef_fit_mat
......@@ -422,21 +422,21 @@ class Gauss_Fit(object):
ub_background.append(item + item * movement_allowance)
# Set up upper and lower bounds:
lb_mat = [lb_a, # amplitude
coef_guess_mat[:, 1] - position_range, # x position
coef_guess_mat[:, 2] - position_range, # y position
lb_mat = [lb_a, # amplitude
coef_guess_mat[:, 1] - position_range, # x position
coef_guess_mat[:, 2] - position_range, # y position
[np.max([0, value - value * movement_allowance]) for value in coef_guess_mat[:, 3]], # sigma x
[np.max([0, value - value * movement_allowance]) for value in coef_guess_mat[:, 4]], # sigma y
coef_guess_mat[:, 5] - 2 * 3.14159, # theta
lb_background] # background
coef_guess_mat[:, 5] - 2 * 3.14159, # theta
lb_background] # background
ub_mat = [ub_a, # amplitude
coef_guess_mat[:, 1] + position_range, # x position
coef_guess_mat[:, 2] + position_range, # y position
ub_mat = [ub_a, # amplitude
coef_guess_mat[:, 1] + position_range, # x position
coef_guess_mat[:, 2] + position_range, # y position
coef_guess_mat[:, 3] + coef_guess_mat[:, 3] * movement_allowance, # sigma x
coef_guess_mat[:, 4] + coef_guess_mat[:, 4] * movement_allowance, # sigma y
coef_guess_mat[:, 5] + 2 * 3.14159, # theta
ub_background] # background
coef_guess_mat[:, 5] + 2 * 3.14159, # theta
ub_background] # background
lb_mat = np.transpose(lb_mat)
ub_mat = np.transpose(ub_mat)
......@@ -608,46 +608,3 @@ class Gauss_Fit(object):
fig.show()
return self.motif_converged_dataset
# if __name__=='__main__':
# file_name = r"C:\Users\o2d\Documents\pycroscopy\\test_scripts\\testing_gauss_fit\image 04.h5"
# folder_path, file_path = os.path.split(file_name)
#
# file_base_name, file_extension = file_name.rsplit('.')
# h5_file = h5py.File(file_name, mode='r+')
# # look at the data tree in the h5
# '''
# # define a small function called 'print_tree' to look at the folder tree structure
# def print_tree(parent):
# print(parent.name)
# if isinstance(parent, h5py.Group):
# for child in parent:
# print_tree(parent[child])
# '''
#
# #print('Datasets and datagroups within the file:')
# file_handle = h5_file
# #print_tree(file_handle)
#
# cropped_clean_image = h5_file['/Measurement_000/Channel_000/Raw_Data-Windowing_000/Image_Windows-SVD_000/U-Cluster_000/Labels-Atom_Finding_000/Cropped_Clean_Image']
# atom_grp = cropped_clean_image.parent
# guess_params = atom_grp['Guess_Positions']
#
# num_nearest_neighbors = 4
# psf_width = atom_grp.attrs['psf_width']
# win_size = atom_grp.attrs['motif_win_size']
#
# fitting_parms = {'fit_region_size': win_size * 0.5, # region to consider when fitting
# 'num_nearest_neighbors': num_nearest_neighbors,
# 'sigma_guess': 3, # starting guess for gaussian standard deviation
# 'position_range': win_size / 4,# range that the fitted position can go from initial guess position[pixels]
# 'max_function_evals': 100,
# 'fitting_tolerance': 1E-4,
# 'symmetric': True,
# 'background': True,
# 'movement_allowance': 5.0} # percent of movement allowed (on some parameters)
#
# foo = Gauss_Fit(atom_grp, fitting_parms)
#
#
pycroscopy/analysis/utils/be_loop.py
View file @ 831cbad2
......@@ -57,25 +57,25 @@ def calculate_loop_centroid(vdc, loop_vals):
vdc = np.squeeze(np.array(vdc))
num_steps = vdc.size
x_vals = np.zeros(num_steps - 1)
y_vals = np.zeros(num_steps - 1)
area_vals = np.zeros(num_steps - 1)
for index in range(num_steps - 1):
x_i = vdc[index]
x_i1 = vdc[index + 1]
y_i = loop_vals[index]
y_i1 = loop_vals[index + 1]
x_vals[index] = (x_i + x_i1)*(x_i*y_i1 - x_i1*y_i)
y_vals[index] = (y_i + y_i1)*(x_i*y_i1 - x_i1*y_i)
area_vals[index] = (x_i*y_i1 - x_i1*y_i)
x_vals[index] = (x_i + x_i1) * (x_i * y_i1 - x_i1 * y_i)
y_vals[index] = (y_i + y_i1) * (x_i * y_i1 - x_i1 * y_i)
area_vals[index] = (x_i * y_i1 - x_i1 * y_i)
area = 0.50 * np.sum(area_vals)
cent_x = (1.0/(6.0*area)) * np.sum(x_vals)
cent_y = (1.0/(6.0*area)) * np.sum(y_vals)
cent_x = (1.0 / (6.0 * area)) * np.sum(x_vals)
cent_y = (1.0 / (6.0 * area)) * np.sum(y_vals)
return (cent_x, cent_y), area
......@@ -185,7 +185,7 @@ def projectLoop(vdc, amp_vec, phase_vec):
pdist_vals = np.zeros(shape=len(xdat_fit))
for point in range(len(xdat_fit)):
pdist_vals[point] = np.linalg.norm(np.array([0, 0]) -
np.array([xdat_fit[point], ydat_fit[point]]))
np.array([xdat_fit[point], ydat_fit[point]]))
min_point_ind = np.where(pdist_vals == pdist_vals.min())[0][0]
min_point = np.array([xdat_fit[min_point_ind], ydat_fit[min_point_ind]])
......@@ -238,7 +238,8 @@ def projectLoop(vdc, amp_vec, phase_vec):
'Centroid': centroid, 'Geometric Area': geo_area} # Dictionary of Results from projecting
return results
###############################################################################
......@@ -262,19 +263,19 @@ def loop_fit_function(vdc, coef_vec):
a = coef_vec[:5]
b = coef_vec[5:]
d = 1000
v1 = np.asarray(vdc[:int(len(vdc) / 2)])
v2 = np.asarray(vdc[int(len(vdc) / 2):])
g1 = (b[1]-b[0])/2*(erf((v1-a[2])*d)+1)+b[0]
g2 = (b[3]-b[2])/2*(erf((v2-a[3])*d)+1)+b[2]
y1 = (g1 * erf((v1-a[2])/g1) + b[0])/(b[0]+b[1])
y2 = (g2 * erf((v2-a[3])/g2) + b[2])/(b[2]+b[3])
g1 = (b[1] - b[0]) / 2 * (erf((v1 - a[2]) * d) + 1) + b[0]
g2 = (b[3] - b[2]) / 2 * (erf((v2 - a[3]) * d) + 1) + b[2]
y1 = (g1 * erf((v1 - a[2]) / g1) + b[0]) / (b[0] + b[1])
y2 = (g2 * erf((v2 - a[3]) / g2) + b[2]) / (b[2] + b[3])
f1 = a[0] + a[1] * y1 + a[4] * v1
f2 = a[0] + a[1] * y2 + a[4] * v2
f1 = a[0] + a[1]*y1 + a[4]*v1
f2 = a[0] + a[1]*y2 + a[4]*v2
loop_eval = np.hstack((f1, f2))
return loop_eval
......@@ -313,7 +314,7 @@ def loop_fit_jacobian(vdc, coef_vec):
g2 = (b[3] - b[2]) / 2 * (erf((V2 - a[3]) * d) + 1) + b[2]
# Some useful fractions
oosqpi = 1.0/np.sqrt(np.pi)
oosqpi = 1.0 / np.sqrt(np.pi)
oob01 = 1.0 / (b[0] + b[1])
oob23 = 1.0 / (b[2] + b[3])
......@@ -325,16 +326,16 @@ def loop_fit_jacobian(vdc, coef_vec):
dg2b2 = -0.5 * (erf(V2 - a[3]) * d + 1)
dg2b3 = -dg2b2
Y1 = oob01 * (g1 * erf((V1-a[2])/g1) + b[0])
Y2 = oob23 * (g2 * erf((V2-a[3])/g2) + b[2])
Y1 = oob01 * (g1 * erf((V1 - a[2]) / g1) + b[0])
Y2 = oob23 * (g2 * erf((V2 - a[3]) / g2) + b[2])
# Derivatives of Y1 and Y2
dY1g1 = oob01 * (erf((V1 - a[2]) / g1) - 2 * oosqpi * np.exp(-(V1 - a[2])**2 / g1**2))
dY2g2 = oob23 * (erf((V2 - a[4]) / g2) - 2 * oosqpi * np.exp(-(V2 - a[4])**2 / g2**2))
dY1g1 = oob01 * (erf((V1 - a[2]) / g1) - 2 * oosqpi * np.exp(-(V1 - a[2]) ** 2 / g1 ** 2))
dY2g2 = oob23 * (erf((V2 - a[4]) / g2) - 2 * oosqpi * np.exp(-(V2 - a[4]) ** 2 / g2 ** 2))
dY1a2 = -2 * oosqpi * oob01 * np.exp(-(V1 - a[2])**2 / g1**2)
dY1a2 = -2 * oosqpi * oob01 * np.exp(-(V1 - a[2]) ** 2 / g1 ** 2)
dY2a3 = -2 * oosqpi * oob01 * np.exp(-(V2 - a[3])**2 / g2**2)
dY2a3 = -2 * oosqpi * oob01 * np.exp(-(V2 - a[3]) ** 2 / g2 ** 2)
dY1b0 = oob01 * (1 - Y1)
dY1b1 = -oob01 * Y1
......@@ -387,7 +388,7 @@ def loop_fit_jacobian(vdc, coef_vec):
return J
###############################################################################
......@@ -409,50 +410,50 @@ def get_switching_coefs(loop_centroid, loop_coef_vec):
anv = loop_coef_vec[0:5]
bnv = loop_coef_vec[5:]
nuc_threshold = .97
lc_x = loop_centroid[0]
# Upper Branch
B = bnv[2]
nuc_v01a = B*erfinv((nuc_threshold*bnv[2] + nuc_threshold*bnv[3] - bnv[2]) /B) - abs(anv[2])
nuc_v01a = B * erfinv((nuc_threshold * bnv[2] + nuc_threshold * bnv[3] - bnv[2]) / B) - abs(anv[2])
B = bnv[3]
nuc_v01b = B*erfinv((nuc_threshold*bnv[2] + nuc_threshold*bnv[3] - bnv[2]) /B) -abs(anv[2])
nuc_v01b = B * erfinv((nuc_threshold * bnv[2] + nuc_threshold * bnv[3] - bnv[2]) / B) - abs(anv[2])
# Choose the voltage which is closer to the centroid x value.
if abs(nuc_v01a-lc_x) <= abs(nuc_v01b-lc_x):
if abs(nuc_v01a - lc_x) <= abs(nuc_v01b - lc_x):
nuc_1 = nuc_v01a
else:
nuc_1 = nuc_v01b
if np.isinf(nuc_1):
nuc_1 = 0
# Lower Branch
B = bnv[0]
nuc_v02a = B*erfinv(((1-nuc_threshold)*bnv[0] + (1-nuc_threshold)*bnv[1] - bnv[0])/B) + abs(anv[3])
nuc_v02a = B * erfinv(((1 - nuc_threshold) * bnv[0] + (1 - nuc_threshold) * bnv[1] - bnv[0]) / B) + abs(anv[3])
B = bnv[1]
nuc_v02b = B*erfinv(((1-nuc_threshold)*bnv[0] + (1-nuc_threshold)*bnv[1] - bnv[0])/B) + abs(anv[3])
if abs(nuc_v02a-lc_x) <= abs(nuc_v02b-lc_x):
nuc_v02b = B * erfinv(((1 - nuc_threshold) * bnv[0] + (1 - nuc_threshold) * bnv[1] - bnv[0]) / B) + abs(anv[3])
if abs(nuc_v02a - lc_x) <= abs(nuc_v02b - lc_x):
nuc_2 = nuc_v02a
else:
nuc_2 = nuc_v02b
if np.isinf(nuc_2):
nuc_2 = 0
# calculate switching parameters
switching_coef_vec = np.zeros(shape=(9, 1))
switching_coef_vec[0] = loop_coef_vec[2]
switching_coef_vec[1] = loop_coef_vec[3]
switching_coef_vec[2] = (loop_coef_vec[3]+loop_coef_vec[2])/2
switching_coef_vec[3] = loop_coef_vec[0]+loop_coef_vec[1]
switching_coef_vec[2] = (loop_coef_vec[3] + loop_coef_vec[2]) / 2
switching_coef_vec[3] = loop_coef_vec[0] + loop_coef_vec[1]
switching_coef_vec[4] = loop_coef_vec[0]
switching_coef_vec[5] = loop_coef_vec[1]
switching_coef_vec[6] = abs((loop_coef_vec[2]-loop_coef_vec[3])*abs(loop_coef_vec[1]))
switching_coef_vec[6] = abs((loop_coef_vec[2] - loop_coef_vec[3]) * abs(loop_coef_vec[1]))
switching_coef_vec[7] = nuc_1
switching_coef_vec[8] = nuc_2
......@@ -460,10 +461,11 @@ def get_switching_coefs(loop_centroid, loop_coef_vec):
'Imprint': switching_coef_vec[2], 'R+': switching_coef_vec[3],
'R-': switching_coef_vec[4], 'Switchable Polarization': switching_coef_vec[5],
'Work of Switching': switching_coef_vec[6],
'Nucleation Bias 1': nuc_1, 'Nucleation Bias 2': nuc_2}
'Nucleation Bias 1': nuc_1, 'Nucleation Bias 2': nuc_2}
return switching_labels_dict
##############################################################################
......@@ -497,7 +499,7 @@ def calc_switching_coef_vec(loop_coef_vec, nuc_threshold):
nuc_v02a = bnv[0] * erfinv(((1 - nuc_threshold) * bnv[0] + (1 - nuc_threshold) * bnv[1] - bnv[0]) / bnv[0]) + anv[2]
nuc_v02b = bnv[1] * erfinv(((1 - nuc_threshold) * bnv[0] + (1 - nuc_threshold) * bnv[1] - bnv[0]) / bnv[1]) + anv[2]
nuc_v02 = np.where(np.isfinite(nuc_v02a),nuc_v02a,
nuc_v02 = np.where(np.isfinite(nuc_v02a), nuc_v02a,
np.where(np.isfinite(nuc_v02b), nuc_v02b, -1E-10))
switching_coef_vec['V+'] = loop_coef_vec[:, 3]
......@@ -608,8 +610,8 @@ def generate_guess(vdc, pr_vec, show_plots=False):
y_intersections = []
for pair in range(outline_1.shape[0]):
x_pt = find_intersection(outline_1[pair], outline_2[pair],
[geom_centroid[0], hull.min_bound[1]],
[geom_centroid[0], hull.max_bound[1]])
[geom_centroid[0], hull.min_bound[1]],
[geom_centroid[0], hull.max_bound[1]])
if type(x_pt) != type(None):
y_intersections.append(x_pt)
......@@ -620,8 +622,8 @@ def generate_guess(vdc, pr_vec, show_plots=False):
x_intersections = []
for pair in range(outline_1.shape[0]):
x_pt = find_intersection(outline_1[pair], outline_2[pair],
[hull.min_bound[0], geom_centroid[1]],
[hull.max_bound[0], geom_centroid[1]])
[hull.min_bound[0], geom_centroid[1]],
[hull.max_bound[0], geom_centroid[1]])
if type(x_pt) != type(None):
x_intersections.append(x_pt)
......@@ -636,8 +638,8 @@ def generate_guess(vdc, pr_vec, show_plots=False):
max_y_intercept = max(y_intersections[0][1], y_intersections[1][1])
if len(x_intersections) == 0:
min_x_intercept = min(vdc)/2.0
max_x_intercept = max(vdc)/2.0
min_x_intercept = min(vdc) / 2.0
max_x_intercept = max(vdc) / 2.0
else:
min_x_intercept = min(x_intersections[0][0], x_intersections[1][0])
max_x_intercept = max(x_intersections[0][0], x_intersections[1][0])
......@@ -672,6 +674,7 @@ def generate_guess(vdc, pr_vec, show_plots=False):
return init_guess_coef_vec
###############################################################################
......
pycroscopy/analysis/utils/be_sho.py
View file @ 831cbad2
......@@ -22,7 +22,7 @@ def SHOfunc(parms, w_vec):
Vector of frequency values
"""
return parms[0] * exp(1j * parms[3]) * parms[1] ** 2 / \
(w_vec ** 2 - 1j * w_vec * parms[1] / parms[2] - parms[1] ** 2)
(w_vec ** 2 - 1j * w_vec * parms[1] / parms[2] - parms[1] ** 2)
def SHOestimateGuess(w_vec, resp_vec, num_points=5):
......@@ -43,69 +43,71 @@ def SHOestimateGuess(w_vec, resp_vec, num_points=5):
retval : tuple
SHO fit parameters arranged as amplitude, frequency, quality factor, phase
"""
ii = np.argsort(abs(resp_vec))[::-1]
a_mat = np.array([])
e_vec = np.array([])
for c1 in range(num_points):
for c2 in range(c1+1,num_points):
for c2 in range(c1 + 1, num_points):
w1 = w_vec[ii[c1]]
w2 = w_vec[ii[c2]]
X1 = real(resp_vec[ii[c1]])
X2 = real(resp_vec[ii[c2]])
Y1 = imag(resp_vec[ii[c1]])
Y2 = imag(resp_vec[ii[c2]])
denom = (w1*(X1**2 - X1*X2 + Y1*(Y1 - Y2)) + w2*(-X1*X2 + X2**2 - Y1*Y2 + Y2**2))
denom = (w1 * (X1 ** 2 - X1 * X2 + Y1 * (Y1 - Y2)) + w2 * (-X1 * X2 + X2 ** 2 - Y1 * Y2 + Y2 ** 2))
if denom > 0:
a = ((w1**2 - w2**2)*(w1*X2*(X1**2 + Y1**2) - w2*X1*(X2**2 + Y2**2)))/denom
b = ((w1**2 - w2**2)*(w1*Y2*(X1**2 + Y1**2) - w2*Y1*(X2**2 + Y2**2)))/denom
c = ((w1**2 - w2**2)*(X2*Y1 - X1*Y2))/denom
d = (w1**3*(X1**2 + Y1**2) -
w1**2*w2*(X1*X2 + Y1*Y2) -
w1*w2**2*(X1*X2 + Y1*Y2) +
w2**3*(X2**2 + Y2**2))/denom
if d>0:
a = ((w1 ** 2 - w2 ** 2) * (w1 * X2 * (X1 ** 2 + Y1 ** 2) - w2 * X1 * (X2 ** 2 + Y2 ** 2))) / denom
b = ((w1 ** 2 - w2 ** 2) * (w1 * Y2 * (X1 ** 2 + Y1 ** 2) - w2 * Y1 * (X2 ** 2 + Y2 ** 2))) / denom
c = ((w1 ** 2 - w2 ** 2) * (X2 * Y1 - X1 * Y2)) / denom
d = (w1 ** 3 * (X1 ** 2 + Y1 ** 2) -
w1 ** 2 * w2 * (X1 * X2 + Y1 * Y2) -
w1 * w2 ** 2 * (X1 * X2 + Y1 * Y2) +
w2 ** 3 * (X2 ** 2 + Y2 ** 2)) / denom
if d > 0:
a_mat = append(a_mat, [a, b, c, d])
A_fit=abs(a+1j*b)/d
w0_fit=sqrt(d)
Q_fit =-sqrt(d)/c
phi_fit=arctan2(-b,-a)
H_fit = A_fit*w0_fit**2*exp(1j*phi_fit)/(w_vec ** 2 - 1j * w_vec * w0_fit / Q_fit - w0_fit ** 2)
A_fit = abs(a + 1j * b) / d
w0_fit = sqrt(d)
Q_fit = -sqrt(d) / c
phi_fit = arctan2(-b, -a)
H_fit = A_fit * w0_fit ** 2 * exp(1j * phi_fit) / (
w_vec ** 2 - 1j * w_vec * w0_fit / Q_fit - w0_fit ** 2)
e_vec = append(e_vec,
sum((real(H_fit) - real(resp_vec)) ** 2) +
sum((imag(H_fit) - imag(resp_vec)) ** 2))
if a_mat.size>0:
a_mat=a_mat.reshape(-1,4)
weight_vec = (1/e_vec)**4
w_sum=sum(weight_vec)
a_w = sum(weight_vec*a_mat[:,0])/w_sum
b_w = sum(weight_vec*a_mat[:,1])/w_sum
c_w = sum(weight_vec*a_mat[:,2])/w_sum
d_w = sum(weight_vec*a_mat[:,3])/w_sum
A_fit = abs(a_w+1j*b_w)/d_w
if a_mat.size > 0:
a_mat = a_mat.reshape(-1, 4)
weight_vec = (1 / e_vec) ** 4
w_sum = sum(weight_vec)
a_w = sum(weight_vec * a_mat[:, 0]) / w_sum
b_w = sum(weight_vec * a_mat[:, 1]) / w_sum
c_w = sum(weight_vec * a_mat[:, 2]) / w_sum
d_w = sum(weight_vec * a_mat[:, 3]) / w_sum
A_fit = abs(a_w + 1j * b_w) / d_w
w0_fit = sqrt(d_w)
Q_fit = -sqrt(d_w)/c_w
phi_fit = np.arctan2(-b_w,-a_w)
H_fit = A_fit*w0_fit**2*exp(1j*phi_fit)/(w_vec ** 2 - 1j * w_vec * w0_fit / Q_fit - w0_fit ** 2)
if np.std(abs(resp_vec))/np.std(abs(resp_vec-H_fit))<1.2 or w0_fit<np.min(w_vec) or w0_fit>np.max(w_vec):
p0=SHOfastGuess(w_vec, resp_vec)
else:
p0=np.array([A_fit, w0_fit, Q_fit, phi_fit])
Q_fit = -sqrt(d_w) / c_w
phi_fit = np.arctan2(-b_w, -a_w)
H_fit = A_fit * w0_fit ** 2 * exp(1j * phi_fit) / (w_vec ** 2 - 1j * w_vec * w0_fit / Q_fit - w0_fit ** 2)
if np.std(abs(resp_vec)) / np.std(abs(resp_vec - H_fit)) < 1.2 or w0_fit < np.min(w_vec) or w0_fit > np.max(
w_vec):
p0 = SHOfastGuess(w_vec, resp_vec)
else:
p0 = np.array([A_fit, w0_fit, Q_fit, phi_fit])
else:
p0=SHOfastGuess(w_vec, resp_vec)
p0 = SHOfastGuess(w_vec, resp_vec)
return p0
......@@ -128,9 +130,10 @@ def SHOfastGuess(w_vec, resp_vec, qual_factor=200):
SHO fit parameters arranged as [amplitude, frequency, quality factor, phase]
"""
amp_vec = abs(resp_vec)
i_max = int(len(resp_vec)/2)
i_max = int(len(resp_vec) / 2)
return np.array([np.mean(amp_vec) / qual_factor, w_vec[i_max], qual_factor, np.angle(resp_vec[i_max])])
def SHOlowerBound(w_vec):
"""
Provides the lower bound for the SHO fitting function
......
pycroscopy/io/microdata.py
View file @ 831cbad2
......@@ -19,7 +19,7 @@ class MicroData(object):
"""
Generic class that is extended by the MicroDataGroup and MicroDataset objects
"""
def __init__(self, name, parent):
'''
Parameters
......@@ -35,12 +35,14 @@ class MicroData(object):
self.parent = parent
self.indexed = False
class MicroDataGroup(MicroData):
"""
Holds data that will be converted to a h5.Group by io.ioHDF5
Note that it can also hold information (e.g. attributes) of an h5.File.