In [8]:
%load_ext autoreload
%autoreload 2
In [9]:
import IPython
from scipy import io
import pandas as pd
import h5py
import matplotlib.pyplot as plt
from matplotlib import cm, ticker, gridspec
import gen_matrix_2d as gm2
import numpy as np
import cmocean
%matplotlib inline
In [10]:
# Define figure size
figsize=[10,8]
# Scientific notation format
# From: http://stackoverflow.com/questions/25983218/scientific-notation-colorbar-in-matplotlib
def fmt(x, pos):
a, b = '{:.2e}'.format(x).split('e')
b = int(b)
return r'${} \times 10^{{{}}}$'.format(a, b)
In [11]:
# Load data
data = pd.read_hdf('conv_full_80/conv.hdf')
# Remove failed methods
data = data[data['info']==0]
# Extract classes
krylov = data[data['class']=='krylov']
direct = data[data['class']=='direct']
# Import solution vectors
sol = h5py.File('conv_full_80/conv_sol.hdf5','r')
# Determine number of spatial & angular resolutions used in trials
nsg = len(np.unique([key.split('x')[0] for key in sol.keys()]))
nag = len(np.unique([key.split('x')[1] for key in sol.keys()]))
In [12]:
ls
In [13]:
np.unique([key.split('x')[0] for key in sol.keys()])
Out[13]:
In [14]:
np.unique([key.split('x')[1] for key in sol.keys()])
Out[14]:
In [15]:
# List of matrix sizes
list(sol.keys())
Out[15]:
In [16]:
# List of methods
list(sol['10x10'].keys())
Out[16]:
In [17]:
# Show a subset of data table
ns=40
nth=40
whichrows = np.logical_and.reduce([
data['info']==0,
data['nx']==ns,
data['ny']==ns,
data['nth']==nth
])
data[whichrows]
Out[17]:
In [18]:
# Compare irradiance plots for two different methods
ns = 40
nth = 40
method_list = ['lu','lgmres']
for method in method_list:
plt.figure(figsize=figsize)
gm2.plot_irrad(
sol['{}x{}'.format(ns,nth)][method][:],
var_lengths=[ns,ns,nth],
var_order=[0,1,2],
title=method,
cmap=cmocean.cm.thermal
)
In [19]:
list(sol['10x10'].keys())
Out[19]:
In [20]:
#sol.close()
In [21]:
# fig = plt.figure(figsize=[10,6])
# ax = plt.gca()
# srtd = stat_it.sort_values('mat_size')
# srtd.groupby('method').plot(x='mat_size',y='dt',ax=ax,marker='o')
# ax.set_yscale('log')
# plt.xlabel('Matrix size (m=n)')
# plt.ylabel('Log time to converge (log s)')
# plt.title('Convergence times: Stationary iterative methods')
# plt.legend([v[0] for v in srtd.groupby('method')['method']])
# plt.savefig('plots/stat_it_dt_log.png')
#
# fig = plt.figure(figsize=[10,6])
# ax = plt.gca()
# srtd.groupby('method').plot(x='mat_size',y='dt',ax=ax,marker='o')
# plt.xlabel('Matrix size (m=n)')
# plt.ylabel('Time to converge (s)')
# plt.title('Convergence times: Stationary iterative methods')
# plt.legend([v[0] for v in srtd.groupby('method')['method']])
# plt.savefig('plots/stat_it_dt.png')
In [22]:
fig = plt.figure(figsize=figsize)
ax = plt.gca()
srtd = krylov.sort_values('mat_size')
srtd.groupby('method').plot(x='mat_size',y='dt',ax=ax,marker='o')
ax.set_yscale('log')
plt.xlabel('Matrix size (m=n)')
plt.ylabel('Log time to converge (log s)')
plt.title('Convergence times: Krylov subspace methods')
plt.legend([v[0] for v in srtd.groupby('method')['method']])
plt.savefig('plots/krylov_dt_log.png')
fig = plt.figure(figsize=figsize)
ax = plt.gca()
srtd.groupby('method').plot(x='mat_size',y='dt',ax=ax,marker='o')
plt.xlabel('Matrix size (m=n)')
plt.ylabel('Time to converge (s)')
plt.title('Convergence times: Krylov subspace methods')
plt.legend([v[0] for v in srtd.groupby('method')['method']])
plt.savefig('plots/krylov_dt.png')
In [23]:
# Isolate LGMRES results
(
krylov
[krylov['method']=='lgmres']
.sort_values(['nx','nth'])
)
Out[23]:
In [24]:
# Extract dt values for lgmres to 2d array
# Indexed by [ns,na]
dt_grid = (
krylov
[krylov['method']=='lgmres']
.sort_values(['nx','nth'])
['dt']
.values
.reshape([nsg,nag])
)
# Create meshgrids for ns, na
ns_list = np.sort(np.unique(krylov['nx']))
na_list = np.sort(np.unique(krylov['nth']))
ns_grid,na_grid = np.meshgrid(ns_list,na_list,indexing='ij')
In [25]:
cmap = cm.viridis_r
plt.figure(figsize=figsize)
#plt.contour(ns_grid,na_grid,np.log(dt_grid),cmap=cmap)
plt.scatter(ns_grid,na_grid,c=np.log(dt_grid),
cmap=cmap,edgecolor='k',s=100)
plt.colorbar(label='Solution time')
plt.xlabel('Spatial grid points')
plt.ylabel('Angular grid points')
plt.title('LGMRES solution times for various grids')
# Add text above dots
vshift = max(na_list)/25
for ii in range(nsg):
ns = ns_list[ii]
for jj in range(nag):
na = na_list[jj]
dt = dt_grid[ii][jj]
plt.text(ns,na+vshift,'{:.2e}'.format(dt),ha='center',va='bottom')
plt.tight_layout()
plot_margin = vshift
x0, x1, y0, y1 = plt.axis()
plt.axis((x0 - plot_margin,
x1 + plot_margin,
y0 - plot_margin,
y1 + plot_margin * 2))
Out[25]:
In [26]:
# Create list of lists from solutions for method
# method_name is the string identifying the method
# sol is a dict-like object containing the solution
# sol is indexed by matrix name, then method name
def method_to_grid(method_name,sol,ns_list,na_list):
# Determine nsg, nag
nsg = len(ns_list)
nag = len(na_list)
grid = []
for ii in range(nsg):
grid.append([])
for jj in range(nag):
grid[ii].append(sol['{:.0f}x{:.0f}'
.format(ns_list[ii],na_list[jj])]
[method_name]
[:]
)
return grid
In [27]:
# Create grid of plots for a particular method
# qq = radiance quantity to plot (2D array of shape (ns,na))
def irrad_grid(qq,ns_list,na_list,method,var_order=[0,1,2],
cmap=cmocean.cm.thermal,title=None):
# Determine nsg, nag
nsg = len(ns_list)
nag = len(na_list)
# Axis grid
# Use GridSpec to allow margins to be adjusted nicely
# without title overlapping column labels
# See: http://stackoverflow.com/questions/8248467/matplotlib-tight-layout-doesnt-take-into-account-figure-subtitle
fig = plt.figure(figsize=[12,12.5])
gs = gridspec.GridSpec(nsg,nag)
ax_list = [fig.add_subplot(ss) for ss in gs]
# Font size for row & column labels
fs = 20
# Axes counter
kk = 0
for ii in range(nsg):
ns = int(ns_list[ii])
for jj in range(nag):
na = int(na_list[jj])
var_lengths = [ns,ns,na]
# Row & Column labels
if jj == 0:
ylabel = 'NS = {}'.format(ns)
else:
ylabel = ''
if ii == nsg-1:
xlabel = 'NA = {}'.format(na)
else:
xlabel = ''
title = ''
plt.sca(ax_list[kk])
gm2.plot_irrad(qq[ii][jj],var_lengths,var_order,cmap=cmap,
cbar_kwargs=dict(
format=ticker.FuncFormatter(fmt)))
plt.colorbar()
ax_list[kk].set_xlabel(xlabel,fontsize=fs)
ax_list[kk].set_ylabel(ylabel,fontsize=fs)
ax_list[kk].set_title(title,fontsize=fs)
ax_list[kk].set_xticks([])
ax_list[kk].set_yticks([])
# Increment axes counter
kk += 1
# Figure title
plt.suptitle('{} Convergence Grid'.format(method),fontsize=fs+4)
plt.subplots_adjust()
# Reduce margins
#gs.tight_layout(fig, rect=[0,0,1,1])
gs.tight_layout(fig, rect=[0,0,1,0.95])
In [28]:
g=method_to_grid('lgmres',sol,ns_list,na_list)
irrad_grid(g,ns_list,na_list,'LGMRES')
In [30]:
g1=method_to_grid('lgmres',sol,ns_list,na_list)
g2=method_to_grid('lu',sol,ns_list,na_list)
In [31]:
gdiff = [[np.abs(g1[ii][jj]-g2[ii][jj][:,0]) for jj in range(len(g1[0]))] for ii in range(len(g1))]
In [32]:
irrad_grid(gdiff,ns_list,na_list,'LGMRES - LU')
In [33]:
maxdiff = np.abs(s1-s2[:,0]).max()
gm2.plot_irrad(s1-s2[:,0],[40,40,40],[0,1,2],
cmap=cmocean.cm.thermal,
cbar_kwargs=dict(
format=ticker.FuncFormatter(fmt)),
title='LGMRES - LU'
)
In [34]:
# Contour plots of convergence time
# Reshape times into 2D array
grid =
fig = plt.figure(figsize=[10,6])
ax = plt.gca()
srtd = krylov.sort_values('mat_size')
srtd.groupby('method').plot(x='mat_size',y='dt',ax=ax,marker='o')
ax.set_yscale('log')
plt.xlabel('Matrix size (m=n)')
plt.ylabel('Log time to converge (log s)')
plt.title('Convergence times: Krylov subspace methods')
plt.legend([v[0] for v in srtd.groupby('method')['method']])
plt.savefig('plots/krylov_dt_log.png')
In [35]:
fig = plt.figure(figsize=[10,6])
ax = plt.gca()
srtd = direct.sort_values('mat_size')
srtd.groupby('method').plot(x='mat_size',y='dt',ax=ax,marker='o')
ax.set_yscale('log')
plt.xlabel('Matrix size (m=n)')
plt.ylabel('Log time to converge (log s)')
plt.title('Convergence times: Direct methods')
plt.legend([v[0] for v in srtd.groupby('method')['method']])
plt.savefig('plots/direct_dt_log.png')
fig = plt.figure(figsize=[10,6])
ax = plt.gca()
srtd.groupby('method').plot(x='mat_size',y='dt',ax=ax,marker='o')
plt.xlabel('Matrix size (m=n)')
plt.ylabel('Time to converge (s)')
plt.title('Convergence times: Direct methods')
plt.legend([v[0] for v in srtd.groupby('method')['method']])
plt.savefig('plots/direct_dt.png')
In [ ]: