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

added override to LESolver

tempestpy_adaptions
Stefan Pranger 4 years ago
parent
9a8e8fec5f
commit
4a018c0ca7
1 changed files with 120 additions and 114 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. 234
      src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp

234
src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp
View File

@ -5,6 +5,7 @@
#include "storm/environment/solver/MinMaxSolverEnvironment.h"
#include "storm/environment/solver/OviSolverEnvironment.h"
#include "storm/environment/solver/MultiplierEnvironment.h"
#include "storm/utility/ConstantsComparator.h"
#include "storm/utility/KwekMehlhorn.h"
@ -19,27 +20,27 @@
namespace storm {
namespace solver {
template<typename ValueType>
IterativeMinMaxLinearEquationSolver<ValueType>::IterativeMinMaxLinearEquationSolver(std::unique_ptr<LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory) : linearEquationSolverFactory(std::move(linearEquationSolverFactory)) {
// Intentionally left empty
}
template<typename ValueType>
IterativeMinMaxLinearEquationSolver<ValueType>::IterativeMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, std::unique_ptr<LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory) : StandardMinMaxLinearEquationSolver<ValueType>(A), linearEquationSolverFactory(std::move(linearEquationSolverFactory)) {
// Intentionally left empty.
}
template<typename ValueType>
IterativeMinMaxLinearEquationSolver<ValueType>::IterativeMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType>&& A, std::unique_ptr<LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory) : StandardMinMaxLinearEquationSolver<ValueType>(std::move(A)), linearEquationSolverFactory(std::move(linearEquationSolverFactory)) {
// Intentionally left empty.
}
template<typename ValueType>
MinMaxMethod IterativeMinMaxLinearEquationSolver<ValueType>::getMethod(Environment const& env, bool isExactMode) const {
// Adjust the method if none was specified and we want exact or sound computations.
auto method = env.solver().minMax().getMethod();
if (isExactMode && method != MinMaxMethod::PolicyIteration && method != MinMaxMethod::RationalSearch && method != MinMaxMethod::ViToPi) {
if (env.solver().minMax().isMethodSetFromDefault()) {
STORM_LOG_INFO("Selecting 'Policy iteration' as the solution technique to guarantee exact results. If you want to override this, please explicitly specify a different method.");
@ -58,7 +59,7 @@ namespace storm {
STORM_LOG_THROW(method == MinMaxMethod::ValueIteration || method == MinMaxMethod::PolicyIteration || method == MinMaxMethod::RationalSearch || method == MinMaxMethod::SoundValueIteration || method == MinMaxMethod::IntervalIteration || method == MinMaxMethod::OptimisticValueIteration || method == MinMaxMethod::ViToPi, storm::exceptions::InvalidEnvironmentException, "This solver does not support the selected method.");
return method;
}
template<typename ValueType>
bool IterativeMinMaxLinearEquationSolver<ValueType>::internalSolveEquations(Environment const& env, OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
bool result = false;
@ -87,14 +88,14 @@ namespace storm {
default:
STORM_LOG_THROW(false, storm::exceptions::InvalidEnvironmentException, "This solver does not implement the selected solution method");
}
return result;
}
template<typename ValueType>
bool IterativeMinMaxLinearEquationSolver<ValueType>::solveInducedEquationSystem(Environment const& env, std::unique_ptr<LinearEquationSolver<ValueType>>& linearEquationSolver, std::vector<uint64_t> const& scheduler, std::vector<ValueType>& x, std::vector<ValueType>& subB, std::vector<ValueType> const& originalB) const {
assert(subB.size() == x.size());
// Resolve the nondeterminism according to the given scheduler.
bool convertToEquationSystem = this->linearEquationSolverFactory->getEquationProblemFormat(env) == LinearEquationSolverProblemFormat::EquationSystem;
storm::storage::SparseMatrix<ValueType> submatrix = this->A->selectRowsFromRowGroups(scheduler, convertToEquationSystem);
@ -102,7 +103,7 @@ namespace storm {
submatrix.convertToEquationSystem();
}
storm::utility::vector::selectVectorValues<ValueType>(subB, scheduler, this->A->getRowGroupIndices(), originalB);
// Check whether the linear equation solver is already initialized
if (!linearEquationSolver) {
// Initialize the equation solver
@ -116,14 +117,14 @@ namespace storm {
// Solve the equation system for the 'DTMC' and return true upon success
return linearEquationSolver->solveEquations(env, x, subB);
}
template<typename ValueType>
bool IterativeMinMaxLinearEquationSolver<ValueType>::solveEquationsPolicyIteration(Environment const& env, OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
// Create the initial scheduler.
std::vector<storm::storage::sparse::state_type> scheduler = this->hasInitialScheduler() ? this->getInitialScheduler() : std::vector<storm::storage::sparse::state_type>(this->A->getRowGroupCount());
return performPolicyIteration(env, dir, x, b, std::move(scheduler));
}
template<typename ValueType>
bool IterativeMinMaxLinearEquationSolver<ValueType>::performPolicyIteration(Environment const& env, OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<storm::storage::sparse::state_type>&& initialPolicy) const {
std::vector<storm::storage::sparse::state_type> scheduler = std::move(initialPolicy);
@ -162,7 +163,7 @@ namespace storm {
do {
// Solve the equation system for the 'DTMC'.
solveInducedEquationSystem(environmentOfSolver, solver, scheduler, x, subB, b);
// Go through the multiplication result and see whether we can improve any of the choices.
bool schedulerImproved = false;
for (uint_fast64_t group = 0; group < this->A->getRowGroupCount(); ++group) {
@ -172,14 +173,14 @@ namespace storm {
if (choice - this->A->getRowGroupIndices()[group] == currentChoice) {
continue;
}
// Create the value of the choice.
ValueType choiceValue = storm::utility::zero<ValueType>();
for (auto const& entry : this->A->getRow(choice)) {
choiceValue += entry.getValue() * x[entry.getColumn()];
}
choiceValue += b[choice];
// If the value is strictly better than the solution of the inner system, we need to improve the scheduler.
// TODO: If the underlying solver is not precise, this might run forever (i.e. when a state has two choices where the (exact) values are equal).
// only changing the scheduler if the values are not equal (modulo precision) would make this unsound.
@ -190,12 +191,12 @@ namespace storm {
}
}
}
// If the scheduler did not improve, we are done.
if (!schedulerImproved) {
status = SolverStatus::Converged;
}
// Update environment variables.
++iterations;
status = this->updateStatus(status, x, dir == storm::OptimizationDirection::Minimize ? SolverGuarantee::GreaterOrEqual : SolverGuarantee::LessOrEqual, iterations, env.solver().minMax().getMaximalNumberOfIterations());
@ -203,21 +204,21 @@ namespace storm {
// Potentially show progress.
this->showProgressIterative(iterations);
} while (status == SolverStatus::InProgress);
this->reportStatus(status, iterations);
// If requested, we store the scheduler for retrieval.
if (this->isTrackSchedulerSet()) {
this->schedulerChoices = std::move(scheduler);
}
if (!this->isCachingEnabled()) {
clearCache();
}
return status == SolverStatus::Converged || status == SolverStatus::TerminatedEarly;
}
template<typename ValueType>
bool IterativeMinMaxLinearEquationSolver<ValueType>::valueImproved(OptimizationDirection dir, ValueType const& value1, ValueType const& value2) const {
if (dir == OptimizationDirection::Minimize) {
@ -230,7 +231,7 @@ namespace storm {
template<typename ValueType>
MinMaxLinearEquationSolverRequirements IterativeMinMaxLinearEquationSolver<ValueType>::getRequirements(Environment const& env, boost::optional<storm::solver::OptimizationDirection> const& direction, bool const& hasInitialScheduler) const {
auto method = getMethod(env, storm::NumberTraits<ValueType>::IsExact || env.solver().isForceExact());
// Check whether a linear equation solver is needed and potentially start with its requirements
bool needsLinEqSolver = false;
needsLinEqSolver |= method == MinMaxMethod::PolicyIteration;
@ -259,7 +260,7 @@ namespace storm {
requirements.requireUniqueSolution();
}
requirements.requireLowerBounds();
} else if (method == MinMaxMethod::IntervalIteration) {
// Interval iteration requires a unique solution and lower+upper bounds
if (!this->hasUniqueSolution()) {
@ -296,18 +297,18 @@ namespace storm {
template<typename ValueType>
typename IterativeMinMaxLinearEquationSolver<ValueType>::ValueIterationResult IterativeMinMaxLinearEquationSolver<ValueType>::performValueIteration(Environment const& env, OptimizationDirection dir, std::vector<ValueType>*& currentX, std::vector<ValueType>*& newX, std::vector<ValueType> const& b, ValueType const& precision, bool relative, SolverGuarantee const& guarantee, uint64_t currentIterations, uint64_t maximalNumberOfIterations, storm::solver::MultiplicationStyle const& multiplicationStyle) const {
STORM_LOG_ASSERT(currentX != newX, "Vectors must not be aliased.");
// Get handle to multiplier.
storm::solver::Multiplier<ValueType> const& multiplier = *this->multiplierA;
// Allow aliased multiplications.
bool useGaussSeidelMultiplication = multiplicationStyle == storm::solver::MultiplicationStyle::GaussSeidel;
// Proceed with the iterations as long as the method did not converge or reach the maximum number of iterations.
uint64_t iterations = currentIterations;
SolverStatus status = SolverStatus::InProgress;
while (status == SolverStatus::InProgress) {
// Compute x' = min/max(A*x + b).
@ -318,12 +319,12 @@ namespace storm {
} else {
multiplier.multiplyAndReduce(env, dir, *currentX, &b, *newX);
}
// Determine whether the method converged.
if (storm::utility::vector::equalModuloPrecision<ValueType>(*currentX, *newX, precision, relative)) {
status = SolverStatus::Converged;
}
// Update environment variables.
std::swap(currentX, newX);
++iterations;
@ -332,7 +333,7 @@ namespace storm {
// Potentially show progress.
this->showProgressIterative(iterations);
}
return ValueIterationResult(iterations - currentIterations, status);
}
@ -367,7 +368,7 @@ namespace storm {
}
return true;
}
if (!auxiliaryRowGroupVector) {
auxiliaryRowGroupVector = std::make_unique<std::vector<ValueType>>(this->A->getRowGroupCount());
}
@ -376,14 +377,14 @@ namespace storm {
}
storm::solver::helper::OptimisticValueIterationHelper<ValueType> helper(*this->A);
// x has to start with a lower bound.
this->createLowerBoundsVector(x);
std::vector<ValueType>* lowerX = &x;
std::vector<ValueType>* upperX = auxiliaryRowGroupVector.get();
auto statusIters = helper.solveEquations(env, lowerX, upperX, b,
env.solver().minMax().getRelativeTerminationCriterion(),
storm::utility::convertNumber<ValueType>(env.solver().minMax().getPrecision()),
@ -413,14 +414,19 @@ namespace storm {
if (!this->multiplierA) {
this->multiplierA = storm::solver::MultiplierFactory<ValueType>().create(env, *this->A);
}
// TODO cleanup
if(env.solver().multiplier().getOptimizationDirectionOverride().is_initialized()) {
multiplierA->setOptimizationDirectionOverride(env.solver().multiplier().getOptimizationDirectionOverride().get());
}
if (!auxiliaryRowGroupVector) {
auxiliaryRowGroupVector = std::make_unique<std::vector<ValueType>>(this->A->getRowGroupCount());
}
// By default, we can not provide any guarantee
SolverGuarantee guarantee = SolverGuarantee::None;
if (this->hasInitialScheduler()) {
// Solve the equation system induced by the initial scheduler.
std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> linEqSolver;
@ -469,7 +475,7 @@ namespace storm {
std::vector<ValueType>* newX = auxiliaryRowGroupVector.get();
std::vector<ValueType>* currentX = &x;
this->startMeasureProgress();
ValueIterationResult result = performValueIteration(env, dir, currentX, newX, b, storm::utility::convertNumber<ValueType>(env.solver().minMax().getPrecision()), env.solver().minMax().getRelativeTerminationCriterion(), guarantee, 0, env.solver().minMax().getMaximalNumberOfIterations(), env.solver().minMax().getMultiplicationStyle());
@ -477,27 +483,27 @@ namespace storm {
if (currentX == auxiliaryRowGroupVector.get()) {
std::swap(x, *currentX);
}
this->reportStatus(result.status, result.iterations);
// If requested, we store the scheduler for retrieval.
if (this->isTrackSchedulerSet()) {
this->schedulerChoices = std::vector<uint_fast64_t>(this->A->getRowGroupCount());
this->multiplierA->multiplyAndReduce(env, dir, x, &b, *auxiliaryRowGroupVector.get(), &this->schedulerChoices.get());
}
if (!this->isCachingEnabled()) {
clearCache();
}
return result.status == SolverStatus::Converged || result.status == SolverStatus::TerminatedEarly;
}
template<typename ValueType>
void preserveOldRelevantValues(std::vector<ValueType> const& allValues, storm::storage::BitVector const& relevantValues, std::vector<ValueType>& oldValues) {
storm::utility::vector::selectVectorValues(oldValues, relevantValues, allValues);
}
/*!
* This version of value iteration is sound, because it approaches the solution from below and above. This
* technique is due to Haddad and Monmege (Interval iteration algorithm for MDPs and IMDPs, TCS 2017) and was
@ -511,28 +517,28 @@ namespace storm {
if (!this->multiplierA) {
this->multiplierA = storm::solver::MultiplierFactory<ValueType>().create(env, *this->A);
}
if (!auxiliaryRowGroupVector) {
auxiliaryRowGroupVector = std::make_unique<std::vector<ValueType>>(this->A->getRowGroupCount());
}
// Allow aliased multiplications.
bool useGaussSeidelMultiplication = env.solver().minMax().getMultiplicationStyle() == storm::solver::MultiplicationStyle::GaussSeidel;
std::vector<ValueType>* lowerX = &x;
this->createLowerBoundsVector(*lowerX);
this->createUpperBoundsVector(this->auxiliaryRowGroupVector, this->A->getRowGroupCount());
std::vector<ValueType>* upperX = this->auxiliaryRowGroupVector.get();
std::vector<ValueType>* tmp = nullptr;
if (!useGaussSeidelMultiplication) {
auxiliaryRowGroupVector2 = std::make_unique<std::vector<ValueType>>(lowerX->size());
tmp = auxiliaryRowGroupVector2.get();
}
// Proceed with the iterations as long as the method did not converge or reach the maximum number of iterations.
uint64_t iterations = 0;
SolverStatus status = SolverStatus::InProgress;
bool doConvergenceCheck = true;
bool useDiffs = this->hasRelevantValues() && !env.solver().minMax().isSymmetricUpdatesSet();
@ -636,7 +642,7 @@ namespace storm {
status = storm::utility::vector::equalModuloPrecision<ValueType>(*lowerX, *upperX, precision, relative) ? SolverStatus::Converged : status;
}
}
// Update environment variables.
++iterations;
doConvergenceCheck = !doConvergenceCheck;
@ -650,33 +656,33 @@ namespace storm {
// Potentially show progress.
this->showProgressIterative(iterations);
}
this->reportStatus(status, iterations);
// We take the means of the lower and upper bound so we guarantee the desired precision.
ValueType two = storm::utility::convertNumber<ValueType>(2.0);
storm::utility::vector::applyPointwise<ValueType, ValueType, ValueType>(*lowerX, *upperX, *lowerX, [&two] (ValueType const& a, ValueType const& b) -> ValueType { return (a + b) / two; });
// Since we shuffled the pointer around, we need to write the actual results to the input/output vector x.
if (&x == tmp) {
std::swap(x, *tmp);
} else if (&x == this->auxiliaryRowGroupVector.get()) {
std::swap(x, *this->auxiliaryRowGroupVector);
}
// If requested, we store the scheduler for retrieval.
if (this->isTrackSchedulerSet()) {
this->schedulerChoices = std::vector<uint_fast64_t>(this->A->getRowGroupCount());
this->multiplierA->multiplyAndReduce(env, dir, x, &b, *this->auxiliaryRowGroupVector, &this->schedulerChoices.get());
}
if (!this->isCachingEnabled()) {
clearCache();
}
return status == SolverStatus::Converged;
}
template<typename ValueType>
bool IterativeMinMaxLinearEquationSolver<ValueType>::solveEquationsSoundValueIteration(Environment const& env, OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
@ -690,7 +696,7 @@ namespace storm {
} else {
this->soundValueIterationHelper = std::make_unique<storm::solver::helper::SoundValueIterationHelper<ValueType>>(std::move(*this->soundValueIterationHelper), x, *this->auxiliaryRowGroupVector, env.solver().minMax().getRelativeTerminationCriterion(), storm::utility::convertNumber<ValueType>(env.solver().minMax().getPrecision()));
}
// Prepare initial bounds for the solution (if given)
if (this->hasLowerBound()) {
this->soundValueIterationHelper->setLowerBound(this->getLowerBound(true));
@ -698,16 +704,16 @@ namespace storm {
if (this->hasUpperBound()) {
this->soundValueIterationHelper->setUpperBound(this->getUpperBound(true));
}
storm::storage::BitVector const* relevantValuesPtr = nullptr;
if (this->hasRelevantValues()) {
relevantValuesPtr = &this->getRelevantValues();
}
SolverStatus status = SolverStatus::InProgress;
this->startMeasureProgress();
uint64_t iterations = 0;
while (status == SolverStatus::InProgress && iterations < env.solver().minMax().getMaximalNumberOfIterations()) {
++iterations;
this->soundValueIterationHelper->performIterationStep(dir, b);
@ -716,12 +722,12 @@ namespace storm {
} else {
status = this->updateStatus(status, this->hasCustomTerminationCondition() && this->soundValueIterationHelper->checkCustomTerminationCondition(this->getTerminationCondition()), iterations, env.solver().minMax().getMaximalNumberOfIterations());
}
// Potentially show progress.
this->showProgressIterative(iterations);
}
this->soundValueIterationHelper->setSolutionVector();
// If requested, we store the scheduler for retrieval.
if (this->isTrackSchedulerSet()) {
this->schedulerChoices = std::vector<uint_fast64_t>(this->A->getRowGroupCount());
@ -729,14 +735,14 @@ namespace storm {
}
this->reportStatus(status, iterations);
if (!this->isCachingEnabled()) {
clearCache();
}
return status == SolverStatus::Converged;
}
template<typename ValueType>
bool IterativeMinMaxLinearEquationSolver<ValueType>::solveEquationsViToPi(Environment const& env, OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
// First create an (inprecise) vi solver to get a good initial strategy for the (potentially precise) policy iteration solver.
@ -760,11 +766,11 @@ namespace storm {
STORM_LOG_INFO("Found initial policy using Value Iteration. Starting Policy iteration now.");
return performPolicyIteration(env, dir, x, b, std::move(initialSched));
}
template<typename ValueType>
bool IterativeMinMaxLinearEquationSolver<ValueType>::isSolution(storm::OptimizationDirection dir, storm::storage::SparseMatrix<ValueType> const& matrix, std::vector<ValueType> const& values, std::vector<ValueType> const& b) {
storm::utility::ConstantsComparator<ValueType> comparator;
auto valueIt = values.begin();
auto bIt = b.begin();
for (uint64_t group = 0; group < matrix.getRowGroupCount(); ++group, ++valueIt) {
@ -778,18 +784,18 @@ namespace storm {
for (auto endRow = matrix.getRowGroupIndices()[group + 1]; row < endRow; ++row, ++bIt) {
ValueType newValue = *bIt;
newValue += matrix.multiplyRowWithVector(row, values);
if ((dir == storm::OptimizationDirection::Minimize && newValue < groupValue) || (dir == storm::OptimizationDirection::Maximize && newValue > groupValue)) {
groupValue = newValue;
}
}
// If the value does not match the one in the values vector, the given vector is not a solution.
if (!comparator.isEqual(groupValue, *valueIt)) {
return false;
}
}
// Checked all values at this point.
return true;
}
@ -797,7 +803,7 @@ namespace storm {
template<typename ValueType>
template<typename RationalType, typename ImpreciseType>
bool IterativeMinMaxLinearEquationSolver<ValueType>::sharpen(storm::OptimizationDirection dir, uint64_t precision, storm::storage::SparseMatrix<RationalType> const& A, std::vector<ImpreciseType> const& x, std::vector<RationalType> const& b, std::vector<RationalType>& tmp) {
for (uint64_t p = 0; p <= precision; ++p) {
storm::utility::kwek_mehlhorn::sharpen(p, x, tmp);
@ -817,18 +823,18 @@ namespace storm {
storm::storage::SparseMatrix<storm::RationalNumber> rationalA = this->A->template toValueType<storm::RationalNumber>();
std::vector<storm::RationalNumber> rationalX(x.size());
std::vector<storm::RationalNumber> rationalB = storm::utility::vector::convertNumericVector<storm::RationalNumber>(b);
if (!this->multiplierA) {
this->multiplierA = storm::solver::MultiplierFactory<ValueType>().create(env, *this->A);
}
if (!auxiliaryRowGroupVector) {
auxiliaryRowGroupVector = std::make_unique<std::vector<ValueType>>(this->A->getRowGroupCount());
}
// Forward the call to the core rational search routine.
bool converged = solveEquationsRationalSearchHelper<storm::RationalNumber, ImpreciseType>(env, dir, *this, rationalA, rationalX, rationalB, *this->A, x, b, *auxiliaryRowGroupVector);
// Translate back rational result to imprecise result.
auto targetIt = x.begin();
for (auto it = rationalX.begin(), ite = rationalX.end(); it != ite; ++it, ++targetIt) {
@ -838,30 +844,30 @@ namespace storm {
if (!this->isCachingEnabled()) {
this->clearCache();
}
return converged;
}
template<typename ValueType>
template<typename ImpreciseType>
typename std::enable_if<std::is_same<ValueType, ImpreciseType>::value && NumberTraits<ValueType>::IsExact, bool>::type IterativeMinMaxLinearEquationSolver<ValueType>::solveEquationsRationalSearchHelper(Environment const& env, OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
// Version for when the overall value type is exact and the same type is to be used for the imprecise part.
if (!this->multiplierA) {
this->multiplierA = storm::solver::MultiplierFactory<ValueType>().create(env, *this->A);
}
if (!auxiliaryRowGroupVector) {
auxiliaryRowGroupVector = std::make_unique<std::vector<ValueType>>(this->A->getRowGroupCount());
}
// Forward the call to the core rational search routine.
bool converged = solveEquationsRationalSearchHelper<ValueType, ImpreciseType>(env, dir, *this, *this->A, x, b, *this->A, *auxiliaryRowGroupVector, b, x);
if (!this->isCachingEnabled()) {
this->clearCache();
}
return converged;
}
@ -870,10 +876,10 @@ namespace storm {
typename std::enable_if<!std::is_same<ValueType, ImpreciseType>::value, bool>::type IterativeMinMaxLinearEquationSolver<ValueType>::solveEquationsRationalSearchHelper(Environment const& env, OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
// Version for when the overall value type is exact and the imprecise one is not. We first try to solve the
// problem using the imprecise data type and fall back to the exact type as needed.
// Translate A to its imprecise version.
storm::storage::SparseMatrix<ImpreciseType> impreciseA = this->A->template toValueType<ImpreciseType>();
// Translate x to its imprecise version.
std::vector<ImpreciseType> impreciseX(x.size());
{
@ -884,23 +890,23 @@ namespace storm {
*targetIt = storm::utility::convertNumber<ImpreciseType, ValueType>(*sourceIt);
}
}
// Create temporary storage for an imprecise x.
std::vector<ImpreciseType> impreciseTmpX(x.size());
// Translate b to its imprecise version.
std::vector<ImpreciseType> impreciseB(b.size());
auto targetIt = impreciseB.begin();
for (auto sourceIt = b.begin(); targetIt != impreciseB.end(); ++targetIt, ++sourceIt) {
*targetIt = storm::utility::convertNumber<ImpreciseType, ValueType>(*sourceIt);
}
// Create imprecise solver from the imprecise data.
IterativeMinMaxLinearEquationSolver<ImpreciseType> impreciseSolver(std::make_unique<storm::solver::GeneralLinearEquationSolverFactory<ImpreciseType>>());
impreciseSolver.setMatrix(impreciseA);
impreciseSolver.setCachingEnabled(true);
impreciseSolver.multiplierA = storm::solver::MultiplierFactory<ImpreciseType>().create(env, impreciseA);
bool converged = false;
try {
// Forward the call to the core rational search routine.
@ -908,7 +914,7 @@ namespace storm {
impreciseSolver.clearCache();
} catch (storm::exceptions::PrecisionExceededException const& e) {
STORM_LOG_WARN("Precision of value type was exceeded, trying to recover by switching to rational arithmetic.");
if (!auxiliaryRowGroupVector) {
auxiliaryRowGroupVector = std::make_unique<std::vector<ValueType>>(this->A->getRowGroupCount());
}
@ -918,7 +924,7 @@ namespace storm {
for (auto it = impreciseX.begin(), ite = impreciseX.end(); it != ite; ++it, ++targetIt) {
*targetIt = storm::utility::convertNumber<ValueType>(*it);
}
// Get rid of the superfluous data structures.
impreciseX = std::vector<ImpreciseType>();
impreciseTmpX = std::vector<ImpreciseType>();
@ -928,15 +934,15 @@ namespace storm {
if (!this->multiplierA) {
this->multiplierA = storm::solver::MultiplierFactory<ValueType>().create(env, *this->A);
}
// Forward the call to the core rational search routine, but now with our value type as the imprecise value type.
converged = solveEquationsRationalSearchHelper<ValueType, ValueType>(env, dir, *this, *this->A, x, b, *this->A, *auxiliaryRowGroupVector, b, x);
}
if (!this->isCachingEnabled()) {
this->clearCache();
}
return converged;
}
@ -945,12 +951,12 @@ namespace storm {
static std::vector<RationalType>* getTemporary(std::vector<RationalType>& rationalX, std::vector<ImpreciseType>*& currentX, std::vector<ImpreciseType>*& newX) {
return &rationalX;
}
static void swapSolutions(std::vector<RationalType>& rationalX, std::vector<RationalType>*& rationalSolution, std::vector<ImpreciseType>& x, std::vector<ImpreciseType>*& currentX, std::vector<ImpreciseType>*& newX) {
// Nothing to do.
}
};
template<typename RationalType>
struct TemporaryHelper<RationalType, RationalType> {
static std::vector<RationalType>* getTemporary(std::vector<RationalType>& rationalX, std::vector<RationalType>*& currentX, std::vector<RationalType>*& newX) {
@ -960,7 +966,7 @@ namespace storm {
static void swapSolutions(std::vector<RationalType>& rationalX, std::vector<RationalType>*& rationalSolution, std::vector<RationalType>& x, std::vector<RationalType>*& currentX, std::vector<RationalType>*& newX) {
if (&rationalX == rationalSolution) {
// In this case, the rational solution is in place.
// However, since the rational solution is no alias to current x, the imprecise solution is stored
// in current x and and rational x is not an alias to x, we can swap the contents of currentX to x.
std::swap(x, *currentX);
@ -971,7 +977,7 @@ namespace storm {
}
}
};
template<typename ValueType>
template<typename RationalType, typename ImpreciseType>
bool IterativeMinMaxLinearEquationSolver<ValueType>::solveEquationsRationalSearchHelper(Environment const& env, OptimizationDirection dir, IterativeMinMaxLinearEquationSolver<ImpreciseType> const& impreciseSolver, storm::storage::SparseMatrix<RationalType> const& rationalA, std::vector<RationalType>& rationalX, std::vector<RationalType> const& rationalB, storm::storage::SparseMatrix<ImpreciseType> const& A, std::vector<ImpreciseType>& x, std::vector<ImpreciseType> const& b, std::vector<ImpreciseType>& tmpX) const {
@ -989,29 +995,29 @@ namespace storm {
while (status == SolverStatus::InProgress && overallIterations < env.solver().minMax().getMaximalNumberOfIterations()) {
// Perform value iteration with the current precision.
typename IterativeMinMaxLinearEquationSolver<ImpreciseType>::ValueIterationResult result = impreciseSolver.performValueIteration(env, dir, currentX, newX, b, storm::utility::convertNumber<ImpreciseType, ValueType>(precision), env.solver().minMax().getRelativeTerminationCriterion(), SolverGuarantee::LessOrEqual, overallIterations, env.solver().minMax().getMaximalNumberOfIterations(), env.solver().minMax().getMultiplicationStyle());
// At this point, the result of the imprecise value iteration is stored in the (imprecise) current x.
++valueIterationInvocations;
STORM_LOG_TRACE("Completed " << valueIterationInvocations << " value iteration invocations, the last one with precision " << precision << " completed in " << result.iterations << " iterations.");
// Count the iterations.
overallIterations += result.iterations;
// Compute maximal precision until which to sharpen.
uint64_t p = storm::utility::convertNumber<uint64_t>(storm::utility::ceil(storm::utility::log10<ValueType>(storm::utility::one<ValueType>() / precision)));
// Make sure that currentX and rationalX are not aliased.
std::vector<RationalType>* temporaryRational = TemporaryHelper<RationalType, ImpreciseType>::getTemporary(rationalX, currentX, newX);
// Sharpen solution and place it in the temporary rational.
bool foundSolution = sharpen(dir, p, rationalA, *currentX, rationalB, *temporaryRational);
// After sharpen, if a solution was found, it is contained in the free rational.
if (foundSolution) {
status = SolverStatus::Converged;
TemporaryHelper<RationalType, ImpreciseType>::swapSolutions(rationalX, temporaryRational, x, currentX, newX);
} else {
// Increase the precision.
@ -1020,14 +1026,14 @@ namespace storm {
status = this->updateStatus(status, false, overallIterations, env.solver().minMax().getMaximalNumberOfIterations());
}
// Swap the two vectors if the current result is not in the original x.
if (currentX != originalX) {
std::swap(x, tmpX);
}
this->reportStatus(status, overallIterations);
return status == SolverStatus::Converged || status == SolverStatus::TerminatedEarly;
}
@ -1035,18 +1041,18 @@ namespace storm {
bool IterativeMinMaxLinearEquationSolver<ValueType>::solveEquationsRationalSearch(Environment const& env, OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
return solveEquationsRationalSearchHelper<double>(env, dir, x, b);
}
template<typename ValueType>
void IterativeMinMaxLinearEquationSolver<ValueType>::computeOptimalValueForRowGroup(uint_fast64_t group, OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, uint_fast64_t* choice) const {
uint64_t row = this->A->getRowGroupIndices()[group];
uint64_t groupEnd = this->A->getRowGroupIndices()[group + 1];
assert(row != groupEnd);
auto bIt = b.begin() + row;
ValueType& xi = x[group];
xi = this->A->multiplyRowWithVector(row, x) + *bIt;
uint64_t optimalRow = row;
for (++row, ++bIt; row < groupEnd; ++row, ++bIt) {
ValueType choiceVal = this->A->multiplyRowWithVector(row, x) + *bIt;
if (minimize(dir)) {
@ -1065,7 +1071,7 @@ namespace storm {
*choice = optimalRow - this->A->getRowGroupIndices()[group];
}
}
template<typename ValueType>
void IterativeMinMaxLinearEquationSolver<ValueType>::clearCache() const {
multiplierA.reset();
@ -1075,9 +1081,9 @@ namespace storm {
optimisticValueIterationHelper.reset();
StandardMinMaxLinearEquationSolver<ValueType>::clearCache();
}
template class IterativeMinMaxLinearEquationSolver<double>;
#ifdef STORM_HAVE_CARL
template class IterativeMinMaxLinearEquationSolver<storm::RationalNumber>;
#endif

Write Preview
Loading…
Cancel
Save