@@ -175,9 +175,16 @@ evolve_result evolve(Info *info, double duration, int64_t *state, double *rates)
175175 // Find the total for all propensities
176176 total = 0.0 ;
177177 for (reaction = 0 ; reaction < reactions_count ; reaction ++ ) {
178+ if (propensities [reaction ] < 0 ) {
179+ status = 3 ; // a negative propensity
180+ }
178181 total += propensities [reaction ];
179182 }
180183
184+ if (status > 0 ) {
185+ break ;
186+ }
187+
181188 if (isnan (total )) {
182189 printf ("failed simulation: total propensity is NaN\n" );
183190 int max_reaction = 0 ;
@@ -208,6 +215,11 @@ evolve_result evolve(Info *info, double duration, int64_t *state, double *rates)
208215 // `interval` from an exponential distribution.
209216 interval = sample_exponential (random_state , total );
210217 point = sample_uniform (random_state ) * total ;
218+ if (point > total ) {
219+ // If roundoff made point > total, the `progress` loop would go past the
220+ // end of the array.
221+ point = total ;
222+ }
211223
212224 // If we have surpassed the provided duration we can exit now
213225 if (now + interval > duration ) {
@@ -219,6 +231,8 @@ evolve_result evolve(Info *info, double duration, int64_t *state, double *rates)
219231 choice = 0 ;
220232 progress = 0.0 ;
221233
234+ // Note: Even if `point` happens to be 0, this needs to skip 0 propensity
235+ // choices to avoid computing negative counts.
222236 while (progress + propensities [choice ] < point || propensities [choice ] == 0 ) {
223237 progress += propensities [choice ];
224238 choice += 1 ;
0 commit comments