Storm 1.10.0.1
A Modern Probabilistic Model Checker
Loading...
Searching...
No Matches
ValueIterationOperator.h
Go to the documentation of this file.
1#pragma once
2#include <algorithm>
3#include <functional>
4#include <optional>
5#include <utility>
6#include <vector>
7
8#include <boost/range/adaptor/reversed.hpp>
9#include <boost/range/irange.hpp>
10
18#include "storm/utility/vector.h" // TODO
19
20namespace storm {
21class Environment;
22
23namespace storage {
24template<typename T>
25class SparseMatrix;
26}
27
28namespace solver::helper {
29
38template<typename ValueType, bool TrivialRowGrouping, typename SolutionType>
40 public:
42
50 template<bool Backward = true>
51 void setMatrix(storm::storage::SparseMatrix<ValueType> const& matrix, std::vector<IndexType> const* rowGroupIndices = nullptr);
52
59 void setMatrixForwards(storm::storage::SparseMatrix<ValueType> const& matrix, std::vector<IndexType> const* rowGroupIndices = nullptr);
60
67 void setMatrixBackwards(storm::storage::SparseMatrix<ValueType> const& matrix, std::vector<IndexType> const* rowGroupIndices = nullptr);
68
98 template<typename OperandType, typename OffsetType, typename BackendType>
99 bool apply(OperandType const& operandIn, OperandType& operandOut, OffsetType const& offsets, BackendType& backend) const {
100 // We assume in the normal apply case that we can just maximize.
101 return applyRobust<OptimizationDirection::Maximize>(operandIn, operandOut, offsets, backend);
102 }
103
104 template<OptimizationDirection RobustDir, typename OperandType, typename OffsetType, typename BackendType>
105 bool applyRobust(OperandType const& operandIn, OperandType& operandOut, OffsetType const& offsets, BackendType& backend) const {
106 if (hasSkippedRows) {
107 if (backwards) {
108 return apply<OperandType, OffsetType, BackendType, true, true, RobustDir>(operandOut, operandIn, offsets, backend);
109 } else {
110 return apply<OperandType, OffsetType, BackendType, false, true, RobustDir>(operandOut, operandIn, offsets, backend);
111 }
112 } else {
113 if (backwards) {
114 return apply<OperandType, OffsetType, BackendType, true, false, RobustDir>(operandOut, operandIn, offsets, backend);
115 } else {
116 return apply<OperandType, OffsetType, BackendType, false, false, RobustDir>(operandOut, operandIn, offsets, backend);
117 }
118 }
119 }
120
124 template<typename OperandType, typename OffsetType, typename BackendType>
125 bool applyInPlace(OperandType& operand, OffsetType const& offsets, BackendType& backend) const {
126 return apply(operand, operand, offsets, backend);
127 }
128
129 template<OptimizationDirection RobustDir, typename OperandType, typename OffsetType, typename BackendType>
130 bool applyInPlaceRobust(OperandType& operand, OffsetType const& offsets, BackendType& backend) const {
131 return applyRobust<RobustDir>(operand, operand, offsets, backend);
132 }
133
141 void setIgnoredRows(bool useLocalRowIndices, std::function<bool(IndexType, IndexType)> const& ignore);
142
146 void unsetIgnoredRows();
147
151 std::vector<IndexType> const& getRowGroupIndices() const;
152
159 std::vector<SolutionType>& allocateAuxiliaryVector(uint64_t size, std::optional<SolutionType> const& initialValue = {});
160
164 void freeAuxiliaryVector();
165
166 private:
171 template<typename OperandType, typename OffsetType, typename BackendType, bool Backward, bool SkipIgnoredRows, OptimizationDirection RobustDirection>
172 bool apply(OperandType& operandOut, OperandType const& operandIn, OffsetType const& offsets, BackendType& backend) const {
173 auto const outSize = TrivialRowGrouping ? getSize(operandOut) : rowGroupIndices->size() - 1;
174 STORM_LOG_ASSERT(TrivialRowGrouping || getSize(operandOut) >= outSize, "Dimension mismatch");
175 backend.startNewIteration();
176 auto matrixValueIt = matrixValues.cbegin();
177 auto matrixColumnIt = matrixColumns.cbegin();
178 for (auto groupIndex : indexRange<Backward>(0, outSize)) {
179 STORM_LOG_ASSERT(matrixColumnIt != matrixColumns.end(), "VI Operator in invalid state.");
180 STORM_LOG_ASSERT(*matrixColumnIt >= StartOfRowIndicator, "VI Operator in invalid state.");
181 if constexpr (TrivialRowGrouping) {
182 // Ugly special case
183 if constexpr (std::is_same<BackendType, RobustSchedulerTrackingBackend<double, RobustDirection, TrivialRowGrouping>>::value) {
184 // Intentionally different method name
185 backend.processRobustRow(applyRow<RobustDirection>(matrixColumnIt, matrixValueIt, operandIn, offsets, groupIndex), groupIndex,
186 applyCache.robustOrder);
187 } else {
188 // Generic nextRow interface
189 backend.firstRow(applyRow<RobustDirection>(matrixColumnIt, matrixValueIt, operandIn, offsets, groupIndex), groupIndex, groupIndex);
190 }
191 } else {
192 IndexType rowIndex = (*rowGroupIndices)[groupIndex];
193 if constexpr (SkipIgnoredRows) {
194 rowIndex += skipMultipleIgnoredRows(matrixColumnIt, matrixValueIt);
195 }
196 backend.firstRow(applyRow<RobustDirection>(matrixColumnIt, matrixValueIt, operandIn, offsets, rowIndex), groupIndex, rowIndex);
197 while (*matrixColumnIt < StartOfRowGroupIndicator) {
198 ++rowIndex;
199 if (!SkipIgnoredRows || !skipIgnoredRow(matrixColumnIt, matrixValueIt)) {
200 backend.nextRow(applyRow<RobustDirection>(matrixColumnIt, matrixValueIt, operandIn, offsets, rowIndex), groupIndex, rowIndex);
201 }
202 }
203 }
204 if constexpr (isPair<OperandType>::value) {
205 backend.applyUpdate(operandOut.first[groupIndex], operandOut.second[groupIndex], groupIndex);
206 } else {
207 backend.applyUpdate(operandOut[groupIndex], groupIndex);
208 }
209 if (backend.abort()) {
210 return backend.converged();
211 }
212 }
213 STORM_LOG_ASSERT(matrixColumnIt + 1 == matrixColumns.cend(), "Unexpected position of matrix column iterator.");
214 STORM_LOG_ASSERT(matrixValueIt == matrixValues.cend(), "Unexpected position of matrix column iterator.");
215 backend.endOfIteration();
216 return backend.converged();
217 }
218
219 // Auxiliary methods to deal with various OperandTypes and OffsetTypes
220
221 template<typename OpT>
222 OpT initializeRowRes(std::vector<OpT> const&, OpT const& offsets, uint64_t offsetIndex) const {
223 return offsets;
224 }
225
226 template<typename OpT, typename OffT>
227 OpT initializeRowRes(std::vector<OpT> const&, std::vector<OffT> const& offsets, uint64_t offsetIndex) const {
228 return offsets[offsetIndex];
229 }
230
231 template<typename OpT1, typename OpT2, typename OffT>
232 std::pair<OpT1, OpT2> initializeRowRes(std::pair<std::vector<OpT1>, std::vector<OpT2>> const&, std::vector<OffT> const& offsets,
233 uint64_t offsetIndex) const {
234 return {offsets[offsetIndex], offsets[offsetIndex]};
235 }
236
237 template<typename OpT1, typename OpT2, typename OffT1, typename OffT2>
238 std::pair<OpT1, OpT2> initializeRowRes(std::pair<std::vector<OpT1>, std::vector<OpT2>> const&, std::pair<std::vector<OffT1> const*, OffT2> const& offsets,
239 uint64_t offsetIndex) const {
240 return {(*offsets.first)[offsetIndex], offsets.second};
241 }
242
243 template<OptimizationDirection RobustDirection, typename OpT, typename OffT>
244 OpT robustInitializeRowRes(std::vector<OpT> const&, std::vector<OffT> const& offsets, uint64_t offsetIndex) const {
245 if constexpr (RobustDirection == OptimizationDirection::Maximize) {
246 return offsets[offsetIndex].upper();
247 } else {
248 return offsets[offsetIndex].lower();
249 }
250 }
251
252 template<OptimizationDirection RobustDirection, typename OpT1, typename OpT2, typename OffT>
253 std::pair<OpT1, OpT2> robustInitializeRowRes(std::pair<std::vector<OpT1>, std::vector<OpT2>> const&, std::vector<OffT> const& offsets,
254 uint64_t offsetIndex) const {
255 STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "Value Iteration is not implemented with pairs and interval-models.");
256
257 return {offsets[offsetIndex].upper(), offsets[offsetIndex].upper()};
258 }
259
260 template<OptimizationDirection RobustDirection, typename OpT1, typename OpT2, typename OffT1, typename OffT2>
261 std::pair<OpT1, OpT2> robustInitializeRowRes(std::pair<std::vector<OpT1>, std::vector<OpT2>> const&,
262 std::pair<std::vector<OffT1> const*, OffT2> const& offsets, uint64_t offsetIndex) const {
263 STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "Value Iteration is not implemented with pairs and interval-models.");
264
265 return {(*offsets.first)[offsetIndex], offsets.second};
266 }
267
271 template<OptimizationDirection RobustDirection, typename OperandType, typename OffsetType>
272 auto applyRow(std::vector<IndexType>::const_iterator& matrixColumnIt, typename std::vector<ValueType>::const_iterator& matrixValueIt,
273 OperandType const& operand, OffsetType const& offsets, uint64_t offsetIndex) const {
274 if constexpr (std::is_same_v<ValueType, storm::Interval>) {
275 return applyRowRobust<RobustDirection>(matrixColumnIt, matrixValueIt, operand, offsets, offsetIndex);
276 } else {
277 return applyRowStandard(matrixColumnIt, matrixValueIt, operand, offsets, offsetIndex);
278 }
279 }
280
281 template<typename OperandType, typename OffsetType>
282 auto applyRowStandard(std::vector<IndexType>::const_iterator& matrixColumnIt, typename std::vector<ValueType>::const_iterator& matrixValueIt,
283 OperandType const& operand, OffsetType const& offsets, uint64_t offsetIndex) const {
284 STORM_LOG_ASSERT(*matrixColumnIt >= StartOfRowIndicator, "VI Operator in invalid state.");
285 auto result{initializeRowRes(operand, offsets, offsetIndex)};
286 for (++matrixColumnIt; *matrixColumnIt < StartOfRowIndicator; ++matrixColumnIt, ++matrixValueIt) {
287 if constexpr (isPair<OperandType>::value) {
288 result.first += operand.first[*matrixColumnIt] * (*matrixValueIt);
289 result.second += operand.second[*matrixColumnIt] * (*matrixValueIt);
290 } else {
291 result += operand[*matrixColumnIt] * (*matrixValueIt);
292 }
293 }
294 return result;
295 }
296
297 // Aux function for applyRowRobust
298 template<OptimizationDirection RobustDirection>
299 struct AuxCompare {
300 bool operator()(const std::pair<SolutionType, std::pair<SolutionType, uint64_t>>& a,
301 const std::pair<SolutionType, std::pair<SolutionType, uint64_t>>& b) const {
302 if constexpr (RobustDirection == OptimizationDirection::Maximize) {
303 return a.first > b.first;
304 } else {
305 return a.first < b.first;
306 }
307 }
308 };
309
310 template<OptimizationDirection RobustDirection, typename OperandType, typename OffsetType>
311 auto applyRowRobust(std::vector<IndexType>::const_iterator& matrixColumnIt, typename std::vector<ValueType>::const_iterator& matrixValueIt,
312 OperandType const& operand, OffsetType const& offsets, uint64_t offsetIndex) const {
313 STORM_LOG_ASSERT(*matrixColumnIt >= StartOfRowIndicator, "VI Operator in invalid state.");
314 auto result{robustInitializeRowRes<RobustDirection>(operand, offsets, offsetIndex)};
315
316 applyCache.robustOrder.clear();
317
318 if (applyCache.hasOnlyConstants.size() > 0 && applyCache.hasOnlyConstants.get(offsetIndex)) {
319 for (++matrixColumnIt; *matrixColumnIt < StartOfRowIndicator; ++matrixColumnIt, ++matrixValueIt) {
320 auto const lower = matrixValueIt->lower();
321 if constexpr (isPair<OperandType>::value) {
322 STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "Value Iteration is not implemented with pairs and interval-models.");
323 // Notice the unclear semantics here in terms of how to order things.
324 } else {
325 result += operand[*matrixColumnIt] * lower;
326 }
327 }
328 return result;
329 }
330
331 SolutionType remainingValue{storm::utility::one<SolutionType>()};
332 uint64_t orderCounter = 0;
333 for (++matrixColumnIt; *matrixColumnIt < StartOfRowIndicator; ++matrixColumnIt, ++matrixValueIt, ++orderCounter) {
334 auto const lower = matrixValueIt->lower();
335 if constexpr (isPair<OperandType>::value) {
336 STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "Value Iteration is not implemented with pairs and interval-models.");
337 // Notice the unclear semantics here in terms of how to order things.
338 } else {
339 result += operand[*matrixColumnIt] * lower;
340 }
341 remainingValue -= lower;
342 auto const diameter = matrixValueIt->upper() - lower;
343 if (!storm::utility::isZero(diameter)) {
344 applyCache.robustOrder.emplace_back(operand[*matrixColumnIt], std::make_pair(diameter, orderCounter));
345 }
346 }
347 if (storm::utility::isZero(remainingValue)) {
348 return result;
349 }
350
351 static AuxCompare<RobustDirection> cmp;
352 std::sort(applyCache.robustOrder.begin(), applyCache.robustOrder.end(), cmp);
353
354 for (auto const& pair : applyCache.robustOrder) {
355 auto availableMass = std::min(pair.second.first, remainingValue);
356 result += availableMass * pair.first;
357 remainingValue -= availableMass;
358 if (storm::utility::isZero(remainingValue)) {
359 return result;
360 }
361 }
363 // sad states allowed (they're having a bummer summer)
364 (storm::utility::isOne(remainingValue) && applyCache.robustOrder.size() == 0),
365 "Remaining value should be zero (all prob mass taken) or it should be a sad state, but is " << remainingValue);
366 return result;
367 }
368
369 // Auxiliary helpers used for metaprogramming
370 template<bool Backward>
371 auto indexRange(IndexType start, IndexType end) const {
372 if constexpr (Backward) {
373 return boost::adaptors::reverse(boost::irange(start, end));
374 } else {
375 return boost::irange(start, end);
376 }
377 }
378
379 template<typename T>
380 uint64_t getSize(std::vector<T> const& vec) const {
381 return vec.size();
382 }
383
384 template<typename T1, typename T2>
385 uint64_t getSize(std::pair<T1, T2> const& pairOfVec) const {
386 return pairOfVec.first.size();
387 }
388
389 template<typename>
390 struct isPair : std::false_type {};
391
392 template<typename T1, typename T2>
393 struct isPair<std::pair<T1, T2>> : std::true_type {};
394
398 template<bool Backward = true>
399 void setIgnoredRows(bool useLocalRowIndices, std::function<bool(IndexType, IndexType)> const& ignore);
400
404 void moveToEndOfRow(std::vector<IndexType>::iterator& matrixColumnIt) const;
405
409 bool skipIgnoredRow(std::vector<IndexType>::const_iterator& matrixColumnIt, typename std::vector<ValueType>::const_iterator& matrixValueIt) const;
410
414 uint64_t skipMultipleIgnoredRows(std::vector<IndexType>::const_iterator& matrixColumnIt,
415 typename std::vector<ValueType>::const_iterator& matrixValueIt) const;
416
420 std::vector<ValueType> matrixValues;
421
426 std::vector<IndexType> matrixColumns;
427
431 std::vector<IndexType> const* rowGroupIndices;
432
436 bool backwards{true};
437
441 bool hasSkippedRows{false};
442
446 std::vector<SolutionType> auxiliaryVector;
447
451 bool auxiliaryVectorUsedExternally{false};
452
453 // Due to a GCC bug we have to add this dummy template type here
454 // https://stackoverflow.com/questions/49707184/explicit-specialization-in-non-namespace-scope-does-not-compile-in-gcc
455 template<typename ApplyValueType, typename Dummy>
456 struct ApplyCache {};
457
458 template<typename Dummy>
459 struct ApplyCache<storm::Interval, Dummy> {
460 mutable std::vector<std::pair<SolutionType, std::pair<SolutionType, uint64_t>>> robustOrder;
461 storage::BitVector hasOnlyConstants;
462 };
463
467 ApplyCache<ValueType, int> applyCache;
468
472 static constexpr IndexType StartOfRowIndicator = 1ull << 63; // 10000..0
473
477 static constexpr IndexType StartOfRowGroupIndicator = StartOfRowIndicator + (1ull << 62); // 11000..0
478
482 static constexpr IndexType SkipNumEntriesMask = ~StartOfRowGroupIndicator; // 00111..1
483};
484
485} // namespace solver::helper
486} // namespace storm
This class represents the Value Iteration Operator (also known as Bellman operator).
std::vector< SolutionType > & allocateAuxiliaryVector(uint64_t size, std::optional< SolutionType > const &initialValue={})
Allocates additional storage that can be used e.g.
void setIgnoredRows(bool useLocalRowIndices, std::function< bool(IndexType, IndexType)> const &ignore)
Sets rows that will be skipped when applying the operator.
bool applyRobust(OperandType const &operandIn, OperandType &operandOut, OffsetType const &offsets, BackendType &backend) const
std::vector< IndexType > const & getRowGroupIndices() const
void setMatrix(storm::storage::SparseMatrix< ValueType > const &matrix, std::vector< IndexType > const *rowGroupIndices=nullptr)
Initializes this operator with the given data.
void freeAuxiliaryVector()
Clears the auxiliary vector, invalidating any references to it.
bool applyInPlaceRobust(OperandType &operand, OffsetType const &offsets, BackendType &backend) const
bool apply(OperandType const &operandIn, OperandType &operandOut, OffsetType const &offsets, BackendType &backend) const
Applies the operator with the given operands, offsets, and backend.
bool applyInPlace(OperandType &operand, OffsetType const &offsets, BackendType &backend) const
Same as apply but with operandOut==operandIn.
void setMatrixForwards(storm::storage::SparseMatrix< ValueType > const &matrix, std::vector< IndexType > const *rowGroupIndices=nullptr)
Initializes this operator with the given data for forward iterations (starting with the smallest row ...
void setMatrixBackwards(storm::storage::SparseMatrix< ValueType > const &matrix, std::vector< IndexType > const *rowGroupIndices=nullptr)
Initializes this operator with the given data for backward iterations (starting with the largest row ...
A class that holds a possibly non-square matrix in the compressed row storage format.
#define STORM_LOG_ASSERT(cond, message)
Definition macros.h:11
#define STORM_LOG_THROW(cond, exception, message)
Definition macros.h:30
bool isOne(ValueType const &a)
Definition constants.cpp:36
bool isAlmostZero(ValueType const &a)
Definition constants.cpp:56
bool isZero(ValueType const &a)
Definition constants.cpp:41
LabParser.cpp.
Definition cli.cpp:18
carl::Interval< double > Interval