diff --git a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp index d72eeee13..529ce05b0 100644 --- a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp +++ b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp @@ -219,20 +219,6 @@ namespace storm { } } - template - ValueType SparseMarkovAutomatonCslHelper::poisson(ValueType lambda, uint64_t i) { - ValueType res = pow(lambda, i); - ValueType fac = 1; - for (uint64_t j = i ; j>0 ; j--){ - fac = fac *j; - } - res = res / fac ; - res = res * exp(-lambda); - return res; - } - - - template ::SupportsExponential, int>::type> void SparseMarkovAutomatonCslHelper::calculateVu(Environment const& env, std::vector> const& relativeReachability, OptimizationDirection dir, uint64_t k, uint64_t node, uint64_t const kind, ValueType lambda, uint64_t probSize, std::vector>>& unifVectors, storm::storage::SparseMatrix const& fullTransitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, std::unique_ptr> const& solver, std::ofstream& logfile, std::vector const& poisson){ if (unifVectors[1][k][node]!=-1){return;} //dynamic programming. avoiding multiple calculation. @@ -245,7 +231,9 @@ namespace storm { calculateUnifPlusVector(env, N-1-(i-k),node,2,lambda,probSize,relativeReachability,dir,unifVectors,fullTransitionMatrix, markovianStates,psiStates,solver, logfile, poisson); //old: relativeReachability, dir, (N-1-(i-k)),node,lambda,wu,fullTransitionMatrix,markovianStates,psiStates, solver); } - res+=poisson[i]*unifVectors[2][N-1-(i-k)][node]; + if (i(); for (uint64_t i = k ; i1e-300){ + if (i> foxGlynnResult = storm::utility::numerical::getFoxGlynnCutoff( - T * lambda, 1e+300, epsilon * epsilon * kappa); + storm::utility::numerical::FoxGlynnResult foxGlynnResult = storm::utility::numerical::foxGlynn(lambda*T, epsilon*kappa); // 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; } -*/ - std::vector poisson = foxGlynnProb(lambda*T, N+1, epsilon*kappa); // (4) define vectors/matrices std::vector init(numberOfStates, -1); @@ -693,13 +678,13 @@ namespace storm { for (uint64_t k = N; k <= N; k--) { calculateUnifPlusVector(env, k, i, 0, lambda, probSize, relReachability, dir, unifVectors, fullTransitionMatrix, markovianStates, psiStates, solver, logfile, - poisson); + foxGlynnResult.weights); calculateUnifPlusVector(env, k, i, 2, lambda, probSize, relReachability, dir, unifVectors, fullTransitionMatrix, markovianStates, psiStates, solver, logfile, - poisson); + foxGlynnResult.weights); calculateVu(env, relReachability, dir, k, i, 1, lambda, probSize, unifVectors, fullTransitionMatrix, markovianStates, psiStates, solver, logfile, - poisson); + foxGlynnResult.weights); //also use iteration to keep maxNorm of vd and vup to date, so the loop-condition is easy to prove ValueType diff = std::abs(unifVectors[0][k][i] - unifVectors[1][k][i]); maxNorm = std::max(maxNorm, diff); diff --git a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h index bff568487..2b62d244f 100644 --- a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h +++ b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h @@ -88,330 +88,6 @@ namespace storm { //TODO: move this - typedef struct FoxGlynn - { - int left; - int right; - double total_weight; - double *weights; - } FoxGlynn; - - static std::vector foxGlynnProb(double lambdaT, int N, double precision){ - - FoxGlynn* fg = NULL; - if(!fox_glynn(lambdaT, DBL_MIN, DBL_MAX, precision, &fg)) { - printf("ERROR: fox-glynn failed\n"); - return std::vector{}; - } - - long double sumOfPoissonProbs = 0.0; - std::vector poisson_p(N,0.0); - unsigned long iter_num; - - std::cout << "fg left " << fg->left << " fh right " << fg->right <<"\n"; - //for(int i=fg->left; i<=fg->right; i++) { - for (int i = 0; iweights[i-fg->left]/fg->total_weight; - sumOfPoissonProbs+=poisson_p[i]; - } - - for(int i=fg->left-1; i>=0; i--) { - poisson_p[i] = poisson_p[i+1]*((i+1)/(lambdaT)); - sumOfPoissonProbs+=poisson_p[i]; - } - - iter_num = fg->right; - - freeFG(fg); - - return poisson_p; - } - - static bool finder(const int m, const double lambda, const double tau, const double omega, - const double epsilon, double * pw_m, FoxGlynn *pFG) - { - /*The pi constant*/ - static const double pi = 3.14159265358979323846264; - static const double lambda_25 = 25.0; - static const double lambda_400 = 40; - - const double sqrt_2_pi = sqrt( 2.0 * pi ); - const double sqrt_2 = sqrt(2.0); - const double sqrt_lambda = sqrt(lambda); - double lambda_max, k, k_rtp = HUGE_VAL, k_prime, c_m_inf, result, al, dkl, bl; - - /*Simple bad cases, when we quit*/ - if( lambda == 0.0 ) - { - printf("ERROR: Fox-Glynn: lambda = 0, terminating the algorithm\n"); - return false; - } - /* The requested error level must not be smaller than the minimum machine precision - (needed to guarantee the convergence of the error conditions) */ - if( epsilon < tau) - { - printf("ERROR: Fox-Glynn: epsilon < tau, invalid error level, terminating the algorithm\n"); - printf("epsilon %f, tau %f\n",epsilon,tau); - return false; - } - /* zero is used as left truncation point for lambda <= 25 */ - pFG->left = 0; - lambda_max = lambda; - - /* for lambda below 25 the exponential can be smaller than tau */ - /* if that is the case we expect underflows and warn the user */ - if( 0.0 < lambda && lambda <= lambda_25 ) - { - if( exp( -lambda ) <= tau ) - { - printf("ERROR: Fox-Glynn: 0 < lambda < 25, underflow. The results are UNRELIABLE.\n"); - } - } - - bl = (1.0 + 1.0/lambda) * exp(1.0 / (8.0 * lambda)); - - /****Compute pFG->right truncation point****/ - /*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(lambda < lambda_400) - lambda_max = lambda_400; - k = 4; - al = (1.0+1.0/lambda_max) * exp(1.0/16.0) * sqrt_2; - dkl = exp(-2.0/9.0 * (k * sqrt(2.0 * lambda_max) + 1.5 )); - dkl = 1.0 / (1.0 - dkl); - /* 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 neccesary.*/ - while((epsilon/2.0) < ((al * dkl * exp(-(k*k)/2.0))/(k*sqrt_2_pi))) - { - k++; - dkl = exp(-2.0/9.0 * (k * sqrt_2 * sqrt(lambda_max) + 1.5 )); - dkl = 1.0 / (1.0 - dkl); - } - k_rtp = k; - pFG->right = (int)ceil(m + k_rtp * sqrt_2 * sqrt(lambda_max) + 1.5 ); - - - /****Compute pFG->left truncation point****/ - /* compute the left truncation point for lambda > 25 */ - /* for lambda <= 25 we use zero as left truncation point */ - if(lambda > lambda_25) - { - /*Start looking for the left truncation point*/ - /* start search at k=4 (taken from original Fox-Glynn paper */ - k = 4; - /* increase the left truncation point as long as we fulfill the error condition */ - - /* This loop is a modification to the original Fox-Glynn paper. - The search for the left 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 neccesary.*/ - while((epsilon/2.0) < ((bl * exp(-(k*k)/2.0))/(k * sqrt_2_pi))) - k++; - /*Finally the left truncation point is found*/ - pFG->left = (int)floor(m - k*sqrt(lambda)- 1.5 ); - /* for small lambda the above calculation can yield negative truncation points, crop them here */ - if(pFG->left < 0) - pFG->left = 0; - /* perform underflow check */ - k_prime = k + 3.0 / (2.0 * sqrt_lambda); - /*We take the c_m_inf = 0.02935 / sqrt( m ), as for lambda >= 25 - c_m = 1 / ( sqrt( 2.0 * pi * m ) ) * exp( m - lambda - 1 / ( 12.0 * m ) ) => c_m_inf*/ - c_m_inf = 0.02935 / sqrt((double) m); - result = 0.0; - if( 0.0 < k_prime && k_prime <= sqrt_lambda / 2.0 ) - { - result = c_m_inf * exp( - pow(k_prime,2.0) / 2.0 - pow(k_prime, 3.0) / (3.0 * sqrt_lambda) ); - } - else - { - if( k_prime <= sqrt( m + 1.0 ) / m ) - { - double result_1 = c_m_inf * pow( - 1.0 - k_prime / sqrt((double) (m + 1)), - k_prime * sqrt((double) (m + 1))); - double result_2 = exp( - lambda ); - /*Take the maximum*/ - result = ( result_1 > result_2 ? result_1 : result_2); - } - else - { - /*NOTE: It will be an underflow error*/; - printf("ERROR: Fox-Glynn: lambda >= 25, underflow. The results are UNRELIABLE.\n"); - } - } - if ( result * omega / ( 1.0e+10 * ( pFG->right - pFG->left ) ) <= tau ) - { - printf("ERROR: Fox-Glynn: lambda >= 25, underflow. The results are UNRELIABLE.\n"); - } - } - - - - /*We still have to perform an underflow check for the right truncation point when lambda >= 400*/ - if( lambda >= lambda_400 ) - { - k_prime = k_rtp * sqrt_2 + 3.0 / (2.0 * sqrt_lambda); - /*We take the c_m_inf = 0.02935 / sqrt( m ), as for lambda >= 25 - c_m = 1 / ( sqrt( 2.0 * pi * m ) ) * exp( m - lambda - 1 / ( 12.0 * m ) ) => c_m_inf*/ - c_m_inf = 0.02935 / sqrt((double) m); - result = c_m_inf * exp( - pow( k_prime + 1.0 , 2.0 ) / 2.0 ); - if( result * omega / ( 1.0e+10 * ( pFG->right - pFG->left ) ) <= tau) - { - printf("ERROR: Fox-Glynn: lambda >= 400, underflow. The results are UNRELIABLE.\n"); - } - } - /*Time to set the initial value for weights*/ - *pw_m = omega / ( 1.0e+10 * ( pFG->right - pFG->left ) ); - - return true; - } - -/***************************************************************************** -Name : weighter -Role : The WEIGHTER function from the Fox-Glynn algorithm -@param : double lambda: (rate of uniformization)*(mission time) -@param : double tau: underflow -@param : double omega: overflow -@param : double epsilon: error bound -@param : FoxGlynn *: return by reference -@return : TRUE if everything is fine, otherwise FALSE. - This is the F parameter of Fox-Glynn finder function. -remark : -******************************************************************************/ - static bool weighter(const double lambda, const double tau, const double omega, const double epsilon, FoxGlynn *pFG) - { - static const double pi = 3.14159265358979323846264; - static const double lambda_25 = 25.0; - static const double lambda_400 = 40; - /*The magic m point*/ - const int m = (int)floor(lambda); - double w_m = 0; - int j, s, t; - - if( ! finder( m, lambda, tau, omega, epsilon, &w_m, pFG ) ) - return false; - - /*Allocate space for weights*/ - pFG->weights = (double *) calloc((size_t) (pFG->right - pFG->left + 1), - sizeof(double)); - /*Set an initial weight*/ - pFG->weights[ m - pFG->left ] = w_m; - - /*Fill the left side of the array*/ - for( j = m; j > pFG->left; j-- ) - pFG->weights[ ( j - pFG->left ) - 1 ] = ( j / lambda ) * pFG->weights[ j - pFG->left ]; - - /*Fill the right side of the array, have two cases lambda < 400 & lambda >= 400*/ - if( lambda < lambda_400 ) - { - /*Perform the underflow check, according to Fox-Glynn*/ - if( pFG->right > 600 ) - { - printf("ERROR: Fox-Glynn: pFG->right > 600, underflow is possible\n"); - return false; - } - /*Compute weights*/ - for( j = m; j < pFG->right; j++ ) - { - double q = lambda / ( j + 1 ); - if( pFG->weights[ j - pFG->left ] > tau / q ) - { - pFG->weights[ ( j - pFG->left ) + 1 ] = q * pFG->weights[ j - pFG->left ]; - }else{ - pFG->right = j; - break; /*It's time to compute W*/ - } - } - }else{ - /*Compute weights*/ - for( j = m; j < pFG->right; j++ ) - pFG->weights[ ( j - pFG->left ) + 1 ] = ( lambda / ( j + 1 ) ) * pFG->weights[ j - pFG->left ]; - } - - /*It is time to compute the normalization weight W*/ - pFG->total_weight = 0.0; - s = pFG->left; - t = pFG->right; - - while( s < t ) - { - if( pFG->weights[ s - pFG->left ] <= pFG->weights[ t - pFG->left ] ) - { - pFG->total_weight += pFG->weights[ s - pFG->left ]; - s++; - }else{ - pFG->total_weight += pFG->weights[ t - pFG->left ]; - t--; - } - } - pFG->total_weight += pFG->weights[ s - pFG->left ]; - - /* printf("Fox-Glynn: ltp = %d, rtp = %d, w = %10.15le \n", pFG->left, pFG->right, pFG->total_weight); */ - - return true; - } - -/***************************************************************************** -Name : fox_glynn -Role : get poisson probabilities. -@param : double lambda: (rate of uniformization)*(mission time) -@param : double tau: underflow -@param : double omega: overflow -@param : double epsilon: error bound -@param : FoxGlynn **: return a new FoxGlynn structure by reference -@return : TRUE if it worked fine, otherwise false -remark : -******************************************************************************/ - static bool fox_glynn(const double lambda, const double tau, const double omega, const double epsilon, FoxGlynn **ppFG) - { - /* printf("Fox-Glynn: lambda = %3.3le, epsilon = %1.8le\n",lambda, epsilon); */ - - *ppFG = (FoxGlynn *) calloc((size_t) 1, sizeof(FoxGlynn)); - (*ppFG)->weights = NULL; - - return weighter(lambda, tau, omega, epsilon, *ppFG); - } - - -/** -* Fries the memory allocated for the FoxGlynn structure -* @param fg the structure to free -*/ - static void freeFG(FoxGlynn * fg) - { - if( fg ){ - if( fg->weights ) - free(fg->weights); - free(fg); - } - } - - /*! - * Computes the poission-distribution - * - * - * @param parameter lambda to use - * @param point i - * TODO: replace with Fox-Lynn - * @return the probability - */ - - - /*! - * Computes the poission-distribution - * - * - * @param parameter lambda to use - * @param point i - * TODO: replace with Fox-glynn - * @return the probability - */ - template - static ValueType poisson(ValueType lambda, uint64_t i); template ::SupportsExponential, int>::type=0> static uint64_t trajans(storm::storage::SparseMatrix const& TransitionMatrix, uint64_t node, std::vector& disc, std::vector& finish, uint64_t * counter);