Method development for geometric functions

26/02/20

Aims:

  • Develop \(\beta_{L,M}\) formalism.
  • Develop corresponding numerical methods.
  • Speed things up (see low-level benchmarking notebook).
  • Analyse geometric terms.

Setup

[50]:

# Imports
import numpy as np
import pandas as pd
import xarray as xr
from functools import lru_cache  # For function result caching

# Special functions
# from scipy.special import sph_harm
import spherical_functions as sf
import quaternion

# Performance & benchmarking libraries
# from joblib import Memory
# import xyzpy as xyz
import numba as nb

# Timings with ttictoc
# https://github.com/hector-sab/ttictoc
# from ttictoc import TicToc

# Package fns.
# For module testing, include path to module here
import sys
import os
modPath = r'D:\code\github\ePSproc'  # Win test machine
# modPath = r'/home/femtolab/github/ePSproc/'  # Linux test machine
sys.path.append(modPath)
import epsproc as ep
# TODO: tidy this up!
from epsproc.util import matEleSelector
from epsproc.geomFunc import geomCalc

Exploring Wigner 3js

In photoionization calculations, there is a lot of angular momentum coupling to deal with. Typically, 4 to 6 Wigner 3j terms appear (depending on the formalism), and/or higher-order terms in cases where couplings are included.

\begin{equation} W = \left(\begin{array}{ccc} l & l' & L\\ m & m' & M \end{array}\right) \end{equation}

Since this is, effectively, a 6D space, dimensions \((l_{max}, l_{max}, 2l_{max}, 2l_{max}+1, 2l_{max}+1, 4l_{max}+1)\) things can get large quickly. For small \(l_{max}\) it’s easy to look at some values directly…

(For more details on 3j symbols, see Wikipedia and Wolfram Mathworld; for more on the numerics, see the test notebook, and benchmarks.)

TODO: test numerically some of the symmetry properties here, by summation.

[ ]:



[109]:

# Calculate some values.
# w3jTable will output all values up to l=lp=Lmax (hence L=2Lmax)
lmax = 1
w3jlist = geomCalc.w3jTable(Lmax = lmax, form = '2d')  # For form = '2d', the function will output only valid entries as a coordinate table

print(w3jlist.shape)
print(f'Max value: {w3jlist[:,-1].max()}, min value: {w3jlist[:,-1].min()}\n')

# Print the table - output format has rows (l, lp, L, m, mp, M, 3j)
print(w3jlist)

(26, 7)
Max value: 1.0, min value: -0.5773502691896257

[[ 0.          0.          0.          0.          0.          0.
   1.        ]
 [ 0.          1.          1.          0.         -1.          1.
   0.57735027]
 [ 0.          1.          1.          0.          0.          0.
  -0.57735027]
 [ 0.          1.          1.          0.          1.         -1.
   0.57735027]
 [ 1.          0.          1.         -1.          0.          1.
   0.57735027]
 [ 1.          0.          1.          0.          0.          0.
  -0.57735027]
 [ 1.          0.          1.          1.          0.         -1.
   0.57735027]
 [ 1.          1.          2.         -1.         -1.          2.
   0.4472136 ]
 [ 1.          1.          1.         -1.          0.          1.
   0.40824829]
 [ 1.          1.          2.         -1.          0.          1.
  -0.31622777]
 [ 1.          1.          0.         -1.          1.          0.
   0.57735027]
 [ 1.          1.          1.         -1.          1.          0.
  -0.40824829]
 [ 1.          1.          2.         -1.          1.          0.
   0.18257419]
 [ 1.          1.          1.          0.         -1.          1.
  -0.40824829]
 [ 1.          1.          2.          0.         -1.          1.
  -0.31622777]
 [ 1.          1.          0.          0.          0.          0.
  -0.57735027]
 [ 1.          1.          1.          0.          0.          0.
   0.        ]
 [ 1.          1.          2.          0.          0.          0.
   0.36514837]
 [ 1.          1.          1.          0.          1.         -1.
   0.40824829]
 [ 1.          1.          2.          0.          1.         -1.
  -0.31622777]
 [ 1.          1.          0.          1.         -1.          0.
   0.57735027]
 [ 1.          1.          1.          1.         -1.          0.
   0.40824829]
 [ 1.          1.          2.          1.         -1.          0.
   0.18257419]
 [ 1.          1.          1.          1.          0.         -1.
  -0.40824829]
 [ 1.          1.          2.          1.          0.         -1.
  -0.31622777]
 [ 1.          1.          2.          1.          1.         -2.
   0.4472136 ]]

[110]:

# Recalculate and set to Xarray output format, then plot with ep.lmPlot()

w3j = geomCalc.w3jTable(Lmax = lmax, form = 'xdaLM')

# Check number of valid entries matches basic table above
print(f'Number of valid (non-NaN) elements: {w3j.count()}')

# Set parameters to restack the Xarray into (L,M) pairs
plotDimsRed = ['l', 'm', 'lp', 'mp']
xDim = {'LM':['L','M']}

# Plot with ep.lmPlot(), real values
daPlot, daPlotpd, legendList, gFig = ep.lmPlot(w3j, plotDims=plotDimsRed, xDim=xDim, pType = 'r')

Number of valid (non-NaN) elements: <xarray.DataArray 'w3jStacked' ()>
array(26)
Plotting data (No filename), pType=r, thres=0.01, with Seaborn

C:\Program Files (x86)\Microsoft Visual Studio\Shared\Anaconda3_64\envs\ePSdev\lib\site-packages\xarray\core\nputils.py:215: RuntimeWarning:

All-NaN slice encountered
../_images/methods_geometric_method_dev_260220_tidy_160320_7_2.png 
[111]:

# A complementary visulization is to call directly the sns.clustermap plot, use clustering and plot by category labels - see https://seaborn.pydata.org/index.html
# (ep.lpPlot uses a modified version of this routine.)
ep.snsMatMod.clustermap(daPlotpd.fillna(0), center=0, cmap="vlag", row_cluster=True, col_cluster=True)

[111]:

<epsproc._sns_matrixMod.ClusterGrid at 0x221eedaca20>
../_images/methods_geometric_method_dev_260220_tidy_160320_8_1.png 
[53]:

# Print out values by QNs (Pandas table)
daPlotpd

[53]:
L 0 1 2
M 0 -1 0 1 -2 -1 0 1 2
l m lp mp
0 0 0 0 1.00000 NaN NaN NaN NaN NaN NaN NaN NaN
1 -1 NaN NaN NaN 0.577350 NaN NaN NaN NaN NaN
0 NaN NaN -0.577350 NaN NaN NaN NaN NaN NaN
1 NaN 0.577350 NaN NaN NaN NaN NaN NaN NaN
1 -1 0 0 NaN NaN NaN 0.577350 NaN NaN NaN NaN NaN
1 -1 NaN NaN NaN NaN NaN NaN NaN NaN 0.447214
0 NaN NaN NaN 0.408248 NaN NaN NaN -0.316228 NaN
1 0.57735 NaN -0.408248 NaN NaN NaN 0.182574 NaN NaN
0 0 0 NaN NaN -0.577350 NaN NaN NaN NaN NaN NaN
1 -1 NaN NaN NaN -0.408248 NaN NaN NaN -0.316228 NaN
0 -0.57735 NaN 0.000000 NaN NaN NaN 0.365148 NaN NaN
1 NaN 0.408248 NaN NaN NaN -0.316228 NaN NaN NaN
1 0 0 NaN 0.577350 NaN NaN NaN NaN NaN NaN NaN
1 -1 0.57735 NaN 0.408248 NaN NaN NaN 0.182574 NaN NaN
0 NaN -0.408248 NaN NaN NaN -0.316228 NaN NaN NaN
1 NaN NaN NaN NaN 0.447214 NaN NaN NaN NaN

This ends up as a relatively sparse array, since many combinations are invalid (do not follow angular momentum selection rules), hence there are many NaN terms.

The results can also be output as a 6D sparse array, using the Sparse library.

[54]:

# Calculate and output in Sparse array format
w3jSparse = geomCalc.w3jTable(Lmax = lmax, form = 'ndsparse')
w3jSparse

[54]:
Formatcoo
Data Typefloat64
Shape(2, 2, 3, 2, 2, 3)
nnz25
Density0.1736111111111111
Read-onlyTrue
Size1.4K
Storage ratio1.2

Here nnz is the number of non-zero elements.

[55]:

# Try a larger Lmax and plot only.
lmax = 3
w3j = geomCalc.w3jTable(Lmax = lmax, form = 'xdaLM')

# Check number of valid entries matches basic table above
print(f'Number of valid (non-NaN) elements: {w3j.count()}')

# Set parameters to restack the Xarray into (L,M) pairs
plotDimsRed = ['l', 'm', 'lp', 'mp']
xDim = {'LM':['L','M']}

# Plot with ep.lmPlot(), real values
daPlot, daPlotpd, legendList, gFig = ep.lmPlot(w3j, plotDims=plotDimsRed, xDim=xDim, pType = 'r')

Number of valid (non-NaN) elements: <xarray.DataArray 'w3jStacked' ()>
array(820)
Plotting data (No filename), pType=r, thres=0.01, with Seaborn

C:\Program Files (x86)\Microsoft Visual Studio\Shared\Anaconda3_64\envs\ePSdev\lib\site-packages\xarray\core\nputils.py:215: RuntimeWarning:

All-NaN slice encountered
../_images/methods_geometric_method_dev_260220_tidy_160320_13_2.png 
[56]:

# Resort axis by (l,lp)
plotDimsRed = ['l', 'lp', 'm', 'mp']
xDim = {'LM':['L','M']}

# Plot with ep.lmPlot(), real values
daPlot, daPlotpd, legendList, gFig = ep.lmPlot(w3j, plotDims=plotDimsRed, xDim=xDim, pType = 'r')

C:\Program Files (x86)\Microsoft Visual Studio\Shared\Anaconda3_64\envs\ePSdev\lib\site-packages\xarray\core\nputils.py:215: RuntimeWarning:

All-NaN slice encountered


Plotting data (No filename), pType=r, thres=0.01, with Seaborn
../_images/methods_geometric_method_dev_260220_tidy_160320_14_2.png 
[57]:

# A complementary visulization is to call directly the sns.clustermap plot, use clustering and plot by category labels - see https://seaborn.pydata.org/index.html
# (ep.lpPlot uses a modified version of this routine.)
ep.snsMatMod.clustermap(daPlotpd.fillna(0), center=0, cmap="vlag", row_cluster=True, col_cluster=True)

[57]:

<epsproc._sns_matrixMod.ClusterGrid at 0x22183e4dfd0>
../_images/methods_geometric_method_dev_260220_tidy_160320_15_1.png

This clearly shows that the valid terms become sparser at higher \(l\), and the couplings become smaller.

Structure can also be examined using other methods, e.g. correlation functions (see, for example, Seaborn Discovering structure in heatmap data. The example here shows Panda’s standard Pearson correlation coefficient, which may (or may not) be particularly meaningful here… but does show structures.

[58]:

ep.snsMatMod.clustermap(daPlotpd.fillna(0).T.corr(), center=0, cmap="vlag", row_cluster=True, col_cluster=True)

[58]:

<epsproc._sns_matrixMod.ClusterGrid at 0x22183f526a0>
../_images/methods_geometric_method_dev_260220_tidy_160320_17_1.png 
[59]:

# Test some other Seaborn methods... these likely won't scale well for large lmax!
import seaborn as sns

# Recalculate for small lmax
lmax = 1
w3j = geomCalc.w3jTable(Lmax = lmax, form = 'xdaLM')

# Set parameters to restack the Xarray into (L,M) pairs
plotDimsRed = ['l', 'm', 'lp', 'mp']
xDim = {'LM':['L','M']}

# Plot with ep.lmPlot(), real values
daPlot, daPlotpd, legendList, gFig = ep.lmPlot(w3j, plotDims=plotDimsRed, xDim=xDim, pType = 'r')

# Try sns pairplot
# sns.pairplot(daPlotpd.fillna(0).T)  # Big grids!
# sns.pairplot(daPlotpd.fillna(0))  # OK, not particularly informative
# sns.pairplot(daPlotpd.fillna(0), hue = 'l')  # Doesn't work - multindex issue?

C:\Program Files (x86)\Microsoft Visual Studio\Shared\Anaconda3_64\envs\ePSdev\lib\site-packages\xarray\core\nputils.py:215: RuntimeWarning:

All-NaN slice encountered


Plotting data (No filename), pType=r, thres=0.01, with Seaborn
../_images/methods_geometric_method_dev_260220_tidy_160320_18_2.png

\(E_{P,R}\) tensor

The coupling of two 1-photon terms can be written as a tensor contraction:

\begin{equation} E_{PR}(\hat{e})=[e\otimes e^{*}]_{R}^{P}=[P]^{\frac{1}{2}}\sum_{p}(-1)^{R}\left(\begin{array}{ccc} 1 & 1 & P\\ p & R-p & -R \end{array}\right)e_{p}e_{R-p}^{*}\label{eq:EPR-defn-1} \end{equation}

Where \(e_{p}\) and \(e_{R-p}\) define the field strengths for the polarizations \(p\) and \(R-p\), which are coupled into the spherical tensor \(E_{PR}\).

[60]:

# Calculate EPR terms, all QNs, with field strengths e = 1
EPRtable = geomCalc.EPR(form = '2d')

# Output values as list, [l, lp, P, p, R-p, R, EPR]
print(EPRtable)

[[ 1.          1.          0.         -1.          1.          0.
   0.57735027]
 [ 1.          1.          1.         -1.          0.          1.
  -0.70710678]
 [ 1.          1.          1.         -1.          1.          0.
  -0.70710678]
 [ 1.          1.          2.         -1.         -1.          2.
   1.        ]
 [ 1.          1.          2.         -1.          0.          1.
   0.70710678]
 [ 1.          1.          2.         -1.          1.          0.
   0.40824829]
 [ 1.          1.          0.          0.          0.          0.
  -0.57735027]
 [ 1.          1.          1.          0.         -1.          1.
   0.70710678]
 [ 1.          1.          1.          0.          1.         -1.
  -0.70710678]
 [ 1.          1.          2.          0.         -1.          1.
   0.70710678]
 [ 1.          1.          2.          0.          0.          0.
   0.81649658]
 [ 1.          1.          2.          0.          1.         -1.
   0.70710678]
 [ 1.          1.          0.          1.         -1.          0.
   0.57735027]
 [ 1.          1.          1.          1.         -1.          0.
   0.70710678]
 [ 1.          1.          1.          1.          0.         -1.
   0.70710678]
 [ 1.          1.          2.          1.         -1.          0.
   0.40824829]
 [ 1.          1.          2.          1.          0.         -1.
   0.70710678]
 [ 1.          1.          2.          1.          1.         -2.
   1.        ]]

As before, we can visualise these values…

[61]:

EPRX = geomCalc.EPR(form = 'xarray')

# Set parameters to restack the Xarray into (L,M) pairs
# plotDimsRed = ['l', 'p', 'lp', 'R-p']
plotDimsRed = ['p', 'R-p']
xDim = {'PR':['P','R']}

# Plot with ep.lmPlot(), real values
daPlot, daPlotpd, legendList, gFig = ep.lmPlot(EPRX, plotDims=plotDimsRed, xDim=xDim, pType = 'r')

# Seem to have some all-NaN cols persisting here, not sure why...

Plotting data (No filename), pType=r, thres=0.01, with Seaborn

No handles with labels found to put in legend.
../_images/methods_geometric_method_dev_260220_tidy_160320_23_2.png 
[62]:

# Tabulate
daPlotpd.dropna(axis = 1, how = 'all')

[62]:
P 0 1 2
R 0 -1 0 1 -2 -1 0 1 2
p R-p
-1 -1 NaN NaN NaN NaN NaN NaN NaN NaN 1.0
0 NaN NaN NaN -0.707107 NaN NaN NaN 0.707107 NaN
1 0.57735 NaN -0.707107 NaN NaN NaN 0.408248 NaN NaN
0 -1 NaN NaN NaN 0.707107 NaN NaN NaN 0.707107 NaN
0 -0.57735 NaN NaN NaN NaN NaN 0.816497 NaN NaN
1 NaN -0.707107 NaN NaN NaN 0.707107 NaN NaN NaN
1 -1 0.57735 NaN 0.707107 NaN NaN NaN 0.408248 NaN NaN
0 NaN 0.707107 NaN NaN NaN 0.707107 NaN NaN NaN
1 NaN NaN NaN NaN 1.0 NaN NaN NaN NaN 
[63]:

# A complementary visulization is to call directly the sns.clustermap plot, use clustering and plot by category labels - see https://seaborn.pydata.org/index.html
# (ep.lpPlot uses a modified version of this routine.)
ep.snsMatMod.clustermap(daPlotpd.fillna(0), center=0, cmap="vlag", row_cluster=True, col_cluster=True)

[63]:

<epsproc._sns_matrixMod.ClusterGrid at 0x221812b0ac8>
../_images/methods_geometric_method_dev_260220_tidy_160320_25_1.png 
[64]:

ep.snsMatMod.clustermap(daPlotpd.fillna(0).T.corr(), center=0, cmap="vlag", row_cluster=True, col_cluster=True)

[64]:

<epsproc._sns_matrixMod.ClusterGrid at 0x2218055fdd8>
../_images/methods_geometric_method_dev_260220_tidy_160320_26_1.png

\(B_{L,M}\) term

The coupling of the partial wave pairs, \(|l,m\rangle\) and \(|l',m'\rangle\), into the observable set of \(\{L,M\}\) is defined by a tensor contraction with two 3j terms.

\begin{equation} B_{L,M}=(-1)^{m}\left(\frac{(2l+1)(2l'+1)(2L+1)}{4\pi}\right)^{1/2}\left(\begin{array}{ccc} l & l' & L\\ 0 & 0 & 0 \end{array}\right)\left(\begin{array}{ccc} l & l' & L\\ -m & m' & M \end{array}\right) \end{equation}

(See notebook ePSproc_wigner3j_dataStructures_260220.ipynb on Bemo for additonal dev details)

Note that this term is equivalent, effectively, as the triple intergral over spherical harmonics (e.g. eqn. 3.119 in Zare):

\begin{equation} \intop_{0}^{2\pi}\intop_{0}^{\pi}Y_{J_{3}M_{3}}(\theta,\phi)Y_{J_{2}M_{2}}(\theta,\phi)Y_{J_{1}M_{1}}(\theta,\phi)\sin\theta d\theta d\phi=\left(\frac{(2J_{1}+1)(2J_{2}+1)(2J_{3}+1)}{4\pi}\right)^{1/2}\left(\begin{array}{ccc} J_{1} & J_{2} & J_{3}\\ 0 & 0 & 0 \end{array}\right)\left(\begin{array}{ccc} J_{1} & J_{2} & J_{3}\\ M_{1} & M_{2} & M_{3} \end{array}\right) \end{equation}

… or contraction over a pair of harmonics into a resultant harmonic (e.g. eqns. C.21, C.22 in Blum), and this is how the term arises in the derivation.

\begin{equation} Y_{J_{1}M_{1}}(\theta,\phi)Y_{J_{2}M_{2}}(\theta,\phi)=\sum_{J_{3}M_{3}}\left(\frac{(2J_{1}+1)(2J_{2}+1)(2J_{3}+1)}{4\pi}\right)^{1/2}\left(\begin{array}{ccc} J_{1} & J_{2} & J_{3}\\ 0 & 0 & 0 \end{array}\right)\left(\begin{array}{ccc} J_{1} & J_{2} & J_{3}\\ M_{1} & M_{2} & M_{3} \end{array}\right)Y^*_{J_{3}M_{3}}(\theta,\phi) \end{equation}

(Again, for more details on see Wikipedia and Wolfram Mathworld.)

Note also some defns. use conjugate spherical harmonics, can convert as (Blum C.21):

\begin{equation} \beta_{L,M}^{\mu_{i},\mu_{f}}Y_{LM}^{*}(\theta_{\hat{k}},\phi_{\hat{k}})=\beta_{L,-M}^{\mu_{i},\mu_{f}}(-1)^{M}Y_{L,-M}(\theta_{\hat{k}},\phi_{\hat{k}})\label{eq:sph-conj-conv} \end{equation}
[65]:

Lmax = 1
BLMtable = geomCalc.betaTerm(Lmax = Lmax, form = 'xdaLM')  # Output as stacked Xarray

[66]:

BLMtable

[66]:

<xarray.DataArray 'w3jStacked' (lSet: 6, mSet: 9)>
array([[        nan,         nan,         nan,         nan,  0.28209479,
                nan,         nan,         nan,         nan],
       [        nan,         nan,         nan, -0.28209479,  0.28209479,
        -0.28209479,         nan,         nan,         nan],
       [        nan,  0.28209479,         nan,         nan,  0.28209479,
                nan,         nan,  0.28209479,         nan],
       [        nan,         nan,  0.28209479,         nan,  0.28209479,
                nan,  0.28209479,         nan,         nan],
       [        nan,         nan,         nan,         nan,         nan,
                nan,         nan,         nan,         nan],
       [-0.30901936,  0.21850969, -0.12615663, -0.21850969,  0.25231325,
        -0.21850969, -0.12615663,  0.21850969, -0.30901936]])
Coordinates:
  * mSet     (mSet) MultiIndex
  - m        (mSet) int64 -1 -1 -1 0 0 0 1 1 1
  - mp       (mSet) int64 -1 0 1 -1 0 1 -1 0 1
  - M        (mSet) int64 2 1 0 1 0 -1 0 -1 -2
  * lSet     (lSet) MultiIndex
  - l        (lSet) int64 0 0 1 1 1 1
  - lp       (lSet) int64 0 1 0 1 1 1
  - L        (lSet) int64 0 1 1 0 1 2
Attributes:
    dataType:  Wigner3j
[67]:

plotDimsRed = ['l', 'm', 'lp', 'mp']
xDim = {'LM':['L','M']}
# daPlot = BLMtable.unstack().stack({'LM':['L','M']})
# ep.lmPlot(w3jXcombMult.unstack().stack({'LM':['l','m']}), plotDims=['lp', 'L', 'mp', 'M'], xDim='LM', SFflag = False)
# daPlot, daPlotpd, legendList, gFig = ep.lmPlot(BLMtable.unstack().stack(xDim), plotDims=plotDimsRed, xDim=xDim, SFflag = False, squeeze = False)
# daPlot, daPlotpd, legendList, gFig = ep.lmPlot(BLMtable, plotDims=plotDimsRed, xDim=xDim)
# ep.lmPlot(w3jXcombMult.unstack(), plotDims=['lp', 'L', 'mp', 'M'], xDim='L', SFflag = False)

# daPlotpd = daPlot.unstack().stack(plotDim = plotDimsRed).to_pandas().dropna(axis = 1).T

[68]:

# daPlot, daPlotpd, legendList, gFig = ep.lmPlot(BLMtable.unstack().stack(xDim), plotDims=plotDimsRed, xDim=xDim, pType = 'r')
daPlot, daPlotpd, legendList, gFig = ep.lmPlot(BLMtable, plotDims=plotDimsRed, xDim=xDim, pType = 'r')

C:\Program Files (x86)\Microsoft Visual Studio\Shared\Anaconda3_64\envs\ePSdev\lib\site-packages\xarray\core\nputils.py:215: RuntimeWarning:

All-NaN slice encountered


Plotting data (No filename), pType=r, thres=0.01, with Seaborn
../_images/methods_geometric_method_dev_260220_tidy_160320_32_2.png 
[69]:

daPlotpd

[69]:
L 0 1 2
M 0 -1 0 1 -2 -1 0 1 2
l m lp mp
0 0 0 0 0.282095 NaN NaN NaN NaN NaN NaN NaN NaN
1 -1 NaN NaN NaN -0.282095 NaN NaN NaN NaN NaN
0 NaN NaN 0.282095 NaN NaN NaN NaN NaN NaN
1 NaN -0.282095 NaN NaN NaN NaN NaN NaN NaN
1 -1 0 0 NaN NaN NaN 0.282095 NaN NaN NaN NaN NaN
1 -1 NaN NaN NaN NaN NaN NaN NaN NaN -0.309019
0 NaN NaN NaN NaN NaN NaN NaN 0.21851 NaN
1 0.282095 NaN NaN NaN NaN NaN -0.126157 NaN NaN
0 0 0 NaN NaN 0.282095 NaN NaN NaN NaN NaN NaN
1 -1 NaN NaN NaN NaN NaN NaN NaN -0.21851 NaN
0 0.282095 NaN NaN NaN NaN NaN 0.252313 NaN NaN
1 NaN NaN NaN NaN NaN -0.21851 NaN NaN NaN
1 0 0 NaN 0.282095 NaN NaN NaN NaN NaN NaN NaN
1 -1 0.282095 NaN NaN NaN NaN NaN -0.126157 NaN NaN
0 NaN NaN NaN NaN NaN 0.21851 NaN NaN NaN
1 NaN NaN NaN NaN -0.309019 NaN NaN NaN NaN 
[70]:

# A complementary visulization is to call directly the sns.clustermap plot, use clustering and plot by category labels - see https://seaborn.pydata.org/index.html
# (ep.lpPlot uses a modified version of this routine.)
ep.snsMatMod.clustermap(daPlotpd.fillna(0).T, center=0, cmap="vlag", row_cluster=True, col_cluster=True)

[70]:

<epsproc._sns_matrixMod.ClusterGrid at 0x221eed13a20>
../_images/methods_geometric_method_dev_260220_tidy_160320_34_1.png 
[71]:

ep.snsMatMod.clustermap(daPlotpd.fillna(0).T.corr(), center=0, cmap="vlag", row_cluster=True, col_cluster=True)

[71]:

<epsproc._sns_matrixMod.ClusterGrid at 0x22183ac17b8>
../_images/methods_geometric_method_dev_260220_tidy_160320_35_1.png 
[72]:

ep.snsMatMod.clustermap(daPlotpd.fillna(0).corr(), center=0, cmap="vlag", row_cluster=True, col_cluster=True)

[72]:

<epsproc._sns_matrixMod.ClusterGrid at 0x221839db400>
../_images/methods_geometric_method_dev_260220_tidy_160320_36_1.png

What do we learn here…? This is a nice way to visualize the selection rules into the observable: for instance, only terms \(l=l'\) and \(m=-m'\) contribute to the overall photoinoization cross-section term (\(L=0, M=0\)). However, since these terms are fairly simply followed algebraically in this case, via the rules inherent in the 3j product, this is not particularly insightful. These visualizations will become more useful when dealing with real sets of matrix elements, and specific polarization geometries, which will further modulate the \(B_{L,M}\) terms.

[73]:

# PLOT AS A FN. of dl, dm... should just need to recalc dims or add a separate label for this?
# dl = BLMtable.l - BLMtable.lp
# dm = BLMtable.m - BLMtable.mp

BLMdeltas = BLMtable.unstack().stack({'LM':['L','M'], 'dl':['l','lp'], 'dm':['m','mp']})
dl = BLMdeltas.l - BLMdeltas.lp
dm = BLMdeltas.m - BLMdeltas.mp

# Set new dims & plot
BLMdeltas = BLMdeltas.assign_coords(dl=dl, dm=dm)  # This works, but not plotting correctly for some reason - maybe dim names?
plotDimsRed = ['dl','dm']
# BLMdeltas = BLMdeltas.assign_coords(dl=dl, dm=dm).rename({'dl':'l', 'dm':'m'})  # Hmm, same result - values present in daPlotpd not showing up here, or losing correlations...?
# plotDimsRed = ['l','m']
xDim = {'LM':['L','M']}
daPlot, daPlotpd, legendList, gFig = ep.lmPlot(BLMdeltas, plotDims=plotDimsRed, xDim=xDim, pType = 'r')

C:\Program Files (x86)\Microsoft Visual Studio\Shared\Anaconda3_64\envs\ePSdev\lib\site-packages\xarray\core\nputils.py:215: RuntimeWarning:

All-NaN slice encountered


Plotting data (No filename), pType=r, thres=0.01, with Seaborn

No handles with labels found to put in legend.
../_images/methods_geometric_method_dev_260220_tidy_160320_38_3.png 
[74]:

daPlotpd

[74]:
L 0 1 2
M 0 -1 0 1 -2 -1 0 1 2
dl dm
0 0 0.282095 NaN NaN NaN NaN NaN NaN NaN NaN
-1 1 NaN NaN NaN -0.282095 NaN NaN NaN NaN NaN
0 NaN NaN 0.282095 NaN NaN NaN NaN NaN NaN
-1 NaN -0.282095 NaN NaN NaN NaN NaN NaN NaN
1 -1 NaN NaN NaN 0.282095 NaN NaN NaN NaN NaN
0 NaN NaN 0.282095 NaN NaN NaN NaN NaN NaN
1 NaN 0.282095 NaN NaN NaN NaN NaN NaN NaN
0 0 NaN NaN NaN NaN NaN NaN NaN NaN -0.309019
-1 NaN NaN NaN NaN NaN NaN NaN 0.21851 NaN
-2 0.282095 NaN NaN NaN NaN NaN -0.126157 NaN NaN
1 NaN NaN NaN NaN NaN NaN NaN -0.21851 NaN
0 0.282095 NaN NaN NaN NaN NaN 0.252313 NaN NaN
-1 NaN NaN NaN NaN NaN -0.21851 NaN NaN NaN
2 0.282095 NaN NaN NaN NaN NaN -0.126157 NaN NaN
1 NaN NaN NaN NaN NaN 0.21851 NaN NaN NaN
0 NaN NaN NaN NaN -0.309019 NaN NaN NaN NaN 
[75]:

ep.snsMatMod.clustermap(daPlotpd.fillna(0))

[75]:

<epsproc._sns_matrixMod.ClusterGrid at 0x221840ae048>
../_images/methods_geometric_method_dev_260220_tidy_160320_40_1.png

Note correlation of terms, with final \(L\) and \(M\) odd or even terms grouped according to odd or even \(\Delta l\) and \(\Delta m\) respectively.

[76]:

ep.snsMatMod.clustermap(daPlotpd.fillna(0).T.corr())

[76]:

<epsproc._sns_matrixMod.ClusterGrid at 0x22181071d68>
../_images/methods_geometric_method_dev_260220_tidy_160320_42_1.png 
[79]:

# Test correlation fns.
# This fails with NaNs present it seems
# ep.snsMatMod.clustermap(daPlot.dropna(dim='plotDim', how='all').dropna(dim='LM',how='all').fillna(0).to_pandas().T.corr())
# ep.snsMatMod.clustermap(daPlot.dropna(dim='plotDim', how='all').dropna(dim='LM',how='all').fillna(0).to_pandas().corr())
ep.snsMatMod.clustermap(daPlotpd.fillna(0).T.corr())
ep.snsMatMod.clustermap(daPlotpd.fillna(0).corr())

[79]:

<epsproc._sns_matrixMod.ClusterGrid at 0x221817b2828>
../_images/methods_geometric_method_dev_260220_tidy_160320_43_1.png
../_images/methods_geometric_method_dev_260220_tidy_160320_43_2.png 
[80]:

# Switch plotting dims
plotDimsRed = ['L','M']
xDim = {'llpmmp':['l', 'm', 'lp', 'mp']}
daPlot, daPlotpd, legendList, gFig = ep.lmPlot(BLMtable, plotDims=plotDimsRed, xDim=xDim, mMax = 2, pType = 'r')

C:\Program Files (x86)\Microsoft Visual Studio\Shared\Anaconda3_64\envs\ePSdev\lib\site-packages\xarray\core\nputils.py:215: RuntimeWarning:

All-NaN slice encountered


Plotting data (No filename), pType=r, thres=0.01, with Seaborn

No handles with labels found to put in legend.
../_images/methods_geometric_method_dev_260220_tidy_160320_44_3.png

\(\Lambda\) Term

Define MF projection term, \(\Lambda_{R',R}(R_{\hat{n}})\):

\begin{equation} \Lambda_{R',R}(R_{\hat{n}})=(-1)^{(R')}\left(\begin{array}{ccc} 1 & 1 & P\\ \mu & -\mu' & R' \end{array}\right)D_{-R',-R}^{P}(R_{\hat{n}}) \end{equation}

This is similar to the \(E_{PR}\) term, and essentially rotates it into the MF.

[81]:

# Calculate MF projection terms, all QNs, for Euler angles corresponding to (z,x,y) cases
lambdaTerm, lambdaTable, lambdaD, QNs = geomCalc.MFproj(form = '2d')

# Output values as list, ['l', 'lp', 'P', 'mu', 'mup', 'Rp', 'R']
# print(lambdaTerm.shape)
# print(np.c_[QNs[0:10,:], lambdaTerm[0:10,:]])

[82]:

# Xarray form
lambdaTerm, lambdaTable, lambdaD, QNs = geomCalc.MFproj(form = 'xarray')

[83]:

# Check values Xarray vs. 2d case for consistency
# TO DO

[84]:

# In this case the QNs are left unstacked to facilitate multiplication across different terms - this fails for stacked (multindex) cases here.
lambdaTerm.dims

[84]:

('Rp', 'l', 'lp', 'P', 'mu', 'mup', 'R', 'Euler')

[85]:

# Plot all QNs vs. Euler angles (by label)
plotDimsRed = ['l', 'lp', 'P', 'mu', 'mup', 'R', 'Rp']
xDim = 'Labels'
daPlot, daPlotpd, legendList, gFig = ep.lmPlot(lambdaTerm, plotDims=plotDimsRed, xDim=xDim, pType = 'a')   # Plot abs values
daPlotPhase, daPlotpdPhase, legendList, gFig = ep.lmPlot(lambdaTerm, plotDims=plotDimsRed, xDim=xDim, pType = 'phase')   # Plot phase

Plotting data (No filename), pType=a, thres=0.01, with Seaborn
Plotting data (No filename), pType=phase, thres=0.01, with Seaborn
../_images/methods_geometric_method_dev_260220_tidy_160320_50_1.png
../_images/methods_geometric_method_dev_260220_tidy_160320_50_2.png 
[86]:

daPlotpd

[86]:
Labels z x y
l lp P mu mup R Rp
1 1 0 -1 1 0 0 0.577350 0.577350 0.577350
0 0 0 0 0.577350 0.577350 0.577350
1 -1 0 0 0.577350 0.577350 0.577350
1 -1 0 -1 1 0.000000 0.204124 0.204124
0 1 0.000000 0.288675 0.288675
1 1 0.408248 0.204124 0.204124
1 -1 0 0.000000 0.288675 0.288675
0 0 0.408248 0.000000 0.000000
1 0 0.000000 0.288675 0.288675
0 -1 -1 1 0.000000 0.204124 0.204124
0 1 0.000000 0.288675 0.288675
1 1 0.408248 0.204124 0.204124
1 -1 -1 0.408248 0.204124 0.204124
0 -1 0.000000 0.288675 0.288675
1 -1 0.000000 0.204124 0.204124
1 -1 -1 0 0.000000 0.288675 0.288675
0 0 0.408248 0.000000 0.000000
1 0 0.000000 0.288675 0.288675
0 -1 -1 0.408248 0.204124 0.204124
0 -1 0.000000 0.288675 0.288675
1 -1 0.000000 0.204124 0.204124
2 -1 -1 -2 2 0.000000 0.111803 0.111803
-1 2 0.000000 0.223607 0.223607
0 2 0.000000 0.273861 0.273861
1 2 0.000000 0.223607 0.223607
2 2 0.447214 0.111803 0.111803
0 -2 1 0.000000 0.158114 0.158114
-1 1 0.000000 0.158114 0.158114
1 1 0.316228 0.158114 0.158114
2 1 0.000000 0.158114 0.158114
1 -2 0 0.000000 0.111803 0.111803
0 0 0.182574 0.091287 0.091287
2 0 0.000000 0.111803 0.111803
0 -1 -2 1 0.000000 0.158114 0.158114
-1 1 0.000000 0.158114 0.158114
1 1 0.316228 0.158114 0.158114
2 1 0.000000 0.158114 0.158114
0 -2 0 0.000000 0.223607 0.223607
0 0 0.365148 0.182574 0.182574
2 0 0.000000 0.223607 0.223607
1 -2 -1 0.000000 0.158114 0.158114
-1 -1 0.316228 0.158114 0.158114
1 -1 0.000000 0.158114 0.158114
2 -1 0.000000 0.158114 0.158114
1 -1 -2 0 0.000000 0.111803 0.111803
0 0 0.182574 0.091287 0.091287
2 0 0.000000 0.111803 0.111803
0 -2 -1 0.000000 0.158114 0.158114
-1 -1 0.316228 0.158114 0.158114
1 -1 0.000000 0.158114 0.158114
2 -1 0.000000 0.158114 0.158114
1 -2 -2 0.447214 0.111803 0.111803
-1 -2 0.000000 0.223607 0.223607
0 -2 0.000000 0.273861 0.273861
1 -2 0.000000 0.223607 0.223607
2 -2 0.000000 0.111803 0.111803 
[87]:

# Test correlation fns.
# ep.snsMatMod.clustermap(daPlotpd.fillna(0).T.corr())  # This currently failing with NaNs or inf?
ep.snsMatMod.clustermap(daPlotpd.fillna(0).corr())

[87]:

<epsproc._sns_matrixMod.ClusterGrid at 0x221eeaef1d0>
../_images/methods_geometric_method_dev_260220_tidy_160320_52_1.png

\(\beta_{L,M}^{MF}\) rewrite

The various terms defined above can be used to redefine the full MF observables, expressed as a set of \(\beta_{L,M}\) coefficients.

The original (full) form for the MF equation, as implemented in ePSproc (note this is defined for an expansion in non-conjugate harmonics, \(Y_{LM}(\theta_{\hat{k}},\phi_{\hat{k}})\)):

\begin{eqnarray} \beta_{L,-M}^{\mu_{i},\mu_{f}} & = & \sum_{l,m,\mu}\sum_{l',m',\mu'}(-1)^{M}(-1)^{m}(-1)^{(\mu'-\mu_{0})}\left(\frac{(2l+1)(2l'+1)(2L+1)}{4\pi}\right)^{1/2}\left(\begin{array}{ccc} l & l' & L\\ 0 & 0 & 0 \end{array}\right)\left(\begin{array}{ccc} l & l' & L\\ -m & m' & -M \end{array}\right)\nonumber \\ & \times & \sum_{P,R',R}(2P+1)(-1)^{(R'-R)}\left(\begin{array}{ccc} 1 & 1 & P\\ \mu & -\mu' & R' \end{array}\right)\left(\begin{array}{ccc} 1 & 1 & P\\ \mu_{0} & -\mu_{0} & R \end{array}\right)D_{-R',-R}^{P}(R_{\hat{n}})I_{l,m,\mu}^{p_{i}\mu_{i},p_{f}\mu_{f}}(E)I_{l',m',\mu'}^{p_{i}\mu_{i},p_{f}\mu_{f}*}(E) \end{eqnarray}

Where \(I_{l,m,\mu}^{p_{i}\mu_{i},p_{f}\mu_{f}}(E)\) are the energy-dependent dipole matrix elements.

Substituting the terms above (plus some tidying up with phases etc. - see other notes for full derivation):

\begin{eqnarray} \beta_{L,-M}^{\mu_{i},\mu_{f}} & = & \sum_{l,m,\mu}\sum_{l',m',\mu'}(-1)^{(\mu'-\mu_{0})}{\color{red}B_{L,-M}}\nonumber \\ & \times & \sum_{P,R',R}{\color{red}E_{P-R}(\hat{e})\Lambda_{R',R}(R_{\hat{n}})}I_{l,m,\mu}^{p_{i}\mu_{i},p_{f}\mu_{f}}(E)I_{l',m',\mu'}^{p_{i}\mu_{i},p_{f}\mu_{f}*}(E) \end{eqnarray}

And can also reorder…

\begin{eqnarray*} \beta_{L,-M}^{\mu_{i},\mu_{f}} & = & \sum_{P,R',R}{\color{red}E_{P-R}(\hat{e})}\sum_{l,m,\mu}\sum_{l',m',\mu'}(-1)^{(\mu'-\mu_{0})}{\color{red}B_{L,-M}\Lambda_{R',R}(R_{\hat{n}})}I_{l,m,\mu}^{p_{i}\mu_{i},p_{f}\mu_{f}}(E)I_{l',m',\mu'}^{p_{i}\mu_{i},p_{f}\mu_{f}*}(E) \end{eqnarray*}

Add additional (function) labels to track QNs:

\begin{eqnarray*} \beta_{L,-M}^{\mu_{i},\mu_{f}} & = & \sum_{P,R',R}{\color{red}E_{P-R}(\hat{e};\mu_{0})}\sum_{l,m,\mu}\sum_{l',m',\mu'}(-1)^{(\mu'-\mu_{0})}{\color{red}\Lambda_{R',R}(R_{\hat{n}};\mu,P,R,R')B_{L,-M}(l,l',m,m')}I_{l,m,\mu}^{p_{i}\mu_{i},p_{f}\mu_{f}}(E)I_{l',m',\mu'}^{p_{i}\mu_{i},p_{f}\mu_{f}*}(E) \end{eqnarray*}

Cf. R&U rotational wavepacket case, see also Chpt. 12 in Quantum Met book for similar formalism separating geometric channels:

\begin{equation} \varUpsilon_{L,M}^{u}=\sum_{ll'}\sum_{\lambda\lambda'}\sum_{mm'}(-1)^{m}\sqrt{\frac{(2l+1)(2l'+1)(2L+1)}{4\pi}}\left(\begin{array}{ccc} l & l' & L\\ m & -m' & M \end{array}\right)\left(\begin{array}{ccc} l & l' & L\\ 0 & 0 & 0 \end{array}\right)\gamma_{\alpha_{+}l\lambda ml'\lambda'm'}^{u}\label{eq:geom-channels} \end{equation}
\begin{equation} \Lambda^{\zeta\zeta'}=\sum_{\Gamma,\Gamma'}\sum_{\mu,\mu'}\sum_{h,h'}b_{hl\lambda}^{\Gamma\mu*}b_{h'l'\lambda'}^{\Gamma'\mu'}\boldsymbol{D}_{hl}^{\Gamma\mu*}(q,\,k)\boldsymbol{D}_{h'l'}^{\Gamma'\mu'}(q',\,k) \end{equation}
\begin{equation} \beta_{L,M}^{u}=\sum_{\zeta,\zeta'}\varUpsilon_{L,M}^{u,\zeta\zeta'}\Lambda^{\zeta\zeta'} \end{equation}

Testing multiplications…

[88]:

lmax = 2
p = [0]  # Set single pol state only - NOTE EPR() currently sets all ep=1, need to fix this!

# Calculate various tensors...

# *** EPR
EPRX = geomCalc.EPR(form = 'xarray', p = p)

# Set parameters to restack the Xarray into (L,M) pairs
# plotDimsRed = ['l', 'p', 'lp', 'R-p']
plotDimsRed = ['p', 'R-p']
xDim = {'PR':['P','R']}

# Plot with ep.lmPlot(), real values
daPlot, daPlotpd, legendList, gFig = ep.lmPlot(EPRX, plotDims=plotDimsRed, xDim=xDim, pType = 'r')


# *** Blm term to lmax
BLMtable = geomCalc.betaTerm(Lmax = lmax, form = 'xdaLM')

plotDimsRed = ['l', 'm', 'lp', 'mp']
xDim = {'LM':['L','M']}

daPlot, daPlotpd, legendList, gFig = ep.lmPlot(BLMtable, plotDims=plotDimsRed, xDim=xDim, pType = 'r')


# *** Lambda term
lambdaTerm, lambdaTable, lambdaD, QNs = geomCalc.MFproj(form = 'xarray')

plotDimsRed = ['l', 'lp', 'P', 'mu', 'mup', 'R', 'Rp']
xDim = 'Labels'
daPlot, daPlotpd, legendList, gFig = ep.lmPlot(lambdaTerm, plotDims=plotDimsRed, xDim=xDim, pType = 'a')   # Plot abs values


#%% Basic multiplication (no matrix elements)

# Fix/remove QNs for multiplication

# EPRXresort = EPRX.unstack().squeeze().drop('R-p')  # Not quite right - need to flatten on R-p?
# EPRXresort = EPRX.unstack().squeeze()  # This removes singleton dims from data, but leaves labels.
EPRXresort = EPRX.unstack().squeeze().drop('l').drop('lp')  # This removes (l,lp) dims fully.

lambdaTermResort = lambdaTerm.squeeze().drop('l').drop('lp')

BLMtableResort = BLMtable.unstack()

# Apply additional phase conventions
Rphase = True
Mphase = True
mupPhase = True

# Set phase conventions +ve <> -ve terms by in-place coord multiplications
# This should allow for simple selection of values by coords.
if Rphase:
    EPRXresort['R'] *= -1

if Mphase:
    BLMtableResort['M'] *= -1

# Test big mult...
mTerm = EPRXresort * lambdaTermResort * BLMtableResort  # Multiplication works OK... but might be an ugly result...

# Set additional phase term, (-1)^(mup-p)
if mupPhase:
    mupPhaseTerm = np.power(-1, np.abs(mTerm.mup - mTerm.p))
    mTerm = mTerm * mupPhaseTerm

mTerm.attrs['dataType'] = 'MulTest'
mTerm.attrs['file'] = 'MulTest'  # Temporarily adding this, not sure why this is an issue here however (not an issue for other cases...)

Plotting data (No filename), pType=r, thres=0.01, with Seaborn

No handles with labels found to put in legend.
C:\Program Files (x86)\Microsoft Visual Studio\Shared\Anaconda3_64\envs\ePSdev\lib\site-packages\xarray\core\nputils.py:215: RuntimeWarning:

All-NaN slice encountered


Plotting data (No filename), pType=r, thres=0.01, with Seaborn
Plotting data (No filename), pType=a, thres=0.01, with Seaborn
../_images/methods_geometric_method_dev_260220_tidy_160320_58_3.png
../_images/methods_geometric_method_dev_260220_tidy_160320_58_4.png
../_images/methods_geometric_method_dev_260220_tidy_160320_58_5.png 
[89]:

mTerm.dims

[89]:

('P', 'R-p', 'R', 'Rp', 'mu', 'mup', 'Euler', 'l', 'lp', 'L', 'm', 'mp', 'M')

[90]:

mTerm.shape

[90]:

(3, 3, 3, 5, 3, 3, 3, 3, 3, 5, 5, 5, 9)

[91]:

# Try making sense of this...
# NOTE throwing values to console is very slow here due to high number of dims... calcs and plotting are OK however.

# sumDims = ['l', 'lp', 'P', 'mu', 'mup', 'Rp', 'm', 'mp']  # Define dims to sum over
sumDims = ['P', 'mu', 'mup', 'Rp', ]  # Define dims to sum over
# selDims = {'p':0, 'R':0, 'R-p':0}  # Use selDims to select pol state
selDims = {'R':0, 'R-p':0}  # Use selDims to select pol state
# plotDimsRed = ['L','M']
plotDimsRed = ['L','M', 'l', 'lp', 'm', 'mp']
xDim = 'Labels'

daPlot, daPlotpd, legendList, gFig = ep.lmPlot(mTerm, plotDims=plotDimsRed, xDim=xDim, sumDims=sumDims, selDims=selDims, pType = 'a')   # Plot abs values

Plotting data MulTest, pType=a, thres=0.01, with Seaborn
../_images/methods_geometric_method_dev_260220_tidy_160320_61_1.png 
[92]:

# Single pol state and geometry....

sumDims = ['P', 'mu', 'mup', 'Rp', ]  # Define dims to sum over
# selDims = {'p':0, 'R':0, 'R-p':0,'Labels':'z'}  # Use selDims to select pol state
selDims = {'R':0, 'R-p':0,'Labels':'z'}  # Use selDims to select pol state
# plotDimsRed = ['L','M']
plotDimsRed = ['l', 'lp', 'm', 'mp']
xDim = {'LM':['L','M']}

daPlot, daPlotpd, legendList, gFig = ep.lmPlot(mTerm, plotDims=plotDimsRed, xDim=xDim, sumDims=sumDims, selDims=selDims, pType = 'a')   # Plot abs values
# This currently includes many non-valid (L,M) - need to further reindex before plotting?  They are not present in full plot above.

Plotting data MulTest, pType=a, thres=0.01, with Seaborn
../_images/methods_geometric_method_dev_260220_tidy_160320_62_1.png 
[93]:

# daPlotpd  # .dropna(how='any',axis=1)  - many invalid null terms here??? (Hence can't drop with dropna)
daPlotpd[daPlotpd.abs() >= 0.01].dropna(how='all', axis=1)  # Slightly better!

[93]:
L 0 1 2 3 4
M 0 1 0 -1 2 1 0 -1 -2 3 ... -3 4 3 2 1 0 -1 -2 -3 -4
l lp m mp
0 0 0 0 0.282095 NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
1 0 -1 NaN NaN NaN 0.282095 NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
0 NaN NaN 0.282095 NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
1 NaN 0.282095 NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2 0 -2 NaN NaN NaN NaN NaN NaN NaN NaN 0.282095 NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
2 2 2 -2 0.282095 NaN NaN NaN NaN NaN 0.180224 NaN NaN NaN ... NaN NaN NaN NaN NaN 0.040299 NaN NaN NaN NaN
-1 NaN NaN NaN NaN NaN 0.220728 NaN NaN NaN NaN ... NaN NaN NaN NaN 0.090112 NaN NaN NaN NaN NaN
0 NaN NaN NaN NaN 0.180224 NaN NaN NaN NaN NaN ... NaN NaN NaN 0.156078 NaN NaN NaN NaN NaN NaN
1 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN 0.238414 NaN NaN NaN NaN NaN NaN NaN
2 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN 0.337168 NaN NaN NaN NaN NaN NaN NaN NaN

81 rows × 25 columns

[94]:

# With cluster sorting
ep.snsMatMod.clustermap(daPlotpd[daPlotpd.abs() >= 0.01].dropna(how='all', axis=1).fillna(0))

[94]:

<epsproc._sns_matrixMod.ClusterGrid at 0x2218048de48>
../_images/methods_geometric_method_dev_260220_tidy_160320_64_1.png

Calculations with matrix elements to be explored in the next installment…