33template<
typename ValueType,
typename SolutionType>
36 : linearEquationSolverFactory(
std::move(linearEquationSolverFactory)) {
40template<
typename ValueType,
typename SolutionType>
47template<
typename ValueType,
typename SolutionType>
54template<
typename ValueType,
typename SolutionType>
59 if (isExactMode && method != MinMaxMethod::PolicyIteration && method != MinMaxMethod::RationalSearch && method != MinMaxMethod::ViToPi) {
62 "Selecting 'Policy iteration' as the solution technique to guarantee exact results. If you want to override this, please explicitly specify a "
64 method = MinMaxMethod::PolicyIteration;
66 STORM_LOG_WARN(
"The selected solution method " <<
toString(method) <<
" does not guarantee exact results.");
68 }
else if (env.
solver().
isForceSoundness() && method != MinMaxMethod::SoundValueIteration && method != MinMaxMethod::IntervalIteration &&
69 method != MinMaxMethod::PolicyIteration && method != MinMaxMethod::RationalSearch && method != MinMaxMethod::OptimisticValueIteration &&
70 method != MinMaxMethod::GuessingValueIteration) {
72 method = MinMaxMethod::OptimisticValueIteration;
76 <<
"' as the solution technique to guarantee sound results. If you want to override this, please explicitly specify a different method.");
78 STORM_LOG_WARN(
"The selected solution method does not guarantee sound results.");
81 STORM_LOG_THROW(method == MinMaxMethod::ValueIteration || method == MinMaxMethod::PolicyIteration || method == MinMaxMethod::RationalSearch ||
82 method == MinMaxMethod::SoundValueIteration || method == MinMaxMethod::IntervalIteration ||
83 method == MinMaxMethod::OptimisticValueIteration || method == MinMaxMethod::GuessingValueIteration || method == MinMaxMethod::ViToPi,
84 storm::exceptions::InvalidEnvironmentException,
"This solver does not support the selected method '" <<
toString(method) <<
"'.");
88template<
typename ValueType,
typename SolutionType>
90 std::vector<SolutionType>& x, std::vector<ValueType>
const& b)
const {
93 case MinMaxMethod::ValueIteration:
94 result = solveEquationsValueIteration(env, dir, x, b);
96 case MinMaxMethod::OptimisticValueIteration:
97 result = solveEquationsOptimisticValueIteration(env, dir, x, b);
99 case MinMaxMethod::GuessingValueIteration:
100 result = solveEquationsGuessingValueIteration(env, dir, x, b);
102 case MinMaxMethod::PolicyIteration:
103 result = solveEquationsPolicyIteration(env, dir, x, b);
105 case MinMaxMethod::RationalSearch:
106 result = solveEquationsRationalSearch(env, dir, x, b);
108 case MinMaxMethod::IntervalIteration:
109 result = solveEquationsIntervalIteration(env, dir, x, b);
111 case MinMaxMethod::SoundValueIteration:
112 result = solveEquationsSoundValueIteration(env, dir, x, b);
114 case MinMaxMethod::ViToPi:
115 result = solveEquationsViToPi(env, dir, x, b);
118 STORM_LOG_THROW(
false, storm::exceptions::InvalidEnvironmentException,
"This solver does not implement the selected solution method");
124template<
typename ValueType,
typename SolutionType>
126 if (!viOperatorTriv && !viOperatorNontriv) {
127 if (this->A->hasTrivialRowGrouping()) {
129 viOperatorTriv = std::make_shared<helper::ValueIterationOperator<ValueType, true, SolutionType>>();
130 viOperatorTriv->setMatrixBackwards(*this->A);
131 if constexpr (!std::is_same_v<ValueType, storm::Interval>) {
136 viOperatorNontriv = std::make_shared<helper::ValueIterationOperator<ValueType, false, SolutionType>>();
137 viOperatorNontriv->setMatrixBackwards(*this->A);
141 viOperatorNontriv = std::make_shared<helper::ValueIterationOperator<ValueType, false, SolutionType>>();
142 viOperatorNontriv->setMatrixBackwards(*this->A);
145 if (this->choiceFixedForRowGroup) {
147 assert(this->initialScheduler);
148 auto callback = [&](uint64_t groupIndex, uint64_t localRowIndex) {
149 return this->choiceFixedForRowGroup->get(groupIndex) && this->initialScheduler->at(groupIndex) != localRowIndex;
151 if (viOperatorTriv) {
152 viOperatorTriv->setIgnoredRows(
true, callback);
154 if (viOperatorNontriv) {
155 viOperatorNontriv->setIgnoredRows(
true, callback);
160template<
typename ValueType,
typename SolutionType>
161void IterativeMinMaxLinearEquationSolver<ValueType, SolutionType>::extractScheduler(std::vector<SolutionType>& x, std::vector<ValueType>
const& b,
163 if (std::is_same_v<ValueType, storm::Interval> && this->A->hasTrivialRowGrouping()) {
165 if (!this->robustSchedulerIndex) {
166 this->robustSchedulerIndex = std::vector<uint64_t>(x.size(), 0);
168 this->robustSchedulerIndex->resize(x.size(), 0);
170 uint64_t numSchedulerChoices = 0;
171 for (uint64_t row = 0; row < this->A->getRowCount(); ++row) {
172 this->robustSchedulerIndex->at(row) = numSchedulerChoices;
173 numSchedulerChoices += this->A->getRow(row).getNumberOfEntries();
176 if (!this->schedulerChoices) {
177 this->schedulerChoices = std::vector<uint64_t>(numSchedulerChoices, 0);
179 this->schedulerChoices->resize(numSchedulerChoices, 0);
183 if (!this->schedulerChoices) {
184 this->schedulerChoices = std::vector<uint64_t>(x.size(), 0);
186 this->schedulerChoices->resize(x.size(), 0);
192 "Expected VI operator to be initialized for scheduler extraction. Initializing now, but this is inefficient.");
193 if (!viOperatorTriv && !viOperatorNontriv) {
196 if (viOperatorTriv) {
197 if constexpr (std::is_same<ValueType, storm::Interval>() && std::is_same<SolutionType, double>()) {
199 schedHelper.computeScheduler(x, b, dir, *this->schedulerChoices, robust, updateX ? &x : nullptr, this->robustSchedulerIndex);
201 STORM_LOG_ERROR(
"SchedulerTrackingHelper not implemented for this setting (trivial row grouping but not Interval->double).");
204 if (viOperatorNontriv) {
206 schedHelper.computeScheduler(x, b, dir, *this->schedulerChoices, robust, updateX ? &x : nullptr, this->robustSchedulerIndex);
210template<
typename ValueType,
typename SolutionType>
211bool IterativeMinMaxLinearEquationSolver<ValueType, SolutionType>::solveInducedEquationSystem(
212 Environment
const& env, std::unique_ptr<LinearEquationSolver<SolutionType>>& linearEquationSolver, std::vector<uint64_t>
const& scheduler,
213 std::vector<SolutionType>& x, std::vector<ValueType>& subB, std::vector<ValueType>
const& originalB,
OptimizationDirection dir)
const {
214 if constexpr (std::is_same_v<ValueType, storm::Interval>) {
215 if constexpr (std::is_same_v<SolutionType, storm::Interval>) {
217 "We did not implement solving induced equation systems for interval-based models outside of robust VI.");
221 STORM_LOG_ASSERT(subB.size() == x.size(),
"Sizes of subB and x do not coincide.");
222 STORM_LOG_ASSERT(this->linearEquationSolverFactory !=
nullptr,
"Wrong constructor was called.");
229 auto schedulerIterator = scheduler.begin();
230 for (uint64_t rowIndex = 0; rowIndex < this->A->getRowCount(); rowIndex++) {
231 auto const& row = this->A->getRow(rowIndex);
233 static std::vector<SolutionType> tmp;
236 SolutionType probLeft = storm::utility::one<SolutionType>();
238 for (
auto const& entry : row) {
239 tmp.push_back(entry.getValue().lower());
240 probLeft -= entry.getValue().lower();
243 auto const& rowIter = row.begin();
244 for (uint64_t i = 0;
i < row.getNumberOfEntries();
i++, schedulerIterator++) {
246 auto const& entry = rowIter[*schedulerIterator];
247 auto const diameter = entry.getValue().upper() - entry.getValue().lower();
249 tmp[*schedulerIterator] += value;
256 for (uint64_t i = 0;
i < row.getNumberOfEntries();
i++) {
257 auto const& entry = rowIter[
i];
258 newMatrixBuilder.addNextValue(rowIndex, entry.getColumn(), tmp[i]);
261 STORM_LOG_ASSERT(schedulerIterator == scheduler.end(),
"Offset issue in scheduler?");
265 std::vector<SolutionType> b;
266 for (
auto const& entry : subB) {
267 if (dir == OptimizationDirection::Maximize) {
268 b.push_back(entry.upper());
270 b.push_back(entry.lower());
274 auto const submatrix = newMatrixBuilder.build();
277 if (!linearEquationSolver) {
279 linearEquationSolver = this->linearEquationSolverFactory->create(env, std::move(submatrix));
280 linearEquationSolver->setBoundsFromOtherSolver(*
this);
281 linearEquationSolver->setCachingEnabled(
true);
284 linearEquationSolver->setMatrix(std::move(submatrix));
287 return linearEquationSolver->solveEquations(env, x, b);
289 STORM_LOG_ASSERT(subB.size() == x.size(),
"Sizes of subB and x do not coincide.");
290 STORM_LOG_ASSERT(this->linearEquationSolverFactory !=
nullptr,
"Wrong constructor was called.");
297 if (convertToEquationSystem) {
300 storm::utility::vector::selectVectorValues<ValueType>(subB, scheduler, this->A->getRowGroupIndices(), originalB);
303 if (!linearEquationSolver) {
305 linearEquationSolver = this->linearEquationSolverFactory->create(env, std::move(submatrix));
306 linearEquationSolver->setBoundsFromOtherSolver(*
this);
307 linearEquationSolver->setCachingEnabled(
true);
310 linearEquationSolver->setMatrix(std::move(submatrix));
313 return linearEquationSolver->solveEquations(env, x, subB);
317template<
typename ValueType,
typename SolutionType>
318bool IterativeMinMaxLinearEquationSolver<ValueType, SolutionType>::solveEquationsPolicyIteration(Environment
const& env,
OptimizationDirection dir,
319 std::vector<SolutionType>& x,
320 std::vector<ValueType>
const& b)
const {
321 std::vector<storm::storage::sparse::state_type> scheduler =
322 this->hasInitialScheduler() ? this->getInitialScheduler() :
std::vector<
storm::storage::sparse::
state_type>(this->A->getRowGroupCount());
323 return performPolicyIteration(env, dir, x, b, std::move(scheduler));
326template<
typename ValueType,
typename SolutionType>
327bool IterativeMinMaxLinearEquationSolver<ValueType, SolutionType>::performPolicyIteration(
328 Environment
const& env,
OptimizationDirection dir, std::vector<SolutionType>& x, std::vector<ValueType>
const& b,
329 std::vector<storm::storage::sparse::state_type>&& initialPolicy)
const {
330 if constexpr (std::is_same_v<ValueType, storm::Interval>) {
331 STORM_LOG_THROW(
false, storm::exceptions::NotImplementedException,
"We did not implement policy iteration for interval-based models.");
334 std::vector<storm::storage::sparse::state_type> scheduler = std::move(initialPolicy);
336 if (!auxiliaryRowGroupVector) {
337 auxiliaryRowGroupVector = std::make_unique<std::vector<ValueType>>(this->A->getRowGroupCount());
339 std::vector<ValueType>& subB = *auxiliaryRowGroupVector;
342 std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver;
344 std::unique_ptr<storm::Environment> environmentOfSolverStorage;
345 auto precOfSolver = env.solver().getPrecisionOfLinearEquationSolver(env.solver().getLinearEquationSolverType());
347 bool changePrecision = precOfSolver.first && precOfSolver.first.get() > env.solver().minMax().getPrecision();
348 bool changeRelative = precOfSolver.second && !precOfSolver.second.get() && env.solver().minMax().getRelativeTerminationCriterion();
349 if (changePrecision || changeRelative) {
350 environmentOfSolverStorage = std::make_unique<storm::Environment>(env);
351 boost::optional<storm::RationalNumber> newPrecision;
352 boost::optional<bool> newRelative;
353 if (changePrecision) {
354 newPrecision = env.solver().minMax().getPrecision();
356 if (changeRelative) {
359 environmentOfSolverStorage->solver().setLinearEquationSolverPrecision(newPrecision, newRelative);
362 storm::Environment const& environmentOfSolver = environmentOfSolverStorage ? *environmentOfSolverStorage : env;
365 uint64_t iterations = 0;
366 this->startMeasureProgress();
369 solveInducedEquationSystem(environmentOfSolver, solver, scheduler, x, subB, b, dir);
372 bool schedulerImproved =
false;
374 for (uint_fast64_t group = 0; group < this->A->getRowGroupCount(); ++group) {
375 if (!this->choiceFixedForRowGroup || !this->choiceFixedForRowGroup.get()[group]) {
377 uint_fast64_t currentChoice = scheduler[group];
378 for (uint_fast64_t choice = this->A->getRowGroupIndices()[group]; choice < this->A->getRowGroupIndices()[group + 1]; ++choice) {
380 if (choice - this->A->getRowGroupIndices()[group] == currentChoice) {
385 ValueType choiceValue = storm::utility::zero<ValueType>();
386 for (
auto const& entry : this->A->getRow(choice)) {
387 choiceValue += entry.getValue() * x[entry.getColumn()];
389 choiceValue += b[choice];
394 if (valueImproved(dir, x[group], choiceValue)) {
395 schedulerImproved =
true;
396 scheduler[group] = choice - this->A->getRowGroupIndices()[group];
397 x[group] = std::move(choiceValue);
404 if (!schedulerImproved) {
412 iterations, env.solver().minMax().getMaximalNumberOfIterations());
415 this->showProgressIterative(iterations);
419 this->reportStatus(status, iterations);
422 if (this->isTrackSchedulerSet()) {
423 this->schedulerChoices = std::move(scheduler);
426 if (!this->isCachingEnabled()) {
434template<
typename ValueType,
typename SolutionType>
435bool IterativeMinMaxLinearEquationSolver<ValueType, SolutionType>::valueImproved(
OptimizationDirection dir, ValueType
const& value1,
436 ValueType
const& value2)
const {
437 if (dir == OptimizationDirection::Minimize) {
438 return value2 < value1;
440 return value2 > value1;
444template<
typename ValueType,
typename SolutionType>
446 Environment const& env, boost::optional<storm::solver::OptimizationDirection>
const& direction,
bool const& hasInitialScheduler)
const {
450 bool needsLinEqSolver =
false;
451 needsLinEqSolver |= method == MinMaxMethod::PolicyIteration;
452 needsLinEqSolver |= method == MinMaxMethod::ValueIteration && (this->hasInitialScheduler() || hasInitialScheduler);
453 needsLinEqSolver |= method == MinMaxMethod::ViToPi;
456 if constexpr (std::is_same_v<ValueType, storm::Interval>) {
457 STORM_LOG_ASSERT(!needsLinEqSolver,
"Intervals should not require a linear equation solver.");
459 }
else if (needsLinEqSolver) {
465 if (method == MinMaxMethod::ValueIteration) {
466 if (!this->hasUniqueSolution()) {
472 if (!direction || direction.get() == OptimizationDirection::Maximize) {
475 if (!direction || direction.get() == OptimizationDirection::Minimize) {
480 }
else if (method == MinMaxMethod::OptimisticValueIteration) {
482 if (!this->hasUniqueSolution()) {
487 }
else if (method == MinMaxMethod::GuessingValueIteration) {
489 if (!this->hasUniqueSolution()) {
493 }
else if (method == MinMaxMethod::IntervalIteration) {
495 if (!this->hasUniqueSolution()) {
499 }
else if (method == MinMaxMethod::RationalSearch) {
503 if (!this->hasUniqueSolution()) {
508 }
else if (method == MinMaxMethod::PolicyIteration) {
513 if (!this->hasNoEndComponents() && !this->hasInitialScheduler()) {
516 }
else if (method == MinMaxMethod::SoundValueIteration) {
517 if (!this->hasUniqueSolution()) {
521 }
else if (method == MinMaxMethod::ViToPi) {
523 if (!this->hasUniqueSolution()) {
527 STORM_LOG_THROW(
false, storm::exceptions::InvalidEnvironmentException,
"Unsupported technique for iterative MinMax linear equation solver.");
532template<
typename ValueType,
typename SolutionType>
534 ValueType result = storm::utility::zero<ValueType>();
535 auto oldValueIt = oldValues.begin();
536 for (
auto value : relevantValues) {
537 result = storm::utility::max<ValueType>(result, storm::utility::abs<ValueType>(allValues[value] - *oldValueIt));
543template<
typename ValueType,
typename SolutionType>
544ValueType
computeMaxAbsDiff(std::vector<ValueType>
const& allOldValues, std::vector<ValueType>
const& allNewValues,
546 ValueType result = storm::utility::zero<ValueType>();
547 for (
auto value : relevantValues) {
548 result = storm::utility::max<ValueType>(result, storm::utility::abs<ValueType>(allNewValues[value] - allOldValues[value]));
553template<
typename ValueType,
typename SolutionType>
555 std::vector<SolutionType>& x,
556 std::vector<ValueType>
const& b)
const {
557 if constexpr (std::is_same_v<ValueType, storm::Interval>) {
558 STORM_LOG_THROW(
false, storm::exceptions::NotImplementedException,
"We did not implement optimistic value iteration for interval-based models.");
563 x.assign(x.size(), storm::utility::zero<SolutionType>());
564 if (this->isTrackSchedulerSet()) {
565 this->schedulerChoices = std::vector<uint_fast64_t>(x.size(), 0);
572 helper::OptimisticValueIterationHelper<ValueType, false> oviHelper(viOperatorNontriv);
574 std::optional<ValueType> lowerBound, upperBound;
575 if (this->hasLowerBound()) {
576 lowerBound = this->getLowerBound(
true);
578 if (this->hasUpperBound()) {
579 upperBound = this->getUpperBound(
true);
581 uint64_t numIterations{0};
582 auto oviCallback = [&](
SolverStatus const& current, std::vector<ValueType>
const& v) {
583 this->showProgressIterative(numIterations);
586 this->createLowerBoundsVector(x);
587 std::optional<ValueType> guessingFactor;
591 this->startMeasureProgress();
593 upperBound, oviCallback);
594 this->reportStatus(status, numIterations);
597 if (this->isTrackSchedulerSet()) {
598 this->extractScheduler(x, b, dir, this->isUncertaintyRobust());
601 if (!this->isCachingEnabled()) {
609template<
typename ValueType,
typename SolutionType>
610bool IterativeMinMaxLinearEquationSolver<ValueType, SolutionType>::solveEquationsGuessingValueIteration(Environment
const& env,
OptimizationDirection dir,
611 std::vector<SolutionType>& x,
612 std::vector<ValueType>
const& b)
const {
613 if constexpr (std::is_same_v<ValueType, storm::Interval>) {
614 STORM_LOG_THROW(
false, storm::exceptions::NotImplementedException,
"We did not implement guessing value iteration for interval-based models.");
620 auto upperX = std::make_unique<std::vector<SolutionType>>(x.size());
624 uint64_t numIterations{0};
626 auto gviCallback = [&](helper::GVIData<ValueType>
const& data) {
627 this->showProgressIterative(numIterations);
628 bool terminateEarly = this->hasCustomTerminationCondition() && this->getTerminationCondition().terminateNow(data.x,
SolverGuarantee::LessOrEqual) &&
630 return this->updateStatus(data.status, terminateEarly, numIterations, env.solver().minMax().getMaximalNumberOfIterations());
633 this->createLowerBoundsVector(lowerX);
634 this->createUpperBoundsVector(*upperX);
636 this->startMeasureProgress();
637 auto statusIters = helper.solveEquations(lowerX, *upperX, b, numIterations,
638 storm::utility::convertNumber<ValueType>(env.solver().minMax().getPrecision()), dir, gviCallback);
639 auto two = storm::utility::convertNumber<ValueType>(2.0);
640 storm::utility::vector::applyPointwise<ValueType, ValueType, ValueType>(
641 lowerX, *upperX, x, [&two](ValueType
const& a, ValueType
const& b) -> ValueType {
return (a + b) / two; });
643 this->reportStatus(statusIters, numIterations);
646 if (this->isTrackSchedulerSet()) {
647 this->extractScheduler(x, b, dir, this->isUncertaintyRobust());
650 if (!this->isCachingEnabled()) {
658template<
typename ValueType,
typename SolutionType>
659bool IterativeMinMaxLinearEquationSolver<ValueType, SolutionType>::solveEquationsValueIteration(Environment
const& env,
OptimizationDirection dir,
660 std::vector<SolutionType>& x,
661 std::vector<ValueType>
const& b)
const {
666 if (this->hasInitialScheduler()) {
667 if (!auxiliaryRowGroupVector) {
668 auxiliaryRowGroupVector = std::make_unique<std::vector<ValueType>>(this->A->getRowGroupCount());
671 std::unique_ptr<storm::solver::LinearEquationSolver<SolutionType>> linEqSolver;
673 std::unique_ptr<storm::Environment> environmentOfSolverStorage;
674 auto precOfSolver = env.solver().getPrecisionOfLinearEquationSolver(env.solver().getLinearEquationSolverType());
676 bool changePrecision = precOfSolver.first && precOfSolver.first.get() > env.solver().minMax().getPrecision();
677 bool changeRelative = precOfSolver.second && !precOfSolver.second.get() && env.solver().minMax().getRelativeTerminationCriterion();
678 if (changePrecision || changeRelative) {
679 environmentOfSolverStorage = std::make_unique<storm::Environment>(env);
680 boost::optional<storm::RationalNumber> newPrecision;
681 boost::optional<bool> newRelative;
682 if (changePrecision) {
683 newPrecision = env.solver().minMax().getPrecision();
685 if (changeRelative) {
688 environmentOfSolverStorage->solver().setLinearEquationSolverPrecision(newPrecision, newRelative);
691 storm::Environment const& environmentOfSolver = environmentOfSolverStorage ? *environmentOfSolverStorage : env;
693 solveInducedEquationSystem(environmentOfSolver, linEqSolver, this->getInitialScheduler(), x, *auxiliaryRowGroupVector, b, dir);
697 }
else if (!this->hasUniqueSolution()) {
699 this->createLowerBoundsVector(x);
702 this->createUpperBoundsVector(x);
705 }
else if (this->hasCustomTerminationCondition()) {
707 this->createLowerBoundsVector(x);
710 this->createUpperBoundsVector(x);
715 uint64_t numIterations{0};
717 this->showProgressIterative(numIterations);
718 return this->updateStatus(current, x, guarantee, numIterations, env.solver().minMax().getMaximalNumberOfIterations());
720 this->startMeasureProgress();
722 if (this->A->hasTrivialRowGrouping()) {
725 auto status = viHelper.VI(x, b, numIterations, env.solver().minMax().getRelativeTerminationCriterion(),
726 storm::utility::convertNumber<SolutionType>(env.solver().minMax().getPrecision()), dir, viCallback,
727 env.solver().minMax().getMultiplicationStyle(), this->isUncertaintyRobust());
728 this->reportStatus(status, numIterations);
731 if (this->isTrackSchedulerSet()) {
732 this->extractScheduler(x, b, dir, this->isUncertaintyRobust());
735 if (!this->isCachingEnabled()) {
743 auto status = viHelper.VI(x, b, numIterations, env.solver().minMax().getRelativeTerminationCriterion(),
744 storm::utility::convertNumber<SolutionType>(env.solver().minMax().getPrecision()), dir, viCallback,
745 env.solver().minMax().getMultiplicationStyle(), this->isUncertaintyRobust());
746 this->reportStatus(status, numIterations);
749 if (this->isTrackSchedulerSet()) {
750 this->extractScheduler(x, b, dir, this->isUncertaintyRobust());
753 if (!this->isCachingEnabled()) {
761template<
typename ValueType,
typename SolutionType>
772template<
typename ValueType,
typename SolutionType>
774 std::vector<SolutionType>& x,
775 std::vector<ValueType>
const& b)
const {
776 if constexpr (std::is_same_v<ValueType, storm::Interval>) {
777 STORM_LOG_THROW(
false, storm::exceptions::NotImplementedException,
"We did not implement intervaliteration for interval-based models");
781 helper::IntervalIterationHelper<ValueType, false> iiHelper(viOperatorNontriv);
783 auto lowerBoundsCallback = [&](std::vector<SolutionType>& vector) { this->createLowerBoundsVector(vector); };
784 auto upperBoundsCallback = [&](std::vector<SolutionType>& vector) { this->createUpperBoundsVector(vector); };
786 uint64_t numIterations{0};
787 auto iiCallback = [&](helper::IIData<ValueType>
const& data) {
788 this->showProgressIterative(numIterations);
789 bool terminateEarly = this->hasCustomTerminationCondition() && this->getTerminationCondition().terminateNow(data.x,
SolverGuarantee::LessOrEqual) &&
793 std::optional<storm::storage::BitVector> optionalRelevantValues;
794 if (this->hasRelevantValues()) {
795 optionalRelevantValues = this->getRelevantValues();
797 this->startMeasureProgress();
799 dir, iiCallback, optionalRelevantValues);
800 this->reportStatus(status, numIterations);
803 if (this->isTrackSchedulerSet()) {
804 this->extractScheduler(x, b, dir, this->isUncertaintyRobust());
807 if (!this->isCachingEnabled()) {
815template<
typename ValueType,
typename SolutionType>
816bool IterativeMinMaxLinearEquationSolver<ValueType, SolutionType>::solveEquationsSoundValueIteration(Environment
const& env,
OptimizationDirection dir,
817 std::vector<SolutionType>& x,
818 std::vector<ValueType>
const& b)
const {
819 if constexpr (std::is_same_v<ValueType, storm::Interval>) {
820 STORM_LOG_THROW(
false, storm::exceptions::NotImplementedException,
"SoundVI does not handle interval-based models");
824 assert(x.size() == this->A->getRowGroupCount());
826 std::optional<ValueType> lowerBound, upperBound;
827 if (this->hasLowerBound()) {
828 lowerBound = this->getLowerBound(
true);
830 if (this->hasUpperBound()) {
831 upperBound = this->getUpperBound(
true);
836 auto precision = storm::utility::convertNumber<ValueType>(env.solver().minMax().getPrecision());
837 uint64_t numIterations{0};
838 auto sviCallback = [&](
typename helper::SoundValueIterationHelper<ValueType, false>::SVIData
const& current) {
839 this->showProgressIterative(numIterations);
840 return this->updateStatus(current.status,
841 this->hasCustomTerminationCondition() && current.checkCustomTerminationCondition(this->getTerminationCondition()),
842 numIterations, env.solver().minMax().getMaximalNumberOfIterations());
844 this->startMeasureProgress();
845 helper::SoundValueIterationHelper<ValueType, false> sviHelper(viOperatorNontriv);
846 std::optional<storm::storage::BitVector> optionalRelevantValues;
847 if (this->hasRelevantValues()) {
848 optionalRelevantValues = this->getRelevantValues();
850 auto status = sviHelper.SVI(x, b, numIterations, env.solver().minMax().getRelativeTerminationCriterion(), precision, dir, lowerBound, upperBound,
851 sviCallback, optionalRelevantValues);
854 if (this->isTrackSchedulerSet()) {
855 this->extractScheduler(x, b, dir, this->isUncertaintyRobust());
858 this->reportStatus(status, numIterations);
860 if (!this->isCachingEnabled()) {
868template<
typename ValueType,
typename SolutionType>
869bool IterativeMinMaxLinearEquationSolver<ValueType, SolutionType>::solveEquationsViToPi(Environment
const& env,
OptimizationDirection dir,
870 std::vector<SolutionType>& x, std::vector<ValueType>
const& b)
const {
871 if constexpr (std::is_same_v<ValueType, storm::Interval>) {
872 STORM_LOG_THROW(
false, storm::exceptions::NotImplementedException,
"ViToPi does not handle interval-based models");
876 std::vector<storm::storage::sparse::state_type> initialSched;
878 Environment viEnv = env;
879 viEnv.solver().minMax().setMethod(MinMaxMethod::ValueIteration);
880 viEnv.solver().setForceExact(
false);
881 viEnv.solver().setForceSoundness(
false);
882 auto impreciseSolver = GeneralMinMaxLinearEquationSolverFactory<double>().create(viEnv, this->A->template toValueType<double>());
883 impreciseSolver->setHasUniqueSolution(this->hasUniqueSolution());
884 impreciseSolver->setTrackScheduler(
true);
885 if (this->hasInitialScheduler()) {
886 auto initSched = this->getInitialScheduler();
887 impreciseSolver->setInitialScheduler(std::move(initSched));
889 auto impreciseSolverReq = impreciseSolver->getRequirements(viEnv, dir);
890 STORM_LOG_THROW(!impreciseSolverReq.hasEnabledCriticalRequirement(), storm::exceptions::UnmetRequirementException,
891 "The value-iteration based solver has an unmet requirement: " << impreciseSolverReq.getEnabledRequirementsAsString());
892 impreciseSolver->setRequirementsChecked(
true);
893 auto xVi = storm::utility::vector::convertNumericVector<double>(x);
894 auto bVi = storm::utility::vector::convertNumericVector<double>(b);
895 impreciseSolver->solveEquations(viEnv, dir, xVi, bVi);
896 initialSched = impreciseSolver->getSchedulerChoices();
898 STORM_LOG_INFO(
"Found initial policy using Value Iteration. Starting Policy iteration now.");
899 return performPolicyIteration(env, dir, x, b, std::move(initialSched));
902template<
typename ValueType,
typename SolutionType>
903bool IterativeMinMaxLinearEquationSolver<ValueType, SolutionType>::solveEquationsRationalSearch(Environment
const& env,
OptimizationDirection dir,
904 std::vector<SolutionType>& x,
905 std::vector<ValueType>
const& b)
const {
906 if constexpr (std::is_same_v<ValueType, storm::Interval>) {
907 STORM_LOG_THROW(
false, storm::exceptions::NotImplementedException,
"Rational search does not handle interval-based models");
912 std::shared_ptr<helper::ValueIterationOperator<storm::RationalNumber, false>> exactOp;
913 std::shared_ptr<helper::ValueIterationOperator<double, false>> impreciseOp;
914 std::function<bool(uint64_t, uint64_t)> fixedChoicesCallback;
915 if (this->choiceFixedForRowGroup) {
917 assert(this->initialScheduler);
918 fixedChoicesCallback = [&](uint64_t groupIndex, uint64_t localRowIndex) {
919 return this->choiceFixedForRowGroup->get(groupIndex) && this->initialScheduler->at(groupIndex) != localRowIndex;
923 if constexpr (std::is_same_v<ValueType, storm::RationalNumber>) {
924 exactOp = viOperatorNontriv;
925 impreciseOp = std::make_shared<helper::ValueIterationOperator<double, false>>();
926 impreciseOp->setMatrixBackwards(this->A->template toValueType<double>(), &this->A->getRowGroupIndices());
927 if (this->choiceFixedForRowGroup) {
928 impreciseOp->setIgnoredRows(
true, fixedChoicesCallback);
930 }
else if constexpr (std::is_same_v<ValueType, double>) {
931 impreciseOp = viOperatorNontriv;
932 exactOp = std::make_shared<helper::ValueIterationOperator<storm::RationalNumber, false>>();
933 exactOp->setMatrixBackwards(this->A->template toValueType<storm::RationalNumber>(), &this->A->getRowGroupIndices());
934 if (this->choiceFixedForRowGroup) {
935 exactOp->setIgnoredRows(
true, fixedChoicesCallback);
940 uint64_t numIterations{0};
942 this->showProgressIterative(numIterations);
943 return this->updateStatus(current, x,
SolverGuarantee::None, numIterations, env.solver().minMax().getMaximalNumberOfIterations());
945 this->startMeasureProgress();
946 auto status = rsHelper.RS(x, b, numIterations, storm::utility::convertNumber<ValueType>(env.solver().minMax().getPrecision()), dir, rsCallback);
948 this->reportStatus(status, numIterations);
951 if (this->isTrackSchedulerSet()) {
952 this->extractScheduler(x, b, dir, this->isUncertaintyRobust());
955 if (!this->isCachingEnabled()) {
963template<
typename ValueType,
typename SolutionType>
965 auxiliaryRowGroupVector.reset();
966 if (viOperatorTriv) {
967 viOperatorTriv.reset();
969 if (viOperatorNontriv) {
970 viOperatorNontriv.reset();
SolverEnvironment & solver()
uint64_t const & getMaximalNumberOfIterations() const
storm::RationalNumber const & getPrecision() const
bool const & isMethodSetFromDefault() const
storm::solver::MinMaxMethod const & getMethod() const
bool isForceRequireUnique() const
bool const & getRelativeTerminationCriterion() const
std::optional< storm::RationalNumber > const & getUpperBoundGuessingFactor() const
OviSolverEnvironment const & ovi() const
MinMaxSolverEnvironment & minMax()
bool isForceExact() const
bool isForceSoundness() const
virtual void clearCache() const override
Clears the currently cached data that has been stored during previous calls of the solver.
IterativeMinMaxLinearEquationSolver(std::unique_ptr< LinearEquationSolverFactory< SolutionType > > &&linearEquationSolverFactory)
virtual bool internalSolveEquations(Environment const &env, OptimizationDirection dir, std::vector< SolutionType > &x, std::vector< ValueType > const &b) const override
virtual MinMaxLinearEquationSolverRequirements getRequirements(Environment const &env, boost::optional< storm::solver::OptimizationDirection > const &direction=boost::none, bool const &hasInitialScheduler=false) const override
Retrieves the requirements of this solver for solving equations with the current settings.
virtual void clearCache() const
Clears the currently cached data that has been stored during previous calls of the solver.
MinMaxLinearEquationSolverRequirements & requireBounds(bool critical=true)
MinMaxLinearEquationSolverRequirements & requireUniqueSolution(bool critical=true)
MinMaxLinearEquationSolverRequirements & requireLowerBounds(bool critical=true)
MinMaxLinearEquationSolverRequirements & requireValidInitialScheduler(bool critical=true)
MinMaxLinearEquationSolverRequirements & requireUpperBounds(bool critical=true)
Implements guessing value iteration.
Implements rational search.
Helper class to extract optimal scheduler choices from a MinMax equation system solution.
A bit vector that is internally represented as a vector of 64-bit values.
A class that can be used to build a sparse matrix by adding value by value.
A class that holds a possibly non-square matrix in the compressed row storage format.
void convertToEquationSystem()
Transforms the matrix into an equation system.
SparseMatrix selectRowsFromRowGroups(std::vector< index_type > const &rowGroupToRowIndexMapping, bool insertDiagonalEntries=true) const
Selects exactly one row from each row group of this matrix and returns the resulting matrix.
#define STORM_LOG_INFO(message)
#define STORM_LOG_WARN(message)
#define STORM_LOG_ERROR(message)
#define STORM_LOG_ASSERT(cond, message)
#define STORM_LOG_WARN_COND(cond, message)
#define STORM_LOG_THROW(cond, exception, message)
SFTBDDChecker::ValueType ValueType
bool constexpr maximize(OptimizationDirection d)
void preserveOldRelevantValues(std::vector< ValueType > const &allValues, storm::storage::BitVector const &relevantValues, std::vector< ValueType > &oldValues)
std::string toString(GurobiSolverMethod const &method)
Yields a string representation of the GurobiSolverMethod.
ValueType computeMaxAbsDiff(std::vector< ValueType > const &allValues, storm::storage::BitVector const &relevantValues, std::vector< ValueType > const &oldValues)
void selectVectorValues(std::vector< T > &vector, storm::storage::BitVector const &positions, std::vector< T > const &values)
Selects the elements from a vector at the specified positions and writes them consecutively into anot...
bool hasNonZeroEntry(std::vector< T > const &v)
ValueType min(ValueType const &first, ValueType const &second)
bool isZero(ValueType const &a)