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

made qvi code more readable

tempestpy_adaptions
TimQu 7 years ago
parent
674a30c154
commit
f168df139d
1 changed files with 96 additions and 105 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. 153
      src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp

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

@ -800,21 +800,11 @@ namespace storm {
}
}
bool checkRestartCriterion() {
return false;
// iterations <= restartMaxIterations && (minimize(dir) ? restartThreshold * improvedPrimaryBound > primaryBound : restartThreshold * primaryBound > improvedPrimaryBound
}
bool isPreciseEnough(ValueType const& xi, ValueType const& yi, ValueType const& lb, ValueType const& ub) {
return yi * (ub - lb) <= storm::utility::abs<ValueType>((relative ? (precision * xi) : (precision * storm::utility::convertNumber<ValueType>(2.0))));
}
template<OptimizationDirection dir>
bool checkConvergenceUpdateBounds(uint64_t const& iterations, storm::storage::BitVector const* relevantValues = nullptr) {
bool checkConvergenceUpdateBounds(storm::storage::BitVector const* relevantValues = nullptr) {
if (convergencePhase1) {
if (checkConvergencePhase1()) {
STORM_LOG_INFO("Quick Value Iteration took " << iterations << " iterations for first convergence phase.");
firstIndexViolatingConvergence = 0;
if (relevantValues != nullptr) {
firstIndexViolatingConvergence = relevantValues->getNextSetIndex(firstIndexViolatingConvergence);
@ -829,9 +819,60 @@ namespace storm {
// The difference between lower and upper bound has to be < precision at every (relevant) value
// For efficiency reasons we first check whether it is worth to compute the actual bounds. We do so by considering possibly too tight bounds
ValueType lowerBoundCandidate, upperBoundCandidate;
if (preliminaryConvergenceCheck<dir>(lowerBoundCandidate, upperBoundCandidate)) {
updateLowerUpperBound<dir>(lowerBoundCandidate, upperBoundCandidate);
checkIfDecisionValueBlocks<dir>();
return checkConvergencePhase2<dir>(relevantValues);
}
return false;
}
void setSolutionVector() {
STORM_LOG_WARN_COND(hasLowerBound && hasUpperBound, "No lower or upper result bound could be computed within the given number of Iterations.");
ValueType meanBound = (upperBound + lowerBound) / storm::utility::convertNumber<ValueType>(2.0);
storm::utility::vector::applyPointwise(x, y, x, [&meanBound] (ValueType const& xi, ValueType const& yi) { return xi + yi * meanBound; });
STORM_LOG_INFO("Quick Value Iteration terminated with lower value bound "
<< (hasLowerBound ? lowerBound : storm::utility::zero<ValueType>()) << (hasLowerBound ? "" : "(none)")
<< " and upper value bound "
<< (hasUpperBound ? upperBound : storm::utility::zero<ValueType>()) << (hasUpperBound ? "" : "(none)")
<< ". Decision value is "
<< (hasDecisionValue ? decisionValue : storm::utility::zero<ValueType>()) << (hasDecisionValue ? "" : "(none)")
<< ".");
}
private:
ValueType lowerBoundCandidate = x[minIndex] / (storm::utility::one<ValueType>() - y[minIndex]);
ValueType upperBoundCandidate = x[maxIndex] / (storm::utility::one<ValueType>() - y[maxIndex]);
bool checkConvergencePhase1() {
// Return true if y ('the probability to stay within the matrix') is < 1 at every entry
for (; firstIndexViolatingConvergence != y.size(); ++firstIndexViolatingConvergence) {
static_assert(NumberTraits<ValueType>::IsExact || std::is_same<ValueType, double>::value, "Considered ValueType not handled.");
if (NumberTraits<ValueType>::IsExact) {
if (storm::utility::isOne(y[firstIndexViolatingConvergence])) {
return false;
}
} else {
if (storm::utility::isAlmostOne(storm::utility::convertNumber<double>(y[firstIndexViolatingConvergence]))) {
return false;
}
}
}
convergencePhase1 = false;
return true;
}
bool isPreciseEnough(ValueType const& xi, ValueType const& yi, ValueType const& lb, ValueType const& ub) {
return yi * (ub - lb) <= storm::utility::abs<ValueType>((relative ? (precision * xi) : (precision * storm::utility::convertNumber<ValueType>(2.0))));
}
template<OptimizationDirection dir>
bool preliminaryConvergenceCheck(ValueType& lowerBoundCandidate, ValueType& upperBoundCandidate) {
lowerBoundCandidate = x[minIndex] / (storm::utility::one<ValueType>() - y[minIndex]);
upperBoundCandidate = x[maxIndex] / (storm::utility::one<ValueType>() - y[maxIndex]);
// Make sure that these candidates are at least as tight as the already known bounds
if (hasLowerBound && lowerBoundCandidate < lowerBound) {
lowerBoundCandidate = lowerBound;
@ -839,31 +880,23 @@ namespace storm {
if (hasUpperBound && upperBoundCandidate > upperBound) {
upperBoundCandidate = upperBound;
}
bool computeActualBounds = isPreciseEnough(x[firstIndexViolatingConvergence], y[firstIndexViolatingConvergence], lowerBoundCandidate, upperBoundCandidate);
if (!computeActualBounds) {
if (decisionValueBlocks) {
ValueType improvedPrimaryBound = x[getPrimaryIndex<dir>()] + getPrimaryBound<dir>() * y[getPrimaryIndex<dir>()];
assert(better<dir>(getPrimaryBound<dir>(), improvedPrimaryBound));
computeActualBounds = checkRestartCriterion();
} else {
computeActualBounds = hasDecisionValue && better<dir>(decisionValue, getPrimaryBound<dir>());
if (isPreciseEnough(x[firstIndexViolatingConvergence], y[firstIndexViolatingConvergence], lowerBoundCandidate, upperBoundCandidate)) {
return true;
}
if (!decisionValueBlocks) {
return hasDecisionValue && better<dir>(decisionValue, getPrimaryBound<dir>());
}
return false;
}
if (computeActualBounds) {
template<OptimizationDirection dir>
void updateLowerUpperBound(ValueType& lowerBoundCandidate, ValueType& upperBoundCandidate) {
auto xIt = x.begin();
auto xIte = x.end();
auto yIt = y.begin();
ValueType improvedPrimaryBound;
bool computedImprovedPrimaryBound = false;
for (uint64_t index = 0; xIt != xIte; ++xIt, ++yIt, ++index) {
ValueType currentBound = *xIt / (storm::utility::one<ValueType>() - *yIt);
if (decisionValueBlocks) {
ValueType currentImprovedBound = *xIt + getPrimaryBound<dir>() * (*yIt);
if (!computedImprovedPrimaryBound || better<dir>(currentImprovedBound, improvedPrimaryBound)) {
computedImprovedPrimaryBound = true;
getPrimaryIndex<dir>() = index;
improvedPrimaryBound = std::move(currentImprovedBound);
}
if (better<dir>(getSecondaryBound<dir>(), currentBound)) {
getSecondaryIndex<dir>() = index;
getSecondaryBound<dir>() = std::move(currentBound);
@ -884,13 +917,19 @@ namespace storm {
if ((minimize(dir) || !decisionValueBlocks) && (!hasUpperBound || upperBoundCandidate < upperBound)) {
setUpperBound(upperBoundCandidate);
}
}
template<OptimizationDirection dir>
void checkIfDecisionValueBlocks() {
// Check whether the decision value blocks now (i.e. further improvement of the primary bound would lead to a non-optimal scheduler).
if (!decisionValueBlocks && hasDecisionValue && better<dir>(decisionValue, getPrimaryBound<dir>())) {
getPrimaryBound<dir>() = decisionValue;
decisionValueBlocks = true;
}
}
template<OptimizationDirection dir>
bool checkConvergencePhase2(storm::storage::BitVector const* relevantValues = nullptr) {
// Check whether the desired precision is reached
if (isPreciseEnough(x[firstIndexViolatingConvergence], y[firstIndexViolatingConvergence], lowerBound, upperBound)) {
// The current index satisfies the desired bound. We now move to the next index that violates it
@ -905,55 +944,12 @@ namespace storm {
} else {
if (!isPreciseEnough(x[firstIndexViolatingConvergence], y[firstIndexViolatingConvergence], lowerBound, upperBound)) {
// not converged yet
break;
}
}
}
}
// Check whether we should restart
if (computedImprovedPrimaryBound && checkRestartCriterion()) {
STORM_LOG_INFO("Restarting QVI after " << iterations << " iterations. Improved bound from " << getPrimaryBound<dir>() << " to " << improvedPrimaryBound << ".");
getPrimaryBound<dir>() = improvedPrimaryBound;
restart();
}
}
return false;
}
void setSolutionVector() {
STORM_LOG_WARN_COND(hasLowerBound && hasUpperBound, "No lower or upper result bound could be computed within the given number of Iterations.");
ValueType meanBound = (upperBound + lowerBound) / storm::utility::convertNumber<ValueType>(2.0);
storm::utility::vector::applyPointwise(x, y, x, [&meanBound] (ValueType const& xi, ValueType const& yi) { return xi + yi * meanBound; });
STORM_LOG_INFO("Quick Value Iteration terminated with lower value bound "
<< (hasLowerBound ? lowerBound : storm::utility::zero<ValueType>()) << (hasLowerBound ? "" : "(none)")
<< " and upper value bound "
<< (hasUpperBound ? upperBound : storm::utility::zero<ValueType>()) << (hasUpperBound ? "" : "(none)")
<< ". Decision value is "
<< (hasDecisionValue ? decisionValue : storm::utility::zero<ValueType>()) << (hasDecisionValue ? "" : "(none)")
<< ".");
}
private:
bool checkConvergencePhase1() {
// Return true if y ('the probability to stay within the matrix') is < 1 at every entry
for (; firstIndexViolatingConvergence != y.size(); ++firstIndexViolatingConvergence) {
static_assert(NumberTraits<ValueType>::IsExact || std::is_same<ValueType, double>::value, "Considered ValueType not handled.");
if (NumberTraits<ValueType>::IsExact) {
if (storm::utility::isOne(y[firstIndexViolatingConvergence])) {
return false;
}
} else {
if (storm::utility::isAlmostOne(storm::utility::convertNumber<double>(y[firstIndexViolatingConvergence]))) {
return false;
}
}
}
convergencePhase1 = false;
return true;
return false;
}
std::vector<ValueType>& x;
@ -982,11 +978,6 @@ namespace storm {
QuickValueIterationHelper<ValueType> helper(x, *this->auxiliaryRowGroupVector, env.solver().minMax().getRelativeTerminationCriterion(), storm::utility::convertNumber<ValueType>(env.solver().minMax().getPrecision()), this->A->getSizeOfLargestRowGroup());
// Get the precision
uint64_t restartMaxIterations = env.solver().minMax().getQviRestartMaxIterations();
ValueType restartThreshold = storm::utility::convertNumber<ValueType>(env.solver().minMax().getQviRestartThreshold());
// Prepare initial bounds for the solution (if given)
if (this->hasLowerBound()) {
helper.setLowerBound(this->getLowerBound(true));
@ -1007,13 +998,13 @@ namespace storm {
while (status == SolverStatus::InProgress && iterations < env.solver().minMax().getMaximalNumberOfIterations()) {
if (minimize(dir)) {
helper.template performIterationStep<OptimizationDirection::Minimize>(*this->A, b);
if (helper.template checkConvergenceUpdateBounds<OptimizationDirection::Minimize>(iterations, relevantValuesPtr)) {
if (helper.template checkConvergenceUpdateBounds<OptimizationDirection::Minimize>(relevantValuesPtr)) {
status = SolverStatus::Converged;
}
} else {
assert(maximize(dir));
helper.template performIterationStep<OptimizationDirection::Maximize>(*this->A, b);
if (helper.template checkConvergenceUpdateBounds<OptimizationDirection::Maximize>(iterations, relevantValuesPtr)) {
if (helper.template checkConvergenceUpdateBounds<OptimizationDirection::Maximize>(relevantValuesPtr)) {
status = SolverStatus::Converged;
}
}

Write Preview
Loading…
Cancel
Save