import numpy as np
import lepm.lattice_elasticity as le
import lepm.le_geometry as leg
import lepm.structure as lestructure
import lepm.plotting.colormaps as cmaps
import lepm.plotting.science_plot_style as sps
##################
import os
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib.collections import LineCollection
from matplotlib.collections import PatchCollection
import matplotlib.gridspec as gridspec
import argparse
import matplotlib.path as mplpath
import matplotlib.patches as mpatches
from mpl_toolkits.axes_grid1 import make_axes_locatable
import glob
try:
import cPickle as pickle
except ImportError:
import pickle
import copy
from scipy.ndimage import imread
from scipy import fftpack
import scipy
from scipy.spatial import Delaunay
import scipy.optimize as opt
import sys
import plotting.colormaps as cmaps
import math
import cmath
import shapely.geometry as sg
'''
Generate lattices using the lattice class: hexagonal, deformed kagome, triangular, square, hexagonalperBC, etc.
Differences between my code and Lisa's
*Ni --> NL
*Nk --> KL
KL : Correponds to NL matrix. 1 corresponds to a true connection, 0 signifies there is not a connection, -1 is periodic connection
'''
[docs]class Lattice:
"""Create a lattice instance.
Attributes of the lattice are:
xy :
BL, NL, KL, LVUC, LV, UC, Pvxydict,
LL : tuple of 2 floats
BBox :
lp : dict
polygons : list of lists of ints
list of lists of indices of each polygon
lp is a dictionary ('lattice parameters') which can contain:
LatticeTop, shape, NH, NV,
rootdir, meshfn, lattice_exten, phi, delta, theta, eta, x1, x2, x3, z, source, loadlattice_number,
check, Ndefects, Bvec, dislocation_xy, target_z, make_slit, cutz_method, cutLfrac
lattice.gxy : tuple of NX x NY arrays
"""
def __init__(self, lp={}, xy=np.array([]), NL=np.array([]), KL=np.array([]), BL=np.array([]), polygons=None):
"""Create a lattice instance."""
self.xy = xy
self.NL = NL
self.KL = KL
self.BL = BL
self.lp = lp
self.PVx = None
self.PVy = None
self.PVxydict = None
self.polygons = polygons
# Non-essential attributes
self.NLNNN = None
self.KLNNN = None
self.NLNNNangles = None
self.gxy = None
self.gr = None
def build(self):
from lepm.make_lattice import build_lattice
build_lattice(self)
def __hash__(self):
return hash((self.xy, self.NL, self.KL, self.BL, self.PVx, self.PVy))
def __eq__(self, other):
return hasattr(other, 'xy') and hasattr(other, 'NL') and hasattr(other, 'KL') and hasattr(other, 'BL') and \
self.xy == other.xy and self.NL == other.NL and self.KL == other.KL and self.BL == other.BL and \
self.PVx == other.PVx and self.PVy == other.PVy
def __ne__(self, other):
# Not strictly necessary, but to avoid having both x==y and x!=y
# True at the same time
return not(self == other)
def get_meshfn(self, attribute=True):
if 'meshfn' in self.lp:
return self.lp['meshfn']
else:
meshfn = le.find_meshfn(self.lp)
print 'output of le.find_meshfn() = ', meshfn
if attribute:
self.lp['meshfn'] = meshfn
[docs] def load(self, meshfn='auto', load_polygons=False, load_gxy=False):
"""Load a saved lattice into the lattice instance. If meshfn is specified, loads that lattice.
Otherwise, attempts to load lattice based on parameter self.lp['meshfn']. If that is also unavailable,
loads from lp[rootdir]/networks/self.LatticeTop/self.lp[lattice_exten]_NH_x_NV.
"""
print '\n\n\n\nLoading network: meshfn == ', meshfn
if meshfn == 'auto':
try:
fn = self.lp['meshfn']
print 'glob.glob(fn) = ', glob.glob(fn)
if not glob.glob(fn):
print '\n\n\n\nWARNING: Replacing rootdir path with midways in lattice params dict\n\n\n'
fn = fn.replace('/Users/npmitchell/Dropbox/Soft_Matter/GPU/', '/home/npmitchell/scratch-midway/')
self.lp['meshfn'] = fn
except KeyError:
print 'No meshfn specified for lattice.load(), attempting to find a match...'
self.get_meshfn()
fn = self.lp['meshfn']
# raise RuntimeError('There is no meshfn in the lp (params) dictionary for this lattice instance')
else:
fnglob = sorted(glob.glob(meshfn))
is_a_dir = np.where(np.array([os.path.isdir(ii) for ii in fnglob]))[0]
print 'meshfn = ', meshfn
fn = fnglob[is_a_dir[0]]
# print 'is_a_dir = ', is_a_dir
# print 'np.size(is_a_dir) = ', np.size(is_a_dir)
# print 'fn = ', fn
if np.size(is_a_dir) > 1:
print 'Found multiple lattices matching meshfn in lattice.load(). Using the first matching lattice.'
fn = fn[0]
self.lp['meshfn'] = fn
if fn[-1] == '/':
fn = fn[:-1]
print 'fn = ', fn
self.lp = le.load_params(fn+'/', 'lattice_params')
if not glob.glob(self.lp['meshfn']):
print "\n\n\n\nWARNING: Replacing rootdir path with midway's in lattice params dict:\n" +\
self.lp['meshfn'] + '\n--->\n' + fn + '\n\n'
self.lp['meshfn'] = fn
self.xy = np.loadtxt(fn + '_xy.txt', delimiter=',')
self.NL = np.loadtxt(fn + '_NL.txt', delimiter=',', dtype=int)
self.KL = np.loadtxt(fn + '_KL.txt', delimiter=',', dtype=int)
try:
self.BL = np.loadtxt(fn + '_BL.txt', delimiter=',', dtype=int)
except IOError:
self.BL = le.NL2BL(self.NL, self.KL)
name = fn.split('/')[-1]
try:
self.PVxydict = le.load_evaled_dict(fn + '/', filename=name + '_PVxydict.txt')
self.PVx = np.loadtxt(fn + '/'+name+'_PVx.txt', delimiter=',', dtype=float)
self.PVy = np.loadtxt(fn + '/' + name + '_PVy.txt', delimiter=',', dtype=float)
if self.PVx == self.PVy:
self.PVx, self.PVy = le.PVxydict2PVxPVy(self.PVxydict, self.NL)
header = 'PVx: ijth element of PVx are the x-components of the vector taking NL[i,j] to ' +\
'its image as seen by particle i'
np.savetxt(fn + '/' + name + '_PVx.txt', self.PVx, delimiter=',', header=header)
header = 'PVy: ijth element of PVy are the y-components of the vector taking NL[i,j] to ' +\
'its image as seen by particle i'
np.savetxt(fn + '/' + name + '_PVy.txt', self.PVy, delimiter=',', header=header)
except:
try:
self.PVxydict = le.load_evaled_dict(fn + '/', filename=name + '_PVxydict.txt')
self.PVx, self.PVy = le.PVxydict2PVxPVy(self.PVxydict, self.NL)
header = 'PVx: ijth element of PVx are the x-components of the vector taking NL[i,j] to ' +\
'its image as seen by particle i'
np.savetxt(fn + '/' + name + '_PVx.txt', self.PVx, delimiter=',', header=header)
header = 'PVy: ijth element of PVy are the y-components of the vector taking NL[i,j] to ' +\
'its image as seen by particle i'
np.savetxt(fn + '/' + name + '_PVy.txt', self.PVy, delimiter=',', header=header)
except IOError:
print 'No periodic vectors stored for this lattice, assuming not periodic.'
if load_polygons:
self.load_polygons()
if load_gxy:
self.load_gxy()
def load_polygons(self):
# Attempt to load polygons
if glob.glob(le.prepdir(self.lp['meshfn']) + "polygons.pkl"):
self.polygons = pickle.load(open(le.prepdir(self.lp['meshfn']) + "polygons.pkl", "r"))
else:
print 'No polygons pickle stored for this lattice. Creating it! <-- meshfn = ', self.lp['meshfn']
self.save_polygons(attribute=True)
return self.polygons
def load_gxy(self):
fn = le.prepder(self.lp['meshfn']) + "gxy.txt"
if glob.glob(fn):
gxy = np.loadtxt(fn, delimiter=',')
raise RuntimeError('Need to cast loaded array gxy as xygrid for saving function...')
lestructure.save_gxy()
[docs] def save(self, skip_polygons=False):
"""Save lattice to hard disk, outputting txt files (and images) to meshfn, which is located at lp[rootdir]/networks/LatticeTop/lattice_exten_NH_x_NV
"""
print('Saving... '+self.lp['lattice_exten'])
fbase = le.prepdir(self.lp['rootdir']) + 'networks/'+self.lp['LatticeTop']+'/'
le.ensure_dir(fbase)
if self.lp['NP_load'] == 0:
fmain = self.lp['lattice_exten']+'_'+'{0:06d}'.format(self.lp['NH'])+'_x_'+'{0:06d}'.format(self.lp['NV'])+ self.lp['slit_exten']
else:
fmain = self.lp['lattice_exten'] + '_' + 'NP{0:06d}'.format(self.lp['NP_load']) + self.lp['slit_exten']
meshfn = fbase + fmain
self.lp['meshfn'] = meshfn
fextn = '.txt'
print 'Saving ', fbase+fmain+'_xy'+fextn
np.savetxt(fbase+fmain+'_xy'+fextn, self.xy, fmt='%.18e', delimiter=',', header='x,y')
np.savetxt(fbase+fmain+'_BL'+fextn, self.BL, fmt='%i', delimiter=',', header='BL')
np.savetxt(fbase+fmain+'_NL'+fextn, self.NL, fmt='%i', delimiter=',', header='NL')
np.savetxt(fbase+fmain+'_KL'+fextn, self.KL, fmt='%i', delimiter=',', header='KL')
try:
np.savetxt(meshfn+'/'+fmain+'_LVUC'+fextn, self.LVUC, fmt='%i', delimiter=',',
header='LVUC : lattice vectors and unit cell vector identification')
except:
print 'Could not output LVUC...'
# try:
# np.savetxt(meshfn+'/'+fmain+'_LV'+fextn, self.LV, fmt='%.18e', delimiter=',', header='LV : lattice vectors')
# except:
# print 'Could not output LV...'
# try:
# np.savetxt(meshfn+'/'+fmain+'_UC'+fextn, self.UC, fmt='%.18e',
# delimiter=',', header='UC : unit cell vectors --> vectors to points in repeated macrocell')
# except:
# print 'Could not output UC...'
# try:
# np.savetxt(meshfn+'/'+fmain+'_LL'+fextn, self.lp['LL'], fmt='%.18e', delimiter=',',
# header='LL : real-space dimensions in x,y (for use in sampling k-space, for ex.)')
# except:
# print 'Could not output LL...'
# try:
# np.savetxt(meshfn+'/'+fmain+'_polygon'+fextn, self.polygon, fmt='%.18e', delimiter=',',
# header='polygon : real-space bounding polygon, for cropping lattice')
# except:
# print 'Could not output polygon...'
# Save everything else as lattice_params.txt
header = 'lattice parameters dictionary'
# self.lp['LatticeTop'] = self.LatticeTop
# self.lp['shape'] = self.shape
# self.lp['NH'] = self.NH
# self.lp['NV'] = self.NV
# self.lp['BBox'] = self.lp['BBox']
le.ensure_dir(meshfn+'/')
le.save_dict(self.lp, meshfn+'/lattice_params.txt', header, keyfmt = 'auto', valfmt = 'auto', padding_var = 7)
if 'periodicBC' in self.lp:
if self.lp['periodicBC']:
# PVxydict
header = 'PVxydict : '
filename = meshfn+'/'+fmain+'_PVxydict'+fextn
le.save_dict(self.PVxydict, filename, header )
header = 'PVx: ijth element of PVx are the x-components of the vector taking NL[i,j] to its image' +\
' as seen by particle i'
filename = meshfn+'/'+fmain+'_PVx'+fextn
np.savetxt(filename, self.PVx, delimiter=',', header=header)
header = 'PVy: ijth element of PVy are the y-components of the vector taking NL[i,j] to its image' +\
' as seen by particle i'
filename = meshfn+'/'+fmain+'_PVy'+fextn
np.savetxt(filename, self.PVy, delimiter=',', header=header)
# Save polygons
if self.lp['LatticeTop'] != 'linear' and skip_polygons is False:
self.save_polygons()
# Check lattice
print('Plotting and printing image...')
fname = fbase+fmain+'.png'
title = self.lp['lattice_exten'] + ' ' + str(self.lp['NH']) + ' x ' + str(self.lp['NV'])
xlimv = max( max(self.lp['NH'] * 0.5 + 5, self.lp['NV'] * 0.5 + 5),
max((np.max(self.xy[:, 0]) * 1.05, np.max(self.xy[:, 1]) * 1.05)))
ylimv = xlimv
climv = 0.1
# Register cmap
if not 'BlueBlackRed' in plt.colormaps():
plt.register_cmap(name='BlueBlackRed', cmap=cmaps.BlueBlackRed)
if self.BL.size > 0:
le.movie_plot_2D(self.xy, self.BL, 0*(self.BL[:,0]), fname, title, NL=self.NL, KL=self.KL,\
PVx=self.PVx, PVy=self.PVy, xlimv=xlimv, ylimv=ylimv, climv=climv,
colormap='BlueBlackRed', bgcolor='#FFFFFF' ) #,colorpoly=True)
else:
'''There are no bonds'''
le.movie_plot_2D(self.xy, self.BL, np.zeros((len(self.BL),1)), fname, title, \
PVx=self.PVx, PVy=self.PVy, xlimv=xlimv,ylimv=ylimv,climv=climv,
colormap='BlueBlackRed', bgcolor='#FFFFFF')
[docs] def get_polygons(self, attribute=False, save_if_missing=True, check=False):
"""Obtain the polygons comprising the network
Returns
-------
polygons : list of lists of ints
list of lists of indices of each polygon
"""
if self.polygons is not None:
return self.polygons
elif glob.glob(le.prepdir(self.lp['meshfn']) + "polygons.pkl"):
return self.load_polygons()
else:
polygons = le.extract_polygons_lattice(self.xy, self.BL, self.NL, self.KL, PVx=self.PVx, PVy=self.PVy,
viewmethod=False, check=check)
if attribute:
self.polygons = polygons
if save_if_missing:
self.save_polygons()
return polygons
[docs] def save_polygons(self, attribute=True, check=False):
"""Obtain the polygons comprising the network and save it as pickle"""
if self.polygons is not None:
polygons = self.polygons
else:
polygons = le.extract_polygons_lattice(self.xy, self.BL, self.NL, self.KL, PVx=self.PVx, PVy=self.PVy,
viewmethod=False, check=check)
if attribute:
self.polygons = polygons
print 'dumping polygons: ', le.prepdir(self.lp['meshfn']) + "polygons.pkl"
pickle.dump(polygons, open(le.prepdir(self.lp['meshfn']) + "polygons.pkl", "wb"))
def get_NNN_info(self, attribute=False):
if self.NLNNN is not None and self.KLNNN is not None and self.NLNNNangles is not None:
return self.NLNNN, self.KLNNN, self.NLNNNangles
else:
NLNNN, KLNNN = le.calc_NLNNN_and_KLNNN(self.xy, self.BL, self.NL, self.KL, self.PVx, self.PVy)
NLNNNangles = le.NNN_bond_angles(self.xy, self.NL, self.KL, NLNNN, KLNNN)
if attribute:
self.NLNNN = NLNNN
self.KLNNN = KLNNN
self.NLNNNangles = self.NLNNNangles
return NLNNN, KLNNN, NLNNNangles
[docs] def bond_length_histogram(self, **kwargs):
"""
Keyword arguments
-----------------
fig : figure instance or None
ax : axis instance or None
outdir : str or None
check : bool
"""
print 'creating bond length histogram (in lattice_class.py)'
lestructure.bond_length_histogram(self.xy, self.NL, self.KL, self.BL, PVx=self.PVx, PVy=self.PVy, **kwargs)
[docs] def number_variance(self, **kwargs): # outdir=None,check=False):
"""Computes number variance = (<N^2> - <N>^2) / <N> as a function of radius of circles sampling the pointset.
"""
radV, varN, a_regr = lestructure.calc_number_variance(self.xy, self.lp['LL'], **kwargs)
return np.dstack((radV, varN))[0], a_regr
[docs] def calc_gxy_gr(self, **kwargs):
"""
Parameters
----------
**kwargs : outdir=None, dr=0.1, eps=1e-7, maxgr=0.04, check=False
Returns
----------
gxy_info : xc, yc, gxy
gr_info : rc, grV
"""
xgrid, ygrid, gxy_grid, rcenters, grV = lestructure.calc_gxy_gr(self.xy, self.lp['BBox'], **kwargs)
return (xgrid, ygrid, gxy_grid), np.dstack((rcenters, grV))[0]
def calc_fancy_gr(self, **kwargs):
return lestructure.calc_fancy_gr(self.xy, self.lp['BBox'], **kwargs)
def pointset_fft(self, **kwargs):
return lestructure.pointset_fft(self.xy, **kwargs)
def pointset(self, **kwargs):
return lestructure.pointset(self.xy, **kwargs)
def boundary(self):
return le.extract_boundary(self.xy, self.NL, self.KL, self.BL, check=self.lp['check'])
[docs] def structure_factor(self, **kwargs):
"""
Returns
--------
Skxy : tuple of 3 float arrays
Skr : (NX*NY) x 2 float array
"""
kx, ky, Skmesh, kr, Skr = lestructure.calc_structure_factor(self.xy, self.lp['LL'], **kwargs)
return (kx, ky, Skmesh), np.dstack((kr, Skr))[0]
def plot_NNNangle_hist(self, outdir=None, title=r'$2 \theta_{nml}$, next-nearest neighbor bond angle distribution',
fname='figNNNangles', FSFS=12, show=True):
""""""
if outdir is not None and outdir != 'none':
le.ensure_dir(le.prepdir(outdir))
NLNNN, KLNNN, NLNNNangles = self.get_NNN_info()
# print 'NLNNNangles = ', NLNNNangles
thetaH = np.mod(2.*NLNNNangles, 2.*np.pi)
# print 'thetaH= ', thetaH
thetaH[thetaH > np.pi] = - (2. * np.pi - thetaH[thetaH > np.pi])
# print 'thetaH= ', thetaH
fig, ax = plt.subplots(1, 2, figsize=(15, 7))
gs = gridspec.GridSpec(7, 7)
ax[0].set_position(gs[1:6, 0:2].get_position(fig))
ax[1].set_position(gs[1:6, 3:].get_position(fig))
le.movie_plot_2D(self.xy, self.BL, 0*(self.BL[:, 0]), None, '', NL=self.NL, KL=self.KL,
PVx=self.PVx, PVy=self.PVy, axcb=None,
colorz=False, colormap='BlueBlackRed', bgcolor='#FFFFFF', axis_off=True, ax=ax[0])
ax[1].hist(thetaH[np.where(KLNNN > 0)].ravel() / np.pi, bins=100)
ax[1].set_title(title, fontsize=FSFS)
ax[1].set_ylabel('Occurence', fontsize=FSFS)
ax[1].set_xlabel(r'$2 \theta_{nml}/\pi$', fontsize=FSFS)
ax[1].set_xlim([-1, 1])
if outdir is not None:
# pickle.dump(fig, file(outdir+fname+'.pkl','w'))
plt.savefig(outdir+fname+'.png')
if show:
print 'Displaying figure...'
plt.show()
else:
return fig, ax
[docs] def plot_BW_lat(self, fig=None, ax=None, meshfn='./', exten='.pdf', save=True, close=True, axis_off=True,
title='auto', **kwargs):
"""Plot (and save if desired) a black and white image of the lattice with connectivity.
If save is False, displays the lattice if show is True. Otherwise, adds lattice to axis.
Parameters
----------
kwargs : keyword arguments for le.movie_plot_2D()
"""
# Register cmap
if 'BlueBlackRed' not in plt.colormaps():
plt.register_cmap(name='BlueBlackRed', cmap=cmaps.BlueBlackRed)
print('Plotting as black and white...')
if save:
if meshfn != './' and meshfn != 'none':
le.ensure_dir(le.prepdir(meshfn))
fname = meshfn + '/' + meshfn.split('/')[-1] + '_BW' + exten
else:
fname = 'none'
if title == 'auto':
title = self.lp['lattice_exten'] + ' ' + str(self.lp['NH']) + ' x ' + str(self.lp['NV'])
# xlimv = max(max(self.lp['NH'] * 0.5 + 5, self.lp['NV'] * 0.5 + 5),
# xlimv = max((np.max(self.xy[:, 0]) * 1.05, np.max(self.xy[:, 1]) * 1.05))
# ylimv = xlimv
climv = 0.1
if self.BL.size > 0:
# if save:
[ax, axcb] = le.movie_plot_2D(self.xy, self.BL, 0*(self.BL[:, 0]), fname, title, fig=fig, ax=ax, NL=self.NL,
KL=self.KL, PVx=self.PVx, PVy=self.PVy, climv=climv,
axcb=None, colorz=False, colormap='BlueBlackRed', bgcolor='#FFFFFF',
axis_off=axis_off, **kwargs)
# else:
# le.display_lattice_2D(self.xy, self.BL, title=title, NL=self.NL, KL=self.KL,
# PVx=self.PVx, PVy=self.PVy, xlimv=xlimv, ylimv=ylimv, climv=climv,
# colorz=False, colormap='BlueBlackRed', bgcolor='#FFFFFF', axis_off=axis_off,
# close=close)
else:
raise RuntimeWarning('Could not save BW plot since BL is empty (no bonds)!')
return [ax, axcb]
def get_gxy_gr(self, attribute=False):
if self.gxy is not None and self.gr is not None:
gxy = self.gxy
gr = self.gr
else:
gxy, gr = self.calc_gxy_gr()
if attribute:
self.gxy = gxy
self.gr = gr
return gxy, gr
def save_gxy_gr(self, outdir=None):
self.get_gxy_gr()
raise RuntimeError('Lattice.save_gxy_gr() is not written!')
[docs] def plot_numbered(self, **kwargs):
"""
Parameters
----------
**kwargs : keyword arguments for self.plot_BW_lat()
"""
[ax, axcb] = self.plot_BW_lat(meshfn=None, save=False, close=True, **kwargs)
for ii in range(len(self.xy)):
ax.text(self.xy[ii, 0] + 0.3, self.xy[ii, 1], str(ii))
return [ax, axcb]
[docs] def summarize_structure(self):
"""If customization in the plots characterizing structure are desired,
first use methods of the lattice passed to this method. As an example, to use
custom arguments in the gxy_gr calculation, run:
lattice.calc_gxy_gr( custom arguments here... )
lattice.summarize_structure()
"""
plt.close('all')
cmaps.ensure_cmaps()
if 'meshfn' in self.lp:
outdir = le.prepdir(self.lp['meshfn'])
le.ensure_dir(outdir)
fname = outdir + 'structure_summary.png' #self.lp['meshfn'].split('/')[-1] + '.pdf'
else:
fname = './structure_summary.png'
# Set up figure
Wfig = 180
x0s = 10
y0s = 10
ws = Wfig * 0.22
hs = 3./4.*ws
# space between figures
vspace = 10
# space above top figure
tspace = 10
wB = ws*0.8
hB = hs*0.8
wss = wB*.3
hss = wss*3./4.
# label space
lbs = ws*0.15
Hfig = y0s + ws*2 + vspace + tspace
# create figure
fig = sps.figure_in_mm(Wfig, Hfig)
label_params = dict(size=12, fontweight='normal')
# a: scatter, b: bond hist, c: gxy, d: gr, dINSET: grinset e: variance, f: Skxy, g: Skr
axes = [
sps.axes_in_mm(x0, y0, width, height, label=part, label_params=label_params)
for x0, y0, width, height, part in (
[x0s , y0s+hs+vspace , wB , wB , 'a'], # snippet
[x0s+ws , y0s+hs+vspace+lbs , wB , hB , 'b'], # bond histogram
[x0s+ws+ws , y0s+hs+vspace+lbs , hB , hB , 'c'], # g(x,y)
[x0s+ws+ws+ws , y0s+hs+vspace+lbs , wB , hB , 'd'], # g(r)
[x0s+ws+ws+ws+wss*2.5 , y0s+hs+vspace+hss*3.0, wss , hss, '' ], # g(r) inset
[x0s , y0s , wB , hB , 'e'], # nvariance
[x0s+ws+ws , y0s , hB , hB , 'f'], # S(kx,ky)
[x0s+ws+ws+ws , y0s , wB , hB , 'g'], # S(k)
)
]
try:
gxyGrid = self.gxy
grVec = self.gr
except:
gxyGrid, grVec = self.calc_gxy_gr() #outdir=None, dr=0.1, eps=1e-7, check=False)
try:
(varNvec, a_regr) = self.nvariance
except:
(varNvec, a_regr) = self.number_variance() #outdir=None, check=False)
try:
SkxyGrid, SkrVec = self.Skxy, self.Skr
except:
SkxyGrid, SkrVec = self.structure_factor()
# Plot snippet of lattice
keep = np.where(np.logical_and(self.xy[:,0]<15, self.xy[:,1]<15))[0]
xysnip, NLsnip, KLsnip, BLsnip = le.remove_pts(keep, self.xy, self.BL)
xlimv = 10
ylimv = xlimv
if BLsnip.size >0:
'''Plot snippet of lattice'''
if len(keep) == len(self.xy):
le.movie_plot_2D(self.xy, self.BL, 0 * (self.BL[:, 0]), fname='none', title='', NL=self.NL, KL=self.KL, \
PVx=self.PVx, PVy=self.PVy, ax=axes[0], fig=fig, axcb='none', xlimv=xlimv, ylimv=ylimv, \
colorz=False, colormap='BlueBlackRed', bgcolor='#FFFFFF', axis_off=True)
else:
le.movie_plot_2D(xysnip, np.abs(BLsnip), 0 * (BLsnip[:, 0]), fname='none', title='', NL=NLsnip,
KL=KLsnip, \
PVx=None, PVy=None, ax=axes[0], fig=fig, axcb='none', xlimv=xlimv, ylimv=ylimv, \
colorz=False, colormap='BlueBlackRed', bgcolor='#FFFFFF', axis_off=True)
ind = 1
# Histogram bond list
self.bond_length_histogram(fig=fig, ax=axes[ind], outdir=None, check=False, savetxt=False)
ind += 1
# Plot gxy as heatmap
# axes[ind].pcolor( gxyGrid[0], gxyGrid[1], gxyGrid[2], cmap='viridis', vmin=0.0) #vmax=4.0)
imgxy = axes[ind].pcolormesh( gxyGrid[0], gxyGrid[1], gxyGrid[2], cmap='viridis', vmin=0.0)
divider = make_axes_locatable(axes[ind])
# Append axes to the right of ax3, with 20% width of ax3
cdiv = divider.append_axes("right", size="05%", pad=0.05)
# Create colorbar in the appended axes
# Tick locations can be set with the kwarg `ticks` and the format of the ticklabels with kwarg `format`
cb = plt.colorbar(imgxy, cax=cdiv, format="%.2f")
# cb.set_label(r'$ g(x,y)$')
axes[ind].axis('image')
titlestr = r'$g(\mathbf{x})$' # with system size=({0:0.3f}'.format(2.* np.min(np.abs(self.lp['BBox'])))+')'
# print 'title = ', titlestr
axes[ind].set_title(titlestr)
axes[ind].set_xlim(-5, 5)
axes[ind].set_ylim(-5, 5)
ind += 1
# Plot g(r)
axes[ind].plot(grVec[:, 0], grVec[:, 1], 'r.-')
axes[ind].set_xlim(0, 5)
axes[ind].set_xlabel(r'$r$')
axes[ind].set_ylabel(r'$g(r)$')
ind += 1
# Plot g(r) inset
axes[ind].plot(grVec[:, 0], grVec[:, 1], 'r-')
ind += 1
# Plot number variance
axes[ind].plot(varNvec[:, 0], varNvec[:, 1], 'r.-')
axes[ind].set_xlabel(r'$R$')
axes[ind].set_ylabel(r'$\sigma^2(R)$')
ind += 1
# Plot S(kx,ky) as heatmap
imgSk = axes[ind].pcolor( SkxyGrid[0], SkxyGrid[1], SkxyGrid[2], cmap='viridis', vmin=0.0)
divider = make_axes_locatable(axes[ind])
cdiv = divider.append_axes("right", size="05%", pad=0.05)
# Create colorbar in the appended axes
# Tick locations can be set with the kwarg `ticks` and the format of the ticklabels with kwarg `format`
cb = plt.colorbar(imgSk, cax=cdiv, format="%.0f")
# cb.set_label(r'$ g(x,y)$')
axes[ind].axis('image')
titlestr = r'$S(\mathbf{k})$' # with system size='+'({0:0.3f}'.format(2.* np.min(np.abs(self.lp['BBox'])))+')'
# print 'title = ', titlestr
axes[ind].set_title(titlestr)
ind += 1
# Plot S(k)
# print 'SkrVec = ', SkrVec
SkrVec = SkrVec[~np.isnan(SkrVec[:,1]),:]
axes[ind].plot(SkrVec[:,0], SkrVec[:,1], 'r.-')
axes[ind].axis('tight')
axes[ind].set_xlim(0, np.max(SkrVec[:,0]))
axes[ind].set_ylim(0, 5)
axes[ind].set_xlabel(r'$k$')
axes[ind].set_ylabel(r'$S(k)$')
plt.savefig(fname)
if self.lp['check']:
plt.show()
if __name__ == '__main__':
'''Perform an example of using the lattice class'''
# check input arguments for timestamp (name of simulation is timestamp)
parser = argparse.ArgumentParser(description='Specify time string (timestr) for gyro simulation.')
# Task
parser.add_argument('-NNNanglehist', '--NNNanglehist', help='Make a sequence of NNN angle histograms',
action='store_true')
parser.add_argument('-NNNvary_param', '--NNNvary_param', help='Parameter to vary in NNN angle visualization',
type=str, default='delta')
parser.add_argument('-loadlat', '--loadlat', help='Demonstrate loading a lattice from meshfn string',
action='store_true')
parser.add_argument('-view_lattice', '--view_lattice', help='Preview a saved lattice by loading it',
action='store_true')
parser.add_argument('-view_numbered', '--view_numbered', help='Preview a saved lattice and label particles',
action='store_true')
parser.add_argument('-summarize', '--summarize', help='Summarize the structure of a lattice',
action='store_true')
parser.add_argument('-redo_gxy', '--redo_gxy', help='Redo calculation of g(x,y)',
action='store_true')
parser.add_argument('-maxgxy', '--maxgxy', help='Max plotted value for g(x,y)', type=float, default=0.04)
parser.add_argument('-load_save_gxy', '--load_save_gxy', help='Load and plot g(x,y)', action='store_true')
parser.add_argument('-pointset', '--pointset', help='Plot the lattice point set', action='store_true')
# Lattice Geometry Parameters
parser.add_argument('-N', '--N', help='Mesh width AND height, in number of lattice spacings (leave blank to specify separate dims)', type=int, default=-1)
parser.add_argument('-NH', '--NH', help='Mesh width, in number of lattice spacings', type=int, default=20)
parser.add_argument('-NV', '--NV', help='Mesh height, in number of lattice spacings', type=int, default=20)
parser.add_argument('-NP', '--NP_load', help='Number of particles in mesh, overwrites N, NH, and NV.',
type=int, default=0)
parser.add_argument('-LT', '--LatticeTop', help='Lattice topology: linear, hexagonal, triangular, deformed_kagome, hyperuniform, circlebonds, penroserhombTri',
type=str, default='hexagonal')
parser.add_argument('-shape', '--shape', help='Shape of the overall mesh geometry', type=str, default='square')
parser.add_argument('-Nhist', '--Nhist', help='Number of bins for approximating S(k)', type=int, default=50)
parser.add_argument('-kr_max', '--kr_max', help='Upper bound for magnitude of k in calculating S(k)', type=int, default=30)
parser.add_argument('-check', '--check', help='Check outputs during computation of lattice', action='store_true')
parser.add_argument('-nice_plot', '--nice_plot', help='Output nice pdf plots of lattice', action='store_true')
# For loading and coordination
parser.add_argument('-LLID', '--loadlattice_number',
help='If LT=hyperuniform/isostatic, selects which lattice to use', type=str, default='01')
parser.add_argument('-LLz', '--loadlattice_z', help='If LT=hyperuniform/isostatic, selects what z index to use',
type=str, default='001')
parser.add_argument('-source', '--source',
help='Selects who made the lattice to load, if loaded from source (ulrich, hexner, etc)',
type=str, default='hexner')
parser.add_argument('-cut_z', '--cut_z',
help='Declare whether or not to cut bonds to obtain target coordination number z', type=bool,
default=False)
parser.add_argument('-cutz_method', '--cutz_method',
help='Method for cutting z from initial loaded-lattice value to target_z (highest or random)',
type=str, default='none')
parser.add_argument('-z', '--target_z', help='Coordination number to enforce', type=float, default=-1)
parser.add_argument('-Ndefects', '--Ndefects', help='Number of defects to introduce', type=int, default=1)
parser.add_argument('-Bvec', '--Bvec', help='Direction of burgers vectors of dislocations (random, W, SW, etc)',
type=str, default='W')
parser.add_argument('-dislocxy', '--dislocation_xy',
help='Position of single dislocation, if not centered at (0,0), as strings sep by / (ex: 1/4.4)',
type=str, default='none')
parser.add_argument('-thres', '--thres', help='Threshold distance (from letters in uofc, for ex)',
type=float, default=1.0)
# Global geometric params
parser.add_argument('-periodic', '--periodicBC', help='Enforce periodic boundary conditions', action='store_true')
parser.add_argument('-slit', '--make_slit', help='Make a slit in the mesh', action='store_true')
parser.add_argument('-phi', '--phi', help='Shear angle for hexagonal (honeycomb) lattice in radians/pi',
type=float, default=0.0)
parser.add_argument('-delta', '--delta', help='Deformation angle for hexagonal (honeycomb) lattice in radians/pi',
type=float, default=120./180.)
parser.add_argument('-eta', '--eta', help='Randomization/jitter in lattice', type=float, default=0.000)
parser.add_argument('-theta', '--theta', help='Overall rotation of lattice', type=float, default=0.000)
parser.add_argument('-alph', '--alph', help='Twist angle for twisted_kagome, max is pi/3', type=float, default=0.000)
parser.add_argument('-huno', '--hyperuniform_number', help='Hyperuniform realization number', type=str, default='01')
parser.add_argument('-skip_gr', '--skip_gr', help='Skip calculation of g(r) correlation function for the lattice',
action='store_true')
parser.add_argument('-skip_gxy', '--skip_gxy',
help='Skip calculation of g(x,y) 2D correlation function for the lattice', action='store_true')
parser.add_argument('-skip_sigN', '--skip_sigN', help='Skip calculation of variance_N(R)', action='store_true')
parser.add_argument('-fancy_gr', '--fancy_gr',
help='Perform careful calculation of g(r) correlation function for the ENTIRE lattice',
action='store_true')
args = parser.parse_args()
if args.N >0:
NH = args.N
NV = args.N
else:
NH = args.NH
NV = args.NV
phi = np.pi* args.phi
delta = np.pi* args.delta
strain =0.00 #initial
# z = 4.0 #target z
if args.LatticeTop == 'linear':
shape = 'line'
else:
shape = args.shape
theta = args.theta
eta = args.eta
transpose_lattice=0
make_slit = args.make_slit
# deformed kagome params
x1 = 0.0
x2 = 0.0
x3 = 0.0
z = 0.0
rootdir = '/Users/npmitchell/Dropbox/Soft_Matter/GPU/'
lp = {'LatticeTop': args.LatticeTop,
'shape': shape,
'NH': NH,
'NV': NV,
'NP_load': args.NP_load,
'rootdir': '/Users/npmitchell/Dropbox/Soft_Matter/GPU/',
'phi': phi,
'delta': delta,
'theta': theta,
'eta': eta,
'x1': x1,
'x2': x2,
'x3': x3,
'z': z,
'source': args.source,
'loadlattice_number': args.loadlattice_number,
'check': args.check,
'Ndefects': args.Ndefects,
'Bvec': args.Bvec,
'dislocation_xy': args.dislocation_xy,
'target_z': args.target_z,
'make_slit': args.make_slit,
'cutz_method': args.cutz_method,
'cutLfrac': 0.0,
'conf': 01,
'periodicBC': args.periodicBC,
'loadlattice_z': args.loadlattice_z,
'alph': args.alph,
'origin': np.array([0., 0.]),
'thres': args.thres,
}
if args.NNNanglehist:
cmaps.register_colormaps()
outroot = '/Users/npmitchell/Dropbox/Soft_Matter/GPU/experiments/NNNangles/'
outdir = outroot + lp['LatticeTop'] + '/'
if lp['LatticeTop'] == 'hexagonal':
le.ensure_dir(outdir)
if args.NNNvary_param == 'delta':
for delta in np.arange(0.66666, 1.23, 0.02):
lp['delta'] = delta * np.pi
lat = Lattice(lp)
lat.build()
print 'lat.xy = ', lat.xy
# lat.plot_BW_lat(meshfn= outroot + lp['LatticeTop'], exten='_delta'+le.float2pstr(delta)+'pi.png')
fname = 'N'+'{0:03d}'.format(max(lp['NH'], lp['NP_load'])) + \
'NNNangles_phi'+le.float2pstr(lp['phi']/np.pi)+'pi'+le.float2pstr(delta/np.pi)+'pi.png'
lat.plot_NNNangle_hist(outdir=outdir, fname=fname, show=False)
elif args.NNNvary_param == 'phi':
for phi in np.arange(0.0, 0.5, 0.0075):
lp['phi'] = phi*np.pi
lat = Lattice(lp)
lat.build()
print 'lat.xy = ', lat.xy
# lat.plot_BW_lat(meshfn= outroot + lp['LatticeTop'], exten='_delta'+le.float2pstr(delta)+'pi.png')
fname = 'N'+'{0:03d}'.format(max(lp['NH'], lp['NP_load'])) + \
'NNNangles_phi'+le.float2pstr(phi)+'pi'+'_delta'+le.float2pstr(lp['delta']/np.pi)+'pi.png'
lat.plot_NNNangle_hist(outdir=outdir, fname=fname, show=False)
else:
le.ensure_dir(outdir)
lat = Lattice(lp)
lat.build()
print 'lat.xy = ', lat.xy
lat.plot_BW_lat(meshfn=outroot + lp['LatticeTop'], exten='.pdf')
fname = 'N'+'{0:03d}'.format(max(lp['NH'], lp['NP_load'])) + 'NNNangles_' +lp['LatticeTop']+'.png'
lat.plot_NNNangle_hist(outdir=outdir, fname=fname, show=False)
if args.loadlat:
# Check loading
meshfn = '/Users/npmitchell/Dropbox/Soft_Matter/GPU/networks/hucentroid/hucentroid_square_periodic_d01_NP000020'
lat = Lattice()
lat.load(meshfn=meshfn)
lat.PVx, lat.PVy = le.PVxydict2PVxPVy(lat.PVxydict, lat.NL)
BM = le.NL2BM(lat.xy, lat.NL, lat.KL, PVx=lat.PVx, PVy=lat.PVy)
bL = le.BM2bL(lat.NL, BM, lat.BL)
plt.hist(bL)
plt.show()
if args.view_lattice:
# pointset
lat = Lattice(lp)
lat.load()
lat.plot_BW_lat(fig=None, ax=None, meshfn=None, exten='.pdf', save=False, close=True, axis_off=False,
title='auto')
plt.show()
if args.view_numbered:
lat = Lattice(lp)
lat.load()
lat.plot_numbered()
plt.show()
if args.summarize:
lat.summarize_structure()
if args.redo_gxy:
lat = Lattice(lp)
lat.load(load_polygons=False)
lat.calc_gxy_gr(outdir=le.prepdir(lp['meshfn']), dr=0.1, eps=1e-7, check=False, maxgxy=args.maxgxy)
if args.load_save_gxy:
lat = Lattice(lp)
lat.load(load_polygons=False)
lat.save_gxy_gr()
if args.nice_plot:
meshfn = le.find_meshfn(lp)
lp['meshfn'] = meshfn
lat = Lattice(lp)
lat.load()
print 'Saving nice BW plot...'
lat.plot_BW_lat(meshfn=lat.lp['meshfn'])
if args.pointset:
meshfn = le.find_meshfn(lp)
lp['meshfn'] = meshfn
lat = Lattice(lp)
lat.load()
print '\nPerforming image of points...'
lat.pointset(outdir=lat.lp['meshfn'], wsfrac=0.8)
