interpolate_simulation_result does not yield consistent timepoints

[nijawa]

Bug description

The function mio::interpolate_simulation_result does not always yield the expected number of timepoints. If the simulation goes from 0 to days, interpolate_simulation_result does not always result in interpolated timepoints from 0 to days.

Version

Windows

To reproduce

1. For me this was easiest seen in pycode\examples\simulation\sde_sirs_simple.py but any example should work really
2. Set days = 5 and run, printed time points are 0, 1, 2, 3, 4, 5, expected 0, 1, 2, 3, 4, 5
3. Set days = 4 and run, printed timepoints are 0, 1, 2, 3, expected 0, 1, 2, 3, 4

Relevant log output

Example results:

For t = 5:

Time             Susceptible      Infected         Recovered       
         0.00000       8900.00000        100.00000       1000.00000
         1.00000          1.60501       8355.46295       1642.93204
         2.00000          3.15233       7586.63010       2410.21757
         3.00000          5.47813       6887.20686       3107.31502
         4.00000          8.72828       6250.83330       3740.43843
         5.00000         13.44976       5719.91700       4266.63324

For t = 4

Time             Susceptible      Infected         Recovered       
         0.00000       8900.00000        100.00000       1000.00000
         1.00000          1.93701       8383.50728       1614.55571
         2.00000          2.54839       7612.70909       2384.74252
         3.00000          7.00763       6935.67804       3057.31434

Add any relevant information, e.g. used compiler, screenshots.

Checklist

• Attached labels, especially loc:: or model:: labels.
• Linked to project

[mknaranja]

Thanks for this important report. However, I assume to know why this happens. It solves the ODE from 0 to day 4 in double arithmetic with a potentially adaptive integrator which then yields a last evaluation point at 4-1e-14 for instance. I think we already at some point discussed tolerances for this functions and integrators but don't know about the current state. Could you check on that?

[nijawa]

The example in questions is a SDE model and the integrator should use a fixed width of 1e-3. I will take a deeper look into this.

[nijawa]

@mknaranja This indeed seems to be a floating point issue. On my local machine, the last timepoint for days = 4 and dt = 1e-3 seems to be 3.99999999999967048580629125353880226612091064453125. This plus the tolerance of 1e-14 does not exceed 4 and as such the timepoint 4 is not interpolated. For days = 5 I see the last point being hit perfectly 5.00000000000000000000000000000000000000000000000000. I am not familiar enough with memilios FP handling and how to handle cumulative error of those floating point operations to suggest a good fix but I would assume (dynamically?) lowering the tolerance could be enough.

[mknaranja]

Actually, the tolerance is really 1e-14 (was my guess before...)

memilio/cpp/memilio/data/analyze_result.h

Line 45 in 590381f

TimeSeries<FP> interpolate_simulation_result(const TimeSeries<FP>& simulation_result, const FP abs_tol = 1e-14);

and the critical lines are here:

memilio/cpp/memilio/data/analyze_result.h

Lines 188 to 190 in 590381f

	const auto day0 = (t0 - abs_tol < ceil(t0) - 1) ? floor(t0) : ceil(t0);
	// add another day if the last time point is equal to day_max up to absolute tolerance tol
	const auto day_max = (t_max + abs_tol > floor(t_max) + 1) ? ceil(t_max) : floor(t_max);

1. I am actually unsure how the (stochastic) numerical integrator handles reaching the maximum time t_max and if there, a tolerance would be needed. @nijawa Could you check there on your own?
2. A good (without any negative impact) solution would, at least for the interpolation, just be to increase the tolerance in the interpolate function to 1e-10 or am I wrong @reneSchm ?

[reneSchm]

I am actually unsure how the (stochastic) numerical integrator handles reaching the maximum time t_max and if there, a tolerance would be needed

I can answer that, it is the same as the ODE integrator (just that it uses 2 DerivFunctions). So it uses the same relative tolerance for tmax, with a constant 1e-10.

A good (without any negative impact) solution would, at least for the interpolation, just be to increase the tolerance in the interpolate function to 1e-10

I think that would fix it mostly, but maybe we should use our zero_tolerance (1e-12 for double, 1e-6 for float) constant for both the integrator as well as the interpolation default? We could also add a relaxing factor (times 2 or 5, maybe?) for the interpolation, and switch to a relative error on the day0/day_max computation?