Storm 1.10.0.1
A Modern Probabilistic Model Checker
Loading...
Searching...
No Matches
GuessingValueIterationHelper.cpp
Go to the documentation of this file.
2
4
8
10
12
13namespace storm::solver::helper {
14
15namespace gviinternal {
16
17template<typename ValueType>
19
20template<typename ValueType>
22 const auto& rowGroupIndices = matrix.getRowGroupIndices();
23 IndexType i = weights.size();
24 while (i > 0) {
25 --i;
26
27 auto den = storm::utility::convertNumber<ValueType>(rowGroupIndices[i + 1] - rowGroupIndices[i]);
28 auto weightOnGroup = weights[i] / den;
29 weights[i] = 0;
30
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();
34 }
35 }
36 }
37}
38
39template<typename ValueType>
40IndexType IterationHelper<ValueType>::selectRowGroupToGuess(const std::vector<ValueType>& lowerX, const std::vector<ValueType>& upperX) {
41 weights = upperX;
42 // upper - lower
43 storm::utility::vector::addScaledVector<ValueType, ValueType>(weights, lowerX, -1);
44
45 std::vector<ValueType> additiveWeights(weights.size(), 0);
46 for (int i = 0; i < 50; i++) {
47 swipeWeights();
48 storm::utility::vector::addScaledVector(additiveWeights, weights, 1);
49 }
50
51 return static_cast<IndexType>(std::max_element(additiveWeights.begin(), additiveWeights.end()) - additiveWeights.begin());
52}
53
54template<typename ValueType>
55ValueType IterationHelper<ValueType>::getMaxLength(std::vector<ValueType>& lower, std::vector<ValueType>& upper) const {
56 auto itL = lower.begin();
57 const auto itLEnd = lower.end();
58 auto itU = upper.begin();
59
60 ValueType result = *itU - *itL;
61 for (++itL, ++itU; itL != itLEnd; ++itL, ++itU) {
62 ValueType length = *itU - *itL;
63 result = std::max(result, length);
64 }
65 return result;
66}
67
68template<typename ValueType>
69ValueType IterationHelper<ValueType>::getSumLength(std::vector<ValueType>& lower, std::vector<ValueType>& upper) const {
70 ValueType sum = 0;
71 auto itL = lower.begin();
72 const auto itLEnd = lower.end();
73 auto itU = upper.begin();
74
75 for (; itL != itLEnd; ++itL, ++itU) {
76 sum += *itU - *itL;
77 }
78 return sum;
79}
80
81} // namespace gviinternal
82
83template<typename ValueType, storm::OptimizationDirection Dir>
85 public:
86 ValueType xGuessed, yGuessed;
87
89 GVIBackend(uint64_t guessedRowGroup) : guessedRowGroup(guessedRowGroup) {}
90
92
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);
96 }
97
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);
101 }
102
103 void applyUpdate(ValueType& xCurr, ValueType& yCurr, uint64_t rowGroup) {
104 if (guessedRowGroup && rowGroup == *guessedRowGroup) {
105 xGuessed = *xBest;
106 yGuessed = *yBest;
107 return;
108 }
109 xCurr = std::max(xCurr, *xBest);
110 yCurr = std::min(yCurr, *yBest);
111 }
112
113 void endOfIteration() const {}
114
115 bool constexpr converged() const {
116 return false;
117 }
118
119 bool constexpr abort() const {
120 return false;
121 }
122
123 private:
125 std::optional<uint64_t> guessedRowGroup;
126};
127
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);
132}
133
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);
144}
145
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,
150 gviinternal::IndexType rowGroupToGuess, ValueType guessValue, ValueType precision,
151 std::function<SolverStatus(GVIData<ValueType> const&)> const& iterationCallback) {
152 guessLower = lowerX;
153 guessUpper = upperX;
154 guessLower[rowGroupToGuess] = guessUpper[rowGroupToGuess] = guessValue;
155 ValueType sumLengthBefore = 0, maxLengthBefore = 0;
156 GVIBackend<ValueType, Dir> guessingBackend(rowGroupToGuess);
157 GVIBackend<ValueType, Dir> iiBackend;
158
159 auto status = SolverStatus::InProgress;
160 while (status == SolverStatus::InProgress) {
161 ++numIterations;
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) {
167 lowerX = guessLower;
168 lowerX[rowGroupToGuess] = guessedNewLower;
170 }
171 if (guessedNewUpper <= guessValue) {
172 upperX = guessUpper;
173 upperX[rowGroupToGuess] = guessedNewUpper;
175 }
176
177 auto sumLength = iterationHelper.getSumLength(guessLower, guessUpper);
178 if (sumLength == sumLengthBefore) {
179 // nothing changed. abort the verification
181 }
182 sumLengthBefore = sumLength;
183
184 if (iterationHelper.getMaxLength(lowerX, upperX) < 2 * precision) {
186 }
187
188 if (iterationCallback) {
189 status = iterationCallback({lowerX, upperX, status});
190 }
191 }
193}
194
195template<typename ValueType, bool TrivialRowGrouping>
197 std::vector<ValueType>& lowerX, std::vector<ValueType>& upperX, const std::vector<ValueType>& b, uint64_t& numIterations, ValueType precision,
198 std::optional<storm::solver::OptimizationDirection> dir, std::function<SolverStatus(GVIData<ValueType> const&)> const& iterationCallback) {
199 if (!dir.has_value() || minimize(dir.value()))
200 return solveEquations<OptimizationDirection::Minimize>(lowerX, upperX, b, numIterations, precision, iterationCallback);
201 else
202 return solveEquations<OptimizationDirection::Maximize>(lowerX, upperX, b, numIterations, precision, iterationCallback);
203}
204
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,
209 std::function<SolverStatus(GVIData<ValueType> const&)> const& iterationCallback) {
210 // do n iterations first
211 auto status = doIterations<Dir>(lowerX, upperX, b, numIterations, lowerX.size(), precision, iterationCallback);
212 if (status != SolverStatus::InProgress)
213 return status;
214 while (status == SolverStatus::InProgress) {
215 auto rowGroupToGuess = selectRowGroupToGuess(lowerX, upperX);
216 bool didVerify = false;
217 // try verification using different fractions of the interval
218 for (int den = 2; den < 30 && !didVerify; ++den) {
219 for (int num = 1; num < den; num++) {
220 if (std::gcd(num, den) != 1)
221 continue;
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;
225 if (verifyResult == gviinternal::VerifyResult::Verified || verifyResult == gviinternal::VerifyResult::Converged) {
226 didVerify = true;
227 break;
228 }
229 }
230 }
231 if (didVerify) {
232 if (iterationHelper.getMaxLength(lowerX, upperX) < 2 * precision)
234 } else {
235 break;
236 }
237 }
238
239 // verification has failed. iterate until convergence
240 return doIterations<Dir>(lowerX, upperX, b, numIterations, {}, precision, iterationCallback);
241}
242
243template<typename ValueType, bool TrivialRowGrouping>
246 : viOperator(viOperator), iterationHelper(matrix) {
247 // Intentionally left empty.
248}
249
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;
258 auto status = SolverStatus::InProgress;
259 while (status == SolverStatus::InProgress && (!maxIterations.has_value() || localIterations < maxIterations)) {
260 applyInPlace<Dir>(lowerX, upperX, b, iiBackend);
261 ++localIterations;
262 ++numIterations;
263
264 auto sumLength = iterationHelper.getSumLength(lowerX, upperX);
265 if (sumLengthBefore == sumLength) {
266 // nothing changed. abort the iterations.
268 }
269 sumLengthBefore = sumLength;
270
271 if (iterationHelper.getMaxLength(lowerX, upperX) < 2 * precision) {
273 }
274 if (iterationCallback) {
275 status = iterationCallback(GVIData<ValueType>{lowerX, upperX, status});
276 }
277 }
278
279 if (iterationHelper.getMaxLength(lowerX, upperX) < 2 * precision)
281 return status;
282}
283
284template class GuessingValueIterationHelper<double, true>;
285template class GuessingValueIterationHelper<double, false>;
286template class GuessingValueIterationHelper<storm::RationalNumber, true>;
287template class GuessingValueIterationHelper<storm::RationalNumber, false>;
288} // namespace storm::solver::helper
void nextRow(std::pair< ValueType, ValueType > &&value, uint64_t rowGroup, uint64_t row)
void applyUpdate(ValueType &xCurr, ValueType &yCurr, uint64_t rowGroup)
void firstRow(std::pair< ValueType, ValueType > &&value, uint64_t rowGroup, uint64_t row)
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.
Definition Extremum.h:15
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...
Definition vector.h:464