From 0d1de8aba9e929e0b6dc71270066e9cff72db018 Mon Sep 17 00:00:00 2001 From: Timo Philipp Gros Date: Sat, 6 Jan 2018 13:23:20 +0100 Subject: [PATCH] restructured code, SCC missing --- .../helper/SparseMarkovAutomatonCslHelper.cpp | 303 +++++------------- .../helper/SparseMarkovAutomatonCslHelper.h | 28 +- 2 files changed, 76 insertions(+), 255 deletions(-) diff --git a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp index 4b419d92e..ad95dee20 100644 --- a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp +++ b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp @@ -158,70 +158,13 @@ namespace storm { STORM_LOG_THROW(false, storm::exceptions::InvalidOperationException, "Computing bounded reachability probabilities is unsupported for this value type."); } - template::SupportsExponential, int>::type> - void SparseMarkovAutomatonCslHelper::printTransitions(const uint64_t N, ValueType const diff, - storm::storage::SparseMatrix const &fullTransitionMatrix, - std::vector const &exitRateVector, storm::storage::BitVector const &markovianStates, - storm::storage::BitVector const &psiStates, std::vector> relReachability, - const storage::BitVector &cycleStates, const storage::BitVector &cycleGoalStates, - std::vector>> &unifVectors, std::ofstream& logfile) { - - auto const& rowGroupIndices = fullTransitionMatrix.getRowGroupIndices(); - auto numberOfStates = fullTransitionMatrix.getRowGroupCount(); - - //Transition Matrix - logfile << "number of states = num of row group count " << numberOfStates << "\n"; - for (uint_fast64_t i = 0; i < fullTransitionMatrix.getRowGroupCount(); i++) { - logfile << " from node " << i << " "; - auto from = rowGroupIndices[i]; - auto to = rowGroupIndices[i+1]; - for (auto j = from ; j < to; j++){ - for (auto &v : fullTransitionMatrix.getRow(j)) { - if (markovianStates[i]){ - logfile << v.getValue() *exitRateVector[i] << " -> "<< v.getColumn() << "\t"; - } else { - logfile << v.getValue() << " -> "<< v.getColumn() << "\t"; - } - } - logfile << "\n"; - } - } - logfile << "\n"; - - logfile << "probStates\tmarkovianStates\tgoalStates\tcycleStates\tcycleGoalStates\n"; - for (int i =0 ; i< markovianStates.size() ; i++){ - logfile << (~markovianStates)[i] << "\t\t" << markovianStates[i] << "\t\t" << psiStates[i] << "\t\t" << cycleStates[i] << "\t\t" << cycleGoalStates[i] << "\n"; - } - - logfile << "Iteration for N = " << N << " maximal difference was " << diff << "\n"; - - logfile << "vd: \n"; - for (uint64_t i =0 ; i::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, storm::utility::numerical::FoxGlynnResult const & poisson){ + 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, + storm::utility::numerical::FoxGlynnResult const & poisson, bool cycleFree){ if (unifVectors[1][k][node]!=-1){return;} //dynamic programming. avoiding multiple calculation. uint64_t N = unifVectors[1].size()-1; auto const& rowGroupIndices = fullTransitionMatrix.getRowGroupIndices(); @@ -229,7 +172,7 @@ namespace storm { ValueType res =0; for (uint64_t i = k ; i < N ; i++ ){ if (unifVectors[2][N-1-(i-k)][node]==-1){ - calculateUnifPlusVector(env, N-1-(i-k),node,2,lambda,probSize,relativeReachability,dir,unifVectors,fullTransitionMatrix, markovianStates,psiStates,solver, logfile, poisson); + calculateUnifPlusVector(env, N-1-(i-k),node,2,lambda,probSize,relativeReachability,dir,unifVectors,fullTransitionMatrix, markovianStates,psiStates,solver, logfile, poisson, cycleFree); //old: relativeReachability, dir, (N-1-(i-k)),node,lambda,wu,fullTransitionMatrix,markovianStates,psiStates, solver); } if (i>=poisson.left && i<=poisson.right){ @@ -250,7 +193,8 @@ 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, storm::utility::numerical::FoxGlynnResult const & poisson) { + std::unique_ptr> const &solver, std::ofstream& logfile, + storm::utility::numerical::FoxGlynnResult const & poisson, bool cycleFree) { if (unifVectors[kind][k][node]!=-1){ @@ -300,7 +244,7 @@ namespace storm { for (auto &element : line){ uint64_t to = element.getColumn(); if (unifVectors[kind][k+1][to]==-1){ - calculateUnifPlusVector(env, k+1,to,kind,lambda,probSize,relativeReachability,dir,unifVectors,fullTransitionMatrix,markovianStates,psiStates,solver, logfile, poisson); + calculateUnifPlusVector(env, k+1,to,kind,lambda,probSize,relativeReachability,dir,unifVectors,fullTransitionMatrix,markovianStates,psiStates,solver, logfile, poisson, cycleFree); } res+=element.getValue()*unifVectors[kind][k+1][to]; } @@ -310,164 +254,80 @@ namespace storm { } //probabilistic non-goal State - res = -1; - uint64_t rowStart = rowGroupIndices[node]; - uint64_t rowEnd = rowGroupIndices[node+1]; - for (uint64_t i = rowStart; i< rowEnd; i++){ - auto line = fullTransitionMatrix.getRow(i); - ValueType between = 0; - for (auto& element: line){ - uint64_t to = element.getColumn(); - if (to==node){ - continue; - } - if (unifVectors[kind][k][to]==-1){ - calculateUnifPlusVector(env, k,to,kind,lambda,probSize,relativeReachability,dir,unifVectors,fullTransitionMatrix,markovianStates,psiStates,solver, logfile, poisson); - } - between+=element.getValue()*unifVectors[kind][k][to]; - } - if (maximize(dir)){ - res = std::max(res,between); - } else { - if (res!=-1){ - res = std::min(res,between); - } else { - res = between; - } - } - } - unifVectors[kind][k][node]=res; - //logfile << print << "probabilistic state: "<< " res = " << unifVectors[kind][k][node] << " but calculated more \n"; - - //end probabilistic states - } - - - template ::SupportsExponential, int>::type> - uint64_t SparseMarkovAutomatonCslHelper::trajans(storm::storage::SparseMatrix const& transitionMatrix, uint64_t node, std::vector& disc, std::vector& finish, uint64_t* counter) { - auto const& rowGroupIndice = transitionMatrix.getRowGroupIndices(); - - disc[node] = *counter; - finish[node] = *counter; - (*counter)+=1; - - auto from = rowGroupIndice[node]; - auto to = rowGroupIndice[node+1]; - - for(uint64_t i =from; i::SupportsExponential, int>::type> - 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> - storm::storage::BitVector SparseMarkovAutomatonCslHelper::identifyProbCyclesGoalStates(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& cycleStates) { - - storm::storage::BitVector goalStates(cycleStates.size(), false); - auto const& rowGroupIndices = transitionMatrix.getRowGroupIndices(); - - for (uint64_t i = 0 ; i < transitionMatrix.getRowGroupCount() ; i++){ - if (!cycleStates[i]){ + //not cycle free - use solver technique, calling SVI per default + //solving all sub-MDP's in one iteration + std::vector b(probSize, 0), x(numberOfProbStates,0); + //calculate b + uint64_t lineCounter=0; + for (int i =0; i::SupportsExponential, int>::type> - storm::storage::BitVector SparseMarkovAutomatonCslHelper::identifyProbCycles(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates){ - - storm::storage::BitVector const& probabilisticStates = ~markovianStates; - storm::storage::BitVector const& probabilisticNonGoalStates = ~markovianStates & ~psiStates; - - storm::storage::SparseMatrix const& probMatrix = transitionMatrix.getSubmatrix(true, probabilisticNonGoalStates, probabilisticNonGoalStates); - uint64_t probSize = probMatrix.getRowGroupCount(); - std::vector disc(probSize, 0), finish(probSize, 0); - - uint64_t counter =1; - - for (uint64_t i =0; isolveEquations(env, dir, x, b); - storm::storage::BitVector cycleStates(markovianStates.size(), false); - for (int i = 0 ; i< finish.size() ; i++){ - auto f = finish[i]; - for (int j =i+1; j::SupportsExponential, int>::type> void SparseMarkovAutomatonCslHelper::deleteProbDiagonals(storm::storage::SparseMatrix& transitionMatrix, storm::storage::BitVector const& markovianStates){ auto const& rowGroupIndices = transitionMatrix.getRowGroupIndices(); @@ -512,6 +372,7 @@ namespace storm { storm::solver::MinMaxLinearEquationSolverFactory const &minMaxLinearEquationSolverFactory) { STORM_LOG_TRACE("Using UnifPlus to compute bounded until probabilities."); + bool cycleFree; std::ofstream logfile("U+logfile.txt", std::ios::app); //logfile << "Using U+\n"; @@ -571,7 +432,6 @@ namespace storm { //calculate relative reachability - /* for (uint64_t i = 0; i < numberOfStates; i++) { if (markovianStates[i]) { continue; @@ -593,15 +453,14 @@ namespace storm { requirements.clearBounds(); STORM_LOG_THROW(requirements.empty(), storm::exceptions::UncheckedRequirementException, "Cannot establish requirements for solver."); - */ std::unique_ptr> solver; - /*if (probSize != 0) { + if (probSize != 0) { solver = minMaxLinearEquationSolverFactory.create(env, probMatrix); solver->setHasUniqueSolution(); solver->setBounds(storm::utility::zero(), storm::utility::one()); solver->setRequirementsChecked(); solver->setCachingEnabled(true); - } */ + } // while not close enough to precision: do { //logfile << "starting iteration\n"; @@ -661,15 +520,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, - foxGlynnResult); + foxGlynnResult, cycleFree); calculateUnifPlusVector(env, k, i, 2, lambda, probSize, relReachability, dir, unifVectors, fullTransitionMatrix, markovianStates, psiStates, solver, logfile, - foxGlynnResult); + foxGlynnResult, cycleFree); calculateVu(env, relReachability, dir, k, i, 1, lambda, probSize, unifVectors, fullTransitionMatrix, markovianStates, psiStates, solver, logfile, - 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]); + foxGlynnResult, cycleFree); } } @@ -678,23 +535,11 @@ namespace storm { ValueType diff = std::abs(unifVectors[0][0][i]-unifVectors[1][0][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; - std::cout << "Finished Iteration for N = " << N << " with difference " << maxNorm << "\n"; - } while (maxNorm > epsilon /* * (1 - kappa)*/); + } while (maxNorm > epsilon*(1 - kappa)); logfile.close(); return unifVectors[0][0]; diff --git a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h index ac005ce38..82492aa6e 100644 --- a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h +++ b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h @@ -60,10 +60,7 @@ namespace storm { 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(Environment const& env, 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, storm::utility::numerical::FoxGlynnResult const & poisson); + static void calculateUnifPlusVector(Environment const& env, 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, storm::utility::numerical::FoxGlynnResult const & poisson, bool cycleFree); template ::SupportsExponential, int>::type=0> static void deleteProbDiagonals(storm::storage::SparseMatrix& transitionMatrix, storm::storage::BitVector const& markovianStates); @@ -80,28 +77,12 @@ namespace storm { return id-1; } - template ::SupportsExponential, int>::type=0> - static storm::storage::BitVector identifyProbCyclesGoalStates(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& cycleStates); - - - 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 - - - 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); - /* * Computes vu vector according to UnifPlus * */ template ::SupportsExponential, int>::type=0> - static void 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, storm::utility::numerical::FoxGlynnResult const & poisson); - - - + static void 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, storm::utility::numerical::FoxGlynnResult const & poisson, bool cycleFree); /*! * Prints the TransitionMatrix and the vectors vd, vu, wu to the logfile @@ -109,11 +90,6 @@ namespace storm { * */ - template ::SupportsExponential, int>::type=0> - static void printTransitions(const uint64_t N, ValueType const diff, storm::storage::SparseMatrix const& fullTransitionMatrix, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, - storm::storage::BitVector const& psiStates, std::vector> relReachability, - storm::storage::BitVector const& cycleStates , storm::storage::BitVector const& cycleGoalStates ,std::vector>>& unifVectors, std::ofstream& logfile); - template ::SupportsExponential, int>::type = 0> static void computeBoundedReachabilityProbabilities(Environment const& env, OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRates, storm::storage::BitVector const& goalStates, storm::storage::BitVector const& markovianNonGoalStates, storm::storage::BitVector const& probabilisticNonGoalStates, std::vector& markovianNonGoalValues, std::vector& probabilisticNonGoalValues, ValueType delta, uint64_t numberOfSteps, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory);