diff --git a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp index 4264a9b28..b4c7b0526 100644 --- a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp +++ b/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,8 +267,10 @@ namespace storm { // Vd res = storm::utility::zero(); for (uint64_t i = k ; i::SupportsExponential, int>::type> - std::vector 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."); + template::SupportsExponential, int>::type> + std::vector 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(); + + //bitvectors to identify different kind of states + storm::storage::BitVector markovianStates = markovStates; + storm::storage::BitVector allStates(markovianStates.size(), true); + storm::storage::BitVector probabilisticStates = ~markovianStates; + + + //vectors to save calculation + std::vector> vd, vu, wu; + std::vector>> unifVectors{}; + + + //transitions from goalStates will be ignored. still: they are not allowed to be probabilistic! + for (uint64_t i = 0; i < psiStates.size(); i++) { + if (psiStates[i]) { + markovianStates.set(i, true); + probabilisticStates.set(i, false); + } + } + + //transition matrix with diagonal entries. The values can be changed during uniformisation + std::vector exitRate{exitRateVector}; + typename storm::storage::SparseMatrix fullTransitionMatrix = transitionMatrix.getSubmatrix( + true, allStates, allStates, true); + // delete diagonals + deleteProbDiagonals(fullTransitionMatrix, markovianStates); + typename storm::storage::SparseMatrix 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); + probSize = probMatrix.getRowCount(); + } + + auto &rowGroupIndices = fullTransitionMatrix.getRowGroupIndices(); + + + //(1) define horizon, epsilon, kappa , N, lambda, + uint64_t numberOfStates = fullTransitionMatrix.getRowGroupCount(); + double T = boundsPair.second; + ValueType kappa = storm::utility::one() / 10; // would be better as option-parameter + ValueType epsilon = storm::settings::getModule().getPrecision(); + ValueType lambda = exitRateVector[0]; + for (ValueType act: exitRateVector) { + lambda = std::max(act, lambda); + } + uint64_t N; + - std::ofstream logfile("U+logfile.txt", std::ios::app); - ValueType maxNorm = storm::utility::zero(); - ValueType oldDiff = -storm::utility::zero(); - //bitvectors to identify different kind of states - storm::storage::BitVector markovianStates = markovStates; - storm::storage::BitVector allStates(markovianStates.size(), true); - storm::storage::BitVector probabilisticStates = ~markovianStates; + //calculate relative ReachabilityVectors + std::vector in(numberOfStates, 0); + std::vector> relReachability(fullTransitionMatrix.getRowCount(), in); - //vectors to save calculation - std::vector> vd,vu,wu; - std::vector>> unifVectors{}; + //calculate relative reachability - //transitions from goalStates will be ignored. still: they are not allowed to be probabilistic! - for (uint64_t i =0 ; i &act = relReachability[j]; + for (auto element: fullTransitionMatrix.getRow(j)) { + if (markovianStates[element.getColumn()]) { + act[element.getColumn()] = element.getValue(); } } - - //transition matrix with diagonal entries. The values can be changed during uniformisation - std::vector exitRate{exitRateVector}; - typename storm::storage::SparseMatrix fullTransitionMatrix = transitionMatrix.getSubmatrix(true, allStates , allStates , true); - // delete diagonals - deleteProbDiagonals(fullTransitionMatrix, markovianStates); - typename storm::storage::SparseMatrix 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); - probSize = probMatrix.getRowCount(); } + } - auto& rowGroupIndices = fullTransitionMatrix.getRowGroupIndices(); + //create equitation solver + storm::solver::MinMaxLinearEquationSolverRequirements requirements = minMaxLinearEquationSolverFactory.getRequirements( + true, dir); + requirements.clearBounds(); + STORM_LOG_THROW(requirements.empty(), storm::exceptions::UncheckedRequirementException, + "Cannot establish requirements for solver."); + + std::unique_ptr> solver; + if (probSize != 0) { + solver = minMaxLinearEquationSolverFactory.create(probMatrix); + solver->setHasUniqueSolution(); + solver->setBounds(storm::utility::zero(), storm::utility::one()); + solver->setRequirementsChecked(); + solver->setCachingEnabled(true); + } + // while not close enough to precision: + do { + maxNorm = storm::utility::zero(); + // (2) update parameter + N = ceil(lambda * T * exp(2) - log(kappa * epsilon)); + + // (3) uniform - just applied to markovian states + for (uint_fast64_t i = 0; i < fullTransitionMatrix.getRowGroupCount(); i++) { + if (!markovianStates[i]) { + continue; + } + uint64_t from = rowGroupIndices[i]; //markovian state -> no Nondeterminism -> only one row + if (exitRate[i] == lambda) { + continue; //already unified + } - //(1) define horizon, epsilon, kappa , N, lambda, - uint64_t numberOfStates = fullTransitionMatrix.getRowGroupCount(); - double T = boundsPair.second; - ValueType kappa = storm::utility::one() /10; // would be better as option-parameter - ValueType epsilon = storm::settings::getModule().getPrecision(); - ValueType lambda = exitRateVector[0]; - for (ValueType act: exitRateVector) { - lambda = std::max(act, lambda); + auto line = fullTransitionMatrix.getRow(from); + ValueType exitOld = exitRate[i]; + ValueType exitNew = lambda; + for (auto &v : line) { + if (v.getColumn() == i) { //diagonal element + ValueType newSelfLoop = exitNew - exitOld + v.getValue(); + ValueType newRate = newSelfLoop / exitNew; + v.setValue(newRate); + } else { //modify probability + ValueType propOld = v.getValue(); + ValueType propNew = propOld * exitOld / exitNew; + v.setValue(propNew); + } + } + exitRate[i] = exitNew; } - uint64_t N; + // calculate poisson distribution + std::tuple> foxGlynnResult = storm::utility::numerical::getFoxGlynnCutoff( + T * lambda, 1e+300, epsilon * epsilon * kappa); - //calculate relative ReachabilityVectors - std::vector in(numberOfStates, 0); - std::vector> relReachability(fullTransitionMatrix.getRowCount(),in); + // 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 init(numberOfStates, -1); + vd = std::vector>(N + 1, init); + vu = std::vector>(N + 1, init); + wu = std::vector>(N + 1, init); + + unifVectors.clear(); + unifVectors.push_back(vd); + unifVectors.push_back(vd); + unifVectors.push_back(vd); + + //define 0=vd 1=vu 2=wu + // (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)); + //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 - //calculate relative reachability + // (6) double lambda - for(uint64_t i=0 ; i& act = relReachability[j]; - for(auto element: fullTransitionMatrix.getRow(j)){ - if (markovianStates[element.getColumn()]){ - act[element.getColumn()]=element.getValue(); - } - } - } - } + lambda = 2 * lambda; - //create equitation solver - storm::solver::MinMaxLinearEquationSolverRequirements requirements = minMaxLinearEquationSolverFactory.getRequirements(true, dir); - requirements.clearBounds(); - STORM_LOG_THROW(requirements.empty(), storm::exceptions::UncheckedRequirementException, "Cannot establish requirements for solver."); - - std::unique_ptr> solver; - if (probSize!=0){ - solver = minMaxLinearEquationSolverFactory.create(probMatrix); - solver->setHasUniqueSolution(); - solver->setBounds(storm::utility::zero(), storm::utility::one()); - solver->setRequirementsChecked(); - solver->setCachingEnabled(true); - } - // while not close enough to precision: - do { - maxNorm = storm::utility::zero(); - // (2) update parameter - N = ceil(lambda*T*exp(2)-log(kappa*epsilon)); - - // (3) uniform - just applied to markovian states - for (uint_fast64_t i = 0; i < fullTransitionMatrix.getRowGroupCount(); i++) { - if (!markovianStates[i]) { - continue; - } - uint64_t from = rowGroupIndices[i]; //markovian state -> no Nondeterminism -> only one row - - if (exitRate[i] == lambda) { - continue; //already unified - } - - auto line = fullTransitionMatrix.getRow(from); - ValueType exitOld = exitRate[i]; - ValueType exitNew = lambda; - for (auto &v : line) { - if (v.getColumn() == i) { //diagonal element - ValueType newSelfLoop = exitNew - exitOld + v.getValue(); - ValueType newRate = newSelfLoop / exitNew; - v.setValue(newRate); - } else { //modify probability - ValueType propOld = v.getValue(); - ValueType propNew = propOld * exitOld / exitNew; - v.setValue(propNew); - } - } - exitRate[i] = exitNew; - } - - // calculate poisson distribution - std::tuple> foxGlynnResult = storm::utility::numerical::getFoxGlynnCutoff(T*lambda, 1e+300, epsilon*kappa/ 8.0); - - // 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 init(numberOfStates, -1); - vd = std::vector> (N + 1, init); - vu = std::vector> (N + 1, init); - wu = std::vector> (N + 1, init); - - unifVectors.clear(); - unifVectors.push_back(vd); - unifVectors.push_back(vd); - unifVectors.push_back(vd); - - //define 0=vd 1=vu 2=wu - // (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)); - //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 - - // (6) double lambda - - lambda=2*lambda; - - // (7) escape if not coming closer to solution - if (oldDiff!=-1){ - if (oldDiff==maxNorm){ - std::cout << "Not coming closer to solution as " << maxNorm << "/n"; - break; - } - } - oldDiff = maxNorm; - } while (maxNorm>epsilon*(1-kappa)); - - logfile.close(); - return unifVectors[0][0]; + // (7) escape if not coming closer to solution + if (oldDiff != -1) { + if (oldDiff == maxNorm) { + std::cout << "Not coming closer to solution as " << maxNorm << "/n"; + break; + } + } + oldDiff = maxNorm; + } while (maxNorm > epsilon * (1 - kappa)); + + logfile.close(); + return unifVectors[0][0]; }