Commit 181619f5 authored by Somnath, Suhas's avatar Somnath, Suhas Committed by GitHub
Browse files

Merge pull request #117 from pycroscopy/cades_dev

GIV updates
parents 973516fa b8202f08
......@@ -77,10 +77,6 @@ class Model(object):
verbose : Boolean (Optional)
Whether or not to print log statements
Returns
-------
None
"""
if self._parallel:
......@@ -139,15 +135,12 @@ class Model(object):
def _get_data_chunk(self, verbose=False):
"""
Returns the next chunk of data for the guess or the fit
Reads the next chunk of data for the guess or the fit into memory
Parameters
-----
None
Returns
--------
verbose : bool, optional
Whether or not to print log statements
"""
if self._start_pos < self.h5_main.shape[0]:
self._end_pos = int(min(self.h5_main.shape[0], self._start_pos + self._max_pos_per_read))
......@@ -182,16 +175,19 @@ class Model(object):
else:
self.guess = self.h5_guess[self._start_pos:self._end_pos, :]
def _set_results(self, is_guess=False):
def _set_results(self, is_guess=False, verbose=False):
"""
Writes the provided guess or fit results into appropriate datasets.
Given that the guess and fit datasets are relatively small, we should be able to hold them in memory just fine
Parameters
---------
is_guess : Boolean
is_guess : bool, optional
Default - False
Flag that differentiates the guess from the fit
verbose : bool, optional
Default - False
Whether or not to print log statements
"""
statement = 'guess'
......@@ -203,13 +199,14 @@ class Model(object):
targ_dset = self.h5_fit
source_dset = self.fit
"""print('Writing data to positions: {} to {}'.format(self.__start_pos, self._end_pos))
targ_dset[self._start_pos:self._end_pos, :] = source_dset"""
if verbose:
print('Writing data to positions: {} to {}'.format(self.__start_pos, self._end_pos))
targ_dset[:, :] = source_dset
# flush the file
self.hdf.flush()
print('Finished writing ' + statement + ' results to file!')
if verbose:
print('Finished writing ' + statement + ' results to file!')
def _create_guess_datasets(self):
"""
......@@ -229,9 +226,8 @@ class Model(object):
None
"""
warn('Please override the _create_guess_datasets specific to your model')
self.guess = None # replace with actual h5 dataset
pass
raise NotImplementedError('Please override the _create_guess_datasets specific to your model')
def _create_fit_datasets(self):
"""
......@@ -251,9 +247,8 @@ class Model(object):
None
"""
warn('Please override the _create_fit_datasets specific to your model')
self.fit = None # replace with actual h5 dataset
pass
raise NotImplementedError('Please override the _create_fit_datasets specific to your model')
def do_guess(self, processors=None, strategy='wavelet_peaks',
options={"peak_widths": np.array([10, 200]), "peak_step": 20}):
......
......@@ -18,6 +18,7 @@ from ..io.io_hdf5 import ioHDF5
from ..io.io_utils import recommendCores
from ..io.microdata import MicroDataGroup, MicroDataset
from ..io.hdf_utils import getH5DsetRefs, getAuxData, link_as_main, copyAttributes, linkRefAsAlias
from ..viz.plot_utils import set_tick_font_size
def do_bayesian_inference(V, IV_point, freq, num_x_steps=251, gam=0.03, e=10.0, sigma=10., sigmaC=1.,
......@@ -179,10 +180,10 @@ def do_bayesian_inference(V, IV_point, freq, num_x_steps=251, gam=0.03, e=10.0,
return results_dict
def plot_bayesian_spot_from_h5(h5_bayesian_grp, h5_resh, pix_ind):
def plot_bayesian_spot_from_h5(h5_bayesian_grp, h5_resh, pix_ind, r_extra_override=None):
"""
Plots the basic Bayesian Inference results for a specific pixel
Parameters
----------
h5_bayesian_grp : h5py.Datagroup reference
......@@ -191,6 +192,8 @@ def plot_bayesian_spot_from_h5(h5_bayesian_grp, h5_resh, pix_ind):
Dataset containing the raw / filtered measured current split by pixel
pix_ind : unsigned int
Integer index of the desired pixel
r_extra_override : float, Optional
Default - will not override
Returns
-------
......@@ -202,6 +205,15 @@ def plot_bayesian_spot_from_h5(h5_bayesian_grp, h5_resh, pix_ind):
h5_vr = h5_bayesian_grp['vr']
h5_irec = h5_bayesian_grp['irec']
h5_cap = h5_bayesian_grp['capacitance']
freq = h5_bayesian_grp.attrs['freq']
try:
r_extra = h5_bayesian_grp['Rextra'] # should be 220 Ohms according to calibration
except KeyError:
# Old / incorrect inference model
r_extra = 0
if r_extra_override is not None:
r_extra = 220
possibly_rolled_bias = getAuxData(h5_irec, auxDataName=['Spectroscopic_Values'])[0]
split_directions = h5_bayesian_grp.attrs['split_directions']
......@@ -215,11 +227,11 @@ def plot_bayesian_spot_from_h5(h5_bayesian_grp, h5_resh, pix_ind):
cap_val = h5_cap[pix_ind]
return plot_bayesian_results(orig_bias, possibly_rolled_bias, i_meas, bias_interp, mr_vec, i_recon, vr_vec,
split_directions, cap_val, pix_pos=h5_pos[pix_ind])
split_directions, cap_val, freq, r_extra, pix_pos=h5_pos[pix_ind])
def plot_bayesian_results(orig_bias, possibly_rolled_bias, i_meas, bias_interp, mr_vec, i_recon, vr_vec,
split_directions, cap_val, pix_pos=[0, 0]):
split_directions, cap_val, freq, r_extra, pix_pos=[0, 0]):
"""
Plots the basic Bayesian Inference results for a specific pixel
......@@ -242,7 +254,11 @@ def plot_bayesian_results(orig_bias, possibly_rolled_bias, i_meas, bias_interp,
split_directions : Boolean
Whether or not to compute the forward and reverse portions of the loop separately
cap_val : float
Calculated capacitance
Inferred capacitance in nF
freq : float
Excitation frequency
r_extra : float
Resistance of extra resistor [Ohms] necessary to get correct resistance values
pix_pos : list of two numbers
Pixel row and column positions or values
......@@ -251,19 +267,31 @@ def plot_bayesian_results(orig_bias, possibly_rolled_bias, i_meas, bias_interp,
fig : matplotlib.pyplot figure handle
Handle to figure
"""
font_size_1 = 14
font_size_2 = 16
half_x_ind = int(0.5 * bias_interp.size)
ex_amp = np.max(bias_interp)
colors = [['b', 'r', 'g'], ['violet', 'r', 'g']]
colors = [['red', 'orange'], ['blue', 'cyan']]
syms = [['-', '--', '--'], ['-', ':', ':']]
names = ['Forw', 'Rev']
names = ['Forward', 'Reverse']
if not split_directions:
colors = colors[0]
syms = syms[0]
names = ['']
# Need to calculate the correct resistance in the original domain:
omega = 2 * np.pi * freq
cos_omega_t = np.roll(orig_bias, int(-0.25 * orig_bias.size))
mean_cap_val = 0.5 * np.mean(cap_val) # Not sure why we need the 0.5
i_cap = mean_cap_val * omega * cos_omega_t # * nF -> nA
i_extra = r_extra * mean_cap_val * orig_bias # ohms * nF * V -> nA
i_correct = i_meas - i_cap - i_extra
i_correct_rolled = np.roll(i_correct, int(-0.25 * orig_bias.size))
orig_half_pt = int(0.5 * orig_bias.size)
st_dev = np.sqrt(vr_vec)
good_pts = np.where(st_dev < 10)[0]
good_forw = good_pts[np.where(good_pts < half_x_ind)[0]]
......@@ -271,44 +299,54 @@ def plot_bayesian_results(orig_bias, possibly_rolled_bias, i_meas, bias_interp,
pos_limits = mr_vec + st_dev
neg_limits = mr_vec - st_dev
fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(10, 5))
fig, axes = plt.subplots(ncols=3, figsize=(15, 5))
# fig.subplots_adjust(wspace=3.5)
axes[0].set_xlabel('Voltage (V)')
axes[0].set_ylabel('Resistance (GOhm)')
axes[0].set_title('R(V), mean C = ' + str(round(np.mean(cap_val) * 1E3, 2)) + 'pF, ' + 'at row ' +
str(pix_pos[0]) + ', col ' + str(pix_pos[1]))
axes[0].set_ylabel('Resistance (G\Omega)', fontsize=font_size_2)
if split_directions:
pts_to_plot = [good_forw, good_rev]
else:
pts_to_plot = [good_pts]
for type_ind, pts_list, cols_set, sym_set, set_name in zip(range(len(names)), pts_to_plot, colors, syms, names):
axes[0].plot(bias_interp[pts_list], mr_vec[pts_list], cols_set[0], linestyle=sym_set[0], linewidth=2,
label='R(V) {}'.format(set_name))
axes[0].plot(bias_interp[pts_list], pos_limits[pts_list], cols_set[1], linestyle=sym_set[1],
label='R(V)+$\sigma$ {}'.format(set_name))
axes[0].plot(bias_interp[pts_list], neg_limits[pts_list], cols_set[2], linestyle=sym_set[2],
label='R(V)-$\sigma$ {}'.format(set_name))
for type_ind, axis, pts_list, cols_set, sym_set, set_name in zip(range(len(names)),
axes[:2], pts_to_plot,
colors, syms, names):
axis.set_title('$R(V)$ ' + set_name + ' at Row = ' + str(pix_pos[0]) +
' Col =' + str(pix_pos[1]), fontsize=font_size_2)
axis.plot(bias_interp[pts_list], mr_vec[pts_list], cols_set[0],
linestyle=sym_set[0], linewidth=3, label='R(V)')
axis.fill_between(bias_interp[pts_list], pos_limits[pts_list], neg_limits[pts_list],
alpha=0.25, color=cols_set[1],
label='R(V)+-$\sigma$')
axis.set_xlabel('Voltage (V)', fontsize=font_size_2)
axes[0].legend(loc='best')
axes[0].set_ylim((0, 20))
axes[0].set_xlim((-ex_amp, ex_amp))
axis.legend(loc='upper left', fontsize=font_size_1)
axis.set_xlim((-ex_amp, ex_amp))
axes[1].plot(orig_bias, i_meas, 'b', label='I$_{meas}$')
axes[1].plot(possibly_rolled_bias, i_recon, 'cyan', linestyle='--', label='I$_{rec}$')
# ################### CURRENT PLOT ##########################
axes[2].plot(orig_bias, i_meas, 'b', linewidth=2, label='I$_{meas}$')
axes[2].plot(possibly_rolled_bias, i_recon, 'cyan', linewidth=2, linestyle='--', label='I$_{rec}$')
if split_directions:
axes[1].plot(bias_interp[:half_x_ind], (bias_interp / mr_vec)[:half_x_ind], 'r', label='I$_{Bayes} Forw$')
axes[1].plot(bias_interp[half_x_ind:], (bias_interp / mr_vec)[half_x_ind:], 'orange', label='I$_{Bayes} Rev$')
axes[2].plot(cos_omega_t[:orig_half_pt], i_correct_rolled[:orig_half_pt],
'r', linewidth=3, label='I$_{Bayes} Forw$')
axes[2].plot(cos_omega_t[orig_half_pt:], i_correct_rolled[orig_half_pt:],
'orange', linewidth=3, label='I$_{Bayes} Rev$')
else:
axes[1].plot(bias_interp, bias_interp / mr_vec, 'g', label='I$_{Bayes}$')
axes[1].plot(orig_bias, i_correct, 'g', label='I$_{Bayes}$')
# axes[2].legend(loc='upper right', bbox_to_anchor=(-.1, 0.30), fontsize=font_size_1)
axes[2].legend(loc='best', fontsize=font_size_1)
axes[2].set_xlabel('Voltage(V)', fontsize=font_size_2)
axes[2].set_title('$I(V)$ at row ' + str(pix_pos[0]) + ', col ' + str(pix_pos[1]),
fontsize=font_size_2)
axes[2].set_ylabel('Current (nA)', fontsize=font_size_2)
axes[1].legend(loc='best')
axes[1].set_xlabel('Voltage(V)')
axes[1].set_title('Bayesian Inference at row ' + str(pix_pos[0]) + ', col ' + str(pix_pos[1]))
set_tick_font_size(axes, font_size_1)
axes[1].set_ylabel('Current (nA)')
fig.tight_layout()
return fig
......
......@@ -19,7 +19,7 @@ from ..io.io_hdf5 import ioHDF5
from ..io.hdf_utils import getH5DsetRefs, linkRefs, getAuxData, link_as_main, copyAttributes, copy_main_attributes
from ..io.io_utils import getTimeStamp
from ..io.microdata import MicroDataGroup, MicroDataset
from ..viz.plot_utils import rainbow_plot
from ..viz.plot_utils import rainbow_plot, set_tick_font_size
from ..io.translators.utils import build_ind_val_dsets
......@@ -87,13 +87,19 @@ def test_filter(resp_wfm, filter_parms, samp_rate, show_plots=True, use_rainbow_
Get parameters from the dictionary.
'''
noise_band_filter = filter_parms.get('band_filt_[Hz]', 1)
if noise_band_filter is None:
noise_band_filter = 1
if isinstance(noise_band_filter, Iterable):
noise_band_filter = noiseBandFilter(num_pts, samp_rate, noise_band_filter[0],
noise_band_filter[1])
if verbose and isinstance(noise_band_filter, Iterable):
print('Calculated valid noise_band_filter')
else:
noise_band_filter = 1
low_pass_filter = filter_parms.get('LPF_cutOff_[Hz]', -1)
if low_pass_filter is None:
low_pass_filter = 1
if low_pass_filter > 0:
low_pass_filter = makeLPF(num_pts, samp_rate, low_pass_filter)
if verbose and isinstance(low_pass_filter, Iterable):
......@@ -107,17 +113,22 @@ def test_filter(resp_wfm, filter_parms, samp_rate, show_plots=True, use_rainbow_
harmonic_filter[1], harmonic_filter[2])
if verbose and isinstance(harmonic_filter, Iterable):
print('Calculated valid harmonic filter')
else:
harmonic_filter = 1
composite_filter = noise_band_filter * low_pass_filter * harmonic_filter
noise_floor = filter_parms.get('noise_threshold', None)
if type(noise_floor) not in [float, np.float, np.float16, np.float32, np.float64]:
noise_floor = None
fft_pix_data = np.fft.fftshift(np.fft.fft(resp_wfm))
if 0 < noise_floor < 1:
noise_floor = getNoiseFloor(fft_pix_data, noise_floor)[0]
if show_plots:
l_ind = int(0.5 * num_pts)
if type(composite_filter) == np.ndarray:
if isinstance(composite_filter, np.ndarray):
r_ind = np.max(np.where(composite_filter > 0)[0])
else:
r_ind = num_pts
......@@ -139,17 +150,23 @@ def test_filter(resp_wfm, filter_parms, samp_rate, show_plots=True, use_rainbow_
ax_raw = plt.subplot2grid((2, 4), (0, 0), colspan=lhs_colspan)
ax_filt = plt.subplot2grid((2, 4), (1, 0), colspan=lhs_colspan)
axes = [ax_raw, ax_filt]
set_tick_font_size(axes, 14)
else:
fig = None
axes = None
if show_plots:
amp = np.abs(fft_pix_data)
ax_raw.semilogy(w_vec[l_ind:r_ind], amp[l_ind:r_ind])
ax_raw.semilogy(w_vec[l_ind:r_ind], (composite_filter[l_ind:r_ind] + np.min(amp)) * (np.max(amp) - np.min(amp)))
ax_raw.semilogy(w_vec[l_ind:r_ind], amp[l_ind:r_ind], label='Raw')
ax_raw.semilogy(w_vec[l_ind:r_ind],
(composite_filter[l_ind:r_ind] + np.min(amp)) * (np.max(amp) - np.min(amp)),
linewidth=3, color='orange', label='Composite Filter')
if noise_floor is not None:
ax_raw.semilogy(w_vec[l_ind:r_ind], np.ones(r_ind - l_ind) * noise_floor)
ax_raw.set_title('Raw Signal')
ax_raw.semilogy(w_vec[l_ind:r_ind], np.ones(r_ind - l_ind) * noise_floor,
linewidth=2, color='r', label='Noise Threshold')
ax_raw.legend(loc='best', fontsize=14)
ax_raw.set_title('Raw Signal', fontsize=16)
ax_raw.set_ylabel('Magnitude (a. u.)', fontsize=14)
fft_pix_data *= composite_filter
......@@ -158,8 +175,9 @@ def test_filter(resp_wfm, filter_parms, samp_rate, show_plots=True, use_rainbow_
if show_plots:
ax_filt.semilogy(w_vec[l_ind:r_ind], np.abs(fft_pix_data[l_ind:r_ind]))
ax_filt.set_title('Filtered Signal')
ax_filt.set_xlabel('Frequency(kHz)')
ax_filt.set_title('Filtered Signal', fontsize=16)
ax_filt.set_xlabel('Frequency(kHz)', fontsize=14)
ax_filt.set_ylabel('Magnitude (a. u.)', fontsize=14)
if noise_floor is not None:
ax_filt.set_ylim(bottom=noise_floor) # prevents the noise threshold from messing up plots
......@@ -170,9 +188,10 @@ def test_filter(resp_wfm, filter_parms, samp_rate, show_plots=True, use_rainbow_
rainbow_plot(ax_loops, excit_wfm[l_resp_ind:r_resp_ind], filt_data[l_resp_ind:r_resp_ind] * 1E+3)
else:
ax_loops.plot(excit_wfm[l_resp_ind:r_resp_ind], filt_data[l_resp_ind:r_resp_ind] * 1E+3)
ax_loops.set_title('AI vs AO')
ax_loops.set_xlabel('Input Bias (V)')
ax_loops.set_ylabel('Deflection (mV)')
ax_loops.set_title('AI vs AO', fontsize=16)
ax_loops.set_xlabel('Input Bias (V)', fontsize=14)
ax_loops.set_ylabel('Deflection (mV)', fontsize=14)
set_tick_font_size(ax_loops, 14)
axes.append(ax_loops)
fig.tight_layout()
return filt_data, fig, axes
......
......@@ -4,114 +4,133 @@ Created on 7/17/16 10:08 AM
"""
from __future__ import division, print_function, absolute_import
from warnings import warn
import sys
import numpy as np
import psutil
import scipy
import itertools
import multiprocessing as mp
from ..io.hdf_utils import checkIfMain
from ..io.io_hdf5 import ioHDF5
import multiprocessing as mp
from ..io.io_utils import recommendCores
class Process(object):
"""
Encapsulates the typical steps performed when applying a processing function to a dataset.
Parameters
----------
h5_main : h5py.Dataset instance
The dataset over which the analysis will be performed. This dataset should be linked to the spectroscopic
indices and values, and position indices and values datasets.
"""
def __init__(self, h5_main):
def __init__(self, h5_main, cores=None):
"""
Parameters
----------
h5_main : h5py.Dataset instance
The dataset over which the analysis will be performed. This dataset should be linked to the spectroscopic
indices and values, and position indices and values datasets.
cores : uint, optional
Default - None
How many cores to use for the computation
"""
# Checking if dataset is "Main"
if checkIfMain(h5_main):
self.h5_main = h5_main
self.hdf = ioHDF5(self.h5_main.file)
else:
warn('Provided dataset is not a "Main" dataset with necessary ancillary datasets')
return
raise ValueError('Provided dataset is not a "Main" dataset with necessary ancillary datasets')
# Determining the max size of the data that can be put into memory
self._setMemoryAndCPUs()
self._set_memory_and_cores(cores=cores)
self._start_pos = 0
self._end_pos = self.h5_main.shape[0]
def _setMemoryAndCPUs(self):
self._results = None
self.h5_results_grp = None
def _set_memory_and_cores(self, cores=None, verbose=False):
"""
Checks hardware limitations such as memory, # cpus and sets the recommended datachunk sizes and the
number of cores to be used by analysis methods.
Returns
-------
None
Parameters
----------
cores : uint, optional
Default - None
How many cores to use for the computation
verbose : bool, optional
Whether or not to print log statements
"""
if self._parallel:
self._maxCpus = psutil.cpu_count() - 2
if cores is None:
self._cores = psutil.cpu_count() - 2
else:
self._maxCpus = 1
self._maxMemoryMB = psutil.virtual_memory().available / 1e6 # in MB
cores = int(abs(cores))
self._cores = min(psutil.cpu_count(), max(1, cores))
self._max_mem_mb = psutil.virtual_memory().available / 1e6 # in MB
self._maxDataChunk = self._maxMemoryMB / self._maxCpus
max_data_chunk = self._max_mem_mb / self._cores
# Now calculate the number of positions that can be stored in memory in one go.
mb_per_position = self.h5_main.dtype.itemsize * self.h5_main.shape[1] / 1e6
self._max_pos_per_read = int(np.floor(self._maxDataChunk / mb_per_position))
print('Allowed to read {} pixels per chunk'.format(self._max_pos_per_read))
self._max_pos_per_read = int(np.floor(max_data_chunk / mb_per_position))
def _get_data_chunk(self):
"""
Returns a chunk of data for the guess or the fit
if verbose:
print('Allowed to read {} pixels per chunk'.format(self._max_pos_per_read))
print('Allowed to use up to', str(self._cores), 'cores and', str(self._max_mem_mb), 'MB of memory')
Parameters
----------
None
def _unit_function(*args):
raise NotImplementedError('Please override the _create_results_datasets specific to your model')
Returns
-------
def _read_data_chunk(self, verbose=False):
"""
Reads a chunk of data for the intended computation into memory
Parameters
-----
verbose : bool, optional
Whether or not to print log statements
"""
if self._start_pos < self.h5_main.shape[0]:
self._end_pos = int(min(self.h5_main.shape[0], self._start_pos + self._max_pos_per_read))
self.data = self.h5_main[self._start_pos:self._end_pos, :]
print('Reading pixels {} to {} of {}'.format(self._start_pos, self._end_pos, self.h5_main.shape[0]))
if verbose:
print('Reading pixels {} to {} of {}'.format(self._start_pos, self._end_pos, self.h5_main.shape[0]))
# DON'T update the start position
else:
print('Finished reading all data!')
if verbose:
print('Finished reading all data!')
self.data = None
def _set_results(self):
def _write_results_chunk(self):
"""
Writes the provided guess or fit results into appropriate datasets.
Given that the guess and fit datasets are relatively small, we should be able to hold them in memory just fine
Writes the computed results into appropriate datasets.
This needs to be rewritten since the processed data is expected to be at least as large as the dataset
"""
# Now update the start position
self._start_pos = self._end_pos
warn('Please override the _set_results specific to your model')
pass
raise NotImplementedError('Please override the _set_results specific to your model')
def _create_results_datasets(self):
"""
Model specific call that will write the h5 group, guess dataset, corresponding spectroscopic datasets and also
Process specific call that will write the h5 group, guess dataset, corresponding spectroscopic datasets and also
link the guess dataset to the spectroscopic datasets. It is recommended that the ancillary datasets be populated
within this function.
"""
warn('Please override the _create_results_datasets specific to your model')
pass
raise NotImplementedError('Please override the _create_results_datasets specific to your model')
def compute(self, func, **kwargs):
def compute(self, *args, **kwargs):
"""
Parameters
----------
func
kwargs:
processors: int
number of processors to use. Default all processors on the system except for 1.
......@@ -124,40 +143,78 @@ class Process(object):
self._create_results_datasets()
self._start_pos = 0
processors = kwargs.get("processors", self._maxCpus)
results = list()
if self._parallel:
# start pool of workers
print('Computing in parallel ... launching %i kernels...' % processors)
pool = mp.Pool(processors)
self._get_data_chunk()
while self.data is not None: # as long as we have not reached the end of this data set:
# apply guess to this data chunk:
tasks = [vector for vector in self.data]
chunk = int(self.data.shape[0] / processors)
jobs = pool.imap(func, tasks, chunksize=chunk)
# get Results from different processes
print('Extracting Results...')
temp = [j for j in jobs]
# Reformat the data to the appropriate type and or do additional computation now
# Write to file
self._set_results()
# read the next chunk
self._get_data_chunk()
# Finished reading the entire data set
print('closing %i kernels...' % processors)
pool.close()
else:
print("Computing Guesses In Serial ...")
self._get_data_chunk()
while self.data is not None: # as long as we have not reached the end of this data set:
temp = [func(vector) for vector in self.data]
results.append(self._reformatResults(temp))
# read the next chunk
self._get_data_chunk()
# reorder to get one numpy array out