sp
/
tempest
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

started c++ifying David Jansen's implementation of Fox-Glynn

tempestpy_adaptions
dehnert 7 years ago
parent
48b0a40d8a
commit
27558e2140
2 changed files with 298 additions and 264 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 290
      src/storm/utility/numerical.cpp
  2. 272
      src/storm/utility/numerical.h

290
src/storm/utility/numerical.cpp
View File

@ -0,0 +1,290 @@
#include <vector>
#include <tuple>
#include <cmath>
#include <boost/math/constants/constants.hpp>
#include "storm/utility/macros.h"
#include "storm/utility/constants.h"
#include "storm/exceptions/InvalidArgumentException.h"
#include "storm/exceptions/OutOfRangeException.h"
#include "storm/exceptions/PrecisionExceededException.h"
namespace storm {
namespace utility {
namespace numerical {
template<typename ValueType>
FoxGlynnResult<ValueType>::FoxGlynnResult() : left(0), right(0), totalWeight(0) {
// Intentionally left empty.
}
template<typename ValueType>
FoxGlynnResult<ValueType> foxGlynnFinder(ValueType lambda, ValueType epsilon) {
int left, right;
FoxGlynn * pFG;
/* tau is only used in underflow checks, which we are going to do in the
logarithm domain. */
tau = log(tau);
/* In error bound comparisons, we always compare with epsilon*sqrt_2_pi.*/
epsilon *= sqrt_2_pi;
/****Compute pFG->left truncation point****/
if ( m < 25 )
{
/* for lambda below 25 the exponential can be smaller than tau */
/* if that is the case we expect underflows and warn the user */
if ( -lambda <= tau )
{
printf("ERROR: Fox-Glynn: 0 < lambda < 25, underflow near Poi(%g,"
"0) = %.2e. The results are UNRELIABLE.\n",
lambda, exp(-lambda));
}
/* zero is used as left truncation point for lambda <= 25 */
left = 0;
}
/* compute the left truncation point for lambda >= 25 */
/* for lambda < 25 we use zero as left truncation point */
else
{
const double bl = (1 + 1 / lambda) * exp((1/lambda) * 0.125);
const double sqrt_lambda = sqrt(lambda);
int k;
/*Start looking for the left truncation point*/
/* start search at k=4 (taken from original Fox-Glynn paper) */
/* increase the left truncation point until we fulfil the error
condition */
for ( k = 4 ; TRUE ; k++ ) {
double max_err;
left = m - (int) ceil(k*sqrt_lambda + 0.5);
/* for small lambda the above calculation can yield negative truncation points, crop them here */
if ( left <= 0 ) {
left = 0;
break;
}
/* Note that Propositions 2-4 in Fox--Glynn mix up notation: they
write Phi where they mean 1 - Phi. (In Corollaries 1 and 2, phi is
used correctly again.) */
max_err = bl * exp(-0.5 * (k*k)) / k;
if ( max_err * 2 <= epsilon ) {
/* If the error on the left hand side is smaller, we can be
more lenient on the right hand side. To this end, we now set
epsilon to the part of the error that has not yet been eaten
up by the left-hand truncation. */
epsilon -= max_err;
break;
}
}
/*Finally the left truncation point is found*/
}
/****Compute pFG->right truncation point****/
{
double lambda_max;
int m_max, k;
/*According to Fox-Glynn, if lambda < 400 we should take lambda = 400,
otherwise use the original value. This is for computing the right truncation point*/
if ( m < 400 ) {
lambda_max = lambda_400;
m_max = 400;
epsilon *= 0.662608824988162441697980;
/* i.e. al = (1+1/400) * exp(1/16) * sqrt_2; epsilon /= al; */
} else {
lambda_max = lambda;
m_max = m;
epsilon *= (1 - 1 / (lambda + 1)) * 0.664265347050632847802225;
/* i.e. al = (1+1/lambda) * exp(1/16) * sqrt_2; epsilon /= al; */
}
/* find right truncation point */
/* This loop is a modification to the original Fox-Glynn paper.
The search for the right truncation point is only terminated by
the error condition and not by the stop index from the FG paper.
This can yield more accurate results if necessary.*/
for ( k = 4 ; TRUE ; k++ )
{
/* dkl_inv is between 1 - 1e-33 and 1 if lambda_max >= 400 and
k >= 4; this will always be rounded to 1.0. We therefore leave the
factor out.
double dkl_inv=1 - exp(-266/401.0 * (k*sqrt(2*lambda_max) + 1.5));
*/
if ( k * epsilon /* actually: "k * (dkl_inv*epsilon/al)", which
has been calculated above */ >= exp(-0.5*(k*k)) )
break;
}
right = m_max + (int) ceil(k * sqrt(2 * lambda_max) + 0.5);
if ( right > m_max + (int) ceil((lambda_max + 1) * 0.5) ) {
printf("ERROR: Fox-Glynn: right = %d >> lambda = %g, cannot "
"bound the right tail. The results are "
"UNRELIABLE.\n", right, lambda_max);
}
}
/*Time to set the initial value for weights*/
pFG = calloc((size_t) 1, sizeof(FoxGlynn) +
(right - left) * sizeof(pFG->weights[0]));
if ( NULL == pFG ) {
err_msg_3(err_MEMORY, "finder(%d,%g,_,%g,_)", m, lambda, omega, NULL);
}
pFG->right = right;
pFG->left = left;
pFG->weights[m - left] = omega / ( 1.0e+10 * (right - left) );
if ( m >= 25 )
{
/* perform underflow check */
double result, log_c_m_inf;
int i;
/* we are going to compare with tau - log(w[m]) */
tau -= log(pFG->weights[m - left]);
/*We take the c_m_inf = 0.14627 / sqrt( m ), as for lambda >= 25
c_m = 1 / ( sqrt( 2.0 * pi * m ) ) * exp( m - lambda - 1 / ( 12.0 * m ) ) => c_m_inf*/
/* Note that m-lambda is in the interval (-1,0],
and -1/(12*m) is in [-1/(12*25),0).
So, exp(m-lambda - 1/(12*m)) is in (exp(-1-1/(12*25)),exp(0)).
Therefore, we can improve the lower bound on c_m to
exp(-1-1/(12*25)) / sqrt(2*pi) = ~0.14627. Its logarithm is
-1 - 1/(12*25) - log(2*pi) * 0.5 = ~ -1.922272 (rounded towards
-infinity). */
log_c_m_inf = -1.922272 - log((double) m) * 0.5;
/* We use FG's Proposition 6 directly (and not Corollary 4 i and ii),
as k_prime may be too large if pFG->left == 0. */
i = m - left;
if ( i <= left /* equivalent to 2*i <= m,
equivalent to i <= lambda/2 */ )
{
/* Use Proposition 6 (i). Note that Fox--Glynn are off by one in
the proof of this proposition; they sum up to i-1, but should have
summed up to i. */
result = log_c_m_inf
- i * (i+1) * (0.5 + (2*i+1)/(6*lambda)) / lambda;
}
else
{
/* Use Corollary 4 (iii). Note that k_prime <= sqrt(m+1)/m is a
misprint for k_prime <= m/sqrt(m+1), which is equivalent to
left >= 0, which holds trivially. */
result = -lambda;
if ( 0 != left ) {
/* also use Proposition 6 (ii) */
double result_1 = log_c_m_inf + i * log(1 - i/(double) (m+1));
/*Take the maximum*/
if ( result_1 > result )
result = result_1;
}
}
if ( result <= tau )
{
const int log10_result = (int) floor(result * log10_e);
printf("ERROR: Fox-Glynn: lambda >= 25, underflow near Poi(%g,%d)"
" <= %.2fe%+d. The results are UNRELIABLE.\n",
lambda, left, exp(result - log10_result/log10_e),
log10_result);
}
/*We still have to perform an underflow check for the right truncation point when lambda >= 400*/
if ( m >= 400 )
{
/* Use Proposition 5 of Fox--Glynn */
i = right - m;
result = log_c_m_inf - i * (i + 1) / (2 * lambda);
if ( result <= tau )
{
const int log10_result = (int) floor(result * log10_e);
printf("ERROR: Fox-Glynn: lambda >= 400, underflow near "
"Poi(%g,%d) <= %.2fe%+d. The results are "
"UNRELIABLE.\n", lambda, right,
exp(result - log10_result / log10_e),
log10_result);
}
}
}
return pFG;
}
template<typename ValueType>
FoxGlynnResult<ValueType> foxGlynnWeighter(ValueType lambda, ValueType epsilon) {
/*The magic m point*/
const uint64_t m = (int) floor(lambda);
int j, t;
FoxGlynn * pFG;
pFG = finder(m, lambda, tau, omega, epsilon);
if ( NULL == pFG ) {
err_msg_4(err_CALLBY, "weighter(%g,%g,%g,%g)", lambda,
tau, omega, epsilon, NULL);
}
/*Fill the left side of the array*/
for ( j = m - pFG->left ; j > 0 ; j-- )
pFG->weights[j-1] = (j+pFG->left) / lambda * pFG->weights[j];
t = pFG->right - pFG->left;
/*Fill the right side of the array, have two cases lambda < 400 & lambda >= 400*/
if ( m < 400 )
{
/*Perform the underflow check, according to Fox-Glynn*/
if( pFG->right > 600 )
{
printf("ERROR: Fox-Glynn: pFG->right > 600, underflow is possible\n");
freeFG(pFG);
return NULL;
}
/*Compute weights*/
for ( j = m - pFG->left ; j < t ; j++ )
{
double q = lambda / (j + 1 + pFG->left);
if ( pFG->weights[j] > tau / q )
{
pFG->weights[j + 1] = q * pFG->weights[j];
}else{
t = j;
pFG->right = j + pFG->left;
break; /*It's time to compute W*/
}
}
}else{
/*Compute weights*/
for ( j = m - pFG->left ; j < t ; j++ )
pFG->weights[j+1] = lambda / (j+1 + pFG->left) * pFG->weights[j];
}
/*It is time to compute the normalization weight W*/
pFG->total_weight = 0.0;
j = 0;
/* t was set above */
while( j < t )
{
if ( pFG->weights[j] <= pFG->weights[t] )
{
pFG->total_weight += pFG->weights[j];
j++;
}else{
pFG->total_weight += pFG->weights[t];
t--;
}
}
pFG->total_weight += pFG->weights[j];
/* printf("Fox-Glynn: ltp = %d, rtp = %d, w = %10.15le \n", pFG->left, pFG->right, pFG->total_weight); */
return pFG;
}
template<typename ValueType>
FoxGlynnResult<ValueType> foxGlynn(ValueType lambda, ValueType epsilon) {
STORM_LOG_THROW(lambda > 0, storm::exceptions::InvalidArgumentException, "Fox-Glynn requires positive lambda.");
return foxGlynnWeighter(lambda, epsilon);
}
template FoxGlynnResult<double> foxGlynn(double lambda, double epsilon);
}
}
}

272
src/storm/utility/numerical.h
View File

@ -1,278 +1,22 @@
#include <vector>
#include <tuple>
#include <cmath>
#include <boost/math/constants/constants.hpp>
#include "storm/utility/macros.h"
#include "storm/utility/constants.h"
#include "storm/exceptions/InvalidArgumentException.h"
#include "storm/exceptions/OutOfRangeException.h"
#include "storm/exceptions/PrecisionExceededException.h"
#include <cstdint>
namespace storm {
namespace utility {
namespace numerical {
template<typename ValueType>
std::tuple<uint_fast64_t, uint_fast64_t, ValueType, std::vector<ValueType>> getFoxGlynnCutoff(ValueType lambda, ValueType overflow, ValueType accuracy) {
STORM_LOG_THROW(lambda != storm::utility::zero<ValueType>(), storm::exceptions::InvalidArgumentException, "Error in Fox-Glynn algorithm: lambda must not be zero.");
// This code is a modified version of the one in PRISM. According to their implementation, for lambda
// smaller than 400, we compute the result using the naive method.
if (lambda < 400) {
ValueType eToPowerMinusLambda = std::exp(-lambda);
ValueType targetValue = (1 - accuracy) / eToPowerMinusLambda;
std::vector<ValueType> weights;
ValueType exactlyKEvents = 1;
ValueType atMostKEvents = exactlyKEvents;
weights.push_back(exactlyKEvents * eToPowerMinusLambda);
uint_fast64_t k = 1;
do {
exactlyKEvents *= lambda / k;
atMostKEvents += exactlyKEvents;
weights.push_back(exactlyKEvents * eToPowerMinusLambda);
++k;
} while (atMostKEvents < targetValue);
return std::make_tuple(0, k - 1, 1.0, weights);
} else {
STORM_LOG_THROW(accuracy >= 1e-10, storm::exceptions::InvalidArgumentException, "Error in Fox-Glynn algorithm: the accuracy must not be below 1e-10.");
// Factor from Fox&Glynn's paper. The paper does not explain where it comes from.
ValueType factor = 1e+10;
// Now start the Finder algorithm to find the truncation points.
ValueType m = std::floor(lambda);
uint_fast64_t leftTruncationPoint = 0, rightTruncationPoint = 0;
{
// Factors used by the corollaries explained in Fox & Glynns paper.
// Square root of pi.
ValueType sqrtpi = 1.77245385090551602729;
// Square root of 2.
ValueType sqrt2 = 1.41421356237309504880;
// Set up a_\lambda, b_\lambda, and the square root of lambda.
ValueType aLambda = 0, bLambda = 0, sqrtLambda = 0;
if (m < 400) {
sqrtLambda = std::sqrt(400.0);
aLambda = (1.0 + 1.0 / 400.0) * std::exp(0.0625) * sqrt2;
bLambda = (1.0 + 1.0 / 400.0) * std::exp(0.125 / 400.0);
} else {
sqrtLambda = std::sqrt(lambda);
aLambda = (1.0 + 1.0 / lambda) * std::exp(0.0625) * sqrt2;
bLambda = (1.0 + 1.0 / lambda) * std::exp(0.125 / lambda);
}
// Use Corollary 1 from the paper to find the right truncation point.
uint_fast64_t k = 4;
ValueType dkl = 1.0 / (1 - std::exp(-(2.0 / 9.0) * (k * sqrt2 * sqrtLambda + 1.5)));
// According to David Jansen the upper bound can be ignored to achieve more accurate results.
// Right hand side of the equation in Corollary 1.
while ((accuracy / 2.0) < (aLambda * dkl * std::exp(-(k*k / 2.0)) / (k * sqrt2 * sqrtpi))) {
++k;
// d(k,Lambda) from the paper.
dkl = 1.0 / (1 - std::exp(-(2.0 / 9.0)*(k * sqrt2 * sqrtLambda + 1.5)));
}
// Left hand side of the equation in Corollary 1.
rightTruncationPoint = static_cast<uint_fast64_t>(std::ceil((m + k * sqrt2 * sqrtLambda + 1.5)));
// Use Corollary 2 to find left truncation point.
k = 4;
// Right hand side of the equation in Corollary 2.
while ((accuracy / 2.0) < ((bLambda * std::exp(-(k*k / 2.0))) / (k * sqrt2 * sqrtpi))) {
++k;
}
// Left hand side of the equation in Corollary 2.
leftTruncationPoint = static_cast<uint_fast64_t>(std::max(std::trunc(m - k * sqrtLambda - 1.5), storm::utility::zero<ValueType>()));
// Check for underflow.
STORM_LOG_THROW(std::trunc((m - k * sqrtLambda - 1.5)) >= 0, storm::exceptions::OutOfRangeException, "Error in Fox-Glynn algorithm: Underflow of left truncation point.");
}
std::vector<ValueType> weights(rightTruncationPoint - leftTruncationPoint + 1);
weights[m - leftTruncationPoint] = overflow / (factor * (rightTruncationPoint - leftTruncationPoint));
for (uint_fast64_t j = m; j > leftTruncationPoint; --j) {
weights[j - 1 - leftTruncationPoint] = (j / lambda) * weights[j - leftTruncationPoint];
}
for (uint_fast64_t j = m; j < rightTruncationPoint; ++j) {
weights[j + 1 - leftTruncationPoint] = (lambda / (j + 1)) * weights[j - leftTruncationPoint];
}
// Compute the total weight and start with smallest to avoid roundoff errors.
ValueType totalWeight = storm::utility::zero<ValueType>();
uint_fast64_t s = leftTruncationPoint;
uint_fast64_t t = rightTruncationPoint;
while (s < t) {
if (weights[s - leftTruncationPoint] <= weights[t - leftTruncationPoint]) {
totalWeight += weights[s - leftTruncationPoint];
++s;
} else {
totalWeight += weights[t - leftTruncationPoint];
--t;
}
}
totalWeight += weights[s - leftTruncationPoint];
return std::make_tuple(leftTruncationPoint, rightTruncationPoint, totalWeight, weights);
}
}
template<typename ValueType>
std::tuple<uint64_t, uint64_t, ValueType> foxGlynnFinder(ValueType lambda, ValueType underflow, ValueType overflow, ValueType accuracy) {
uint64_t m = static_cast<uint64_t>(lambda);
int64_t L = 0;
uint64_t R = 0;
uint64_t k = 4;
ValueType pi = boost::math::constants::pi<ValueType>();
if (m < 25) {
STORM_LOG_THROW(std::exp(-static_cast<double>(m)) >= underflow, storm::exceptions::PrecisionExceededException, "Underflow in Fox-Glynn.");
L = 0;
} else {
ValueType rhsb = (accuracy / 2) * std::sqrt(2 * pi) / ((1 + 1/lambda) * std::exp(1 / (8 * lambda)));
k = 4;
do {
L = m - std::ceil(k * std::sqrt(lambda) + 0.5);
if (L <= 0) {
break;
}
if (std::exp(-(k*k) / 2) / k <= rhsb) {
break;
}
++k;
} while (true);
}
uint64_t mmax;
uint64_t maxr;
ValueType r2l;
ValueType rhsa;
if (m < 400) {
mmax = 400;
r2l = std::sqrt(2 * 400);
rhsa = (accuracy / 2) * std::sqrt(2 * pi) / ((1 + 1.0/400) * std::sqrt(2.0) * std::exp(1.0 / 16));
maxr = 400 + static_cast<uint64_t>(std::ceil((400 + 1) / 2.0));
} else {
mmax = m;
r2l = std::sqrt(2 * lambda);
rhsa = (accuracy / 2) * std::sqrt(2 * pi) / ((1 + 1.0/lambda) * std::sqrt(2.0) * std::exp(1.0 / 16));
maxr = m + static_cast<uint64_t>(std::ceil((lambda + 1) / 2.0));
}
k = 4;
do {
R = mmax + static_cast<uint64_t>(std::ceil(k * r2l + 0.5));
STORM_LOG_THROW(R <= maxr, storm::exceptions::PrecisionExceededException, "Fox-Glynn: cannot bound right tail.");
ValueType dkl = storm::utility::one<ValueType>();
if (dkl * std::exp(-(k * k / 2.0)) / k <= rhsa) {
break;
}
++k;
} while (true);
ValueType wm = std::pow<ValueType>(10, -10) * overflow * (R - L + 1);
if (m >= 25) {
ValueType lcm = -1 - (1.0/12*25) - std::log(2 * pi) - 0.5 * std::log(static_cast<ValueType>(m));
uint64_t i = m - L;
ValueType llb;
if (i <= static_cast<uint64_t>(L)) {
llb = -static_cast<ValueType>(i)*(i + 1)*(((2 * i + 1) / (6 * lambda)) + 0.5) * (1 / lambda) + lcm;
} else {
llb = std::max(i * std::log(1 - (i / (m + 1.0))), -lambda);
}
STORM_LOG_THROW(llb >= std::log(underflow) - std::log(wm), storm::exceptions::PrecisionExceededException, "Fox-Glynn: Underflow at lower bound.");
if (m >= 400) {
llb = lcm - std::pow(R - m + 1, 2) / (2 * lambda);
STORM_LOG_THROW(llb >= std::log(underflow) - std::log(wm), storm::exceptions::PrecisionExceededException, "Fox-Glynn: Underflow at upper bound.");
}
}
return std::make_tuple(static_cast<uint64_t>(L), R, wm);
}
template<typename ValueType>
std::tuple<uint64_t, uint64_t, ValueType, std::vector<ValueType>> foxGlynnWeighter(ValueType lambda, ValueType accuracy) {
STORM_LOG_THROW(lambda > storm::utility::zero<ValueType>(), storm::exceptions::InvalidArgumentException, "Fox-Glynn: Lambda must be positive.");
ValueType underflow = std::numeric_limits<ValueType>::min();
ValueType overflow = std::numeric_limits<ValueType>::max();
struct FoxGlynnResult {
FoxGlynnResult();
// Initialize.
uint64_t m = static_cast<uint64_t>(lambda);
std::tuple<uint64_t, uint64_t, ValueType> finderResult = foxGlynnFinder(lambda, underflow, overflow, accuracy);
uint64_t L = std::get<0>(finderResult);
uint64_t R = std::get<1>(finderResult);
ValueType wm = std::get<2>(finderResult);
std::vector<ValueType> w(R - L + 1);
w.back() = wm;
// Down.
for (uint64_t j = m; j > L; --j) {
w[j - L - 1] = (j / lambda) * w[j - L];
}
// Up.
if (lambda < 400) {
// Special.
STORM_LOG_THROW(R <= 600, storm::exceptions::PrecisionExceededException, "Potential underflow in Fox-Glynn.");
for (uint64_t j = m; j < R; ++j) {
auto q = lambda / (j + 1);
if (w[j - L] > underflow / q) {
w[j - L + 1] = q * w[j - L];
} else {
R = j;
break;
}
}
} else {
for (uint64_t j = m; j < R; ++j) {
w[j - L + 1] = (lambda / (j + 1)) * w[j - L];
}
}
// Compute W.
ValueType W = storm::utility::zero<ValueType>();
uint64_t s = L;
uint64_t t = R;
while (s < t) {
if (w[s - L] <= w[t - L]) {
W += w[s - L];
++s;
} else {
W += w[t - L];
--t;
}
}
W += w[s - L];
// Return.
return std::make_tuple(L, R, W, w);
uint64_t left;
uint64_t right;
ValueType totalWeight;
std::vector<ValueType> weights;
}
template<typename ValueType>
std::tuple<uint_fast64_t, uint_fast64_t, ValueType, std::vector<ValueType>> foxGlynn(ValueType lambda, ValueType accuracy) {
return foxGlynnWeighter(lambda, accuracy);
}
FoxGlynnResult<ValueType> foxGlynn(ValueType lambda, ValueType epsilon);
}
}

Write Preview
Loading…
Cancel
Save