memilio
memilio copied to clipboard
SciCompMod
NAN in multivariant stochastic SEIR Model
Bug description
NAN in compartments in SEIR2V Model
Version
Windows
To reproduce
- build code
- execute sde_seir2v_example
(https://github.com/SciCompMod/memilio/blob/1009-nan-in-multivariant-seir-model/cpp/examples/sde_seir2v.cpp)
Relevant log output
No response
Add any relevant information, e.g. used compiler, screenshots.
No response
Checklist
- [x] Attached labels, especially loc:: or model:: labels.
- [x] Linked to project
I use std::min and std::max instead of std::clamp to ensure that flows are in the correct bonds. Usage of std::clamp gives a "invalid bounds arguments passed to std::clamp" error which is consistent with the NAN output
@nijawa : Yes, nan is probably an incorrect input but clamp should prevent any appearance of nan beforehand. So there's potentially another problem? My questions still is "How does the nan appear?" Can you do a break point before the first nan and check what the actual miscomputation is?
@mknaranja At t = 21.5 we have the compartment InfectedV1V2 at 0.014175747933863577 with an ingoing flow of 0 and an outgoing flow of 0.14175747933863578. With step size = 0.1 that would be 0.014175747933863577 - 0.014175747933863578 in the next step which would be negative. It is displayed as -0.0000 but with higher precision the value seems to be -1.73472347597681e-18. A negative compartment will of course result in nan.
I suspect that the issue is that we dont minimize the flow to the compartment but rather we minimize it to compartment / step_size. compartment / step_size * step_size can differ from compartment in floating point calc.
I dont have a good way to resolve this right now
I think that compartment / step_size * step_size can differ from compartment is not the problem, that basically is the thing for all flops. The problem, as you state, rather is that we do not have a tolerance like eps here. Can you make a suggestion based on @reneSchm's already suggested changes? You can now branch from main.
my seir2v model is not yet in main yet so if I branch from main that model will be missing I think? There also still is a bit I have to tidy up before it should get merged into main.
Is there a prefered eps value for memilio? Else I would propose 1e-6
Yes, the change would just address the functionality not merge the model yet.
As we are in double precision, eps (or I propose tol) should be much smaller. Minimum 1e-10 or rather 1e-14.
Clamping with a tolerance also is problematic if compartments are 0. To fix this I could either use a multiplicative tolerance (something like 1-1e-10) or use std::max(std::min to enforce non-negative flows over clamp errors
I suggest clamping to zero exactly, and only use tolerances on the upper bound. Then you don't need a relative tolerance, at least I would've just subtracted tol. You can test which works better for you.
You should avoid getting NaNs in the first place, min/max will probably not reliably recover the results if there are NaNs involved.
One other thing you could consider is post-processing the results of get_flows in the advance function. There you can just take the max of each compartment and 0. If you've clamped the flows before, the error from that should be at most around 1e-16.
my seir2v model is not yet in main yet so if I branch from main that model will be missing I think?
Depends. If you have nothing commited, then you can just pull again to update all unrelated files. If there are edits that would be overwritten, git will complain. If you did make commits to your branch, you can try rebasing on main, or if that fails merge it. Then there is also the good old method of opening a new branch and copying files manually. Can be easier, occasionally.
I suggest clamping to zero exactly, and only use tolerances on the upper bound. Then you don't need a relative tolerance, at least I would've just subtracted
tol. You can test which works better for you.You should avoid getting NaNs in the first place, min/max will probably not reliably recover the results if there are NaNs involved.
One other thing you could consider is post-processing the results of get_flows in the advance function. There you can just take the max of each compartment and 0. If you've clamped the flows before, the error from that should be at most around 1e-16.
The tolerance works on the upper limit of the clamp. std::clamp(flow, 0, compartment / step_size - tol) doesnt seem to work if compartment = 0
I am not sure. If the problem still persists, it would be possible to catch these errors in the simulation, as the SDE models use their own simulations anyways.