From 7db58c6374662f61ffb829494cde5f0023dad097 Mon Sep 17 00:00:00 2001 From: Timo Philipp Gros Date: Tue, 28 Nov 2017 17:22:30 +0100 Subject: [PATCH] using existing fox glynn now --- .../helper/SparseMarkovAutomatonCslHelper.cpp | 340 ++++++++++-------- .../helper/SparseMarkovAutomatonCslHelper.h | 3 + 2 files changed, 196 insertions(+), 147 deletions(-) diff --git a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp index d83286338..4264a9b28 100644 --- a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp +++ b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp @@ -369,6 +369,45 @@ namespace storm { return finish[node]; } + template::SupportsExponential, int>::type= 0> + void SparseMarkovAutomatonCslHelper::identify( + storm::storage::SparseMatrix const &fullTransitionMatrix, + storm::storage::BitVector const &markovianStates, storm::storage::BitVector const& psiStates) { + auto indices = fullTransitionMatrix.getRowGroupIndices(); + bool realProb = false; + bool NDM = false; + bool Alternating = true; + bool probStates = false; + bool markStates = false; + + for (uint64_t i=0; i::SupportsExponential, int>::type=0> storm::storage::BitVector SparseMarkovAutomatonCslHelper::identifyProbCyclesGoalStates(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& cycleStates) { @@ -463,174 +502,181 @@ namespace storm { 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 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(); + 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; - //(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; + //vectors to save calculation + std::vector> vd,vu,wu; + std::vector>> unifVectors{}; - //calculate relative ReachabilityVectors - std::vector in(numberOfStates, 0); - std::vector> relReachability(fullTransitionMatrix.getRowCount(),in); + //transitions from goalStates will be ignored. still: they are not allowed to be probabilistic! + for (uint64_t i =0 ; i 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(); - //calculate relative reachability - for(uint64_t i=0 ; i() /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 from = rowGroupIndices[i]; - auto to = rowGroupIndices[i+1]; - for (auto j=from ; j& act = relReachability[j]; - for(auto element: fullTransitionMatrix.getRow(j)){ - if (markovianStates[element.getColumn()]){ - act[element.getColumn()]=element.getValue(); - } - } - } - } + uint64_t N; - //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 relative ReachabilityVectors + std::vector in(numberOfStates, 0); + std::vector> relReachability(fullTransitionMatrix.getRowCount(),in); - // 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); - 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, 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); - } - } - printTransitions(N, maxNorm, fullTransitionMatrix,exitRate,markovianStates,psiStates,relReachability,psiStates, psiStates,unifVectors, logfile); //TODO remove - // (6) double lambda - lambda=2*lambda; + //calculate relative reachability - // (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)); + for(uint64_t i=0 ; i& act = relReachability[j]; + for(auto element: fullTransitionMatrix.getRow(j)){ + if (markovianStates[element.getColumn()]){ + act[element.getColumn()]=element.getValue(); + } + } + } + } + + //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]; - logfile.close(); - return unifVectors[0][0]; } template ::SupportsExponential, int>::type> diff --git a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h index ff536e71f..ab036c8c3 100644 --- a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h +++ b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h @@ -56,6 +56,9 @@ namespace storm { static std::vector computeReachabilityTimes(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); private: + template ::SupportsExponential, int>::type=0> + static void identify(storm::storage::SparseMatrix const& fullTransitionMatrix, storm::storage::BitVector const& markovianStates,storm::storage::BitVector const& psiStates); + 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, std::vector const& poisson);