diff --git a/CHANGELOG.md b/CHANGELOG.md index 2f9274625..118c4c8c5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -18,6 +18,7 @@ Version 1.3.x - New binary `storm-conv` that handles conversions between model files - New binary `storm-pomdp` that handles the translation of POMDPs to pMCs. - Maximal progress assumption is now applied while building Markov Automata (sparse engine). +- Improved Unif+ implementation for Markov Automata, significantly reduced memory consumption. - Added support for expected time properties for discrete time models - Bug fix in the parser for DRN (MDPs and MAs might have been affected). - `storm-gspn`: Improved .pnpro parser diff --git a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp index c78948727..b2fcad350 100644 --- a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp +++ b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp @@ -34,56 +34,96 @@ namespace storm { namespace modelchecker { namespace helper { + + /** + * Data structure holding result vectors (vLower, vUpper, wUpper) for Unif+. + */ + template + struct UnifPlusVectors { + UnifPlusVectors() { + // Intentionally empty + } + + /** + * Initialize results vectors. vLowerOld, vUpperOld and wUpper[k=N] are initialized with zeros. + */ + UnifPlusVectors(uint64_t steps, uint64_t noStates) : numberOfStates(noStates), steps(steps), resLowerOld(numberOfStates, storm::utility::zero()), resLowerNew(numberOfStates, -1), resUpper(numberOfStates, storm::utility::zero()), wUpperOld(numberOfStates, storm::utility::zero()), wUpperNew(numberOfStates, -1) { + // Intentionally left empty + } + + /** + * Prepare new iteration by setting the new result vectors as old result vectors, and initializing the new result vectors with -1 again. + */ + void prepareNewIteration() { + resLowerOld.swap(resLowerNew); + std::fill(resLowerNew.begin(), resLowerNew.end(), -1); + wUpperOld.swap(wUpperNew); + std::fill(wUpperNew.begin(), wUpperNew.end(), -1); + } + + uint64_t numberOfStates; + uint64_t steps; + std::vector resLowerOld; + std::vector resLowerNew; + std::vector resUpper; + std::vector wUpperOld; + std::vector wUpperNew; + }; template - void calculateUnifPlusVector(Environment const& env, uint64_t k, uint64_t state, uint64_t const kind, ValueType lambda, uint64_t numberOfProbabilisticChoices, 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, storm::utility::numerical::FoxGlynnResult const& poisson, bool cycleFree) { - - if (unifVectors[kind][k][state] != -1) { + void calculateUnifPlusVector(Environment const& env, uint64_t k, uint64_t state, bool calcLower, ValueType lambda, uint64_t numberOfProbabilisticChoices, std::vector> const & relativeReachability, OptimizationDirection dir, UnifPlusVectors& unifVectors, storm::storage::SparseMatrix const& fullTransitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, std::unique_ptr> const& solver, storm::utility::numerical::FoxGlynnResult const& poisson, bool cycleFree) { + // Set reference to acutal vector + std::vector& resVectorOld = calcLower ? unifVectors.resLowerOld : unifVectors.wUpperOld; + std::vector& resVectorNew = calcLower ? unifVectors.resLowerNew : unifVectors.wUpperNew; + + if (resVectorNew[state] != -1) { // Result already calculated. return; } auto numberOfStates = fullTransitionMatrix.getRowGroupCount(); - uint64_t N = unifVectors[kind].size() - 1; + uint64_t N = unifVectors.steps; auto const& rowGroupIndices = fullTransitionMatrix.getRowGroupIndices(); ValueType res; // First case, k==N, independent from kind of state. if (k == N) { - unifVectors[kind][k][state] = storm::utility::zero(); + STORM_LOG_ASSERT(false, "Result for k=N was already calculated."); + resVectorNew[state] = storm::utility::zero(); return; } // Goal state, independent from kind of state. if (psiStates[state]) { - if (kind == 0) { - // Vd + if (calcLower) { + // v lower res = storm::utility::zero(); for (uint64_t i = k; i < N; ++i){ if (i >= poisson.left && i <= poisson.right) { - ValueType between = poisson.weights[i - poisson.left]; - res += between; + res += poisson.weights[i - poisson.left]; } } - unifVectors[kind][k][state] = res; + resVectorNew[state] = res; } else { - // WU - unifVectors[kind][k][state] = storm::utility::one(); + // w upper + resVectorNew[state] = storm::utility::one(); } return; } - + + // Markovian non-goal state. if (markovianStates[state]) { res = storm::utility::zero(); for (auto const& element : fullTransitionMatrix.getRow(rowGroupIndices[state])) { - uint64_t to = element.getColumn(); - if (unifVectors[kind][k+1][to] == -1) { - calculateUnifPlusVector(env, k+1, to, kind, lambda, numberOfProbabilisticChoices, relativeReachability, dir, unifVectors, fullTransitionMatrix, markovianStates, psiStates, solver, poisson, cycleFree); + uint64_t successor = element.getColumn(); + if (resVectorOld[successor] == -1) { + STORM_LOG_ASSERT(false, "Need to calculate previous result."); + calculateUnifPlusVector(env, k+1, successor, calcLower, lambda, numberOfProbabilisticChoices, relativeReachability, dir, unifVectors, fullTransitionMatrix, markovianStates, psiStates, solver, poisson, cycleFree); } - res += element.getValue()*unifVectors[kind][k+1][to]; + res += element.getValue() * resVectorOld[successor]; } - unifVectors[kind][k][state]=res; + resVectorNew[state]=res; return; } @@ -93,7 +133,7 @@ namespace storm { res = -1; for (uint64_t i = rowGroupIndices[state]; i < rowGroupIndices[state + 1]; ++i) { auto row = fullTransitionMatrix.getRow(i); - ValueType between = 0; + ValueType between = storm::utility::zero(); for (auto const& element : row) { uint64_t successor = element.getColumn(); @@ -102,10 +142,10 @@ namespace storm { continue; } - if (unifVectors[kind][k][successor] == -1) { - calculateUnifPlusVector(env, k, successor, kind, lambda, numberOfProbabilisticChoices, relativeReachability, dir, unifVectors, fullTransitionMatrix, markovianStates, psiStates, solver, poisson, cycleFree); + if (resVectorNew[successor] == -1) { + calculateUnifPlusVector(env, k, successor, calcLower, lambda, numberOfProbabilisticChoices, relativeReachability, dir, unifVectors, fullTransitionMatrix, markovianStates, psiStates, solver, poisson, cycleFree); } - between += element.getValue() * unifVectors[kind][k][successor]; + between += element.getValue() * resVectorNew[successor]; } if (maximize(dir)) { res = storm::utility::max(res, between); @@ -117,11 +157,11 @@ namespace storm { } } } - unifVectors[kind][k][state] = res; + resVectorNew[state] = res; return; } - // If we arrived at this point, the model is not cycle free. Use the solver to solve the unterlying equation system. + // If we arrived at this point, the model is not cycle free. Use the solver to solve the underlying equation system. uint64_t numberOfProbabilisticStates = numberOfStates - markovianStates.getNumberOfSetBits(); std::vector b(numberOfProbabilisticChoices, storm::utility::zero()); std::vector x(numberOfProbabilisticStates, storm::utility::zero()); @@ -142,10 +182,10 @@ namespace storm { continue; } - if (unifVectors[kind][k][successor] == -1) { - calculateUnifPlusVector(env, k, successor, kind, lambda, numberOfProbabilisticStates, relativeReachability, dir, unifVectors, fullTransitionMatrix, markovianStates, psiStates, solver, poisson, cycleFree); + if (resVectorNew[successor] == -1) { + calculateUnifPlusVector(env, k, successor, calcLower, lambda, numberOfProbabilisticStates, relativeReachability, dir, unifVectors, fullTransitionMatrix, markovianStates, psiStates, solver, poisson, cycleFree); } - res = res + relativeReachability[j][stateCount] * unifVectors[kind][k][successor]; + res += relativeReachability[j][stateCount] * resVectorNew[successor]; ++stateCount; } @@ -158,28 +198,7 @@ namespace storm { solver->solveEquations(env, dir, x, b); // Expand the solution for the probabilistic states to all states. - storm::utility::vector::setVectorValues(unifVectors[kind][k], ~markovianStates, x); - } - - template - void calculateVu(Environment const& env, std::vector> const& relativeReachability, OptimizationDirection dir, uint64_t k, uint64_t state, uint64_t const kind, ValueType lambda, uint64_t numberOfProbabilisticStates, std::vector>>& unifVectors, storm::storage::SparseMatrix const& fullTransitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, std::unique_ptr> const& solver, storm::utility::numerical::FoxGlynnResult const & poisson, bool cycleFree) { - - // Avoiding multiple computation of the same value. - if (unifVectors[1][k][state] != -1) { - return; - } - uint64_t N = unifVectors[1].size() - 1; - - ValueType res = storm::utility::zero(); - for (uint64_t i = k; i < N; ++i) { - if (unifVectors[2][N-1-(i-k)][state] == -1) { - calculateUnifPlusVector(env, N-1-(i-k), state, 2, lambda, numberOfProbabilisticStates, relativeReachability, dir, unifVectors, fullTransitionMatrix, markovianStates, psiStates, solver, poisson, cycleFree); - } - if (i >= poisson.left && i <= poisson.right) { - res += poisson.weights[i - poisson.left] * unifVectors[2][N-1-(i-k)][state]; - } - } - unifVectors[1][k][state] = res; + storm::utility::vector::setVectorValues(resVectorNew, ~markovianStates, x); } template @@ -227,9 +246,9 @@ namespace storm { // Searching for SCCs in probabilistic fragment to decide which algorithm is applied. storm::storage::StronglyConnectedComponentDecomposition sccDecomposition(transitionMatrix, probabilisticStates, true, false); bool cycleFree = sccDecomposition.empty(); - + // Vectors to store computed vectors. - std::vector>> unifVectors(3); + UnifPlusVectors unifVectors; // Transitions from goal states will be ignored. However, we mark them as non-probabilistic to make sure // we do not apply the MDP algorithm to them. @@ -255,17 +274,22 @@ namespace storm { // (1) define/declare horizon, epsilon, kappa, N, lambda, maxNorm uint64_t numberOfStates = fullTransitionMatrix.getRowGroupCount(); - double T = boundsPair.second; + // 'Unpack' the bounds to make them more easily accessible. + double lowerBound = boundsPair.first; + double upperBound = boundsPair.second; + // Lower bound > 0 is not implemented! + STORM_LOG_THROW(lowerBound == 0, storm::exceptions::NotImplementedException, "Support for lower bound > 0 not implemented in Unif+."); + // Truncation error // TODO: make kappa a parameter. ValueType kappa = storm::utility::one() / 10; + // Approximation error ValueType epsilon = storm::settings::getModule().getPrecision(); + // Lambda is largest exit rate ValueType lambda = exitRateVector[0]; for (ValueType const& rate : exitRateVector) { lambda = std::max(rate, lambda); } - STORM_LOG_TRACE("Initial lambda is " << lambda << "."); - uint64_t N; - ValueType maxNorm = storm::utility::zero(); + STORM_LOG_DEBUG("Initial lambda is " << lambda << "."); // Compute the relative reachability vectors and create solver for models with SCCs. std::vector> relativeReachabilities(transitionMatrix.getRowCount()); @@ -300,16 +324,19 @@ namespace storm { } } + ValueType maxNorm = storm::utility::zero(); + // Maximal step size + uint64_t N; + storm::utility::ProgressMeasurement progressIterations("iterations"); + size_t iteration = 0; + progressIterations.startNewMeasurement(iteration); // Loop until result is within precision bound. - std::vector init(numberOfStates, -1); do { - maxNorm = storm::utility::zero(); - // (2) update parameter - N = storm::utility::ceil(lambda * T * std::exp(2) - storm::utility::log(kappa * epsilon)); + N = storm::utility::ceil(lambda * upperBound * std::exp(2) - storm::utility::log(kappa * epsilon)); // (3) uniform - just applied to Markovian states. - for (uint64_t i = 0; i < fullTransitionMatrix.getRowGroupCount(); i++) { + for (uint64_t i = 0; i < numberOfStates; i++) { if (!markovianAndGoalStates[i] || psiStates[i]) { continue; } @@ -340,7 +367,7 @@ namespace storm { } // Compute poisson distribution. - storm::utility::numerical::FoxGlynnResult foxGlynnResult = storm::utility::numerical::foxGlynn(lambda * T, epsilon * kappa / 100); + storm::utility::numerical::FoxGlynnResult foxGlynnResult = storm::utility::numerical::foxGlynn(lambda * upperBound, epsilon * kappa / 100); // Scale the weights so they sum to one. for (auto& element : foxGlynnResult.weights) { @@ -348,35 +375,48 @@ namespace storm { } // (4) Define vectors/matrices. - std::vector> v = std::vector>(N + 1, init); - - unifVectors[0] = v; - unifVectors[1] = v; - unifVectors[2] = v; + // Initialize result vectors and already insert zeros for iteration N + unifVectors = UnifPlusVectors(N, numberOfStates); - // Define 0=vd, 1=vu, 2=wu. // (5) Compute vectors and maxNorm. - for (uint64_t i = 0; i < numberOfStates; ++i) { - for (uint64_t k = N; k <= N; --k) { - calculateUnifPlusVector(env, k, i, 0, lambda, numberOfProbabilisticChoices, relativeReachabilities, dir, unifVectors, fullTransitionMatrix, markovianAndGoalStates, psiStates, solver, foxGlynnResult, cycleFree); - calculateUnifPlusVector(env, k, i, 2, lambda, numberOfProbabilisticChoices, relativeReachabilities, dir, unifVectors, fullTransitionMatrix, markovianAndGoalStates, psiStates, solver, foxGlynnResult, cycleFree); - calculateVu(env, relativeReachabilities, dir, k, i, 1, lambda, numberOfProbabilisticChoices, unifVectors, fullTransitionMatrix, markovianAndGoalStates, psiStates, solver, foxGlynnResult, cycleFree); + // Iteration k = N was already performed by initializing with zeros. + + // Iterations k < N + storm::utility::ProgressMeasurement progressSteps("steps in iteration " + std::to_string(iteration)); + progressSteps.setMaxCount(N); + progressSteps.startNewMeasurement(0); + for (int64_t k = N-1; k >= 0; --k) { + if (k < (int64_t)(N-1)) { + unifVectors.prepareNewIteration(); + } + for (uint64_t state = 0; state < numberOfStates; ++state) { + // Calculate results for lower bound and wUpper + calculateUnifPlusVector(env, k, state, true, lambda, numberOfProbabilisticChoices, relativeReachabilities, dir, unifVectors, fullTransitionMatrix, markovianAndGoalStates, psiStates, solver, foxGlynnResult, cycleFree); + calculateUnifPlusVector(env, k, state, false, lambda, numberOfProbabilisticChoices, relativeReachabilities, dir, unifVectors, fullTransitionMatrix, markovianAndGoalStates, psiStates, solver, foxGlynnResult, cycleFree); + // Calculate result for upper bound + uint64_t index = N-1-k; + if (index >= foxGlynnResult.left && index <= foxGlynnResult.right) { + STORM_LOG_ASSERT(unifVectors.wUpperNew[state] != -1, "wUpper was not computed before."); + unifVectors.resUpper[state] += foxGlynnResult.weights[index - foxGlynnResult.left] * unifVectors.wUpperNew[state]; + } } + progressSteps.updateProgress(N-k); } // Only iterate over result vector, as the results can only get more precise. + maxNorm = storm::utility::zero(); for (uint64_t i = 0; i < numberOfStates; i++){ - ValueType diff = storm::utility::abs(unifVectors[0][0][i] - unifVectors[1][0][i]); + ValueType diff = storm::utility::abs(unifVectors.resUpper[i] - unifVectors.resLowerNew[i]); maxNorm = std::max(maxNorm, diff); } // (6) Double lambda. lambda *= 2; - STORM_LOG_TRACE("Increased lambda to " << lambda << ", max diff is " << maxNorm << "."); - + STORM_LOG_DEBUG("Increased lambda to " << lambda << ", max diff is " << maxNorm << "."); + progressIterations.updateProgress(++iteration); } while (maxNorm > epsilon * (1 - kappa)); - return unifVectors[0][0]; + return unifVectors.resLowerNew; } template @@ -569,7 +609,12 @@ namespace storm { return computeBoundedUntilProbabilitiesImca(env, dir, transitionMatrix, exitRateVector, markovianStates, psiStates, boundsPair); } else { STORM_LOG_ASSERT(settings.getMarkovAutomatonBoundedReachabilityMethod() == storm::settings::modules::MinMaxEquationSolverSettings::MarkovAutomatonBoundedReachabilityMethod::UnifPlus, "Unknown solution method."); - return computeBoundedUntilProbabilitiesUnifPlus(env, dir, transitionMatrix, exitRateVector, markovianStates, psiStates, boundsPair); + if (!storm::utility::isZero(boundsPair.first)) { + STORM_LOG_WARN("Using IMCA method because Unif+ does not support a lower bound > 0."); + return computeBoundedUntilProbabilitiesImca(env, dir, transitionMatrix, exitRateVector, markovianStates, psiStates, boundsPair); + } else { + return computeBoundedUntilProbabilitiesUnifPlus(env, dir, transitionMatrix, exitRateVector, markovianStates, psiStates, boundsPair); + } } }