15namespace gviinternal {
17template<
typename ValueType>
20template<
typename ValueType>
22 const auto& rowGroupIndices = matrix.getRowGroupIndices();
27 auto den = storm::utility::convertNumber<ValueType>(rowGroupIndices[i + 1] - rowGroupIndices[i]);
28 auto weightOnGroup = weights[i] / den;
31 for (
auto row = rowGroupIndices[i]; row < rowGroupIndices[i + 1]; ++row) {
32 for (
auto const& entry : matrix.getRow(row)) {
33 weights[entry.getColumn()] += weightOnGroup * entry.getValue();
39template<
typename ValueType>
43 storm::utility::vector::addScaledVector<ValueType, ValueType>(weights, lowerX, -1);
45 std::vector<ValueType> additiveWeights(weights.size(), 0);
46 for (
int i = 0; i < 50; i++) {
51 return static_cast<IndexType>(std::max_element(additiveWeights.begin(), additiveWeights.end()) - additiveWeights.begin());
54template<
typename ValueType>
56 auto itL = lower.begin();
57 const auto itLEnd = lower.end();
58 auto itU = upper.begin();
60 ValueType result = *itU - *itL;
61 for (++itL, ++itU; itL != itLEnd; ++itL, ++itU) {
62 ValueType length = *itU - *itL;
63 result = std::max(result, length);
68template<
typename ValueType>
71 auto itL = lower.begin();
72 const auto itLEnd = lower.end();
73 auto itU = upper.begin();
75 for (; itL != itLEnd; ++itL, ++itU) {
83template<
typename ValueType, storm::OptimizationDirection Dir>
89 GVIBackend(uint64_t guessedRowGroup) : guessedRowGroup(guessedRowGroup) {}
93 void firstRow(std::pair<ValueType, ValueType>&& value, [[maybe_unused]] uint64_t rowGroup, [[maybe_unused]] uint64_t row) {
94 xBest = std::move(value.first);
95 yBest = std::move(value.second);
98 void nextRow(std::pair<ValueType, ValueType>&& value, [[maybe_unused]] uint64_t rowGroup, [[maybe_unused]] uint64_t row) {
99 xBest &= std::move(value.first);
100 yBest &= std::move(value.second);
103 void applyUpdate(ValueType& xCurr, ValueType& yCurr, uint64_t rowGroup) {
104 if (guessedRowGroup && rowGroup == *guessedRowGroup) {
109 xCurr = std::max(xCurr, *xBest);
110 yCurr = std::min(yCurr, *yBest);
125 std::optional<uint64_t> guessedRowGroup;
128template<
typename ValueType,
bool TrivialRowGrouping>
129gviinternal::IndexType GuessingValueIterationHelper<ValueType, TrivialRowGrouping>::selectRowGroupToGuess(std::vector<ValueType>& lowerX,
130 std::vector<ValueType>& upperX) {
131 return iterationHelper.selectRowGroupToGuess(lowerX, upperX);
134template<
typename ValueType,
bool TrivialRowGrouping>
135template<OptimizationDirection Dir>
136void GuessingValueIterationHelper<ValueType, TrivialRowGrouping>::applyInPlace(std::vector<ValueType>& lowerX, std::vector<ValueType>& upperX,
137 const std::vector<ValueType>& b, GVIBackend<ValueType, Dir>& backend) {
138 std::pair<std::vector<ValueType>, std::vector<ValueType>> xy;
139 xy.first.swap(lowerX);
140 xy.second.swap(upperX);
141 viOperator->applyInPlace(xy, b, backend);
142 lowerX.swap(xy.first);
143 upperX.swap(xy.second);
146template<
typename ValueType,
bool TrivialRowGrouping>
147template<OptimizationDirection Dir>
148std::pair<gviinternal::VerifyResult, SolverStatus> GuessingValueIterationHelper<ValueType, TrivialRowGrouping>::tryVerify(
149 std::vector<ValueType>& lowerX, std::vector<ValueType>& upperX,
const std::vector<ValueType>& b, uint64_t& numIterations,
151 std::function<
SolverStatus(GVIData<ValueType>
const&)>
const& iterationCallback) {
154 guessLower[rowGroupToGuess] = guessUpper[rowGroupToGuess] = guessValue;
155 ValueType sumLengthBefore = 0, maxLengthBefore = 0;
156 GVIBackend<ValueType, Dir> guessingBackend(rowGroupToGuess);
157 GVIBackend<ValueType, Dir> iiBackend;
162 applyInPlace<Dir>(lowerX, upperX, b, iiBackend);
163 applyInPlace<Dir>(guessLower, guessUpper, b, guessingBackend);
164 auto guessedNewLower = guessingBackend.xGuessed;
165 auto guessedNewUpper = guessingBackend.yGuessed;
166 if (guessValue <= guessedNewLower) {
168 lowerX[rowGroupToGuess] = guessedNewLower;
171 if (guessedNewUpper <= guessValue) {
173 upperX[rowGroupToGuess] = guessedNewUpper;
177 auto sumLength = iterationHelper.getSumLength(guessLower, guessUpper);
178 if (sumLength == sumLengthBefore) {
182 sumLengthBefore = sumLength;
184 if (iterationHelper.getMaxLength(lowerX, upperX) < 2 * precision) {
188 if (iterationCallback) {
189 status = iterationCallback({lowerX, upperX, status});
195template<
typename ValueType,
bool TrivialRowGrouping>
197 std::vector<ValueType>& lowerX, std::vector<ValueType>& upperX,
const std::vector<ValueType>& b, uint64_t& numIterations, ValueType precision,
199 if (!dir.has_value() ||
minimize(dir.value()))
200 return solveEquations<OptimizationDirection::Minimize>(lowerX, upperX, b, numIterations, precision, iterationCallback);
202 return solveEquations<OptimizationDirection::Maximize>(lowerX, upperX, b, numIterations, precision, iterationCallback);
205template<
typename ValueType,
bool TrivialRowGrouping>
206template<OptimizationDirection Dir>
208 std::vector<ValueType>& lowerX, std::vector<ValueType>& upperX,
const std::vector<ValueType>& b, uint64_t& numIterations, ValueType precision,
211 auto status = doIterations<Dir>(lowerX, upperX, b, numIterations, lowerX.size(), precision, iterationCallback);
215 auto rowGroupToGuess = selectRowGroupToGuess(lowerX, upperX);
216 bool didVerify =
false;
218 for (
int den = 2; den < 30 && !didVerify; ++den) {
219 for (
int num = 1; num < den; num++) {
220 if (std::gcd(num, den) != 1)
222 auto guessValue = (lowerX[rowGroupToGuess] * num + upperX[rowGroupToGuess] * (den - num)) / den;
223 auto [verifyResult, verifyStatus] = tryVerify<Dir>(lowerX, upperX, b, numIterations, rowGroupToGuess, guessValue, precision, iterationCallback);
224 status = verifyStatus;
232 if (iterationHelper.getMaxLength(lowerX, upperX) < 2 * precision)
240 return doIterations<Dir>(lowerX, upperX, b, numIterations, {}, precision, iterationCallback);
243template<
typename ValueType,
bool TrivialRowGrouping>
246 : viOperator(viOperator), iterationHelper(matrix) {
250template<
typename ValueType,
bool TrivialRowGrouping>
251template<OptimizationDirection Dir>
253 std::vector<ValueType>& lowerX, std::vector<ValueType>& upperX,
const std::vector<ValueType>& b, uint64_t& numIterations,
254 std::optional<uint64_t> maxIterations,
const ValueType& precision, std::function<
SolverStatus(
GVIData<ValueType> const&)>
const& iterationCallback) {
255 ValueType sumLengthBefore = 0, maxLengthBefore = 0;
256 uint64_t localIterations = 0;
260 applyInPlace<Dir>(lowerX, upperX, b, iiBackend);
264 auto sumLength = iterationHelper.getSumLength(lowerX, upperX);
265 if (sumLengthBefore == sumLength) {
269 sumLengthBefore = sumLength;
271 if (iterationHelper.getMaxLength(lowerX, upperX) < 2 * precision) {
274 if (iterationCallback) {
275 status = iterationCallback(GVIData<ValueType>{lowerX, upperX, status});
279 if (iterationHelper.getMaxLength(lowerX, upperX) < 2 * precision)
284template class GuessingValueIterationHelper<double, true>;
285template class GuessingValueIterationHelper<double, false>;
286template class GuessingValueIterationHelper<storm::RationalNumber, true>;
287template class GuessingValueIterationHelper<storm::RationalNumber, false>;
void nextRow(std::pair< ValueType, ValueType > &&value, uint64_t rowGroup, uint64_t row)
bool constexpr converged() const
void applyUpdate(ValueType &xCurr, ValueType &yCurr, uint64_t rowGroup)
void firstRow(std::pair< ValueType, ValueType > &&value, uint64_t rowGroup, uint64_t row)
GVIBackend(uint64_t guessedRowGroup)
bool constexpr abort() const
void endOfIteration() const
Implements guessing value iteration.
GuessingValueIterationHelper(std::shared_ptr< ValueIterationOperator< ValueType, TrivialRowGrouping > > viOperator, const storage::SparseMatrix< ValueType > &matrix)
SolverStatus solveEquations(std::vector< ValueType > &lowerX, std::vector< ValueType > &upperX, const std::vector< ValueType > &b, uint64_t &numIterations, ValueType precision, std::optional< storm::solver::OptimizationDirection > dir, std::function< SolverStatus(GVIData< ValueType > const &)> const &iterationCallback)
This class represents the Value Iteration Operator (also known as Bellman operator).
ValueType getSumLength(std::vector< ValueType > &lower, std::vector< ValueType > &upper) const
ValueType getMaxLength(std::vector< ValueType > &lower, std::vector< ValueType > &upper) const
IndexType selectRowGroupToGuess(const std::vector< ValueType > &lowerX, const std::vector< ValueType > &upperX)
IterationHelper(storm::storage::SparseMatrix< ValueType > const &matrix_)
A class that holds a possibly non-square matrix in the compressed row storage format.
Stores and manages an extremal (maximal or minimal) value.
bool constexpr minimize(OptimizationDirection d)
void addScaledVector(std::vector< InValueType1 > &firstOperand, std::vector< InValueType2 > const &secondOperand, InValueType3 const &factor)
Computes x:= x + a*y, i.e., adds each element of the first vector and (the corresponding element of t...