Generate Synthetic Data for the Multi-Cycle Fitting Example#

Produces the data.csv / energy.csv / time.csv used by ../example.ipynb.
The data is a single 2D dataset: a GLP peak on a linear background (the same system as 01_basic_fitting). After t = 0 the peak position alternates between two repeating subcycles — a negative exponential shift and a mirrored positive one with a slower relaxation — convolved with a Gaussian instrument response function (IRF).

Ground truth model

LinBack (linear background) + GLP peak — see models_energy_truth.yaml
Multi-cycle time dependence on the peak position GLP_01_x0: ['IRF', 'MonoExpNeg', 'MonoExpPos'] at frequency=0.25 — see models_time_truth.yaml

Re-run this notebook to regenerate the data in place.

from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import trspecfit

1. Define the Ground Truth Model#

Build a 2D model on the same energy/time axes the example fits over. The *_truth.yaml files hold the ground truth parameters. frequency=0.25 means the dynamics repeat with period 1/0.25 = 4 time units; the two subcycle models split each period equally, so each subcycle lasts 2 time units.

project = trspecfit.Project(path=Path.cwd(), config_file=None)

# Save the axes first and reload them, so the simulation runs on bit-identical
# axes to what ../example.ipynb later loads. Subcycle masks switch discretely
# at period boundaries, so even 1e-7 axis differences from CSV rounding would
# flip the subcycle assignment of samples that land exactly on a boundary.
data_fmt = '%.6e'
np.savetxt('energy.csv', np.arange(0, 20, 0.05), fmt=data_fmt)    # 400 points
np.savetxt('time.csv', np.arange(-4, 16, 0.05), fmt=data_fmt)     # 400 points
energy = np.loadtxt('energy.csv')
time = np.loadtxt('time.csv')

file = trspecfit.File(
    parent_project=project,
    energy=energy,
    time=time,
)

# Load the ground truth energy model and add the multi-cycle time dependence.
file.load_model(model_yaml='models_energy_truth.yaml', model_info='base')
file.add_time_dependence(
    target_model='base',
    target_parameter='GLP_01_x0',
    dynamics_yaml='models_time_truth.yaml',
    dynamics_model=['IRF', 'MonoExpNeg', 'MonoExpPos'],
    frequency=0.25,
)

file.describe_model()

Loading model to describe energy- (and time-)dependent data: base
Loading model to describe time-dependence of a model parameter: GLP_01_x0 of IRF_MonoExpNeg_MonoExpPos model
model name: base
function: LinBack from trspecfit.functions.energy
function: GLP from trspecfit.functions.energy
all lmfit.Parameters() [flattened and sorted alphabetically]:
Name                        Value      Min      Max   Stderr     Vary     Expr Brute_Step
GLP_01_A                       12     -inf      inf     None     True     None     None
GLP_01_F                        1     -inf      inf     None     True     None     None
GLP_01_m                     0.25     -inf      inf     None     True     None     None
GLP_01_x0                       9     -inf      inf     None     True     None     None
GLP_01_x0_expFun_01_A          -2     -inf      inf     None     True     None     None
GLP_01_x0_expFun_01_t0          0     -inf      inf     None     True     None     None
GLP_01_x0_expFun_01_tau       0.4     -inf      inf     None     True     None     None
GLP_01_x0_expFun_02_A           2     -inf      inf     None    False -GLP_01_x0_expFun_01_A     None
GLP_01_x0_expFun_02_t0          0     -inf      inf     None     True     None     None
GLP_01_x0_expFun_02_tau       0.8     -inf      inf     None     True     None     None
GLP_01_x0_gaussCONV_SD       0.15    1e-06        1     None     True     None     None
LinBack_b                       2     -inf      inf     None     True     None     None
LinBack_m                    0.08     -inf      inf     None     True     None     None
LinBack_xStart                  0     -inf      inf     None     True     None     None
LinBack_xStop                  20     -inf      inf     None     True     None     None

<site-packages>/trspecfit/mcp.py:2883: UserWarning: 7 time sample(s) lie within a relative tolerance of 1e-06 of a subcycle boundary (t = 2, 4, 6, 8, 10, ...). Subcycle assignment switches discretely there, so the model prediction for these samples can silently change with the floating-point representation of the time axis (e.g. after saving it to text and reloading). Shift these samples off the boundaries or sample strictly inside subcycles.
  self.normalize_time()  # update the normalization of the time axis

2. Simulate and Save#

Photon-counting detection draws each pixel from a Poisson distribution, so the noise is signal-dependent (shot noise) — the realistic regime for counting detectors. Lower counts_per_delay ⇒ noisier data.

sim = trspecfit.Simulator(
    model=file.model_active,
    detection='photon_counting',
    counts_per_delay=5000,
    seed=42,
)
clean, noisy, noise = sim.simulate_2d()

# Save noisy data (data is [time, energy], matching File's expectation).
# The axes were already saved (and reloaded) above.
np.savetxt('data.csv', noisy, delimiter=',', fmt=data_fmt)

emask = (energy >= 5) & (energy <= 18)
snr = clean[:, emask].max() / np.std(noise[:, emask])
print(f'data.csv  |  shape {noisy.shape}  |  peak SNR ~ {snr:.1f}')

data.csv  |  shape (400, 400)  |  peak SNR ~ 14.4

3. Quick Visual Check#

The 2D map shows the peak position alternating between the two subcycles after t = 0; the 1D slice makes the noise level visible. In the 2D map, look for the sawtooth: sharp (IRF-blurred) jumps at each subcycle boundary (t = 0, 2, 4, …) relaxing with two different time constants.

fig, axes = plt.subplots(1, 2, figsize=(10, 4))

im = axes[0].pcolormesh(energy, time, noisy, shading='auto')
axes[0].set_xlabel('Energy')
axes[0].set_ylabel('Time')
axes[0].set_title('Simulated 2D data')
plt.colorbar(im, ax=axes[0])

t_ind = np.abs(time - 0.2).argmin()
axes[1].plot(energy, noisy[t_ind], lw=0.8, label='noisy')
axes[1].plot(energy, clean[t_ind], lw=2.0, label='clean (truth)')
axes[1].set_xlabel('Energy')
axes[1].set_ylabel('Intensity')
axes[1].set_title(f'Slice at t = {time[t_ind]:.1f} (early subcycle 1)')
axes[1].legend()

plt.tight_layout()
plt.show()

../../../../_images/1c0bd1e47de4864606952bc5f3e4a057e6cdcd8726028e662d4d6ffe69b9043e.png

Ground Truth Parameters#

Parameter	Value	Description
`LinBack m`	0.08	Background slope
`LinBack b`	2.0	Background intercept (at `xStart`)
`LinBack xStart` / `xStop`	0 / 20	Background clamp range
`GLP A`	12	Peak amplitude
`GLP x0`	9	Peak center (static part)
`GLP F`	1.0	Peak width
`GLP m`	0.25	Gaussian/Lorentzian mixing
`frequency`	0.25	Repetition rate: period 4, two 2-unit subcycles
`gaussCONV SD`	0.15	IRF width (global element 0)
`expFun_01 A`	-2	Subcycle-1 shift amplitude
`expFun_01 tau`	0.4	Subcycle-1 relaxation time
`expFun_02 A`	`-expFun_01_A` = +2	Subcycle-2 amplitude (linked)
`expFun_02 tau`	0.8	Subcycle-2 relaxation time (independent)

Model files#

The YAML files next to this notebook, as shipped with the example.

`models_energy_truth.yaml`#

# ground truth energy model for data generation
# (same structure as the fit model in ../models_energy.yaml)

base:
    LinBack:
      m: [8E-2, True]
      b: [2.0, True]
      xStart: [0, True]
      xStop: [20, True]
    GLP:
      A: [12, True]
      x0: [9, True]
      F: [1.0, True]
      m: [0.25, True]

`models_time_truth.yaml`#

# ground truth time models for data generation
# Multi-cycle dynamics: ['IRF', 'MonoExpNeg', 'MonoExpPos'] at frequency=0.25,
# i.e. a 4-unit period split into two 2-unit subcycles (after t = 0).

# global element 0: Gaussian instrument response, convolves all time traces
IRF:
    gaussCONV:
        SD: [0.15, True, 1.0E-6, 1]

# subcycle 1: negative shift, decays back to 0 on its local clock
MonoExpNeg:
    expFun:
        A: [-2, True]
        tau: [0.4, True]
        t0: [0, True]

# subcycle 2: amplitude mirrors subcycle 1 via expression, tau independent
MonoExpPos:
    expFun:
        A: ["-expFun_01_A"]
        tau: [0.8, True]
        t0: [0, True]

Generate Synthetic Data for the Multi-Cycle Fitting Example

Contents