diff --git a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp index c176330c3..dbd60cd39 100644 --- a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp +++ b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp @@ -37,28 +37,47 @@ namespace storm { template struct UnifPlusVectors { - std::vector> resLower; - std::vector> resUpper; + UnifPlusVectors(uint64_t steps, uint64_t noStates) : numberOfStates(noStates), resLowerOld(numberOfStates, -1), resLowerNew(numberOfStates, -1), resUpperOld(numberOfStates, -1), resUpperNew(numberOfStates, -1) { + // Intentionally empty. + wUpper = std::vector>(steps, std::vector(numberOfStates, -1)); + } + + /** + * 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); + resUpperOld.swap(resUpperNew); + std::fill(resUpperNew.begin(), resUpperNew.end(), -1); + } + + uint64_t numberOfStates; + std::vector resLowerOld; + std::vector resLowerNew; + std::vector resUpperOld; + std::vector resUpperNew; std::vector> wUpper; }; template 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) { - std::vector>& resVector = calcLower ? unifVectors.resLower : unifVectors.wUpper; + std::vector& resVectorOld = calcLower ? unifVectors.resLowerOld : unifVectors.wUpper[k+1]; + std::vector& resVectorNew = calcLower ? unifVectors.resLowerNew : unifVectors.wUpper[k]; - if (resVector[k][state] != -1) { + if (resVectorNew[state] != -1) { // Result already calculated. return; } auto numberOfStates = fullTransitionMatrix.getRowGroupCount(); - uint64_t N = resVector.size() - 1; + uint64_t N = unifVectors.wUpper.size() - 1; auto const& rowGroupIndices = fullTransitionMatrix.getRowGroupIndices(); ValueType res; // First case, k==N, independent from kind of state. if (k == N) { - resVector[k][state] = storm::utility::zero(); + resVectorNew[state] = storm::utility::zero(); return; } @@ -72,10 +91,10 @@ namespace storm { res += poisson.weights[i - poisson.left]; } } - resVector[k][state] = res; + resVectorNew[state] = res; } else { // w upper - resVector[k][state] = storm::utility::one(); + resVectorNew[state] = storm::utility::one(); } return; } @@ -86,12 +105,13 @@ namespace storm { res = storm::utility::zero(); for (auto const& element : fullTransitionMatrix.getRow(rowGroupIndices[state])) { uint64_t successor = element.getColumn(); - if (resVector[k+1][successor] == -1) { + 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() * resVector[k+1][successor]; + res += element.getValue() * resVectorOld[successor]; } - resVector[k][state]=res; + resVectorNew[state]=res; return; } @@ -110,10 +130,10 @@ namespace storm { continue; } - if (resVector[k][successor] == -1) { + if (resVectorNew[successor] == -1) { calculateUnifPlusVector(env, k, successor, calcLower, lambda, numberOfProbabilisticChoices, relativeReachability, dir, unifVectors, fullTransitionMatrix, markovianStates, psiStates, solver, poisson, cycleFree); } - between += element.getValue() * resVector[k][successor]; + between += element.getValue() * resVectorNew[successor]; } if (maximize(dir)) { res = storm::utility::max(res, between); @@ -125,7 +145,7 @@ namespace storm { } } } - resVector[k][state] = res; + resVectorNew[state] = res; return; } @@ -150,10 +170,10 @@ namespace storm { continue; } - if (resVector[k][successor] == -1) { + if (resVectorNew[successor] == -1) { calculateUnifPlusVector(env, k, successor, calcLower, lambda, numberOfProbabilisticStates, relativeReachability, dir, unifVectors, fullTransitionMatrix, markovianStates, psiStates, solver, poisson, cycleFree); } - res += relativeReachability[j][stateCount] * resVector[k][successor]; + res += relativeReachability[j][stateCount] * resVectorNew[successor]; ++stateCount; } @@ -166,28 +186,30 @@ namespace storm { solver->solveEquations(env, dir, x, b); // Expand the solution for the probabilistic states to all states. - storm::utility::vector::setVectorValues(resVector[k], ~markovianStates, x); + storm::utility::vector::setVectorValues(resVectorNew, ~markovianStates, x); } template void calculateResUpper(Environment const& env, std::vector> const& relativeReachability, OptimizationDirection dir, uint64_t k, uint64_t state, ValueType lambda, uint64_t numberOfProbabilisticStates, 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) { // Avoiding multiple computation of the same value. - if (unifVectors.resUpper[k][state] != -1) { + if (unifVectors.resUpperNew[state] != -1) { + STORM_LOG_ASSERT(false, "Result was already calculated."); return; } - uint64_t N = unifVectors.resUpper.size() - 1; + uint64_t N = unifVectors.wUpper.size() - 1; ValueType res = storm::utility::zero(); for (uint64_t i = k; i < N; ++i) { if (unifVectors.wUpper[N-1-(i-k)][state] == -1) { - calculateUnifPlusVector(env, N-1-(i-k), state, false, lambda, numberOfProbabilisticStates, relativeReachability, dir, unifVectors, fullTransitionMatrix, markovianStates, psiStates, solver, poisson, cycleFree); + STORM_LOG_ASSERT(false, "Need to calculate previous result."); + calculateUnifPlusVector(env, N-1-(i-k), state, false, 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.wUpper[N-1-(i-k)][state]; } } - unifVectors.resUpper[k][state] = res; + unifVectors.resUpperNew[state] = res; } template @@ -237,7 +259,7 @@ namespace storm { bool cycleFree = sccDecomposition.empty(); // Vectors to store computed vectors. - UnifPlusVectors unifVectors; + UnifPlusVectors unifVectors(0, 0); // 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. @@ -313,7 +335,6 @@ namespace storm { } } - std::vector init(numberOfStates, -1); ValueType maxNorm = storm::utility::zero(); // Maximal step size uint64_t N; @@ -362,10 +383,7 @@ namespace storm { } // (4) Define vectors/matrices. - std::vector> v = std::vector>(N + 1, init); - unifVectors.resLower = v; - unifVectors.resUpper = v; - unifVectors.wUpper = v; + unifVectors = UnifPlusVectors(N+1, numberOfStates); // (5) Compute vectors and maxNorm. for (int64_t k = N; k >= 0; --k) { @@ -374,12 +392,13 @@ namespace storm { calculateUnifPlusVector(env, k, i, false, lambda, numberOfProbabilisticChoices, relativeReachabilities, dir, unifVectors, fullTransitionMatrix, markovianAndGoalStates, psiStates, solver, foxGlynnResult, cycleFree); calculateResUpper(env, relativeReachabilities, dir, k, i, lambda, numberOfProbabilisticChoices, unifVectors, fullTransitionMatrix, markovianAndGoalStates, psiStates, solver, foxGlynnResult, cycleFree); } + unifVectors.prepareNewIteration(); } // 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.resUpper[0][i] - unifVectors.resLower[0][i]); + ValueType diff = storm::utility::abs(unifVectors.resUpperOld[i] - unifVectors.resLowerOld[i]); maxNorm = std::max(maxNorm, diff); } @@ -389,7 +408,7 @@ namespace storm { } while (maxNorm > epsilon * (1 - kappa)); - return unifVectors.resLower[0]; + return unifVectors.resLowerOld; } template