From 39b9fcbe65a968846cdb54685d36107ffab9aef6 Mon Sep 17 00:00:00 2001 From: Joachim Klein Date: Sat, 19 May 2018 19:29:22 +0200 Subject: [PATCH] Fox-Glynn computation: Apply accuracy check also for the 'naive' case (avoid infinite loop) (addresses #72) If the desired accuracy is to small, the naive weight computation can run into an infinite loop, as floating point precision / rounding leads to non-termination of the loop. To partially address this, we move the accuracy check up, so that it is also applied in the case of the 'naive' computation. This should make it much less likely to run into this in practice. For a full fix, we'd either need to check for non-progress in the loop or do an analysis that the floating-point precision of double always suffices for the remaining allowed input values. --- prism/src/explicit/FoxGlynn.java | 13 +++++++++---- prism/src/prism/prism.cc | 17 +++++++++++------ 2 files changed, 20 insertions(+), 10 deletions(-) diff --git a/prism/src/explicit/FoxGlynn.java b/prism/src/explicit/FoxGlynn.java index 98df6add..96725180 100644 --- a/prism/src/explicit/FoxGlynn.java +++ b/prism/src/explicit/FoxGlynn.java @@ -74,7 +74,12 @@ public final class FoxGlynn if (q_tmax == 0.0) { throw new PrismException("Overflow: TA parameter qtmax = time * maxExitRate = 0."); } - else if (q_tmax < 400) + + if (accuracy < 1e-10) { + throw new PrismException("Overflow: Accuracy is smaller than Fox Glynn can handle (must be at least 1e-10)."); + } + + if (q_tmax < 400) { //here naive approach should have better performance than Fox Glynn final double expcoef = Math.exp(-q_tmax); //the "e^-lambda" part of p.m.f. of Poisson dist. int k; //denotes that we work with event "k steps occur" @@ -91,6 +96,9 @@ public final class FoxGlynn //add further steps until you have accumulated enough k = 1; do { + // TODO: catch case where lastval gets so small that + // accnum never reaches desval due to rounding/floating point precision errors (infinite loop) + lastval *= q_tmax / k; // invariant: lastval = q_tmax^k / k! accum += lastval; w.add(lastval * expcoef); @@ -112,9 +120,6 @@ public final class FoxGlynn } else { //use actual Fox Glynn for q_tmax>400 - if (accuracy < 1e-10) { - throw new PrismException("Overflow: Accuracy is smaller than Fox Glynn can handle (must be at least 1e-10)."); - } final double factor = 1e+10; //factor from the paper, it has no real explanation there final int m = (int) q_tmax; //mode //run FINDER to get left, right and weight[m] diff --git a/prism/src/prism/prism.cc b/prism/src/prism/prism.cc index 9de8ddd5..e4aa89d8 100644 --- a/prism/src/prism/prism.cc +++ b/prism/src/prism/prism.cc @@ -116,7 +116,14 @@ EXPORT FoxGlynnWeights fox_glynn(double q_tmax, double underflow, double overflo fgw.right = -1; return fgw; } - else if (q_tmax < 400) + + if (accuracy < 1e-10) { + printf("Overflow: Accuracy is smaller than Fox Glynn can handle (must be at least 1e-10)."); + fgw.right = -1; + return fgw; + } + + if (q_tmax < 400) { //here naive approach should have better performance than Fox Glynn const double expcoef = exp(-q_tmax); //the "e^-lambda" part of p.m.f. of Poisson dist. int k; //denotes that we work with event "k steps occur" @@ -133,6 +140,9 @@ EXPORT FoxGlynnWeights fox_glynn(double q_tmax, double underflow, double overflo //add further steps until you have accumulated enough k = 1; do { + // TODO: catch case where lastval gets so small that + // accnum never reaches desval due to rounding/floating point precision errors (infinite loop) + lastval *= q_tmax / k; // invariant: lastval = q_tmax^k / k! accum += lastval; w.push_back(lastval * expcoef); @@ -154,11 +164,6 @@ EXPORT FoxGlynnWeights fox_glynn(double q_tmax, double underflow, double overflo } else { //use actual Fox Glynn for q_tmax>400 - if (accuracy < 1e-10) { - printf("Overflow: Accuracy is smaller than Fox Glynn can handle (must be at least 1e-10)."); - fgw.right = -1; - return fgw; - } const double factor = 1e+10; const long m = (long) q_tmax; //mode