Browse Source

Avoided duplicated code for sound value iteration

tempestpy_adaptions
TimQu 7 years ago
parent
commit
a24de86ce1
  1. 406
      src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp
  2. 213
      src/storm/solver/NativeLinearEquationSolver.cpp
  3. 418
      src/storm/solver/helper/SoundValueIterationHelper.cpp
  4. 147
      src/storm/solver/helper/SoundValueIterationHelper.h

406
src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp

@ -10,6 +10,7 @@
#include "storm/utility/KwekMehlhorn.h"
#include "storm/utility/NumberTraits.h"
#include "storm/solver/helper/SoundValueIterationHelper.h"
#include "storm/utility/Stopwatch.h"
#include "storm/utility/vector.h"
#include "storm/utility/macros.h"
@ -604,395 +605,6 @@ namespace storm {
return status == SolverStatus::Converged;
}
template<typename ValueType>
class SoundValueIterationHelper {
public:
typedef uint32_t IndexType;
SoundValueIterationHelper(storm::storage::SparseMatrix<ValueType> const& matrix, std::vector<ValueType>& x, std::vector<ValueType>& y, bool relative, ValueType const& precision) : x(x), y(y), hasLowerBound(false), hasUpperBound(false), minIndex(0), maxIndex(0), relative(relative), precision(precision), rowGroupIndices(matrix.getRowGroupIndices()) {
STORM_LOG_THROW(matrix.getEntryCount() < std::numeric_limits<IndexType>::max(), storm::exceptions::NotSupportedException, "The number of matrix entries is too large for the selected index type.");
uint64_t sizeOfLargestRowGroup = matrix.getSizeOfLargestRowGroup();
xTmp.resize(sizeOfLargestRowGroup);
yTmp.resize(sizeOfLargestRowGroup);
x.assign(x.size(), storm::utility::zero<ValueType>());
y.assign(x.size(), storm::utility::one<ValueType>());
hasDecisionValue = false;
decisionValueBlocks = false;
convergencePhase1 = true;
firstIndexViolatingConvergence = 0;
numRows = matrix.getRowCount();
matrixValues.clear();
matrixColumns.clear();
rowIndications.clear();
matrixValues.reserve(matrix.getNonzeroEntryCount());
matrixColumns.reserve(matrix.getColumnCount());
rowIndications.reserve(numRows + 1);
rowIndications.push_back(0);
for (IndexType r = 0; r < numRows; ++r) {
for (auto const& entry : matrix.getRow(r)) {
matrixValues.push_back(entry.getValue());
matrixColumns.push_back(entry.getColumn());
}
rowIndications.push_back(matrixValues.size());
}
}
inline void setLowerBound(ValueType const& value) {
hasLowerBound = true;
lowerBound = value;
}
inline void setUpperBound(ValueType const& value) {
hasUpperBound = true;
upperBound = value;
}
template<OptimizationDirection dir>
inline bool better(ValueType const& val1, ValueType const& val2) {
return maximize(dir) ? val1 > val2 : val1 < val2;
}
template<OptimizationDirection dir>
inline ValueType& getPrimaryBound() {
return maximize(dir) ? upperBound : lowerBound;
}
template<OptimizationDirection dir>
inline bool& hasPrimaryBound() {
return maximize(dir) ? hasUpperBound : hasLowerBound;
}
template<OptimizationDirection dir>
inline ValueType& getSecondaryBound() {
return maximize(dir) ? lowerBound : upperBound;
}
template<OptimizationDirection dir>
inline uint64_t& getPrimaryIndex() {
return maximize(dir) ? maxIndex : minIndex;
}
template<OptimizationDirection dir>
inline uint64_t& getSecondaryIndex() {
return maximize(dir) ? minIndex : maxIndex;
}
void multiplyRow(IndexType const& rowIndex, ValueType const& bi, ValueType& xi, ValueType& yi) {
assert(rowIndex < numRows);
ValueType xRes = bi;
ValueType yRes = storm::utility::zero<ValueType>();
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];
yRes += *entryIt * y[*colIt];
}
xi = std::move(xRes);
yi = std::move(yRes);
}
template<OptimizationDirection dir>
void performIterationStep(std::vector<ValueType> const& b) {
if (!decisionValueBlocks) {
performIterationStepUpdateDecisionValue<dir>(b);
} else {
assert(decisionValue == getPrimaryBound<dir>());
auto xIt = x.rbegin();
auto yIt = y.rbegin();
auto groupStartIt = rowGroupIndices.rbegin();
uint64_t groupEnd = *groupStartIt;
++groupStartIt;
for (auto groupStartIte = rowGroupIndices.rend(); groupStartIt != groupStartIte; groupEnd = *(groupStartIt++), ++xIt, ++yIt) {
// Perform the iteration for the first row in the group
IndexType row = *groupStartIt;
ValueType xBest, yBest;
multiplyRow(row, b[row], xBest, yBest);
++row;
// Only do more work if there are still rows in this row group
if (row != groupEnd) {
ValueType xi, yi;
ValueType bestValue = xBest + yBest * getPrimaryBound<dir>();
for (;row < groupEnd; ++row) {
// Get the multiplication results
multiplyRow(row, b[row], xi, yi);
ValueType currentValue = xi + yi * getPrimaryBound<dir>();
// Check if the current row is better then the previously found one
if (better<dir>(currentValue, bestValue)) {
xBest = std::move(xi);
yBest = std::move(yi);
bestValue = std::move(currentValue);
} else if (currentValue == bestValue && yBest > yi) {
// If the value for this row is not strictly better, it might still be equal and have a better y value
xBest = std::move(xi);
yBest = std::move(yi);
}
}
}
*xIt = std::move(xBest);
*yIt = std::move(yBest);
}
}
}
template<OptimizationDirection dir>
void performIterationStepUpdateDecisionValue(std::vector<ValueType> const& b) {
auto xIt = x.rbegin();
auto yIt = y.rbegin();
auto groupStartIt = rowGroupIndices.rbegin();
uint64_t groupEnd = *groupStartIt;
++groupStartIt;
for (auto groupStartIte = rowGroupIndices.rend(); groupStartIt != groupStartIte; groupEnd = *(groupStartIt++), ++xIt, ++yIt) {
// Perform the iteration for the first row in the group
uint64_t row = *groupStartIt;
ValueType xBest, yBest;
multiplyRow(row, b[row], xBest, yBest);
++row;
// Only do more work if there are still rows in this row group
if (row != groupEnd) {
ValueType xi, yi;
uint64_t xyTmpIndex = 0;
if (hasPrimaryBound<dir>()) {
ValueType bestValue = xBest + yBest * getPrimaryBound<dir>();
for (;row < groupEnd; ++row) {
// Get the multiplication results
multiplyRow(row, b[row], xi, yi);
ValueType currentValue = xi + yi * getPrimaryBound<dir>();
// Check if the current row is better then the previously found one
if (better<dir>(currentValue, bestValue)) {
if (yBest < yi) {
// We need to store the 'old' best value as it might be relevant for the decision value
xTmp[xyTmpIndex] = std::move(xBest);
yTmp[xyTmpIndex] = std::move(yBest);
++xyTmpIndex;
}
xBest = std::move(xi);
yBest = std::move(yi);
bestValue = std::move(currentValue);
} else if (yBest > yi) {
// If the value for this row is not strictly better, it might still be equal and have a better y value
if (currentValue == bestValue) {
xBest = std::move(xi);
yBest = std::move(yi);
} else {
xTmp[xyTmpIndex] = std::move(xi);
yTmp[xyTmpIndex] = std::move(yi);
++xyTmpIndex;
}
}
}
} else {
for (;row < groupEnd; ++row) {
multiplyRow(row, b[row], xi, yi);
// Update the best choice
if (yi > yBest || (yi == yBest && better<dir>(xi, xBest))) {
xTmp[xyTmpIndex] = std::move(xBest);
yTmp[xyTmpIndex] = std::move(yBest);
++xyTmpIndex;
xBest = std::move(xi);
yBest = std::move(yi);
} else {
xTmp[xyTmpIndex] = std::move(xi);
yTmp[xyTmpIndex] = std::move(yi);
++xyTmpIndex;
}
}
}
// Update the decision value
for (uint64_t i = 0; i < xyTmpIndex; ++i) {
ValueType deltaY = yBest - yTmp[i];
if (deltaY > storm::utility::zero<ValueType>()) {
ValueType newDecisionValue = (xTmp[i] - xBest) / deltaY;
if (!hasDecisionValue || better<dir>(newDecisionValue, decisionValue)) {
decisionValue = std::move(newDecisionValue);
hasDecisionValue = true;
}
}
}
}
*xIt = std::move(xBest);
*yIt = std::move(yBest);
}
}
template<OptimizationDirection dir>
bool checkConvergenceUpdateBounds(storm::storage::BitVector const* relevantValues = nullptr) {
if (convergencePhase1) {
if (checkConvergencePhase1()) {
firstIndexViolatingConvergence = 0;
if (relevantValues != nullptr) {
firstIndexViolatingConvergence = relevantValues->getNextSetIndex(firstIndexViolatingConvergence);
}
} else {
return false;
}
}
STORM_LOG_ASSERT(!std::any_of(y.begin(), y.end(), [](ValueType value){return storm::utility::isOne(value);}), "Did not expect staying-probability 1 at this point.");
// Reaching this point means that we are in Phase 2:
// 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("Sound 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;
}
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;
}
if (hasUpperBound && upperBoundCandidate > upperBound) {
upperBoundCandidate = upperBound;
}
if (isPreciseEnough(x[firstIndexViolatingConvergence], y[firstIndexViolatingConvergence], lowerBoundCandidate, upperBoundCandidate)) {
return true;
}
if (!decisionValueBlocks) {
return hasDecisionValue && better<dir>(decisionValue, getPrimaryBound<dir>());
}
return false;
}
template<OptimizationDirection dir>
void updateLowerUpperBound(ValueType& lowerBoundCandidate, ValueType& upperBoundCandidate) {
auto xIt = x.begin();
auto xIte = x.end();
auto yIt = y.begin();
for (uint64_t index = 0; xIt != xIte; ++xIt, ++yIt, ++index) {
ValueType currentBound = *xIt / (storm::utility::one<ValueType>() - *yIt);
if (decisionValueBlocks) {
if (better<dir>(getSecondaryBound<dir>(), currentBound)) {
getSecondaryIndex<dir>() = index;
getSecondaryBound<dir>() = std::move(currentBound);
}
} else {
if (currentBound < lowerBoundCandidate) {
minIndex = index;
lowerBoundCandidate = std::move(currentBound);
} else if (currentBound > upperBoundCandidate) {
maxIndex = index;
upperBoundCandidate = std::move(currentBound);
}
}
}
if ((maximize(dir) || !decisionValueBlocks) && (!hasLowerBound || lowerBoundCandidate > lowerBound)) {
setLowerBound(lowerBoundCandidate);
}
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
while (true) {
++firstIndexViolatingConvergence;
if (relevantValues != nullptr) {
firstIndexViolatingConvergence = relevantValues->getNextSetIndex(firstIndexViolatingConvergence);
}
if (firstIndexViolatingConvergence == x.size()) {
// Converged!
return true;
} else {
if (!isPreciseEnough(x[firstIndexViolatingConvergence], y[firstIndexViolatingConvergence], lowerBound, upperBound)) {
// not converged yet
return false;
}
}
}
}
return false;
}
std::vector<ValueType>& x;
std::vector<ValueType>& y;
std::vector<ValueType> xTmp, yTmp;
ValueType lowerBound, upperBound, decisionValue;
bool hasLowerBound, hasUpperBound, hasDecisionValue;
uint64_t minIndex, maxIndex;
bool convergencePhase1;
bool decisionValueBlocks;
uint64_t firstIndexViolatingConvergence;
bool relative;
ValueType precision;
std::vector<ValueType> matrixValues;
std::vector<IndexType> matrixColumns;
std::vector<IndexType> rowIndications;
std::vector<typename storm::storage::SparseMatrix<ValueType>::index_type> const& rowGroupIndices;
IndexType numRows;
};
template<typename ValueType>
bool IterativeMinMaxLinearEquationSolver<ValueType>::solveEquationsSoundValueIteration(Environment const& env, OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
@ -1003,7 +615,7 @@ namespace storm {
}
// TODO: implement caching for the helper
SoundValueIterationHelper<ValueType> helper(*this->A, x, *this->auxiliaryRowGroupVector, env.solver().minMax().getRelativeTerminationCriterion(), storm::utility::convertNumber<ValueType>(env.solver().minMax().getPrecision()));
storm::solver::helper::SoundValueIterationHelper<ValueType> helper(*this->A, 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()) {
@ -1023,17 +635,9 @@ namespace storm {
uint64_t iterations = 0;
while (status == SolverStatus::InProgress && iterations < env.solver().minMax().getMaximalNumberOfIterations()) {
if (minimize(dir)) {
helper.template performIterationStep<OptimizationDirection::Minimize>(b);
if (helper.template checkConvergenceUpdateBounds<OptimizationDirection::Minimize>(relevantValuesPtr)) {
status = SolverStatus::Converged;
}
} else {
assert(maximize(dir));
helper.template performIterationStep<OptimizationDirection::Maximize>(b);
if (helper.template checkConvergenceUpdateBounds<OptimizationDirection::Maximize>(relevantValuesPtr)) {
status = SolverStatus::Converged;
}
helper.performIterationStep(dir, b);
if (helper.checkConvergenceUpdateBounds(dir, relevantValuesPtr)) {
status = SolverStatus::Converged;
}
// Update environment variables.

213
src/storm/solver/NativeLinearEquationSolver.cpp

@ -9,6 +9,7 @@
#include "storm/utility/NumberTraits.h"
#include "storm/utility/constants.h"
#include "storm/utility/vector.h"
#include "storm/solver/helper/SoundValueIterationHelper.h"
#include "storm/solver/Multiplier.h"
#include "storm/exceptions/InvalidStateException.h"
#include "storm/exceptions/InvalidEnvironmentException.h"
@ -563,216 +564,6 @@ namespace storm {
return converged;
}
template<typename ValueType>
class SoundPowerHelper {
public:
typedef uint32_t IndexType;
SoundPowerHelper(storm::storage::SparseMatrix<ValueType> const& matrix, std::vector<ValueType>& x, std::vector<ValueType>& y, bool relative, ValueType const& precision) : x(x), y(y), hasLowerBound(false), hasUpperBound(false), minIndex(0), maxIndex(0), relative(relative), precision(precision) {
STORM_LOG_THROW(matrix.getEntryCount() < std::numeric_limits<IndexType>::max(), storm::exceptions::NotSupportedException, "The number of matrix entries is too large for the selected index type.");
x.assign(x.size(), storm::utility::zero<ValueType>());
y.assign(x.size(), storm::utility::one<ValueType>());
convergencePhase1 = true;
firstIndexViolatingConvergence = 0;
numRows = matrix.getRowCount();
matrixValues.clear();
matrixColumns.clear();
rowIndications.clear();
matrixValues.reserve(matrix.getNonzeroEntryCount());
matrixColumns.reserve(matrix.getColumnCount());
rowIndications.reserve(numRows + 1);
rowIndications.push_back(0);
for (IndexType r = 0; r < numRows; ++r) {
for (auto const& entry : matrix.getRow(r)) {
matrixValues.push_back(entry.getValue());
matrixColumns.push_back(entry.getColumn());
}
rowIndications.push_back(matrixValues.size());
}
}
inline void setLowerBound(ValueType const& value) {
hasLowerBound = true;
lowerBound = value;
}
inline void setUpperBound(ValueType const& value) {
hasUpperBound = true;
upperBound = value;
}
void multiplyRow(IndexType const& rowIndex, ValueType const& bi, ValueType& xi, ValueType& yi) {
assert(rowIndex < numRows);
ValueType xRes = bi;
ValueType yRes = storm::utility::zero<ValueType>();
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];
yRes += *entryIt * y[*colIt];
}
xi = std::move(xRes);
yi = std::move(yRes);
}
void performIterationStep(std::vector<ValueType> const& b) {
auto xIt = x.rbegin();
auto yIt = y.rbegin();
IndexType row = numRows;
while (row > 0) {
--row;
multiplyRow(row, b[row], *xIt, *yIt);
++xIt;
++yIt;
}
}
bool checkConvergenceUpdateBounds(storm::storage::BitVector const* relevantValues = nullptr) {
if (convergencePhase1) {
if (checkConvergencePhase1()) {
firstIndexViolatingConvergence = 0;
if (relevantValues != nullptr) {
firstIndexViolatingConvergence = relevantValues->getNextSetIndex(firstIndexViolatingConvergence);
}
} else {
return false;
}
}
STORM_LOG_ASSERT(!std::any_of(y.begin(), y.end(), [](ValueType value){return storm::utility::isOne(value);}), "Did not expect staying-probability 1 at this point.");
// Reaching this point means that we are in Phase 2:
// 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(lowerBoundCandidate, upperBoundCandidate)) {
updateLowerUpperBound(lowerBoundCandidate, upperBoundCandidate);
return checkConvergencePhase2(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("Sound Power 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)")
<< ".");
}
private:
bool checkConvergencePhase1() {
// Return true if y ('the probability to stay within the 'maybestates') 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))));
}
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;
}
if (hasUpperBound && upperBoundCandidate > upperBound) {
upperBoundCandidate = upperBound;
}
if (isPreciseEnough(x[firstIndexViolatingConvergence], y[firstIndexViolatingConvergence], lowerBoundCandidate, upperBoundCandidate)) {
return true;
}
return false;
}
void updateLowerUpperBound(ValueType& lowerBoundCandidate, ValueType& upperBoundCandidate) {
auto xIt = x.begin();
auto xIte = x.end();
auto yIt = y.begin();
for (uint64_t index = 0; xIt != xIte; ++xIt, ++yIt, ++index) {
ValueType currentBound = *xIt / (storm::utility::one<ValueType>() - *yIt);
if (currentBound < lowerBoundCandidate) {
minIndex = index;
lowerBoundCandidate = std::move(currentBound);
} else if (currentBound > upperBoundCandidate) {
maxIndex = index;
upperBoundCandidate = std::move(currentBound);
}
}
if (!hasLowerBound || lowerBoundCandidate > lowerBound) {
setLowerBound(lowerBoundCandidate);
}
if (!hasUpperBound || upperBoundCandidate < upperBound) {
setUpperBound(upperBoundCandidate);
}
}
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
while (true) {
++firstIndexViolatingConvergence;
if (relevantValues != nullptr) {
firstIndexViolatingConvergence = relevantValues->getNextSetIndex(firstIndexViolatingConvergence);
}
if (firstIndexViolatingConvergence == x.size()) {
// Converged!
return true;
} else {
if (!isPreciseEnough(x[firstIndexViolatingConvergence], y[firstIndexViolatingConvergence], lowerBound, upperBound)) {
// not converged yet
return false;
}
}
}
}
return false;
}
std::vector<ValueType>& x;
std::vector<ValueType>& y;
ValueType lowerBound, upperBound, decisionValue;
bool hasLowerBound, hasUpperBound, hasDecisionValue;
uint64_t minIndex, maxIndex;
bool convergencePhase1;
uint64_t firstIndexViolatingConvergence;
std::vector<ValueType> matrixValues;
std::vector<IndexType> matrixColumns;
std::vector<IndexType> rowIndications;
IndexType numRows;
bool relative;
ValueType precision;
};
template<typename ValueType>
bool NativeLinearEquationSolver<ValueType>::solveEquationsSoundPower(Environment const& env, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
@ -784,7 +575,7 @@ namespace storm {
}
// TODO: implement caching for the helper
SoundPowerHelper<ValueType> helper(*this->A, x, *this->cachedRowVector, env.solver().native().getRelativeTerminationCriterion(), storm::utility::convertNumber<ValueType>(env.solver().native().getPrecision()));
storm::solver::helper::SoundValueIterationHelper<ValueType> helper(*this->A, x, *this->cachedRowVector, env.solver().native().getRelativeTerminationCriterion(), storm::utility::convertNumber<ValueType>(env.solver().native().getPrecision()));
// Prepare initial bounds for the solution (if given)
if (this->hasLowerBound()) {

418
src/storm/solver/helper/SoundValueIterationHelper.cpp

@ -0,0 +1,418 @@
#include "storm/solver/helper/SoundValueIterationHelper.h"
#include "storm/storage/SparseMatrix.h"
#include "storm/storage/BitVector.h"
#include "storm/utility/vector.h"
#include "storm/utility/macros.h"
#include "storm/utility/NumberTraits.h"
#include "storm/exceptions/NotSupportedException.h"
namespace storm {
namespace solver {
namespace helper {
template<typename ValueType>
SoundValueIterationHelper<ValueType>::SoundValueIterationHelper(storm::storage::SparseMatrix<ValueType> const& matrix, std::vector<ValueType>& x, std::vector<ValueType>& y, bool relative, ValueType const& precision) : x(x), y(y), hasLowerBound(false), hasUpperBound(false), hasDecisionValue(false), convergencePhase1(true), decisionValueBlocks(false), firstIndexViolatingConvergence(0), minIndex(0), maxIndex(0), relative(relative), precision(precision), rowGroupIndices(nullptr) {
STORM_LOG_THROW(matrix.getEntryCount() < std::numeric_limits<IndexType>::max(), storm::exceptions::NotSupportedException, "The number of matrix entries is too large for the selected index type.");
if (!matrix.hasTrivialRowGrouping()) {
rowGroupIndices = &matrix.getRowGroupIndices();
uint64_t sizeOfLargestRowGroup = matrix.getSizeOfLargestRowGroup();
xTmp.resize(sizeOfLargestRowGroup);
yTmp.resize(sizeOfLargestRowGroup);
}
x.assign(x.size(), storm::utility::zero<ValueType>());
y.assign(x.size(), storm::utility::one<ValueType>());
numRows = matrix.getRowCount();
matrixValues.clear();
matrixColumns.clear();
rowIndications.clear();
matrixValues.reserve(matrix.getNonzeroEntryCount());
matrixColumns.reserve(matrix.getColumnCount());
rowIndications.reserve(numRows + 1);
rowIndications.push_back(0);
for (IndexType r = 0; r < numRows; ++r) {
for (auto const& entry : matrix.getRow(r)) {
matrixValues.push_back(entry.getValue());
matrixColumns.push_back(entry.getColumn());
}
rowIndications.push_back(matrixValues.size());
}
}
template<typename ValueType>
SoundValueIterationHelper<ValueType>::SoundValueIterationHelper(SoundValueIterationHelper<ValueType>&& oldHelper, std::vector<ValueType>& x, std::vector<ValueType>& y, bool relative, ValueType const& precision) : x(x), y(y), xTmp(std::move(oldHelper.xTmp)), yTmp(std::move(oldHelper.yTmp)), hasLowerBound(false), hasUpperBound(false), hasDecisionValue(false), convergencePhase1(true), decisionValueBlocks(false), firstIndexViolatingConvergence(0), minIndex(0), maxIndex(0), relative(relative), precision(precision), numRows(std::move(oldHelper.numRows)), matrixValues(std::move(oldHelper.matrixValues)), matrixColumns(std::move(oldHelper.matrixColumns)), rowIndications(std::move(oldHelper.rowIndications)), rowGroupIndices(oldHelper.rowGroupIndices) {
x.assign(x.size(), storm::utility::zero<ValueType>());
y.assign(x.size(), storm::utility::one<ValueType>());
}
template<typename ValueType>
void SoundValueIterationHelper<ValueType>::setLowerBound(ValueType const& value) {
hasLowerBound = true;
lowerBound = value;
}
template<typename ValueType>
void SoundValueIterationHelper<ValueType>::setUpperBound(ValueType const& value) {
hasUpperBound = true;
upperBound = value;
}
template<typename ValueType>
void SoundValueIterationHelper<ValueType>::multiplyRow(IndexType const& rowIndex, ValueType const& bi, ValueType& xi, ValueType& yi) {
assert(rowIndex < numRows);
ValueType xRes = bi;
ValueType yRes = storm::utility::zero<ValueType>();
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];
yRes += *entryIt * y[*colIt];
}
xi = std::move(xRes);
yi = std::move(yRes);
}
template<typename ValueType>
void SoundValueIterationHelper<ValueType>::performIterationStep(OptimizationDirection const& dir, std::vector<ValueType> const& b) {
if (rowGroupIndices) {
if (minimize(dir)) {
performIterationStep<InternalOptimizationDirection::Minimize>(b);
} else {
performIterationStep<InternalOptimizationDirection::Maximize>(b);
}
} else {
performIterationStep(b);
}
}
template<typename ValueType>
void SoundValueIterationHelper<ValueType>::performIterationStep(std::vector<ValueType> const& b) {
auto xIt = x.rbegin();
auto yIt = y.rbegin();
IndexType row = numRows;
while (row > 0) {
--row;
multiplyRow(row, b[row], *xIt, *yIt);
++xIt;
++yIt;
}
}
template<typename ValueType>
template<typename SoundValueIterationHelper<ValueType>::InternalOptimizationDirection dir>
void SoundValueIterationHelper<ValueType>::performIterationStep(std::vector<ValueType> const& b) {
if (!decisionValueBlocks) {
performIterationStepUpdateDecisionValue<dir>(b);
} else {
assert(decisionValue == getPrimaryBound<dir>());
auto xIt = x.rbegin();
auto yIt = y.rbegin();
auto groupStartIt = rowGroupIndices->rbegin();
uint64_t groupEnd = *groupStartIt;
++groupStartIt;
for (auto groupStartIte = rowGroupIndices->rend(); groupStartIt != groupStartIte; groupEnd = *(groupStartIt++), ++xIt, ++yIt) {
// Perform the iteration for the first row in the group
IndexType row = *groupStartIt;
ValueType xBest, yBest;
multiplyRow(row, b[row], xBest, yBest);
++row;
// Only do more work if there are still rows in this row group
if (row != groupEnd) {
ValueType xi, yi;
ValueType bestValue = xBest + yBest * getPrimaryBound<dir>();
for (;row < groupEnd; ++row) {
// Get the multiplication results
multiplyRow(row, b[row], xi, yi);
ValueType currentValue = xi + yi * getPrimaryBound<dir>();
// Check if the current row is better then the previously found one
if (better<dir>(currentValue, bestValue)) {
xBest = std::move(xi);
yBest = std::move(yi);
bestValue = std::move(currentValue);
} else if (currentValue == bestValue && yBest > yi) {
// If the value for this row is not strictly better, it might still be equal and have a better y value
xBest = std::move(xi);
yBest = std::move(yi);
}
}
}
*xIt = std::move(xBest);
*yIt = std::move(yBest);
}
}
}
template<typename ValueType>
template<typename SoundValueIterationHelper<ValueType>::InternalOptimizationDirection dir>
void SoundValueIterationHelper<ValueType>::performIterationStepUpdateDecisionValue(std::vector<ValueType> const& b) {
auto xIt = x.rbegin();
auto yIt = y.rbegin();
auto groupStartIt = rowGroupIndices->rbegin();
uint64_t groupEnd = *groupStartIt;
++groupStartIt;
for (auto groupStartIte = rowGroupIndices->rend(); groupStartIt != groupStartIte; groupEnd = *(groupStartIt++), ++xIt, ++yIt) {
// Perform the iteration for the first row in the group
uint64_t row = *groupStartIt;
ValueType xBest, yBest;
multiplyRow(row, b[row], xBest, yBest);
++row;
// Only do more work if there are still rows in this row group
if (row != groupEnd) {
ValueType xi, yi;
uint64_t xyTmpIndex = 0;
if (hasPrimaryBound<dir>()) {
ValueType bestValue = xBest + yBest * getPrimaryBound<dir>();
for (;row < groupEnd; ++row) {
// Get the multiplication results
multiplyRow(row, b[row], xi, yi);
ValueType currentValue = xi + yi * getPrimaryBound<dir>();
// Check if the current row is better then the previously found one
if (better<dir>(currentValue, bestValue)) {
if (yBest < yi) {
// We need to store the 'old' best value as it might be relevant for the decision value
xTmp[xyTmpIndex] = std::move(xBest);
yTmp[xyTmpIndex] = std::move(yBest);
++xyTmpIndex;
}
xBest = std::move(xi);
yBest = std::move(yi);
bestValue = std::move(currentValue);
} else if (yBest > yi) {
// If the value for this row is not strictly better, it might still be equal and have a better y value
if (currentValue == bestValue) {
xBest = std::move(xi);
yBest = std::move(yi);
} else {
xTmp[xyTmpIndex] = std::move(xi);
yTmp[xyTmpIndex] = std::move(yi);
++xyTmpIndex;
}
}
}
} else {
for (;row < groupEnd; ++row) {
multiplyRow(row, b[row], xi, yi);
// Update the best choice
if (yi > yBest || (yi == yBest && better<dir>(xi, xBest))) {
xTmp[xyTmpIndex] = std::move(xBest);
yTmp[xyTmpIndex] = std::move(yBest);
++xyTmpIndex;
xBest = std::move(xi);
yBest = std::move(yi);
} else {
xTmp[xyTmpIndex] = std::move(xi);
yTmp[xyTmpIndex] = std::move(yi);
++xyTmpIndex;
}
}
}
// Update the decision value
for (uint64_t i = 0; i < xyTmpIndex; ++i) {
ValueType deltaY = yBest - yTmp[i];
if (deltaY > storm::utility::zero<ValueType>()) {
ValueType newDecisionValue = (xTmp[i] - xBest) / deltaY;
if (!hasDecisionValue || better<dir>(newDecisionValue, decisionValue)) {
decisionValue = std::move(newDecisionValue);
hasDecisionValue = true;
}
}
}
}
*xIt = std::move(xBest);
*yIt = std::move(yBest);
}
}
template<typename ValueType>
bool SoundValueIterationHelper<ValueType>::checkConvergenceUpdateBounds(OptimizationDirection const& dir, storm::storage::BitVector const* relevantValues) {
if (rowGroupIndices) {
if (minimize(dir)) {
return checkConvergenceUpdateBounds<InternalOptimizationDirection::Minimize>(relevantValues);
} else {
return checkConvergenceUpdateBounds<InternalOptimizationDirection::Maximize>(relevantValues);
}
} else {
return checkConvergenceUpdateBounds(relevantValues);
}
}
template<typename ValueType>
bool SoundValueIterationHelper<ValueType>::checkConvergenceUpdateBounds(storm::storage::BitVector const* relevantValues) {
return checkConvergenceUpdateBounds<InternalOptimizationDirection::None>(relevantValues);
}
template<typename ValueType>
template<typename SoundValueIterationHelper<ValueType>::InternalOptimizationDirection dir>
bool SoundValueIterationHelper<ValueType>::checkConvergenceUpdateBounds(storm::storage::BitVector const* relevantValues) {
if (convergencePhase1) {
if (checkConvergencePhase1()) {
firstIndexViolatingConvergence = 0;
if (relevantValues != nullptr) {
firstIndexViolatingConvergence = relevantValues->getNextSetIndex(firstIndexViolatingConvergence);
}
} else {
return false;
}
}
STORM_LOG_ASSERT(!std::any_of(y.begin(), y.end(), [](ValueType value){return storm::utility::isOne(value);}), "Did not expect staying-probability 1 at this point.");
// Reaching this point means that we are in Phase 2:
// 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);
if (dir != InternalOptimizationDirection::None) {
checkIfDecisionValueBlocks<dir>();
}
return checkConvergencePhase2(relevantValues);
}
return false;
}
template<typename ValueType>
void SoundValueIterationHelper<ValueType>::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("Sound 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)")
<< ".");
}
template<typename ValueType>
bool SoundValueIterationHelper<ValueType>::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;
}
template<typename ValueType>
bool SoundValueIterationHelper<ValueType>::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<typename ValueType>
template<typename SoundValueIterationHelper<ValueType>::InternalOptimizationDirection dir>
bool SoundValueIterationHelper<ValueType>::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;
}
if (hasUpperBound && upperBoundCandidate > upperBound) {
upperBoundCandidate = upperBound;
}
if (isPreciseEnough(x[firstIndexViolatingConvergence], y[firstIndexViolatingConvergence], lowerBoundCandidate, upperBoundCandidate)) {
return true;
}
if (dir != InternalOptimizationDirection::None && !decisionValueBlocks) {
return hasDecisionValue && better<dir>(decisionValue, getPrimaryBound<dir>());
}
return false;
}
template<typename ValueType>
template<typename SoundValueIterationHelper<ValueType>::InternalOptimizationDirection dir>
void SoundValueIterationHelper<ValueType>::updateLowerUpperBound(ValueType& lowerBoundCandidate, ValueType& upperBoundCandidate) {
auto xIt = x.begin();
auto xIte = x.end();
auto yIt = y.begin();
for (uint64_t index = 0; xIt != xIte; ++xIt, ++yIt, ++index) {
ValueType currentBound = *xIt / (storm::utility::one<ValueType>() - *yIt);
if (dir != InternalOptimizationDirection::None && decisionValueBlocks) {
if (better<dir>(getSecondaryBound<dir>(), currentBound)) {
getSecondaryIndex<dir>() = index;
getSecondaryBound<dir>() = std::move(currentBound);
}
} else {
if (currentBound < lowerBoundCandidate) {
minIndex = index;
lowerBoundCandidate = std::move(currentBound);
} else if (currentBound > upperBoundCandidate) {
maxIndex = index;
upperBoundCandidate = std::move(currentBound);
}
}
}
if ((dir != InternalOptimizationDirection::Minimize || !decisionValueBlocks) && (!hasLowerBound || lowerBoundCandidate > lowerBound)) {
setLowerBound(lowerBoundCandidate);
}
if ((dir != InternalOptimizationDirection::Maximize || !decisionValueBlocks) && (!hasUpperBound || upperBoundCandidate < upperBound)) {
setUpperBound(upperBoundCandidate);
}
}
template<typename ValueType>
template<typename SoundValueIterationHelper<ValueType>::InternalOptimizationDirection dir>
void SoundValueIterationHelper<ValueType>::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<typename ValueType>
bool SoundValueIterationHelper<ValueType>::checkConvergencePhase2(storm::storage::BitVector const* relevantValues) {
// 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
while (true) {
++firstIndexViolatingConvergence;
if (relevantValues != nullptr) {
firstIndexViolatingConvergence = relevantValues->getNextSetIndex(firstIndexViolatingConvergence);
}
if (firstIndexViolatingConvergence == x.size()) {
// Converged!
return true;
} else {
if (!isPreciseEnough(x[firstIndexViolatingConvergence], y[firstIndexViolatingConvergence], lowerBound, upperBound)) {
// not converged yet
return false;
}
}
}
}
return false;
}
template class SoundValueIterationHelper<double>;
template class SoundValueIterationHelper<storm::RationalNumber>;
}
}
}

147
src/storm/solver/helper/SoundValueIterationHelper.h

@ -0,0 +1,147 @@
#pragma once
#include <vector>
#include "storm/solver/OptimizationDirection.h"
namespace storm {
namespace storage {
template<typename ValueType>
class SparseMatrix;
class BitVector;
}
namespace solver {
namespace helper {
template<typename ValueType>
class SoundValueIterationHelper {
public:
typedef uint32_t IndexType;
/*!
* Creates a new helper from the given data
*/
SoundValueIterationHelper(storm::storage::SparseMatrix<ValueType> const& matrix, std::vector<ValueType>& x, std::vector<ValueType>& y, bool relative, ValueType const& precision);
/*!
* Creates a helper from the given data, considering the same matrix as the given old helper
*/
SoundValueIterationHelper(SoundValueIterationHelper<ValueType>&& oldHelper, std::vector<ValueType>& x, std::vector<ValueType>& y, bool relative, ValueType const& precision);
/*!
* Sets the currently known lower / upper bound
*/
void setLowerBound(ValueType const& value);
void setUpperBound(ValueType const& value);
void setSolutionVector();
/*!
* Performs one iteration step with respect to the given optimization direction.
*/
void performIterationStep(OptimizationDirection const& dir, std::vector<ValueType> const& b);
/*!
* Performs one iteration step, assuming that the row grouping of the initial matrix is trivial.
*/
void performIterationStep(std::vector<ValueType> const& b);
/*!
* Checks for convergence and updates the known lower/upper bounds.
*/
bool checkConvergenceUpdateBounds(OptimizationDirection const& dir, storm::storage::BitVector const* relevantValues = nullptr);
/*!
* Checks for convergence and updates the known lower/upper bounds, assuming that the row grouping of the initial matrix is trivial.
*/
bool checkConvergenceUpdateBounds(storm::storage::BitVector const* relevantValues = nullptr);
private:
enum class InternalOptimizationDirection {
None, Minimize, Maximize
};
template<InternalOptimizationDirection dir>
void performIterationStep(std::vector<ValueType> const& b);
template<InternalOptimizationDirection dir>
void performIterationStepUpdateDecisionValue(std::vector<ValueType> const& b);
void multiplyRow(IndexType const& rowIndex, ValueType const& bi, ValueType& xi, ValueType& yi);
template<InternalOptimizationDirection dir>
bool checkConvergenceUpdateBounds(storm::storage::BitVector const* relevantValues = nullptr);
bool checkConvergencePhase1();
bool checkConvergencePhase2(storm::storage::BitVector const* relevantValues = nullptr);
bool isPreciseEnough(ValueType const& xi, ValueType const& yi, ValueType const& lb, ValueType const& ub);
template<InternalOptimizationDirection dir>
bool preliminaryConvergenceCheck(ValueType& lowerBoundCandidate, ValueType& upperBoundCandidate);
template<InternalOptimizationDirection dir>
void updateLowerUpperBound(ValueType& lowerBoundCandidate, ValueType& upperBoundCandidate);
template<InternalOptimizationDirection dir>
void checkIfDecisionValueBlocks();
// Auxiliary helper functions to avoid case distinctions due to different optimization directions
template<InternalOptimizationDirection dir>
inline bool better(ValueType const& val1, ValueType const& val2) {
return (dir == InternalOptimizationDirection::Maximize) ? val1 > val2 : val1 < val2;
}
template<InternalOptimizationDirection dir>
inline ValueType& getPrimaryBound() {
return (dir == InternalOptimizationDirection::Maximize) ? upperBound : lowerBound;
}
template<InternalOptimizationDirection dir>
inline bool& hasPrimaryBound() {
return (dir == InternalOptimizationDirection::Maximize) ? hasUpperBound : hasLowerBound;
}
template<InternalOptimizationDirection dir>
inline ValueType& getSecondaryBound() {
return (dir == InternalOptimizationDirection::Maximize) ? lowerBound : upperBound;
}
template<InternalOptimizationDirection dir>
inline uint64_t& getPrimaryIndex() {
return (dir == InternalOptimizationDirection::Maximize) ? maxIndex : minIndex;
}
template<InternalOptimizationDirection dir>
inline uint64_t& getSecondaryIndex() {
return (dir == InternalOptimizationDirection::Maximize) ? minIndex : maxIndex;
}
std::vector<ValueType>& x;
std::vector<ValueType>& y;
std::vector<ValueType> xTmp, yTmp;
ValueType lowerBound, upperBound, decisionValue;
bool hasLowerBound, hasUpperBound, hasDecisionValue;
bool convergencePhase1;
bool decisionValueBlocks;
uint64_t firstIndexViolatingConvergence;
uint64_t minIndex, maxIndex;
bool relative;
ValueType precision;
IndexType numRows;
std::vector<ValueType> matrixValues;
std::vector<IndexType> matrixColumns;
std::vector<IndexType> rowIndications;
std::vector<uint_fast64_t> const* rowGroupIndices;
};
}
}
}
Loading…
Cancel
Save