sp
/
tempest
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

First working version of time-bounded reachability for Markov automata. 

Former-commit-id: 6501cbfca4
tempestpy_adaptions
dehnert 11 years ago
parent
2ef0405104
commit
18711c01a3
6 changed files with 207 additions and 134 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 203
      src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.h
  2. 77
      src/solver/AbstractNondeterministicLinearEquationSolver.h
  3. 39
      src/solver/GmmxxNondeterministicLinearEquationSolver.h
  4. 2
      src/storage/SparseMatrix.cpp
  5. 12
      src/storm.cpp
  6. 8
      src/utility/solver.h

203
src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.h
View File

@ -22,7 +22,7 @@ namespace storm {
class SparseMarkovAutomatonCslModelChecker : public AbstractModelChecker<ValueType> {
public:
explicit SparseMarkovAutomatonCslModelChecker(storm::models::MarkovAutomaton<ValueType> const& model, storm::solver::AbstractNondeterministicLinearEquationSolver<ValueType>* linearEquationSolver) : AbstractModelChecker<ValueType>(model), minimumOperatorStack(), linearEquationSolver(linearEquationSolver) {
explicit SparseMarkovAutomatonCslModelChecker(storm::models::MarkovAutomaton<ValueType> const& model) : AbstractModelChecker<ValueType>(model), minimumOperatorStack() {
// Intentionally left empty.
}
@ -70,36 +70,20 @@ namespace storm {
return result;
}
std::vector<ValueType> checkTimeBoundedEventually(bool min, storm::storage::BitVector const& goalStates, uint_fast64_t lowerBound, uint_fast64_t upperBound) const {
if (!this->getModel().isClosed()) {
throw storm::exceptions::InvalidArgumentException() << "Unable to compute time-bounded reachability on non-closed Markov automaton.";
}
// (1) Compute the number of steps we need to take.
ValueType lambda = this->getModel().getMaximalExitRate();
ValueType delta = (2 * storm::settings::Settings::getInstance()->getOptionByLongName("precision").getArgument(0).getValueAsDouble()) / (upperBound * lambda * lambda);
LOG4CPLUS_INFO(logger, "Determined delta to be " << delta << ".");
// Get some data fields for convenient access.
std::vector<uint_fast64_t> const& nondeterministicChoiceIndices = this->getModel().getNondeterministicChoiceIndices();
std::vector<ValueType> const& exitRates = this->getModel().getExitRates();
typename storm::storage::SparseMatrix<ValueType> const& transitionMatrix = this->getModel().getTransitionMatrix();
// (2) Compute four sparse matrices:
// * a matrix A_MS with all (discretized) transitions from Markovian non-goal states to all Markovian non-goal states.
// * a matrix A_MStoPS with all (discretized) transitions from Markovian non-goal states to all probabilistic non-goal states.
// * a matrix A_PS with all (non-discretized) transitions from probabilistic non-goal states to other probabilistic non-goal states. This matrix has more rows than columns.
// * a matrix A_PStoMS with all (non-discretized) transitions from probabilistic non-goal states to all Markovian non-goal states. This matrix may have any shape.
storm::storage::BitVector const& markovianNonGoalStates = this->getModel().getMarkovianStates() & ~goalStates;
storm::storage::BitVector const& probabilisticNonGoalStates = ~this->getModel().getMarkovianStates() & ~goalStates;
typename storm::storage::SparseMatrix<ValueType> aMarkovian = this->getModel().getTransitionMatrix().getSubmatrix(markovianNonGoalStates, this->getModel().getNondeterministicChoiceIndices(), true);
typename storm::storage::SparseMatrix<ValueType> aMarkovianToProbabilistic = this->getModel().getTransitionMatrix().getSubmatrix(markovianNonGoalStates, probabilisticNonGoalStates, nondeterministicChoiceIndices);
typename storm::storage::SparseMatrix<ValueType> aProbabilistic = this->getModel().getTransitionMatrix().getSubmatrix(probabilisticNonGoalStates, this->getModel().getNondeterministicChoiceIndices());
typename storm::storage::SparseMatrix<ValueType> aProbabilisticToMarkovian = this->getModel().getTransitionMatrix().getSubmatrix(probabilisticNonGoalStates, markovianNonGoalStates, nondeterministicChoiceIndices);
static void computeBoundedReachability(bool min, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, std::vector<ValueType> const& exitRates, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& goalStates, storm::storage::BitVector const& markovianNonGoalStates, storm::storage::BitVector const& probabilisticNonGoalStates, std::vector<ValueType>& markovianNonGoalValues, std::vector<ValueType>& probabilisticNonGoalValues, ValueType delta, uint_fast64_t numberOfSteps) {
// Start by computing four sparse matrices:
// * a matrix aMarkovian with all (discretized) transitions from Markovian non-goal states to all Markovian non-goal states.
// * a matrix aMarkovianToProbabilistic with all (discretized) transitions from Markovian non-goal states to all probabilistic non-goal states.
// * a matrix aProbabilistic with all (non-discretized) transitions from probabilistic non-goal states to other probabilistic non-goal states.
// * a matrix aProbabilisticToMarkovian with all (non-discretized) transitions from probabilistic non-goal states to all Markovian non-goal states.
typename storm::storage::SparseMatrix<ValueType> aMarkovian = transitionMatrix.getSubmatrix(markovianNonGoalStates, nondeterministicChoiceIndices, true);
typename storm::storage::SparseMatrix<ValueType> aMarkovianToProbabilistic = transitionMatrix.getSubmatrix(markovianNonGoalStates, probabilisticNonGoalStates, nondeterministicChoiceIndices);
std::vector<uint_fast64_t> markovianNondeterministicChoiceIndices = computeNondeterministicChoiceIndicesForConstraint(nondeterministicChoiceIndices, markovianNonGoalStates);
typename storm::storage::SparseMatrix<ValueType> aProbabilistic = transitionMatrix.getSubmatrix(probabilisticNonGoalStates, nondeterministicChoiceIndices);
typename storm::storage::SparseMatrix<ValueType> aProbabilisticToMarkovian = transitionMatrix.getSubmatrix(probabilisticNonGoalStates, markovianNonGoalStates, nondeterministicChoiceIndices);
std::vector<uint_fast64_t> probabilisticNondeterministicChoiceIndices = computeNondeterministicChoiceIndicesForConstraint(nondeterministicChoiceIndices, probabilisticNonGoalStates);
// The matrices with transitions from Markovian states need to be digitized.
// Digitize aMarkovian. Based on whether the transition is a self-loop or not, we apply the two digitization rules.
uint_fast64_t rowIndex = 0;
for (auto state : markovianNonGoalStates) {
@ -125,53 +109,128 @@ namespace storm {
++rowIndex;
}
// (3) Initialize two vectors:
// * v_PS holds the probability values of probabilistic non-goal states.
// * v_MS holds the probability values of Markovian non-goal states.
std::vector<ValueType> vProbabilistic(probabilisticNonGoalStates.getNumberOfSetBits());
std::vector<ValueType> vMarkovian(markovianNonGoalStates.getNumberOfSetBits());
// Initialize the two vectors that hold the variable one-step probabilities to all target states for probabilistic and Markovian (non-goal) states.
std::vector<ValueType> bProbabilistic(aProbabilistic.getRowCount());
std::vector<ValueType> bMarkovian(markovianNonGoalStates.getNumberOfSetBits());
// (4) Compute the two fixed right-hand side vectors, one for Markovian states and one for the probabilistic ones.
std::vector<ValueType> bProbabilistic = this->getModel().getTransitionMatrix().getConstrainedRowSumVector(probabilisticNonGoalStates, nondeterministicChoiceIndices, ~this->getModel().getMarkovianStates() & goalStates, aProbabilistic.getRowCount());
storm::storage::BitVector probabilisticGoalStates = ~this->getModel().getMarkovianStates() & goalStates;
std::vector<ValueType> bMarkovian;
bMarkovian.reserve(markovianNonGoalStates.getNumberOfSetBits());
// Compute the two fixed right-hand side vectors, one for Markovian states and one for the probabilistic ones.
std::vector<ValueType> bProbabilisticFixed = transitionMatrix.getConstrainedRowSumVector(probabilisticNonGoalStates, nondeterministicChoiceIndices, goalStates, aProbabilistic.getRowCount());
std::vector<ValueType> bMarkovianFixed;
bMarkovianFixed.reserve(markovianNonGoalStates.getNumberOfSetBits());
for (auto state : markovianNonGoalStates) {
bMarkovian.push_back(storm::utility::constGetZero<ValueType>());
bMarkovianFixed.push_back(storm::utility::constGetZero<ValueType>());
typename storm::storage::SparseMatrix<ValueType>::Rows row = transitionMatrix.getRow(rowIndex);
typename storm::storage::SparseMatrix<ValueType>::Rows row = transitionMatrix.getRow(nondeterministicChoiceIndices[state]);
for (auto element : row) {
bMarkovian.back() += (1 - std::exp(-exitRates[state] * delta)) * element.value();
if (goalStates.get(element.column())) {
bMarkovianFixed.back() += (1 - std::exp(-exitRates[state] * delta)) * element.value();
}
}
}
std::cout << delta << std::endl;
std::cout << markovianNonGoalStates.toString() << std::endl;
std::cout << aMarkovian.toString() << std::endl;
std::cout << aMarkovianToProbabilistic.toString() << std::endl;
std::cout << probabilisticNonGoalStates.toString() << std::endl;
std::cout << aProbabilistic.toString() << std::endl;
std::cout << aProbabilisticToMarkovian.toString() << std::endl;
std::cout << bProbabilistic << std::endl;
std::cout << bMarkovian << std::endl;
std::unique_ptr<storm::solver::AbstractNondeterministicLinearEquationSolver<ValueType>> nondeterministiclinearEquationSolver = storm::utility::solver::getNondeterministicLinearEquationSolver<ValueType>();
// (3) Perform value iteration
// Perform the actual value iteration
// * loop until the step bound has been reached
// * in the loop:
// * perform value iteration using A_PSwG, v_PS and the vector b where b = (A * 1_G)|PS + A_PStoMS * v_MS
// and 1_G being the characteristic vector for all goal states.
// * perform one timed-step using v_MS := A_MSwG * v_MS + A_MStoPS * v_PS + (A * 1_G)|MS
//
// After the loop, perform one more step of the value iteration for PS states.
std::vector<ValueType> markovianNonGoalValuesSwap(markovianNonGoalValues);
std::vector<ValueType> multiplicationResultScratchMemory(aProbabilistic.getRowCount());
std::vector<ValueType> aProbabilisticScratchMemory(probabilisticNonGoalValues.size());
for (uint_fast64_t currentStep = 0; currentStep < numberOfSteps; ++currentStep) {
// Start by (re-)computing bProbabilistic = bProbabilisticFixed + aProbabilisticToMarkovian * vMarkovian.
aProbabilisticToMarkovian.multiplyWithVector(markovianNonGoalValues, bProbabilistic);
storm::utility::vector::addVectorsInPlace(bProbabilistic, bProbabilisticFixed);
// Now perform the inner value iteration for probabilistic states.
nondeterministiclinearEquationSolver->solveEquationSystem(min, aProbabilistic, probabilisticNonGoalValues, bProbabilistic, probabilisticNondeterministicChoiceIndices, &multiplicationResultScratchMemory, &aProbabilisticScratchMemory);
// (Re-)compute bMarkovian = bMarkovianFixed + aMarkovianToProbabilistic * vProbabilistic.
aMarkovianToProbabilistic.multiplyWithVector(probabilisticNonGoalValues, bMarkovian);
storm::utility::vector::addVectorsInPlace(bMarkovian, bMarkovianFixed);
aMarkovian.multiplyWithVector(markovianNonGoalValues, markovianNonGoalValuesSwap);
std::swap(markovianNonGoalValues, markovianNonGoalValuesSwap);
storm::utility::vector::addVectorsInPlace(markovianNonGoalValues, bMarkovian);
}
// After the loop, perform one more step of the value iteration for PS states.
aProbabilisticToMarkovian.multiplyWithVector(probabilisticNonGoalValues, bProbabilistic);
storm::utility::vector::addVectorsInPlace(bProbabilistic, bProbabilisticFixed);
nondeterministiclinearEquationSolver->solveEquationSystem(min, aProbabilistic, probabilisticNonGoalValues, bProbabilistic, probabilisticNondeterministicChoiceIndices);
}
std::vector<ValueType> checkTimeBoundedEventually(bool min, storm::storage::BitVector const& goalStates, ValueType lowerBound, ValueType upperBound) const {
// Check whether the automaton is closed.
if (!this->getModel().isClosed()) {
throw storm::exceptions::InvalidArgumentException() << "Unable to compute time-bounded reachability on non-closed Markov automaton.";
}
// (4) Finally, create the result vector out of 1_G and v_all.
// Check whether the given bounds were valid.
if (lowerBound < storm::utility::constGetZero<ValueType>() || upperBound < storm::utility::constGetZero<ValueType>() || upperBound < lowerBound) {
throw storm::exceptions::InvalidArgumentException() << "Illegal interval [";
}
// Get some data fields for convenient access.
std::vector<uint_fast64_t> const& nondeterministicChoiceIndices = this->getModel().getNondeterministicChoiceIndices();
typename storm::storage::SparseMatrix<ValueType> const& transitionMatrix = this->getModel().getTransitionMatrix();
std::vector<ValueType> const& exitRates = this->getModel().getExitRates();
storm::storage::BitVector const& markovianStates = this->getModel().getMarkovianStates();
// (1) Compute the accuracy we need to achieve the required error bound.
ValueType maxExitRate = this->getModel().getMaximalExitRate();
ValueType delta = (2 * storm::settings::Settings::getInstance()->getOptionByLongName("precision").getArgument(0).getValueAsDouble()) / (upperBound * maxExitRate * maxExitRate);
// (2) Compute the number of steps we need to make for the interval.
uint_fast64_t numberOfSteps = static_cast<uint_fast64_t>(std::ceil((upperBound - lowerBound) / delta));
std::cout << "Performing " << numberOfSteps << " iterations (delta=" << delta << ") for interval [" << lowerBound << ", " << upperBound << "]." << std::endl;
// (3) Compute the non-goal states and initialize two vectors
// * vProbabilistic holds the probability values of probabilistic non-goal states.
// * vMarkovian holds the probability values of Markovian non-goal states.
storm::storage::BitVector const& markovianNonGoalStates = markovianStates & ~goalStates;
storm::storage::BitVector const& probabilisticNonGoalStates = ~markovianStates & ~goalStates;
std::vector<ValueType> vProbabilistic(probabilisticNonGoalStates.getNumberOfSetBits());
std::vector<ValueType> vMarkovian(markovianNonGoalStates.getNumberOfSetBits());
// Return dummy vector for the time being.
return std::vector<ValueType>();
computeBoundedReachability(min, transitionMatrix, nondeterministicChoiceIndices, exitRates, markovianStates, goalStates, markovianNonGoalStates, probabilisticNonGoalStates, vMarkovian, vProbabilistic, delta, numberOfSteps);
// (4) If the lower bound of interval was non-zero, we need to take the current values as the starting values for a subsequent value iteration.
if (lowerBound != storm::utility::constGetZero<ValueType>()) {
std::vector<ValueType> vAllProbabilistic((~markovianStates).getNumberOfSetBits());
std::vector<ValueType> vAllMarkovian(markovianStates.getNumberOfSetBits());
// Create the starting value vectors for the next value iteration based on the results of the previous one.
storm::utility::vector::setVectorValues<ValueType>(vAllProbabilistic, ~markovianStates % goalStates, storm::utility::constGetOne<ValueType>());
storm::utility::vector::setVectorValues<ValueType>(vAllProbabilistic, ~markovianStates % ~goalStates, vProbabilistic);
std::cout << vAllProbabilistic << std::endl;
storm::utility::vector::setVectorValues<ValueType>(vAllMarkovian, markovianStates % goalStates, storm::utility::constGetOne<ValueType>());
storm::utility::vector::setVectorValues<ValueType>(vAllMarkovian, markovianStates % ~goalStates, vMarkovian);
std::cout << vAllMarkovian << std::endl;
// Compute the number of steps to reach the target interval.
numberOfSteps = static_cast<uint_fast64_t>(std::ceil(lowerBound / delta));
std::cout << "Performing " << numberOfSteps << " iterations (delta=" << delta << ") for interval [0, " << lowerBound << "]." << std::endl;
// FIXME: According to IMCA, the value of delta needs to be recomputed here, but using the equations from the source this doesn't make sense, because they always evaluate to the original value of delta.
// Compute the bounded reachability for interval [0, b-a].
computeBoundedReachability(min, transitionMatrix, nondeterministicChoiceIndices, exitRates, markovianStates, storm::storage::BitVector(this->getModel().getNumberOfStates()), markovianStates, ~markovianStates, vAllMarkovian, vAllProbabilistic, delta, numberOfSteps);
// Create the result vector out of vAllProbabilistic and vAllMarkovian and return it.
std::vector<ValueType> result(this->getModel().getNumberOfStates());
storm::utility::vector::setVectorValues(result, ~markovianStates, vAllProbabilistic);
storm::utility::vector::setVectorValues(result, markovianStates, vAllMarkovian);
return result;
} else {
// Create the result vector out of 1_G, vProbabilistic and vMarkovian and return it.
std::vector<ValueType> result(this->getModel().getNumberOfStates());
storm::utility::vector::setVectorValues<ValueType>(result, goalStates, storm::utility::constGetOne<ValueType>());
storm::utility::vector::setVectorValues(result, probabilisticNonGoalStates, vProbabilistic);
storm::utility::vector::setVectorValues(result, markovianNonGoalStates, vMarkovian);
return result;
}
}
std::vector<ValueType> checkLongRunAverage(bool min, storm::storage::BitVector const& goalStates) const {
@ -370,11 +429,8 @@ namespace storm {
// Finalize the matrix and solve the corresponding system of equations.
sspMatrix.finalize();
std::vector<ValueType> x(numberOfStatesNotInMecs + mecDecomposition.size());
if (linearEquationSolver != nullptr) {
this->linearEquationSolver->solveEquationSystem(min, sspMatrix, x, b, sspNondeterministicChoiceIndices);
} else {
throw storm::exceptions::InvalidStateException() << "No valid linear equation solver available.";
}
std::unique_ptr<storm::solver::AbstractNondeterministicLinearEquationSolver<ValueType>> nondeterministiclinearEquationSolver = storm::utility::solver::getNondeterministicLinearEquationSolver<ValueType>();
nondeterministiclinearEquationSolver->solveEquationSystem(min, sspMatrix, x, b, sspNondeterministicChoiceIndices);
// Prepare result vector.
std::vector<ValueType> result(this->getModel().getNumberOfStates());
@ -451,7 +507,7 @@ namespace storm {
// Then, we can eliminate the rows and columns for all states whose values are already known to be 0.
std::vector<ValueType> x(maybeStates.getNumberOfSetBits());
std::vector<uint_fast64_t> subNondeterministicChoiceIndices = this->computeNondeterministicChoiceIndicesForConstraint(maybeStates);
std::vector<uint_fast64_t> subNondeterministicChoiceIndices = computeNondeterministicChoiceIndicesForConstraint(this->getModel().getNondeterministicChoiceIndices(), maybeStates);
storm::storage::SparseMatrix<ValueType> submatrix = this->getModel().getTransitionMatrix().getSubmatrix(maybeStates, this->getModel().getNondeterministicChoiceIndices());
// Now prepare the mean sojourn times for all states so they can be used as the right-hand side of the equation system.
@ -465,11 +521,8 @@ namespace storm {
storm::utility::vector::selectVectorValuesRepeatedly(b, maybeStates, this->getModel().getNondeterministicChoiceIndices(), meanSojournTimes);
// Solve the corresponding system of equations.
if (linearEquationSolver != nullptr) {
this->linearEquationSolver->solveEquationSystem(min, submatrix, x, b, subNondeterministicChoiceIndices);
} else {
throw storm::exceptions::InvalidStateException() << "No valid linear equation solver available.";
}
std::unique_ptr<storm::solver::AbstractNondeterministicLinearEquationSolver<ValueType>> nondeterministiclinearEquationSolver = storm::utility::solver::getNondeterministicLinearEquationSolver<ValueType>();
nondeterministiclinearEquationSolver->solveEquationSystem(min, submatrix, x, b, subNondeterministicChoiceIndices);
// Create resulting vector.
std::vector<ValueType> result(this->getModel().getNumberOfStates());
@ -497,10 +550,7 @@ namespace storm {
* @param constraint A bit vector specifying which states are kept.
* @returns A vector of the nondeterministic choice indices of the subsystem induced by the given constraint.
*/
std::vector<uint_fast64_t> computeNondeterministicChoiceIndicesForConstraint(storm::storage::BitVector const& constraint) const {
// First, get a reference to the full nondeterministic choice indices.
std::vector<uint_fast64_t> const& nondeterministicChoiceIndices = this->getModel().getNondeterministicChoiceIndices();
static std::vector<uint_fast64_t> computeNondeterministicChoiceIndicesForConstraint(std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, storm::storage::BitVector const& constraint) {
// Reserve the known amount of slots for the resulting vector.
std::vector<uint_fast64_t> subNondeterministicChoiceIndices(constraint.getNumberOfSetBits() + 1);
uint_fast64_t currentRowCount = 0;
@ -522,9 +572,6 @@ namespace storm {
return subNondeterministicChoiceIndices;
}
// An object that is used for solving linear equations and performing matrix-vector multiplication.
std::unique_ptr<storm::solver::AbstractNondeterministicLinearEquationSolver<ValueType>> linearEquationSolver;
};
}
}

77
src/solver/AbstractNondeterministicLinearEquationSolver.h
View File

@ -13,9 +13,19 @@ namespace storm {
template<class Type>
class AbstractNondeterministicLinearEquationSolver {
public:
AbstractNondeterministicLinearEquationSolver() {
storm::settings::Settings* s = storm::settings::Settings::getInstance();
precision = s->getOptionByLongName("precision").getArgument(0).getValueAsDouble();
maxIterations = s->getOptionByLongName("maxIterations").getArgument(0).getValueAsUnsignedInteger();
relative = s->getOptionByLongName("relative").getArgument(0).getValueAsBoolean();
}
AbstractNondeterministicLinearEquationSolver(double precision, uint_fast64_t maxIterations, bool relative) : precision(precision), maxIterations(maxIterations), relative(relative) {
// Intentionally left empty.
}
virtual AbstractNondeterministicLinearEquationSolver<Type>* clone() const {
return new AbstractNondeterministicLinearEquationSolver<Type>();
return new AbstractNondeterministicLinearEquationSolver<Type>(this->precision, this->maxIterations, this->relative);
}
/*!
@ -69,38 +79,39 @@ namespace storm {
* as there are states in the MDP.
* @returns The solution vector x of the system of linear equations as the content of the parameter x.
*/
virtual void solveEquationSystem(bool minimize, storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type> const& b, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices) const {
LOG4CPLUS_INFO(logger, "Starting iterative solver.");
// Get the settings object to customize solving.
storm::settings::Settings* s = storm::settings::Settings::getInstance();
// Get relevant user-defined settings for solving the equations.
double precision = s->getOptionByLongName("precision").getArgument(0).getValueAsDouble();
uint_fast64_t maxIterations = s->getOptionByLongName("maxIterations").getArgument(0).getValueAsUnsignedInteger();
bool relative = s->getOptionByLongName("relative").getArgument(0).getValueAsBoolean();
virtual void solveEquationSystem(bool minimize, storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type> const& b, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, std::vector<Type>* multiplyResult = nullptr, std::vector<Type>* newX = nullptr) const {
// LOG4CPLUS_INFO(logger, "Starting iterative solver.");
// Set up the environment for the power method.
std::vector<Type> multiplyResult(A.getRowCount());
// auto startTime = std::chrono::high_resolution_clock::now();
// std::chrono::nanoseconds totalTime(0);
if (multiplyResult == nullptr) {
multiplyResult = new std::vector<Type>(A.getRowCount());
}
std::vector<Type>* currentX = &x;
std::vector<Type>* newX = new std::vector<Type>(x.size());
bool xMemoryProvided = true;
if (newX == nullptr) {
newX = new std::vector<Type>(x.size());
xMemoryProvided = false;
}
std::vector<Type>* swap = nullptr;
uint_fast64_t iterations = 0;
bool converged = false;
// Proceed with the iterations as long as the method did not converge or reach the
// user-specified maximum number of iterations.
while (!converged && iterations < maxIterations) {
// Compute x' = A*x + b.
A.multiplyWithVector(*currentX, multiplyResult);
storm::utility::vector::addVectorsInPlace(multiplyResult, b);
A.multiplyWithVector(*currentX, *multiplyResult);
storm::utility::vector::addVectorsInPlace(*multiplyResult, b);
// Reduce the vector x' by applying min/max for all non-deterministic choices as given by the topmost
// element of the min/max operator stack.
if (minimize) {
storm::utility::vector::reduceVectorMin(multiplyResult, *newX, nondeterministicChoiceIndices);
storm::utility::vector::reduceVectorMin(*multiplyResult, *newX, nondeterministicChoiceIndices);
} else {
storm::utility::vector::reduceVectorMax(multiplyResult, *newX, nondeterministicChoiceIndices);
storm::utility::vector::reduceVectorMax(*multiplyResult, *newX, nondeterministicChoiceIndices);
}
// Determine whether the method converged.
@ -112,24 +123,32 @@ namespace storm {
newX = swap;
++iterations;
}
// If we performed an odd number of iterations, we need to swap the x and currentX, because the newest result
// is currently stored in currentX, but x is the output vector.
if (iterations % 2 == 1) {
std::swap(x, *currentX);
delete currentX;
} else {
if (!xMemoryProvided) {
delete currentX;
}
} else if (!xMemoryProvided) {
delete newX;
}
// Check if the solver converged and issue a warning otherwise.
if (converged) {
LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iterations << " iterations.");
} else {
LOG4CPLUS_WARN(logger, "Iterative solver did not converge.");
}
// // Check if the solver converged and issue a warning otherwise.
// if (converged) {
// LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iterations << " iterations.");
// } else {
// LOG4CPLUS_WARN(logger, "Iterative solver did not converge.");
// }
}
protected:
double precision;
uint_fast64_t maxIterations;
bool relative;
};
} // namespace solver

39
src/solver/GmmxxNondeterministicLinearEquationSolver.h
View File

@ -47,42 +47,41 @@ namespace storm {
delete gmmxxMatrix;
}
virtual void solveEquationSystem(bool minimize, storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type> const& b, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices) const override {
// Get the settings object to customize solving.
storm::settings::Settings* s = storm::settings::Settings::getInstance();
// Get relevant user-defined settings for solving the equations.
double precision = s->getOptionByLongName("precision").getArgument(0).getValueAsDouble();
uint_fast64_t maxIterations = s->getOptionByLongName("maxIterations").getArgument(0).getValueAsUnsignedInteger();
bool relative = s->getOptionByLongName("relative").getArgument(0).getValueAsBoolean();
virtual void solveEquationSystem(bool minimize, storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type> const& b, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, std::vector<Type>* multiplyResult = nullptr, std::vector<Type>* newX = nullptr) const override {
// Transform the transition probability matrix to the gmm++ format to use its arithmetic.
gmm::csr_matrix<Type>* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(A);
// Set up the environment for the power method.
std::vector<Type> multiplyResult(A.getRowCount());
if (multiplyResult == nullptr) {
multiplyResult = new std::vector<Type>(A.getRowCount());
}
std::vector<Type>* currentX = &x;
std::vector<Type>* newX = new std::vector<Type>(x.size());
bool xMemoryProvided = true;
if (newX == nullptr) {
newX = new std::vector<Type>(x.size());
xMemoryProvided = false;
}
std::vector<Type>* swap = nullptr;
uint_fast64_t iterations = 0;
bool converged = false;
// Proceed with the iterations as long as the method did not converge or reach the user-specified maximum number
// of iterations.
while (!converged && iterations < maxIterations) {
while (!converged && iterations < this->maxIterations) {
// Compute x' = A*x + b.
gmm::mult(*gmmxxMatrix, *currentX, multiplyResult);
gmm::add(b, multiplyResult);
gmm::mult(*gmmxxMatrix, *currentX, *multiplyResult);
gmm::add(b, *multiplyResult);
// Reduce the vector x' by applying min/max for all non-deterministic choices.
if (minimize) {
storm::utility::vector::reduceVectorMin(multiplyResult, *newX, nondeterministicChoiceIndices);
storm::utility::vector::reduceVectorMin(*multiplyResult, *newX, nondeterministicChoiceIndices);
} else {
storm::utility::vector::reduceVectorMax(multiplyResult, *newX, nondeterministicChoiceIndices);
storm::utility::vector::reduceVectorMax(*multiplyResult, *newX, nondeterministicChoiceIndices);
}
// Determine whether the method converged.
converged = storm::utility::vector::equalModuloPrecision(*currentX, *newX, precision, relative);
converged = storm::utility::vector::equalModuloPrecision(*currentX, *newX, this->precision, this->relative);
// Update environment variables.
swap = currentX;
@ -95,8 +94,10 @@ namespace storm {
// is currently stored in currentX, but x is the output vector.
if (iterations % 2 == 1) {
std::swap(x, *currentX);
delete currentX;
} else {
if (!xMemoryProvided) {
delete currentX;
}
} else if (!xMemoryProvided) {
delete newX;
}
delete gmmxxMatrix;

2
src/storage/SparseMatrix.cpp
View File

@ -577,7 +577,7 @@ namespace storage {
LOG4CPLUS_DEBUG(logger, "Determined size of submatrix to be " << subRowCount << "x" << rowGroupConstraint.getNumberOfSetBits() << ".");
// Create and initialize resulting matrix.
SparseMatrix result(subRowCount, rowGroupConstraint.getNumberOfSetBits());
SparseMatrix result(subRowCount, columnConstraint.getNumberOfSetBits());
result.initialize(subNonZeroEntries);
// Create a temporary vector that stores for each index whose bit is set

12
src/storm.cpp
View File

@ -473,13 +473,11 @@ int main(const int argc, const char* argv[]) {
LOG4CPLUS_INFO(logger, "Model is a Markov automaton.");
std::shared_ptr<storm::models::MarkovAutomaton<double>> markovAutomaton = parser.getModel<storm::models::MarkovAutomaton<double>>();
markovAutomaton->close();
storm::modelchecker::csl::SparseMarkovAutomatonCslModelChecker<double> mc(*markovAutomaton, new storm::solver::AbstractNondeterministicLinearEquationSolver<double>());
std::cout << mc.checkExpectedTime(true, markovAutomaton->getLabeledStates("goal")) << std::endl;
std::cout << mc.checkExpectedTime(false, markovAutomaton->getLabeledStates("goal")) << std::endl;
std::cout << mc.checkLongRunAverage(true, markovAutomaton->getLabeledStates("goal")) << std::endl;
std::cout << mc.checkLongRunAverage(false, markovAutomaton->getLabeledStates("goal")) << std::endl;
std::cout << mc.checkTimeBoundedEventually(true, markovAutomaton->getLabeledStates("goal"), 0, 1) << std::endl;
std::cout << mc.checkTimeBoundedEventually(false, markovAutomaton->getLabeledStates("goal"), 0, 1) << std::endl;
storm::modelchecker::csl::SparseMarkovAutomatonCslModelChecker<double> mc(*markovAutomaton);
// std::cout << mc.checkExpectedTime(true, markovAutomaton->getLabeledStates("goal")) << std::endl;
// std::cout << mc.checkLongRunAverage(true, markovAutomaton->getLabeledStates("goal")) << std::endl;
// std::cout << mc.checkTimeBoundedEventually(true, markovAutomaton->getLabeledStates("goal"), 0, 1) << std::endl;
std::cout << mc.checkTimeBoundedEventually(true, markovAutomaton->getLabeledStates("goal"), 1, 2) << std::endl;
break;
}
case storm::models::Unknown:

8
src/utility/solver.h
View File

@ -1,6 +1,8 @@
#ifndef STORM_UTILITY_SOLVER_H_
#define STORM_UTILITY_SOLVER_H_
#include "src/solver/AbstractNondeterministicLinearEquationSolver.h"
#include "src/solver/GmmxxNondeterministicLinearEquationSolver.h"
#include "src/solver/GurobiLpSolver.h"
namespace storm {
@ -10,6 +12,12 @@ namespace storm {
std::unique_ptr<storm::solver::LpSolver> getLpSolver(std::string const& name) {
return std::unique_ptr<storm::solver::LpSolver>(new storm::solver::GurobiLpSolver(name));
}
template<typename ValueType>
std::unique_ptr<storm::solver::AbstractNondeterministicLinearEquationSolver<ValueType>> getNondeterministicLinearEquationSolver() {
return std::unique_ptr<storm::solver::AbstractNondeterministicLinearEquationSolver<ValueType>>(new storm::solver::AbstractNondeterministicLinearEquationSolver<ValueType>());
// return std::unique_ptr<storm::solver::AbstractNondeterministicLinearEquationSolver<ValueType>>(new storm::solver::GmmxxNondeterministicLinearEquationSolver<ValueType>());
}
}
}
}

Write Preview
Loading…
Cancel
Save