diff --git a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp index 67eec89ee..3ef29e743 100644 --- a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp +++ b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp @@ -213,7 +213,7 @@ namespace storm { } template ::SupportsExponential, int>::type> - void SparseMarkovAutomatonCslHelper::calculateVu(OptimizationDirection dir, uint64_t k, uint64_t node, ValueType lambda, std::vector>& vu, std::vector>& wu, storm::storage::SparseMatrix const& fullTransitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates){ + void SparseMarkovAutomatonCslHelper::calculateVu(std::vector> const& relativeReachability, OptimizationDirection dir, uint64_t k, uint64_t node, ValueType lambda, std::vector>& vu, std::vector>& wu, storm::storage::SparseMatrix const& fullTransitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory){ if (vu[k][node]!=-1){return;} //dynamic programming. avoiding multiple calculation. uint64_t N = vu.size()-1; auto rowGroupIndices = fullTransitionMatrix.getRowGroupIndices(); @@ -221,7 +221,7 @@ namespace storm { ValueType res =0; for (uint64_t i = k ; i < N ; i++ ){ if (wu[N-1-(i-k)][node]==-1){ - calculateWu(dir, (N-1-(i-k)),node,lambda,wu,fullTransitionMatrix,markovianStates,psiStates); + calculateWu(relativeReachability, dir, (N-1-(i-k)),node,lambda,wu,fullTransitionMatrix,markovianStates,psiStates, minMaxLinearEquationSolverFactory); } res+=poisson(lambda, i)*wu[N-1-(i-k)][node]; } @@ -229,7 +229,7 @@ namespace storm { } template ::SupportsExponential, int>::type> - void SparseMarkovAutomatonCslHelper::calculateWu(OptimizationDirection dir, uint64_t k, uint64_t node, ValueType lambda, std::vector>& wu, storm::storage::SparseMatrix const& fullTransitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates){ + void SparseMarkovAutomatonCslHelper::calculateWu(std::vector> const& relativeReachability, OptimizationDirection dir, uint64_t k, uint64_t node, ValueType lambda, std::vector>& wu, storm::storage::SparseMatrix const& fullTransitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory){ if (wu[k][node]!=-1){return;} //dynamic programming. avoiding multiple calculation. uint64_t N = wu.size()-1; auto const& rowGroupIndices = fullTransitionMatrix.getRowGroupIndices( ); @@ -251,11 +251,72 @@ namespace storm { for (auto &element : line){ uint64_t to = element.getColumn(); if (wu[k+1][to]==-1){ - calculateWu(dir, k+1,to,lambda,wu,fullTransitionMatrix,markovianStates,psiStates); + calculateWu(relativeReachability, dir, k+1,to,lambda,wu,fullTransitionMatrix,markovianStates,psiStates, minMaxLinearEquationSolverFactory); } res+=element.getValue()*wu[k+1][to]; } } else { + + // applying simpleMDP stuff + auto numberOfStates=relativeReachability[0].size(); + typename storm::storage::SparseMatrix probMatrix = fullTransitionMatrix.getSubmatrix(true, ~markovianStates , ~markovianStates, true); + auto probSize = probMatrix.getRowGroupCount(); + std::vector b(probSize, 0), x(probSize,0); + + //calculate b + for (uint64_t i =0 ; i> solver = minMaxLinearEquationSolverFactory.create(probMatrix); + + solver->setHasUniqueSolution(); + solver->setBounds(storm::utility::zero(), storm::utility::one()); + solver->setRequirementsChecked(); + solver->setCachingEnabled(true); + + + solver->solveEquations(dir, x, b); + + for (uint64_t i =0 ; i::SupportsExponential, int>::type> - void SparseMarkovAutomatonCslHelper::calculateVd(OptimizationDirection dir, uint64_t k, uint64_t node, ValueType lambda, std::vector>& vd, storm::storage::SparseMatrix const& fullTransitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates){ + void SparseMarkovAutomatonCslHelper::calculateVd(std::vector> const& relativeReachability, OptimizationDirection dir, uint64_t k, uint64_t node, ValueType lambda, std::vector>& vd, storm::storage::SparseMatrix const& fullTransitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory){ std::ofstream logfile("U+logfile.txt", std::ios::app); @@ -324,11 +384,83 @@ namespace storm { for (auto &element : line){ uint64_t to = element.getColumn(); if (vd[k+1][to]==-1){ - calculateVd(dir,k+1,to,lambda,vd, fullTransitionMatrix, markovianStates,psiStates); + calculateVd(relativeReachability, dir,k+1,to,lambda,vd, fullTransitionMatrix, markovianStates,psiStates, minMaxLinearEquationSolverFactory); } res+=element.getValue()*vd[k+1][to]; } + + vd[k][node]=res; + logfile << " res = " << res << "\n"; } else { //no-goal prob state + // applying simpleMDP stuff + logfile << "probState, applying MDP stuff "; + auto numberOfStates=relativeReachability[0].size(); + typename storm::storage::SparseMatrix probMatrix = fullTransitionMatrix.getSubmatrix(true, ~markovianStates , ~markovianStates, true); + auto probSize = probMatrix.getRowGroupCount(); + std::vector b(probSize, 0), x(probSize,0); + + //calculate b + for (uint64_t i =0 ; i> solver = minMaxLinearEquationSolverFactory.create(probMatrix); + + solver->setHasUniqueSolution(); + solver->setBounds(storm::utility::zero(), storm::utility::one()); + solver->setRequirementsChecked(); + solver->setCachingEnabled(true); + + + solver->solveEquations(dir, x, b); + logfile << "goal vector x is \n"; + for (int i =0 ; i::SupportsExponential, int>::type=0> @@ -482,7 +613,7 @@ 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& markovianStates, storm::storage::BitVector const& psiStates){ + std::vector SparseMarkovAutomatonCslHelper::unifPlus(OptimizationDirection dir, std::pair const& boundsPair, std::vector const& exitRateVector, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& markovianStates, 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); @@ -516,24 +647,32 @@ namespace storm { uint64_t N; //calculate relative ReachabilityVectors - std::vector in(fullTransitionMatrix.getRowCount(), 0); - std::vector> relReachability(numberOfStates,in); + std::vector in(numberOfStates, 0); + std::vector> relReachability(fullTransitionMatrix.getRowCount(),in); printTransitions(relReachability, exitRateVector, fullTransitionMatrix, markovianStates, psiStates, psiStates, psiStates, vd,vu,wu); // TODO: delete when develepmont is finished - for(int i=0 ; i& act = relReachability[j]; for(auto element: fullTransitionMatrix.getRow(j)){ if (markovianStates[element.getColumn()]){ - act[element.getColumn]=element.getValue(); + act[element.getColumn()]=element.getValue(); } } - } } + + + printTransitions(relReachability, exitRateVector, fullTransitionMatrix, markovianStates, psiStates, psiStates, psiStates, vd,vu,wu); // TODO: delete when develepmont is finished // while not close enough to precision: do { // (2) update parameter @@ -577,15 +716,17 @@ namespace storm { // (5) calculate vectors and maxNorm for (uint64_t i = 0; i < numberOfStates; i++) { for (uint64_t k = N; k <= N; k--) { - calculateVd(dir, k, i, T*lambda, vd, fullTransitionMatrix, markovianStates, psiStates); - calculateWu(dir, k, i, T*lambda, wu, fullTransitionMatrix, markovianStates, psiStates); - calculateVu(dir, k, i, T*lambda, vu, wu, fullTransitionMatrix, markovianStates, psiStates); + calculateVd(relReachability, dir, k, i, T*lambda, vd, fullTransitionMatrix, markovianStates, psiStates, minMaxLinearEquationSolverFactory); + calculateWu(relReachability, dir, k, i, T*lambda, wu, fullTransitionMatrix, markovianStates, psiStates, minMaxLinearEquationSolverFactory); + calculateVu(relReachability, dir, k, i, T*lambda, vu, wu, fullTransitionMatrix, markovianStates, psiStates, minMaxLinearEquationSolverFactory); //also use iteration to keep maxNorm of vd and vu up to date, so the loop-condition is easy to prove ValueType diff = std::abs(vd[k][i]-vu[k][i]); maxNorm = std::max(maxNorm, diff); } } + printTransitions(relReachability, exitRateVector, fullTransitionMatrix, markovianStates, psiStates, psiStates, psiStates, vd,vu,wu); // TODO: delete when develepmont is finished + // (6) double lambda lambda=2*lambda; @@ -668,8 +809,7 @@ namespace storm { } else { STORM_LOG_ASSERT(markovAutomatonSettings.getTechnique() == storm::settings::modules::MarkovAutomatonSettings::BoundedReachabilityTechnique::UnifPlus, "Unknown solution technique."); - // Why is optimization direction not passed? - return unifPlus(dir, boundsPair, exitRateVector, transitionMatrix, markovianStates, psiStates); + return unifPlus(dir, boundsPair, exitRateVector, transitionMatrix, markovianStates, psiStates, minMaxLinearEquationSolverFactory); } } diff --git a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h index 376b4e88a..7dac4cf6f 100644 --- a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h +++ b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h @@ -28,7 +28,7 @@ namespace storm { * */ template ::SupportsExponential, int>::type=0> - static std::vector unifPlus(OptimizationDirection dir, std::pair const& boundsPair, std::vector const& exitRateVector, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates); + static std::vector unifPlus(OptimizationDirection dir, std::pair const& boundsPair, std::vector const& exitRateVector, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); template ::SupportsExponential, int>::type = 0> static std::vector computeBoundedUntilProbabilities(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, std::pair const& boundsPair, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); @@ -98,14 +98,14 @@ namespace storm { * */ template ::SupportsExponential, int>::type=0> - static void calculateVd(OptimizationDirection dir, uint64_t k, uint64_t node, ValueType lambda, std::vector>& vd, storm::storage::SparseMatrix const& fullTransitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates); + static void calculateVd(std::vector> const& relativeReachability, OptimizationDirection dir, uint64_t k, uint64_t node, ValueType lambda, std::vector>& vd, storm::storage::SparseMatrix const& fullTransitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); /*! * Computes vu vector according to UnifPlus * */ template ::SupportsExponential, int>::type=0> - static void calculateVu(OptimizationDirection dir, uint64_t k, uint64_t node, ValueType lambda, std::vector>& vu, std::vector>& wu, storm::storage::SparseMatrix const& fullTransitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates); + static void calculateVu(std::vector> const& relativeReachability, OptimizationDirection dir, uint64_t k, uint64_t node, ValueType lambda, std::vector>& vu, std::vector>& wu, storm::storage::SparseMatrix const& fullTransitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); /*! @@ -113,7 +113,7 @@ namespace storm { * */ template ::SupportsExponential, int>::type=0> - static void calculateWu(OptimizationDirection dir, uint64_t k, uint64_t node, ValueType lambda, std::vector>& wu, storm::storage::SparseMatrix const& fullTransitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates); + static void calculateWu(std::vector> const& relativeReachability, OptimizationDirection dir, uint64_t k, uint64_t node, ValueType lambda, std::vector>& wu, storm::storage::SparseMatrix const& fullTransitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); /*! * Prints the TransitionMatrix and the vectors vd, vu, wu to the logfile