fixed use of FoxLynn after CutOff

src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp

@@ -242,7 +242,7 @@ namespace storm {
 
 
                 if (unifVectors[kind][k][node]!=-1){
-                    logfile << "already calculated for k = " << k << " node = " << node << "\n";
+                    //logfile << "already calculated for k = " << k << " node = " << node << "\n";
                     return;
                 }
                 std::string print = std::string("calculating vector ") + std::to_string(kind)   +  " for k = " + std::to_string(k) + " node " + std::to_string(node) +" \t";
@@ -267,9 +267,11 @@ namespace storm {
                         // Vd
                         res = storm::utility::zero<ValueType>();
                         for (uint64_t i = k ; i<N ; i++){
+                            if (i < poisson.size()){
                                 ValueType between = poisson[i];
                                 res+=between;
                             }
+                        }
                         unifVectors[kind][k][node]=res;
                     } else {
                         // WU
@@ -501,7 +503,13 @@ namespace storm {
             }
 
             template<typename ValueType, typename std::enable_if<storm::NumberTraits<ValueType>::SupportsExponential, int>::type>
-            std::vector<ValueType> SparseMarkovAutomatonCslHelper::unifPlus(OptimizationDirection dir, std::pair<double, double> const& boundsPair, std::vector<ValueType> const& exitRateVector, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& markovStates, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory){
+            std::vector<ValueType> SparseMarkovAutomatonCslHelper::unifPlus(OptimizationDirection dir,
+                                                                            std::pair<double, double> const &boundsPair,
+                                                                            std::vector<ValueType> const &exitRateVector,
+                                                                            storm::storage::SparseMatrix<ValueType> const &transitionMatrix,
+                                                                            storm::storage::BitVector const &markovStates,
+                                                                            storm::storage::BitVector const &psiStates,
+                                                                            storm::solver::MinMaxLinearEquationSolverFactory<ValueType> const &minMaxLinearEquationSolverFactory) {
                 STORM_LOG_TRACE("Using UnifPlus to compute bounded until probabilities.");
 
 
@@ -530,13 +538,15 @@ namespace storm {
 
                 //transition matrix with diagonal entries. The values can be changed during uniformisation
                 std::vector<ValueType> exitRate{exitRateVector};
-                    typename storm::storage::SparseMatrix<ValueType> fullTransitionMatrix = transitionMatrix.getSubmatrix(true, allStates , allStates , true);
+                typename storm::storage::SparseMatrix<ValueType> fullTransitionMatrix = transitionMatrix.getSubmatrix(
+                        true, allStates, allStates, true);
                 // delete diagonals
                 deleteProbDiagonals(fullTransitionMatrix, markovianStates);
                 typename storm::storage::SparseMatrix<ValueType> probMatrix{};
                 uint64_t probSize = 0;
                 if (probabilisticStates.getNumberOfSetBits() != 0) { //work around in case there are no prob states
-                        probMatrix = fullTransitionMatrix.getSubmatrix(true, probabilisticStates , probabilisticStates, true);
+                    probMatrix = fullTransitionMatrix.getSubmatrix(true, probabilisticStates, probabilisticStates,
+                                                                   true);
                     probSize = probMatrix.getRowCount();
                 }
 
@@ -582,9 +592,11 @@ namespace storm {
                 }
 
                 //create equitation solver
-                                    storm::solver::MinMaxLinearEquationSolverRequirements requirements = minMaxLinearEquationSolverFactory.getRequirements(true, dir);
+                storm::solver::MinMaxLinearEquationSolverRequirements requirements = minMaxLinearEquationSolverFactory.getRequirements(
+                        true, dir);
                 requirements.clearBounds();
-                                    STORM_LOG_THROW(requirements.empty(), storm::exceptions::UncheckedRequirementException, "Cannot establish requirements for solver.");
+                STORM_LOG_THROW(requirements.empty(), storm::exceptions::UncheckedRequirementException,
+                                "Cannot establish requirements for solver.");
 
                 std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver;
                 if (probSize != 0) {
@@ -629,12 +641,17 @@ namespace storm {
                     }
 
                     // calculate poisson distribution
-                                        std::tuple<uint_fast64_t, uint_fast64_t, ValueType, std::vector<ValueType>> foxGlynnResult = storm::utility::numerical::getFoxGlynnCutoff(T*lambda, 1e+300, epsilon*kappa/ 8.0);
+
+                    std::tuple<uint_fast64_t, uint_fast64_t, ValueType, std::vector<ValueType>> foxGlynnResult = storm::utility::numerical::getFoxGlynnCutoff(
+                                T * lambda, 1e+300, epsilon * epsilon * kappa);
+
 
                     // Scale the weights so they add up to one.
                     for (auto &element : std::get<3>(foxGlynnResult)) {
                         element /= std::get<2>(foxGlynnResult);
                     }
+
+
                     // (4) define vectors/matrices
                     std::vector<ValueType> init(numberOfStates, -1);
                     vd = std::vector<std::vector<ValueType>>(N + 1, init);
@@ -650,15 +667,22 @@ 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, std::get<3>(foxGlynnResult));
-                                                    calculateUnifPlusVector(k,i,2,lambda,probSize,relReachability,dir,unifVectors,fullTransitionMatrix,markovianStates,psiStates,solver, logfile, std::get<3>(foxGlynnResult));
-                                                    calculateVu(relReachability,dir,k,i,1,lambda,probSize,unifVectors,fullTransitionMatrix,markovianStates,psiStates,solver, logfile, std::get<3>(foxGlynnResult));
+                            calculateUnifPlusVector(k, i, 0, lambda, probSize, relReachability, dir, unifVectors,
+                                                    fullTransitionMatrix, markovianStates, psiStates, solver, logfile,
+                                                    std::get<3>(foxGlynnResult));
+                            calculateUnifPlusVector(k, i, 2, lambda, probSize, relReachability, dir, unifVectors,
+                                                    fullTransitionMatrix, markovianStates, psiStates, solver, logfile,
+                                                    std::get<3>(foxGlynnResult));
+                            calculateVu(relReachability, dir, k, i, 1, lambda, probSize, unifVectors,
+                                        fullTransitionMatrix, markovianStates, psiStates, solver, logfile,
+                                        std::get<3>(foxGlynnResult));
                             //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);
                         }
                     }
-                                        printTransitions(N, maxNorm, fullTransitionMatrix,exitRate,markovianStates,psiStates,relReachability,psiStates, psiStates,unifVectors, logfile); //TODO remove
+                    printTransitions(N, maxNorm, fullTransitionMatrix, exitRate, markovianStates, psiStates,
+                                     relReachability, psiStates, psiStates, unifVectors, logfile); //TODO remove
 
                     // (6) double lambda