From 6066ecd5909cceffca33e57a8a1c8a1f97f11735 Mon Sep 17 00:00:00 2001
From: Matthias Volk <matthias.volk@cs.rwth-aachen.de>
Date: Wed, 5 Dec 2018 18:35:54 +0100
Subject: [PATCH 01/12] Added struct for Unif+ vectors
---
.../helper/SparseMarkovAutomatonCslHelper.cpp | 86 ++++++++++---------
1 file changed, 46 insertions(+), 40 deletions(-)
diff --git a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
index a0c8829bc..4d0aeb271 100644
--- a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
+++ b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
@@ -34,29 +34,37 @@
namespace storm {
namespace modelchecker {
namespace helper {
+
+ template<typename ValueType>
+ struct UnifPlusVectors {
+ std::vector<std::vector<ValueType>> resLower;
+ std::vector<std::vector<ValueType>> resUpper;
+ std::vector<std::vector<ValueType>> wUpper;
+ };
template<typename ValueType>
- void calculateUnifPlusVector(Environment const& env, uint64_t k, uint64_t state, uint64_t const kind, ValueType lambda, uint64_t numberOfProbabilisticChoices, std::vector<std::vector<ValueType>> const & relativeReachability, OptimizationDirection dir, std::vector<std::vector<std::vector<ValueType>>>& unifVectors, storm::storage::SparseMatrix<ValueType> const& fullTransitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> const& solver, storm::utility::numerical::FoxGlynnResult<ValueType> 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<std::vector<ValueType>> const & relativeReachability, OptimizationDirection dir, UnifPlusVectors<ValueType>& unifVectors, storm::storage::SparseMatrix<ValueType> const& fullTransitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> const& solver, storm::utility::numerical::FoxGlynnResult<ValueType> const& poisson, bool cycleFree) {
+ std::vector<std::vector<ValueType>>& resVector = calcLower ? unifVectors.resLower : unifVectors.wUpper;
+
+ if (resVector[k][state] != -1) {
// Result already calculated.
return;
}
auto numberOfStates = fullTransitionMatrix.getRowGroupCount();
- uint64_t N = unifVectors[kind].size() - 1;
+ uint64_t N = resVector.size() - 1;
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<ValueType>();
+ resVector[k][state] = storm::utility::zero<ValueType>();
return;
}
// Goal state, independent from kind of state.
if (psiStates[state]) {
- if (kind == 0) {
+ if (calcLower) {
// Vd
res = storm::utility::zero<ValueType>();
for (uint64_t i = k; i < N; ++i){
@@ -65,10 +73,10 @@ namespace storm {
res += between;
}
}
- unifVectors[kind][k][state] = res;
+ resVector[k][state] = res;
} else {
// WU
- unifVectors[kind][k][state] = storm::utility::one<ValueType>();
+ resVector[k][state] = storm::utility::one<ValueType>();
}
return;
}
@@ -78,12 +86,12 @@ namespace storm {
res = storm::utility::zero<ValueType>();
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);
+ if (resVector[k+1][to] == -1) {
+ calculateUnifPlusVector(env, k+1, to, calcLower, lambda, numberOfProbabilisticChoices, relativeReachability, dir, unifVectors, fullTransitionMatrix, markovianStates, psiStates, solver, poisson, cycleFree);
}
- res += element.getValue()*unifVectors[kind][k+1][to];
+ res += element.getValue()*resVector[k+1][to];
}
- unifVectors[kind][k][state]=res;
+ resVector[k][state]=res;
return;
}
@@ -102,10 +110,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 (resVector[k][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() * resVector[k][successor];
}
if (maximize(dir)) {
res = storm::utility::max(res, between);
@@ -117,7 +125,7 @@ namespace storm {
}
}
}
- unifVectors[kind][k][state] = res;
+ resVector[k][state] = res;
return;
}
@@ -142,10 +150,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 (resVector[k][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 = res + relativeReachability[j][stateCount] * resVector[k][successor];
++stateCount;
}
@@ -158,28 +166,28 @@ 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);
+ storm::utility::vector::setVectorValues(resVector[k], ~markovianStates, x);
}
template <typename ValueType>
- void calculateVu(Environment const& env, std::vector<std::vector<ValueType>> const& relativeReachability, OptimizationDirection dir, uint64_t k, uint64_t state, uint64_t const kind, ValueType lambda, uint64_t numberOfProbabilisticStates, std::vector<std::vector<std::vector<ValueType>>>& unifVectors, storm::storage::SparseMatrix<ValueType> const& fullTransitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> const& solver, storm::utility::numerical::FoxGlynnResult<ValueType> const & poisson, bool cycleFree) {
+ void calculateResUpper(Environment const& env, std::vector<std::vector<ValueType>> const& relativeReachability, OptimizationDirection dir, uint64_t k, uint64_t state, ValueType lambda, uint64_t numberOfProbabilisticStates, UnifPlusVectors<ValueType>& unifVectors, storm::storage::SparseMatrix<ValueType> const& fullTransitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> const& solver, storm::utility::numerical::FoxGlynnResult<ValueType> const & poisson, bool cycleFree) {
// Avoiding multiple computation of the same value.
- if (unifVectors[1][k][state] != -1) {
+ if (unifVectors.resUpper[k][state] != -1) {
return;
}
- uint64_t N = unifVectors[1].size() - 1;
+ uint64_t N = unifVectors.resUpper.size() - 1;
ValueType res = storm::utility::zero<ValueType>();
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 (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);
}
if (i >= poisson.left && i <= poisson.right) {
- res += poisson.weights[i - poisson.left] * unifVectors[2][N-1-(i-k)][state];
+ res += poisson.weights[i - poisson.left] * unifVectors.wUpper[N-1-(i-k)][state];
}
}
- unifVectors[1][k][state] = res;
+ unifVectors.resUpper[k][state] = res;
}
template <typename ValueType>
@@ -227,9 +235,9 @@ namespace storm {
// Searching for SCCs in probabilistic fragment to decide which algorithm is applied.
storm::storage::StronglyConnectedComponentDecomposition<ValueType> sccDecomposition(transitionMatrix, probabilisticStates, true, false);
bool cycleFree = sccDecomposition.empty();
-
+
// Vectors to store computed vectors.
- std::vector<std::vector<std::vector<ValueType>>> unifVectors(3);
+ UnifPlusVectors<ValueType> 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.
@@ -349,24 +357,22 @@ namespace storm {
// (4) Define vectors/matrices.
std::vector<std::vector<ValueType>> v = std::vector<std::vector<ValueType>>(N + 1, init);
+ unifVectors.resLower = v;
+ unifVectors.resUpper = v;
+ unifVectors.wUpper = v;
- unifVectors[0] = v;
- unifVectors[1] = v;
- unifVectors[2] = v;
-
- // 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);
+ for (int64_t k = N; k >= 0; --k) {
+ calculateUnifPlusVector(env, k, i, true, lambda, numberOfProbabilisticChoices, relativeReachabilities, dir, unifVectors, fullTransitionMatrix, markovianAndGoalStates, psiStates, solver, foxGlynnResult, cycleFree);
+ 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);
}
}
// Only iterate over result vector, as the results can only get more precise.
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[0][i] - unifVectors.resLower[0][i]);
maxNorm = std::max(maxNorm, diff);
}
@@ -376,7 +382,7 @@ namespace storm {
} while (maxNorm > epsilon * (1 - kappa));
- return unifVectors[0][0];
+ return unifVectors.resLower[0];
}
template <typename ValueType>
From b8b2c58dabcc939eceb7b84cd68a54b058f3a81b Mon Sep 17 00:00:00 2001
From: Matthias Volk <matthias.volk@cs.rwth-aachen.de>
Date: Wed, 5 Dec 2018 18:36:45 +0100
Subject: [PATCH 02/12] Started on some refactoring in Unif+
---
.../helper/SparseMarkovAutomatonCslHelper.cpp | 25 ++++++++++++-------
1 file changed, 16 insertions(+), 9 deletions(-)
diff --git a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
index 4d0aeb271..03a445f36 100644
--- a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
+++ b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
@@ -263,17 +263,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<ValueType>() / 10;
+ // Approximation error
ValueType epsilon = storm::settings::getModule<storm::settings::modules::GeneralSettings>().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<ValueType>();
// Compute the relative reachability vectors and create solver for models with SCCs.
std::vector<std::vector<ValueType>> relativeReachabilities(transitionMatrix.getRowCount());
@@ -308,16 +313,17 @@ namespace storm {
}
}
- // Loop until result is within precision bound.
std::vector<ValueType> init(numberOfStates, -1);
+ ValueType maxNorm = storm::utility::zero<ValueType>();
+ // Maximal step size
+ uint64_t N;
+ // Loop until result is within precision bound.
do {
- maxNorm = storm::utility::zero<ValueType>();
-
// (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;
}
@@ -348,7 +354,7 @@ namespace storm {
}
// Compute poisson distribution.
- storm::utility::numerical::FoxGlynnResult<ValueType> foxGlynnResult = storm::utility::numerical::foxGlynn(lambda * T, epsilon * kappa / 100);
+ storm::utility::numerical::FoxGlynnResult<ValueType> foxGlynnResult = storm::utility::numerical::foxGlynn(lambda * upperBound, epsilon * kappa / 100);
// Scale the weights so they sum to one.
for (auto& element : foxGlynnResult.weights) {
@@ -371,6 +377,7 @@ namespace storm {
}
// Only iterate over result vector, as the results can only get more precise.
+ maxNorm = storm::utility::zero<ValueType>();
for (uint64_t i = 0; i < numberOfStates; i++){
ValueType diff = storm::utility::abs(unifVectors.resUpper[0][i] - unifVectors.resLower[0][i]);
maxNorm = std::max(maxNorm, diff);
From cbd61396139d7394dadc3944c5606014cc06f3f6 Mon Sep 17 00:00:00 2001
From: Matthias Volk <matthias.volk@cs.rwth-aachen.de>
Date: Thu, 6 Dec 2018 11:28:00 +0100
Subject: [PATCH 03/12] Small changes
---
.../helper/SparseMarkovAutomatonCslHelper.cpp | 24 +++++++++----------
1 file changed, 12 insertions(+), 12 deletions(-)
diff --git a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
index 03a445f36..6f948d59e 100644
--- a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
+++ b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
@@ -65,31 +65,31 @@ namespace storm {
// Goal state, independent from kind of state.
if (psiStates[state]) {
if (calcLower) {
- // Vd
+ // v lower
res = storm::utility::zero<ValueType>();
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];
}
}
resVector[k][state] = res;
} else {
- // WU
+ // w upper
resVector[k][state] = storm::utility::one<ValueType>();
}
return;
}
-
+
+
// Markovian non-goal state.
if (markovianStates[state]) {
res = storm::utility::zero<ValueType>();
for (auto const& element : fullTransitionMatrix.getRow(rowGroupIndices[state])) {
- uint64_t to = element.getColumn();
- if (resVector[k+1][to] == -1) {
- calculateUnifPlusVector(env, k+1, to, calcLower, lambda, numberOfProbabilisticChoices, relativeReachability, dir, unifVectors, fullTransitionMatrix, markovianStates, psiStates, solver, poisson, cycleFree);
+ uint64_t successor = element.getColumn();
+ if (resVector[k+1][successor] == -1) {
+ calculateUnifPlusVector(env, k+1, successor, calcLower, lambda, numberOfProbabilisticChoices, relativeReachability, dir, unifVectors, fullTransitionMatrix, markovianStates, psiStates, solver, poisson, cycleFree);
}
- res += element.getValue()*resVector[k+1][to];
+ res += element.getValue() * resVector[k+1][successor];
}
resVector[k][state]=res;
return;
@@ -101,7 +101,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<ValueType>();
for (auto const& element : row) {
uint64_t successor = element.getColumn();
@@ -129,7 +129,7 @@ namespace storm {
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<ValueType> b(numberOfProbabilisticChoices, storm::utility::zero<ValueType>());
std::vector<ValueType> x(numberOfProbabilisticStates, storm::utility::zero<ValueType>());
@@ -153,7 +153,7 @@ namespace storm {
if (resVector[k][successor] == -1) {
calculateUnifPlusVector(env, k, successor, calcLower, lambda, numberOfProbabilisticStates, relativeReachability, dir, unifVectors, fullTransitionMatrix, markovianStates, psiStates, solver, poisson, cycleFree);
}
- res = res + relativeReachability[j][stateCount] * resVector[k][successor];
+ res += relativeReachability[j][stateCount] * resVector[k][successor];
++stateCount;
}
From 8ae800b13057f11fd7a1077da5a0628743c03912 Mon Sep 17 00:00:00 2001
From: Matthias Volk <matthias.volk@cs.rwth-aachen.de>
Date: Thu, 6 Dec 2018 11:42:53 +0100
Subject: [PATCH 04/12] Changed iteration order to iterate over stepsize in
outer loop
---
.../csl/helper/SparseMarkovAutomatonCslHelper.cpp | 4 ++--
1 file changed, 2 insertions(+), 2 deletions(-)
diff --git a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
index 6f948d59e..c176330c3 100644
--- a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
+++ b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
@@ -368,8 +368,8 @@ namespace storm {
unifVectors.wUpper = v;
// (5) Compute vectors and maxNorm.
- for (uint64_t i = 0; i < numberOfStates; ++i) {
- for (int64_t k = N; k >= 0; --k) {
+ for (int64_t k = N; k >= 0; --k) {
+ for (uint64_t i = 0; i < numberOfStates; ++i) {
calculateUnifPlusVector(env, k, i, true, lambda, numberOfProbabilisticChoices, relativeReachabilities, dir, unifVectors, fullTransitionMatrix, markovianAndGoalStates, psiStates, solver, foxGlynnResult, cycleFree);
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);
From f9fb90499d4b29144345216eaf656dcfbee5cfe6 Mon Sep 17 00:00:00 2001
From: Matthias Volk <matthias.volk@cs.rwth-aachen.de>
Date: Thu, 6 Dec 2018 16:13:24 +0100
Subject: [PATCH 05/12] Only keep track of results from the last iteration
(instead of all iterations) for 2 of the 3 vectors
---
.../helper/SparseMarkovAutomatonCslHelper.cpp | 77 ++++++++++++-------
1 file changed, 48 insertions(+), 29 deletions(-)
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<typename ValueType>
struct UnifPlusVectors {
- std::vector<std::vector<ValueType>> resLower;
- std::vector<std::vector<ValueType>> 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<std::vector<ValueType>>(steps, std::vector<ValueType>(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<ValueType> resLowerOld;
+ std::vector<ValueType> resLowerNew;
+ std::vector<ValueType> resUpperOld;
+ std::vector<ValueType> resUpperNew;
std::vector<std::vector<ValueType>> wUpper;
};
template<typename ValueType>
void calculateUnifPlusVector(Environment const& env, uint64_t k, uint64_t state, bool calcLower, ValueType lambda, uint64_t numberOfProbabilisticChoices, std::vector<std::vector<ValueType>> const & relativeReachability, OptimizationDirection dir, UnifPlusVectors<ValueType>& unifVectors, storm::storage::SparseMatrix<ValueType> const& fullTransitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> const& solver, storm::utility::numerical::FoxGlynnResult<ValueType> const& poisson, bool cycleFree) {
- std::vector<std::vector<ValueType>>& resVector = calcLower ? unifVectors.resLower : unifVectors.wUpper;
+ std::vector<ValueType>& resVectorOld = calcLower ? unifVectors.resLowerOld : unifVectors.wUpper[k+1];
+ std::vector<ValueType>& 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<ValueType>();
+ resVectorNew[state] = storm::utility::zero<ValueType>();
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<ValueType>();
+ resVectorNew[state] = storm::utility::one<ValueType>();
}
return;
}
@@ -86,12 +105,13 @@ namespace storm {
res = storm::utility::zero<ValueType>();
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 <typename ValueType>
void calculateResUpper(Environment const& env, std::vector<std::vector<ValueType>> const& relativeReachability, OptimizationDirection dir, uint64_t k, uint64_t state, ValueType lambda, uint64_t numberOfProbabilisticStates, UnifPlusVectors<ValueType>& unifVectors, storm::storage::SparseMatrix<ValueType> const& fullTransitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> const& solver, storm::utility::numerical::FoxGlynnResult<ValueType> 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<ValueType>();
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 <typename ValueType>
@@ -237,7 +259,7 @@ namespace storm {
bool cycleFree = sccDecomposition.empty();
// Vectors to store computed vectors.
- UnifPlusVectors<ValueType> unifVectors;
+ UnifPlusVectors<ValueType> 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<ValueType> init(numberOfStates, -1);
ValueType maxNorm = storm::utility::zero<ValueType>();
// Maximal step size
uint64_t N;
@@ -362,10 +383,7 @@ namespace storm {
}
// (4) Define vectors/matrices.
- std::vector<std::vector<ValueType>> v = std::vector<std::vector<ValueType>>(N + 1, init);
- unifVectors.resLower = v;
- unifVectors.resUpper = v;
- unifVectors.wUpper = v;
+ unifVectors = UnifPlusVectors<ValueType>(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<ValueType>();
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 <typename ValueType>
From 6220b114b52a5ef09101f10a5587811662c4fb12 Mon Sep 17 00:00:00 2001
From: Matthias Volk <matthias.volk@cs.rwth-aachen.de>
Date: Thu, 6 Dec 2018 20:11:57 +0100
Subject: [PATCH 06/12] Small simplifications
---
.../helper/SparseMarkovAutomatonCslHelper.cpp | 81 ++++++++++---------
1 file changed, 44 insertions(+), 37 deletions(-)
diff --git a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
index dbd60cd39..3af9714b8 100644
--- a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
+++ b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
@@ -35,11 +35,23 @@ namespace storm {
namespace modelchecker {
namespace helper {
+ /**
+ * Data structure holding result vectors (vLower, vUpper, wUpper) for Unif+.
+ */
template<typename ValueType>
struct UnifPlusVectors {
- 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<std::vector<ValueType>>(steps, std::vector<ValueType>(numberOfStates, -1));
+ 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), resLowerOld(numberOfStates, storm::utility::zero<ValueType>()), resLowerNew(numberOfStates, -1), resUpperOld(numberOfStates, storm::utility::zero<ValueType>()), resUpperNew(numberOfStates, storm::utility::zero<ValueType>()) {
+ // For wUpper we have to keep track of all previous results
+ wUpper = std::vector<std::vector<ValueType>>(steps+1, std::vector<ValueType>(numberOfStates, -1));
+ // Initialize entries for step N with zeros
+ std::fill(wUpper[steps].begin(), wUpper[steps].end(), storm::utility::zero<ValueType>());
}
/**
@@ -49,7 +61,7 @@ namespace storm {
resLowerOld.swap(resLowerNew);
std::fill(resLowerNew.begin(), resLowerNew.end(), -1);
resUpperOld.swap(resUpperNew);
- std::fill(resUpperNew.begin(), resUpperNew.end(), -1);
+ std::fill(resUpperNew.begin(), resUpperNew.end(), storm::utility::zero<ValueType>());
}
uint64_t numberOfStates;
@@ -62,6 +74,7 @@ namespace storm {
template<typename ValueType>
void calculateUnifPlusVector(Environment const& env, uint64_t k, uint64_t state, bool calcLower, ValueType lambda, uint64_t numberOfProbabilisticChoices, std::vector<std::vector<ValueType>> const & relativeReachability, OptimizationDirection dir, UnifPlusVectors<ValueType>& unifVectors, storm::storage::SparseMatrix<ValueType> const& fullTransitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> const& solver, storm::utility::numerical::FoxGlynnResult<ValueType> const& poisson, bool cycleFree) {
+ // Set reference to acutal vector
std::vector<ValueType>& resVectorOld = calcLower ? unifVectors.resLowerOld : unifVectors.wUpper[k+1];
std::vector<ValueType>& resVectorNew = calcLower ? unifVectors.resLowerNew : unifVectors.wUpper[k];
@@ -188,29 +201,6 @@ namespace storm {
// Expand the solution for the probabilistic states to all states.
storm::utility::vector::setVectorValues(resVectorNew, ~markovianStates, x);
}
-
- template <typename ValueType>
- void calculateResUpper(Environment const& env, std::vector<std::vector<ValueType>> const& relativeReachability, OptimizationDirection dir, uint64_t k, uint64_t state, ValueType lambda, uint64_t numberOfProbabilisticStates, UnifPlusVectors<ValueType>& unifVectors, storm::storage::SparseMatrix<ValueType> const& fullTransitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> const& solver, storm::utility::numerical::FoxGlynnResult<ValueType> const & poisson, bool cycleFree) {
-
- // Avoiding multiple computation of the same value.
- if (unifVectors.resUpperNew[state] != -1) {
- STORM_LOG_ASSERT(false, "Result was already calculated.");
- return;
- }
- uint64_t N = unifVectors.wUpper.size() - 1;
-
- ValueType res = storm::utility::zero<ValueType>();
- for (uint64_t i = k; i < N; ++i) {
- if (unifVectors.wUpper[N-1-(i-k)][state] == -1) {
- 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.resUpperNew[state] = res;
- }
template <typename ValueType>
void eliminateProbabilisticSelfLoops(storm::storage::SparseMatrix<ValueType>& transitionMatrix, storm::storage::BitVector const& markovianStates) {
@@ -259,7 +249,7 @@ namespace storm {
bool cycleFree = sccDecomposition.empty();
// Vectors to store computed vectors.
- UnifPlusVectors<ValueType> unifVectors(0, 0);
+ UnifPlusVectors<ValueType> 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.
@@ -383,22 +373,39 @@ namespace storm {
}
// (4) Define vectors/matrices.
- unifVectors = UnifPlusVectors<ValueType>(N+1, numberOfStates);
+ // Initialize result vectors and already insert zeros for iteration N
+ unifVectors = UnifPlusVectors<ValueType>(N, numberOfStates);
// (5) Compute vectors and maxNorm.
- for (int64_t k = N; k >= 0; --k) {
- for (uint64_t i = 0; i < numberOfStates; ++i) {
- calculateUnifPlusVector(env, k, i, true, lambda, numberOfProbabilisticChoices, relativeReachabilities, dir, unifVectors, fullTransitionMatrix, markovianAndGoalStates, psiStates, solver, foxGlynnResult, cycleFree);
- 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);
+ // Iteration k = N was already performed by initializing with zeros.
+
+ // Iterations k < N
+ 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
+ // resUpperNew was already initialized with zeros
+ uint64_t left = std::max(foxGlynnResult.left, (uint64_t)(k));
+ uint64_t right = std::min(foxGlynnResult.right, N-1);
+ for (uint64_t state = 0; state < numberOfStates; ++state) {
+ for (uint64_t i = left; i <= right; ++i) {
+ STORM_LOG_ASSERT(unifVectors.wUpper[N-1-(i-k)][state] != -1, "wUpper was not computed before.");
+ unifVectors.resUpperNew[state] += foxGlynnResult.weights[i - foxGlynnResult.left] * unifVectors.wUpper[N-1-(i-k)][state];
+ }
}
- unifVectors.prepareNewIteration();
}
// Only iterate over result vector, as the results can only get more precise.
maxNorm = storm::utility::zero<ValueType>();
for (uint64_t i = 0; i < numberOfStates; i++){
- ValueType diff = storm::utility::abs(unifVectors.resUpperOld[i] - unifVectors.resLowerOld[i]);
+ ValueType diff = storm::utility::abs(unifVectors.resUpperNew[i] - unifVectors.resLowerNew[i]);
maxNorm = std::max(maxNorm, diff);
}
@@ -408,7 +415,7 @@ namespace storm {
} while (maxNorm > epsilon * (1 - kappa));
- return unifVectors.resLowerOld;
+ return unifVectors.resLowerNew;
}
template <typename ValueType>
From 4c5b041340bb2c5fdb1abf5351fc5c6431c6c262 Mon Sep 17 00:00:00 2001
From: Matthias Volk <matthias.volk@cs.rwth-aachen.de>
Date: Thu, 6 Dec 2018 20:12:19 +0100
Subject: [PATCH 07/12] Debug output
---
.../csl/helper/SparseMarkovAutomatonCslHelper.cpp | 5 +++--
1 file changed, 3 insertions(+), 2 deletions(-)
diff --git a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
index 3af9714b8..1671f6245 100644
--- a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
+++ b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
@@ -90,6 +90,7 @@ namespace storm {
// First case, k==N, independent from kind of state.
if (k == N) {
+ STORM_LOG_ASSERT(false, "Result for k=N was already calculated.");
resVectorNew[state] = storm::utility::zero<ValueType>();
return;
}
@@ -290,7 +291,7 @@ namespace storm {
for (ValueType const& rate : exitRateVector) {
lambda = std::max(rate, lambda);
}
- STORM_LOG_TRACE("Initial lambda is " << lambda << ".");
+ STORM_LOG_DEBUG("Initial lambda is " << lambda << ".");
// Compute the relative reachability vectors and create solver for models with SCCs.
std::vector<std::vector<ValueType>> relativeReachabilities(transitionMatrix.getRowCount());
@@ -411,7 +412,7 @@ namespace storm {
// (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 << ".");
} while (maxNorm > epsilon * (1 - kappa));
From d054f3c64a5909100c5d6f2cf6b27b98cd5a7553 Mon Sep 17 00:00:00 2001
From: Matthias Volk <matthias.volk@cs.rwth-aachen.de>
Date: Fri, 7 Dec 2018 00:06:58 +0100
Subject: [PATCH 08/12] Result for upper bounds needs only be calculated for
k=0
---
.../helper/SparseMarkovAutomatonCslHelper.cpp | 18 ++++++++++--------
1 file changed, 10 insertions(+), 8 deletions(-)
diff --git a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
index 1671f6245..0b7d03b95 100644
--- a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
+++ b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
@@ -391,14 +391,16 @@ namespace storm {
calculateUnifPlusVector(env, k, state, false, lambda, numberOfProbabilisticChoices, relativeReachabilities, dir, unifVectors, fullTransitionMatrix, markovianAndGoalStates, psiStates, solver, foxGlynnResult, cycleFree);
}
- // Calculate result for upper bound
- // resUpperNew was already initialized with zeros
- uint64_t left = std::max(foxGlynnResult.left, (uint64_t)(k));
- uint64_t right = std::min(foxGlynnResult.right, N-1);
- for (uint64_t state = 0; state < numberOfStates; ++state) {
- for (uint64_t i = left; i <= right; ++i) {
- STORM_LOG_ASSERT(unifVectors.wUpper[N-1-(i-k)][state] != -1, "wUpper was not computed before.");
- unifVectors.resUpperNew[state] += foxGlynnResult.weights[i - foxGlynnResult.left] * unifVectors.wUpper[N-1-(i-k)][state];
+ if (k == 0) {
+ // Calculate result for upper bound
+ // resUpperNew was already initialized with zeros
+ uint64_t left = std::max(foxGlynnResult.left, (uint64_t)(k));
+ uint64_t right = std::min(foxGlynnResult.right, N-1);
+ for (uint64_t state = 0; state < numberOfStates; ++state) {
+ for (uint64_t i = left; i <= right; ++i) {
+ STORM_LOG_ASSERT(unifVectors.wUpper[N-1-(i-k)][state] != -1, "wUpper was not computed before.");
+ unifVectors.resUpperNew[state] += foxGlynnResult.weights[i - foxGlynnResult.left] * unifVectors.wUpper[N-1-(i-k)][state];
+ }
}
}
}
From 01a23856a8e3c34ffdc6ac9d5b4b3f409ea77068 Mon Sep 17 00:00:00 2001
From: Matthias Volk <matthias.volk@cs.rwth-aachen.de>
Date: Fri, 7 Dec 2018 00:26:33 +0100
Subject: [PATCH 09/12] Drastically decreased memory consumption of Unif+
---
.../helper/SparseMarkovAutomatonCslHelper.cpp | 42 ++++++++-----------
1 file changed, 17 insertions(+), 25 deletions(-)
diff --git a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
index 0b7d03b95..427eec25f 100644
--- a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
+++ b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
@@ -47,11 +47,8 @@ namespace storm {
/**
* Initialize results vectors. vLowerOld, vUpperOld and wUpper[k=N] are initialized with zeros.
*/
- UnifPlusVectors(uint64_t steps, uint64_t noStates) : numberOfStates(noStates), resLowerOld(numberOfStates, storm::utility::zero<ValueType>()), resLowerNew(numberOfStates, -1), resUpperOld(numberOfStates, storm::utility::zero<ValueType>()), resUpperNew(numberOfStates, storm::utility::zero<ValueType>()) {
- // For wUpper we have to keep track of all previous results
- wUpper = std::vector<std::vector<ValueType>>(steps+1, std::vector<ValueType>(numberOfStates, -1));
- // Initialize entries for step N with zeros
- std::fill(wUpper[steps].begin(), wUpper[steps].end(), storm::utility::zero<ValueType>());
+ UnifPlusVectors(uint64_t steps, uint64_t noStates) : numberOfStates(noStates), steps(steps), resLowerOld(numberOfStates, storm::utility::zero<ValueType>()), resLowerNew(numberOfStates, -1), resUpper(numberOfStates, storm::utility::zero<ValueType>()), wUpperOld(numberOfStates, storm::utility::zero<ValueType>()), wUpperNew(numberOfStates, -1) {
+ // Intentionally left empty
}
/**
@@ -60,23 +57,24 @@ namespace storm {
void prepareNewIteration() {
resLowerOld.swap(resLowerNew);
std::fill(resLowerNew.begin(), resLowerNew.end(), -1);
- resUpperOld.swap(resUpperNew);
- std::fill(resUpperNew.begin(), resUpperNew.end(), storm::utility::zero<ValueType>());
+ wUpperOld.swap(wUpperNew);
+ std::fill(wUpperNew.begin(), wUpperNew.end(), -1);
}
uint64_t numberOfStates;
+ uint64_t steps;
std::vector<ValueType> resLowerOld;
std::vector<ValueType> resLowerNew;
- std::vector<ValueType> resUpperOld;
- std::vector<ValueType> resUpperNew;
- std::vector<std::vector<ValueType>> wUpper;
+ std::vector<ValueType> resUpper;
+ std::vector<ValueType> wUpperOld;
+ std::vector<ValueType> wUpperNew;
};
template<typename ValueType>
void calculateUnifPlusVector(Environment const& env, uint64_t k, uint64_t state, bool calcLower, ValueType lambda, uint64_t numberOfProbabilisticChoices, std::vector<std::vector<ValueType>> const & relativeReachability, OptimizationDirection dir, UnifPlusVectors<ValueType>& unifVectors, storm::storage::SparseMatrix<ValueType> const& fullTransitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> const& solver, storm::utility::numerical::FoxGlynnResult<ValueType> const& poisson, bool cycleFree) {
// Set reference to acutal vector
- std::vector<ValueType>& resVectorOld = calcLower ? unifVectors.resLowerOld : unifVectors.wUpper[k+1];
- std::vector<ValueType>& resVectorNew = calcLower ? unifVectors.resLowerNew : unifVectors.wUpper[k];
+ std::vector<ValueType>& resVectorOld = calcLower ? unifVectors.resLowerOld : unifVectors.wUpperOld;
+ std::vector<ValueType>& resVectorNew = calcLower ? unifVectors.resLowerNew : unifVectors.wUpperNew;
if (resVectorNew[state] != -1) {
// Result already calculated.
@@ -84,7 +82,7 @@ namespace storm {
}
auto numberOfStates = fullTransitionMatrix.getRowGroupCount();
- uint64_t N = unifVectors.wUpper.size() - 1;
+ uint64_t N = unifVectors.steps;
auto const& rowGroupIndices = fullTransitionMatrix.getRowGroupIndices();
ValueType res;
@@ -389,26 +387,20 @@ namespace storm {
// 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);
- }
-
- if (k == 0) {
// Calculate result for upper bound
- // resUpperNew was already initialized with zeros
- uint64_t left = std::max(foxGlynnResult.left, (uint64_t)(k));
- uint64_t right = std::min(foxGlynnResult.right, N-1);
- for (uint64_t state = 0; state < numberOfStates; ++state) {
- for (uint64_t i = left; i <= right; ++i) {
- STORM_LOG_ASSERT(unifVectors.wUpper[N-1-(i-k)][state] != -1, "wUpper was not computed before.");
- unifVectors.resUpperNew[state] += foxGlynnResult.weights[i - foxGlynnResult.left] * unifVectors.wUpper[N-1-(i-k)][state];
- }
+ 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];
}
}
+
}
// Only iterate over result vector, as the results can only get more precise.
maxNorm = storm::utility::zero<ValueType>();
for (uint64_t i = 0; i < numberOfStates; i++){
- ValueType diff = storm::utility::abs(unifVectors.resUpperNew[i] - unifVectors.resLowerNew[i]);
+ ValueType diff = storm::utility::abs(unifVectors.resUpper[i] - unifVectors.resLowerNew[i]);
maxNorm = std::max(maxNorm, diff);
}
From 94a1d4710347cb2cbcfcbe06f969a2a117774cec Mon Sep 17 00:00:00 2001
From: Matthias Volk <matthias.volk@cs.rwth-aachen.de>
Date: Fri, 7 Dec 2018 11:32:36 +0100
Subject: [PATCH 10/12] Use IMCA method for bounded until with lower bound > 0
---
.../csl/helper/SparseMarkovAutomatonCslHelper.cpp | 7 ++++++-
1 file changed, 6 insertions(+), 1 deletion(-)
diff --git a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
index 427eec25f..32e95fd36 100644
--- a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
+++ b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
@@ -603,7 +603,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);
+ }
}
}
From 230ac204807619277381c366e5eae29bd33e647d Mon Sep 17 00:00:00 2001
From: Matthias Volk <matthias.volk@cs.rwth-aachen.de>
Date: Fri, 7 Dec 2018 17:02:40 +0100
Subject: [PATCH 11/12] Added progress measurements for Unif+ iterations and
steps
---
.../csl/helper/SparseMarkovAutomatonCslHelper.cpp | 10 ++++++++--
1 file changed, 8 insertions(+), 2 deletions(-)
diff --git a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
index 32e95fd36..ccc0cc324 100644
--- a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
+++ b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
@@ -327,6 +327,9 @@ namespace storm {
ValueType maxNorm = storm::utility::zero<ValueType>();
// 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.
do {
// (2) update parameter
@@ -379,6 +382,9 @@ namespace storm {
// 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();
@@ -394,7 +400,7 @@ namespace storm {
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.
@@ -407,7 +413,7 @@ namespace storm {
// (6) Double lambda.
lambda *= 2;
STORM_LOG_DEBUG("Increased lambda to " << lambda << ", max diff is " << maxNorm << ".");
-
+ progressIterations.updateProgress(++iteration);
} while (maxNorm > epsilon * (1 - kappa));
return unifVectors.resLowerNew;
From 99808240bff7e27d7b11586c26e1d6eaa837bc34 Mon Sep 17 00:00:00 2001
From: Matthias Volk <matthias.volk@cs.rwth-aachen.de>
Date: Sat, 8 Dec 2018 00:25:10 +0100
Subject: [PATCH 12/12] Updated changelog
---
CHANGELOG.md | 1 +
1 file changed, 1 insertion(+)
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