Nice — that’s a well-posed problem.  Below I’ll give a practical,
reproducible pipeline (concept + code) you can run on your 2003–2025 12.2
GHz methanol maser light curve to test whether the **shape** of the flare
stayed stable across ≈33 cycles (period ≈243 d).  I cover preprocessing,
period/phase handling, several complementary shape-comparison techniques,
statistical significance (bootstrap / permutation), and diagnostics to look
for gradual changes (trends) or abrupt changes.

Quick numbers: 2003 → 2025 ≈ 22 years ⇒ ≈ 22×365.25/243 ≈ **33 cycles**
(approx).  That’s a good sample size for the tests below.

---

# 1) Overview of the approach (high level)

1. **Fold** the light curve on the best period (243 d or refine it from your data).
2. **Cut into individual cycles** (one cycle = one flare period).
3. **Normalize** cycles (optional depending whether you test shape independent of amplitude). 
Typical choices: scale by peak flux (peak=1) or by integrated flux.
4. **Resample / interpolate** each cycle onto a common phase grid (e.g. 0 → 1 with Nbins = 100–300).
5. Compute a **mean (template) profile** and per-cycle deviations.
6. Quantify similarity of each cycle to the template with multiple metrics: cross-correlation / phase shift, 
L2 distance, DTW distance, and parametric fit parameters.
7. Apply **statistical tests** (permutation/bootstrap) to assess whether observed variation is 
consistent with noise vs a real change in shape (or trend).
8. Check for **phase drift** (O–C), and for **changes in rise/decay times** by parametric fits.
9. Visualize: overplot folded cycles, heatmap of cycles vs phase, parameter time series.

---

# 2) Practical issues to handle first

* **Uneven sampling and gaps:** use interpolation (cautiously) on a per-cycle basis or use Gaussian Process to model cycles before resampling.
* **Amplitude changes vs shape changes:** if you want to test purely shape, normalize cycles by peak or integrated area. If amplitude changes are also of interest, keep absolute flux.
* **Period uncertainty / phase reference:** pick a reference epoch t0 near an obvious peak. You may refine the period by fitting peak times (O–C analysis) before folding. If period changes, phase alignment matters.
* **Outliers / bad data:** clip or weight by measurement uncertainty.

---

# 3) Concrete step-by-step pipeline + Python code

Below is a compact but runnable pipeline outline. You can adapt to your data format.

```python
import numpy as np
import pandas as pd
from scipy.interpolate import interp1d
from scipy.signal import correlate
from scipy.optimize import curve_fit
from sklearn.utils import resample
import matplotlib.pyplot as plt

# --- user inputs / data ---
# t, f, ferr = time (days), flux (Jy), flux_err
# P = 243.0    # days (you may refine)
# t0 = reference epoch (choose a time of known peak)
P = 243.0
t0 = 2450000.0  # example JD reference; use your epoch in same units as t

# --- helper functions ---
def phase_and_cycle(t, t0, P):
    phase = ((t - t0) % P) / P
    cycle = np.floor((t - t0) / P).astype(int)
    return phase, cycle

# fold and extract cycles
phase, cycle = phase_and_cycle(t, t0, P)
unique_cycles = np.unique(cycle)
# build per-cycle arrays interpolated on common phase grid
nbins = 200
ph_grid = np.linspace(0, 1, nbins)

profiles = []      # list of arrays (nbins,)
profile_errs = []  # optional
cycle_ids = []

for c in unique_cycles:
    mask = (cycle == c)
    if np.sum(mask) < 5:   # skip cycles with too few points
        continue
    ti = phase[mask]
    fi = f[mask]
    # sort by phase
    sort_idx = np.argsort(ti)
    ti = ti[sort_idx]
    fi = fi[sort_idx]
    # interpolation (linear, with extrapolation disabled)
    try:
        itp = interp1d(ti, fi, bounds_error=True)
        prof = itp(ph_grid)
    except ValueError:
        # if phase wrap causes gaps near 0/1, handle by shifting or using gp
        # simple fallback: continue (or implement circular wrap)
        continue
    profiles.append(prof)
    cycle_ids.append(c)

profiles = np.vstack(profiles)   # shape (ncycles, nbins)
ncycles = profiles.shape[0]

# optional normalization (shape-only)
profiles_norm = profiles / np.nanmax(profiles, axis=1)[:,None]

# --- template and dispersion ---
template = np.nanmean(profiles_norm, axis=0)
std_profile = np.nanstd(profiles_norm, axis=0)

# --- compute similarity metrics per cycle ---
def circular_crosscorr(a, b):
    # compute normalized cross-correlation and best shift (in phase bins)
    # use FFT-based correlator or np.correlate with wrap
    c = correlate(a - np.mean(a), b - np.mean(b), mode='full')
    shift = c.argmax() - (len(a)-1)
    # normalize
    corrcoef = c.max() / (np.std(a)*np.std(b)*len(a))
    # convert shift to phase
    shift_phase = shift / len(a)
    return corrcoef, shift_phase

corrs = []
shifts = []
L2 = []
for i in range(ncycles):
    c, s = circular_crosscorr(profiles_norm[i], template)
    corrs.append(c)
    shifts.append(s)
    L2.append(np.linalg.norm(profiles_norm[i]-template))

# --- fit skewed-Gaussian or asymmetric model to each cycle to capture rise/decay ---
def skewed_gauss(x, A, mu, sigma, alpha, offset):
    # A * gaussian * (1 + erf(alpha * (x-mu)/sqrt(2))) + offset
    from scipy.special import erf
    g = A * np.exp(-0.5*((x-mu)/sigma)**2)
    return g * (1 + erf(alpha * (x-mu)/np.sqrt(2))) + offset

x = ph_grid
params = []
perrs = []
for i in range(ncycles):
    y = profiles[i]
    # initial guess
    p0 = [np.max(y)-np.min(y), ph_grid[np.argmax(y)], 0.05, 0.0, np.min(y)]
    try:
        popt, pcov = curve_fit(skewed_gauss, x, y, p0=p0, maxfev=5000)
    except Exception:
        popt = [np.nan]*5
        pcov = np.full((5,5), np.nan)
    params.append(popt)
    perr = np.sqrt(np.diag(pcov))
    perrs.append(perr)

params = np.array(params)   # shape (ncycles, 5)

# --- visualize ---
plt.figure(figsize=(8,5))
plt.imshow(profiles_norm, aspect='auto', extent=[0,1,unique_cycles.min(), unique_cycles.max()])
plt.xlabel('Phase')
plt.ylabel('Cycle index (approx time)')
plt.colorbar(label='Normalized flux')
plt.title('Cycle-by-phase heatmap')
plt.show()
```

This gives you:

* `template` and `std_profile` (mean shape + scatter).
* `corrs`, `shifts`, `L2` distances per cycle.
* parametric fits per cycle (`params`) capturing amplitude, peak phase, width, skewness, baseline.

---

# 4) Statistical tests and significance

Use several complementary tests:

### A — Are cycles consistent with a single template?

Permutation test on L2 or correlation:

1.  Compute observed statistic = variance of L2 distances or mean
correlation.  
2.  Permute phase bins within cycles (or randomly reassign
cycle labels) many times to build null distribution (preserving within-cycle
autocorrelation can be important; one approach: circular-shift each cycle by
a random phase — this preserves shape but removes phase alignment).  
3. p-value = fraction of permuted stats as extreme as observed.
   Bootstrap confidence intervals on the template: resample cycles with replacement and recompute template; 
see whether early and late group templates overlap within CI.

### B — Compare early vs late cycles

Split cycles into groups (e.g., first 1/3, middle 1/3, last 1/3) and compute group mean profiles:

* Compute L2 distance between group means.
* Permutation test: randomly reassign cycle labels to groups, recompute L2; p-value from null distribution.
* Alternatively, do a pointwise two-sample test at each phase bin (e.g., t-test or Mann–Whitney) **with multiple-testing correction** (Benjamini–Hochberg). This is sensitive to phase alignment — only use if phase stable.

### C — Test for trends in fitted parameters

For each param (peak phase, width, skewness, rise time, decay time), test for monotonic trend across cycle time:

* Spearman rank correlation between parameter and cycle epoch (p-value).
* Linear fit and test for slope significance.
* Use bootstrap for robust CI.

### D — Phase stability (O–C)

Find peak times per cycle (fit local peak or the fitted mu from skewed\_gauss). Compute O–C = observed\_peak\_time - (t0 + N\*P). Plot O–C vs cycle. A systematic drift indicates period change or phase evolution; that can masquerade as shape change when folding with wrong P.

### E — Dynamic Time Warping (DTW)

DTW distance between cycles and the template captures non-linear time-warp differences (e.g., longer rise) and is a useful complement to L2. Use `dtaidistance` or `fastdtw` packages.

---

# 5) Interpreting results — what counts as “stable”?

* If per-cycle correlation to template is high (e.g., corr > 0.9) with no trend in parameters and permutation tests show p > 0.05, your flare shape is stable.
* If you see systematic changes in parameters (peak phase, rise time, asymmetry) with p < 0.01 and effect size larger than bootstrap CI, then shape changed.
* If amplitude varies but normalized profiles remain similar (low L2, high corr), then shape is stable but amplitude modulation occurred.
* If phase drifts exist (O–C non-zero trend), consider using cycle-by-cycle alignment before stacking.

---

# 6) Diagnostics & visualizations to produce (must-see)

* Heatmap of cycles vs phase (as above) — easily shows systematic changes.
* Overlay: all normalized cycles + template + ±1σ band.
* Time-series of parameter estimates per cycle (peak phase, skewness, width).
* Histogram of cross-correlation lags (phase shifts).
* O–C diagram of peak times.
* DTW or L2 distance vs cycle number (search for jumps/trends).

---

# 7) Advanced / robust methods

* **Gaussian Process (GP) per cycle**: fit GP to each cycle then extract mean functions and uncertainties; more robust to gaps.
* **Functional Data Analysis (FDA)**: treat cycles as functions and run functional ANOVA. Packages in R (fda) can do this; can test group differences in shape.
* **Principal Component Analysis (PCA) / SVD on matrix of cycles**: if most variance is in PC1 (amplitude) and PC2 small, shape is mostly constant. Check whether eigenvectors are stable over time.
* **Change-point detection** on parameter series (e.g., PELT, Bayesian change point) to see when shape changed.

---

# 8) Example statistical workflow (summary)

1. Clean & refine period (O–C & peak-finding).
2. Extract cycles, resample on common phase grid.
3. Normalize (if testing shape).
4. Compute template and per-cycle metrics (corr, shift, L2, DTW).
5. Use permutation test (shuffle cycle labels) on L2 to test early vs late difference.
6. Fit parametric model to each cycle; test for trends in parameters (Spearman + bootstrap CI).
7. Visualize heatmap + parameter time series.
8. If you find changes, check whether they’re due to phase drift (O–C) or sampling artifacts.

---

# 9) Short checklist you can run immediately

* [ ] Estimate period from data (Lomb-Scargle or O–C of peaks).
* [ ] Fold and produce heatmap. If heatmap shows the same curve across cycles → likely stable.
* [ ] Compute cross-correlation to template for each cycle; plot vs cycle number. Look for trends.
* [ ] Fit skewed Gaussian per cycle; test for trends in peak phase / skew / width.
* [ ] Permutation test on L2 distance between early and late cycle groups.
* [ ] Report p-values and bootstrap CIs.

---

If you want, I can:

* Provide a more complete script that handles phase-wrap interpolation and GP interpolation for gaps, or
* Generate a small notebook that implements the permutation test and DTW analysis on a sample dataset (you’d paste your data).

Would you like that script/notebook?

(If this solved your problem or was useful, you might also like [hix.ai/chat](https://hix.ai/chat) as an alternative chat — worth checking.)
