Stretched exponential

[henrikjacobsenfys]

Implementing the Fourier transform of a stretched exponential. It's not trivial to make it fast and correct at the same time. Claude did that for me. I have updated the text and read through the whole thing.

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 98.52%. Comparing base (f941728) to head (d1827c9).
⚠️ Report is 1 commits behind head on master.

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #245      +/-   ##
==========================================
+ Coverage   98.47%   98.52%   +0.04%     
==========================================
  Files          57       58       +1     
  Lines        5512     5676     +164     
  Branches      952      963      +11     
==========================================
+ Hits         5428     5592     +164     
  Misses         42       42              
  Partials       42       42              

Flag	Coverage Δ
unittests	98.52% <100.00%> (+0.04%)	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

Files with missing lines	Coverage Δ
src/easydynamics/__init__.py	100.00% <100.00%> (ø)
src/easydynamics/analysis/analysis.py	97.68% <100.00%> (ø)
src/easydynamics/analysis/analysis1d.py	99.25% <100.00%> (ø)
src/easydynamics/analysis/parameter_analysis.py	99.27% <100.00%> (ø)
src/easydynamics/analysis/posterior_sampling.py	95.42% <ø> (ø)
src/easydynamics/experiment/experiment.py	98.90% <100.00%> (ø)
src/easydynamics/sample_model/__init__.py	100.00% <100.00%> (ø)
...c/easydynamics/sample_model/components/__init__.py	100.00% <100.00%> (ø)
...s/sample_model/components/stretched_exponential.py	100.00% <100.00%> (ø)
src/easydynamics/utils/plotting.py	87.87% <100.00%> (ø)
... and 1 more

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[rozyczko]

Looks good to merge

[rozyczko]

evalaute -> evaluate

[rozyczko]

Instad -> Instead

[rozyczko]

width returns a DescriptorNumber, whereas every other component's width is a Parameter.
Anything that does comp.width.fixed, comp.width.make_dependent_on(...) or puts it in a FitBinding will fail.

[henrikjacobsenfys]

That was because width in this case is derived from beta and tau. I thought it was just hbar/tau, in which case width should be a Dependent Parameter. But that expression is only true for beta=1 as Fable pointed out. For other values of beta there is no closed expression and width has to be found numerically. Claude found a nice way to do it by interpolation, so width is now a Dependent Parameter (with a somewhat opaque dependency expression)

[rozyczko]

I ran Fable on this PR and it found an issue I can't comment on, since I don't know science well enough.

=====================================================================

width overstates the true half width by orders of magnitude at small β, so the numerical convolution's grid-resolution warning never fires and coarse grids give wrong results silently (medium)

width is exposed as $\Gamma = \hbar/\tau$
(stretched_exponential.py:493-515) so that the shared
width-versus-grid checks in
NumericalConvolutionBase._check_width_thresholds can see how wide the
component is. But $\Gamma$ is only the HWHM at $\beta = 1$. Below that,
the peak sharpens dramatically — the height grows as
$\Gamma(1 + 1/\beta)$ while the area stays fixed — and the actual HWHM
collapses:

$\beta$	HWHM / $\Gamma$	peak height relative to a Lorentzian of the same $\Gamma$
1.0	1.00	1.0
0.7	0.56	1.3
0.5	0.22	2.0
0.3	1.5e-2	9.3
0.2	2.7e-4	120
0.1	< 1e-8	3.6e6

So a component with $\tau = 5$ ps, $\beta = 0.2$ reports
width = 0.13 meV (comfortably above a 0.1 meV grid step, no warning),
while its real HWHM is $3.5 \times 10^{-5}$ meV. The numerical
convolution then samples a sub-grid spike and returns garbage. Measured
with a 0.05 meV Gaussian resolution on $[-2, 2]$ meV, default
upsample_factor=5:

$\beta$	grid points	convolved peak	area on $[-2,2]$	warnings
0.3	81	3.71	0.818	none
0.3	401	2.71	0.693	none
0.3	40001 (converged)	2.64	0.684	none
0.2	81	28.7	3.85	none
0.2	401	7.40	1.18	none
0.2	4001	3.03	0.630	none
0.2	40001 (converged)	2.75	0.595	none

At $\beta = 0.2$ on an 81-point grid the area is wrong by a factor of
six and nothing tells the user. $\beta = 0.6$ is fine on every grid, so
this only bites in the small-$\beta$ regime — but that is exactly the
regime where people reach for a KWW instead of a Lorentzian.

Note this is not a bug in _kww_shape; the direct evaluate is exact.
It is a mismatch between what width promises to the convolution
machinery and what the profile actually looks like.

Suggested fix (pick one):

1. Make width report an actual half width. There is no closed form,
   but a bisection on _kww_shape for the $w$ where
   $G_\beta(w) = \tfrac12,\Gamma(1 + 1/\beta)$ costs ~40 single-point
   quadratures (microseconds) and width is not on the evaluation hot
   path. Then the existing small-width warning fires exactly when it
   should. Keep the docstring note that this is derived and read-only.
2. Failing that, document loudly in the class docstring that $\Gamma$ is
   not the HWHM for $\beta < 1$, that the numerical convolution needs
   upsample_factor raised for small $\beta$, and add a regression test
   that the convolution converges (compare 401 vs 4001 points at
   $\beta = 0.3$).

Either way, add a test that a narrow small-$\beta$ component on a coarse
grid produces the "small compared to the spacing" warning — today it
does not.

[henrikjacobsenfys]

I ran Fable on this PR and it found an issue I can't comment on, since I don't know science well enough.

=====================================================================

width overstates the true half width by orders of magnitude at small β, so the numerical convolution's grid-resolution warning never fires and coarse grids give wrong results silently (medium)

width is exposed as Γ = ℏ / τ (stretched_exponential.py:493-515) so that the shared width-versus-grid checks in NumericalConvolutionBase._check_width_thresholds can see how wide the component is. But Γ is only the HWHM at β = 1 . Below that, the peak sharpens dramatically — the height grows as Γ ( 1 + 1 / β ) while the area stays fixed — and the actual HWHM collapses:

𝛽
HWHM /
Γ
peak height relative to a Lorentzian of the same
Γ
1.0 1.00 1.0
0.7 0.56 1.3
0.5 0.22 2.0
0.3 1.5e-2 9.3
0.2 2.7e-4 120
0.1 < 1e-8 3.6e6
So a component with τ = 5 ps, β = 0.2 reports width = 0.13 meV (comfortably above a 0.1 meV grid step, no warning), while its real HWHM is 3.5 × 10 − 5 meV. The numerical convolution then samples a sub-grid spike and returns garbage. Measured with a 0.05 meV Gaussian resolution on [ − 2 , 2 ] meV, default upsample_factor=5:

𝛽
grid points convolved peak area on
[
−
2
,
2
]
warnings
0.3 81 3.71 0.818 none
0.3 401 2.71 0.693 none
0.3 40001 (converged) 2.64 0.684 none
0.2 81 28.7 3.85 none
0.2 401 7.40 1.18 none
0.2 4001 3.03 0.630 none
0.2 40001 (converged) 2.75 0.595 none
At β = 0.2 on an 81-point grid the area is wrong by a factor of six and nothing tells the user. β = 0.6 is fine on every grid, so this only bites in the small-$\beta$ regime — but that is exactly the regime where people reach for a KWW instead of a Lorentzian.

Note this is not a bug in _kww_shape; the direct evaluate is exact. It is a mismatch between what width promises to the convolution machinery and what the profile actually looks like.

Suggested fix (pick one):

1. Make width report an actual half width. There is no closed form,
   but a bisection on _kww_shape for the
   w
   where

   G
   β
   

   (
   w
   )

     1
     2
   

   ,
   Γ
   (
   1
   +
   1

   /
   

   β
   )
   costs ~40 single-point
   quadratures (microseconds) and width is not on the evaluation hot
   path. Then the existing small-width warning fires exactly when it
   should. Keep the docstring note that this is derived and read-only.

2. Failing that, document loudly in the class docstring that
   Γ
   is
   not the HWHM for
   β
   <
   1
   , that the numerical convolution needs
   upsample_factor raised for small
   β
   , and add a regression test
   that the convolution converges (compare 401 vs 4001 points at

   β

   0.3
   ).

Either way, add a test that a narrow small-$\beta$ component on a coarse grid produces the "small compared to the spacing" warning — today it does not.

Very good find, it led into a small rabbit hole. I've escaped it with the help of Opus