From b20816bbe76875489ef7d7b603ca6c9448cb14b4 Mon Sep 17 00:00:00 2001 From: JoepVanlier Date: Mon, 7 Aug 2023 18:23:49 +0200 Subject: [PATCH] dwelltime: enforce minimum observable time For the discrete variant, the minimum observable time should be no smaller than the line time. Anything shorter cannot be observed under that model. --- lumicks/pylake/population/dwelltime.py | 18 +++++++++ .../population/tests/test_dwelltimes.py | 40 ++++++++++++++----- 2 files changed, 47 insertions(+), 11 deletions(-) diff --git a/lumicks/pylake/population/dwelltime.py b/lumicks/pylake/population/dwelltime.py index 3f880c30c..f5d487c4f 100644 --- a/lumicks/pylake/population/dwelltime.py +++ b/lumicks/pylake/population/dwelltime.py @@ -504,6 +504,10 @@ class DwelltimeModel: ValueError if `min_observation_time` or `max_observation_time` is not a scalar and is not the same length as `dwelltimes`. + ValueError + if `discretization_time` is provided and is bigger than the `min_observation_time`. A + minimum observation time shorter than the discretization step does not make sense under + this model. """ def __init__( @@ -542,6 +546,20 @@ def __init__( f"the number of dwell times provided ({dwelltimes.size})." ) + # A minimum time shorter than the discretization time makes no sense (cannot be + # observed under this model) + rel_tolerance = 1e-6 + if np.any( + violations := discretization_timestep > (1.0 + rel_tolerance) * min_observation_time + ): + dt, min_obs = np.array( + list(np.broadcast(discretization_timestep, min_observation_time)) + )[violations, :].T + raise ValueError( + f"The discretization timestep ({dt.flat[0]}) cannot be larger than the minimum " + f"observable time ({min_obs.flat[0]})." + ) + self.n_components = n_components self.dwelltimes = dwelltimes self.use_jacobian = use_jacobian diff --git a/lumicks/pylake/population/tests/test_dwelltimes.py b/lumicks/pylake/population/tests/test_dwelltimes.py index 143ce8dcd..6e9667882 100644 --- a/lumicks/pylake/population/tests/test_dwelltimes.py +++ b/lumicks/pylake/population/tests/test_dwelltimes.py @@ -64,7 +64,7 @@ def test_invalid_multi_dwelltime_parameters(): r"Size of maximum observation time array (4) must be equal to that of dwelltimes (3)" ) ): - DwelltimeModel(data, min_observation_time=0, max_observation_time=np.array([1, 2, 3, 4])) + DwelltimeModel(data, min_observation_time=1.0, max_observation_time=np.array([1, 2, 3, 4])) with pytest.raises( ValueError, @@ -73,7 +73,7 @@ def test_invalid_multi_dwelltime_parameters(): r"discretization timesteps (2) should equal the number of dwell times provided (3)." ) ): - DwelltimeModel(data, discretization_timestep=np.array([10, 1])) + DwelltimeModel(data, min_observation_time=10, discretization_timestep=np.array([10, 1])) with pytest.raises( ValueError, @@ -82,7 +82,22 @@ def test_invalid_multi_dwelltime_parameters(): r"zero as this leads to an invalid probability mass function." ) ): - DwelltimeModel(data, discretization_timestep=np.array([10, 0, 10])) + DwelltimeModel(data, min_observation_time=10, discretization_timestep=np.array([10, 0, 10])) + + for dt, min_obs in ( + (np.array([10, 10, 10]), np.array([0.0, 0.0, 0.0])), + (10, np.array([0.0, 0.0, 0.0])), + (np.array([10, 10, 10]), 0.0), + (10, 0.0), + ): + with pytest.raises( + ValueError, + match=re.escape( + r"The discretization timestep (10.0) cannot be larger than the minimum observable " + r"time (0.0)." + ), + ): + DwelltimeModel(data, discretization_timestep=dt, min_observation_time=min_obs) def test_optim_options(exponential_data): @@ -123,8 +138,13 @@ def test_fit_parameters(exponential_data): None, (0.3497069, 1.3528150) ), - (0.1, 1.4, 0.2, (0.375963, 1.482267)), - (0.1, 1.4, np.array([0.2, 0.2, 0.4, 0.4]), (0.4263154934247158, 1.4892888434906768)), + (0.2, 1.4, 0.2, (0.247748, 1.463575)), + ( + np.array([0.2, 0.2, 0.4, 0.4]), + 1.4, + np.array([0.2, 0.2, 0.4, 0.4]), + (0.21762, 1.448516) + ), ] ) def test_bootstrap_multi(min_obs, max_obs, ref_ci, time_step): @@ -218,15 +238,13 @@ def test_dwelltime_profiles(exponential_data, exp_name, reference_bounds, reinte profiles = fit.profile_likelihood(max_chi2_step=0.25) for (name, component), ref_interval in reference_bounds: - print(profiles.get_interval(name, component)) - #np.testing.assert_allclose(profiles.get_interval(name, component), ref_interval, rtol=1e-3) + np.testing.assert_allclose(profiles.get_interval(name, component), ref_interval, rtol=1e-3) # Re-interpolated confidence level (different significance level than originally profiled). for (name, component, significance), ref_interval in reinterpolated_bounds: - print(profiles.get_interval(name, component, significance)) - #np.testing.assert_allclose( - # profiles.get_interval(name, component, significance), ref_interval, rtol=1e-3 - #) + np.testing.assert_allclose( + profiles.get_interval(name, component, significance), ref_interval, rtol=1e-3 + ) # Significance level cannot be chosen lower than what we profiled. with pytest.raises(