Negative values in state matrix when stochastic

[ethho]

Problem

On d444bea, negative values appear in E and Ia compartments of state array, only when running model with stochastic draws.

$ tox -- tests/test_models.py
# ...
>           raise ValueError(err)
E           ValueError: Found negative value(s) in DataArray: <xarray.DataArray (vertex: 13, compt: 2, age: 5, risk: 2)>
E           array([[[[nan, nan],
E                    [-2., nan],
E                    [nan, nan],
E                    [nan, nan],
E                    [nan, nan]],
E           
E                   [[nan, nan],
E                    [nan, nan],
E                    [nan, nan],
E                    [nan, nan],
E                    [nan, nan]]],
E           
E           
E                  [[[nan, -1.],
E                    [nan, nan],
E                    [nan, nan],
E                    [nan, nan],
E                    [nan, nan]],
E           
E                   [[nan, nan],
E           ...
E                    [nan, nan]],
E           
E                   [[-1., nan],
E                    [nan, nan],
E                    [nan, nan],
E                    [nan, nan],
E                    [nan, nan]]],
E           
E           
E                  [[[nan, nan],
E                    [-1., nan],
E                    [nan, nan],
E                    [nan, nan],
E                    [nan, nan]],
E           
E                   [[nan, nan],
E                    [nan, nan],
E                    [nan, nan],
E                    [nan, nan],
E                    [nan, nan]]]])
E           Coordinates:
E             * age      (age) <U5 '<5' '5-17' '18-49' '50-64' '65+'
E             * risk     (risk) <U4 'low' 'high'
E             * compt    (compt) <U5 'E' 'Ia'
E             * vertex   (vertex) int32 78644 78648 78652 78654 ... 78741 78747 78752 78753

episimlab/utils/variable.py:64: ValueError
# ...

Root Cause

It seems that some compartments (vertices in the compartment graph) are being depleted beyond zero. The key is that this error occurs only when running the model non-deterministically. I suspected that this was an issue where weights of lower priority edges (e.g. E2Pa has lower priority than S2E) were being applied regardless of the scaling factor k calculated at:

https://github.com/eho-tacc/episimlab/blob/d444bea66ad92a57bc5b5044e165fea0880c51ce/episimlab/compt_model.py#L62-L75

calc_k should account for fully depleted u compartments, since k goes to zero when state of u goes to zero. The problem is that the edge_weight method, in which the stochastic draw is performed, is called during calculation of k, and again when determining the applied edge weight:

https://github.com/eho-tacc/episimlab/blob/d444bea66ad92a57bc5b5044e165fea0880c51ce/episimlab/compt_model.py#L83-L83

When edge_weight is non-deterministic, the returned weight is almost certainly different when the same function is called with same u, v args multiple times. Consequently, the result of the first stochastic draw is sometimes greater than the result of the second draw, resulting in a k that is not sufficiently small to prevent depletion when the edge weight returned by the second draw is applied.

[ethho]

Solution

Implemented on 6493908. Simply cache calls to edge_weight (now called get_edge_weight) in a dictionary that clears at every timestep:

https://github.com/eho-tacc/episimlab/blob/5fd7bc1e0511a2804955078be10a4d0756b40942/episimlab/compt_model.py#L97-L102

This ensures that we draw at most once per edge per timestep.