N2O AF calcs v1¶
22/03/24
Alignment frame calcs preliminary:
- Load and plot alignment data (ADMs) for all temperatures.
- Compute AF results for given temp (with t downsampling), note default case assumes linear pol, $z$-axis alignment.
- Plot XS and betas, interactive heatmaps and line plots.
Note all plots are interactive, mouse-over for data points, and use toolbars for zoom/selection, and widgets for browsing data types (where applicable).
For general methods: https://epsproc.readthedocs.io/en/dev/demos/ePSproc_class_demo_161020.html
In [18]:
# Quick hack to override default HTML template
# NOT required in some JLab versions.
# https://stackoverflow.com/a/63777508
# https://stackoverflow.com/questions/21971449/how-do-i-increase-the-cell-width-of-the-jupyter-ipython-notebook-in-my-browser
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:100% !important; }</style>"))
Load alignment data¶
Now with ADM class.
In [1]:
from epsproc.classes.alignment import ADM
from pathlib import Path
dataPath = Path('~/ePS/N2O/alignment/Wigner D polyatomic normalized')
ADMtype = '.txt'
ADMin = ADM(fileBase=dataPath.expanduser(), ext = ADMtype)
ADMin.loadData(keyType='setADMsT', normType='wignerDpoly')
In [2]:
# import epsproc as ep
# ep.setADMs
In [3]:
ADMin.plot(width=1000)
# ADMin.plot(xlim=(30,50), width=800) # Auto-stack plus pass args & zoom in on specific region (note slider will reset region with interactive zoom)
Load matrix elements¶
In [4]:
# Plotters
# from epsproc.plot import hvPlotters
# Multijob class dev code
from epsproc.classes.multiJob import ePSmultiJob
import warnings
# warnings.filterwarnings('once') # Skip repeated numpy deprecation warnings in current build (xr15 env)
warnings.filterwarnings('ignore') # Skip repeated numpy deprecation warnings in current build (xr15 env)
In [5]:
# # Scan for subdirs, based on existing routine in getFiles()
# fileBase = Path('/home/paul/ePS/OCS/OCS_survey') # Data dir on Stimpy
fileBase = Path('~/fock-mount/globalhome/eps/N2O/N2O_valence').expanduser() # Data dir on Jake
In [6]:
# TODO: fix orb label here, currently relies on (different) fixed format
data = ePSmultiJob(fileBase, verbose = 0)
data.scanFiles()
# data.jobsSummary()
Compute AF¶
Note default case assumes linear pol, $z$-axis alignment.
In [7]:
# Set ADMs to use for calc - set for specific T and downsample
# Subselection, note temperature is set by dataKey here
temp = '10K'
ADMin.subsetADMs(dataKey = temp, trange=[35,45],tStep = 5)
# Plot subselection
ADMin.plot(keys='ADM')
# Set to main data structure for calcs
data.data['ADM'] = ADMin.data['ADM']
In [8]:
# Basic routine (no grouping)
thres = 1e-5
# data.AFBLM(AKQS = data.data['ADM']['ADM'].real, selDims = {'Type':'L'}, thres=1e-2) # OK for 41 t-points (x2 downsample), but missing some orbs and t-points in output!
# data.AFBLM(AKQS = data.data['ADM']['ADM'].real, selDims = {'Type':'L'}, thres=thres) # OK for 41 t-points (x2 downsample), still some missing regions for thres=1e-2
# mostly contiguous for thres = 1e-3,1e-4
data.AFBLM(AKQS = data.data['ADM']['ADM'].real, selDims = {'Type':'L'}, thres=thres) # outputs fairly smooth for thres=1e-5
# Save full class with Pickle
# TODO: GENERAL DEBUG FOR ACCIDENTAL FUNCTION PASSING HERE!!!!
import pickle
outfile = f'N2O_AFBLM_{temp}_{thres}_220323.pickle'
with open(outfile, 'wb') as handle:
# data.data['ADM']['plots']['ADM']['pData'].attrs['pTypeDetails'] = [] # UGH, this breaks pickle unless removed.
pickle.dump(data.data, handle, protocol=pickle.HIGHEST_PROTOCOL)
print(f'Saved data.data with pickle, file: {outfile}')
Plots¶
NOTE: interactive plots. Roll over plots for values and tool options. For heatmaps use the "box select" tool on the histogram at the side to change the colour mapping scale.
Remarks:
- Some unphysical spikes in values <5eV, these are regions with small cross-sections, and are currently filtered out in these plots.
- General patterns seem to be simlar to R-matrix results, as expected, although absolute Blm values seem larger...? Might be a normalisation difference here, TBC.
Cross-sections¶
In [9]:
from epsproc.plot import hvPlotters
import numpy as np
# Set data - use existing plot routine but skip line outs (quite slow, replotted at end of notebook)
data.BLMplot(backend='hv', xDim = 't', renderPlot=False)
# Heatmaps for cross-sections
xrAF = data.plots['BLMplot']['XR'].copy()
Etype = 'Eke'
# Set opts to match sizes - should be able to link plots or set gridspace to handle this?
hvPlotters.setPlotDefaults(fSize = [750,300], imgSize = 600)
cmap = 'coolwarm' # See http://holoviews.org/user_guide/Colormaps.html
# Get cross-sections
XS = np.abs(xrAF.XSrescaled) # Should be in MB
XS.name = 'XS/MB'
hvXS = hvPlotters.hv.Dataset(XS.fillna(0))
# Plot heatmap + ADMs
(hvXS.to(hvPlotters.hv.HeatMap, kdims=['t',Etype]).opts(cmap=cmap).hist() + data.data['ADM']['ADM'].unstack().squeeze().real.hvplot.line(x='t').overlay('K')).cols(1)
Out[9]:
AF-$\beta_{L,M}$¶
In [10]:
# Heatmap for l=2,4,6
# NOTE - clipped cmap to physical values with clim=(-1,2), but some spikes.
# BLM = np.abs(xrAF.XSrescaled) # Should be in MB
# XS.name = 'XS/MB'
# Filter out small XS regions?
xrAFfiltered = xrAF.where(np.abs(xrAF.XSrescaled) > 0.01)
hvDS = hvPlotters.hv.Dataset(xrAFfiltered.unstack('BLM').fillna(0).squeeze(drop=True))
# hvDS = hvPlotters.hv.Dataset(xrAF.unstack('BLM').fillna(0).squeeze(drop=True))
# Plot heatmap + ADMs
(hvDS.select(l=[2,4,6], m=0).to(hvPlotters.hv.HeatMap, kdims=['t',Etype]).opts(cmap=cmap, clim=(-0.5,1)).hist() + data.data['ADM']['ADM'].unstack().squeeze().real.hvplot.line(x='t').overlay('K')).cols(1) #.overlay(['l','t']).layout('Orb').cols(1)
Out[10]:
In [11]:
# Line-outs
# Note this is quite slow to render
data.BLMplot(backend='hv', xDim = 't', ylim=(-1,2), filterXS=0.01)
Out[11]:
Versions¶
Code¶
In [12]:
# print(data.jobInfo['ePolyScat'][0])
# print('Run: ' + jobInfo['Starting'][0].split('at')[1])
In [13]:
import scooby
scooby.Report(additional=['epsproc', 'xarray', 'jupyter', 'holoviews', 'pandas'])
Out[13]:
In [14]:
# Check current Git commit for local ePSproc version
!git -C {Path(ep.__file__).parent} branch
!git -C {Path(ep.__file__).parent} log --format="%H" -n 1
In [15]:
# Check current remote commits
!git ls-remote --heads https://github.com/phockett/ePSproc
Notes¶
22/03/24
v1, quick run based on previous OCS code plus a few recent function updates.
TODO:
- Testing with alternative polarization states.
- More careful tests/checks in low-XS regions (currently just filter out, but see BLM spikes in a few cases without XS filter applied).
- More careful comparison with R-matrix results.
- AF-PAD plots, grids should work here...
Scratch¶
In [16]:
break
In [ ]:
# Compute AFBLMs for aligned case
keys = data._keysCheck(None)
keys = list(set(keys) - set(['ADM']))
#*************** Throw everything in method - breaks for large datasets (RAM issues)
# data.AFBLM(keys = keys, AKQS = ADMs, selDims = {'Type':'L'}) # NOTE THIS USES independently set ADMs!
# data.AFBLM(keys = keys, AKQS = data.data['ADM']['ADM'].real, selDims = {'Type':'L'}, thres=1e-3) # OK for 61 t-points (x2 downsample)
# data.AFBLM(keys = keys, AKQS = data.data['ADM']['ADM'].real, selDims = {'Type':'L'}, thres=1e-3, thresDims = ['Eke','t']) # NEED TO ADD DIM CHECKING TO thresDims
# data.AFBLM(keys = keys, AKQS = data.data['ADM']['ADM'].real, selDims = {'Type':'L'}, thres=1e-4) # 1e-4 working OK for 30 t-points (x4 downsample)
# data.AFBLM(keys = keys, AKQS = data.data['ADM']['ADM'].real, selDims = {'Type':'L'}, thres=1e-5) # 1e-5 working OK for 16 t-points (x8 downsample)
# data.AFBLM(keys = keys, AKQS = data.data['ADM']['ADM'].real, selDims = {'Type':'L'}, thres=1e-6) # 1e-6 working OK for 16 t-points (x8 downsample)
# data.AFBLM(keys = keys, AKQS = data.data['ADM']['ADM'].real, selDims = {'Type':'L'}, thres=None) # FAILS for 16 t-points (x8 downsample) at 60Gb RAM on Jake (Swap file size issue????)
# OK for 8 t-points (x16 downsample) at ~55Gb RAM peak on Jake.
# thres=1e-4) BREAKS KERNEL #, thresDims=['Eke','t']) # TODO: currently thersDims must be in matE only for AFBLM routine.
# 21/07/22 version - above only gives L=2 max, missing L=4,6. But bumping threshold to 1e-4 kills calculation...
# data.AFBLM(keys = keys, selDims = {'Type':'L'}) # Default case... sets/uses ADMs per key
# UPDATE: testing with scale factor for ADMs, and forcing to real - might be better?
#*********** Using grouping to avoid memory issues
# QUITE slow with current code, but more robust
# TODO: make this better with basis set recycle (cf. fitting code)
#
import pickle
def calcAFBLMGroup(data,ADMs,keys,AFBLMdict):
"""Wrap for XR groupby"""
# Class routine - returns to main datastruture
data.AFBLM(keys = keys, AKQS = ADMs, selDims = {'Type':'L'}, thres=None) #.map(standardize)
# Copy to new structure & restack
# keys = data._keysCheck(keys)
# AFBLMdict = {}
# Restack data to lists for merge
for key in keys:
AFBLMdict[key]['AFBLMlist'].append(data.data[key]['AFBLM'].copy())
return AFBLMdict
AFBLMdict = {k:{'AFBLMlist':[]} for k in keys}
# tBins = np.arange(38,39.1,0.5) # Test range
tBins = np.arange(trange[0], trange[1] + 0.1, 0.5)
# trange
import time
for item in data.data['ADM']['ADM'].groupby_bins("t", tBins):
start = time.time()
print(f'Calculating for group: {item[0]}, {item[1].t.size} t-points.')
AFBLMdict = calcAFBLMGroup(data, item[1], keys, AFBLMdict)
end = time.time()
print(f'Delta: {end - start}s')
# Concat XR lists
for key in keys:
AFBLMdict[key]['full'] = xr.concat(AFBLMdict[key]['AFBLMlist'],"t")
AFBLMdict[key]['full'].name = f'{key}_AFBLM'
# Push back to main class
data.data[key]['AFBLM'] = AFBLMdict[key]['full'].copy()
# Save object with Pickle
with open(f'OCS_{key}_AFBLM_220722.pickle', 'wb') as handle:
print(f'Saving {key} with pickle')
pickle.dump(AFBLMdict[key]['full'], handle, protocol=pickle.HIGHEST_PROTOCOL)
# Save full class with Pickle
# NOW GIVING ERRORS ON RERUN (although worked previously), AttributeError: Can't pickle local object 'plotTypeSelector.<locals>.<lambda>'
# NOT SURE WHERE THIS IS COMING FROM AS YET
# UPDATE: IN data.data['ADM']['plots']['ADM']['pData'].attrs['pTypeDetails']
# TODO: GENERAL DEBUG FOR ACCIDENTAL FUNCTION PASSING HERE!!!!
# import pickle
with open(f'OCS_AFBLM_220722.pickle', 'wb') as handle:
data.data['ADM']['plots']['ADM']['pData'].attrs['pTypeDetails'] = [] # UGH, this breaks pickle unless removed.
pickle.dump(data, handle, protocol=pickle.HIGHEST_PROTOCOL)
print(f'Saved data (full class) with pickle')