50 template<
bool Backward = true>
98 template<
typename OperandType,
typename OffsetType,
typename BackendType>
99 bool apply(OperandType
const& operandIn, OperandType& operandOut, OffsetType
const& offsets, BackendType& backend)
const {
101 return applyRobust<OptimizationDirection::Maximize>(operandIn, operandOut, offsets, backend);
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) {
108 return apply<OperandType, OffsetType, BackendType, true, true, RobustDir>(operandOut, operandIn, offsets, backend);
110 return apply<OperandType, OffsetType, BackendType, false, true, RobustDir>(operandOut, operandIn, offsets, backend);
114 return apply<OperandType, OffsetType, BackendType, true, false, RobustDir>(operandOut, operandIn, offsets, backend);
116 return apply<OperandType, OffsetType, BackendType, false, false, RobustDir>(operandOut, operandIn, offsets, backend);
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);
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);
159 std::vector<SolutionType>&
allocateAuxiliaryVector(uint64_t size, std::optional<SolutionType>
const& initialValue = {});
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) {
183 if constexpr (std::is_same<BackendType, RobustSchedulerTrackingBackend<double, RobustDirection, TrivialRowGrouping>>::value) {
185 backend.processRobustRow(applyRow<RobustDirection>(matrixColumnIt, matrixValueIt, operandIn, offsets, groupIndex), groupIndex,
186 applyCache.robustOrder);
189 backend.firstRow(applyRow<RobustDirection>(matrixColumnIt, matrixValueIt, operandIn, offsets, groupIndex), groupIndex, groupIndex);
192 IndexType rowIndex = (*rowGroupIndices)[groupIndex];
193 if constexpr (SkipIgnoredRows) {
194 rowIndex += skipMultipleIgnoredRows(matrixColumnIt, matrixValueIt);
196 backend.firstRow(applyRow<RobustDirection>(matrixColumnIt, matrixValueIt, operandIn, offsets, rowIndex), groupIndex, rowIndex);
197 while (*matrixColumnIt < StartOfRowGroupIndicator) {
199 if (!SkipIgnoredRows || !skipIgnoredRow(matrixColumnIt, matrixValueIt)) {
200 backend.nextRow(applyRow<RobustDirection>(matrixColumnIt, matrixValueIt, operandIn, offsets, rowIndex), groupIndex, rowIndex);
204 if constexpr (isPair<OperandType>::value) {
205 backend.applyUpdate(operandOut.first[groupIndex], operandOut.second[groupIndex], groupIndex);
207 backend.applyUpdate(operandOut[groupIndex], groupIndex);
209 if (backend.abort()) {
210 return backend.converged();
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();
221 template<
typename OpT>
222 OpT initializeRowRes(std::vector<OpT>
const&, OpT
const& offsets, uint64_t offsetIndex)
const {
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];
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]};
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};
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();
248 return offsets[offsetIndex].lower();
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.");
257 return {offsets[offsetIndex].upper(), offsets[offsetIndex].upper()};
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.");
265 return {(*offsets.first)[offsetIndex], offsets.second};
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);
277 return applyRowStandard(matrixColumnIt, matrixValueIt, operand, offsets, offsetIndex);
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);
291 result += operand[*matrixColumnIt] * (*matrixValueIt);
298 template<OptimizationDirection RobustDirection>
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;
305 return a.first < b.first;
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)};
316 applyCache.robustOrder.clear();
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.");
325 result += operand[*matrixColumnIt] * lower;
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.");
339 result += operand[*matrixColumnIt] * lower;
341 remainingValue -= lower;
342 auto const diameter = matrixValueIt->upper() - lower;
344 applyCache.robustOrder.emplace_back(operand[*matrixColumnIt], std::make_pair(diameter, orderCounter));
351 static AuxCompare<RobustDirection> cmp;
352 std::sort(applyCache.robustOrder.begin(), applyCache.robustOrder.end(), cmp);
354 for (
auto const& pair : applyCache.robustOrder) {
355 auto availableMass = std::min(pair.second.first, remainingValue);
356 result += availableMass * pair.first;
357 remainingValue -= availableMass;
365 "Remaining value should be zero (all prob mass taken) or it should be a sad state, but is " << remainingValue);
370 template<
bool Backward>
372 if constexpr (Backward) {
373 return boost::adaptors::reverse(boost::irange(start, end));
375 return boost::irange(start, end);
380 uint64_t getSize(std::vector<T>
const& vec)
const {
384 template<
typename T1,
typename T2>
385 uint64_t getSize(std::pair<T1, T2>
const& pairOfVec)
const {
386 return pairOfVec.first.size();
390 struct isPair : std::false_type {};
392 template<
typename T1,
typename T2>
393 struct isPair<
std::pair<T1, T2>> : std::true_type {};
398 template<
bool Backward = true>
404 void moveToEndOfRow(std::vector<IndexType>::iterator& matrixColumnIt)
const;
409 bool skipIgnoredRow(std::vector<IndexType>::const_iterator& matrixColumnIt,
typename std::vector<ValueType>::const_iterator& matrixValueIt)
const;
414 uint64_t skipMultipleIgnoredRows(std::vector<IndexType>::const_iterator& matrixColumnIt,
415 typename std::vector<ValueType>::const_iterator& matrixValueIt)
const;
420 std::vector<ValueType> matrixValues;
426 std::vector<IndexType> matrixColumns;
431 std::vector<IndexType>
const* rowGroupIndices;
436 bool backwards{
true};
441 bool hasSkippedRows{
false};
446 std::vector<SolutionType> auxiliaryVector;
451 bool auxiliaryVectorUsedExternally{
false};
455 template<
typename ApplyValueType,
typename Dummy>
456 struct ApplyCache {};
458 template<
typename Dummy>
460 mutable std::vector<std::pair<SolutionType, std::pair<SolutionType, uint64_t>>> robustOrder;
461 storage::BitVector hasOnlyConstants;
467 ApplyCache<ValueType, int> applyCache;
472 static constexpr IndexType StartOfRowIndicator = 1ull << 63;
477 static constexpr IndexType StartOfRowGroupIndicator = StartOfRowIndicator + (1ull << 62);
482 static constexpr IndexType SkipNumEntriesMask = ~StartOfRowGroupIndicator;