Nature Communications, 6th November 2017
Link to paper
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
Oak Ridge National Laboratory, Oak Ridge TN 37831, USA
This Jupyter notebook uses pycroscopy to analyze Band Excitation data. We request you to reference the following papers if you use this notebook for your research:
This is a Jupyter Notebook - it contains text and executable code cells
. To learn more about how to use it, please see this video. Please see the image below for some basic tips on using this notebook.
Image courtesy of Jean Bilheux from the neutron imaging GitHub repository.
Note: This notebook was written for the pycroscopy version listed below and is not guaranteed to work on past or future versions of the package
pycroscopy
version: 0.59.4
If you have a different version of pycroscopy
installed, you may consider using the notebook as is and accept the possibility of errors. The cell below will attempt to install the correct versions of the packages. However, if you experience trouble, uninstall the existing version of pycroscopy and install the required version above by executing the following commands in a terminal (Linux / MacOS) / Anaconda prompt (Windows):
pip uninstall pycroscopy
pip install -I pycroscopy==0.59.4
!pip install h5py matplotlib numpy scipy sklearn
!pip install pycroscopy==0.59.4
# 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.''')
# 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
# Miscellaneous pacakges
import os
import sys
import time as tm
# For the main processing, analysis, I/O and visualization functions
import pycroscopy as px
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.plot_map(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.plot_map(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
Also setting up some basic parameters
# Load the raw dataset
h5_path = px.io_utils.uiGetFile('*.h5', 'G-mode IV dataset')
if not os.path.exists(h5_path):
raise FileExistsError('Please check your file path. This file does not exist')
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')
Also loading ancillary datasets and necessary parameters
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]']
Here we will inspect what the IV curves look like for a single scan line.
Note the significant hysteresis in the IV curves.
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)', subtitle_prefix='Pixel', title=None)
fig.savefig(os.path.join(other_figures_folder,
'example_raw_curves_line_' + str(row_ind) +'.png'),
format='png', dpi=150);
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
# 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);
# Set Filter parameters here:
num_spectral_pts = h5_main.shape[1]
frequency_filters = [px.processing.fft.LowPassFilter(num_spectral_pts, samp_rate, 10E+3)]
noise_tol = 1E-6
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
# Test filter on a single line:
row_ind = 50
raw_row = h5_main[row_ind]
filt_line, fig_filt, axes_filt = px.processing.gmode_utils.test_filter(raw_row, frequency_filters, noise_threshold=noise_tol,
show_plots=True)
raw_row = np.reshape(raw_row, (-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)',
subtitle_prefix='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)
Once the optimal filter parameters are chosen, the same parameters can be applied to filter the entire dataset. This can take a few minutes.
sig_filt = px.SignalFilter(h5_main, frequency_filters=frequency_filters,
noise_threshold=noise_tol, write_filtered=True,
write_condensed=False, num_pix=1)
h5_filt_grp = sig_filt.compute()
h5_filt = h5_filt_grp['Filtered_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
# 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]
h5_resh = h5_resh_grp['Reshaped_Data']
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_grp = h5_resh.parent
print('Data was reshaped from shape', h5_filt.shape,
'reshaped to ', h5_resh.shape)
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:
col = 166
row = 116
pix_ind = row * num_cols + col
i_meas = h5_resh[pix_ind]
results = px.analysis.utils.giv_utils.bayesian_inference_on_period(i_meas, single_AO, ex_freq,
r_extra=110, num_x_steps=250,
show_plots=True)
Since Bayesian inference is computationally expensive, we have an option to load existing results in the data file
# Load bayesian inference
amp_gain = 9
i_cleaner = px.GIVBayesian(h5_resh, ex_freq, amp_gain, r_extra=110, num_x_steps=250, cores=3,
max_mem_mb=1024*10, verbose=True)
h5_bayesian_grp = i_cleaner.compute()
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']
#h5_f.close()
We will look at the inference results from five randomly chosen pixels
# Visualize Bayesian Inference results
num_rand_pos = 2
chosen_rows = np.random.randint(0, num_lines, num_rand_pos)
chosen_cols = np.random.randint(0, num_cols, num_rand_pos)
chosen_cols = [105,48]
chosen_rows = [5,18]
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]
i_meas = h5_resh[pix_ind]
results = px.analysis.utils.giv_utils.bayesian_inference_on_period(i_meas, single_AO, ex_freq,
r_extra=110, num_x_steps=250,
show_plots=False)
fig = px.analysis.utils.giv_utils.plot_bayesian_results(single_AO, i_meas, results['IcorrSine'],
results['x'], results['mR'], results['vR'],
i_recon=results['Irec'], broken_resistance=True,
pix_pos=[chosen_rows[index], chosen_cols[index]],
r_max=None, res_scatter=False)
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(other_figures_folder, file_name + '.pdf'), format='pdf', bbox_inches = 'tight',
pad_inches = 2.0)
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="single", evenly_spaced=True,
title=axis_subtitles, heading='Current at different voltages', reverse_dims=True,
pad_mult=(0.15, 0.07), colorbar_label = 'Current (nA)')
fig_name = "current_slices"
fig.savefig(os.path.join(other_figures_folder, fig_name + '.png'), format='png', dpi=300)
fig.savefig(os.path.join(other_figures_folder, fig_name + '.pdf'), format='pdf', dpi=300)
We understand the main trends in the current data via K-Means clustering
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)
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)
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)
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="single", evenly_spaced=True,
title=axis_subtitles, heading='Resistance at different voltages', stdevs=0.05,
colorbar_label = 'Resistance (GOhm)')
fig.subplots_adjust( wspace = 0.8)
fig_name = "resistance_slices"
fig.savefig(os.path.join(other_figures_folder, fig_name + '.png'), format='png', dpi=300)
fig.savefig(os.path.join(other_figures_folder, fig_name + '.pdf'), format='pdf', dpi=300)
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))
_ = px.plot_utils.plot_map(fig, axis, np.reshape(cap_vec, (num_lines, num_cols)),
x_size=scan_width, y_size=scan_width)
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)"""
col = 166
row = 117
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)))
font_size_1 = 14
font_size_2 = 16
fig201,ax201 = plt.subplots(nrows=3, figsize=(6,22))
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')
lg = ax201[0].legend(fontsize=font_size_1, bbox_to_anchor = [1.0,1.0])
lg.get_frame().set_linewidth(0)
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<0.80],yvec_fwd_switched[yvec_fwd_switched<0.80],
'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')
lg = ax201[1].legend(bbox_to_anchor = [1.0,1.0], fontsize=font_size_1)
lg.get_frame().set_linewidth(0)
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)
lg = ax201[2].legend(bbox_to_anchor = [1.0,1.0], fontsize=font_size_1)
lg.get_frame().set_linewidth(0)
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)
Polarization is roughly equal to the area under the 'ears' arising from the switching current
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}
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=1,
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')
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()
for row_ind in range(num_lines):
for col_ind in range(num_cols):
ret_dict = get_polarization(resistances[row_ind,col_ind],bias_interp, ex_freq, neg_bias, pos_bias, 0.5 )
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')
Via the area under the switching current curves
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.plot_map(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)"""
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)"""
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)"""
cap_diameter = 500E-9 # Diameter of cap in m
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)
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, clim=[0,60], show_xy_ticks=True,
cbar_label = 'Polarization ($\mu C$ $cm^{-2}$)')
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)
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)"""
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
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');
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_)
# 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)
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
# 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
At a specific location
# 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)
axis.set_ylim([-0.05, 0.62])
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)
fig, axis = plt.subplots(figsize=(6, 5))
_, _ = px.plot_utils.plot_map(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)
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.plot_map(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)
Just like the analysis with gIV datasets, we begin analysis of sIV datasets by loading the h5 file where the results are stored.
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)
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:')
proc = px.SVD(h5_siv_raw)
h5_siv_svd_grp = proc.compute()
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)
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)
These are statistically the most important (orthogonal) trends in the IV curves
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)
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')
num_clusts = 6
estimator = KMeans(n_clusters=num_clusts)
proc = px.Cluster(h5_siv_raw, estimator)
if proc.duplicate_h5_groups is None:
print('Performing K-means now')
h5_siv_kmeans_grp = proc.compute()
else:
print('Taking previous results')
h5_siv_kmeans_grp = proc.duplicate_h5_groups[-1]
# 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))
_ = px.plot_utils.plot_map(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)
spectral_point = 170
point_siv = [24,39]
fig, axes = plt.subplots(ncols=2, figsize=(11.5, 5))
_ = px.plot_utils.plot_map(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)
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)
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.
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)]
Here we compare one current slice from each technique and a few features common to both images.
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)
While one could use rigorous image transform algorithms to register one image to another, we chose to simply interpolate and crop images because:
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
# 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()
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:]
_ = px.plot_utils.plot_map(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]
_ = px.plot_utils.plot_map(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)
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:]
_ = px.plot_utils.plot_map(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')
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 = px.plot_utils.plot_map(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_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 = px.plot_utils.plot_map(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)
# Done processing the data. Close the file
h5_f.close()
h5_siv_f.close()