diff --git a/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.cpp b/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.cpp index 2d270c720..90d556133 100644 --- a/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.cpp +++ b/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.cpp @@ -634,19 +634,19 @@ namespace storm { // Use Fox-Glynn to get the truncation points and the weights. // std::tuple> foxGlynnResult = storm::utility::numerical::getFoxGlynnCutoff(lambda, 1e+300, storm::settings::getModule().getPrecision() / 8.0); - std::tuple> foxGlynnResult = storm::utility::numerical::foxGlynn(lambda, storm::settings::getModule().getPrecision() / 8.0); - STORM_LOG_DEBUG("Fox-Glynn cutoff points: left=" << std::get<0>(foxGlynnResult) << ", right=" << std::get<1>(foxGlynnResult)); + storm::utility::numerical::FoxGlynnResult foxGlynnResult = storm::utility::numerical::foxGlynn(lambda, storm::settings::getModule().getPrecision() / 8.0); + STORM_LOG_DEBUG("Fox-Glynn cutoff points: left=" << foxGlynnResult.left << ", right=" << foxGlynnResult.right); // Scale the weights so they add up to one. - for (auto& element : std::get<3>(foxGlynnResult)) { - element /= std::get<2>(foxGlynnResult); + for (auto& element : foxGlynnResult.weights) { + element /= foxGlynnResult.totalWeight; } // If the cumulative reward is to be computed, we need to adjust the weights. if (useMixedPoissonProbabilities) { ValueType sum = storm::utility::zero(); - for (auto& element : std::get<3>(foxGlynnResult)) { + for (auto& element : foxGlynnResult.weights) { sum += element; element = (1 - sum) / uniformizationRate; } @@ -656,10 +656,10 @@ namespace storm { // Initialize result. std::vector result; - uint_fast64_t startingIteration = std::get<0>(foxGlynnResult); + uint_fast64_t startingIteration = foxGlynnResult.left; if (startingIteration == 0) { result = values; - storm::utility::vector::scaleVectorInPlace(result, std::get<3>(foxGlynnResult)[0]); + storm::utility::vector::scaleVectorInPlace(result, foxGlynnResult.weights.front()); ++startingIteration; } else { if (useMixedPoissonProbabilities) { @@ -674,9 +674,9 @@ namespace storm { std::unique_ptr> solver = linearEquationSolverFactory.create(env, std::move(uniformizedMatrix), storm::solver::LinearEquationSolverTask::Multiply); solver->setCachingEnabled(true); - if (!useMixedPoissonProbabilities && std::get<0>(foxGlynnResult) > 1) { + if (!useMixedPoissonProbabilities && foxGlynnResult.left > 1) { // Perform the matrix-vector multiplications (without adding). - solver->repeatedMultiply(values, addVector, std::get<0>(foxGlynnResult) - 1); + solver->repeatedMultiply(values, addVector, foxGlynnResult.left - 1); } else if (useMixedPoissonProbabilities) { std::function addAndScale = [&uniformizationRate] (ValueType const& a, ValueType const& b) { return a + b / uniformizationRate; }; @@ -691,10 +691,10 @@ namespace storm { // multiplication, scale and add the result. ValueType weight = 0; std::function addAndScale = [&weight] (ValueType const& a, ValueType const& b) { return a + weight * b; }; - for (uint_fast64_t index = startingIteration; index <= std::get<1>(foxGlynnResult); ++index) { + for (uint_fast64_t index = startingIteration; index <= foxGlynnResult.right; ++index) { solver->repeatedMultiply(values, addVector, 1); - weight = std::get<3>(foxGlynnResult)[index - std::get<0>(foxGlynnResult)]; + weight = foxGlynnResult.weights[index - foxGlynnResult.left]; storm::utility::vector::applyPointwise(result, values, result, addAndScale); } diff --git a/src/storm/utility/numerical.cpp b/src/storm/utility/numerical.cpp index 1143d8bfd..cec474ff7 100644 --- a/src/storm/utility/numerical.cpp +++ b/src/storm/utility/numerical.cpp @@ -1,13 +1,11 @@ -#include -#include -#include +#include "storm/utility/numerical.h" +#include #include #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 - FoxGlynnResult::FoxGlynnResult() : left(0), right(0), totalWeight(0) { + FoxGlynnResult::FoxGlynnResult() : left(0), right(0), totalWeight(storm::utility::zero()) { // 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 FoxGlynnResult foxGlynnFinder(ValueType lambda, ValueType epsilon) { - int left, right; - FoxGlynn * pFG; + ValueType tau = std::numeric_limits::min(); + ValueType omega = std::numeric_limits::max(); + ValueType const sqrt_2_pi = boost::math::constants::root_two_pi(); + ValueType const log10_e = std::log10(boost::math::constants::e()); + + uint64_t m = static_cast(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(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(std::ceil(k * std::sqrt(2 * lambda_max) + 0.5)); + if (right > m_max + static_cast(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 fgresult; + fgresult.left = static_cast(left); + fgresult.right = static_cast(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(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(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 FoxGlynnResult 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::min(); + + // The magic m point. + uint64_t m = static_cast(lambda); + int64_t j, t; + + FoxGlynnResult 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(); 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 diff --git a/src/storm/utility/numerical.h b/src/storm/utility/numerical.h index 233bb8a74..7bd18cebe 100644 --- a/src/storm/utility/numerical.h +++ b/src/storm/utility/numerical.h @@ -13,7 +13,7 @@ namespace storm { uint64_t right; ValueType totalWeight; std::vector weights; - } + }; template FoxGlynnResult foxGlynn(ValueType lambda, ValueType epsilon);