Basic PEMtk fitting class demo

03/06/21

Tidy version for docs, derived from http://127.0.0.1:8888/lab/tree/dev/PEMtk/fitting/fitting_routines_class-demo_v1_110521-tidy-StimpyTest_020621.ipynb

11/05/21

First version of class test/demo.

See dev notebook for background.

Outline of this notebook:

  • Load required packages.

  • Setup pemtkFit object.

    • Set various parameters, either from existing data or new values.

    • Data is handled as a set of dictionaries within the class, self.data[key][dataType], where key is an arbitrary label for, e.g. a specific experiment, calculation etc, and dataType contains a set of values, parameters etc. (Should become clear below!)

    • Methods operate on all self.data items in general, with some special cases: self.data['subset'] contains data to be used in fitting.

  • Simulate data

    • Use ePSproc to simulated aligned-frame measurements.

  • Fit data

Prerequisities

  • Working installation of ePSproc + PEMtk (or local copies of the Git repos, which can be pointed at for setup below).

  • Test/demo data, from ePSproc Github repo.

Setup

Imports

A few standard imports…

[1]:
import sys
import os
from pathlib import Path
# import numpy as np
# import epsproc as ep
# import xarray as xr

from datetime import datetime as dt
timeString = dt.now()

And local module imports. This should work either for installed versions (e.g. via pip install), or for test code via setting the base path below to point at your local copies.

[2]:
# For module testing, include path to module here, otherwise use global installation
if sys.platform == "win32":
    modPath = Path(r'D:\code\github')  # Win test machine
    winFlag = True
else:
    modPath = Path(r'/home/femtolab/github')  # Linux test machine
    winFlag = False

# Append to sys path
sys.path.append((modPath/'ePSproc').as_posix())
sys.path.append((modPath/'PEMtk').as_posix())
[3]:
# ePSproc
import epsproc as ep

# Set data path
# Note this is set here from ep.__path__, but may not be correct in all cases - depends on where the Github repo is.
epDemoDataPath = Path(ep.__path__[0]).parent/'data'
[4]:
# PEMtk
# import pemtk as pm
# from pemtk.data.dataClasses import dataClass

# Import fitting class
from pemtk.fit.fitClass import pemtkFit
[5]:
# Set HTML output style for Xarray in notebooks (optional), may also depend on version of Jupyter notebook or lab, or Xr
# See http://xarray.pydata.org/en/stable/generated/xarray.set_options.html
# if isnotebook():
# xr.set_options(display_style = 'html')
[6]:
# Set some plot options
ep.plot.hvPlotters.setPlotters()

Set & load parameters

There are a few things that need to be configured for a given case…

  • Matrix elements (or (l,m,mu) indicies) to use.

  • Alignment distribution (ADMs).

  • Polarization geometry.

In this demo, a real case will first be simulated with computational values, and then used to test the fitting routines.

(TODO: demo from scratch without known matrix elements.)

For fit testing, start with computational values for the matrix elements. These will be used to simulate data, and also to provide a list of parameters to fit later.

[7]:
# Set for ePSproc test data, available from https://github.com/phockett/ePSproc/tree/master/data
# Note this is set here from ep.__path__, but may not be correct in all cases.

# Multiorb data
dataPath = os.path.join(epDemoDataPath, 'photoionization', 'n2_multiorb')
[8]:
data = pemtkFit(fileBase = dataPath, verbose = 1)
[9]:
# Read data files
data.scanFiles()
data.jobsSummary()

*** Job orb6 details
Key: orb6
Dir D:\code\github\ePSproc\data\photoionization\n2_multiorb, 1 file(s).
{   'batch': 'ePS n2, batch n2_1pu_0.1-50.1eV, orbital A2',
    'event': ' N2 A-state (1piu-1)',
    'orbE': -17.096913836366,
    'orbLabel': '1piu-1'}

*** Job orb5 details
Key: orb5
Dir D:\code\github\ePSproc\data\photoionization\n2_multiorb, 1 file(s).
{   'batch': 'ePS n2, batch n2_3sg_0.1-50.1eV, orbital A2',
    'event': ' N2 X-state (3sg-1)',
    'orbE': -17.341816310545997,
    'orbLabel': '3sg-1'}

*** Job orb6 details
Key: orb6
Dir D:\code\github\ePSproc\data\photoionization\n2_multiorb, 1 file(s).
{   'batch': 'ePS n2, batch n2_1pu_0.1-50.1eV, orbital A2',
    'event': ' N2 A-state (1piu-1)',
    'orbE': -17.096913836366,
    'orbLabel': '1piu-1'}

*** Job orb5 details
Key: orb5
Dir D:\code\github\ePSproc\data\photoionization\n2_multiorb, 1 file(s).
{   'batch': 'ePS n2, batch n2_3sg_0.1-50.1eV, orbital A2',
    'event': ' N2 X-state (3sg-1)',
    'orbE': -17.341816310545997,
    'orbLabel': '3sg-1'}

The class wraps ep.setADMs(). This returns an isotropic distribution by default, or values can be set explicitly from a list. Values are set in self.data['ADM'].

Note: if this is not set, the default value will be used, which is likely not very useful for the fit!

[10]:
# Default case
data.setADMs()
# data.ADM['ADMX']
data.data['ADM']['ADM']
[10]:
Show/Hide data repr Show/Hide attributes
xarray.DataArray
'ADM'
  • ADM: 1
  • t: 1
  • 1
    array([[1]])
    • ADM
      (ADM)
      MultiIndex
      (K, Q, S)
      array([(0, 0, 0)], dtype=object)
    • K
      (ADM)
      int64
      0
      array([0], dtype=int64)
    • Q
      (ADM)
      int64
      0
      array([0], dtype=int64)
    • S
      (ADM)
      int64
      0
      array([0], dtype=int64)
    • t
      (t)
      int32
      0
      units :
      Index
      array([0])
  • dataType :
    ADM
    long_name :
    Axis distribution moments
    units :
    arb
[11]:
# Load time-dependent ADMs for N2 case
# Adapted from ePSproc_AFBLM_testing_010519_300719.m

from scipy.io import loadmat
ADMdataFile = os.path.join(epDemoDataPath, 'alignment', 'N2_ADM_VM_290816.mat')
ADMs = loadmat(ADMdataFile)

# Set tOffset for calcs, 3.76ps!!!
# This is because this is 2-pulse case, and will set t=0 to 2nd pulse (and matches defn. in N2 experimental paper)
tOffset = -3.76
ADMs['time'] = ADMs['time'] + tOffset

data.setADMs(ADMs = ADMs['ADM'], t=ADMs['time'].squeeze(), KQSLabels = ADMs['ADMlist'], addS = True)
data.data['ADM']['ADM']
[11]:
Show/Hide data repr Show/Hide attributes
xarray.DataArray
'ADM'
  • ADM: 4
  • t: 3691
  • (1+0j) (1+0j) ... (0.03781316571382804+4.0131898303059876e-16j)
    array([[ 1.00000000e+00+0.00000000e+00j,  1.00000000e+00+0.00000000e+00j,
             1.00000000e+00+0.00000000e+00j, ...,
             1.00000000e+00+0.00000000e+00j,  1.00000000e+00+0.00000000e+00j,
             1.00000000e+00+0.00000000e+00j],
           [-2.26243113e-17+0.00000000e+00j,  2.43430608e-08+1.04125246e-20j,
             9.80188266e-08+6.89166168e-20j, ...,
             1.05433798e-01-1.62495135e-18j,  1.05433798e-01-1.62495135e-18j,
             1.05433798e-01-1.62495135e-18j],
           [ 1.55724057e-16+0.00000000e+00j, -3.37021111e-10-6.81416260e-20j,
             1.95424253e-10-3.10513374e-19j, ...,
             8.39913132e-02-5.12795441e-17j,  8.39913132e-02-5.12795441e-17j,
             8.39913132e-02-5.12795441e-17j],
           [-7.68430227e-16+0.00000000e+00j, -1.40177466e-11+1.04987400e-19j,
             6.33419102e-10+1.74747003e-18j, ...,
             3.78131657e-02+4.01318983e-16j,  3.78131657e-02+4.01318983e-16j,
             3.78131657e-02+4.01318983e-16j]])
    • ADM
      (ADM)
      MultiIndex
      (K, Q, S)
      array([(0, 0, 0), (2, 0, 0), (4, 0, 0), (6, 0, 0)], dtype=object)
    • K
      (ADM)
      int64
      0 2 4 6
      array([0, 2, 4, 6], dtype=int64)
    • Q
      (ADM)
      int64
      0 0 0 0
      array([0, 0, 0, 0], dtype=int64)
    • S
      (ADM)
      int64
      0 0 0 0
      array([0, 0, 0, 0], dtype=int64)
    • t
      (t)
      float64
      -3.76 -3.76 -3.76 ... 10.1 10.1
      units :
      ps
      array([-3.76    , -3.7598  , -3.7596  , ..., 10.099975, 10.099975, 10.099975])
  • dataType :
    ADM
    long_name :
    Axis distribution moments
    units :
    arb
[12]:
# The ADMplot routine will show a basic line plot, note it needs keys = 'ADM' in the current implementation (otherwise will loop over all keys)
data.ADMplot(keys = 'ADM')
Dataset: ADM, ADM
../_images/fitting_PEMtk_fitting_basic_demo_030621-full_22_1.png
[13]:
# lmPlot is also handy, it minimally needs the correct keys, dataType and xDim specified

# data.lmPlot(keys = 'ADM', dataType = 'ADM', xDim = 't')  # Minimal call

data.lmPlot(keys = 'ADM', dataType = 'ADM', xDim = 't', fillna = True, logFlag = False, cmap = 'vlag')  # Set some additional options
Plotting data (No filename), pType=a, thres=0.01, with Seaborn
No handles with labels found to put in legend.
../_images/fitting_PEMtk_fitting_basic_demo_030621-full_23_2.png

Polarisation geometry/ies

This wraps ep.setPolGeoms. This defaults to (x,y,z) polarization geometries. Values are set in self.data['pol'].

Note: if this is not set, the default value will be used, which is likely not very useful for the fit!

[14]:
data.setPolGeoms()
data.data['pol']['pol']
[14]:
Show/Hide data repr Show/Hide attributes
xarray.DataArray
  • Labels: 3
  • quaternion(1, -0, 0, 0) ... quaternion(0.5, -0.5, 0.5, 0.5)
    array([quaternion(1, -0, 0, 0),
           quaternion(0.707106781186548, -0, 0.707106781186548, 0),
           quaternion(0.5, -0.5, 0.5, 0.5)], dtype=quaternion)
    • Euler
      (Labels)
      object
      (0.0, 0.0, 0.0) ... (1.5707963267948966, 1.5707963267948966, 0.0)
      array([(0.0, 0.0, 0.0), (0.0, 1.5707963267948966, 0.0),
             (1.5707963267948966, 1.5707963267948966, 0.0)], dtype=object)
    • Labels
      (Labels)
      <U18
      'z' 'x' 'y'
      array(['z', 'x', 'y'], dtype='<U18')
  • dataType :
    Euler
[15]:
# # data.setPolGeoms(eulerAngs = [[0,0,0]], labels = ['z'])
# data.setPolGeoms(eulerAngs = [0,0,0], labels = 'z')
# data.data['pol']['pol']  #.swap_dims({'Euler':'Labels'})
[16]:
# data.data['pol']['pol'] = data.data['pol']['pol'].swap_dims({'Euler':'Labels'})
# data.selOpts['pol'] = {'inds': {'Labels': 'z'}}
# data.setSubset(dataKey = 'pol', dataType = 'pol')

Subselect data

Currently handled in the class by setting self.selOpts, this allows for simple reuse of settings as required. Subselected data is set to self.data['subset'][dataType], and is the data the fitting routine will use.

[17]:
# Settings for type subselection are in selOpts[dataType]

# E.g. Matrix element sub-selection
data.selOpts['matE'] = {'thres': 0.01, 'inds': {'Type':'L', 'Eke':1.1}}
data.setSubset(dataKey = 'orb5', dataType = 'matE')  # Subselect from 'orb5' dataset, matrix elements

# Show subselected data
data.data['subset']['matE']
Subselected from dataset 'orb5', dataType 'matE': 36 from 11016 points (0.33%)
[17]:
Show/Hide data repr Show/Hide attributes
xarray.DataArray
'n2_3sg_0.1-50.1eV_A2.inp.out'
  • LM: 6
  • Sym: 2
  • mu: 3
  • it: 1
  • (nan+nanj) (nan+nanj) (nan+nanj) ... (nan+nanj) (nan+nanj)
    array([[[[        nan       +nanj],
             [        nan       +nanj],
             [        nan       +nanj]],
    
            [[        nan       +nanj],
             [        nan       +nanj],
             [ 1.1629411 -1.3536696j ]]],
    
    
           [[[        nan       +nanj],
             [-2.3170597 +1.3587362j ],
             [        nan       +nanj]],
    
            [[        nan       +nanj],
             [        nan       +nanj],
             [        nan       +nanj]]],
    
    
           [[[        nan       +nanj],
             [        nan       +nanj],
             [        nan       +nanj]],
    
            [[ 1.1629411 -1.3536696j ],
             [        nan       +nanj],
             [        nan       +nanj]]],
    
    
           [[[        nan       +nanj],
             [        nan       +nanj],
             [        nan       +nanj]],
    
            [[        nan       +nanj],
             [        nan       +nanj],
             [-0.80272541-0.01697872j]]],
    
    
           [[[        nan       +nanj],
             [ 1.1057219 -0.08717625j],
             [        nan       +nanj]],
    
            [[        nan       +nanj],
             [        nan       +nanj],
             [        nan       +nanj]]],
    
    
           [[[        nan       +nanj],
             [        nan       +nanj],
             [        nan       +nanj]],
    
            [[-0.80272541-0.01697872j],
             [        nan       +nanj],
             [        nan       +nanj]]]])
    • Ehv
      ()
      float64
      18.4
      array(18.4)
    • LM
      (LM)
      MultiIndex
      (l, m)
      array([(1, -1), (1, 0), (1, 1), (3, -1), (3, 0), (3, 1)], dtype=object)
    • l
      (LM)
      int64
      1 1 1 3 3 3
      array([1, 1, 1, 3, 3, 3], dtype=int64)
    • m
      (LM)
      int64
      -1 0 1 -1 0 1
      array([-1,  0,  1, -1,  0,  1], dtype=int64)
    • it
      (it)
      int64
      1
      array([1], dtype=int64)
    • Type
      ()
      <U1
      'L'
      array('L', dtype='<U1')
    • mu
      (mu)
      int64
      -1 0 1
      array([-1,  0,  1], dtype=int64)
    • Eke
      ()
      float64
      1.1
      array(1.1)
    • Sym
      (Sym)
      MultiIndex
      (Cont, Targ, Total)
      array([('SU', 'SG', 'SU'), ('PU', 'SG', 'PU')], dtype=object)
    • Cont
      (Sym)
      object
      'SU' 'PU'
      array(['SU', 'PU'], dtype=object)
    • Targ
      (Sym)
      object
      'SG' 'SG'
      array(['SG', 'SG'], dtype=object)
    • Total
      (Sym)
      object
      'SU' 'PU'
      array(['SU', 'PU'], dtype=object)
    • SF
      ()
      complex128
      (2.2237521+3.6277801j)
      array(2.2237521+3.6277801j)
  • E :
    0.1
    Ehv :
    15.68
    SF :
    (2.1560627+3.741674j)
    Lmax :
    11
    Cont :
    SU
    Targ :
    SG
    Total :
    SU
    QNs :
    ['m', 'l', 'mu', 'ip', 'it', 'Value']
    dataType :
    matE
    file :
    n2_3sg_0.1-50.1eV_A2.inp.out
    fileBase :
    D:\code\github\ePSproc\data\photoionization\n2_multiorb
    fileList :
    n2_3sg_0.1-50.1eV_A2.inp.out
    jobLabel :
    3sg-1
[18]:
# Tabulate the matrix elements
# Not showing as nice table for singleton case - pd.series vs. dataframe?
data.matEtoPD(keys = 'subset', xDim = 'Sym', dropna=False)

# data.data['subset']['matE'].attrs['pd'] # PD set here

*** 3sg-1
Matrix element table, threshold=None, data type=complex128.
Cont  Targ  Total  it  l  m   mu
PU    SG    PU     1   1  -1   1    1.162941-1.353670j
                           1  -1    1.162941-1.353670j
                       3  -1   1   -0.802725-0.016979j
                           1  -1   -0.802725-0.016979j
SU    SG    SU     1   1   0   0   -2.317060+1.358736j
                       3   0   0    1.105722-0.087176j
dtype: complex128
[19]:
# And for the polarisation geometries...
data.selOpts['pol'] = {'inds': {'Labels': 'z'}}
data.setSubset(dataKey = 'pol', dataType = 'pol')
Subselected from dataset 'pol', dataType 'pol': 1 from 3 points (33.33%)
[20]:
# And for the ADMs...

data.selOpts['ADM'] = {}   #{'thres': 0.01, 'inds': {'Type':'L', 'Eke':1.1}}
data.setSubset(dataKey = 'ADM', dataType = 'ADM', sliceParams = {'t':[4, 5, 4]})
data.ADMplot(keys = 'subset')
Subselected from dataset 'ADM', dataType 'ADM': 52 from 14764 points (0.35%)
Dataset: subset, ADM
../_images/fitting_PEMtk_fitting_basic_demo_030621-full_32_1.png
[21]:
# Cusomise plot with return...
# NOT YET IMPLEMENTED
# pltObj = data.ADMplot(keys = 'subset')
# pltObj
[22]:
# Plot from Xarray vs. full dataset
# data.data['subset']['ADM'].where(ADMX['K']>0).real.squeeze().plot.line(x='t');
data.data['subset']['ADM'].real.squeeze().plot.line(x='t', marker = 'x', linestyle='dashed');
data.data['ADM']['ADM'].real.squeeze().plot.line(x='t');
../_images/fitting_PEMtk_fitting_basic_demo_030621-full_34_0.png

Compute AF-\(\beta_{LM}\) and simulate data

With all the components set, some observables can be calculated. For testing, we’ll also use this to simulate an experiemental trace…

Here we’ll use self.afblmMatEfit(), which is also the main fitting routine, and essentially wraps epsproc.afblmXprod() to compute AF-\(\beta_{LM}\)s (for more details, see the ePSproc method development docs).

If called without reference data, the method returns computed AF-\(\beta_{LM}\)s based on the input subsets already created, and also a set of (product) basis functions generated - these can be examined to get a feel for the sensitivity of the geometric part of the problem, and will also be used in fitting to limit repetitive computation.

Compute AF-\(\beta_{LM}\)s

[23]:
# data.afblmMatEfit(data = None)  # OK
BetaNormX, basis = data.afblmMatEfit()  # OK, uses default polarizations & ADMs as set in data['subset']
# BetaNormX, basis = data.afblmMatEfit(ADM = data.data['subset']['ADM'])  # OK, but currently using default polarizations
# BetaNormX, basis = data.afblmMatEfit(ADM = data.data['subset']['ADM'], pol = data.data['pol']['pol'].sel(Labels=['x']))
# BetaNormX, basis = data.afblmMatEfit(ADM = data.data['subset']['ADM'], pol = data.data['pol']['pol'].sel(Labels=['x','y']))  # This fails for a single label...?
# BetaNormX, basis = data.afblmMatEfit(RX=data.data['pol']['pol'])  # This currently fails, need to check for consistency in ep.sphCalc.WDcalc()
                                                                    # - looks like set values and inputs are not consistent in this case? Not passing angs correctly, or overriding?
                                                                    # - See also recently-added sfError flag, which may cause additional problems.

AF-\(\beta_{LM}\)s

The returned objects contain the \(\beta_{LM}\) parameters as an Xarray…

[24]:
BetaNormX
[24]:
Show/Hide data repr Show/Hide attributes
xarray.DataArray
  • Labels: 1
  • t: 13
  • BLM: 15
  • 0j (1.668823121150882+5.788281730622143e-18j) ... 0j
    array([[[ 0.00000000e+00+0.00000000e+00j,
              1.66882312e+00+5.78828173e-18j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              9.25154735e-01-2.04255472e-18j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
             -1.45626494e-01-2.50162312e-18j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              9.05679279e-03+2.69005757e-19j,
              0.00000000e+00+0.00000000e+00j],
            [ 0.00000000e+00+0.00000000e+00j,
              1.65939046e+00-7.23129056e-18j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              9.28209270e-01-1.20305969e-17j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
             -1.35884180e-01+7.60995910e-18j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              7.87360160e-03-5.26232236e-19j,
              0.00000000e+00+0.00000000e+00j],
            [ 0.00000000e+00+0.00000000e+00j,
              1.60819595e+00+5.37946903e-19j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              9.44773203e-01+2.80576809e-17j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
             -8.24684513e-02-1.29331051e-18j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              1.24796733e-03+1.99111244e-19j,
              0.00000000e+00+0.00000000e+00j],
            [ 0.00000000e+00+0.00000000e+00j,
              1.52326324e+00-1.50526506e-18j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              9.72371636e-01+7.27400639e-19j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              4.93224211e-04+1.20678198e-18j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
             -5.69662854e-03-1.12782141e-19j,
              0.00000000e+00+0.00000000e+00j],
            [ 0.00000000e+00+0.00000000e+00j,
              1.43634451e+00+9.52773406e-18j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              1.00104159e+00+2.33569626e-17j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              6.53535497e-02-2.06092931e-18j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              1.13815711e-03-3.50093184e-18j,
              0.00000000e+00+0.00000000e+00j],
            [ 0.00000000e+00+0.00000000e+00j,
              1.38620024e+00+9.20243198e-19j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              1.01780024e+00-1.33980425e-18j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              9.28414140e-02+2.05936950e-18j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              1.15029998e-02-1.86868078e-18j,
              0.00000000e+00+0.00000000e+00j],
            [ 0.00000000e+00+0.00000000e+00j,
              1.39472246e+00-4.45373740e-18j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              1.01490490e+00+1.74312370e-18j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              9.03441438e-02+2.26699795e-18j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              8.28715378e-03+9.94032689e-19j,
              0.00000000e+00+0.00000000e+00j],
            [ 0.00000000e+00+0.00000000e+00j,
              1.45369338e+00+7.91745171e-19j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              9.95369088e-01+2.73033702e-17j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              5.02659008e-02+1.13177636e-18j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              9.75429570e-04-6.51669777e-19j,
              0.00000000e+00+0.00000000e+00j],
            [ 0.00000000e+00+0.00000000e+00j,
              1.53164115e+00+1.13429074e-17j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              9.69964449e-01-3.96161113e-18j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
             -2.24085809e-02-7.42900134e-18j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              5.08531171e-03-3.02631277e-18j,
              0.00000000e+00+0.00000000e+00j],
            [ 0.00000000e+00+0.00000000e+00j,
              1.59428225e+00-1.99642906e-18j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              9.49718737e-01+1.56864870e-17j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
             -8.90643638e-02-5.46842630e-18j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              1.44722472e-02+1.25126193e-18j,
              0.00000000e+00+0.00000000e+00j],
            [ 0.00000000e+00+0.00000000e+00j,
              1.62292950e+00-6.46048690e-18j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              9.40438270e-01+1.76164331e-17j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
             -1.18540934e-01-1.18827450e-18j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              1.80790278e-02+3.96263817e-18j,
              0.00000000e+00+0.00000000e+00j],
            [ 0.00000000e+00+0.00000000e+00j,
              1.61964948e+00-9.65063099e-18j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              9.41433561e-01-1.12629747e-17j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
             -1.11981456e-01+8.46145725e-18j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              1.54232880e-02+1.10625182e-18j,
              0.00000000e+00+0.00000000e+00j],
            [ 0.00000000e+00+0.00000000e+00j,
              1.59961793e+00+3.78243863e-18j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              9.47859045e-01-1.46603781e-18j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
             -8.83338758e-02-1.35167290e-18j,
              0.00000000e+00+0.00000000e+00j,
              0.00000000e+00+0.00000000e+00j,
              1.07442112e-02+1.22430637e-18j,
              0.00000000e+00+0.00000000e+00j]]])
    • Euler
      (Labels)
      object
      (0.0, 0.0, 0.0)
      array([(0.0, 0.0, 0.0)], dtype=object)
    • Labels
      (Labels)
      <U1
      'z'
      array(['z'], dtype='<U1')
    • t
      (t)
      float64
      4.018 4.096 4.175 ... 4.884 4.962
      units :
      ps
      array([4.017657, 4.096371, 4.175086, 4.253801, 4.332516, 4.41123 , 4.489945,
             4.56866 , 4.647375, 4.726089, 4.804804, 4.883519, 4.962233])
    • Ehv
      ()
      float64
      18.4
      array(18.4)
    • it
      ()
      int64
      1
      array(1, dtype=int64)
    • Type
      ()
      <U1
      'L'
      array('L', dtype='<U1')
    • Eke
      ()
      float64
      1.1
      array(1.1)
    • SF
      ()
      complex128
      (2.2237521+3.6277801j)
      array(2.2237521+3.6277801j)
    • XSraw
      (Labels, t)
      complex128
      (1.668823121150882+5.788281730622143e-18j) ... (1.5996179261474746+3.782438628664796e-18j)
      array([[1.66882312+5.78828173e-18j, 1.65939046-7.23129056e-18j,
              1.60819595+5.37946903e-19j, 1.52326324-1.50526506e-18j,
              1.43634451+9.52773406e-18j, 1.38620024+9.20243198e-19j,
              1.39472246-4.45373740e-18j, 1.45369338+7.91745171e-19j,
              1.53164115+1.13429074e-17j, 1.59428225-1.99642906e-18j,
              1.6229295 -6.46048690e-18j, 1.61964948-9.65063099e-18j,
              1.59961793+3.78243863e-18j]])
    • XSrescaled
      (Labels, t)
      complex128
      (5.915823935128086+2.0518924487134524e-17j) ... (5.670497906355172+1.3408395826381391e-17j)
      array([[5.91582394+2.05189245e-17j, 5.88238602-2.56342576e-17j,
              5.70090622+1.90697212e-18j, 5.39982759-5.33602569e-18j,
              5.09170872+3.37749379e-17j, 4.91395189+3.26217720e-18j,
              4.9441624 -1.57880880e-17j, 5.15320886+2.80666355e-18j,
              5.4295265 +4.02095600e-17j, 5.65158343-7.07715676e-18j,
              5.75313528-2.29018298e-17j, 5.74150791-3.42105961e-17j,
              5.67049791+1.34083958e-17j]])
    • XSiso
      ()
      complex128
      (5.36805661023236+0j)
      array(5.36805661+0.j)
    • BLM
      (BLM)
      MultiIndex
      (l, m)
      array([(0, -1), (0, 0), (0, 1), (2, -1), (2, 0), (2, 1), (3, -1), (3, 0),
             (3, 1), (4, -1), (4, 0), (4, 1), (6, -1), (6, 0), (6, 1)], dtype=object)
    • l
      (BLM)
      int64
      0 0 0 2 2 2 3 3 3 4 4 4 6 6 6
      array([0, 0, 0, 2, 2, 2, 3, 3, 3, 4, 4, 4, 6, 6, 6], dtype=int64)
    • m
      (BLM)
      int64
      -1 0 1 -1 0 1 -1 0 1 -1 0 1 -1 0 1
      array([-1,  0,  1, -1,  0,  1, -1,  0,  1, -1,  0,  1, -1,  0,  1], dtype=int64)
  • E :
    0.1
    Ehv :
    15.68
    SF :
    (2.1560627+3.741674j)
    Lmax :
    11
    Cont :
    SU
    Targ :
    SG
    Total :
    SU
    QNs :
    ['m', 'l', 'mu', 'ip', 'it', 'Value']
    dataType :
    BLM
    file :
    n2_3sg_0.1-50.1eV_A2.inp.out
    fileBase :
    D:\code\github\ePSproc\data\photoionization\n2_multiorb
    fileList :
    n2_3sg_0.1-50.1eV_A2.inp.out
    jobLabel :
    3sg-1
    pd :
    Cont Targ Total it l m mu PU SG PU 1 1 -1 1 1.162941-1.353670j 1 -1 1.162941-1.353670j 3 -1 1 -0.802725-0.016979j 1 -1 -0.802725-0.016979j SU SG SU 1 1 0 0 -2.317060+1.358736j 3 0 0 1.105722-0.087176j dtype: complex128
    thres :
    None
[25]:
# ep.BLMplot(BetaNormX, xDim = 't')  # SIGH, not working - issue with Euler/Labels?
# BetaNormX.sel(Labels='z').real.squeeze().plot.line(x='t');
# ep.lmPlot(BetaNormX.sel(Labels='z'), xDim='t', SFflag=False);  #, cmap='vlag');
ep.lmPlot(BetaNormX, xDim='t', SFflag=False);
# data.lmPlot()
Plotting data n2_3sg_0.1-50.1eV_A2.inp.out, pType=a, thres=0.01, with Seaborn
../_images/fitting_PEMtk_fitting_basic_demo_030621-full_41_1.png
[26]:
# Line-plot with Xarray/Matplotlib
# Note there is no filtering here, so this includes some invalid and null terms
BetaNormX.sel(Labels='z').real.squeeze().plot.line(x='t');
../_images/fitting_PEMtk_fitting_basic_demo_030621-full_42_0.png

… and the basis sets as a dictionary.

[27]:
basis.keys()
[27]:
dict_keys(['BLMtableResort', 'polProd', 'phaseConvention', 'BLMRenorm'])

Note that the basis sets here will may be useful for deeper insight into the physics, and fitting routines, and are explored in a separate notebook.

Fitting the data

In order to fit data, and extract matrix elements from an experimental case, we’ll use the lmfit library. This wraps core Scipy fitting routines with additional objects and methods, and is further wrapped for this specific class of problems in pemtkFit class we’re using here.

Set the data to fit

Here we’ll use the values calculated above as our test data. This currently needs to be set as self.data['subset']['AFBLM'] for fitting.

[28]:
# data.data['subset']['AFBLM'] = BetaNormX  # Set manually

data.setData('sim', BetaNormX)  # Set simulated data to master structure as "sim"
data.setSubset('sim','AFBLM')   # Set to 'subset' to use for fitting.

Subselected from dataset 'sim', dataType 'AFBLM': 195 from 195 points (100.00%)
[29]:
# Set basis functions
data.basis = basis

Setting up the fit parameters

In this case, we can work from the existing matrix elements to speed up parameter creation, although in practice this may need to be approached ab initio - nonetheless, the method will be the same, and the ab initio case detailed later.

[30]:
# Input set, as defined earlier
data.data['subset']['matE'].pd
[30]:
Cont  Targ  Total  it  l  m   mu
PU    SG    PU     1   1  -1   1    1.162941-1.353670j
                           1  -1    1.162941-1.353670j
                       3  -1   1   -0.802725-0.016979j
                           1  -1   -0.802725-0.016979j
SU    SG    SU     1   1   0   0   -2.317060+1.358736j
                       3   0   0    1.105722-0.087176j
dtype: complex128
[31]:
# data.setMatEFit()  # Need to fix self.subset usage
data.setMatEFit(data.data['subset']['matE'])  #, Eke=1.1) # Some hard-coded things to fix here! Now roughly working.
Set 6 complex matrix elements to 12 fitting params, see self.params for details.
name value initial value min max vary
m_PU_SG_PU_1_n1_1 1.78461575 1.784615753610107 1.0000e-04 5.00000000 True
m_PU_SG_PU_1_1_n1 1.78461575 1.784615753610107 1.0000e-04 5.00000000 True
m_PU_SG_PU_3_n1_1 0.80290495 0.802904951323892 1.0000e-04 5.00000000 True
m_PU_SG_PU_3_1_n1 0.80290495 0.802904951323892 1.0000e-04 5.00000000 True
m_SU_SG_SU_1_0_0 2.68606212 2.686062120382649 1.0000e-04 5.00000000 True
m_SU_SG_SU_3_0_0 1.10915311 1.109153108617096 1.0000e-04 5.00000000 True
p_PU_SG_PU_1_n1_1 -0.86104140 -0.8610414024232179 -3.14159265 3.14159265 False
p_PU_SG_PU_1_1_n1 -0.86104140 -0.8610414024232179 -3.14159265 3.14159265 True
p_PU_SG_PU_3_n1_1 -3.12044446 -3.1204444620772467 -3.14159265 3.14159265 True
p_PU_SG_PU_3_1_n1 -3.12044446 -3.1204444620772467 -3.14159265 3.14159265 True
p_SU_SG_SU_1_0_0 2.61122920 2.611229196458127 -3.14159265 3.14159265 True
p_SU_SG_SU_3_0_0 -0.07867828 -0.07867827542158025 -3.14159265 3.14159265 True

This sets self.params from the matrix elements, which are a set of (real) parameters for lmfit, as a Parameters object.

Note that:

  • The input matrix elements are converted to magnitude-phase form, hence there are twice the number as the input array, and labelled m or p accordingly, along with a name based on the full set of QNs/indexes set.

  • One phase is set to vary=False, which defines a reference phase. This defaults to the first phase item.

  • Min and max values are defined, by default the ranges are 1e-4<mag<5, -pi<phase<pi.

  • No relationships between the parameters are set by default (apart from the single fixed phase), but can be set manually, see section below.

[32]:
data.params
[32]:
name value initial value min max vary
m_PU_SG_PU_1_n1_1 1.78461575 1.784615753610107 1.0000e-04 5.00000000 True
m_PU_SG_PU_1_1_n1 1.78461575 1.784615753610107 1.0000e-04 5.00000000 True
m_PU_SG_PU_3_n1_1 0.80290495 0.802904951323892 1.0000e-04 5.00000000 True
m_PU_SG_PU_3_1_n1 0.80290495 0.802904951323892 1.0000e-04 5.00000000 True
m_SU_SG_SU_1_0_0 2.68606212 2.686062120382649 1.0000e-04 5.00000000 True
m_SU_SG_SU_3_0_0 1.10915311 1.109153108617096 1.0000e-04 5.00000000 True
p_PU_SG_PU_1_n1_1 -0.86104140 -0.8610414024232179 -3.14159265 3.14159265 False
p_PU_SG_PU_1_1_n1 -0.86104140 -0.8610414024232179 -3.14159265 3.14159265 True
p_PU_SG_PU_3_n1_1 -3.12044446 -3.1204444620772467 -3.14159265 3.14159265 True
p_PU_SG_PU_3_1_n1 -3.12044446 -3.1204444620772467 -3.14159265 3.14159265 True
p_SU_SG_SU_1_0_0 2.61122920 2.611229196458127 -3.14159265 3.14159265 True
p_SU_SG_SU_3_0_0 -0.07867828 -0.07867827542158025 -3.14159265 3.14159265 True

Running a fit…

With the parameters and data set, just call self.fit()!

Statistics and outputs are handled by lmfit, which includes uncertainty estimates and correlations in the fitted parameters.

[33]:
data.fit()
[34]:
# Check fit outputs
data.result
[34]:

Fit Statistics

fitting methodleastsq
# function evals13
# data points195
# variables11
chi-square 5.3662e-31
reduced chi-square 2.9164e-33
Akaike info crit.-14597.7369
Bayesian info crit.-14561.7339

Variables

name value standard error relative error initial value min max vary
m_PU_SG_PU_1_n1_1 1.78461575 1.4711e-09 (0.00%) 1.784615753610107 1.0000e-04 5.00000000 True
m_PU_SG_PU_1_1_n1 1.78461575 1.4711e-09 (0.00%) 1.784615753610107 1.0000e-04 5.00000000 True
m_PU_SG_PU_3_n1_1 0.80290495 2.6141e-09 (0.00%) 0.802904951323892 1.0000e-04 5.00000000 True
m_PU_SG_PU_3_1_n1 0.80290495 2.6141e-09 (0.00%) 0.802904951323892 1.0000e-04 5.00000000 True
m_SU_SG_SU_1_0_0 2.68606212 6.0138e-16 (0.00%) 2.686062120382649 1.0000e-04 5.00000000 True
m_SU_SG_SU_3_0_0 1.10915311 1.4590e-15 (0.00%) 1.109153108617096 1.0000e-04 5.00000000 True
p_PU_SG_PU_1_n1_1 -0.86104140 0.00000000 (0.00%) -0.8610414024232179 -3.14159265 3.14159265 False
p_PU_SG_PU_1_1_n1 -0.86104140 2.5882e-15 (0.00%) -0.8610414024232179 -3.14159265 3.14159265 True
p_PU_SG_PU_3_n1_1 -3.12044446 2.9686e-09 (0.00%) -3.1204444620772467 -3.14159265 3.14159265 True
p_PU_SG_PU_3_1_n1 -3.12044446 2.9686e-09 (0.00%) -3.1204444620772467 -3.14159265 3.14159265 True
p_SU_SG_SU_1_0_0 2.61122920 5.6788e-16 (0.00%) 2.611229196458127 -3.14159265 3.14159265 True
p_SU_SG_SU_3_0_0 -0.07867828 1.5863e-15 (0.00%) -0.07867827542158025 -3.14159265 3.14159265 True

Correlations (unreported correlations are < 0.100)

m_PU_SG_PU_1_n1_1m_PU_SG_PU_1_1_n1-1.0000
m_PU_SG_PU_3_n1_1m_PU_SG_PU_3_1_n1-1.0000
p_PU_SG_PU_3_n1_1p_PU_SG_PU_3_1_n1-1.0000
m_SU_SG_SU_1_0_0m_SU_SG_SU_3_0_0-0.9896
m_SU_SG_SU_3_0_0p_SU_SG_SU_3_0_0-0.8970
m_SU_SG_SU_1_0_0p_SU_SG_SU_3_0_00.8887
p_PU_SG_PU_1_1_n1p_SU_SG_SU_1_0_0-0.7707
p_SU_SG_SU_1_0_0p_SU_SG_SU_3_0_0-0.7543
m_SU_SG_SU_1_0_0p_SU_SG_SU_1_0_0-0.6756
m_SU_SG_SU_3_0_0p_SU_SG_SU_1_0_00.6633
m_PU_SG_PU_1_n1_1p_PU_SG_PU_3_1_n10.4627
m_PU_SG_PU_1_1_n1p_PU_SG_PU_3_1_n1-0.4627
m_PU_SG_PU_1_n1_1p_PU_SG_PU_3_n1_1-0.4627
m_PU_SG_PU_1_1_n1p_PU_SG_PU_3_n1_10.4627
m_PU_SG_PU_1_n1_1p_PU_SG_PU_1_1_n1-0.4471
m_PU_SG_PU_1_1_n1p_PU_SG_PU_1_1_n10.4471
m_PU_SG_PU_1_n1_1m_PU_SG_PU_3_1_n1-0.4414
m_PU_SG_PU_1_1_n1m_PU_SG_PU_3_1_n10.4414
m_PU_SG_PU_1_n1_1m_PU_SG_PU_3_n1_10.4414
m_PU_SG_PU_1_1_n1m_PU_SG_PU_3_n1_1-0.4414
m_PU_SG_PU_3_1_n1p_PU_SG_PU_1_1_n10.3751
m_PU_SG_PU_3_n1_1p_PU_SG_PU_1_1_n1-0.3751
p_PU_SG_PU_1_1_n1p_PU_SG_PU_3_1_n1-0.3431
p_PU_SG_PU_1_1_n1p_PU_SG_PU_3_n1_10.3431
m_PU_SG_PU_3_1_n1m_SU_SG_SU_3_0_0-0.3314
m_PU_SG_PU_3_n1_1m_SU_SG_SU_3_0_00.3314
m_PU_SG_PU_3_1_n1m_SU_SG_SU_1_0_00.3303
m_PU_SG_PU_3_n1_1m_SU_SG_SU_1_0_0-0.3303
p_PU_SG_PU_1_1_n1p_SU_SG_SU_3_0_00.2973
m_SU_SG_SU_1_0_0p_PU_SG_PU_1_1_n10.2958
m_PU_SG_PU_1_n1_1m_SU_SG_SU_3_0_00.2929
m_PU_SG_PU_1_1_n1m_SU_SG_SU_3_0_0-0.2929
m_SU_SG_SU_3_0_0p_PU_SG_PU_1_1_n1-0.2914
m_PU_SG_PU_1_n1_1m_SU_SG_SU_1_0_0-0.2913
m_PU_SG_PU_1_1_n1m_SU_SG_SU_1_0_00.2913
m_SU_SG_SU_3_0_0p_PU_SG_PU_3_1_n10.2865
m_SU_SG_SU_3_0_0p_PU_SG_PU_3_n1_1-0.2865
m_PU_SG_PU_1_n1_1p_SU_SG_SU_1_0_00.2840
m_PU_SG_PU_1_1_n1p_SU_SG_SU_1_0_0-0.2840
m_PU_SG_PU_3_1_n1p_SU_SG_SU_3_0_00.2821
m_PU_SG_PU_3_n1_1p_SU_SG_SU_3_0_0-0.2821
m_PU_SG_PU_3_1_n1p_SU_SG_SU_1_0_0-0.2818
m_PU_SG_PU_3_n1_1p_SU_SG_SU_1_0_00.2818
m_SU_SG_SU_1_0_0p_PU_SG_PU_3_1_n1-0.2814
m_SU_SG_SU_1_0_0p_PU_SG_PU_3_n1_10.2814
m_PU_SG_PU_1_n1_1p_SU_SG_SU_3_0_0-0.2200
m_PU_SG_PU_1_1_n1p_SU_SG_SU_3_0_00.2200
m_PU_SG_PU_3_1_n1p_PU_SG_PU_3_1_n1-0.2023
m_PU_SG_PU_3_n1_1p_PU_SG_PU_3_1_n10.2023
m_PU_SG_PU_3_1_n1p_PU_SG_PU_3_n1_10.2023
m_PU_SG_PU_3_n1_1p_PU_SG_PU_3_n1_1-0.2023
p_PU_SG_PU_3_1_n1p_SU_SG_SU_1_0_00.1343
p_PU_SG_PU_3_n1_1p_SU_SG_SU_1_0_0-0.1343
p_PU_SG_PU_3_1_n1p_SU_SG_SU_3_0_0-0.1313
p_PU_SG_PU_3_n1_1p_SU_SG_SU_3_0_00.1313

Results vs. data can be (crudely) plotted with self.BLMfitPlot() (better plotting routines to follow!).

This will plot all results by default, vs. the subset data used for the fitting routine inputs.

[35]:
# Plot data subset (--x) plus fit (solid lines)
data.BLMfitPlot()
Dataset: subset, AFBLM
Dataset: 0, AFBLM
../_images/fitting_PEMtk_fitting_basic_demo_030621-full_59_1.png
[36]:
# Fit results are currently added to the main data dict by an index number
data.data.keys()
[36]:
dict_keys(['orb6', 'orb5', 'ADM', 'pol', 'subset', 'sim', 0])
[37]:
# Plot results with lmPlot wrapper - this also defaults to show data subset + most recent fit results
data.lmPlotFit()
Plotting data n2_3sg_0.1-50.1eV_A2.inp.out, pType=a, thres=0.01, with Seaborn
Plotting data (No filename), pType=a, thres=0.01, with Seaborn
../_images/fitting_PEMtk_fitting_basic_demo_030621-full_61_1.png
../_images/fitting_PEMtk_fitting_basic_demo_030621-full_61_2.png

Fitting with randomised parameter inputs

Here we might expect some variation in the results, depending on various properites (ionizing channel, dataset size etc.), and also for the fitting to take a bit longer (see benchmarks later).

TODO:

  • More careful analysis routines here, run vs. number of input points & test fiedelity.

  • Parallelize.

  • Variation over runs/general statistical analysis.

[38]:
# Basic randomize routine, [0,1] interval
data.randomizeParams()
[39]:
data.fit()
data.result
[39]:

Fit Statistics

fitting methodleastsq
# function evals687
# data points195
# variables11
chi-square 0.00178259
reduced chi-square 9.6880e-06
Akaike info crit.-2240.52413
Bayesian info crit.-2204.52113

Variables

name value standard error relative error initial value min max vary
m_PU_SG_PU_1_n1_1 1.0000e-04 0.02878885 (28788.85%) 0.02129591573558398 1.0000e-04 5.00000000 True
m_PU_SG_PU_1_1_n1 2.65953073 0.09359366 (3.52%) 0.5003341770940677 1.0000e-04 5.00000000 True
m_PU_SG_PU_3_n1_1 1.77943189 0.83227191 (46.77%) 0.6217616830375992 1.0000e-04 5.00000000 True
m_PU_SG_PU_3_1_n1 2.13431007 0.99444910 (46.59%) 0.5114874598511951 1.0000e-04 5.00000000 True
m_SU_SG_SU_1_0_0 2.71674517 0.04721051 (1.74%) 0.3903702867415594 1.0000e-04 5.00000000 True
m_SU_SG_SU_3_0_0 1.05873885 0.12108910 (11.44%) 0.1518467511430327 1.0000e-04 5.00000000 True
p_PU_SG_PU_1_n1_1 -0.86104140 0.00000000 (0.00%) -0.8610414024232179 -3.14159265 3.14159265 False
p_PU_SG_PU_1_1_n1 1.14805770 0.06955446 (6.06%) 0.045436187967926034 -3.14159265 3.14159265 True
p_PU_SG_PU_3_n1_1 -3.14159265 0.02623102 (0.83%) 0.4122447260821296 -3.14159265 3.14159265 True
p_PU_SG_PU_3_1_n1 3.14159265 0.02336103 (0.74%) 0.09510612061972745 -3.14159265 3.14159265 True
p_SU_SG_SU_1_0_0 -2.02120566 0.05776807 (2.86%) 0.30505876118085085 -3.14159265 3.14159265 True
p_SU_SG_SU_3_0_0 1.77009004 0.15150596 (8.56%) 0.8346487505349893 -3.14159265 3.14159265 True

Correlations (unreported correlations are < 0.100)

m_PU_SG_PU_3_n1_1m_PU_SG_PU_3_1_n1-0.9995
m_SU_SG_SU_1_0_0m_SU_SG_SU_3_0_0-0.9941
m_SU_SG_SU_3_0_0p_SU_SG_SU_3_0_00.9395
m_SU_SG_SU_1_0_0p_SU_SG_SU_3_0_0-0.9376
m_PU_SG_PU_1_1_n1m_PU_SG_PU_3_n1_10.8858
m_PU_SG_PU_1_1_n1m_PU_SG_PU_3_1_n1-0.8829
p_PU_SG_PU_1_1_n1p_PU_SG_PU_3_1_n1-0.7629
m_PU_SG_PU_3_1_n1m_SU_SG_SU_3_0_00.7511
m_PU_SG_PU_3_1_n1m_SU_SG_SU_1_0_0-0.7489
m_PU_SG_PU_3_n1_1m_SU_SG_SU_3_0_0-0.7485
p_PU_SG_PU_3_1_n1p_SU_SG_SU_1_0_0-0.7475
m_PU_SG_PU_3_n1_1m_SU_SG_SU_1_0_00.7465
m_PU_SG_PU_3_1_n1p_SU_SG_SU_3_0_00.6694
m_PU_SG_PU_3_n1_1p_SU_SG_SU_3_0_0-0.6644
m_PU_SG_PU_3_n1_1p_PU_SG_PU_3_1_n1-0.5807
m_PU_SG_PU_3_1_n1p_PU_SG_PU_3_1_n10.5726
m_PU_SG_PU_1_1_n1p_PU_SG_PU_1_1_n10.5688
m_PU_SG_PU_1_1_n1p_PU_SG_PU_3_1_n1-0.5174
p_PU_SG_PU_1_1_n1p_SU_SG_SU_1_0_00.5150
p_PU_SG_PU_1_1_n1p_PU_SG_PU_3_n1_10.4827
m_PU_SG_PU_1_1_n1m_SU_SG_SU_3_0_0-0.4446
p_PU_SG_PU_3_n1_1p_SU_SG_SU_1_0_00.4399
m_PU_SG_PU_1_1_n1m_SU_SG_SU_1_0_00.4394
m_PU_SG_PU_1_n1_1p_PU_SG_PU_1_1_n10.4066
m_PU_SG_PU_1_1_n1p_SU_SG_SU_3_0_0-0.3893
p_SU_SG_SU_1_0_0p_SU_SG_SU_3_0_00.3701
p_PU_SG_PU_3_n1_1p_PU_SG_PU_3_1_n1-0.3492
m_PU_SG_PU_3_n1_1p_PU_SG_PU_1_1_n10.3439
m_PU_SG_PU_3_1_n1p_PU_SG_PU_1_1_n1-0.3308
p_PU_SG_PU_1_1_n1p_SU_SG_SU_3_0_00.3272
m_PU_SG_PU_1_n1_1p_PU_SG_PU_3_1_n1-0.2833
m_SU_SG_SU_3_0_0p_PU_SG_PU_3_1_n10.2446
m_SU_SG_SU_1_0_0p_PU_SG_PU_3_1_n1-0.2388
m_PU_SG_PU_1_n1_1p_SU_SG_SU_1_0_00.1740
m_PU_SG_PU_1_n1_1m_PU_SG_PU_1_1_n10.1679
m_PU_SG_PU_1_n1_1p_PU_SG_PU_3_n1_1-0.1568
m_PU_SG_PU_1_1_n1p_PU_SG_PU_3_n1_10.1397
m_PU_SG_PU_1_1_n1p_SU_SG_SU_1_0_0-0.1345
m_SU_SG_SU_3_0_0p_PU_SG_PU_3_n1_1-0.1326
m_SU_SG_SU_1_0_0p_PU_SG_PU_3_n1_10.1259
m_SU_SG_SU_1_0_0p_PU_SG_PU_1_1_n1-0.1236
m_PU_SG_PU_3_1_n1p_PU_SG_PU_3_n1_1-0.1203
m_PU_SG_PU_1_n1_1m_PU_SG_PU_3_n1_10.1195
m_SU_SG_SU_3_0_0p_PU_SG_PU_1_1_n10.1193
p_PU_SG_PU_3_n1_1p_SU_SG_SU_3_0_00.1174
m_PU_SG_PU_3_n1_1p_PU_SG_PU_3_n1_10.1150
m_SU_SG_SU_3_0_0p_SU_SG_SU_1_0_00.1076
m_SU_SG_SU_1_0_0p_SU_SG_SU_1_0_0-0.1062
[40]:
# Note that if keys is not set, BLMfitPlot will show all fit run results.
data.BLMfitPlot()
Dataset: subset, AFBLM
Dataset: 1, AFBLM
../_images/fitting_PEMtk_fitting_basic_demo_030621-full_65_1.png

Timing for the test fit (next cell) - note this is currently just running with defaults, single core only, and we’re not checking for good fiedelity here. Fits take between 1 and 12 mins, approximately.

Results may vary depending on the inputs… for the current test case:

  • 1st test The slowest run took 12.88 times longer than the fastest. This could mean that an intermediate result is being cached. 50.7 s ± 1min 2s per loop (mean ± std. dev. of 7 runs, 1 loop each)

  • 2nd test The slowest run took 105.38 times longer than the fastest. This could mean that an intermediate result is being cached. 4min 13s ± 7min 6s per loop (mean ± std. dev. of 7 runs, 1 loop each)

TODO:

  • More careful testing & benchmarks.

  • Timing vs. input dataset size.

  • Fitting statistics over fits (fidelity vs. fit vs. dataset size etc.).

[41]:
# %%timeit

# data.randomizeParams()
# data.fit()

Setting parameter relations/constraints

If we know that some of the matrix elements (parameters) are related, this can be set using contraints on the parameters.

In this test case, we know some terms are equal… this should speed up the fitting, and also improve fiedelity.

TODO: automatic relation setting in cases derived from computational matrix elements.

[42]:
# With constraints
# Set param constraints as dict
paramsCons = {}
paramsCons['m_PU_SG_PU_1_n1_1'] = 'm_PU_SG_PU_1_1_n1'
paramsCons['p_PU_SG_PU_1_n1_1'] = 'p_PU_SG_PU_1_1_n1'

paramsCons['m_PU_SG_PU_3_n1_1'] = 'm_PU_SG_PU_3_1_n1'
paramsCons['p_PU_SG_PU_3_n1_1'] = 'p_PU_SG_PU_3_1_n1'

data.setMatEFit(paramsCons = paramsCons)
Set 6 complex matrix elements to 12 fitting params, see self.params for details.
name value initial value min max vary expression
m_PU_SG_PU_1_n1_1 1.78461575 1.784615753610107 1.0000e-04 5.00000000 False m_PU_SG_PU_1_1_n1
m_PU_SG_PU_1_1_n1 1.78461575 1.784615753610107 1.0000e-04 5.00000000 True
m_PU_SG_PU_3_n1_1 0.80290495 0.802904951323892 1.0000e-04 5.00000000 False m_PU_SG_PU_3_1_n1
m_PU_SG_PU_3_1_n1 0.80290495 0.802904951323892 1.0000e-04 5.00000000 True
m_SU_SG_SU_1_0_0 2.68606212 2.686062120382649 1.0000e-04 5.00000000 True
m_SU_SG_SU_3_0_0 1.10915311 1.109153108617096 1.0000e-04 5.00000000 True
p_PU_SG_PU_1_n1_1 -0.86104140 -0.8610414024232179 -3.14159265 3.14159265 False p_PU_SG_PU_1_1_n1
p_PU_SG_PU_1_1_n1 -0.86104140 -0.8610414024232179 -3.14159265 3.14159265 True
p_PU_SG_PU_3_n1_1 -3.12044446 -3.1204444620772467 -3.14159265 3.14159265 False p_PU_SG_PU_3_1_n1
p_PU_SG_PU_3_1_n1 -3.12044446 -3.1204444620772467 -3.14159265 3.14159265 True
p_SU_SG_SU_1_0_0 2.61122920 2.611229196458127 -3.14159265 3.14159265 True
p_SU_SG_SU_3_0_0 -0.07867828 -0.07867827542158025 -3.14159265 3.14159265 True
[43]:
data.randomizeParams()
data.fit()
data.result
[43]:

Fit Statistics

fitting methodleastsq
# function evals161
# data points195
# variables8
chi-square 1.1026e-17
reduced chi-square 5.8961e-20
Akaike info crit.-8626.26451
Bayesian info crit.-8600.08051

Variables

name value standard error relative error initial value min max vary expression
m_PU_SG_PU_1_n1_1 1.78461575 8.1194e-10 (0.00%) 0.6840385724334991 1.0000e-04 5.00000000 False m_PU_SG_PU_1_1_n1
m_PU_SG_PU_1_1_n1 1.78461575 8.1194e-10 (0.00%) 0.6840385724334991 1.0000e-04 5.00000000 True
m_PU_SG_PU_3_n1_1 0.80290495 1.8017e-09 (0.00%) 0.5216275685614492 1.0000e-04 5.00000000 False m_PU_SG_PU_3_1_n1
m_PU_SG_PU_3_1_n1 0.80290495 1.8017e-09 (0.00%) 0.5216275685614492 1.0000e-04 5.00000000 True
m_SU_SG_SU_1_0_0 2.68606212 1.7019e-09 (0.00%) 0.22227352468309725 1.0000e-04 5.00000000 True
m_SU_SG_SU_3_0_0 1.10915311 4.1387e-09 (0.00%) 0.6562471595616396 1.0000e-04 5.00000000 True
p_PU_SG_PU_1_n1_1 -0.65962973 3.2020e-09 (0.00%) 0.6826881478354998 -3.14159265 3.14159265 False p_PU_SG_PU_1_1_n1
p_PU_SG_PU_1_1_n1 -0.65962973 3.2020e-09 (0.00%) 0.6826881478354998 -3.14159265 3.14159265 True
p_PU_SG_PU_3_n1_1 1.59977333 4.1546e-09 (0.00%) 0.8003201107667578 -3.14159265 3.14159265 False p_PU_SG_PU_3_1_n1
p_PU_SG_PU_3_1_n1 1.59977333 4.1546e-09 (0.00%) 0.8003201107667578 -3.14159265 3.14159265 True
p_SU_SG_SU_1_0_0 2.15128498 3.2382e-09 (0.00%) 0.801557156402618 -3.14159265 3.14159265 True
p_SU_SG_SU_3_0_0 -1.44199286 3.9100e-09 (0.00%) 0.901150896749868 -3.14159265 3.14159265 True

Correlations (unreported correlations are < 0.100)

m_SU_SG_SU_1_0_0m_SU_SG_SU_3_0_0-0.9736
m_PU_SG_PU_1_1_n1m_PU_SG_PU_3_1_n1-0.9157
m_SU_SG_SU_1_0_0p_SU_SG_SU_1_0_0-0.7969
p_PU_SG_PU_1_1_n1p_SU_SG_SU_1_0_00.7919
m_SU_SG_SU_3_0_0p_SU_SG_SU_1_0_00.7866
p_PU_SG_PU_1_1_n1p_PU_SG_PU_3_1_n10.7545
p_PU_SG_PU_3_1_n1p_SU_SG_SU_3_0_00.7212
m_PU_SG_PU_3_1_n1p_PU_SG_PU_3_1_n10.6668
m_PU_SG_PU_1_1_n1p_PU_SG_PU_3_1_n1-0.6463
m_PU_SG_PU_1_1_n1p_PU_SG_PU_1_1_n1-0.6318
m_PU_SG_PU_3_1_n1p_PU_SG_PU_1_1_n10.6010
m_SU_SG_SU_3_0_0p_SU_SG_SU_3_0_0-0.5259
m_SU_SG_SU_1_0_0p_SU_SG_SU_3_0_00.4990
m_PU_SG_PU_3_1_n1m_SU_SG_SU_3_0_0-0.3999
m_PU_SG_PU_3_1_n1p_SU_SG_SU_3_0_00.3943
m_PU_SG_PU_1_1_n1m_SU_SG_SU_1_0_0-0.3921
p_PU_SG_PU_3_1_n1p_SU_SG_SU_1_0_00.3847
m_PU_SG_PU_3_1_n1m_SU_SG_SU_1_0_00.3579
m_PU_SG_PU_1_1_n1m_SU_SG_SU_3_0_00.3445
m_PU_SG_PU_1_1_n1p_SU_SG_SU_3_0_0-0.3393
m_SU_SG_SU_3_0_0p_PU_SG_PU_1_1_n10.3276
m_SU_SG_SU_1_0_0p_PU_SG_PU_1_1_n1-0.3107
p_PU_SG_PU_1_1_n1p_SU_SG_SU_3_0_00.2764
m_SU_SG_SU_3_0_0p_PU_SG_PU_3_1_n1-0.1626
m_SU_SG_SU_1_0_0p_PU_SG_PU_3_1_n10.1538
[44]:
# Note that if keys is not set, BLMfitPlot will show only most recent fit run results.
data.BLMfitPlot()
Dataset: subset, AFBLM
Dataset: 2, AFBLM
../_images/fitting_PEMtk_fitting_basic_demo_030621-full_71_1.png

Timing for the test fit (next cell, same method as earlier) - note this is currently just running with defaults, single core only, and we’re not checking for good fiedelity here. Note, also, that these results are in the ~10s range, compared to ~many minutes for the unconstrained case.

  • 1st test 10.7 s ± 2.86 s per loop (mean ± std. dev. of 7 runs, 1 loop each)

  • 2nd test 7.9 s ± 3.5 s per loop (mean ± std. dev. of 7 runs, 1 loop each)

[45]:
%%timeit

data.randomizeParams()
data.fit()
The slowest run took 7.96 times longer than the fastest. This could mean that an intermediate result is being cached.
17.8 s ± 17.1 s per loop (mean ± std. dev. of 7 runs, 1 loop each)
[46]:
# Plot multiple fit sets
# Note that if keys = 'all' is set, BLMfitPlot will show ALL fit run results.
# TODO: use Seaborn/HV for better plotting here!
data.BLMfitPlot(keys = 'all')
Dataset: subset, AFBLM
Dataset: 0, AFBLM
Dataset: 1, AFBLM
Dataset: 2, AFBLM
Dataset: 3, AFBLM
Dataset: 4, AFBLM
Dataset: 5, AFBLM
Dataset: 6, AFBLM
Dataset: 7, AFBLM
Dataset: 8, AFBLM
Dataset: 9, AFBLM
Dataset: 10, AFBLM
../_images/fitting_PEMtk_fitting_basic_demo_030621-full_74_1.png

Versions

[47]:
import scooby
scooby.Report(additional=['epsproc', 'pemtk', 'xarray', 'jupyter'])
[47]:
Sat Jun 05 12:49:56 2021 Eastern Daylight Time
OS Windows CPU(s) 32 Machine AMD64
Architecture 64bit RAM 63.9 GB Environment Jupyter
Python 3.7.3 (default, Apr 24 2019, 15:29:51) [MSC v.1915 64 bit (AMD64)]
epsproc 1.3.0-dev pemtk 0.0.1 xarray 0.15.0
jupyter Version unknown numpy 1.19.2 scipy 1.3.0
IPython 7.12.0 matplotlib 3.3.1 scooby 0.5.6
Intel(R) Math Kernel Library Version 2020.0.0 Product Build 20191125 for Intel(R) 64 architecture applications
[48]:
# 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
* dev
  master
  numba-tests
da12376cc36f640d8974f5ce2c121be3d391caab
[49]:
# Check current remote commits
!git ls-remote --heads git://github.com/phockett/ePSproc
# !git ls-remote --heads git://github.com/phockett/epsman
da12376cc36f640d8974f5ce2c121be3d391caab        refs/heads/dev
82d12cf35b19882d4e9a2cde3d4009fe679cfaee        refs/heads/master
69cd89ce5bc0ad6d465a4bd8df6fba15d3fd1aee        refs/heads/numba-tests
ea30878c842f09d525fbf39fa269fa2302a13b57        refs/heads/revert-9-master
[50]:
# Check current Git commit for local PEMtk version
import pemtk
!git -C {Path(pemtk.__file__).parent} branch
!git -C {Path(pemtk.__file__).parent} log --format="%H" -n 1
* master
fca744ca18b98ecd49fbf17fc79247ebee6b9c3a
[51]:
# Check current remote commits
!git ls-remote --heads git://github.com/phockett/PEMtk
# !git ls-remote --heads git://github.com/phockett/epsman
fca744ca18b98ecd49fbf17fc79247ebee6b9c3a        refs/heads/master