9. Case study: Generalised bootstrapping for a linear heteronuclear scattering system, \(OCS~(C_{\infty v})\)#

In this chapter, the full code and analysis details of the case study for \(OCS\) are given, including obtaining required data, running fits and analysis routines. For more details on the routines, see the PEMtk documentation [20]; for the analysis see particularly the fit fidelity and analysis page, and molecular frame analysis data processing page (full analysis for Ref. [3], illustrating the \(N_2\) case).

9.1. General setup#

In the following code cells (see source notebooks for full details) the general setup routines (as per the outline in Chpt. 7 are executed via a configuration script with presets for the case studies herein.

Additionally, the routines will either run fits, or load existing data if available. Since fitting can be computationally demanding, it is, in general, recommended to approach large fitting problems carefully.

# Configure settings for case study

# Set case study by name
fitSystem='OCS'
fitStem=f"fit_withNoise_orb13"

# Add noise?
addNoise = 'y'
mu, sigma = 0, 0.05  # Up to approx 10% noise (+/- 0.05)

# Batching - number of fits to run between data dumps
batchSize = 10

# Total fits to run
nMax = 10
Hide code cell content 
# Run default config - may need to set full path here

%run '../scripts/setup_notebook_caseStudies_Mod-300723.py'   # Test version with different figure options.
# %run '../scripts/setup_notebook.py'

# Set outputs for notebook or PDF (skips Holoviews plots unless glued)
# Note this is set to default 'pl' in script above
if buildEnv == 'pdf':
    paramPlotBackend = 'sns'    # For category plots with paramPlot
else:
    paramPlotBackend = 'hv'
    
# plotBackend = 'sns'    # For category plots with paramPlot

*** Setting up notebook with standard Quantum Metrology Vol. 3 imports...
For more details see https://pemtk.readthedocs.io/en/latest/fitting/PEMtk_fitting_basic_demo_030621-full.html
To use local source code, pass the parent path to this script at run time, e.g. "setup_fit_demo ~/github"

*** Running: 2023-12-07 10:49:51
Working dir: /home/jovyan/jake-home/buildTmp/_latest_build/html/doc-source/part2
Build env: html
None

* Loading packages...
* sparse not found, sparse matrix forms not available. 
* natsort not found, some sorting functions not available. 
* Setting plotter defaults with epsproc.basicPlotters.setPlotters(). Run directly to modify, or change options in local env.
* Set Holoviews with bokeh.
* pyevtk not found, VTK export not available. 
* Set Holoviews with bokeh.
Jupyter Book      : 0.15.1
External ToC      : 0.3.1
MyST-Parser       : 0.18.1
MyST-NB           : 0.17.2
Sphinx Book Theme : 1.0.1
Jupyter-Cache     : 0.6.1
NbClient          : 0.7.4

# Pull data from web (OCS case)
from epsproc.util.io import getFilesFromGithub

# Set dataName (will be used as download subdir)
dataName = 'OCSfitting'
# OCS matrix elements
fDictMatE, fAllMatE = getFilesFromGithub(subpath='data/photoionization/OCS_multiorb', dataName=dataName, ref='dev')
# OCS alignment data
fDictADM, fAllADM = getFilesFromGithub(subpath='data/alignment/OCS_ADMs_28K_VM_070722', dataName=dataName, ref='dev')

Querying URL: https://api.github.com/repos/phockett/epsproc/contents/data/photoionization/OCS_multiorb?ref=dev
Local file /home/jovyan/jake-home/buildTmp/_latest_build/html/doc-source/part2/OCSfitting/OCS_survey.orb10_E0.1_2.0_30.1eV.inp.out already exists
Local file /home/jovyan/jake-home/buildTmp/_latest_build/html/doc-source/part2/OCSfitting/OCS_survey.orb10_E1.1_2.0_31.1eV.inp.out already exists
Local file /home/jovyan/jake-home/buildTmp/_latest_build/html/doc-source/part2/OCSfitting/OCS_survey.orb11_E0.1_2.0_30.1eV.inp.out already exists
Local file /home/jovyan/jake-home/buildTmp/_latest_build/html/doc-source/part2/OCSfitting/OCS_survey.orb11_E1.1_2.0_31.1eV.inp.out already exists
Querying URL: https://api.github.com/repos/phockett/epsproc/contents/data/alignment/OCS_ADMs_28K_VM_070722?ref=dev
Local file /home/jovyan/jake-home/buildTmp/_latest_build/html/doc-source/part2/OCSfitting/A20_300fs_4p2TW_28K.dat already exists
Local file /home/jovyan/jake-home/buildTmp/_latest_build/html/doc-source/part2/OCSfitting/A40_300fs_4p2TW_28K.dat already exists
Local file /home/jovyan/jake-home/buildTmp/_latest_build/html/doc-source/part2/OCSfitting/A60_300fs_4p2TW_28K.dat already exists
Local file /home/jovyan/jake-home/buildTmp/_latest_build/html/doc-source/part2/OCSfitting/c2t_300fs_4p2TW_28K.dat already exists
Local file /home/jovyan/jake-home/buildTmp/_latest_build/html/doc-source/part2/OCSfitting/time_300fs_4p2TW_28K.dat already exists

# Fitting setup including data generation and parameter creation

# Set datapath, 
dataPath = Path(Path.cwd(),dataName)

# Run general config script with dataPath set above
%run "../scripts/setup_fit_case-studies_270723.py" -d {dataPath} -a {dataPath} -c {fitSystem} -n {addNoise} --sigma {sigma} -a3D 'y'
Hide code cell output 
*** Setting up demo fitting workspace and main `data` class object...
Script: QM3 case studies
For more details see https://pemtk.readthedocs.io/en/latest/fitting/PEMtk_fitting_basic_demo_030621-full.html
To use local source code, pass the parent path to this script at run time, e.g. "setup_fit_demo ~/github"


* Loading packages...

* Loading demo matrix element data from /home/jovyan/jake-home/buildTmp/_latest_build/html/doc-source/part2/OCSfitting
*** Warning: Missing records, expected 64, found 48.
*** Warning: Found 16 blank sets of matrix elements, symmetries ['A2']

*** Job subset details
Key: subset
No 'job' info set for self.data[subset].

*** Job orb13 details
Key: orb13
Dir /home/jovyan/jake-home/buildTmp/_latest_build/html/doc-source/part2/OCSfitting, 1 file(s).
{   'batch': 'ePS OCS, batch OCS_survey, orbital orb10',
    'event': 'orb 10 (A1/S) ionization, basic survey run.',
    'orbE': -17.32548962282056,
    'orbLabel': 'A1/S'}

*** Job orb14 details
Key: orb14
Dir /home/jovyan/jake-home/buildTmp/_latest_build/html/doc-source/part2/OCSfitting, 1 file(s).
{   'batch': 'ePS OCS, batch OCS_survey, orbital orb11',
    'event': ' orb 11 ionization, basic survey run.',
    'orbE': -11.35803261907539,
    'orbLabel': '# OCS, orb 11 ionization, basic survey run.'}
*** Molecular structure

*** Molecular orbital list (from ePS output file)
EH = Energy (Hartrees), E = Energy (eV), NOrbGrp, OrbGrp, GrpDegen = degeneracies and corresponding orbital numbering by group in ePS, NormInt = single centre expansion convergence (should be ~1.0).

*** Warning: some orbital convergences outside single-center expansion convergence tolerance (0.01):
[[1.         0.66056934]
 [2.         0.96318872]
 [4.         0.95118653]
 [6.         0.9785498 ]
 [7.         0.9785498 ]]

*** Job subset details
Key: subset
No 'job' info set for self.data[subset].

*** Job orb13 details
Key: orb13
Dir /home/jovyan/jake-home/buildTmp/_latest_build/html/doc-source/part2/OCSfitting, 1 file(s).
{   'batch': 'ePS OCS, batch OCS_survey, orbital orb10',
    'event': 'orb 10 (A1/S) ionization, basic survey run.',
    'orbE': -17.32548962282056,
    'orbLabel': 'A1/S'}

*** Job orb14 details
Key: orb14
Dir /home/jovyan/jake-home/buildTmp/_latest_build/html/doc-source/part2/OCSfitting, 1 file(s).
{   'batch': 'ePS OCS, batch OCS_survey, orbital orb11',
    'event': ' orb 11 ionization, basic survey run.',
    'orbE': -11.35803261907539,
    'orbLabel': '# OCS, orb 11 ionization, basic survey run.'}
*** Running with default OCS settings.


* Loading demo ADM data from dir /home/jovyan/jake-home/buildTmp/_latest_build/html/doc-source/part2/OCSfitting...

*** Scanning dir
/home/jovyan/jake-home/buildTmp/_latest_build/html/doc-source/part2/OCSfitting
Found 5 .dat file(s)

*** Adding 3D test ADMs...

* Subselecting data...
*** Setting for OCS case study.
Subselected from dataset 'orb13', dataType 'matE': 132 from 7104 points (1.86%)
Subselected from dataset 'pol', dataType 'pol': 1 from 3 points (33.33%)
ADMs: Selecting 51 points from 1775

* Calculating AF-BLMs...

Running for 51 t-points
Subselected from dataset 'sim', dataType 'AFBLM': 2295 from 2295 points (100.00%)

*Setting  up fit parameters (with constraints)...
Set 22 complex matrix elements to 44 fitting params, see self.params for details.
Auto-setting parameters.
Basis set size = 44 params. Fitting with 15 magnitudes and 14 phases floated.

*** Adding Gaussian noise, mu=0, sigma=0.05


*** Setup demo fitting workspace OK.
Dataset: subset, AFBLM
Dataset: sim, AFBLM
../_images/0d6e3a301d636dc771ff0daa694003a49564a825d6b26db4ebf6ce1c15b78a72.png
props Sym EH Occ E NOrbGrp OrbGrp GrpDegen NormInt
orb
1 S -91.990 2.0 -2503.178 1.0 1.0 1.0 0.661
2 S -20.675 2.0 -562.590 1.0 2.0 1.0 0.963
3 S -11.445 2.0 -311.437 1.0 3.0 1.0 1.000
4 S -8.994 2.0 -244.728 1.0 4.0 1.0 0.951
5 S -6.676 2.0 -181.674 1.0 5.0 1.0 0.994
6 P -6.672 2.0 -181.557 1.0 6.0 2.0 0.979
7 P -6.672 2.0 -181.557 2.0 6.0 2.0 0.979
8 S -1.537 2.0 -41.832 1.0 7.0 1.0 0.998
9 S -1.088 2.0 -29.595 1.0 8.0 1.0 0.997
10 S -0.791 2.0 -21.513 1.0 9.0 1.0 0.999
11 P -0.674 2.0 -18.349 1.0 10.0 2.0 1.000
12 P -0.674 2.0 -18.349 2.0 10.0 2.0 1.000
13 S -0.637 2.0 -17.325 1.0 11.0 1.0 0.999
14 P -0.417 2.0 -11.358 1.0 12.0 2.0 0.999
15 P -0.417 2.0 -11.358 2.0 12.0 2.0 0.999
../_images/16ef6a0ecfb5d86dc4f2a2468f6518a1caaa455354599ae9cc72e523a3a3408b.png

9.2. Load existing fit data or run fits#

Note that running fits may be quite time-consuming and computationally intensive, depending on the size of the size of the problem. The default case here will run a small batch for testing if there is no existing data found on the dataPath, otherwise the data is loaded for analysis.

# Look for existing Pickle files on path
dataFiles = list(dataPath.expanduser().glob('*.pickle'))

if not dataFiles:
    print("No data found, executing minimal fitting run...")
    
    # Run fit batch - single
    # data.multiFit(nRange = [n,n+batchSize-1], num_workers=batchSize)

    # Run fit batches with checkpoint files
    for n in np.arange(0,nMax,batchSize):
        print(f'*** Running batch [{n},{n+batchSize-1}], {dt.now().strftime("%d%m%y_%H-%M-%S")}')

        # Run fit batch
        data.multiFit(nRange = [n,n+batchSize-1], num_workers=batchSize)

        # Dump data so far
        data.writeFitData(outStem=f"{fitSystem}_{n+batchSize-1}_{fitStem}")
        
        print(f'Finished batch [{n},{n+batchSize-1}], {dt.now().strftime("%d%m%y_%H-%M-%S")}')
        print(f'Written to file {fitSystem}_{n+batchSize-1}_{fitStem}')

else:
    dataFileIn = dataFiles[-1]   # Add index to select file, although loadFitData will concat multiple files
                                    # Note that concat currently only works for fixed batch sizes however.
    print(f"Set dataFiles: {dataFileIn}")
    data.loadFitData(fList=dataFileIn, dataPath=dataPath)   #.expanduser())
    
    data.BLMfitPlot(keys=['subset','sim'])
    
    # # Check ADMs
    # data.data['subset']['ADM'].unstack().where(data.data['subset']['ADM'].unstack().K>0) \
    #     .real.hvplot.line(x='t').overlay(['K','Q','S'])
    

Set dataFiles: /home/jovyan/jake-home/buildTmp/_latest_build/html/doc-source/part2/OCSfitting/OCS_999_fit_3D-test_withNoise_orb13_200723_05-47-50.pickle
Read data from /home/jovyan/jake-home/buildTmp/_latest_build/html/doc-source/part2/OCSfitting/OCS_999_fit_3D-test_withNoise_orb13_200723_05-47-50.pickle with pickle.
Dataset: subset, AFBLM
Dataset: sim, AFBLM
../_images/46a81ce37ff372ec048d090a0fd63634a441e75dfa992b4e043b02722b9597da.png 
# Check ADMs
# Basic plotter
data.ADMplot(keys = 'subset')

Dataset: subset, ADM
../_images/f535bb88d4225dcd57b58af0f1eb3407102cf1b7c47fcbd171975d898401463d.png
Hide code cell content 
# Check ADMs
# Holoviews
data.data['subset']['ADM'].unstack().where(data.data['subset']['ADM'].unstack().K>0) \
    .real.hvplot.line(x='t').overlay(['K','Q','S'])

# Fits appear as integer indexed items in the main data structure.
data.data.keys()
Hide code cell output 
dict_keys(['subset', 'orb13', 'orb14', 'ADM', 'pol', 'sim', 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223, 224, 225, 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, 236, 237, 238, 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259, 260, 261, 262, 263, 264, 265, 266, 267, 268, 269, 270, 271, 272, 273, 274, 275, 276, 277, 278, 280, 281, 282, 283, 284, 285, 286, 287, 288, 289, 290, 291, 292, 293, 294, 295, 296, 297, 298, 299, 300, 301, 302, 303, 304, 305, 306, 307, 308, 309, 310, 311, 312, 313, 314, 315, 316, 317, 318, 320, 321, 322, 323, 324, 325, 326, 327, 328, 329, 330, 331, 332, 333, 334, 335, 336, 337, 338, 339, 340, 341, 342, 343, 344, 345, 346, 347, 348, 349, 350, 351, 352, 353, 354, 355, 356, 357, 358, 360, 361, 362, 363, 364, 365, 366, 367, 368, 369, 370, 371, 372, 373, 374, 375, 376, 377, 378, 379, 380, 381, 382, 383, 384, 385, 386, 387, 388, 389, 390, 391, 392, 393, 394, 395, 396, 397, 398, 400, 401, 402, 403, 404, 405, 406, 407, 408, 409, 410, 411, 412, 413, 414, 415, 416, 417, 418, 419, 420, 421, 422, 423, 424, 425, 426, 427, 428, 429, 430, 431, 432, 433, 434, 435, 436, 437, 438, 440, 441, 442, 443, 444, 445, 446, 447, 448, 449, 450, 451, 452, 453, 454, 455, 456, 457, 458, 459, 460, 461, 462, 463, 464, 465, 466, 467, 468, 469, 470, 471, 472, 473, 474, 475, 476, 477, 478, 480, 481, 482, 483, 484, 485, 486, 487, 488, 489, 490, 491, 492, 493, 494, 495, 496, 497, 498, 499, 500, 501, 502, 503, 504, 505, 506, 507, 508, 509, 510, 511, 512, 513, 514, 515, 516, 517, 518, 520, 521, 522, 523, 524, 525, 526, 527, 528, 529, 530, 531, 532, 533, 534, 535, 536, 537, 538, 539, 540, 541, 542, 543, 544, 545, 546, 547, 548, 549, 550, 551, 552, 553, 554, 555, 556, 557, 558, 560, 561, 562, 563, 564, 565, 566, 567, 568, 569, 570, 571, 572, 573, 574, 575, 576, 577, 578, 579, 580, 581, 582, 583, 584, 585, 586, 587, 588, 589, 590, 591, 592, 593, 594, 595, 596, 597, 598, 600, 601, 602, 603, 604, 605, 606, 607, 608, 609, 610, 611, 612, 613, 614, 615, 616, 617, 618, 619, 620, 621, 622, 623, 624, 625, 626, 627, 628, 629, 630, 631, 632, 633, 634, 635, 636, 637, 638, 640, 641, 642, 643, 644, 645, 646, 647, 648, 649, 650, 651, 652, 653, 654, 655, 656, 657, 658, 659, 660, 661, 662, 663, 664, 665, 666, 667, 668, 669, 670, 671, 672, 673, 674, 675, 676, 677, 678, 680, 681, 682, 683, 684, 685, 686, 687, 688, 689, 690, 691, 692, 693, 694, 695, 696, 697, 698, 699, 700, 701, 702, 703, 704, 705, 706, 707, 708, 709, 710, 711, 712, 713, 714, 715, 716, 717, 718, 720, 721, 722, 723, 724, 725, 726, 727, 728, 729, 730, 731, 732, 733, 734, 735, 736, 737, 738, 739, 740, 741, 742, 743, 744, 745, 746, 747, 748, 749, 750, 751, 752, 753, 754, 755, 756, 757, 758, 760, 761, 762, 763, 764, 765, 766, 767, 768, 769, 770, 771, 772, 773, 774, 775, 776, 777, 778, 779, 780, 781, 782, 783, 784, 785, 786, 787, 788, 789, 790, 791, 792, 793, 794, 795, 796, 797, 798, 800, 801, 802, 803, 804, 805, 806, 807, 808, 809, 810, 811, 812, 813, 814, 815, 816, 817, 818, 819, 820, 821, 822, 823, 824, 825, 826, 827, 828, 829, 830, 831, 832, 833, 834, 835, 836, 837, 838, 840, 841, 842, 843, 844, 845, 846, 847, 848, 849, 850, 851, 852, 853, 854, 855, 856, 857, 858, 859, 860, 861, 862, 863, 864, 865, 866, 867, 868, 869, 870, 871, 872, 873, 874, 875, 876, 877, 878, 880, 881, 882, 883, 884, 885, 886, 887, 888, 889, 890, 891, 892, 893, 894, 895, 896, 897, 898, 899, 900, 901, 902, 903, 904, 905, 906, 907, 908, 909, 910, 911, 912, 913, 914, 915, 916, 917, 918, 920, 921, 922, 923, 924, 925, 926, 927, 928, 929, 930, 931, 932, 933, 934, 935, 936, 937, 938, 939, 940, 941, 942, 943, 944, 945, 946, 947, 948, 949, 950, 951, 952, 953, 954, 955, 956, 957, 958, 960, 961, 962, 963, 964, 965, 966, 967, 968, 969, 970, 971, 972, 973, 974, 975, 976, 977, 978, 979, 980, 981, 982, 983, 984, 985, 986, 987, 988, 989, 990, 991, 992, 993, 994, 995, 996, 997, 998])

9.3. Post-processing and data overview#

Post-processing involves aggregation of all the fit run results into a single data structure. This can then be analysed statistically and examined for for best-fit results. In the statistical sense, this is essentailly a search for candidate radial matrix elements, based on the assumption that some of the minima found in the \(\chi^2\) hyperspace will be the true results. Even if a clear global minima does not exist, searching for candidate radial matrix elements sets based on clustering of results and multiple local minima is still expected to lead to viable candidates provided that the information content of the dataset is sufficient. However, as discussed elsewhere (see Sect. 4.2), in some cases this may not be the case, and other limitations may apply (e.g. certain parameters may be undefined), or additional data required for unique determination of the radial matrix elements.

For more details on the analysis routines, see the PEMtk documentation [20], particularly the fit fidelity and analysis page, and molecular frame analysis data processing page (full analysis for Ref. [3], illustrating the \(N_2\) case).

# General stats & post-processing to data tables
data.analyseFits()

*** Warning: found MultiIndex for DataFrame data.index - checkDims may have issues with Pandas MultiIndex, but will try anyway.
*** Warning: found MultiIndex for DataFrame data.index - checkDims may have issues with Pandas MultiIndex, but will try anyway.
{ 'Fits': 975,
  'Minima': {'chisqr': 0.2510827613766783, 'redchi': 0.00011080439601795158},
  'Stats': { 'chisqr': min       0.251
mean      0.253
median    0.251
max       2.304
std       0.066
var       0.004
Name: chisqr, dtype: float64,
             'redchi': min       1.108e-04
mean      1.117e-04
median    1.108e-04
max       1.017e-03
std       2.902e-05
var       8.420e-10
Name: redchi, dtype: float64},
  'Success': 971}

# The BLMsetPlot routine will output aggregate fit results.
# Here the spread can be taken as a general indication of the uncertainty of 
# the fitting, and indicate whether the fit is well-characterised/the information 
# content of the data is sufficient.
data.BLMsetPlot(xDim = 't', thres=1e-6)  # With xDim and thres set, for more control over outputs

# Glue plot for later
glue("OCS-fitResultsBLM",data.data['plots']['BLMsetPlot'])
Hide code cell output

Fig. 9.1 Fit overview plot - \(\beta_{L,M}(t)\). Here dashed lines with ‘+’ markers indicates the input data, and bands indicate the mean fit results, where the width is the standard deviation in the fit model results. (See the PEMtk documentation [20] for details, particularly the analysis routines page.)#

# Write aggregate datasets to HDF5 format
# This is more robust than Pickled data, but PEMtk currently only support output for aggregate (post-processed) fit data.

data.processedToHDF5(dataPath = dataPath, outStem = dataFileIn.name, timeStamp=False)

Dumped self.data[fits][dfLong] to /home/jovyan/jake-home/buildTmp/_latest_build/html/doc-source/part2/OCSfitting/OCS_999_fit_3D-test_withNoise_orb13_200723_05-47-50.pickle_dfLong.pdHDF with Pandas .to_hdf() routine.
Dumped data to /home/jovyan/jake-home/buildTmp/_latest_build/html/doc-source/part2/OCSfitting/OCS_999_fit_3D-test_withNoise_orb13_200723_05-47-50.pickle_dfLong.pdHDF with pdHDF.
Dumped self.data[fits][AFxr] to /home/jovyan/jake-home/buildTmp/_latest_build/html/doc-source/part2/OCSfitting/OCS_999_fit_3D-test_withNoise_orb13_200723_05-47-50.pickle_AFxr.pdHDF with Pandas .to_hdf() routine.
Dumped data to /home/jovyan/jake-home/buildTmp/_latest_build/html/doc-source/part2/OCSfitting/OCS_999_fit_3D-test_withNoise_orb13_200723_05-47-50.pickle_AFxr.pdHDF with pdHDF.

# Histogram fit results (reduced chi^2 vs. fit index)
# This may be quite slow for large datasets, setting limited ranges may help

# Use default auto binning
# data.fitHist()

# Example with range set
data.fitHist(thres=1.11e-4, bins=100)

# Glue plot for later
glue("OCS-fitHist",data.data['plots']['fitHistPlot'])

Mask selected 974 results (from 975).

Fig. 9.2 Fit overview plot - \(\chi^2\) vs. fit index. Here bands indicate groupings (local minima) are consistently found.#

Here, Fig. 9.1 shows an overview of the results compared with the input data, and Fig. 9.2 an overview of \(\chi^2\) vs. fit index. Bands in the \(\chi^2\) dimension can indicate groupings (local minima) are consistently found. Assuming each grouping is a viable fit candidate parameter set, these can then be explored in further detail.

9.4. Data exploration#

The general aim in this procedure is to ascertain whether there was a good spread of parameters explored, and a single (or few sets) of best-fit results. There are a few procedures and helper methods for this…

9.4.1. View results#

Single results sets can be viewed in the main data structure, indexed by integers.

# Check keys
fitNumber = 2
data.data[fitNumber].keys()

dict_keys(['AFBLM', 'residual', 'results'])

Here results is an lmFit object, which includes final fit results and information, and AFBLM contains the model (fit) output (i.e. resultant AF-\(\beta_{LM}\) values).

An example is shown below. Of particular note here is which parameters have vary=True - these are included in the fitting - and if there is a column expression, which indicates any parameters defined to have specific relationships (see Chpt. 6). Any correlations found during fitting are also shown, which can also indicate parameters which are related (even if this is not predefined or known a priori).

# Show some results
data.data[fitNumber]['results']
Hide code cell output

Fit Result

9.5. Classify candidate sets#

To probe the minima found, the classifyFits method can be used. This bins results into “candidate” groups, which can then be examined in detail.

# Run with defaults
# data.classifyFits()

# For more control, pass bins
# Here the minima is set at one end, and a %age range used for bins
minVal = data.fitsSummary['Stats']['redchi']['min']    
binRangePC = 1e-5
data.classifyFits(bins = [minVal, minVal + binRangePC*minVal , 20])
success chisqr redchi
count unique top freq count unique top freq count unique top freq
redchiGroup
A 47 1 True 47 47.0 47.0 0.251 1.0 47.0 47.0 1.108e-04 1.0
B 82 1 True 82 82.0 82.0 0.251 1.0 82.0 82.0 1.108e-04 1.0
C 65 1 True 65 65.0 65.0 0.251 1.0 65.0 65.0 1.108e-04 1.0
D 76 1 True 76 76.0 76.0 0.251 1.0 76.0 76.0 1.108e-04 1.0
E 51 1 True 51 51.0 51.0 0.251 1.0 51.0 51.0 1.108e-04 1.0
F 54 1 True 54 54.0 54.0 0.251 1.0 54.0 54.0 1.108e-04 1.0
G 28 1 True 28 28.0 28.0 0.251 1.0 28.0 28.0 1.108e-04 1.0
H 24 1 True 24 24.0 24.0 0.251 1.0 24.0 24.0 1.108e-04 1.0
I 23 1 True 23 23.0 23.0 0.251 1.0 23.0 23.0 1.108e-04 1.0
J 10 1 True 10 10.0 10.0 0.251 1.0 10.0 10.0 1.108e-04 1.0
K 12 1 True 12 12.0 12.0 0.251 1.0 12.0 12.0 1.108e-04 1.0
L 7 1 True 7 7.0 7.0 0.251 1.0 7.0 7.0 1.108e-04 1.0
M 3 1 True 3 3.0 3.0 0.251 1.0 3.0 3.0 1.108e-04 1.0
N 7 1 True 7 7.0 7.0 0.251 1.0 7.0 7.0 1.108e-04 1.0
O 7 1 True 7 7.0 7.0 0.251 1.0 7.0 7.0 1.108e-04 1.0
P 7 1 True 7 7.0 7.0 0.251 1.0 7.0 7.0 1.108e-04 1.0
Q 5 1 True 5 5.0 5.0 0.251 1.0 5.0 5.0 1.108e-04 1.0
R 7 1 True 7 7.0 7.0 0.251 1.0 7.0 7.0 1.108e-04 1.0
S 3 1 True 3 3.0 3.0 0.251 1.0 3.0 3.0 1.108e-04 1.0 
*** Warning: found MultiIndex for DataFrame data.index - checkDims may have issues with Pandas MultiIndex, but will try anyway.
*** Warning: found MultiIndex for DataFrame data.index - checkDims may have issues with Pandas MultiIndex, but will try anyway.
../_images/433498643557ad70cc8ad28008fd103fc9633d838c9bc114f9ff5cb775810caa.png

9.6. Explore candidate result sets#

Drill-down on a candidate set of results, and examine values and spreads. For more details see PEMtk documentation [20], especially the analysis routines page. (See also Sect. 2.3 for details on the plotting libaries implemented here.)

9.6.1. Raw results#

Plot spreads in magnitude and phase parameters. Statistical plots are available for Seaborn and Holoviews backends, with some slightly different options.

# From the candidates, select a group for analysis
selGroup = 'A'

# paramPlot can be used to check the spread on each parameter.
# Plots use Seaborn or Holoviews/Bokeh
# Colour-mapping is controlled by the 'hue' paramter, additionally pass hRound for sig. fig control.
# The remap setting allows for short-hand labels as set in data.lmmu

paramType = 'm' # Set for (m)agnitude or (p)hase parameters
hRound = 14 # Set for cmapping, default may be too small (leads to all grey cmap on points)

data.paramPlot(selectors={'Type':paramType, 'redchiGroup':selGroup}, hue = 'redchi', 
               backend=paramPlotBackend, hvType='violin', 
               returnFlag = True, hRound=hRound, remap = 'lmMap');

*** Warning: found MultiIndex for DataFrame data.index - checkDims may have issues with Pandas MultiIndex, but will try anyway.

paramType = 'p' # Set for (m)agnitude or (p)hase parameters
data.paramPlot(selectors={'Type':paramType, 'redchiGroup':selGroup}, hue = 'redchi', backend=paramPlotBackend, hvType='violin', 
               returnFlag = True, hRound=hRound, remap = 'lmMap'); 

*** Warning: found MultiIndex for DataFrame data.index - checkDims may have issues with Pandas MultiIndex, but will try anyway.

9.6.2. Phases, phase shifts & corrections#

Depending on how the fit was configured, phases may be defined in different ways. To set the phases relative to a speific parameter, and wrap to a specified range, use the phaseCorrection() method. This defaults to using the first parameter as a reference phase, and wraps to \(-\pi:\pi\). The phase-corrected values are output to a new Type, ‘pc’, and a set of normalised magnitudes to ‘n’. Additional settings can be passed for more control, as shown below.

# Run phase correction routine
# Set absFlag=True for unsigned phases (mapped to 0:pi)
# Set useRef=False to set ref phase as 0, otherwise the reference value is set.
phaseCorrParams={'absFlag':True, 'useRef':False}
data.phaseCorrection(**phaseCorrParams)
Hide code cell output 
*** Warning: found MultiIndex for DataFrame data.index - checkDims may have issues with Pandas MultiIndex, but will try anyway.
*** Warning: found MultiIndex for DataFrame data.index - checkDims may have issues with Pandas MultiIndex, but will try anyway.
Set ref param = P_S_P_1_1_n1_1
Param P_S_P_1_1_n1_1 P_S_P_1_n1_1_1 P_S_P_2_1_n1_1 P_S_P_2_n1_1_1 P_S_P_3_1_n1_1 P_S_P_3_n1_1_1 P_S_P_4_1_n1_1 P_S_P_4_n1_1_1 P_S_P_5_1_n1_1 P_S_P_5_n1_1_1 ... P_S_P_7_1_n1_1 P_S_P_7_n1_1_1 S_S_S_0_0_0_1 S_S_S_1_0_0_1 S_S_S_2_0_0_1 S_S_S_3_0_0_1 S_S_S_4_0_0_1 S_S_S_5_0_0_1 S_S_S_6_0_0_1 S_S_S_7_0_0_1
Fit Type redchiGroup
2 m F 5.563e-01 5.563e-01 0.331 0.331 1.015 1.015 0.527 0.527 0.280 0.280 ... 0.773 0.773 0.932 0.398 0.795 0.898 0.735 0.380 0.292 1.029
n F 1.817e-01 1.817e-01 0.108 0.108 0.332 0.332 0.172 0.172 0.091 0.091 ... 0.252 0.252 0.304 0.130 0.260 0.293 0.240 0.124 0.095 0.336
p F -2.200e+00 -2.200e+00 0.945 0.945 -1.391 -1.391 0.167 0.167 0.844 0.844 ... 1.953 1.953 -0.009 0.892 -3.139 2.388 0.197 -0.060 -2.815 1.957
pc F 0.000e+00 0.000e+00 3.138 3.138 0.809 0.809 2.367 2.367 3.044 3.044 ... 2.130 2.130 2.191 3.092 0.938 1.695 2.397 2.140 0.614 2.126
3 m F 1.960e-01 1.960e-01 0.693 0.693 0.979 0.979 0.297 0.297 0.134 0.134 ... 0.824 0.824 0.620 0.183 1.017 0.969 0.521 0.407 0.619 1.040
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
996 pc C 0.000e+00 0.000e+00 3.028 3.028 0.941 0.941 2.804 2.804 3.056 3.056 ... 2.015 2.015 2.681 3.065 0.350 1.702 2.867 1.995 0.101 2.048
998 m J 1.425e-04 1.425e-04 0.759 0.759 0.893 0.893 0.424 0.424 0.161 0.161 ... 0.759 0.759 0.552 0.287 0.820 1.019 0.701 0.530 0.665 0.990
n J 4.655e-05 4.655e-05 0.248 0.248 0.292 0.292 0.139 0.139 0.053 0.053 ... 0.248 0.248 0.180 0.094 0.268 0.333 0.229 0.173 0.217 0.323
p J -2.200e+00 -2.200e+00 2.597 2.597 0.098 0.098 -3.142 -3.142 1.046 1.046 ... -3.142 -3.142 3.142 1.285 -0.396 2.900 2.875 0.582 -0.614 2.933
pc J 0.000e+00 0.000e+00 1.486 1.486 2.298 2.298 0.941 0.941 3.037 3.037 ... 0.941 0.941 0.941 2.798 1.804 1.183 1.208 2.783 1.586 1.150

2072 rows × 22 columns

Examine new data types…

paramType = 'n'
data.paramPlot(selectors={'Type':paramType, 'redchiGroup':selGroup}, hue = 'redchi', 
               backend=paramPlotBackend, hvType='violin', kind='box',
               returnFlag = True, hRound=hRound, remap = 'lmMap');

*** Warning: found MultiIndex for DataFrame data.index - checkDims may have issues with Pandas MultiIndex, but will try anyway.

WARNING:param.Scatter20038: Setting non-parameter attribute kind=box using a mechanism intended only for parameters

paramType = 'pc'
data.paramPlot(selectors={'Type':paramType, 'redchiGroup':selGroup}, hue = 'redchi', 
               backend=paramPlotBackend, hvType='violin', kind='box',
               returnFlag = True, hRound=hRound, remap = 'lmMap');

*** Warning: found MultiIndex for DataFrame data.index - checkDims may have issues with Pandas MultiIndex, but will try anyway.

WARNING:param.Scatter20832: Setting non-parameter attribute kind=box using a mechanism intended only for parameters

9.7. Parameter estimation & fidelity#

For case studies, the fit results can be directly compared to the known input parameters. This should give a feel for how well the data defines the matrix elements (parameters) in this case. In general, probing the correlations and spread of results, and comparing to other (unfitted) results is required to estimate fidelity, see Quantum Metrology Vols. 1 & 2 [4, 9] for further discussion.

9.7.1. Best values and statistics#

To get a final parameter set and associated statistics, based on a subset of the fit results, the paramsReport() method is available. If reference data is available, as for the case studies herein, the paramsCompare() method can also be used to compare with the reference case.

# Parameter summary
data.paramsReport(inds = {'redchiGroup':selGroup})
Hide code cell output 
*** Warning: found MultiIndex for DataFrame data.index - checkDims may have issues with Pandas MultiIndex, but will try anyway.
Set parameter stats to self.paramsSummary.
Type m n p pc
Param Agg
P_S_P_1_1_n1_1 min 2.024e-01 6.611e-02 -2.200 0.000
mean 3.948e-01 1.289e-01 -2.200 0.000
median 4.099e-01 1.339e-01 -2.200 0.000
max 5.163e-01 1.686e-01 -2.200 0.000
std 7.877e-02 2.572e-02 0.000 0.000
var 6.205e-03 6.618e-04 0.000 0.000
P_S_P_1_n1_1_1 min 2.024e-01 6.611e-02 -2.200 0.000
mean 3.948e-01 1.289e-01 -2.200 0.000
median 4.099e-01 1.339e-01 -2.200 0.000
max 5.163e-01 1.686e-01 -2.200 0.000
std 7.877e-02 2.572e-02 0.000 0.000
var 6.205e-03 6.618e-04 0.000 0.000
P_S_P_2_1_n1_1 min 1.000e-04 3.266e-05 -3.142 0.043
mean 3.840e-01 1.254e-01 0.091 1.796
median 3.925e-01 1.282e-01 -0.039 2.012
max 6.380e-01 2.084e-01 3.121 3.106
std 1.374e-01 4.487e-02 1.434 0.830
var 1.888e-02 2.013e-03 2.057 0.688
P_S_P_2_n1_1_1 min 1.000e-04 3.266e-05 -3.142 0.043
mean 3.840e-01 1.254e-01 0.091 1.796
median 3.925e-01 1.282e-01 -0.039 2.012
max 6.380e-01 2.084e-01 3.121 3.106
std 1.374e-01 4.487e-02 1.434 0.830
var 1.888e-02 2.013e-03 2.057 0.688
P_S_P_3_1_n1_1 min 9.460e-01 3.089e-01 -3.142 0.414
mean 1.041e+00 3.401e-01 -0.054 0.888
median 1.038e+00 3.390e-01 -1.311 0.941
max 1.119e+00 3.655e-01 3.142 1.227
std 4.092e-02 1.336e-02 2.709 0.153
var 1.674e-03 1.785e-04 7.341 0.024
P_S_P_3_n1_1_1 min 9.460e-01 3.089e-01 -3.142 0.414
mean 1.041e+00 3.401e-01 -0.054 0.888
median 1.038e+00 3.390e-01 -1.311 0.941
max 1.119e+00 3.655e-01 3.142 1.227
std 4.092e-02 1.336e-02 2.709 0.153
var 1.674e-03 1.785e-04 7.341 0.024
P_S_P_4_1_n1_1 min 2.979e-01 9.728e-02 -3.139 0.642
mean 3.781e-01 1.235e-01 -0.139 1.816
median 3.792e-01 1.238e-01 -0.005 2.135
max 4.374e-01 1.428e-01 3.142 3.074
std 3.503e-02 1.144e-02 1.542 0.654
var 1.227e-03 1.309e-04 2.378 0.428
P_S_P_4_n1_1_1 min 2.979e-01 9.728e-02 -3.139 0.642
mean 3.781e-01 1.235e-01 -0.139 1.816
median 3.792e-01 1.238e-01 -0.005 2.135
max 4.374e-01 1.428e-01 3.142 3.074
std 3.503e-02 1.144e-02 1.542 0.654
var 1.227e-03 1.309e-04 2.378 0.428
P_S_P_5_1_n1_1 min 1.376e-01 4.493e-02 0.597 2.582
mean 2.385e-01 7.790e-02 0.979 2.979
median 2.439e-01 7.966e-02 1.037 2.990
max 2.977e-01 9.721e-02 1.501 3.116
std 4.180e-02 1.365e-02 0.187 0.096
var 1.748e-03 1.864e-04 0.035 0.009
P_S_P_5_n1_1_1 min 1.376e-01 4.493e-02 0.597 2.582
mean 2.385e-01 7.790e-02 0.979 2.979
median 2.439e-01 7.966e-02 1.037 2.990
max 2.977e-01 9.721e-02 1.501 3.116
std 4.180e-02 1.365e-02 0.187 0.096
var 1.748e-03 1.864e-04 0.035 0.009
P_S_P_6_1_n1_1 min 2.518e-01 8.221e-02 -3.142 0.106
mean 4.390e-01 1.434e-01 0.515 1.374
median 4.575e-01 1.494e-01 0.596 0.979
max 5.633e-01 1.840e-01 3.142 2.796
std 7.273e-02 2.375e-02 2.383 0.674
var 5.290e-03 5.641e-04 5.679 0.454
P_S_P_6_n1_1_1 min 2.518e-01 8.221e-02 -3.142 0.106
mean 4.390e-01 1.434e-01 0.515 1.374
median 4.575e-01 1.494e-01 0.596 0.979
max 5.633e-01 1.840e-01 3.142 2.796
std 7.273e-02 2.375e-02 2.383 0.674
var 5.290e-03 5.641e-04 5.679 0.454
P_S_P_7_1_n1_1 min 8.412e-01 2.747e-01 -0.348 1.853
mean 8.528e-01 2.785e-01 0.490 2.125
median 8.540e-01 2.789e-01 -0.090 2.089
max 8.700e-01 2.841e-01 2.129 2.709
std 7.349e-03 2.400e-03 0.935 0.167
var 5.400e-05 5.759e-06 0.875 0.028
P_S_P_7_n1_1_1 min 8.412e-01 2.747e-01 -0.348 1.853
mean 8.528e-01 2.785e-01 0.490 2.125
median 8.540e-01 2.789e-01 -0.090 2.089
max 8.700e-01 2.841e-01 2.129 2.709
std 7.349e-03 2.400e-03 0.935 0.167
var 5.400e-05 5.759e-06 0.875 0.028
S_S_S_0_0_0_1 min 5.400e-01 1.763e-01 -3.142 0.761
mean 7.518e-01 2.455e-01 0.161 1.826
median 7.294e-01 2.382e-01 0.049 2.017
max 9.974e-01 3.257e-01 3.142 2.959
std 1.068e-01 3.488e-02 1.590 0.679
var 1.141e-02 1.217e-03 2.529 0.461
S_S_S_1_0_0_1 min 1.397e-01 4.561e-02 0.546 2.746
mean 3.052e-01 9.966e-02 0.966 2.965
median 3.085e-01 1.007e-01 1.025 2.976
max 3.930e-01 1.284e-01 1.330 3.138
std 5.944e-02 1.941e-02 0.203 0.101
var 3.533e-03 3.767e-04 0.041 0.010
S_S_S_2_0_0_1 min 7.921e-01 2.587e-01 -3.142 0.068
mean 8.599e-01 2.808e-01 0.177 1.348
median 8.508e-01 2.778e-01 0.106 1.160
max 9.820e-01 3.207e-01 3.142 2.398
std 4.777e-02 1.560e-02 2.353 0.635
var 2.282e-03 2.434e-04 5.536 0.403
S_S_S_3_0_0_1 min 8.014e-01 2.617e-01 -0.677 1.524
mean 9.544e-01 3.117e-01 0.380 1.875
median 9.648e-01 3.151e-01 -0.323 1.841
max 1.069e+00 3.491e-01 2.501 2.642
std 5.764e-02 1.882e-02 1.166 0.205
var 3.323e-03 3.543e-04 1.360 0.042
S_S_S_4_0_0_1 min 4.547e-01 1.485e-01 -2.994 0.520
mean 6.150e-01 2.008e-01 0.095 1.810
median 6.241e-01 2.038e-01 0.013 2.095
max 6.911e-01 2.257e-01 3.142 3.137
std 5.920e-02 1.933e-02 1.500 0.642
var 3.504e-03 3.737e-04 2.251 0.412
S_S_S_5_0_0_1 min 2.448e-01 7.994e-02 -3.005 0.784
mean 4.018e-01 1.312e-01 1.237 1.739
median 4.133e-01 1.350e-01 2.223 1.786
max 4.684e-01 1.530e-01 2.540 1.981
std 4.371e-02 1.427e-02 1.553 0.238
var 1.910e-03 2.037e-04 2.412 0.057
S_S_S_6_0_0_1 min 2.452e-01 8.007e-02 -3.142 0.119
mean 4.548e-01 1.485e-01 0.763 1.402
median 4.714e-01 1.540e-01 1.398 1.034
max 5.982e-01 1.953e-01 3.142 2.793
std 8.255e-02 2.696e-02 2.296 0.680
var 6.814e-03 7.267e-04 5.273 0.462
S_S_S_7_0_0_1 min 1.085e+00 3.544e-01 -0.358 1.843
mean 1.098e+00 3.586e-01 0.494 2.156
median 1.098e+00 3.587e-01 -0.054 2.117
max 1.110e+00 3.625e-01 2.140 2.736
std 6.303e-03 2.058e-03 0.905 0.172
var 3.973e-05 4.237e-06 0.819 0.030 
# Parameter comparison
# Note this uses phaseCorrParams as set previously for consistency
data.paramsCompare(phaseCorrParams=phaseCorrParams)
Hide code cell output 
*** Warning: found MultiIndex for DataFrame data.index - checkDims may have issues with Pandas MultiIndex, but will try anyway.
*** Warning: found MultiIndex for DataFrame data.index - checkDims may have issues with Pandas MultiIndex, but will try anyway.
Set ref param = P_S_P_1_1_n1_1
Set parameter comparison to self.paramsSummaryComp.
Param P_S_P_1_1_n1_1 P_S_P_1_n1_1_1 P_S_P_2_1_n1_1 P_S_P_2_n1_1_1 P_S_P_3_1_n1_1 P_S_P_3_n1_1_1 P_S_P_4_1_n1_1 P_S_P_4_n1_1_1 P_S_P_5_1_n1_1 P_S_P_5_n1_1_1 ... P_S_P_7_1_n1_1 P_S_P_7_n1_1_1 S_S_S_0_0_0_1 S_S_S_1_0_0_1 S_S_S_2_0_0_1 S_S_S_3_0_0_1 S_S_S_4_0_0_1 S_S_S_5_0_0_1 S_S_S_6_0_0_1 S_S_S_7_0_0_1
Fit Type
ref m 0.299 0.299 0.735 0.735 0.966 0.966 0.875 0.875 0.541 0.541 ... 0.036 0.036 0.607 0.855 0.967 1.204 0.563 0.535 0.147 0.026
n 0.097 0.097 0.240 0.240 0.315 0.315 0.286 0.286 0.177 0.177 ... 0.012 0.012 0.198 0.279 0.316 0.393 0.184 0.175 0.048 0.008
p -2.200 -2.200 -1.410 -1.410 -2.100 -2.100 -0.800 -0.800 0.131 0.131 ... 2.855 2.855 -2.003 0.983 2.679 -0.535 0.343 2.065 2.594 -1.519
pc 0.000 0.000 0.790 0.790 0.100 0.100 1.401 1.401 2.332 2.332 ... 1.228 1.228 0.197 3.100 1.404 1.665 2.543 2.018 1.489 0.682

4 rows × 22 columns

Param P_S_P_1_1_n1_1 P_S_P_1_n1_1_1 P_S_P_2_1_n1_1 P_S_P_2_n1_1_1 P_S_P_3_1_n1_1 P_S_P_3_n1_1_1 P_S_P_4_1_n1_1 P_S_P_4_n1_1_1 P_S_P_5_1_n1_1 P_S_P_5_n1_1_1 ... P_S_P_7_1_n1_1 P_S_P_7_n1_1_1 S_S_S_0_0_0_1 S_S_S_1_0_0_1 S_S_S_2_0_0_1 S_S_S_3_0_0_1 S_S_S_4_0_0_1 S_S_S_5_0_0_1 S_S_S_6_0_0_1 S_S_S_7_0_0_1
Type Source dType
m mean num 0.395 0.395 0.384 0.384 1.041 1.041 0.378 0.378 0.239 0.239 ... 0.853 0.853 0.752 0.305 0.860 0.954 0.615 0.402 0.455 1.098
ref num 0.299 0.299 0.735 0.735 0.966 0.966 0.875 0.875 0.541 0.541 ... 0.036 0.036 0.607 0.855 0.967 1.204 0.563 0.535 0.147 0.026
diff % 24.383 24.383 91.513 91.513 7.275 7.275 131.508 131.508 126.743 126.743 ... 95.786 95.786 19.306 180.049 12.502 26.162 8.467 33.241 67.712 97.652
num 0.096 0.096 -0.351 -0.351 0.076 0.076 -0.497 -0.497 -0.302 -0.302 ... 0.817 0.817 0.145 -0.549 -0.108 -0.250 0.052 -0.134 0.308 1.072
std % 19.952 19.952 35.782 35.782 3.929 3.929 9.266 9.266 17.524 17.524 ... 0.862 0.862 14.208 19.475 5.556 6.039 9.625 10.877 18.151 0.574
num 0.079 0.079 0.137 0.137 0.041 0.041 0.035 0.035 0.042 0.042 ... 0.007 0.007 0.107 0.059 0.048 0.058 0.059 0.044 0.083 0.006
diff/std % 122.212 122.212 255.754 255.754 185.143 185.143 1419.210 1419.210 723.234 723.234 ... 11115.642 11115.642 135.881 924.509 225.018 433.192 87.969 305.597 373.043 17012.358
n mean num 0.129 0.129 0.125 0.125 0.340 0.340 0.123 0.123 0.078 0.078 ... 0.278 0.278 0.246 0.100 0.281 0.312 0.201 0.131 0.149 0.359
ref num 0.097 0.097 0.240 0.240 0.315 0.315 0.286 0.286 0.177 0.177 ... 0.012 0.012 0.198 0.279 0.316 0.393 0.184 0.175 0.048 0.008
diff % 24.386 24.386 91.506 91.506 7.278 7.278 131.499 131.499 126.734 126.734 ... 95.786 95.786 19.309 180.038 12.498 26.158 8.470 33.236 67.713 97.652
num 0.031 0.031 -0.115 -0.115 0.025 0.025 -0.162 -0.162 -0.099 -0.099 ... 0.267 0.267 0.047 -0.179 -0.035 -0.082 0.017 -0.044 0.101 0.350
std % 19.952 19.952 35.782 35.782 3.929 3.929 9.266 9.266 17.524 17.524 ... 0.862 0.862 14.208 19.475 5.556 6.039 9.625 10.877 18.151 0.574
num 0.026 0.026 0.045 0.045 0.013 0.013 0.011 0.011 0.014 0.014 ... 0.002 0.002 0.035 0.019 0.016 0.019 0.019 0.014 0.027 0.002
diff/std % 122.226 122.226 255.734 255.734 185.231 185.231 1419.117 1419.117 723.186 723.186 ... 11115.592 11115.592 135.902 924.455 224.942 433.115 88.005 305.551 373.049 17012.501
p mean num -2.200 -2.200 0.091 0.091 -0.054 -0.054 -0.139 -0.139 0.979 0.979 ... 0.490 0.490 0.161 0.966 0.177 0.380 0.095 1.237 0.763 0.494
ref num -2.200 -2.200 -1.410 -1.410 -2.100 -2.100 -0.800 -0.800 0.131 0.131 ... 2.855 2.855 -2.003 0.983 2.679 -0.535 0.343 2.065 2.594 -1.519
diff % 0.000 0.000 1646.083 1646.083 3757.652 3757.652 475.859 475.859 86.591 86.591 ... 483.135 483.135 1341.919 1.698 1417.730 240.975 261.343 67.010 239.854 407.175
num 0.000 0.000 1.501 1.501 2.046 2.046 0.661 0.661 0.848 0.848 ... -2.366 -2.366 2.164 -0.016 -2.503 0.915 -0.248 -0.829 -1.830 2.013
std % 0.000 0.000 1572.473 1572.473 4976.382 4976.382 1110.494 1110.494 19.054 19.054 ... 191.036 191.036 986.064 21.018 1332.830 307.217 1581.735 125.579 300.904 183.068
num 0.000 0.000 1.434 1.434 2.709 2.709 1.542 1.542 0.187 0.187 ... 0.935 0.935 1.590 0.203 2.353 1.166 1.500 1.553 2.296 0.905
diff/std % NaN NaN 104.681 104.681 75.510 75.510 42.851 42.851 454.443 454.443 ... 252.903 252.903 136.088 8.079 106.370 78.438 16.523 53.361 79.711 222.417
pc mean num 0.000 0.000 1.796 1.796 0.888 0.888 1.816 1.816 2.979 2.979 ... 2.125 2.125 1.826 2.965 1.348 1.875 1.810 1.739 1.402 2.156
ref num 0.000 0.000 0.790 0.790 0.100 0.100 1.401 1.401 2.332 2.332 ... 1.228 1.228 0.197 3.100 1.404 1.665 2.543 2.018 1.489 0.682
diff % NaN NaN 56.010 56.010 88.748 88.748 22.861 22.861 21.725 21.725 ... 42.228 42.228 89.190 4.549 4.115 11.173 40.498 16.025 6.220 68.373
num 0.000 0.000 1.006 1.006 0.788 0.788 0.415 0.415 0.647 0.647 ... 0.897 0.897 1.629 -0.135 -0.055 0.209 -0.733 -0.279 -0.087 1.474
std % NaN NaN 46.193 46.193 17.279 17.279 36.044 36.044 3.214 3.214 ... 7.845 7.845 37.182 3.392 47.113 10.955 35.471 13.679 48.485 7.998
num 0.000 0.000 0.830 0.830 0.153 0.153 0.654 0.654 0.096 0.096 ... 0.167 0.167 0.679 0.101 0.635 0.205 0.642 0.238 0.680 0.172
diff/std % NaN NaN 121.253 121.253 513.612 513.612 63.425 63.425 676.042 676.042 ... 538.275 538.275 239.873 134.102 8.735 101.990 114.173 117.152 12.828 854.857

28 rows × 22 columns

# Display above results With column name remapping to (l,m) labels only

# With Pandas functionality
data.paramsSummaryComp.rename(columns=data.lmmu['lmMap'])

# With utility method
# summaryRenamed = pemtk.fit._util.renameParams(data.paramsSummaryComp, data.lmmu['lmMap']) 
# summaryRenamed
Param 1,1 1,-1 2,1 2,-1 3,1 3,-1 4,1 4,-1 5,1 5,-1 ... 7,1 7,-1 0,0 1,0 2,0 3,0 4,0 5,0 6,0 7,0
Type Source dType
m mean num 0.395 0.395 0.384 0.384 1.041 1.041 0.378 0.378 0.239 0.239 ... 0.853 0.853 0.752 0.305 0.860 0.954 0.615 0.402 0.455 1.098
ref num 0.299 0.299 0.735 0.735 0.966 0.966 0.875 0.875 0.541 0.541 ... 0.036 0.036 0.607 0.855 0.967 1.204 0.563 0.535 0.147 0.026
diff % 24.383 24.383 91.513 91.513 7.275 7.275 131.508 131.508 126.743 126.743 ... 95.786 95.786 19.306 180.049 12.502 26.162 8.467 33.241 67.712 97.652
num 0.096 0.096 -0.351 -0.351 0.076 0.076 -0.497 -0.497 -0.302 -0.302 ... 0.817 0.817 0.145 -0.549 -0.108 -0.250 0.052 -0.134 0.308 1.072
std % 19.952 19.952 35.782 35.782 3.929 3.929 9.266 9.266 17.524 17.524 ... 0.862 0.862 14.208 19.475 5.556 6.039 9.625 10.877 18.151 0.574
num 0.079 0.079 0.137 0.137 0.041 0.041 0.035 0.035 0.042 0.042 ... 0.007 0.007 0.107 0.059 0.048 0.058 0.059 0.044 0.083 0.006
diff/std % 122.212 122.212 255.754 255.754 185.143 185.143 1419.210 1419.210 723.234 723.234 ... 11115.642 11115.642 135.881 924.509 225.018 433.192 87.969 305.597 373.043 17012.358
n mean num 0.129 0.129 0.125 0.125 0.340 0.340 0.123 0.123 0.078 0.078 ... 0.278 0.278 0.246 0.100 0.281 0.312 0.201 0.131 0.149 0.359
ref num 0.097 0.097 0.240 0.240 0.315 0.315 0.286 0.286 0.177 0.177 ... 0.012 0.012 0.198 0.279 0.316 0.393 0.184 0.175 0.048 0.008
diff % 24.386 24.386 91.506 91.506 7.278 7.278 131.499 131.499 126.734 126.734 ... 95.786 95.786 19.309 180.038 12.498 26.158 8.470 33.236 67.713 97.652
num 0.031 0.031 -0.115 -0.115 0.025 0.025 -0.162 -0.162 -0.099 -0.099 ... 0.267 0.267 0.047 -0.179 -0.035 -0.082 0.017 -0.044 0.101 0.350
std % 19.952 19.952 35.782 35.782 3.929 3.929 9.266 9.266 17.524 17.524 ... 0.862 0.862 14.208 19.475 5.556 6.039 9.625 10.877 18.151 0.574
num 0.026 0.026 0.045 0.045 0.013 0.013 0.011 0.011 0.014 0.014 ... 0.002 0.002 0.035 0.019 0.016 0.019 0.019 0.014 0.027 0.002
diff/std % 122.226 122.226 255.734 255.734 185.231 185.231 1419.117 1419.117 723.186 723.186 ... 11115.592 11115.592 135.902 924.455 224.942 433.115 88.005 305.551 373.049 17012.501
p mean num -2.200 -2.200 0.091 0.091 -0.054 -0.054 -0.139 -0.139 0.979 0.979 ... 0.490 0.490 0.161 0.966 0.177 0.380 0.095 1.237 0.763 0.494
ref num -2.200 -2.200 -1.410 -1.410 -2.100 -2.100 -0.800 -0.800 0.131 0.131 ... 2.855 2.855 -2.003 0.983 2.679 -0.535 0.343 2.065 2.594 -1.519
diff % 0.000 0.000 1646.083 1646.083 3757.652 3757.652 475.859 475.859 86.591 86.591 ... 483.135 483.135 1341.919 1.698 1417.730 240.975 261.343 67.010 239.854 407.175
num 0.000 0.000 1.501 1.501 2.046 2.046 0.661 0.661 0.848 0.848 ... -2.366 -2.366 2.164 -0.016 -2.503 0.915 -0.248 -0.829 -1.830 2.013
std % 0.000 0.000 1572.473 1572.473 4976.382 4976.382 1110.494 1110.494 19.054 19.054 ... 191.036 191.036 986.064 21.018 1332.830 307.217 1581.735 125.579 300.904 183.068
num 0.000 0.000 1.434 1.434 2.709 2.709 1.542 1.542 0.187 0.187 ... 0.935 0.935 1.590 0.203 2.353 1.166 1.500 1.553 2.296 0.905
diff/std % NaN NaN 104.681 104.681 75.510 75.510 42.851 42.851 454.443 454.443 ... 252.903 252.903 136.088 8.079 106.370 78.438 16.523 53.361 79.711 222.417
pc mean num 0.000 0.000 1.796 1.796 0.888 0.888 1.816 1.816 2.979 2.979 ... 2.125 2.125 1.826 2.965 1.348 1.875 1.810 1.739 1.402 2.156
ref num 0.000 0.000 0.790 0.790 0.100 0.100 1.401 1.401 2.332 2.332 ... 1.228 1.228 0.197 3.100 1.404 1.665 2.543 2.018 1.489 0.682
diff % NaN NaN 56.010 56.010 88.748 88.748 22.861 22.861 21.725 21.725 ... 42.228 42.228 89.190 4.549 4.115 11.173 40.498 16.025 6.220 68.373
num 0.000 0.000 1.006 1.006 0.788 0.788 0.415 0.415 0.647 0.647 ... 0.897 0.897 1.629 -0.135 -0.055 0.209 -0.733 -0.279 -0.087 1.474
std % NaN NaN 46.193 46.193 17.279 17.279 36.044 36.044 3.214 3.214 ... 7.845 7.845 37.182 3.392 47.113 10.955 35.471 13.679 48.485 7.998
num 0.000 0.000 0.830 0.830 0.153 0.153 0.654 0.654 0.096 0.096 ... 0.167 0.167 0.679 0.101 0.635 0.205 0.642 0.238 0.680 0.172
diff/std % NaN NaN 121.253 121.253 513.612 513.612 63.425 63.425 676.042 676.042 ... 538.275 538.275 239.873 134.102 8.735 101.990 114.173 117.152 12.828 854.857

28 rows × 22 columns

Hide code cell content 
# Plot values vs. reference cases
# NOTE - experimental code, not yet consolidated and wrapped in PEMtk

paramType = 'm'

# Set new DataFrame including "vary" info (missing in default case)
pDict = 'dfWideTest'
# Try using existing function with extra index set...
data._setWide(indexDims = ['Fit','Type','chisqrGroup','redchiGroup','batch', 'vary'], dataWide='dfWideTest')

# WITH lmMAP remap - good if (l,m) are unique labels
plotData = data.paramPlot(dataDict = pDict, selectors={'vary':True, 'Type':paramType, 'redchiGroup':selGroup}, hue = 'chisqr', 
                          backend='hv', hvType='violin', returnFlag = True, plotScatter=True, hRound=hRound, remap='lmMap') 

# NO REMAP CASE
# plotData = data.paramPlot(dataDict = pDict, selectors={'vary':True, 'Type':paramType, 'redchiGroup':selGroup}, hue = 'chisqr', 
#                           backend='hv', hvType='violin', returnFlag = True, plotScatter=True, hRound=hRound)  #, remap='lmMap') 

p1 = data.data['plots']['paramPlot']

# Plot ref params... CURRENTLY NOT IN paraPlot(), and that also expects fit data so can't reuse directly here.

dataTest = data.data['fits']['dfRef'].copy()
# data.paramPlot(dataDict = 'dfRef')

# Set axis remap
# dataTest.replace({'Param':data.lmmu['lmMap']}, inplace=True)

# Subset
dataTestSub = data._subsetFromXS(selectors = {'Type':paramType}, data = dataTest)  
p2 = dataTestSub.hvplot.scatter(x='Param',y='value', marker='dash', size=500, color='red')

p1*p2

*** Warning: found MultiIndex for DataFrame data.index - checkDims may have issues with Pandas MultiIndex, but will try anyway.
*** Warning: found MultiIndex for DataFrame data.index - checkDims may have issues with Pandas MultiIndex, but will try anyway.
*** Warning: found MultiIndex for DataFrame data.index - checkDims may have issues with Pandas MultiIndex, but will try anyway.
*** Warning: found MultiIndex for DataFrame data.index - checkDims may have issues with Pandas MultiIndex, but will try anyway.

9.8. Using the reconstructed matrix elements#

The results tables are accessible directly, and there are also methods to reformat the best fit results for use in further calculations.

# self.paramsSummary contains the results above as Pandas Dataframe, usual Pandas methods can be applied.
data.paramsSummary['data'].describe()
Param P_S_P_1_1_n1_1 P_S_P_1_n1_1_1 P_S_P_2_1_n1_1 P_S_P_2_n1_1_1 P_S_P_3_1_n1_1 P_S_P_3_n1_1_1 P_S_P_4_1_n1_1 P_S_P_4_n1_1_1 P_S_P_5_1_n1_1 P_S_P_5_n1_1_1 ... P_S_P_7_1_n1_1 P_S_P_7_n1_1_1 S_S_S_0_0_0_1 S_S_S_1_0_0_1 S_S_S_2_0_0_1 S_S_S_3_0_0_1 S_S_S_4_0_0_1 S_S_S_5_0_0_1 S_S_S_6_0_0_1 S_S_S_7_0_0_1
count 188.000 188.000 188.000 188.000 188.000 188.000 188.000 188.000 188.000 188.000 ... 188.000 188.000 188.000 188.000 188.000 188.000 188.000 188.000 188.000 188.000
mean -0.419 -0.419 0.599 0.599 0.554 0.554 0.545 0.545 1.069 1.069 ... 0.937 0.937 0.746 1.084 0.666 0.880 0.680 0.877 0.692 1.027
std 1.042 1.042 1.083 1.083 1.416 1.416 1.125 1.125 1.162 1.162 ... 0.859 0.859 1.087 1.141 1.298 0.860 1.059 1.012 1.276 0.845
min -2.200 -2.200 -3.142 -3.142 -3.142 -3.142 -3.139 -3.139 0.045 0.045 ... -0.348 -0.348 -3.142 0.046 -3.142 -0.677 -2.994 -3.005 -3.142 -0.358
25% -0.550 -0.550 0.113 0.113 0.337 0.337 0.121 0.121 0.128 0.128 ... 0.277 0.277 0.235 0.137 0.274 0.310 0.198 0.138 0.158 0.358
50% 0.033 0.033 0.322 0.322 0.928 0.928 0.359 0.359 0.447 0.447 ... 0.849 0.849 0.700 0.469 0.820 0.941 0.571 0.422 0.479 1.095
75% 0.177 0.177 1.063 1.063 1.039 1.039 0.952 0.952 1.771 1.771 ... 1.956 1.956 1.025 1.684 0.996 1.754 1.150 1.831 1.086 1.950
max 0.516 0.516 3.121 3.121 3.142 3.142 3.142 3.142 3.116 3.116 ... 2.709 2.709 3.142 3.138 3.142 2.642 3.142 2.540 3.142 2.736

8 rows × 22 columns

# To set matrix elements from aggregate fit results, use `seetAggMatE` for Pandas
data.setAggMatE(simpleForm = True)
data.data['agg']['matEpd']

Set reformatted aggregate data to self.data[agg][matEpd].
Type m n p pc comp compC labels
Cont l m mu
P 1 -1 1 0.395 0.129 -2.200 0.000 -0.232-0.319j 0.129+0.000j 1,-1
1 -1 0.395 0.129 -2.200 0.000 -0.232-0.319j 0.129+0.000j 1,1
2 -1 1 0.384 0.125 0.091 1.796 0.382+0.035j -0.028+0.122j 2,-1
1 -1 0.384 0.125 0.091 1.796 0.382+0.035j -0.028+0.122j 2,1
3 -1 1 1.041 0.340 -0.054 0.888 1.040-0.057j 0.215+0.264j 3,-1
1 -1 1.041 0.340 -0.054 0.888 1.040-0.057j 0.215+0.264j 3,1
4 -1 1 0.378 0.123 -0.139 1.816 0.374-0.052j -0.030+0.120j 4,-1
1 -1 0.378 0.123 -0.139 1.816 0.374-0.052j -0.030+0.120j 4,1
5 -1 1 0.239 0.078 0.979 2.979 0.133+0.198j -0.077+0.013j 5,-1
1 -1 0.239 0.078 0.979 2.979 0.133+0.198j -0.077+0.013j 5,1
6 -1 1 0.439 0.143 0.515 1.374 0.382+0.216j 0.028+0.141j 6,-1
1 -1 0.439 0.143 0.515 1.374 0.382+0.216j 0.028+0.141j 6,1
7 -1 1 0.853 0.278 0.490 2.125 0.753+0.401j -0.147+0.237j 7,-1
1 -1 0.853 0.278 0.490 2.125 0.753+0.401j -0.147+0.237j 7,1
S 0 0 0 0.752 0.246 0.161 1.826 0.742+0.121j -0.062+0.238j 0,0
1 0 0 0.305 0.100 0.966 2.965 0.173+0.251j -0.098+0.017j 1,0
2 0 0 0.860 0.281 0.177 1.348 0.847+0.151j 0.062+0.274j 2,0
3 0 0 0.954 0.312 0.380 1.875 0.886+0.354j -0.093+0.297j 3,0
4 0 0 0.615 0.201 0.095 1.810 0.612+0.058j -0.048+0.195j 4,0
5 0 0 0.402 0.131 1.237 1.739 0.132+0.380j -0.022+0.129j 5,0
6 0 0 0.455 0.149 0.763 1.402 0.329+0.314j 0.025+0.146j 6,0
7 0 0 1.098 0.359 0.494 2.156 0.967+0.521j -0.198+0.299j 7,0 
# To set matrix elements from aggregate fit results, use `aggToXR` for Xarray
# data.aggToXR(refKey = 'orb5', returnType = 'ds', conformDims=True)   # use full ref dataset
data.aggToXR(refKey = 'subset', returnType = 'ds', conformDims=True)   # Subselected matE

Added dim Targ
Added dim Total
Added dim Targ
Added dim Total
Set XR dataset for self.data['agg']['matE']
Hide code cell content 
data.data['agg']['matE']

<xarray.Dataset>
Dimensions:  (Type: 2, Sym: 4, LM: 24, mu: 3, it: 1)
Coordinates:
  * Type     (Type) object 'comp' 'compC'
  * Sym      (Sym) MultiIndex
  - Cont     (Sym) object 'P' 'P' 'S' 'S'
  - Targ     (Sym) object 'S' 'U' 'S' 'U'
  - Total    (Sym) object 'P' 'U' 'S' 'U'
  * LM       (LM) MultiIndex
  - l        (LM) int64 0 0 0 1 1 1 2 2 2 3 3 3 4 4 4 5 5 5 6 6 6 7 7 7
  - m        (LM) int64 -1 0 1 -1 0 1 -1 0 1 -1 0 1 -1 0 1 -1 0 1 -1 0 1 -1 0 1
  * mu       (mu) int64 -1 0 1
  * it       (it) int64 1
    Eke      float64 10.1
    Ehv      float64 27.4
    SF       complex128 (2.8514231+2.8292131j)
Data variables:
    comp     (Type, mu, Sym, LM) complex128 (nan+nanj) (nan+nanj) ... (nan+nanj)
    compC    (Type, mu, Sym, LM) complex128 (nan+nanj) (nan+nanj) ... (nan+nanj)
    subset   (LM, Sym, mu, it) complex128 (nan+nanj) (nan+nanj) ... (nan+nanj)

9.8.1. Density matrices#

New (experimental) code for density matrix plots and comparison. See Sect. 3.4 for discussion. Code adapted from the PEMtk documentation [20] MF reconstruction page, original analysis for Ref. [3], illustrating the \(N_2\) case. If the reconstruction is good, the differences (fidelity) should be on the order of the experimental noise level/reconstruction uncertainty, around 10% in the case studies herein; in general the values and patterns of the matrices can also indicate aspects of the retrieval that worked well, or areas where values are poorly defined/recovered from the given dataset.

Hide code cell content 
# Define phase function to test unsigned phases only
def unsignedPhase(da):
    """Convert to unsigned phases."""
    # Set mag, phase
    mag = da.pipe(np.abs)
    phase = da.pipe(np.angle)  # Returns np array only!
    
    # Set unsigned
    magUS = mag.pipe(np.abs)
#     phaseUS = phase.pipe(np.abs)  
    phaseUS = np.abs(phase)
    
    # Set complex
    compFixed = magUS * np.exp(1j* phaseUS)
#     return mag,phase
    return compFixed
Hide code cell content 
# Compute density matrices for retrieved and reference cases, and compare
# v2 - as v1, but differences for unsigned phase case & fix labels
# 26/07/22: messy but working. Some labelling tricks to push back into matPlot() routine

# Import routines
from epsproc.calc import density

# Compose density matrix

# Set dimensions/state vector/representation
# These must be in original data, but will be restacked as necessary to define the effective basis space.
denDims = ['LM', 'mu']
selDims = {'Type':'L'}
pTypes=['r','i']
thres = 1e-2    # 0.2 # Threshold out l>3 terms if using full 'orb5' set.
normME = False
normDen = 'max'
usPhase = True # Use unsigned phases?

# Calculate - Ref case
# matE = data.data['subset']['matE']
# Set data from master class
# k = 'orb5'  # N2 orb5 (SG) dataset
k = 'subset'
matE = data.data[k]['matE']
if normME:
    matE = matE/matE.max()

if usPhase:
    matE = unsignedPhase(matE)
    
daOut, *_ = density.densityCalc(matE, denDims = denDims, selDims = selDims, thres = thres)  # OK

if normDen=='max':
    daOut = daOut/daOut.max()
elif normDen=='trace':
    daOut = daOut/(daOut.sum('Sym').pipe(np.trace)**2)  # Need sym sum here to get 2D trace
    
# daPlot = density.matPlot(daOut.sum('Sym'))
daPlot = density.matPlot(daOut.sum('Sym'), pTypes=pTypes)

# Retrieved
matE = data.data['agg']['matE']['compC']
selDims = {'Type':'compC'}  # For stacked DS case need to set selDims again here to avoid null data selection below.
if normME:
    matE = matE/matE.max()
    
if usPhase:
    matE = unsignedPhase(matE)
    
daOut2, *_ = density.densityCalc(matE, denDims = denDims, selDims = selDims, thres = thres)  # OK

if normDen=='max':
    daOut2 = daOut2/daOut2.max()
elif normDen=='trace':
    daOut2 = daOut2/(daOut2.sum('Sym').pipe(np.trace)**2)
    
daPlot2 = density.matPlot(daOut2.sum('Sym'), pTypes=pTypes)   #.sel(Eke=slice(0.5,1.5,1)))


# Compute difference
if usPhase:
    daDiff = unsignedPhase(daOut.sum('Sym')) - unsignedPhase(daOut2.sum('Sym'))

else:
    daDiff = daOut.sum('Sym') - daOut2.sum('Sym')

daDiff.name = 'Difference'
daPlotDiff = density.matPlot(daDiff, pTypes=pTypes)



#******** Plot
daLayout = (daPlot.redim(pType='Component').layout('Component').relabel('(a) Reference density matrix (unsigned phases)') + daPlot2.opts(show_title=False).layout('pType').opts(show_title=True).relabel('(b) Reconstructed') + 
                daPlotDiff.opts(show_title=False).layout('pType').opts(show_title=True).relabel('(c) Difference')).cols(1)  
daLayout.opts(ep.plot.hvPlotters.opts.HeatMap(width=300, frame_width=300, aspect='square', tools=['hover'], colorbar=True, cmap='coolwarm'))  # .opts(show_title=False)  # .opts(title="Custom Title")  #OK



# Notes on titles... see https://holoviews.org/user_guide/Customizing_Plots.html
#
# .relabel('Test') and .opts(title="Custom Title") OK for whole row titles
#
# daPlot2.opts(show_title=False).layout('pType').opts(show_title=True).relabel('Recon')  Turns off titles per plot, then titles layout
#
# .redim() to modify individual plot group label (from dimension name) 


# Glue figure for later - real part only in this case
# Also clean up axis labels from default state labels ('LM' and 'LM_p' in this case).
glue("OCS-densityComp", daLayout)

Set plot kdims to ['l,m,mu', 'l,m,mu_p']; pass kdims = [dim1,dim2] for more control.
Set plot kdims to ['l,m,mu', 'l,m,mu_p']; pass kdims = [dim1,dim2] for more control.
Set plot kdims to ['l,m,mu', 'l,m,mu_p']; pass kdims = [dim1,dim2] for more control.

WARNING:param.HeatMapPlot24151: Due to internal constraints, when aspect and width/height is set, the bokeh backend uses those values as frame_width/frame_height instead. This ensures the aspect is respected, but means that the plot might be slightly larger than anticipated. Set the frame_width/frame_height explicitly to suppress this warning.
WARNING:param.HeatMapPlot24159: Due to internal constraints, when aspect and width/height is set, the bokeh backend uses those values as frame_width/frame_height instead. This ensures the aspect is respected, but means that the plot might be slightly larger than anticipated. Set the frame_width/frame_height explicitly to suppress this warning.
WARNING:param.HeatMapPlot24188: Due to internal constraints, when aspect and width/height is set, the bokeh backend uses those values as frame_width/frame_height instead. This ensures the aspect is respected, but means that the plot might be slightly larger than anticipated. Set the frame_width/frame_height explicitly to suppress this warning.
WARNING:param.HeatMapPlot24195: Due to internal constraints, when aspect and width/height is set, the bokeh backend uses those values as frame_width/frame_height instead. This ensures the aspect is respected, but means that the plot might be slightly larger than anticipated. Set the frame_width/frame_height explicitly to suppress this warning.
WARNING:param.HeatMapPlot24223: Due to internal constraints, when aspect and width/height is set, the bokeh backend uses those values as frame_width/frame_height instead. This ensures the aspect is respected, but means that the plot might be slightly larger than anticipated. Set the frame_width/frame_height explicitly to suppress this warning.
WARNING:param.HeatMapPlot24230: Due to internal constraints, when aspect and width/height is set, the bokeh backend uses those values as frame_width/frame_height instead. This ensures the aspect is respected, but means that the plot might be slightly larger than anticipated. Set the frame_width/frame_height explicitly to suppress this warning.

Fig. 9.3 Density matrix comparison - rows show (a) reference case (with signs of phases removed), (b) reconstructed case, (c) differences. Columns are (left) imaginary component, (right) real component. If the reconstruction is good, the differences (fidelity) should be on the order of the experimental noise level/reconstruction uncertainty, around 10% in the case studies herein.#

9.8.2. Plot MF PADs#

Routines below adapted from the PEMtk documentation [20] MF reconstruction data processing page (original analysis page for Ref. [3], illustrating the \(N_2\) case). The routines include calls to self.mfpadNumeric() for numerical expansion of the MF-PADs, and self.padPlot() for plotting. Results are illustrated for the retrieved and reference cases in Fig. 9.4 and Fig. 9.5 respectively, and the differential results (reference minus fitted results) in Fig. 9.6.

Hide code cell content 
dataIn = data.data['agg']['matE'].copy()

# Restack for MFPAD calculation and plotter
# Single Eke dim case

# Create empty ePSbase class instance, and set data
# Can then use existing  padPlot() routine for all data
from epsproc.classes.base import ePSbase
dataTest = ePSbase(verbose = 1)

aList = [i for i in dataIn.data_vars]  # List of arrays

# Loop version & propagate attrs
dataType = 'matE'
for item in aList:
    if item.startswith('sub'):
        dataTest.data[item] = {dataType : dataIn[item]}
    else:
        selType = item
        dataTest.data[item] = {dataType : dataIn[item].sel({'Type':selType})}
        
    # Push singleton Eke value to dim for plotter
    dataTest.data[item][dataType] = dataTest.data[item][dataType].expand_dims('Eke')
Hide code cell content 
# Compute MFPADs for a range of cases

# Set Euler angs to include diagonal pol case
pRot = [0, 0, np.pi/2, 0]
tRot = [0, np.pi/2, np.pi/2, np.pi/4]
cRot = [0, 0, 0, 0]
labels = ['z','x','y', 'd']
eulerAngs = np.array([labels, pRot, tRot, cRot]).T   # List form to use later, rows per set of angles

# Should also use MFBLM function below instead of numeric version?
# Numeric version is handy for direct surface and difference case.
R = ep.setPolGeoms(eulerAngs = eulerAngs)

# Comparison and diff
pKey = [i for i in dataIn.data_vars if i!='comp']  # List of arrays
dataTest.mfpadNumeric(keys=pKey, R = R)   # Compute MFPADs for each set of matrix elements using numerical routine

dataTest.data['diff'] = {'TX': dataTest.data['subset']['TX'].sum('Sym')-dataTest.data['compC']['TX'].sum('Sym')}  # Add sum over sym to force matching dims
pKey.extend(['diff'])

# Plot - all cases
# Now run in separate cells below for more stable output
# Erange=[1,2,1]  # Set for a range of Ekes
# Eplot = {'Eke':data.selOpts['matE']['inds']['Eke']}  # Plot for selected Eke (as used for fitting)
# print(f"\n*** Plotting for keys = {pKey}, one per row ***\n")  # Note plot labels could do with some work!
# dataTest.padPlot(keys=pKey, Erange=Erange, backend='pl',returnFlag=True, plotFlag=True) # Generate plotly polar surf plots for each dataset

# Change default plotting config, and then plot in separate cells below
ep.plot.hvPlotters.setPlotters(width=1000, height=500)

* Set Holoviews with bokeh.
Hide code cell content 
# Plot results from reconstructed matE
pKey = 'compC'
print(f"\n*** Plotting for keys = {pKey} ***\n")  # Note plot labels could do with some work!
# dataTest.padPlot(keys=pKey, Erange=Erange, backend='pl',returnFlag=True, plotFlag=True) # Generate plotly polar surf plots for each dataset
Eplot = {'Eke':data.selOpts['matE']['inds']['Eke']}  # Plot for selected Eke (as used for fitting)
dataTest.padPlot(keys=pKey, selDims=Eplot, backend='pl',returnFlag=True, plotFlag=True) # Generate plotly polar surf plots for each dataset

# And GLUE for display later with caption
figObj = dataTest.data[pKey]['plots']['TX']['polar'][0]
glue("OCS-compC", figObj)

*** Plotting for keys = compC ***

Summing over dims: {'Sym'}
Plotting from self.data[compC][TX], facetDims=['Labels', None], pType=a2 with backend=pl.
*** WARNING: plot dataset has min value < 0, min = (-0.3476252085951104+1.2320432095903766j). This may be unphysical and/or result in plotting issues.
Set plot to self.data['compC']['plots']['TX']['polar']

Fig. 9.4 MF-PADs computed from retrieved matrix elements for \((x,y,z,d)\) polarization geometries, where \(d\) is the “diagonal” case with the polarization axis as 45 degrees to the \(z\)-axis.#

Hide code cell content 
# Plot results from reference matE
pKey = 'subset'
print(f"\n*** Plotting for keys = {pKey} ***\n")  # Note plot labels could do with some work!
# dataTest.padPlot(keys=pKey, Erange=Erange, backend='pl',returnFlag=True, plotFlag=True) # Generate plotly polar surf plots for each dataset
Eplot = {'Eke':data.selOpts['matE']['inds']['Eke']}  # Plot for selected Eke (as used for fitting)
dataTest.padPlot(keys=pKey, selDims=Eplot, backend='pl',returnFlag=True, plotFlag=True) # Generate plotly polar surf plots for each dataset

# And GLUE for display later with caption
figObj = dataTest.data[pKey]['plots']['TX']['polar'][0]
glue("OCS-ref", figObj)

*** Plotting for keys = subset ***

Summing over dims: {'Sym'}
Plotting from self.data[subset][TX], facetDims=['Labels', None], pType=a2 with backend=pl.
*** WARNING: plot dataset has min value < 0, min = (-1.0650374116229477+0.052779510753376044j). This may be unphysical and/or result in plotting issues.
Set plot to self.data['subset']['plots']['TX']['polar']

Fig. 9.5 MF-PADs computed from reference ab initio matrix elements for \((x,y,z,d)\) polarization geometries, where \(d\) is the “diagonal” case with the polarization axis as 45 degrees to the \(z\)-axis.#

Hide code cell content 
# Plot normalised differences
pKey = 'diff'
print(f"\n*** Plotting for keys = {pKey} ***\n")  # Note plot labels could do with some work!
# dataTest.padPlot(keys=pKey, Erange=Erange, backend='pl',returnFlag=True, plotFlag=True) # Generate plotly polar surf plots for each dataset
Eplot = {'Eke':data.selOpts['matE']['inds']['Eke']}  # Plot for selected Eke (as used for fitting)
dataTest.padPlot(keys=pKey, selDims=Eplot, backend='pl',returnFlag=True, plotFlag=True) # Generate plotly polar surf plots for each dataset

# And GLUE for display later with caption
figObj = dataTest.data[pKey]['plots']['TX']['polar'][0]
glue("OCS-diff", figObj)

*** Plotting for keys = diff ***

Summing over dims: set()
Plotting from self.data[diff][TX], facetDims=['Labels', None], pType=a2 with backend=pl.
*** WARNING: plot dataset has min value < 0, min = (-1.4260618607121236+0.17721150531570146j). This may be unphysical and/or result in plotting issues.
Set plot to self.data['diff']['plots']['TX']['polar']

Fig. 9.6 MF-PADs differences between retrieved and reference cases for \((x,y,z,d)\) polarization geometries, where \(d\) is the “diagonal” case with the polarization axis as 45 degrees to the \(z\)-axis. Note diffs are normalised to emphasize the shape, but not mangnitudes, of the differences - see the density matrix comparisons for a more rigourous fidelity analysis.#

Hide code cell content 
# Check max differences (abs values)
maxDiff = dataTest.data['diff']['plots']['TX']['pData'].max(dim=['Theta','Phi'])   #.sum(['Theta','Phi']).max()   #.max(dim='Eke')
maxDiff.to_pandas()

Labels
z    2.065
x    3.011
y    3.008
d    2.122
dtype: float64
Hide code cell content 
# Check case without phase correction too - this should indicate poor agreement in general
pKey = 'comp'
dataTest.mfpadNumeric(keys=pKey, R = R) 
dataTest.padPlot(keys=pKey, selDims=Eplot, backend='pl',returnFlag=True, plotFlag=True) # Generate plotly polar surf plots for each dataset

Summing over dims: {'Sym'}
Plotting from self.data[comp][TX], facetDims=['Labels', None], pType=a2 with backend=pl.
*** WARNING: plot dataset has min value < 0, min = (-1.8912066577797655-0.5566477298646336j). This may be unphysical and/or result in plotting issues.
Set plot to self.data['comp']['plots']['TX']['polar']