@ -1,13 +1,11 @@
# include <vector>
# include <tuple>
# include <cmath>
# include "storm/utility/numerical.h"
# 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 {
@ -15,79 +13,93 @@ namespace storm {
namespace numerical {
template < typename ValueType >
FoxGlynnResult < ValueType > : : FoxGlynnResult ( ) : left ( 0 ) , right ( 0 ) , totalWeight ( 0 ) {
FoxGlynnResult < ValueType > : : FoxGlynnResult ( ) : left ( 0 ) , right ( 0 ) , totalWeight ( storm : : utility : : zero < ValueType > ( ) ) {
// Intentionally left empty.
}
/*!
* The following implementation of Fox and Glynn ' s algorithm is taken from David Jansen ' s patched version
* in MRMC , which is based on his paper :
*
* https : //pms.cs.ru.nl/iris-diglib/src/getContent.php?id=2011-Jansen-UnderstandingFoxGlynn
*
* We have only adapted the code to match more of C + + ' s and our coding guidelines .
*/
template < typename ValueType >
FoxGlynnResult < ValueType > foxGlynnFinder ( ValueType lambda , ValueType epsilon ) {
int left , right ;
FoxGlynn * pFG ;
ValueType tau = std : : numeric_limits < ValueType > : : min ( ) ;
ValueType omega = std : : numeric_limits < ValueType > : : max ( ) ;
ValueType const sqrt_2_pi = boost : : math : : constants : : root_two_pi < ValueType > ( ) ;
ValueType const log10_e = std : : log10 ( boost : : math : : constants : : e < ValueType > ( ) ) ;
uint64_t m = static_cast < uint64_t > ( lambda ) ;
int64_t left = 0 ;
int64_t right = 0 ;
/* tau is only used in underflow checks, which we are going to do in the
logarithm domain . */
// 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.*/
// 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 ) ) ;
// Compute 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 ) {
STORM_LOG_WARN ( " Fox-Glynn: 0 < lambda < 25, underflow near Poi( " < < lambda < < " , 0) = " < < std : : exp ( - lambda ) < < " . The results are unreliable. " ) ;
}
/* zero is used as left truncation point for lambda <= 25 */
// 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 */
} else {
// Compute the left truncation point for lambda >= 25 (for lambda < 25 we use zero as left truncation point).
ValueType const bl = ( 1 + 1 / lambda ) * std : : exp ( ( 1 / lambda ) * 0.125 ) ;
ValueType const sqrt_lambda = std : : sqrt ( lambda ) ;
int64_t k ;
for ( k = 4 ; TRUE ; k + + ) {
double max_err ;
// 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 ; ; + + k ) {
ValueType max_err ;
left = m - static_cast < int64_t > ( std : : ceil ( k * sqrt_lambda + 0.5 ) ) ;
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 ) {
// 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 . ) */
// 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 . */
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*/
// Finally the left truncation point is found.
}
/****Compute pFG->right truncation point****/
// Compute 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 ;
ValueType lambda_max ;
int64_t 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 = 400 ;
m_max = 400 ;
epsilon * = 0.662608824988162441697980 ;
/* i.e. al = (1+1/400) * exp(1/16) * sqrt_2; epsilon /= al; */
@ -97,184 +109,157 @@ namespace storm {
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 ) ) )
// 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 ; ; + + 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));
// actually: "k * (dkl_inv*epsilon/al) >= exp(-0.5 * k^2)", but epsilon has been changed appropriately.
if ( k * epsilon > = 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 ) ;
right = m_max + static_cast < int64_t > ( std : : ceil ( k * std : : sqrt ( 2 * lambda_max ) + 0.5 ) ) ;
if ( right > m_max + static_cast < int64_t > ( std : : ceil ( ( lambda_max + 1 ) * 0.5 ) ) ) {
STORM_LOG_WARN ( " Fox-Glynn: right = " < < right < < " >> lambda = " < < lambda_max < < " , cannot bound the right tail. The results are unreliable. " ) ;
}
}
/*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 ) ) ;
// Time to set the initial value for weights.
FoxGlynnResult < ValueType > fgresult ;
fgresult . left = static_cast < uint64_t > ( left ) ;
fgresult . right = static_cast < uint64_t > ( right ) ;
fgresult . weights . resize ( fgresult . right - fgresult . left + 1 ) ;
fgresult . weights [ m - left ] = omega / ( 1.0e+10 * ( right - left ) ) ;
if ( m > = 25 )
{
/* perform underflow check */
double result , log_c_m_inf ;
int i ;
if ( m > = 25 ) {
// Perform underflow check.
ValueType result , log_c_m_inf ;
int64_t 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 ) . */
// we are going to compare with tau - log(w[m]).
tau - = std : : log ( fgresult . 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. */
// 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 . */
// Equivalent to 2*i <= m, equivalent to i <= lambda/2.
if ( i < = left ) {
// 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 . */
} 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) */
if ( left ! = 0 ) {
// 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 )
// 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 ) ;
if ( result < = tau ) {
int64_t const log10_result = static_cast < int64_t > ( std : : floor ( result * log10_e ) ) ;
STORM_LOG_WARN ( " Fox-Glynn: lambda >= 25, underflow near Poi( " < < lambda < < " , " < < left < < " ) <= " < < std : : exp ( result - log10_result / log10_e ) < < log10_result < < " . The results are unreliable. " ) ;
}
/*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 */
// 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 ) ;
if ( result < = tau ) {
int64_t const log10_result = static_cast < int64_t > ( std : : floor ( result * log10_e ) ) ;
STORM_LOG_WARN ( " Fox-Glynn: lambda >= 25, underflow near Poi( " < < lambda < < " , " < < right < < " ) <= " < < std : : exp ( result - log10_result / log10_e ) < < log10_result < < " . The results are unreliable. " ) ;
}
}
}
return pFG ;
return fgresult ;
}
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 ;
ValueType tau = std : : numeric_limits < ValueType > : : min ( ) ;
// The magic m point.
uint64_t m = static_cast < uint64_t > ( lambda ) ;
int64_t j , t ;
FoxGlynnResult < ValueType > result = foxGlynnFinder ( lambda , epsilon ) ;
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 - result . left ; j > 0 ; - - j ) {
result . weights [ j - 1 ] = ( j + result . left ) / lambda * result . weights [ j ] ;
}
/*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 = result . right - result . left ;
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 {
// 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.
STORM_LOG_THROW ( result . right < = 600 , storm : : exceptions : : PrecisionExceededException , " Fox-Glynn: " < < result . right < < " > 600, underflow is possible. " ) ;
// Compute weights.
for ( j = m - result . left ; j < t ; + + j ) {
ValueType q = lambda / ( j + 1 + result . left ) ;
if ( result . weights [ j ] > tau / q ) {
result . weights [ j + 1 ] = q * result . weights [ j ] ;
} else {
t = j ;
pFG - > right = j + pFG - > left ;
break ; /*It's time to compute W*/
result . right = j + result . left ;
// It's time to compute W.
break ;
}
}
} else {
/*Compute weights*/
for ( j = m - pFG - > left ; j < t ; j + + )
pFG - > weights [ j + 1 ] = lambda / ( j + 1 + pFG - > left ) * pFG - > weights [ j ] ;
} else {
// Compute weights.
for ( j = m - result . left ; j < t ; + + j ) {
result . weights [ j + 1 ] = lambda / ( j + 1 + result . left ) * result . weights [ j ] ;
}
}
/*It is time to compute the normalization weight W*/
pFG - > total_weight = 0.0 ;
// It is time to compute the normalization weight W.
result . totalWeight = storm : : utility : : zero < ValueType > ( ) ;
j = 0 ;
/* t was set above */
while ( j < t )
{
if ( pFG - > weights [ j ] < = pFG - > weights [ t ] )
{
pFG - > total_weight + = pFG - > weights [ j ] ;
// t was set above.
while ( j < t ) {
if ( result . weights [ j ] < = result . weights [ t ] ) {
result . totalWeight + = result . weights [ j ] ;
j + + ;
} else {
pFG - > total_weight + = pFG - > weights [ t ] ;
} else {
result . totalWeight + = result . weights [ t ] ;
t - - ;
}
}
pFG - > total_weight + = pFG - > weights [ j ] ;
result . totalWeight + = result . weights [ j ] ;
/* printf("Fox-Glynn: ltp = %d, rtp = %d, w = %10.15le \n", pFG->left, pFG->right, pFG->total_weight); */
STORM_LOG_TRACE ( " Fox-Glynn: ltp = " < < result . left < < " , rtp = " < < result . right < < " , w = " < < result . totalWeight < < " . " ) ;
return pFG ;
return result ;
}
template < typename ValueType >
dehnert
8 years ago