Get error distribution#

Before running this notebook, you must first follow the instruction in README.md to:

  1. Create the Python environment with all dependant Python packages like PyTorch

  2. Activate that python environment

  3. Install ptyrad package into your activated Python environement (only need to install once)

Note: This notebook is designed for visualization of the reconstruction error.

Author: Chia-Hao Lee, cl2696@cornell.edu

00. Setup working directory and imports#

import os
import numpy as np
import torch
import matplotlib.pyplot as plt

# Set this to your desired working directory so you can easily access the data, model, param files
work_dir = "H:/workspace/ptyrad/"

os.chdir(work_dir)
print("Current working dir: ", os.getcwd())
# Note that the output/ directory will be automatically generated under your working directory
Current working dir:  H:\workspace\ptyrad
from ptyrad.core.models import PtychoModel
from ptyrad.init import Initializer
from ptyrad.io.handlers import save_array
from ptyrad.io.load import load_ptyrad
from ptyrad.params import load_params
from ptyrad.plotting.model import plot_forward_pass
from ptyrad.runtime.device import set_gpu_device
from ptyrad.runtime.diagnostics import print_system_info
from ptyrad.runtime.logging import LoggingManager
from ptyrad.solver.reconstruction import make_batches, select_scan_indices

01. Setup file path and parse params#

model_path   = "output/tBL_WSe2/20250908_full_N16384_dp128_flipT100_random32_p6_1obj_1slice_plr1e-4_oalr5e-4_oplr5e-4_slr2e-3_orblur0.4_ozblur1.0_mamp0.03_4.0_oathr0.96_oposc_sng1.0_spr0.1/model_iter0200.hdf5"
params_path  = "output/tBL_WSe2/20250908_full_N16384_dp128_flipT100_random32_p6_1obj_1slice_plr1e-4_oalr5e-4_oplr5e-4_slr2e-3_orblur0.4_ozblur1.0_mamp0.03_4.0_oathr0.96_oposc_sng1.0_spr0.1/tBL_WSe2_reconstruct.yml"
ckpt         = load_ptyrad(model_path)
params       = load_params(params_path, validate=True)

init_params  = params.get('init_params')
model_params = params.get('model_params')
recon_params = params.get('recon_params')
WARNING: Probe aberrations '['probe_defocus', 'probe_c3', 'probe_c5']' in 'init_params' are depracated since PtyRAD v0.1.0b13 and are automatically converted to 'probe_aberrations' dict.
 
\\altas.cac.cornell.edu\cl2696\workspace\ptyrad\src\ptyrad\params\hypertune_params.py:173: UserWarning: hypertune_params.tune_params contains deprecated field 'defocus'. It is currently inactive and ignored, but will be removed in future versions. Please update your params file to remove it.
  warnings.warn(
\\altas.cac.cornell.edu\cl2696\workspace\ptyrad\src\ptyrad\params\hypertune_params.py:173: UserWarning: hypertune_params.tune_params contains deprecated field 'c3'. It is currently inactive and ignored, but will be removed in future versions. Please update your params file to remove it.
  warnings.warn(
\\altas.cac.cornell.edu\cl2696\workspace\ptyrad\src\ptyrad\params\hypertune_params.py:173: UserWarning: hypertune_params.tune_params contains deprecated field 'c5'. It is currently inactive and ignored, but will be removed in future versions. Please update your params file to remove it.
  warnings.warn(
# Modify the init_params so we're loading the reconstructed model
init_params['obj_source'] = 'PtyRAD'
init_params['probe_source'] = 'PtyRAD'
init_params['pos_source'] = 'PtyRAD'
init_params['obj_params'] = model_path
init_params['probe_params'] = model_path
init_params['pos_params'] = model_path

02. Initialize logger, device, initializer, and model#

LoggingManager(log_file='ptyrad_log.txt', log_dir='auto', prefix_time='datetime', show_timestamp=True)

print_system_info()
device = set_gpu_device(gpuid=0) # Pass in `gpuid = None` if you don't have access to a CUDA-compatible GPU. Note that running PtyRAD with CPU would be much slower than on GPU.
init = Initializer(init_params).init_all()
model = PtychoModel(init.init_variables, model_params, device=device)
2026-02-13 22:54:03,404 - ### PtyRAD Logger configuration ###
2026-02-13 22:54:03,404 - log_file       = 'ptyrad_log.txt'. If log_file = None, no log file will be created.
2026-02-13 22:54:03,404 - log_dir        = 'auto'. If log_dir = 'auto', then log will be saved to `output_path` or 'logs/'.
2026-02-13 22:54:03,404 - flush_file     = True. Automatically set to True if `log_file is not None`
2026-02-13 22:54:03,404 - prefix_time    = datetime. If true, preset strings ('date', 'time', 'datetime'), or a string of time format, a datetime str is prefixed to the `log_file`.
2026-02-13 22:54:03,404 - prefix_jobid   = '0'. If not 0, it'll be prefixed to the log file. This is used for hypertune mode with multiple GPUs.
2026-02-13 22:54:03,411 - append_to_file = True. If true, logs will be appended to the existing file. If false, the log file will be overwritten.
2026-02-13 22:54:03,412 - show_timestamp = True. If true, the printed information will contain a timestamp.
2026-02-13 22:54:03,414 -  
2026-02-13 22:54:03,415 - ### System information ###
2026-02-13 22:54:03,416 - Platform: Windows-10-10.0.19045-SP0
2026-02-13 22:54:03,416 - Operating System: Windows 10
2026-02-13 22:54:03,416 - OS Version: 10.0.19045
2026-02-13 22:54:03,416 - Machine: AMD64
2026-02-13 22:54:03,416 - Processor: Intel64 Family 6 Model 85 Stepping 4, GenuineIntel
2026-02-13 22:54:03,416 - Available CPU cores: 12
2026-02-13 22:54:03,416 - Total Memory: 63.49 GB
2026-02-13 22:54:03,416 - Available Memory: 36.90 GB
2026-02-13 22:54:03,416 -  
2026-02-13 22:54:03,416 - ### GPU information ###
2026-02-13 22:54:03,482 - CUDA Available: True
2026-02-13 22:54:03,483 - CUDA Version: 12.8
2026-02-13 22:54:03,488 - Available CUDA GPUs: ['Quadro P5000', 'Quadro P4000']
2026-02-13 22:54:03,488 - CUDA Compute Capability: ['6.1', '6.1']
2026-02-13 22:54:03,488 -   INFO: For torch.compile with Triton, you'll need CUDA GPU with Compute Capability >= 7.0.
2026-02-13 22:54:03,488 -         In addition, Triton does not directly support Windows.
2026-02-13 22:54:03,488 -         For Windows users, please follow the instruction and download `triton-windows` from https://github.com/woct0rdho/triton-windows.
2026-02-13 22:54:03,568 - MIG (Multi-Instance GPU) mode = False
2026-02-13 22:54:03,569 -   INFO: MIG splits a physical GPU into multiple GPU slices, but multiGPU does not support these MIG slices.
2026-02-13 22:54:03,571 -         In addition, multiGPU is currently only available on Linux due to the limited NCCL support.
2026-02-13 22:54:03,572 -       -> If you're doing normal reconstruction/hypertune, you can safely ignore this.
2026-02-13 22:54:03,573 -       -> If you want to do multiGPU, you must provide multiple 'full' GPUs that are not in MIG mode.
2026-02-13 22:54:03,574 -  
2026-02-13 22:54:03,576 - ### Python information ###
2026-02-13 22:54:03,576 - Python Executable: c:\Users\chiahao3\miniforge3\envs\ptyrad\python.exe
2026-02-13 22:54:03,577 - Python Version: 3.12.11 | packaged by Anaconda, Inc. | (main, Jun  5 2025, 12:58:53) [MSC v.1929 64 bit (AMD64)]
2026-02-13 22:54:03,578 -  
2026-02-13 22:54:03,579 - ### Packages information ###
2026-02-13 22:54:03,586 - Numpy Version (metadata): 2.1.2
2026-02-13 22:54:03,588 - PyTorch Version (metadata): 2.7.1+cu128
2026-02-13 22:54:03,592 - Optuna Version (metadata): 4.4.0
2026-02-13 22:54:03,596 - Accelerate Version (metadata): 1.9.0
2026-02-13 22:54:03,601 - PtyRAD Version (ptyrad/__init__.py): 0.1.0b13
2026-02-13 22:54:03,601 - PtyRAD is located at: \\altas.cac.cornell.edu\cl2696\workspace\ptyrad\src\ptyrad\__init__.py
2026-02-13 22:54:03,602 -  
2026-02-13 22:54:03,603 - ### Setting GPU Device ###
2026-02-13 22:54:03,604 - Selected GPU device: cuda:0 (Quadro P5000)
2026-02-13 22:54:03,605 -  
2026-02-13 22:54:03,606 - init_params are displayed below:
2026-02-13 22:54:03,606 -   random_seed: None
2026-02-13 22:54:03,607 -   probe_illum_type: electron
2026-02-13 22:54:03,607 -   probe_kv: 80.0
2026-02-13 22:54:03,607 -   probe_conv_angle: 24.9
2026-02-13 22:54:03,607 -   probe_aberrations: {}
2026-02-13 22:54:03,607 -   beam_kev: None
2026-02-13 22:54:03,607 -   probe_dRn: None
2026-02-13 22:54:03,607 -   probe_Rn: None
2026-02-13 22:54:03,607 -   probe_D_H: None
2026-02-13 22:54:03,607 -   probe_D_FZP: None
2026-02-13 22:54:03,607 -   probe_Ls: None
2026-02-13 22:54:03,607 -   meas_Npix: 128
2026-02-13 22:54:03,607 -   pos_N_scans: 16384
2026-02-13 22:54:03,607 -   pos_N_scan_slow: 128
2026-02-13 22:54:03,607 -   pos_N_scan_fast: 128
2026-02-13 22:54:03,607 -   pos_scan_step_size: 0.429
2026-02-13 22:54:03,607 -   meas_calibration: {'mode': 'fitRBF', 'value': None}
2026-02-13 22:54:03,623 -   probe_pmode_max: 6
2026-02-13 22:54:03,623 -   probe_pmode_init_pows: [0.02]
2026-02-13 22:54:03,623 -   obj_omode_max: 1
2026-02-13 22:54:03,623 -   obj_omode_init_occu: {'occu_type': 'uniform', 'init_occu': None}
2026-02-13 22:54:03,623 -   obj_Nlayer: 1
2026-02-13 22:54:03,623 -   obj_slice_thickness: 12.0
2026-02-13 22:54:03,628 -   meas_permute: None
2026-02-13 22:54:03,629 -   meas_reshape: None
2026-02-13 22:54:03,631 -   meas_flipT: [1, 0, 0]
2026-02-13 22:54:03,631 -   meas_crop: None
2026-02-13 22:54:03,632 -   meas_pad: {'mode': None, 'padding_type': 'power', 'target_Npix': 256, 'value': 0.0, 'threshold': 70}
2026-02-13 22:54:03,634 -   meas_resample: {'mode': None, 'scale_factors': [2.0, 2.0]}
2026-02-13 22:54:03,635 -   meas_add_source_size: None
2026-02-13 22:54:03,636 -   meas_add_detector_blur: None
2026-02-13 22:54:03,636 -   meas_remove_neg_values: {'mode': 'clip_neg', 'value': None, 'force': False}
2026-02-13 22:54:03,636 -   meas_normalization: {'mode': 'max_at_one', 'value': None}
2026-02-13 22:54:03,636 -   meas_add_poisson_noise: None
2026-02-13 22:54:03,636 -   meas_export: None
2026-02-13 22:54:03,636 -   probe_permute: None
2026-02-13 22:54:03,636 -   probe_z_shift: None
2026-02-13 22:54:03,636 -   probe_normalization: {'mode': 'mean_total_ints', 'value': None}
2026-02-13 22:54:03,636 -   pos_scan_flipT: None
2026-02-13 22:54:03,636 -   pos_scan_affine: None
2026-02-13 22:54:03,636 -   pos_scan_rand_std: 0.15
2026-02-13 22:54:03,636 -   obj_z_crop: None
2026-02-13 22:54:03,636 -   obj_z_pad: None
2026-02-13 22:54:03,636 -   obj_z_resample: None
2026-02-13 22:54:03,636 -   meas_source: file
2026-02-13 22:54:03,636 -   meas_params: {'path': 'data\\tBL_WSe2\\Panel_g-h_Themis\\scan_x128_y128.raw', 'key': None, 'shape': None, 'offset': None, 'gap': None}
2026-02-13 22:54:03,636 -   probe_source: PtyRAD
2026-02-13 22:54:03,636 -   probe_params: output/tBL_WSe2/20250908_full_N16384_dp128_flipT100_random32_p6_1obj_1slice_plr1e-4_oalr5e-4_oplr5e-4_slr2e-3_orblur0.4_ozblur1.0_mamp0.03_4.0_oathr0.96_oposc_sng1.0_spr0.1/model_iter0200.hdf5
2026-02-13 22:54:03,652 -   pos_source: PtyRAD
2026-02-13 22:54:03,652 -   pos_params: output/tBL_WSe2/20250908_full_N16384_dp128_flipT100_random32_p6_1obj_1slice_plr1e-4_oalr5e-4_oplr5e-4_slr2e-3_orblur0.4_ozblur1.0_mamp0.03_4.0_oathr0.96_oposc_sng1.0_spr0.1/model_iter0200.hdf5
2026-02-13 22:54:03,652 -   obj_source: PtyRAD
2026-02-13 22:54:03,652 -   obj_params: output/tBL_WSe2/20250908_full_N16384_dp128_flipT100_random32_p6_1obj_1slice_plr1e-4_oalr5e-4_oplr5e-4_slr2e-3_orblur0.4_ozblur1.0_mamp0.03_4.0_oathr0.96_oposc_sng1.0_spr0.1/model_iter0200.hdf5
2026-02-13 22:54:03,652 -   tilt_source: simu
2026-02-13 22:54:03,652 -   tilt_params: {'tilt_type': 'all', 'init_tilts': [[0.0, 0.0]]}
2026-02-13 22:54:03,652 -  
2026-02-13 22:54:03,652 - ### Initializing cache ###
2026-02-13 22:54:03,652 - Loading 'PtyRAD' file from output/tBL_WSe2/20250908_full_N16384_dp128_flipT100_random32_p6_1obj_1slice_plr1e-4_oalr5e-4_oplr5e-4_slr2e-3_orblur0.4_ozblur1.0_mamp0.03_4.0_oathr0.96_oposc_sng1.0_spr0.1/model_iter0200.hdf5 for caching
2026-02-13 22:54:03,754 - use_cached_obj   = True
2026-02-13 22:54:03,755 - use_cached_probe = True
2026-02-13 22:54:03,756 - use_cached_pos   = True
2026-02-13 22:54:03,758 -  
2026-02-13 22:54:03,759 - ### Initializing measurements ###
2026-02-13 22:54:03,761 - Loading measurements from source = 'file'
2026-02-13 22:54:03,761 - Detected measurement file type = '.raw'
2026-02-13 22:54:03,761 - WARNING: Couldn't find the 'shape' in 'meas_params' with file type = '.raw'.
2026-02-13 22:54:03,761 - It is strongly recommended to provide an explicit shape to better load from .raw files
2026-02-13 22:54:03,761 - PtyRAD will still try to load the dataset based on the provided 'init_params', but you may consider setting 'shape': (N_scans, Npix, Npix) inside your 'meas_params' dict.
2026-02-13 22:54:13,628 - Original measurements dtype is float32, casting to float32 (single precision) for computational efficiency.
2026-02-13 22:54:13,628 - Imported meausrements shape / dtype = (16384, 128, 128), dtype = float32
2026-02-13 22:54:14,027 - Imported meausrements int. statistics (min, mean, max) = (-30.8785, 1814.7064, 99441.8047)
2026-02-13 22:54:14,028 - Flipping measurements with [flipud, fliplr, transpose] = [1, 0, 0]
2026-02-13 22:54:14,896 - Removing negative values in measurement with method = clip_neg and value = None
2026-02-13 22:54:15,032 - Minimum value = -30.8785, negative values are clipped to 0 due to the positive px value constraint of measurements
2026-02-13 22:54:17,498 - Normalizing measurements with mode = 'max_at_one' and value = 'None'
2026-02-13 22:54:18,034 - Normalizing by max of the 2D mean pattern intensity: 67097.391
2026-02-13 22:54:18,150 - meausrements shape / dtype = (16384, 128, 128), dtype = float32
2026-02-13 22:54:18,707 - meausrements int. statistics (min, mean, max) = (0.0000, 0.0271, 1.4821)
2026-02-13 22:54:19,040 - No negative values found in measurements. Skipping non-neg correction.
2026-02-13 22:54:19,770 - Pattern total int. statistics (min, mean, max)        = (438.8489, 443.3112, 447.6965), with min/max = 98.0%
2026-02-13 22:54:20,272 - Global meausrements int. statistics (min, mean, max)  = (0.0000, 0.0271, 1.4821)
2026-02-13 22:54:20,272 - measurements                              (N, Ky, Kx) = float32, (16384, 128, 128)
2026-02-13 22:54:20,272 -  
2026-02-13 22:54:20,272 - ### Setting up calibration ###
2026-02-13 22:54:20,272 - meas_calibration mode = 'fitRBF', value = None
2026-02-13 22:54:20,272 - Using loaded raw averaged measurement (before crop/pad/resample) to fit RBF as a part of the meas calibration
2026-02-13 22:54:20,272 - Radius of fitted bright field disk (RBF) = 11.64 px with Npix = 128
2026-02-13 22:54:20,284 - Suggested probe_mask_k radius (RBF*2/Npix) > 0.1819
2026-02-13 22:54:20,284 - Fitting raw averaged measurement with center, radius, and Gaussian blur std as a sanity check
2026-02-13 22:54:20,284 - Note that the fitted Gaussian blur std (detector blur) would be affected by overlapping Bragg disks
2026-02-13 22:54:20,707 - Initial guess: center=(63.78, 63.93), radius=11.64, Gaussian blur std=0.50
2026-02-13 22:54:20,724 - Final fit: center=(63.78, 63.93), radius=11.64, Gaussian blur std=0.54
2026-02-13 22:54:20,725 - WARNING: The 'fitRBF' and 'RBF' calibration methods are highly dependent on the accuracy of user-provided experimental parameters and acquisition conditions,
2026-02-13 22:54:20,726 -          including convergence angle, kV, dose, specimen thickness, and collection angle for the estimation of RBF.
2026-02-13 22:54:20,728 -          For example, a 5-10% error in convergence angle is fairly common.
2026-02-13 22:54:20,728 -          Users are strongly advised to perform proper microscope calibration to ensure accurate results.
2026-02-13 22:54:20,729 -          These method should only be used as a rough estimate and not as a substitute for proper experimental calibration.
2026-02-13 22:54:20,730 - dx (real space pixel size of probe and object) set to 0.1526 Ang with Npix = 128
2026-02-13 22:54:20,731 -  
2026-02-13 22:54:20,731 - ### Setting init_variables dict ###
2026-02-13 22:54:20,733 - Derived values given input init_params:
2026-02-13 22:54:20,733 -   kv          = 80.0 kV
2026-02-13 22:54:20,733 -   wavelength  = 0.0418 Ang
2026-02-13 22:54:20,733 -   conv_angle  = 24.9 mrad
2026-02-13 22:54:20,738 -   Npix        = 128 px
2026-02-13 22:54:20,739 -   dk          = 0.0512 Ang^-1
2026-02-13 22:54:20,740 -   kMax        = 3.2773 Ang^-1
2026-02-13 22:54:20,742 -   da          = 2.1383 mrad
2026-02-13 22:54:20,742 -   angleMax    = 136.8515 mrad
2026-02-13 22:54:20,743 -   RBF         = 11.6447 px (Inferred from the given calibration, NOT necessarily from the loaded measurement data)
2026-02-13 22:54:20,744 -   n_alpha     = 5.4960 (# conv_angle)
2026-02-13 22:54:20,745 -   dx          = 0.1526 Ang, Nyquist-limited dmin = 2*dx = 0.3051 Ang
2026-02-13 22:54:20,748 -   Rayleigh-limited resolution  = 1.0230 Ang (0.61*lambda/alpha for focused probe )
2026-02-13 22:54:20,749 -   Real space probe extent = 19.5282 Ang
2026-02-13 22:54:20,750 -  
2026-02-13 22:54:20,750 - ### Initializing probe ###
2026-02-13 22:54:20,751 - Loading probe from source = 'PtyRAD'
2026-02-13 22:54:20,752 - Loaded probe shape = (6, 128, 128), dtype = complex64
2026-02-13 22:54:20,754 - pmode_now: 6 and pmode_max: 6, leaving the pmode unchanged.
2026-02-13 22:54:20,756 - Orthogonalizing 6 pmodes
2026-02-13 22:54:20,762 - Sorting 6 pmodes by their intensities
2026-02-13 22:54:20,763 - Normalizing probe intensity with mode = 'mean_total_ints' and value = 'None'
2026-02-13 22:54:20,766 - sum(|probe_data|**2) = 443.31, while meas_total_ints (min, mean, max) = (438.8489, 443.3112, 447.6965)
2026-02-13 22:54:20,767 - probe                         (pmode, Ny, Nx) = complex64, (6, 128, 128)
2026-02-13 22:54:20,768 -  
2026-02-13 22:54:20,769 - ### Initializing probe positions ###
2026-02-13 22:54:20,770 - Loading probe positions from source = 'PtyRAD'
2026-02-13 22:54:20,771 - Applying Gaussian distributed random displacement with std = 0.15 px to scan positions
2026-02-13 22:54:20,774 - crop_pos                                (N,2) = int16, (16384, 2)
2026-02-13 22:54:20,775 - crop_pos 1st and last px coords (y,x)         = ([50, 52], [404, 403])
2026-02-13 22:54:20,777 - crop_pos extent (Ang)                         = [54.77040555 55.22809696]
2026-02-13 22:54:20,778 - probe_pos_shifts                        (N,2) = float32, (16384, 2)
2026-02-13 22:54:20,780 -  
2026-02-13 22:54:20,780 - ### Initializing object ###
2026-02-13 22:54:20,782 - Loading object from source = 'PtyRAD'
2026-02-13 22:54:20,804 - omode_now: 1 and omode_max: 1, leaving the omode unchanged.
2026-02-13 22:54:20,807 - object                    (omode, Nz, Ny, Nx) = complex64, (1, 1, 583, 583)
2026-02-13 22:54:20,808 - object extent                 (Z, Y, X) (Ang) = [12.     88.9447 88.9447]
2026-02-13 22:54:20,809 -  
2026-02-13 22:54:20,810 - ### Initializing omode_occu from 'uniform' ###
2026-02-13 22:54:20,811 - omode_occu                            (omode) = float32, (1,)
2026-02-13 22:54:20,812 -  
2026-02-13 22:54:20,813 - ### Initializing H (Fresnel propagator) ###
2026-02-13 22:54:20,814 - Calculating H with probe_shape = [128. 128.] px, dx = 0.1526 Ang, slice_thickness = 12.0000 Ang, lambd = 0.0418 Ang
2026-02-13 22:54:20,819 - H                                    (Ky, Kx) = complex64, (128, 128)
2026-02-13 22:54:20,821 -  
2026-02-13 22:54:20,822 - ### Initializing obj tilts from = 'simu' ###
2026-02-13 22:54:20,823 - Initialized obj_tilts with init_tilts = [[0.0, 0.0]] (theta_y, theta_x) mrad
2026-02-13 22:54:20,824 - obj_tilts                              (N, 2) = float32, (1, 2)
2026-02-13 22:54:20,825 -  
2026-02-13 22:54:20,827 - ### Checking consistency between input params with the initialized variables ###
2026-02-13 22:54:20,828 - Npix, DP measurements, probe, and H shapes are consistent as '128'
2026-02-13 22:54:20,829 - N_scans, len(meas), N_scan_slow*N_scan_fast, len(crop_pos), and len(probe_pos_shifts) are consistent as '16384'
2026-02-13 22:54:20,830 - obj.shape[0] is consistent with len(omode_occu) as '1'
2026-02-13 22:54:20,831 - obj.shape[1] is consistent with Nlayer as '1'
2026-02-13 22:54:20,835 - crop positions (yx_min=[49 47], yx_max=[536 537]) are well contained inside object canvas (Ny,Nx) = (583, 583).
2026-02-13 22:54:20,835 - obj_tilts is consistent with either 1 or N_scans
2026-02-13 22:54:20,836 - Pass the consistency check of initialized variables, initialization is done!
2026-02-13 22:54:20,837 -  
2026-02-13 22:54:20,838 - ### Collecting reconstruction provenance ###
2026-02-13 22:54:21,175 - Reconstruction provenance is collected and initialized.
2026-02-13 22:54:21,175 - ### Initializing PtychoModel model ###
2026-02-13 22:54:22,707 - ### PtychoModel optimizable variables ###
2026-02-13 22:54:22,707 - obja            : torch.Size([1, 1, 583, 583])    , torch.float32   , device:cuda:0, grad:True , lr:5e-04
2026-02-13 22:54:22,707 - objp            : torch.Size([1, 1, 583, 583])    , torch.float32   , device:cuda:0, grad:True , lr:5e-04
2026-02-13 22:54:22,707 - obj_tilts       : torch.Size([1, 2])              , torch.float32   , device:cuda:0, grad:False, lr:0e+00
2026-02-13 22:54:22,707 - slice_thickness : torch.Size([])                  , torch.float32   , device:cuda:0, grad:False, lr:0e+00
2026-02-13 22:54:22,707 - probe           : torch.Size([6, 128, 128, 2])    , torch.float32   , device:cuda:0, grad:True , lr:1e-04
2026-02-13 22:54:22,707 - probe_pos_shifts: torch.Size([16384, 2])          , torch.float32   , device:cuda:0, grad:True , lr:2e-03
2026-02-13 22:54:22,707 -  
2026-02-13 22:54:22,707 - ### Optimizable variables statitsics ###
2026-02-13 22:54:22,722 - Total measurement values  : 268,435,456
2026-02-13 22:54:22,722 - Total optimizing variables: 909,154
2026-02-13 22:54:22,722 - Overdetermined ratio      : 295.26
2026-02-13 22:54:22,722 -  
2026-02-13 22:54:22,722 - ### Model behavior ###
2026-02-13 22:54:22,722 - Obj preblur               : False
2026-02-13 22:54:22,731 - Tilt propagator           : False
2026-02-13 22:54:22,731 - Change slice thickness    : False
2026-02-13 22:54:22,731 - Sub-px probe shift        : True
2026-02-13 22:54:22,731 - Detector blur             : False
2026-02-13 22:54:22,731 - Preload data              : True
2026-02-13 22:54:22,731 - On-the-fly meas padding   : False
2026-02-13 22:54:22,731 - On-the-fly meas resample  : False
2026-02-13 22:54:22,731 -  
2026-02-13 22:54:22,768 - ### Done initializing PtychoModel model ###
2026-02-13 22:54:22,768 -  
# Run a forward pass to check if the reconstructed model (obj, probe, pos) was correctly loaded
indices = np.random.randint(0, init.init_variables["N_scans"], 2)
dp_power = 0.5
plot_forward_pass(model, indices, dp_power)
../_images/2b3484e6b7f2fc14b9662c77c9c854fe70e9c284e844ac2588a4ed3aaa9ad44c.png

03. Run 1 full foward pass for the error distribution#

INDICES_MODE      = recon_params.get('INDICES_MODE')
batch_config      = recon_params.get('BATCH_SIZE', {})
grad_accumulation = batch_config.get("grad_accumulation", 1)
batch_size        = batch_config.get('size') * grad_accumulation
GROUP_MODE        = recon_params.get('GROUP_MODE')

pos = (model.crop_pos + model.opt_probe_pos_shifts).detach().cpu().numpy()

indices = select_scan_indices(
    init.init_variables['N_scan_slow'],
    init.init_variables['N_scan_fast'],
    subscan_slow=INDICES_MODE.get('subscan_slow'),
    subscan_fast=INDICES_MODE.get('subscan_fast'),
    mode=INDICES_MODE.get('mode', 'random'),
)

batches = make_batches(indices, pos, batch_size, mode=GROUP_MODE)
2026-02-13 22:54:25,453 - Selecting indices with the 'full' mode 
2026-02-13 22:54:26,801 - Generated 512 'random' groups of ~32 scan positions in 0.000 sec
# Initialize error_DP to hold the relative error patterns, note that this would take the same amount of GPU VRAM as the 4D-dataset

Npix = model.get_complex_probe_view().shape[-1]
error_DP = torch.zeros((len(indices), Npix, Npix), device=device, dtype=torch.float32)
dp_pow = 0.5

for batch in batches:
    with torch.no_grad():
        model_DP = model(batch)
        measured_DP = model.get_measurements(batch)
        data_mean   = measured_DP.pow(dp_pow).mean(0).clamp(min=1e-8)
        error_DP[batch] = (model_DP.pow(dp_pow) - measured_DP.pow(dp_pow)) / data_mean # error_DP corresponds to the relative error
error_DP_arr = error_DP.detach().cpu().numpy()

# Saving is optional but you can port it into py4DGUI for interactive visualization
save_array(error_DP_arr, file_dir='./', file_name='relative_error', file_format="hdf5", output_shape=None, append_shape=True,)
Saving array with shape = (16384, 128, 128) and dtype = float32
N_scan_slow = init.init_variables['N_scan_slow']
N_scan_fast = init.init_variables['N_scan_fast']

fig, axs = plt.subplots(nrows=1,ncols=2, figsize = (12,6))
axs[0].set_title('Fourier Space Error Distribution')
im0 = axs[0].imshow(error_DP_arr.mean(0), cmap='seismic', vmin=-0.2, vmax=0.2)
fig.colorbar(im0, shrink=0.6)

axs[1].set_title('Real Space Error Distribution')
im1 = axs[1].imshow(error_DP_arr.mean((-1,-2)).reshape(N_scan_slow, N_scan_fast), cmap='seismic', vmin=-0.2, vmax=0.2)
fig.colorbar(im1, shrink=0.6)
plt.show()
../_images/98e3ca50efe2f5a994586c2ba3a31afee73663f9c956050d7e0df85570199e37.png

Note#

  • Since it’s relative error, it naturally contains positive and negative values for each error_DP (N, ky, kx)

  • If we take the average along N, then inevitably some k pixels would have their errors cancel out.

  • Similarly, when we average along (ky, kx), errors for all k pixels would be averaged together and show an overall trend for each probe position.

Fourier Space Error Distribution#

  • Low k region (center) shows low relative error, which is more or less expected as it contains most of the signal and generally SNR ~ 1/sqrt(signal)

  • High k region has higher relative error, indicating model DP predicts stronger intensity than experiment pattern.

Real Space Error Distribution#

  • On atomic sites the relative error tend to be negative values, indicating model DP predicts overall lower intensity than experiment pattern.