diff --git a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp index d88a50cec..a9d8c0158 100644 --- a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp +++ b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp @@ -211,7 +211,7 @@ namespace storm { template ::SupportsExponential, int>::type> - void SparseMarkovAutomatonCslHelper::calculateVu(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){ + void SparseMarkovAutomatonCslHelper::calculateVu(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 poisson){ if (unifVectors[1][k][node]!=-1){return;} //dynamic programming. avoiding multiple calculation. uint64_t N = unifVectors[1].size()-1; auto rowGroupIndices = fullTransitionMatrix.getRowGroupIndices(); @@ -219,10 +219,10 @@ namespace storm { ValueType res =0; for (uint64_t i = k ; i < N ; i++ ){ if (unifVectors[2][N-1-(i-k)][node]==-1){ - calculateUnifPlusVector(N-1-(i-k),node,2,lambda,probSize,relativeReachability,dir,unifVectors,fullTransitionMatrix, markovianStates,psiStates,solver, logfile); + calculateUnifPlusVector(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(lambda, i)*unifVectors[2][N-1-(i-k)][node]; + res+=poisson[i]*unifVectors[2][N-1-(i-k)][node]; } unifVectors[1][k][node]=res; } @@ -238,7 +238,7 @@ namespace storm { storm::storage::SparseMatrix const &fullTransitionMatrix, storm::storage::BitVector const &markovianStates, storm::storage::BitVector const &psiStates, - std::unique_ptr> const &solver, std::ofstream& logfile) { + std::unique_ptr> const &solver, std::ofstream& logfile, std::vector poisson) { if (unifVectors[kind][k][node]!=-1){ @@ -267,7 +267,7 @@ namespace storm { // Vd res = storm::utility::zero(); for (uint64_t i = k ; i SparseMarkovAutomatonCslHelper::unifPlus(OptimizationDirection dir, std::pair const& boundsPair, std::vector const& exitRateVector, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& markovStates, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory){ STORM_LOG_TRACE("Using UnifPlus to compute bounded until probabilities."); - - std::ofstream logfile("U+logfile.txt", std::ios::app); ValueType maxNorm = storm::utility::zero(); ValueType oldDiff = -storm::utility::zero(); @@ -589,6 +587,9 @@ namespace storm { exitRate[i] = exitNew; } + // calculate poisson distribution + std::vector poisson = foxGlynnProb(lambda*T, N, epsilon*kappa); + // (4) define vectors/matrices std::vector init(numberOfStates, -1); vd = std::vector> (N + 1, init); @@ -604,9 +605,9 @@ namespace storm { // (5) calculate vectors and maxNorm for (uint64_t i = 0; i < numberOfStates; i++) { for (uint64_t k = N; k <= N; k--) { - calculateUnifPlusVector(k,i,0,lambda,probSize,relReachability,dir,unifVectors,fullTransitionMatrix,markovianStates,psiStates,solver, logfile); - calculateUnifPlusVector(k,i,2,lambda,probSize,relReachability,dir,unifVectors,fullTransitionMatrix,markovianStates,psiStates,solver, logfile); - calculateVu(relReachability,dir,k,i,1,lambda,probSize,unifVectors,fullTransitionMatrix,markovianStates,psiStates,solver, logfile); + calculateUnifPlusVector(k,i,0,lambda,probSize,relReachability,dir,unifVectors,fullTransitionMatrix,markovianStates,psiStates,solver, logfile, poisson); + calculateUnifPlusVector(k,i,2,lambda,probSize,relReachability,dir,unifVectors,fullTransitionMatrix,markovianStates,psiStates,solver, logfile, poisson); + calculateVu(relReachability,dir,k,i,1,lambda,probSize,unifVectors,fullTransitionMatrix,markovianStates,psiStates,solver, logfile, poisson); //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); @@ -615,6 +616,7 @@ namespace storm { printTransitions(N, maxNorm, fullTransitionMatrix,exitRate,markovianStates,psiStates,relReachability,psiStates, psiStates,unifVectors, logfile); //TODO remove // (6) double lambda + lambda=2*lambda; // (7) escape if not coming closer to solution diff --git a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h index b2ad305a1..5e8c16c33 100644 --- a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h +++ b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h @@ -57,7 +57,7 @@ namespace storm { private: template ::SupportsExponential, int>::type=0> - static void calculateUnifPlusVector(uint64_t k, uint64_t node, uint64_t const kind, ValueType lambda, uint64_t probSize, std::vector> const& relativeReachability, OptimizationDirection dir, 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); + static void calculateUnifPlusVector(uint64_t k, uint64_t node, uint64_t const kind, ValueType lambda, uint64_t probSize, std::vector> const& relativeReachability, OptimizationDirection dir, 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 poisson); template ::SupportsExponential, int>::type=0> static void deleteProbDiagonals(storm::storage::SparseMatrix& transitionMatrix, storm::storage::BitVector const& markovianStates); @@ -80,6 +80,321 @@ namespace storm { template ::SupportsExponential, int>::type=0> static storm::storage::BitVector identifyProbCycles(storm::storage::SparseMatrix const& TransitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates); + + //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; + + //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 * @@ -100,7 +415,7 @@ namespace storm { * */ template ::SupportsExponential, int>::type=0> - static void calculateVu(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); + static void calculateVu(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 poisson);