Sebastian Junges
5 years ago
39 changed files with 759 additions and 518 deletions
-
5CHANGELOG.md
-
5CMakeLists.txt
-
9doc/checklist_new_release.md
-
22resources/3rdparty/CMakeLists.txt
-
7resources/3rdparty/sylvan/src/storm_wrapper.cpp
-
1resources/3rdparty/sylvan/src/sylvan_obj_mtbdd_storm.hpp
-
6resources/3rdparty/sylvan/src/sylvan_obj_storm.cpp
-
7src/CMakeLists.txt
-
40src/storm-cli-utilities/model-handling.h
-
66src/storm-pars/modelchecker/instantiation/SparseDtmcInstantiationModelChecker.cpp
-
82src/storm-pars/modelchecker/instantiation/SparseMdpInstantiationModelChecker.cpp
-
4src/storm-pars/storage/ParameterRegion.cpp
-
27src/storm/builder/DdPrismModelBuilder.cpp
-
5src/storm/environment/solver/OviSolverEnvironment.cpp
-
2src/storm/environment/solver/OviSolverEnvironment.h
-
8src/storm/modelchecker/CheckTask.h
-
2src/storm/modelchecker/helper/infinitehorizon/SparseDeterministicInfiniteHorizonHelper.cpp
-
10src/storm/modelchecker/helper/infinitehorizon/SparseInfiniteHorizonHelper.cpp
-
2src/storm/settings/modules/CoreSettings.cpp
-
9src/storm/settings/modules/OviSolverSettings.cpp
-
3src/storm/settings/modules/OviSolverSettings.h
-
5src/storm/solver/AbstractEquationSolver.cpp
-
1src/storm/solver/AbstractEquationSolver.h
-
63src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp
-
2src/storm/solver/IterativeMinMaxLinearEquationSolver.h
-
57src/storm/solver/NativeLinearEquationSolver.cpp
-
2src/storm/solver/NativeLinearEquationSolver.h
-
13src/storm/solver/TopologicalLinearEquationSolver.cpp
-
11src/storm/solver/TopologicalMinMaxLinearEquationSolver.cpp
-
374src/storm/solver/helper/OptimisticValueIterationHelper.cpp
-
347src/storm/solver/helper/OptimisticValueIterationHelper.h
-
4src/storm/storage/dd/sylvan/InternalSylvanAdd.cpp
-
28src/storm/utility/AutomaticSettings.cpp
-
4src/storm/utility/AutomaticSettings.h
-
8src/storm/utility/Engine.cpp
-
2src/storm/utility/Engine.h
-
17src/test/storm/modelchecker/csl/LraCtmcCslModelCheckerTest.cpp
-
15src/test/storm/modelchecker/prctl/mdp/LraMdpPrctlModelCheckerTest.cpp
-
2version.cmake
@ -0,0 +1,374 @@ |
|||
#include "OptimisticValueIterationHelper.h"
|
|||
|
|||
#include "storm/utility/vector.h"
|
|||
#include "storm/utility/SignalHandler.h"
|
|||
#include "storm/environment/solver/OviSolverEnvironment.h"
|
|||
|
|||
#include "storm/exceptions/NotSupportedException.h"
|
|||
|
|||
#include "storm/utility/macros.h"
|
|||
|
|||
namespace storm { |
|||
|
|||
namespace solver { |
|||
namespace helper { |
|||
namespace oviinternal { |
|||
|
|||
template<typename ValueType> |
|||
ValueType updateIterationPrecision(storm::Environment const& env, ValueType const& diff) { |
|||
auto factor = storm::utility::convertNumber<ValueType>(env.solver().ovi().getPrecisionUpdateFactor()); |
|||
return factor * diff; |
|||
} |
|||
|
|||
template<typename ValueType> |
|||
void guessUpperBoundRelative(std::vector<ValueType> const& x, std::vector<ValueType> &target, ValueType const& relativeBoundGuessingScaler) { |
|||
storm::utility::vector::applyPointwise<ValueType, ValueType>(x, target, [&relativeBoundGuessingScaler] (ValueType const& argument) -> ValueType { return argument * relativeBoundGuessingScaler; }); |
|||
} |
|||
|
|||
template<typename ValueType> |
|||
void guessUpperBoundAbsolute(std::vector<ValueType> const& x, std::vector<ValueType> &target, ValueType const& precision) { |
|||
storm::utility::vector::applyPointwise<ValueType, ValueType>(x, target, [&precision] (ValueType const& argument) -> ValueType { return argument + precision; }); |
|||
} |
|||
|
|||
template <typename ValueType> |
|||
IterationHelper<ValueType>::IterationHelper(storm::storage::SparseMatrix<ValueType> const& matrix) { |
|||
STORM_LOG_THROW(static_cast<uint64_t>(std::numeric_limits<uint32_t>::max()) > matrix.getRowCount() + 1, storm::exceptions::NotSupportedException, "Matrix dimensions too large."); |
|||
STORM_LOG_THROW(static_cast<uint64_t>(std::numeric_limits<uint32_t>::max()) > matrix.getEntryCount(), storm::exceptions::NotSupportedException, "Matrix dimensions too large."); |
|||
matrixValues.reserve(matrix.getNonzeroEntryCount()); |
|||
matrixColumns.reserve(matrix.getColumnCount()); |
|||
rowIndications.reserve(matrix.getRowCount() + 1); |
|||
rowIndications.push_back(0); |
|||
for (IndexType r = 0; r < static_cast<IndexType>(matrix.getRowCount()); ++r) { |
|||
for (auto const& entry : matrix.getRow(r)) { |
|||
matrixValues.push_back(entry.getValue()); |
|||
matrixColumns.push_back(entry.getColumn()); |
|||
} |
|||
rowIndications.push_back(matrixValues.size()); |
|||
} |
|||
if (!matrix.hasTrivialRowGrouping()) { |
|||
rowGroupIndices = &matrix.getRowGroupIndices(); |
|||
} |
|||
} |
|||
|
|||
template <typename ValueType> |
|||
ValueType IterationHelper<ValueType>::singleIterationWithDiff(std::vector<ValueType>& x, std::vector<ValueType> const& b, bool computeRelativeDiff) { |
|||
return singleIterationWithDiffInternal<false, storm::solver::OptimizationDirection::Minimize>(x, b, computeRelativeDiff); |
|||
} |
|||
|
|||
template <typename ValueType> |
|||
ValueType IterationHelper<ValueType>::singleIterationWithDiff(storm::solver::OptimizationDirection const& dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, bool computeRelativeDiff) { |
|||
if (minimize(dir)) { |
|||
return singleIterationWithDiffInternal<true, storm::solver::OptimizationDirection::Minimize>(x, b, computeRelativeDiff); |
|||
} else { |
|||
return singleIterationWithDiffInternal<true, storm::solver::OptimizationDirection::Maximize>(x, b, computeRelativeDiff); |
|||
} |
|||
} |
|||
|
|||
template <typename ValueType> |
|||
template<bool HasRowGroups, storm::solver::OptimizationDirection Dir> |
|||
ValueType IterationHelper<ValueType>::singleIterationWithDiffInternal(std::vector<ValueType>& x, std::vector<ValueType> const& b, bool computeRelativeDiff) { |
|||
STORM_LOG_ASSERT(x.size() > 0, "Empty equation system not expected."); |
|||
ValueType diff = storm::utility::zero<ValueType>(); |
|||
|
|||
IndexType i = x.size(); |
|||
while (i > 0) { |
|||
--i; |
|||
ValueType newXi = HasRowGroups ? multiplyRowGroup<Dir>(i, b, x) : multiplyRow(i, b[i], x); |
|||
ValueType& oldXi = x[i]; |
|||
if (computeRelativeDiff) { |
|||
if (storm::utility::isZero(newXi)) { |
|||
if (storm::utility::isZero(oldXi)) { |
|||
// this operation has no effect:
|
|||
// diff = std::max(diff, storm::utility::zero<ValueType>());
|
|||
} else { |
|||
diff = std::max(diff, storm::utility::one<ValueType>()); |
|||
} |
|||
} else { |
|||
diff = std::max(diff, storm::utility::abs<ValueType>((newXi - oldXi) / newXi)); |
|||
} |
|||
} else { |
|||
diff = std::max(diff, storm::utility::abs<ValueType>(newXi - oldXi)); |
|||
} |
|||
oldXi = std::move(newXi); |
|||
} |
|||
return diff; |
|||
} |
|||
|
|||
template <typename ValueType> |
|||
uint64_t IterationHelper<ValueType>::repeatedIterate(storm::solver::OptimizationDirection const& dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, ValueType precision, bool relative) { |
|||
if (minimize(dir)) { |
|||
return repeatedIterateInternal<true, storm::solver::OptimizationDirection::Minimize>(x, b, precision, relative); |
|||
} else { |
|||
return repeatedIterateInternal<true, storm::solver::OptimizationDirection::Maximize>(x, b, precision, relative); |
|||
} |
|||
} |
|||
|
|||
template <typename ValueType> |
|||
uint64_t IterationHelper<ValueType>::repeatedIterate(std::vector<ValueType>& x, const std::vector<ValueType>& b, ValueType precision, bool relative) { |
|||
return repeatedIterateInternal<false, storm::solver::OptimizationDirection::Minimize>(x, b, precision, relative); |
|||
} |
|||
|
|||
template <typename ValueType> |
|||
template<bool HasRowGroups, storm::solver::OptimizationDirection Dir> |
|||
uint64_t IterationHelper<ValueType>::repeatedIterateInternal(std::vector<ValueType>& x, std::vector<ValueType> const& b, ValueType precision, bool relative) { |
|||
// Do a backwards gauss-seidel style iteration
|
|||
bool convergence = true; |
|||
IndexType i = x.size(); |
|||
while (i > 0) { |
|||
--i; |
|||
ValueType newXi = HasRowGroups ? multiplyRowGroup<Dir>(i, b, x) : multiplyRow(i, b[i], x); |
|||
ValueType& oldXi = x[i]; |
|||
// Check if we converged
|
|||
if (relative) { |
|||
if (storm::utility::isZero(oldXi)) { |
|||
if (!storm::utility::isZero(newXi)) { |
|||
convergence = false; |
|||
break; |
|||
} |
|||
} else if (storm::utility::abs<ValueType>((newXi - oldXi) / oldXi) > precision) { |
|||
convergence = false; |
|||
break; |
|||
} |
|||
} else { |
|||
if (storm::utility::abs<ValueType>((newXi - oldXi)) > precision) { |
|||
convergence = false; |
|||
break; |
|||
} |
|||
} |
|||
} |
|||
if (!convergence) { |
|||
// we now know that we did not converge. We still need to set the remaining values
|
|||
while (i > 0) { |
|||
--i; |
|||
x[i] = HasRowGroups ? multiplyRowGroup<Dir>(i, b, x) : multiplyRow(i, b[i], x); |
|||
} |
|||
} |
|||
return convergence; |
|||
} |
|||
|
|||
template <typename ValueType> |
|||
typename IterationHelper<ValueType>::IterateResult IterationHelper<ValueType>::iterateUpper(storm::solver::OptimizationDirection const& dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, bool takeMinOfOldAndNew) { |
|||
if (minimize(dir)) { |
|||
return iterateUpperInternal<true, storm::solver::OptimizationDirection::Minimize>(x, b, takeMinOfOldAndNew); |
|||
} else { |
|||
return iterateUpperInternal<true, storm::solver::OptimizationDirection::Maximize>(x, b, takeMinOfOldAndNew); |
|||
} |
|||
} |
|||
|
|||
template <typename ValueType> |
|||
typename IterationHelper<ValueType>::IterateResult IterationHelper<ValueType>::iterateUpper(std::vector<ValueType>& x, std::vector<ValueType> const& b, bool takeMinOfOldAndNew) { |
|||
return iterateUpperInternal<false, storm::solver::OptimizationDirection::Minimize>(x, b, takeMinOfOldAndNew); |
|||
} |
|||
|
|||
template <typename ValueType> |
|||
template<bool HasRowGroups, storm::solver::OptimizationDirection Dir> |
|||
typename IterationHelper<ValueType>::IterateResult IterationHelper<ValueType>::iterateUpperInternal(std::vector<ValueType>& x, std::vector<ValueType> const& b, bool takeMinOfOldAndNew) { |
|||
// For each row compare the new upper bound candidate with the old one
|
|||
bool newUpperBoundAlwaysHigherEqual = true; |
|||
bool newUpperBoundAlwaysLowerEqual = true; |
|||
// Do a backwards gauss-seidel style iteration
|
|||
for (IndexType i = x.size(); i > 0;) { |
|||
--i; |
|||
ValueType newXi = HasRowGroups ? multiplyRowGroup<Dir>(i, b, x) : multiplyRow(i, b[i], x); |
|||
ValueType& oldXi = x[i]; |
|||
if (newXi > oldXi) { |
|||
newUpperBoundAlwaysLowerEqual = false; |
|||
if (!takeMinOfOldAndNew) { |
|||
oldXi = newXi; |
|||
} |
|||
} else if (newXi != oldXi) { |
|||
assert(newXi < oldXi); |
|||
newUpperBoundAlwaysHigherEqual = false; |
|||
oldXi = newXi; |
|||
} |
|||
} |
|||
// Return appropriate result
|
|||
if (newUpperBoundAlwaysLowerEqual) { |
|||
if (newUpperBoundAlwaysHigherEqual) { |
|||
return IterateResult::Equal; |
|||
} else { |
|||
return IterateResult::AlwaysLowerOrEqual; |
|||
} |
|||
} else { |
|||
if (newUpperBoundAlwaysHigherEqual) { |
|||
return IterateResult::AlwaysHigherOrEqual; |
|||
} else { |
|||
return IterateResult::Incomparable; |
|||
} |
|||
} |
|||
} |
|||
|
|||
template <typename ValueType> |
|||
ValueType IterationHelper<ValueType>::multiplyRow(IndexType const& rowIndex, ValueType const& bi, std::vector<ValueType> const& x) { |
|||
assert(rowIndex < rowIndications.size()); |
|||
ValueType xRes = bi; |
|||
|
|||
auto entryIt = matrixValues.begin() + rowIndications[rowIndex]; |
|||
auto entryItE = matrixValues.begin() + rowIndications[rowIndex + 1]; |
|||
auto colIt = matrixColumns.begin() + rowIndications[rowIndex]; |
|||
for (; entryIt != entryItE; ++entryIt, ++colIt) { |
|||
xRes += *entryIt * x[*colIt]; |
|||
} |
|||
return xRes; |
|||
} |
|||
|
|||
template <typename ValueType> |
|||
template<storm::solver::OptimizationDirection Dir> |
|||
ValueType IterationHelper<ValueType>::multiplyRowGroup(IndexType const& rowGroupIndex, std::vector<ValueType> const& b, std::vector<ValueType> const& x) { |
|||
STORM_LOG_ASSERT(rowGroupIndices != nullptr, "No row group indices available."); |
|||
auto row = (*rowGroupIndices)[rowGroupIndex]; |
|||
auto const& groupEnd = (*rowGroupIndices)[rowGroupIndex + 1]; |
|||
STORM_LOG_ASSERT(row < groupEnd, "Empty row group not expected."); |
|||
ValueType xRes = multiplyRow(row, b[row], x); |
|||
for (++row; row < groupEnd; ++row) { |
|||
ValueType xCur = multiplyRow(row, b[row], x); |
|||
xRes = minimize(Dir) ? std::min(xRes, xCur) : std::max(xRes, xCur); |
|||
} |
|||
return xRes; |
|||
} |
|||
} |
|||
|
|||
template<typename ValueType> |
|||
OptimisticValueIterationHelper<ValueType>::OptimisticValueIterationHelper(storm::storage::SparseMatrix<ValueType> const& matrix) : iterationHelper(matrix) { |
|||
// Intentionally left empty.
|
|||
} |
|||
|
|||
template<typename ValueType> |
|||
std::pair<SolverStatus, uint64_t> OptimisticValueIterationHelper<ValueType>::solveEquations(Environment const& env, std::vector<ValueType>* lowerX, std::vector<ValueType>* upperX, std::vector<ValueType> const& b, bool relative, ValueType precision, uint64_t maxOverallIterations, boost::optional<storm::solver::OptimizationDirection> dir, boost::optional<storm::storage::BitVector> const& relevantValues) { |
|||
STORM_LOG_ASSERT(lowerX->size() == upperX->size(), "Dimension missmatch."); |
|||
|
|||
// As we will shuffle pointers around, let's store the original positions here.
|
|||
std::vector<ValueType>* initLowerX = lowerX; |
|||
std::vector<ValueType>* initUpperX = upperX; |
|||
|
|||
uint64_t overallIterations = 0; |
|||
uint64_t lastValueIterationIterations = 0; |
|||
uint64_t currentVerificationIterations = 0; |
|||
|
|||
// Get some parameters for the algorithm
|
|||
// 2
|
|||
ValueType two = storm::utility::convertNumber<ValueType>(2.0); |
|||
// Use no termination guaranteed upper bound iteration method
|
|||
bool noTerminationGuarantee = env.solver().ovi().useNoTerminationGuaranteeMinimumMethod(); |
|||
// Desired max difference between upperX and lowerX
|
|||
ValueType doublePrecision = precision * two; |
|||
// Upper bound only iterations
|
|||
uint64_t upperBoundOnlyIterations = env.solver().ovi().getUpperBoundOnlyIterations(); |
|||
ValueType relativeBoundGuessingScaler = (storm::utility::one<ValueType>() + storm::utility::convertNumber<ValueType>(env.solver().ovi().getUpperBoundGuessingFactor()) * precision); |
|||
// Initial precision for the value iteration calls
|
|||
ValueType iterationPrecision = precision; |
|||
|
|||
SolverStatus status = SolverStatus::InProgress; |
|||
|
|||
while (status == SolverStatus::InProgress && overallIterations < maxOverallIterations) { |
|||
// Perform value iteration until convergence
|
|||
lastValueIterationIterations = dir ? iterationHelper.repeatedIterate(dir.get(), *lowerX, b, iterationPrecision, relative) : iterationHelper.repeatedIterate(*lowerX, b, iterationPrecision, relative); |
|||
overallIterations += lastValueIterationIterations; |
|||
|
|||
bool intervalIterationNeeded = false; |
|||
currentVerificationIterations = 0; |
|||
|
|||
if (relative) { |
|||
oviinternal::guessUpperBoundRelative(*lowerX, *upperX, relativeBoundGuessingScaler); |
|||
} else { |
|||
oviinternal::guessUpperBoundAbsolute(*lowerX, *upperX, precision); |
|||
} |
|||
|
|||
bool cancelGuess = false; |
|||
while (status == SolverStatus::InProgress && overallIterations < maxOverallIterations) { |
|||
++overallIterations; |
|||
++currentVerificationIterations; |
|||
// Perform value iteration stepwise for lower bound and guessed upper bound
|
|||
|
|||
// Upper bound iteration
|
|||
auto upperBoundIterResult = dir ? iterationHelper.iterateUpper(dir.get(), *upperX, b, !noTerminationGuarantee) : iterationHelper.iterateUpper(*upperX, b, !noTerminationGuarantee); |
|||
|
|||
if (upperBoundIterResult == oviinternal::IterationHelper<ValueType>::IterateResult::AlwaysHigherOrEqual) { |
|||
// All values moved up (and did not stay the same)
|
|||
// That means the guess for an upper bound is actually a lower bound
|
|||
auto diff = dir ? iterationHelper.singleIterationWithDiff(dir.get(), *upperX, b, relative) : iterationHelper.singleIterationWithDiff(*upperX, b, relative); |
|||
iterationPrecision = oviinternal::updateIterationPrecision(env, diff); |
|||
// We assume to have a single fixed point. We can thus safely set the new lower bound, to the wrongly guessed upper bound
|
|||
// Set lowerX to the upper bound candidate
|
|||
std::swap(lowerX, upperX); |
|||
break; |
|||
} else if (upperBoundIterResult == oviinternal::IterationHelper<ValueType>::IterateResult::AlwaysLowerOrEqual) { |
|||
// All values moved down (and stayed not the same)
|
|||
// This is a valid upper bound. We still need to check the precision.
|
|||
// We can safely use twice the requested precision, as we calculate the center of both vectors
|
|||
bool reachedPrecision; |
|||
if (relevantValues) { |
|||
reachedPrecision = storm::utility::vector::equalModuloPrecision(*lowerX, *upperX, relevantValues.get(), doublePrecision, relative); |
|||
} else { |
|||
reachedPrecision = storm::utility::vector::equalModuloPrecision(*lowerX, *upperX, doublePrecision, relative); |
|||
} |
|||
if (reachedPrecision) { |
|||
status = SolverStatus::Converged; |
|||
break; |
|||
} else { |
|||
// From now on, we keep updating both bounds
|
|||
intervalIterationNeeded = true; |
|||
} |
|||
// The following case below covers that both vectors (old and new) are equal.
|
|||
// Theoretically, this means that the precise fixpoint has been reached. However, numerical instabilities can be tricky and this detection might be incorrect (see the haddad-monmege model).
|
|||
// We therefore disable it. It is very unlikely that we guessed the right fixpoint anyway.
|
|||
//} else if (upperBoundIterResult == oviinternal::IterationHelper<ValueType>::IterateResult::Equal) {
|
|||
// In this case, the guessed upper bound is the precise fixpoint
|
|||
// status = SolverStatus::Converged;
|
|||
// break;
|
|||
} |
|||
|
|||
// Check whether we tried this guess for too long
|
|||
ValueType scaledIterationCount = storm::utility::convertNumber<ValueType>(currentVerificationIterations) * storm::utility::convertNumber<ValueType>(env.solver().ovi().getMaxVerificationIterationFactor()); |
|||
if (!intervalIterationNeeded && scaledIterationCount * iterationPrecision >= storm::utility::one<ValueType>()) { |
|||
cancelGuess = true; |
|||
// In this case we will make one more iteration on the lower bound (mainly to obtain a new iterationPrecision)
|
|||
} |
|||
|
|||
// Lower bound iteration (only if needed)
|
|||
if (cancelGuess || intervalIterationNeeded || currentVerificationIterations > upperBoundOnlyIterations) { |
|||
auto diff = dir ? iterationHelper.singleIterationWithDiff(dir.get(), *lowerX, b, relative) : iterationHelper.singleIterationWithDiff(*lowerX, b, relative); |
|||
|
|||
// Check whether the upper and lower bounds have crossed, i.e., the upper bound is smaller than the lower bound.
|
|||
bool valuesCrossed = false; |
|||
for (uint64_t i = 0; i < lowerX->size(); ++i) { |
|||
if ((*upperX)[i] < (*lowerX)[i]) { |
|||
valuesCrossed = true; |
|||
break; |
|||
} |
|||
} |
|||
|
|||
if (cancelGuess || valuesCrossed) { |
|||
// A new guess is needed.
|
|||
iterationPrecision = oviinternal::updateIterationPrecision(env, diff); |
|||
break; |
|||
} |
|||
} |
|||
} |
|||
if (storm::utility::resources::isTerminate()) { |
|||
status = SolverStatus::Aborted; |
|||
} |
|||
} // end while
|
|||
|
|||
// Swap the results into the output vectors (if necessary).
|
|||
assert(initLowerX != lowerX || (initLowerX == lowerX && initUpperX == upperX)); |
|||
if (initLowerX != lowerX) { |
|||
assert(initUpperX == lowerX); |
|||
assert(initLowerX == upperX); |
|||
lowerX->swap(*upperX); |
|||
} |
|||
|
|||
if (overallIterations > maxOverallIterations) { |
|||
status = SolverStatus::MaximalIterationsExceeded; |
|||
} |
|||
|
|||
return {status, overallIterations}; |
|||
} |
|||
|
|||
|
|||
template class OptimisticValueIterationHelper<double>; |
|||
template class OptimisticValueIterationHelper<storm::RationalNumber>; |
|||
} |
|||
} |
|||
} |
|||
|
|||
@ -1,4 +1,4 @@ |
|||
set(STORM_VERSION_MAJOR 1) |
|||
set(STORM_VERSION_MINOR 6) |
|||
set(STORM_VERSION_PATCH 0) |
|||
set(STORM_VERSION_PATCH 2) |
|||
|
|||
Write
Preview
Loading…
Cancel
Save
Reference in new issue