Ultrafast current imaging by Bayesian inversion¶
Nature Communications, 6th November 2017
S. Somnath1,2,3, K. J. H. Law1,4, A. N. Morozovska5, P. Maksymovych1,2, Y. Kim6, X. Lu7, M. Alexe8, R. Archibald1,4, S. V. Kalinin1,2, S. Jesse1,2 and R. K. Vasudevan1,2
- Institute for Functional Imaging of Materials,
- Center for Nanophase Materials Sciences,
- Advanced Data and Workflows Group, and
Computer Science and Mathematics Division, Oak Ridge National Laboratory, Oak Ridge TN 37831, USA
Institute of Physics, National Academy of Sciences of Ukraine, 46, pr. Nauky, 03028 Kyiv, Ukraine
- School of Advanced Materials Science and Engineering, Sungkyunkwan University (SKKU), Suwon 16419, Republic of Korea
- The State Key Discipline Laboratory of Wide Band Gap Semiconductor Technology, Xidian University, Xi’an, Shaanxi 710071, China
- Department of Physics, University of Warwick, Coventry CV4 7AL, United Kingdom
Installing required packages:¶
In [1]:
!pip install h5py matplotlib numpy scipy sklearn pycroscopy
# Checks Python Version
import sys
if sys.version_info < (3, 5):
print('''This notebook was optimized to work on Python 3.5.
While it may also run on other Python versions,
functionality and performance are not guaranteed
Please consider upgrading your python version.''')
Importing libraries and setting up the notebook¶
In [1]:
# Ensure that this code works on both python 2 and python 3
from __future__ import division, print_function, absolute_import, unicode_literals
# basic numeric computation:
import numpy as np
# Image manipulation
from scipy.ndimage import zoom
# The package used for creating and manipulating HDF5 files:
import h5py
# Plotting and visualization:
import matplotlib.pyplot as plt
# set up notebook to show plots within the notebook
% matplotlib inline
# Data mining
import scipy
from sklearn.cluster import KMeans, DBSCAN
# Miscellaneous pacakges
import os
import sys
import time as tm
# For the main processing, analysis, I/O and visualization functions
import pycroscopy as px
Defining a few handy plot functions:¶
In [24]:
def twin_img_plots(img_list, titles, **kwargs):
if len(img_list) != len(titles):
return None
fig, axes = plt.subplots(ncols=len(img_list), figsize=(11.5, 5))
for axis, img, title in zip(axes.flat, img_list, titles):
im_handle, cbar_handle = px.plot_utils.single_img_cbar_plot(axis, img, **kwargs)
axis.set_title(title, fontsize=16)
axis.set_xlabel('X ($\mu$m)', fontsize=16)
axes[0].set_ylabel('Y ($\mu$m)', fontsize=16)
fig.tight_layout()
return fig, axes
def plot_cluster_results_separately(labels_mat, cluster_centroids, bias_vec, legend_mode=1, **kwargs):
num_clusters = cluster_centroids.shape[0]
fig_lab, axis_lab = plt.subplots(figsize=(5.5,5))
_, _ = px.plot_utils.single_img_cbar_plot(axis_lab, labels_mat,
clim=[0, num_clusters-1],
cmap=plt.get_cmap('viridis', num_clusters),
aspect='auto', show_xy_ticks=True, **kwargs)
axis_lab.set_xlabel('X ($\mu$m)', fontsize=16)
axis_lab.set_ylabel('Y ($\mu$m)', fontsize=16)
axis_lab.set_title('K-Means Cluster Labels', fontsize=16)
fig_lab.tight_layout()
# Plot centroids
fig_width = 5.0
if legend_mode not in [0, 1]:
fig_width = 5.85
fig_centroids, axis_centroids = plt.subplots(figsize=(fig_width, 5))
colors = [ plt.cm.viridis(x) for x in np.linspace(0, 1, cluster_centroids.shape[0]) ]
# print('Number of pixels in each cluster:')
for line_ind in range(cluster_centroids.shape[0]):
cmap=plt.cm.jet
line_color=colors[line_ind]
line_label = 'Cluster ' + str(line_ind)
num_of_cluster_members = len(np.where(labels==line_ind)[0])
# print ("Cluster " + str(line_ind) + ': ' + str(num_of_cluster_members))
if num_of_cluster_members > 10:
axis_centroids.plot(bias_vec, cluster_centroids[line_ind,:],
label=line_label, color=line_color) # marker='o',
axis_centroids.set_xlabel('Voltage (V)', fontsize=16)
axis_centroids.set_ylabel('Current (nA)', fontsize=16)
axis_centroids.set_title('K-Means Cluster Centroids', fontsize=16)
if legend_mode==0:
axis_centroids.legend(loc='lower right', fontsize=14)
elif legend_mode==1:
axis_centroids.legend(loc='upper left', bbox_to_anchor=(1, 1), fontsize=14)
else:
sm = px.plot_utils.make_scalar_mappable(0, num_clusters-1,
cmap=px.plot_utils.discrete_cmap(num_clusters))
plt.colorbar(sm)
px.plot_utils.set_tick_font_size(axis_centroids, 14)
fig_centroids.tight_layout()
return fig_lab, fig_centroids
Selecting the gIV Dataset¶
Also setting up some basic parameters
In [3]:
# Load the raw dataset
h5_path = 'pzt_nanocap_6_bayesian_complete.h5'#'pzt_nanocap_6_ready_for_bayesian-Copy1.h5'
#h5_path = px.io_utils.uiGetFile('*.h5', 'G-mode IV dataset')
print('Working on:\n',h5_path)
# Store some basic parameters in memory
folder_path, filename = os.path.split(h5_path)
try:
os.mkdir(os.path.join(folder_path, 'GIV_Paper_Figures'))
except FileExistsError:
pass
try:
os.mkdir(os.path.join(folder_path, 'GIV_Other_Figures'))
except FileExistsError:
pass
figure_folder = os.path.join(folder_path, 'GIV_Paper_Figures')
other_figures_folder = os.path.join(folder_path, 'GIV_Other_Figures')
Loading the Data¶
Also loading ancillary datasets and necessary parameters
In [4]:
h5_f = h5py.File(h5_path, 'r')
h5_grp = h5_f['Measurement_000/Channel_000']
# This is the reference to the raw data as it was acquired from the microscope
h5_main = h5_grp['Raw_Data']
print('Raw data is of shape:', h5_main.shape)
num_lines = h5_main.shape[0]
num_pts = h5_main.shape[1]
# Reading experimental parameters:
scan_width = h5_grp.attrs['grid_scan_width_[m]']*1E+6
samp_rate = h5_grp.attrs['IO_samp_rate_[Hz]']
ex_freq = h5_grp.attrs['excitation_frequency_[Hz]']
# Getting ancillary information and other parameters
h5_spec_vals = px.hdf_utils.getAuxData(h5_main, auxDataName=['Spectroscopic_Values'])[0]
# Excitation waveform for a single line / row of data
excit_wfm = h5_spec_vals.value
# We expect each pixel to have a single period of the sinusoidal excitation
# Calculating the excitation waveform for a single pixel
pts_per_cycle = int(np.round(samp_rate/ex_freq))
single_AO = excit_wfm[0, :pts_per_cycle]
# Number of pixels within a line:
num_cols = int(h5_main.shape[1] / single_AO.size)
print('Number of pixels:', num_cols, ', each with points:', single_AO.shape)
ex_amp = h5_grp.attrs['excitation_amplitude_[V]']
Visualizing raw data:¶
Here we will inspect what the IV curves look like for a single scan line.
Note the significant hysteresis in the IV curves.
In [4]:
row_ind = 20
# read data for a specific scan line
raw_line_resp = h5_main[row_ind]
# break this up into pixels:
raw_line_mat = np.reshape(raw_line_resp, (-1, single_AO.size))
# Now visualize:
fig, axes = px.plot_utils.plot_loops(single_AO, raw_line_mat, use_rainbow_plots=False, x_label='Bias (V)',
y_label='Current (nA)', subtitles='Pixel', title=None)
fig.savefig(os.path.join(other_figures_folder,
'example_raw_curves_line_' + str(row_ind) +'.png'),
format='png', dpi=150);
Visualizing information in Fourier space¶
Visualizing in the fourier space provides information about the noise floor, frequencies which are noise dominant or signal dominant, etc.
This visualization will guide the design of signal filters to remove the noise
In [5]:
# Preparing the frequency axis:
w_vec = 1E-3*np.linspace(-0.5*samp_rate, 0.5*samp_rate - samp_rate/num_pts, num_pts)
row_ind = 105
F_resp = np.fft.fftshift(np.fft.fft(h5_main[row_ind]))
fig, ax = plt.subplots(figsize=(12, 7))
ax.axvline(x=1E-3*ex_freq, color='r', linewidth=2, label='Excitation')
ax.plot(w_vec[int(0.5*len(w_vec)):], np.log10(np.abs(F_resp[int(0.5*len(w_vec)):])), label='Response')
ax.set_xlabel('Frequency (kHz)', fontsize=16)
ax.set_ylabel('Amplitude (a.u.)', fontsize=16)
ax.legend(fontsize=14)
ax.set_xscale('log')
ax.set_xlim(1E-2, samp_rate*0.5E-3)
ax.set_title('Noise Spectrum for row ' + str(row_ind), fontsize=16)
px.plot_utils.set_tick_font_size(ax, 14)
fig.savefig(os.path.join(other_figures_folder,
'noise_spectrum_line_' + str(row_ind) +'.png'),
format='png', dpi=150);
In [6]:
# Set Filter parameters here:
filter_parms = dict()
filter_parms['noise_threshold'] = 1E-6
filter_parms['comb_[Hz]'] = -1
filter_parms['LPF_cutOff_[Hz]'] = 10E+3
# Noise frequencies - 15.6 kHz ~ 14-17.5, 7.8-8.8, 45-49.9 ~ 48.9414 kHz
filter_parms['band_filt_[Hz]'] = -1 # [[8.3E+3, 15.6E+3, 48.9414E+3], [1E+3, 0.5E+3, 0.1E+3]]
filter_parms['phase_[rad]'] = 0
filter_parms['samp_rate_[Hz]'] = samp_rate
filter_parms['num_pix'] = 1
Visualizing Filtering Results¶
Run the next cell to see what the IV curves look like for the chosen filter parameters. If the current set of filter parameters do not work as well, change the parameters in the previous cell and try again.
Once the results look good, proceed to filter the entire dataset
In [7]:
# Test filter on a single line:
row_ind = 50
filt_line, fig_filt, axes_filt = px.processing.gmode_utils.test_filter(h5_main[row_ind], filter_parms, samp_rate,
show_plots=True, use_rainbow_plots=False)
fig_filt.savefig(os.path.join(other_figures_folder, 'FFT_filter_on_line_{}.png'.format(row_ind)), format='png', dpi=300)
raw_row = np.reshape(h5_main[row_ind], (-1, pts_per_cycle))
filt_row = filt_line.reshape(-1, pts_per_cycle)
fig, axes = px.plot_utils.plot_loops(single_AO, [raw_row, filt_row], dataset_names=['Raw', 'Filtered'],
line_colors=['r', 'b'], x_label='Bias (V)', title='FFT Filtering',
plots_on_side=3, y_label='Current (nA)',
subtitles='Row: ' + str(row_ind) + ' Col:')
fig.savefig(os.path.join(other_figures_folder, 'Example_filtered_loops_from_line_{}.png'.format(row_ind)), format='png', dpi=300)
Filter the entire dataset:¶
Once the optimal filter parameters are chosen, the same parameters can be applied to filter the entire dataset. This can take a few minutes.
In [5]:
h5_filt_grp = px.hdf_utils.findH5group(h5_main, 'FFT_Filtering')
if len(h5_filt_grp) > 0:
print('Taking previously filtered results')
h5_filt_grp = h5_filt_grp[-1]
else:
print('FFT filtering not performed on this dataset. Filtering now:')
h5_filt_grp = px.processing.gmode_utils.fft_filter_dataset(h5_main, filter_parms, write_filtered=True)
h5_filt = h5_filt_grp['Filtered_Data']
Restructuring the data¶
G-mode IV data is acquried continuously for a line of scan, the data therefore is of the form line x time and needs to be restructured to row x col x time where each pixel contains a single period of the sinusoidal excitation waveform
In [6]:
# Check to see if this restructuring was already performed on the dataset:
h5_resh_grp = px.hdf_utils.findH5group(h5_filt, 'Reshape')
if len(h5_resh_grp) > 0:
print('Taking previously reshaped results')
h5_resh_grp = h5_resh_grp[-1]
else:
print('Reshape not performed on this dataset. Reshaping now:')
scan_width = h5_grp.attrs['grid_scan_width_[m]']
h5_resh = px.processing.gmode_utils.reshape_from_lines_to_pixels(h5_filt, pts_per_cycle, scan_width / num_cols)
h5_resh = h5_resh_grp['Reshaped_Data']
print('Data was reshaped from shape', h5_filt.shape,
'reshaped to ', h5_resh.shape)
Bayesian Inference¶
Now that the data has been filtered and reshaped to rows, columns, and time, we can finally perform Bayesian inference on the reshaped dataset.
Try out Bayesian Inference in the cell below:
In [6]:
col = 166
row = 116
pix_ind = row * num_cols + col
i_meas = h5_resh[pix_ind]
results = px.processing.giv_utils.bayesian_inference_on_period(i_meas, single_AO, ex_freq, r_extra=220, num_x_steps=250,
show_plots=True)
Bayesian Inference on dataset¶
Since Bayesian inference is computationally expensive, we have an option to load existing results in the data file
In [7]:
# Load bayesian inference
h5_bayesian_grp = px.hdf_utils.findH5group(h5_resh, 'Bayesian_Inference')
amp_gain = 9
if len(h5_bayesian_grp) == 0:
print('No Bayesian Inference results found - computing now')
t_start = tm.time()
i_cleaner = px.GIVBayesian(h5_resh, ex_freq, amp_gain, r_extra=220, num_x_steps=250, cores=None,
max_mem_mb=1024*10, verbose=True)
h5_bayesian_grp = i_cleaner.compute()
"""h5_bayesian_grp = px.processing.giv_utils.bayesian_inference_dataset(h5_resh, ex_freq, amp_gain,
split_directions=True, Rextra = 52.0,
num_x_steps=125, verbose=False)"""
t_end = tm.time()
print('Took a total of {} hours to apply Bayesian Inference on the entire dataset'.format(np.round((t_end-t_start))))
else:
print('Taking previous results already present in file')
h5_bayesian_grp = h5_bayesian_grp[-1]
bias_interp = np.squeeze(h5_bayesian_grp['Spectroscopic_Values'][()])
h5_resistance = h5_bayesian_grp['Resistance']
h5_capacitance = h5_bayesian_grp['Capacitance']
h5_i_corr_sine = h5_bayesian_grp['Corrected_Current']
Visuaizing Inference Results¶
We will look at the inference results from five randomly chosen pixels
In [7]:
# Visualize Bayesian Inference results
num_rand_pos = 5
chosen_rows = np.random.randint(0, num_lines, num_rand_pos)
chosen_cols = np.random.randint(0, num_cols, num_rand_pos)
for r_mode in [False]:
if r_mode:
suffix = '_scatter'
else:
suffix = '_line'
for index in range(num_rand_pos):
pix_ind = chosen_rows[index] * num_cols + chosen_cols[index]
fig = px.processing.giv_utils.plot_bayesian_spot_from_h5(h5_bayesian_grp, h5_resh, pix_ind, r_max=300,
res_scatter=r_mode)
"""file_name = 'Bayesian_Infer_row_' + str(chosen_rows[index]) + '_col_' + str(chosen_cols[index]) + suffix
fig.savefig(os.path.join(other_figures_folder, file_name + '.png'), format='png', dpi=300)
fig.savefig(os.path.join(figure_folder, file_name + '.pdf'), format='pdf', bbox_inches = 'tight',
pad_inches = 2.0)"""
Visualize current maps at different biases¶
In [24]:
cur_map = currents[:, :, 4]
avg_val = np.mean(cur_map)
std_val = np.std(cur_map)
print(avg_val, std_val)
# im_hand, cbar_hand = px.plot_utils.single_img_cbar_plot(axis, cur_map, show_cbar=True)
fig, axis = plt.subplots()
im_handle = px.plot_utils.plot_map(axis, cur_map, stdevs=0.01)
#im_handle = axis.imshow(cur_map)
cbar = plt.colorbar(im_handle, ax=axis, orientation='vertical',
fraction=0.046, pad=0.04, use_gridspec=True)
In [11]:
"""num_plots = 16
chosen_bias_inds = np.linspace(0, bias_interp.size-1, num_plots, endpoint=True, dtype=np.uint16)
currents = np.reshape(h5_resistance[()], (num_lines, num_cols, -1))
voltages = bias_interp[chosen_bias_inds]
currents = voltages / currents[:, :, chosen_bias_inds]
axis_subtitles = ['V = %.2f V'% cur_val for cur_val in voltages]"""
fig, axes = px.plot_utils.plot_map_stack(currents, num_comps=16, color_bar_mode="each", evenly_spaced=True,
title=axis_subtitles, heading='Current at different voltages', reverse_dims=True,
pad_mult=(0.15, 0.07))
Main trends in the Current¶
We understand the main trends in the current data via K-Means clustering
In [26]:
num_clusters = 7
estimators = KMeans(num_clusters)
rolled_bias = np.roll(single_AO, -1*h5_i_corr_sine.shape[1] // 4)
current_cos = np.roll(h5_i_corr_sine[()], -1*h5_i_corr_sine.shape[1] //4, axis=1)
Forward First:¶
In [28]:
results = estimators.fit(current_cos[:, h5_i_corr_sine.shape[1] //2:])
labels_forw, centroids_forw = px.processing.cluster.reorder_clusters(results.labels_,
results.cluster_centers_)
fig_lab_forw, fig_centr_forw = plot_cluster_results_separately(labels_forw.reshape(num_lines, num_cols),
centroids_forw,
rolled_bias[h5_i_corr_sine.shape[1] //2:],
legend_mode=2,
x_size=scan_width, y_size=scan_width)
file_name = 'k_means_I_corr_forward'
for suffix, fig in zip(['_labels', '_centroids'], [fig_lab_forw, fig_centr_forw]):
new_fig_name = fig_name + suffix
fig.savefig(os.path.join(figure_folder, new_fig_name + '.png'), format='png', dpi=300)
fig.savefig(os.path.join(figure_folder, new_fig_name + '.pdf'), format='pdf', bbox_inches = 'tight',
pad_inches = 2.0)
Reverse next:¶
In [29]:
results = estimators.fit(current_cos[:, :h5_i_corr_sine.shape[1] //2])
labels_rev, centroids_rev = px.processing.cluster.reorder_clusters(results.labels_,
results.cluster_centers_)
fig_lab_rev, fig_centr_rev = plot_cluster_results_separately(labels_rev.reshape(num_lines, num_cols),
centroids_rev,
rolled_bias[:h5_i_corr_sine.shape[1] //2],
legend_mode=2,
x_size=scan_width, y_size=scan_width)
file_name = 'k_means_I_corr_reverse'
for suffix, fig in zip(['_labels', '_centroids'], [fig_lab_rev, fig_centr_rev]):
new_fig_name = fig_name + suffix
fig.savefig(os.path.join(figure_folder, new_fig_name + '.png'), format='png', dpi=300)
fig.savefig(os.path.join(figure_folder, new_fig_name + '.pdf'), format='pdf', bbox_inches = 'tight',
pad_inches = 2.0)
Visualize resistance at different biases¶
In [29]:
resistances = np.reshape(h5_resistance[()], (num_lines, num_cols, -1))
num_plots = 16
voltages = bias_interp[np.linspace(0, bias_interp.size-1, num_plots, endpoint=True, dtype=np.uint16)]
axis_subtitles = ['V = %.2f V'% cur_val for cur_val in voltages]
fig, axes = px.plot_utils.plot_map_stack(resistances, num_comps=16, color_bar_mode="each", evenly_spaced=True,
title=axis_subtitles, heading='Resistance at different voltages', stdevs=0.2)
"""
fig_name = "resistance_slices"
fig.savefig(os.path.join(figure_folder, fig_name + '.png'), format='png', dpi=300)
fig.savefig(os.path.join(figure_folder, fig_name + '.pdf'), format='pdf', dpi=300)"""
Out[29]:
Visualize Capacitance¶
In [33]:
cap_vec = px.io_utils.compound_to_scalar(h5_capacitance[()])
cap_vec = cap_vec.mean(axis=1)* 1E+3 # Convert to pF
fig, axis = plt.subplots(figsize=(6, 5))
single_img_cbar_plot(fig, axis, np.reshape(cap_vec, (num_lines, num_cols)),
x_size=scan_width, y_size=scan_width, stdevs=0.25)
axis.set_title('Capacitance (pF)', fontsize=16)
axis.set_xlabel('X ($\mu$m)', fontsize=16)
axis.set_ylabel('Y ($\mu$m)', fontsize=16)
fig.tight_layout()
"""fig.savefig(os.path.join(figure_folder, 'capacitance_map.png'), format='png', dpi=300)
fig.savefig(os.path.join(figure_folder, 'capacitance_map.pdf'), format='pdf',bbox_inches = 'tight', pad_inches = 2.0)"""
Out[33]:
Extracting Switching Current¶
- The forward switching current is calculated as the difference between the forward and reverse currents for positive biases.
- The reverse switching current is calculated as the difference between the reverse and the forward currents for negative biases
In [17]:
col = 166
row = 116
pix_ind = row * num_cols + col
mr_vec = h5_resistance[pix_ind]
half_ind = mr_vec.size // 2
yvec_fwd = (bias_interp / mr_vec)[:half_ind]
yvec_rev = (bias_interp / mr_vec)[half_ind:]
xvec_fwd = bias_interp[:half_ind]
xvec_rev = bias_interp[half_ind:]
#Set these limits manually
pos_coer_ind = np.argmin(np.abs(xvec_fwd - 2))
neg_coer_ind = np.argmin(np.abs(xvec_rev - -2))
xvec_Pr_fwd = xvec_fwd[pos_coer_ind:]
yvec_fwd_switched = yvec_fwd[pos_coer_ind:]
# remember that reverse is flipped!
yvec_fwd_unswitched = np.flipud(yvec_rev[:(half_ind-pos_coer_ind)])
Pr_forward = yvec_fwd_switched-yvec_fwd_unswitched
xvec_Pr_rev = xvec_rev[neg_coer_ind:]
yvec_rev_switched = yvec_rev[neg_coer_ind:]
# remember that reverse is flipped!
yvec_rev_unswitched = np.flipud(yvec_fwd[:(half_ind-neg_coer_ind)])
Pr_reverse = yvec_rev_switched-yvec_rev_unswitched
integral_fwd_trace = np.trapz(Pr_forward, xvec_Pr_fwd)
integral_rev_trace = np.trapz(Pr_reverse, xvec_Pr_rev)
print ('Areas of fwd switching is: ' + str(abs(round(integral_fwd_trace,4))),
' and rev is: ', str(round(abs(integral_rev_trace),4)))
Plot the switching current for this pixel¶
In [19]:
font_size_1 = 14
font_size_2 = 16
fig201,ax201 = plt.subplots(ncols=3, figsize=(15, 4.5))
ax201[0].plot(xvec_fwd,yvec_fwd , 'r-', linewidth=2, label='Forward')
ax201[0].plot(xvec_rev,yvec_rev , 'b-', linewidth=2, label='Reverse')
ax201[0].plot(xvec_Pr_fwd, yvec_fwd_switched,
'g--', label='Forward NS')
ax201[0].plot(xvec_Pr_rev, yvec_rev_switched,
'k--', label='Reverse NS')
ax201[0].legend(loc='best', fontsize=font_size_1)
ax201[0].set_title('Row = ' + str(col) + ' Col = ' + str(row), fontsize=font_size_2)
ax201[0].set_ylabel('Current (nA)', fontsize=font_size_2)
ax201[0].set_xlabel('Voltage (V)', fontsize=font_size_2)
ax201[0].set_ylim(-0.8, 0.8)
ax201[1].plot(xvec_Pr_fwd,yvec_fwd_switched,
'b-', linewidth=2, label='Forward')
ax201[1].plot(xvec_Pr_fwd,yvec_fwd_unswitched,
'g-', linewidth=2, label='Reverse NS')
ax201[1].plot(xvec_Pr_fwd,yvec_fwd_switched-yvec_fwd_unswitched,
'k-', linewidth=2, label='Difference')
ax201[1].legend(loc='upper right', fontsize=font_size_1)
ax201[1].set_xlabel('Voltage (V)', fontsize=font_size_2)
ax201[1].set_title('Switching on Forward', fontsize=font_size_2)
ax201[2].plot(xvec_Pr_rev,yvec_rev_switched,
'b-', linewidth=2, label='Reverse')
ax201[2].plot(xvec_Pr_rev,yvec_rev_unswitched,
'g-', linewidth=2, label='Forward NS')
ax201[2].plot(xvec_Pr_rev,yvec_rev_switched-yvec_rev_unswitched,
'k-', linewidth=2, label='Difference')
ax201[2].set_ylim(ax201[2].get_ylim()[0], 0.4)
ax201[2].legend(loc='upper right', fontsize=font_size_1)
ax201[2].set_xlabel('Voltage (V)', fontsize=font_size_2)
ax201[2].set_title('Switching on Reverse', fontsize=font_size_2)
px.plot_utils.set_tick_font_size(ax201, font_size_1)
"""
file_name = 'Switched_area_single_point'
fig201.savefig(os.path.join(figure_folder, file_name + '.png'), format='png', dpi=300)
fig201.savefig(os.path.join(figure_folder, file_name + '.pdf'), format='pdf', bbox_inches = 'tight',
pad_inches = 2.0)"""
Out[19]:
Calculate Polarization¶
Polarization is roughly equal to the area under the 'ears' arising from the switching current
In [9]:
def get_polarization(mr_vec, bias_interp, ex_freq, rev_bias_thresh, forw_bias_thresh, thresh_wind):
half_ind = mr_vec.size // 2
yvec_fwd = (bias_interp / mr_vec)[:half_ind]
yvec_rev = (bias_interp / mr_vec)[half_ind:]
xvec_fwd = bias_interp[:half_ind]
xvec_rev = bias_interp[half_ind:]
#Set these limits manually
pos_coer_ind = np.argmin(np.abs(xvec_fwd - forw_bias_thresh))
neg_coer_ind = np.argmin(np.abs(xvec_rev - rev_bias_thresh))
xvec_Pr_fwd = xvec_fwd[pos_coer_ind:]
yvec_fwd_switched = yvec_fwd[pos_coer_ind:]
# remember that reverse is flipped!
yvec_fwd_unswitched = np.flipud(yvec_rev[:(half_ind-pos_coer_ind)])
Pr_forward = yvec_fwd_switched-yvec_fwd_unswitched
xvec_Pr_rev = xvec_rev[neg_coer_ind:]
yvec_rev_switched = yvec_rev[neg_coer_ind:]
# remember that reverse is flipped!
yvec_rev_unswitched = np.flipud(yvec_fwd[:(half_ind-neg_coer_ind)])
Pr_reverse = yvec_rev_switched-yvec_rev_unswitched
integral_fwd_trace = np.trapz(Pr_forward, xvec_Pr_fwd)
integral_rev_trace = np.trapz(Pr_reverse, xvec_Pr_rev)
#Go from the switched area to the polarization with a simple calculation
time_ratio_fwd = float(len(xvec_Pr_fwd))/float(len(xvec_fwd))
time_ratio_rev = float(len(xvec_Pr_rev))/float(len(xvec_rev))
#Calculate the integral for this time step
tstep = (1/ex_freq)*0.5 #time per voltage step (half cycle)
t_fwd = np.linspace(0,tstep*time_ratio_fwd, len(xvec_Pr_fwd))
t_rev = np.linspace(0,tstep*time_ratio_rev, len(xvec_Pr_rev))
Q_forward = np.trapz(Pr_forward*1E-9, t_fwd)
Q_reverse = np.trapz(Pr_reverse*1E-9, t_rev)
ind_fwd = np.where(Pr_forward == np.max(Pr_forward))
ind_rev = np.where(Pr_reverse == np.min(Pr_reverse))
v_fwd = xvec_Pr_fwd[ind_fwd]
v_rev = xvec_Pr_rev[ind_rev]
return {'integrals': [integral_fwd_trace, integral_rev_trace],
'charge': np.abs(np.hstack((Q_forward, Q_reverse))),
'biases': np.hstack((v_fwd,v_rev)),
'indices': np.hstack((ind_fwd,ind_rev)),
'Polarization_Forward': Pr_forward,
'Polarization_Reverse': Pr_reverse}
In [10]:
tstep = (1/ex_freq)*0.5 #time per voltage step (half cycle)
# Manual input here:
pos_bias = 1.8
neg_bias = -1 * pos_bias
h5_resistance = h5_bayesian_grp['Resistance']
bias_interp = np.squeeze(h5_bayesian_grp['Spectroscopic_Values'][()])
half_ind = bias_interp.size // 2
integral_area = np.zeros(shape=(num_lines,num_cols,2)) #For forward and reverse, integrate with respect to V
Q_integrated = np.zeros(shape=(num_lines,num_cols,2)) #For forward and reverse, integrate with respect to t
pos_of_peak = np.zeros(shape=(num_lines,num_cols,2)) #For forward and reverse, voltage where maximum I occurs
peak_index = np.zeros(shape=(num_lines,num_cols,2)) #For forward and reverse, index where maximum I occurs
# We need to determine the size of the vector which will have switching. We can set the limits manually
xvec_fwd = bias_interp[:half_ind]
ps_fwd_size = half_ind - np.argmin(np.abs(xvec_fwd - pos_bias))
xvec_Pr_fwd = xvec_fwd[np.argmin(np.abs(xvec_fwd - pos_bias)):]
time_ratio_fwd = float(len(xvec_Pr_fwd))/float(len(xvec_fwd))
# Same for reverse
xvec_rev = bias_interp[half_ind:]
ps_rev_size = half_ind - np.argmin(np.abs(xvec_rev - neg_bias))
#Now initialize the matrices
Pswitching_fwd = np.zeros(shape=(num_lines,num_cols, ps_fwd_size)) #For forward and reverse,
#index where maximum I occurs
Pswitching_rev = np.zeros(shape=(num_lines,num_cols, ps_rev_size)) #For forward and reverse,
#index where maximum I occurs
import time
t0 = time.time()
results = px.processing.parallel_compute(h5_resistance[()], get_polarization, cores=None,
func_args=[bias_interp, ex_freq, neg_bias, pos_bias, 0.5])
for row_ind in range(num_lines):
for col_ind in range(num_cols):
ret_dict = results[row_ind*num_cols + col_ind ]
integral_area[row_ind,col_ind] = ret_dict['integrals']
Q_integrated[row_ind, col_ind,:] = ret_dict['charge']
pos_of_peak[row_ind,col_ind,:] = ret_dict['biases']
peak_index[row_ind,col_ind,:] = ret_dict['indices']
Pswitching_fwd[row_ind,col_ind,:] = ret_dict['Polarization_Forward']
Pswitching_rev[row_ind,col_ind,:] = ret_dict['Polarization_Reverse']
print('Took', time.time()-t0,'seconds')
Visualize Switching Current¶
Via the area under the switching current curves
In [12]:
num_bins = 1250 #Change as desired for histogram
fig203, ax203 = plt.subplots(nrows=3, ncols=2, figsize=(10,12))
axes = [ax203[0,0], ax203[1,0], ax203[2,0]]
titles = ['Forward Switching Current Area',
'Reverse Switching Current Area',
'Switching current Difference']
imgs = [np.abs(integral_area[:,:,0]),
np.abs(integral_area[:,:,1]),
np.abs(integral_area[:,:,0]) -
np.abs(integral_area[:,:,1])]
for axis, title, img in zip(axes, titles, imgs):
axis.set_title(title, fontsize=16)
axis.set_ylabel('Y ($\mu$m)', fontsize=16)
im_hand, cbar_hand = px.plot_utils.single_img_cbar_plot(axis, img, show_xy_ticks=True,
x_size=scan_width, y_size=scan_width, tick_font_size=14,
clim=[-0.5, 2])
axis.set_xlabel('X ($\mu$m)', fontsize=16)
im20301 = ax203[0,1].hist(np.abs(integral_area[:,:,0]).ravel(), num_bins)
ax203[0,1].set_title('Forward Switching Current Area', fontsize=16)
im20301 = ax203[1,1].hist(np.abs(integral_area[:,:,1]).ravel(), num_bins)
ax203[1,1].set_title('Reverse Switching Current Area', fontsize=16)
im20301 = ax203[2,1].hist((np.abs(integral_area[:,:,0]) -np.abs(integral_area[:,:,1])).ravel(), num_bins)
ax203[2,1].set_title('Difference Area', fontsize=16)
for axis in [ax203[0,1], ax203[1,1], ax203[2,1]]:
axis.set_ylabel('Counts', fontsize=16)
ax203[0,1].set_xticks([])
ax203[1,1].set_xticks([])
ax203[0,1].set_xlim((0,2))
ax203[1,1].set_xlim((0,2))
ax203[2,1].set_xlim((0,2))
px.plot_utils.set_tick_font_size([ax203[0,1], ax203[1,1], ax203[2,1]], font_size=14)
fig203.tight_layout()
"""file_name = 'Switching Area Plots'
fig203.savefig(os.path.join(figure_folder, file_name + '.png'), format='png', dpi=300)
fig203.savefig(os.path.join(figure_folder, file_name + '.pdf'), format='pdf', bbox_inches = 'tight',
pad_inches = 2.0)"""
Out[12]:
Visualize Switched Charge¶
In [22]:
fig, axes = twin_img_plots([Q_integrated[:,:,0]*1E12, Q_integrated[:,:,1]*1E12],
[direction + ' Switching Q (pC)' for direction in ['Forward', 'Reverse']],
x_size=scan_width, y_size=scan_width, stdevs=1.25, show_xy_ticks=True)
"""file_name = 'Switching_Q'
fig.savefig(os.path.join(figure_folder, file_name + '.png'), format='png', dpi=300)
fig.savefig(os.path.join(figure_folder, file_name + '.pdf'), format='pdf', bbox_inches = 'tight',
pad_inches = 2.0)"""
Out[22]:
Visualize where the switching peaks are¶
In [23]:
fig, axes = twin_img_plots([pos_of_peak[:,:,0], pos_of_peak[:,:,1]],
[direction + ' Switching Maximum Bias (V)' for direction in ['Forward', 'Reverse']],
x_size=scan_width, y_size=scan_width, stdevs=2, show_xy_ticks=True)
"""file_name = 'Switching_max_bias'
fig.savefig(os.path.join(figure_folder, file_name + '.png'), format='png', dpi=300)
fig.savefig(os.path.join(figure_folder, file_name + '.pdf'), format='pdf', bbox_inches = 'tight',
pad_inches = 2.0)"""
Out[23]:
In [25]:
cap_diameter = 500E-9 #Diameter of cap
cap_area = np.pi*(cap_diameter/2)**2
P_forward = 0.5*(Q_integrated[:,:,0]/(cap_area)) # Polarization in C/m^2,
P_forward = P_forward*1E6*1E-4 #Polarization in uC (*1E6) /cm^2 (*1E-4)
P_reverse = 0.5*(Q_integrated[:,:,1]/(cap_area)) #Polarization in C/m^2
P_reverse = P_reverse*1E6*1E-4 #Polarization in uC (*1E6) /cm^2 (*1E-4)
Now visualize¶
In [24]:
fig, axes = twin_img_plots([np.abs(P_forward), np.abs(P_reverse)],
['Polarization ' + direction + ' ($\mu$C/cm$^2$)' for direction in ['Forward', 'Reverse']],
x_size=scan_width, y_size=scan_width, stdevs=1, show_xy_ticks=True)
"""file_name = 'Polarization_maps'
fig.savefig(os.path.join(figure_folder, file_name + '.png'), format='png', dpi=300)
fig.savefig(os.path.join(figure_folder, file_name + '.pdf'), format='pdf', bbox_inches = 'tight',
pad_inches = 2.0)"""
Out[24]:
In [26]:
fig101, ax101 = plt.subplots(nrows=1,ncols=2, sharey=True, figsize=(10,5))
good_P_forward = np.abs(P_forward).ravel()
good_P_forward = good_P_forward[np.where(good_P_forward>15)]
good_P_forward = good_P_forward[np.where(good_P_forward<120)]
good_P_reverse = np.abs(P_reverse).ravel()
good_P_reverse = good_P_reverse[np.where(good_P_reverse>10)]
good_P_reverse = good_P_reverse[np.where(good_P_reverse<120)]
ax101[0].hist(good_P_forward, 150)
ax101[1].hist(good_P_reverse, 150)
ax101[0].set_xlabel('Polarization Forward ($\mu$C/cm^2)', fontsize=16)
#ax101[0].set_title('$Polarization Forward(\mu C/cm^2)$', fontsize=16)
ax101[1].set_xlabel('Polarization Reverse ($\mu$C/cm^2)', fontsize=16)
#ax101[1].set_title('$Polarization Reverse(\mu C/cm^2)$', fontsize=16)
#ax101[1].set_ylabel('Counts', fontsize=16)
ax101[0].set_ylabel('Counts', fontsize=16)
ax101[0].set_xlim((15,100))
ax101[1].set_xlim((15,100))
px.plot_utils.set_tick_font_size(ax101, 14)
fig101.tight_layout()
"""file_name = 'P_hist'
fig101.savefig(os.path.join(figure_folder, file_name + '.png'), format='png', dpi=300)
fig101.savefig(os.path.join(figure_folder, file_name + '.pdf'), format='pdf', bbox_inches = 'tight',
pad_inches = 2.0)"""
Out[26]:
Understand nature of the switching events¶
We will do this by using the K-means clustering algorithm to determine the shapes of the curves and their spatial abundance
Before performing K-Means, we need to get rid of outliers
In [27]:
xvec = bias_interp[:half_ind]
P_forward_copy = np.copy(np.abs(P_forward).ravel())
switching_points = np.where(P_forward_copy>20)
Pswitching_fwd_lin = np.copy(Pswitching_fwd).reshape(num_lines*num_cols,-1)
P_switching_fwd = Pswitching_fwd_lin[:,:]
xvec =xvec_Pr_fwd[:]# xvec[70:]
# Getting rid of the handful of outliers:
test = np.where(np.abs(P_switching_fwd) > 3)
cleaned_p_switch_fwd = np.copy(P_switching_fwd)
# For simplicity, set the values of these outlier pixels to 0
cleaned_p_switch_fwd[test[0]] = 0
# See the distribution of these bad pixels
fig, axis = plt.subplots()
bad_pix = np.ones(P_switching_fwd.shape[0], dtype=bool)
bad_pix[test[0]] = False
axis.imshow(bad_pix.reshape(num_lines, num_cols), cmap='gray', origin='lower', alpha=0.5)
axis.set_title('Outlier Pixels');
Performing K-means¶
In [28]:
num_clusters = 12
#Do k-means
estimators = KMeans(num_clusters)
results = estimators.fit(cleaned_p_switch_fwd)
labels, cluster_centroids = px.processing.cluster.reorder_clusters(results.labels_,
results.cluster_centers_)
Visualizing K-means results¶
In [51]:
# Plot the kmeans cluster labels
fig_lab, fig_centroids = plot_cluster_results_separately(labels.reshape(num_lines, num_cols),
cluster_centroids, xvec_Pr_fwd, legend_mode=2,
x_size=scan_width, y_size=scan_width)
"""file_name = 'GIV_Kmeans_Disorder_Labels'
fig_lab.savefig(os.path.join(figure_folder, file_name + '.png'), format='png', dpi=300)
fig_lab.savefig(os.path.join(figure_folder, file_name + '.pdf'), format='pdf', bbox_inches = 'tight',
pad_inches = 2.0)
file_name = 'GIV_Kmeans_Disorder_Centroids'
fig_centroids.savefig(os.path.join(figure_folder, file_name + '.png'), format='png', dpi=300)
fig_centroids.savefig(os.path.join(figure_folder, file_name + '.pdf'), format='pdf', bbox_inches = 'tight',
pad_inches = 2.0)"""
Out[51]:
Disorder Analysis¶
Fit the Gaussian at each point ala Shur paper. Look for point where switched fraction = unswitched fraction. Then fit to a Gaussian, with the constraint that the area equals that of the switched fraction
In [ ]:
#Define some functions
def gauss1(x,*p1):
A1,x1,sigma1= p1
return (A1*np.exp(-((x-x1)**2)/(2*sigma1*sigma1)))
# here you include the penalization factor
def residuals(p,x,y):
#Calcualte the area given the whole xvector
fitted_gaussian = gauss1(xvec, *p) #Gaussian, for whole x
integral = np.trapz(fitted_gaussian,xvec) #Area, for whole x
penalization = abs(integral - area_of_Rayleigh)*1E4 #c.f. with area of Rayleigh
return y - gauss1(x, *p) - penalization
coef_disorder_fit = np.zeros(shape=(P_switching_fwd.shape[0],3))
middle_ind_mat= np.zeros(shape=(P_switching_fwd.shape[0],1))
#Fitting guesses
p0 = [1, 4E-4, 1E-4]
lb = [0,0,0]
ub=[1,1,1]
bounds = (lb,ub)
xvec = np.linspace(0,tstep*time_ratio_fwd, len(xvec_Pr_fwd[:])) #xvec
for pix_ind, yvec in enumerate(P_switching_fwd):
#If it's a switching event
if pix_ind in switching_points[0]:
#Calculate the area
#del area_of_Rayleigh
global area_of_Rayleigh
area_of_Rayleigh = np.trapz(yvec, xvec)
cum_area = scipy.integrate.cumtrapz(yvec,xvec)
#Take the point where the area is half of the total as the cutoff point
middle_index = np.argsort((cum_area - 0.5*area_of_Rayleigh)**2)[0]
#Smaller x,yvecs
new_x, new_y = xvec[middle_index:],yvec[middle_index:]
#Apply the area constraint
if len(new_x)>4:
popt2, pcov2 = scipy.optimize.leastsq(func=residuals, x0=p0, args=(new_x,new_y))
y_fit2 = gauss1(xvec, *popt2)
coef_disorder_fit[pix_ind,:] = popt2
middle_ind_mat[pix_ind,:] = middle_index
Visualize the disorder fitting¶
At a specific location
In [34]:
#row, col = 115, 51
row, col = 110, 134
pix_ind = row * num_cols + col
xvec = np.linspace(0,tstep*time_ratio_fwd, len(xvec_Pr_fwd[:]))
yvec = P_switching_fwd[pix_ind,:] #yvec
new_x = xvec[int(middle_ind_mat[pix_ind][0]):]
new_y = yvec[int(middle_ind_mat[pix_ind][0]):]
y_fit2 = gauss1(xvec, *coef_disorder_fit[pix_ind,:] )
xvec *= 1E+3 # milliseconds
fig, axis = plt.subplots(figsize=(5.25,5))
axis.plot(xvec, yvec, 'bo', markersize=6)
axis.plot(new_x, new_y, 'ro', markersize=6)
axis.fill_between(new_x, new_y, 0, where=new_y >= 0, facecolor='green',
alpha= 0.5, edgecolor = 'w', linewidth=0.0)
axis.fill_between(new_x, new_y, 0, where=new_y <= 0, facecolor='green',
alpha= 0.5, edgecolor = 'w', linewidth=0.0)
axis.plot(xvec, y_fit2, linewidth=3, label='Constrained Fit', color='purple')
axis.legend(loc='upper left', bbox_to_anchor=(1, 1), fontsize=14)
axis.set_ylabel('Current (nA)', fontsize=16)
axis.set_xlabel('Time (ms)', fontsize=16)
axis.set_title('Disorder Fitting Example', fontsize=16)
plt.ticklabel_format(axis='x', style='sci', scilimits=(-2,2))
px.plot_utils.set_tick_font_size(axis, 14)
print (coef_disorder_fit[pix_ind,:])
print (pix_ind in switching_points[0])
file_name = 'Disorder_fit_sp'
fig.savefig(os.path.join(figure_folder, file_name + '.png'), format='png', dpi=300)
fig.savefig(os.path.join(figure_folder, file_name + '.pdf'), format='pdf', bbox_inches = 'tight',
pad_inches = 2.0)
Visualize the spatial distribution of the Disorder¶
In [35]:
fig, axis = plt.subplots(figsize=(6, 5))
_, _ = px.plot_utils.single_img_cbar_plot(axis, (coef_disorder_fit[:,2]*1E3).reshape(num_lines, num_cols),
x_size=scan_width, y_size=scan_width, clim=[0, 0.3],
cbar_label='$\sigma\cdot 10^{-3} (t^{-0.5}$)', aspect='auto')
axis.set_xlabel('X ($\mu$m)', fontsize=16)
axis.set_ylabel('Y ($\mu$m)', fontsize=16)
axis.set_title('Disorder Analysis (Forward)', fontsize=16)
fig.tight_layout()
file_name = 'Disorder_fit_map'
fig.savefig(os.path.join(figure_folder, file_name + '.png'), format='png', dpi=300)
fig.savefig(os.path.join(figure_folder, file_name + '.pdf'), format='pdf', bbox_inches = 'tight',
pad_inches = 2.0)
Dielectric Constant¶
In [37]:
e_r = (cap_vec*1E-12 - 0.480E-12)*d/(e0*cap_area)
e_r_good = e_r[np.where(e_r<200)].ravel()
e_r_good = e_r_good[np.where(e_r_good>0)]
fig, axes = plt.subplots(ncols=2, figsize=(11.5, 5))
_ = px.plot_utils.single_img_cbar_plot(axes[0], e_r.reshape(num_lines, num_cols),
x_size=scan_width, y_size=scan_width, clim=[0, 130])
axes[0].set_xlabel('X ($\mu$m)', fontsize=16)
axes[0].set_ylabel('Y ($\mu$m)', fontsize=16)
axes[0].set_title('Calculated Dielectric Constant', fontsize=16)
axes[1].hist(e_r_good, 150);
axes[1].set_title('Dielectric Constant', fontsize=16)
axes[1].set_xlabel('Dielectric Constant', fontsize=16)
axes[1].set_ylabel('Counts', fontsize=16)
px.plot_utils.set_tick_font_size(axes[1], 14)
fig.tight_layout()
file_name = 'D_hist'
fig.savefig(os.path.join(figure_folder, file_name + '.png'), format='png', dpi=300)
fig.savefig(os.path.join(figure_folder, file_name + '.pdf'), format='pdf', bbox_inches = 'tight',
pad_inches = 2.0)
Standard IV (s-IV)¶
Just like the analysis with gIV datasets, we begin analysis of sIV datasets by loading the h5 file where the results are stored.
In [20]:
h5_siv_path = px.io_utils.uiGetFile('*.h5')
print('Working on SIV file:', h5_siv_path)
siv_folder_path, siv_file_name = os.path.split(h5_siv_path)
h5_siv_f = h5py.File(h5_siv_path, mode='r+')
# Extracting some parameters and handles to the necessary datasets:
h5_siv_grp = h5_siv_f['Measurement_000']
h5_siv_raw = h5_siv_f['Measurement_000/Channel_000/Raw_Data']
siv_num_rows = h5_siv_grp.attrs['num_rows']
siv_num_cols = h5_siv_grp.attrs['num_cols']
siv_num_bias_pts = h5_siv_raw.shape[1]
h5_siv_spec_vals = px.hdf_utils.getAuxData(h5_siv_raw, auxDataName=['Spectroscopic_Values'])[0]
siv_excit_wfm = np.squeeze(h5_siv_spec_vals[()])
siv_y_label = 'Current (nA)'# px.hdf_utils.get_data_descriptor(h5_siv_raw)
SVD on SIV¶
In [21]:
h5_siv_svd_grp = px.hdf_utils.findH5group(h5_siv_raw, 'SVD')
if len(h5_siv_svd_grp) > 0:
print('Taking previous SVD results')
h5_siv_svd_grp = h5_siv_svd_grp[-1]
else:
print('FFT filtering not performed on this dataset. Filtering now:')
h5_siv_svd_grp = px.doSVD(h5_siv_raw)
Visualizing variance of SVD components¶
In [46]:
fig_s, axis_s = px.plot_utils.plotScree(h5_siv_svd_grp['S'])
fig_s.savefig(os.path.join(other_figures_folder, 'SIV_SVD_S.png'), format='png', dpi=150)
Visualizing abundance maps of SVD endmembers¶
In [47]:
loadings = np.reshape(h5_siv_svd_grp['U'][:, :16], (siv_num_rows, siv_num_cols, -1))
fig_u, axes_u = px.plot_utils.plot_map_stack(loadings, num_comps=loadings.shape[-1], heading='Abundance Maps')
fig_u.savefig(os.path.join(other_figures_folder, 'SIV_SVD_U.png'), format='png', dpi=150)
Visualizing endmembers from SVD¶
These are statistically the most important (orthogonal) trends in the IV curves
In [48]:
fig_v, axes_v = px.plot_utils.plot_loops(siv_excit_wfm, h5_siv_svd_grp['V'], evenly_spaced=False, plots_on_side=4,
x_label='Bias (V)', y_label=siv_y_label,
subtitles='Component', title='SVD Endmembers')
fig_v.savefig(os.path.join(other_figures_folder, 'SIV_SVD_V.png'), format='png', dpi=150)
In [32]:
"""h5_siv_filt = px.processing.svd_utils.rebuild_svd(h5_siv_raw, components=16)"""
fig, axes = px.plot_utils.plot_loops(siv_excit_wfm, [h5_siv_raw, h5_siv_filt], line_colors=['red','k'],
dataset_names=['Raw', 'Filtered'], x_label='Bias (V)',
y_label='Current (nA)', title='SVD Filtering of SIV')
k-Means for understanding principle responses and spatial distribution¶
In [27]:
num_clusts = 6
h5_siv_kmeans_grp = px.hdf_utils.check_for_old(h5_siv_raw, 'Cluster',
{'n_clusters':num_clusts, 'n_jobs':1})
if h5_siv_kmeans_grp is None:
print('Performing K-means now')
clusterer = px.processing.Cluster(h5_siv_raw, 'KMeans', n_clusters=num_clusts)
h5_siv_kmeans_grp = clusterer.do_cluster()
else:
print('Taking previous results')
# fig_km, axes_km = px.plot_utils.plot_cluster_h5_group(h5_siv_kmeans_grp, siv_y_label)
fig, axes = plt.subplots(ncols=2, figsize=(11.5,5))
single_img_cbar_plot(fig, axes[0], np.reshape(h5_siv_kmeans_grp['Labels'][()],(siv_num_rows, siv_num_cols)),
x_size=scan_width, y_size=scan_width, clim=[0, num_clusts-1], aspect='auto')
axes[0].set_xlabel('X ($\mu$m)', fontsize=16)
axes[0].set_ylabel('Y ($\mu$m)', fontsize=16)
axes[0].set_title('K-Means Cluster Labels', fontsize=16)
colors = [ plt.cm.viridis(x) for x in np.linspace(0, 1, num_clusts)]
for line_ind in range(num_clusts):
axes[1].plot(siv_excit_wfm, h5_siv_kmeans_grp['Mean_Response'][line_ind], #marker ='o',
label='Cluster ' + str(line_ind), color=colors[line_ind], linewidth=2)
axes[1].set_xlabel('Voltage (V)', fontsize=16)
axes[1].set_ylabel('Current (nA)', fontsize=16)
axes[1].set_title('K-Means Cluster Centroids', fontsize=16)
# axes[1].legend(loc='lower right', fontsize=14)
axes[1].legend(loc='upper left', bbox_to_anchor=(1, 1.12), fontsize=14)
px.plot_utils.set_tick_font_size(axes[1], 14)
fig.tight_layout()
# fig.savefig(os.path.join(other_figures_folder, 'siv_kmeans.png'), format='png', dpi=200)
Visualizing the raw SIV curve at a given pixel¶
In [49]:
spectral_point = 170
point_siv = [24,39]
fig, axes = plt.subplots(ncols=2, figsize=(11.5, 5))
single_img_cbar_plot(fig, axes[0], np.reshape(h5_siv_raw[:,spectral_point],(siv_num_rows, siv_num_cols)),
clim=[-1.5,1.5],
x_size=scan_width, y_size=scan_width) #, cbar_label='Current (nA)')
axes[0].set_title('Current Map, V = ' + str(round(siv_excit_wfm[spectral_point],2)) + 'V',
fontsize=16)
axes[0].set_xlabel('X ($\mu$m)', fontsize=16)
axes[0].set_ylabel('Y ($\mu$m)', fontsize=16)
axes[0].plot(point_siv[1],point_siv[0],marker='o', color='r')#, alpha = 0.8)
siv_pixel = siv_num_cols * point_siv[0] + point_siv[1]
axes[1].plot(siv_excit_wfm, h5_siv_raw[siv_pixel], 'r-', linewidth=2)
axes[1].set_xlabel('Voltage (V)', fontsize=16)
axes[1].set_ylabel ('Current (nA)', fontsize=16)
axes[1].set_title('S I-V Point Measurement', fontsize=16)
fig.tight_layout()
fig.savefig(os.path.join(figure_folder, 'IVPoint.png'), format='png', dpi=200)
fig.savefig(os.path.join(figure_folder, 'IVPoint.eps'), format='eps', dpi=200)
Comparing G-IV with the standard IV at a single location¶
In [51]:
point_fiv = [120,245] #For the fast IV
point_siv = [24,39] #For the standard IV
point_fiv_film = [28,18] #For the fast IV
point_siv_film = [30,18] #For the standard IV
cap_giv_fwd = currents[point_fiv[0],point_fiv[1],:half_ind] #- R2*cap_vec[pix_index,0]
cap_giv_rev = currents[point_fiv[0],point_fiv[1],half_ind:] #- R2*cap_vec[pix_index,1]
film_giv_fwd = currents[point_fiv_film[0],point_fiv_film[1],:half_ind]
film_giv_rev = currents[point_fiv_film[0],point_fiv_film[1],half_ind:]
fig, axis = plt.subplots(figsize=(5.5,5.5))
siv_pixel = siv_num_cols * point_siv[0] + point_siv[1]
axis.plot(siv_excit_wfm, h5_siv_raw[siv_pixel], 'r-', linewidth=2, label = 'Cap S-IV')
axis.plot(bias_interp[:half_ind], cap_giv_fwd, 'b-', linewidth=2, label = 'Cap G-IV')
axis.plot(bias_interp[half_ind:],cap_giv_rev, 'b-', linewidth=2)
siv_pixel = siv_num_cols * point_siv_film[0] + point_siv_film[1]
axis.plot(siv_excit_wfm, h5_siv_raw[siv_pixel], 'g-', linewidth=2, label = 'Film S-IV')
axis.plot(bias_interp[:half_ind], film_giv_fwd, '-', color='orange', linewidth=2, label = 'Film G-IV')
axis.plot(bias_interp[half_ind:], film_giv_rev, '-', color = 'orange', linewidth=2)
axis.set_xlabel('$Voltage (V)$', fontsize=14)
axis.set_ylabel ('$Current (nA)$', fontsize=14)
axis.set_title('I-V Point Comparison', fontsize=16)
axis.legend(loc='best', fontsize=14)
px.plot_utils.set_tick_font_size(axis, 14)
file_name = 'Comparison_IVPoint'
fig.savefig(os.path.join(figure_folder, file_name + '.png'), format='png', dpi=300)
fig.savefig(os.path.join(figure_folder, file_name + '.pdf'), format='pdf', bbox_inches = 'tight',
pad_inches = 2.0)
Register datasets¶
There was some small spatial drift between the gIV and sIV datasets. We therefore need to register these datasets such that we are able to compare the same areas and pixels.
The gIV dataset has 256 rows and 266 columns while the sIV dataset has 40 rows and columns. It is therefore important to either decimate or interpolate the images such that the pixel-to-pixel correspondance is maintained as accurately as possible.
Also, by nature of their implementation, gIV introduces substantially lesser spatial drift within the measurement since the gIV measurement takes around 10-20 mins while sIV takes a few hours.
In [29]:
giv_img = currents[:,:,128]
siv_img = np.reshape(h5_siv_raw[:, 171], (siv_num_rows, siv_num_cols))
#siv_interp = scipy.misc.imresize(siv_img, giv_img.shape)
siv_interp = scipy.ndimage.interpolation.zoom(siv_img,
(giv_img.shape[0]/siv_img.shape[0],
giv_img.shape[1]/siv_img.shape[1]))
giv_points = [(190, 63), (111, 236), (157, 63), (184, 225)]
siv_points = [(190, 17), (111, 182), (157, 19), (190, 175)]
Registration points¶
Here we compare one current slice from each technique and a few features common to both images.
In [31]:
fig, axes = twin_img_plots([giv_img, siv_interp], ['g-IV','s-IV'],
x_size=scan_width, y_size=scan_width, stdevs=1)
colors = ['r','w','orange','b']
for coord, col in zip(giv_points, colors):
axes[0].scatter(coord[1], coord[0], color=col, s=100)
axes[1].imshow(siv_interp, origin='lower')
for coord, col in zip(siv_points, colors):
axes[1].scatter(coord[1], coord[0], color=col, s=100)
Alignment¶
While one could use rigorous image transform algorithms to register one image to another, we chose to simply interpolate and crop images because:
- the sIV image has only moved a little to the left compared to the gIV image.
- The signal to noise ratio in the sIV dataset is too low for registration algorithms
- Only a handful of common capacitors are visible in the s-IV images.
The plots below will show that such simple approaches are sufficent in our case since we can get the two images to overlap very well
In [37]:
# normalizing
siv_norm = (siv_interp - np.amin(siv_interp))/(np.amax(siv_interp)-np.amin(siv_interp))
giv_norm = (giv_img - np.amin(giv_img))/(np.amax(giv_img)-np.amin(giv_img))
# cropping
x_offset = 45
cropped_giv_img = giv_norm[:, x_offset:]
cropped_siv_img = siv_norm[:, :-x_offset]
# Visualizing:
fig, axis = plt.subplots(figsize=(5.5,6))
axis.imshow(cropped_giv_img, origin='lower', cmap='Blues', clim=[0.3, 0.7], aspect='auto')
axis.imshow(cropped_siv_img, cmap = 'gray', origin='lower', alpha = 0.5, clim=[0.85, 0.95], aspect='auto')
axis.set_title('Overlap of aligned sIV and gIV images', fontsize=16)
fig.tight_layout()
Comparing gIV and sIV current slices¶
In [39]:
fig, axes = plt.subplots(nrows=2,ncols=2, figsize=(12,10))
cropped_x_size = scan_width*(currents.shape[1]-45)/currents.shape[1]
for axis, img, title in zip(axes.flat[:2],
[currents[:,:,120], currents[:,:,128]],
['GIV $I$ at $V = $' + str(round(bias_interp[ind],2)) + ' V' for ind in [121, 128]]):
cropped_giv_img = img[:, 45:]
single_img_cbar_plot(fig, axis, cropped_giv_img,
x_size=cropped_x_size, y_size=scan_width, stdevs=2, cbar_label='Current (nA)')
axis.set_title(title, fontsize=16)
for axis, img, title in zip(axes.flat[2:],
[np.reshape(h5_siv_raw[:, ind], (siv_num_rows, siv_num_cols)) for ind in [57, 171]],
['SIV $I$ at $V = $' + str(round(siv_excit_wfm[ind],2)) + ' V' for ind in [57, 171]]):
"""
siv_interp = zoom(img, (currents.shape[0]/img.shape[0], currents.shape[1]/img.shape[1]))
cropped_siv_img = siv_interp[:, :-45]
"""
cropped_siv_img = img[:, :-9]
single_img_cbar_plot(fig, axis, cropped_siv_img,
x_size=cropped_x_size, y_size=scan_width, stdevs=2, cbar_label='Current (nA)')
axis.set_title(title, fontsize=16)
for axis in [axes.flat[0], axes.flat[2]]:
axis.set_ylabel('Y ($\mu$m)', fontsize=16)
for axis in axes.flat[2:]:
axis.set_xlabel('X ($\mu$m)', fontsize=16)
file_name = 'SIV_GIV_current_slices_registered'
fig.savefig(os.path.join(figure_folder, file_name + '.png'), format='png', dpi=300)
fig.savefig(os.path.join(figure_folder, file_name + '.pdf'), format='pdf', bbox_inches = 'tight',
pad_inches = 2.0)
Current maps at switching biases¶
In [65]:
cropped_x_size = scan_width*(currents.shape[1]-45)/currents.shape[1]
targ_volts = [-3.75, 3.75]
neg_inds = np.argwhere(np.abs(bias_interp - targ_volts[0]) < 0.05)
pos_inds = np.argwhere(np.abs(bias_interp - targ_volts[1]) < 0.05)
volt_inds = np.squeeze(np.array([neg_inds[0], pos_inds[0], neg_inds[1], pos_inds[1]]))
giv_slices = [currents[:,:,ind] for ind in volt_inds]
chosen_data = np.array(giv_slices)
chosen_avg = np.mean(chosen_data)
chosen_std = np.std(chosen_data)
stdevs = 1.75
clims = [chosen_avg - stdevs * chosen_std,
chosen_avg + stdevs * chosen_std]
for cmap in ['viridis', 'bwr']:
fig, axes = plt.subplots(nrows=2,ncols=2, figsize=(12,10))
for axis, ind, img, direction in zip(axes.flat, volt_inds, giv_slices,
['Forw','Forw','Rev','Rev']):
cropped_giv_img = img[:, 45:]
single_img_cbar_plot(fig, axis, cropped_giv_img,
x_size=cropped_x_size, y_size=scan_width,
clim=clims, cbar_label='Current (nA)',
cmap=cmap)
axis.set_title('GIV $I_{' + direction + '}$ at $V = $' + str(round(bias_interp[ind],2)) +
' V', fontsize=16)
file_name = 'GIV_current_slices_forw_rev_' + cmap
fig.savefig(os.path.join(figure_folder, file_name + '.png'), format='png', dpi=300)
fig.savefig(os.path.join(figure_folder, file_name + '.pdf'), format='pdf', bbox_inches = 'tight')
gIV current maps for all biases¶
In [160]:
cur_mean = np.mean(currents)
cur_std = np.std(currents)
stdevs=2
clims = [cur_mean - stdevs * cur_std,
cur_mean + stdevs * cur_std]
cropped_giv_img = currents[:, 45:, 0]
fig, axis = plt.subplots(figsize=(6,5))
im_handle, cbar_handle = single_img_cbar_plot(fig, axis, cropped_giv_img,
x_size=cropped_x_size, y_size=scan_width, clim=clims, cbar_label='Current (nA)')
axis.set_ylabel('Y ($\mu$m)', fontsize=16)
axis.set_xlabel('X ($\mu$m)', fontsize=16)
for v_ind in range(currents.shape[2]):
cropped_giv_img = currents[:, 45:, v_ind]
im_handle.set_data(cropped_giv_img)
bias_val = str(round(bias_interp[v_ind],2))
axis.set_title('GIV $I$ at $V = $' + bias_val + ' V', fontsize=16)
file_name = 'giv_cropped_current_slice_' + str(v_ind) + '_' + bias_val
fig.savefig(os.path.join(other_figures_folder, file_name + '.png'), format='png', dpi=300)
SIV current maps for all biases¶
In [74]:
siv_currents = np.reshape(h5_siv_raw[()], (siv_num_rows, siv_num_cols, h5_siv_raw.shape[1]))
cur_mean = np.mean(siv_currents)
cur_std = np.std(siv_currents)
stdevs=1.5
clims = [cur_mean - stdevs * cur_std,
cur_mean + stdevs * cur_std]
cropped_siv_img = currents[:, :-9, 0]
fig, axis = plt.subplots(figsize=(6,5))
im_handle, cbar_handle = single_img_cbar_plot(fig, axis, cropped_siv_img,
x_size=cropped_x_size, y_size=scan_width, clim=clims, cbar_label='Current (nA)')
axis.set_ylabel('Y ($\mu$m)', fontsize=16)
axis.set_xlabel('X ($\mu$m)', fontsize=16)
for v_ind in range(siv_currents.shape[2]):
cropped_siv_img = siv_currents[:, :-9, v_ind]
im_handle.set_data(cropped_siv_img)
bias_val = str(round(siv_excit_wfm[v_ind],2))
axis.set_title('SIV $I$ at $V = $' + bias_val + ' V', fontsize=16)
file_name = 'siv_cropped_current_slice_' + str(v_ind) + '_' + bias_val
fig.savefig(os.path.join(other_figures_folder, file_name + '.png'), format='png', dpi=300)
In [305]:
# Done processing the data. Close the file
h5_f.close()
h5_siv_f.close()
In [ ]: