diff --git a/src/storm/modelchecker/multiobjective/pcaa/StandardPcaaWeightVectorChecker.cpp b/src/storm/modelchecker/multiobjective/pcaa/StandardPcaaWeightVectorChecker.cpp index beaae1814..bd93a44ef 100644 --- a/src/storm/modelchecker/multiobjective/pcaa/StandardPcaaWeightVectorChecker.cpp +++ b/src/storm/modelchecker/multiobjective/pcaa/StandardPcaaWeightVectorChecker.cpp @@ -8,6 +8,8 @@ #include "storm/models/sparse/MarkovAutomaton.h" #include "storm/models/sparse/StandardRewardModel.h" #include "storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.h" +#include "storm/modelchecker/prctl/helper/DsMpiUpperRewardBoundsComputer.h" +#include "storm/modelchecker/prctl/helper/BaierUpperRewardBoundsComputer.h" #include "storm/solver/MinMaxLinearEquationSolver.h" #include "storm/utility/graph.h" #include "storm/utility/macros.h" @@ -176,20 +178,17 @@ namespace storm { std::unique_ptr> solver = solverFactory.create(env, ecQuotient->matrix); solver->setTrackScheduler(true); solver->setHasUniqueSolution(true); + solver->setOptimizationDirection(storm::solver::OptimizationDirection::Maximize); auto req = solver->getRequirements(env, storm::solver::OptimizationDirection::Maximize, true); - boost::optional lowerBound = this->computeWeightedResultBound(true, weightVector, objectivesWithNoUpperTimeBound); - if (lowerBound) { - solver->setLowerBound(lowerBound.get()); + setBoundsToSolver(*solver, req.requiresLowerBounds(), req.requiresUpperBounds(), weightVector, objectivesWithNoUpperTimeBound, ecQuotient->matrix, ecQuotient->rowsWithSumLessOne, ecQuotient->auxChoiceValues); + if (solver->hasLowerBound()) { req.clearLowerBounds(); } - boost::optional upperBound = this->computeWeightedResultBound(false, weightVector, objectivesWithNoUpperTimeBound); - if (upperBound) { - solver->setUpperBound(upperBound.get()); + if (solver->hasUpperBound()) { req.clearUpperBounds(); } STORM_LOG_THROW(req.empty(), storm::exceptions::UncheckedRequirementException, "At least one requirement was not checked."); solver->setRequirementsChecked(true); - solver->setOptimizationDirection(storm::solver::OptimizationDirection::Maximize); // Use the (0...0) vector as initial guess for the solution. std::fill(ecQuotient->auxStateValues.begin(), ecQuotient->auxStateValues.end(), storm::utility::zero()); @@ -209,11 +208,11 @@ namespace storm { if (storm::solver::minimize(this->objectives[objIndex].formula->getOptimalityType())) { storm::utility::vector::scaleVectorInPlace(objectiveResults[objIndex], -storm::utility::one()); } - for (uint_fast64_t objIndex2 = 0; objIndex2 < this->objectives.size(); ++objIndex2) { - if (objIndex != objIndex2) { - objectiveResults[objIndex2] = std::vector(transitionMatrix.getRowGroupCount(), storm::utility::zero()); - } - } + for (uint_fast64_t objIndex2 = 0; objIndex2 < this->objectives.size(); ++objIndex2) { + if (objIndex != objIndex2) { + objectiveResults[objIndex2] = std::vector(transitionMatrix.getRowGroupCount(), storm::utility::zero()); + } + } } else { storm::storage::SparseMatrix deterministicMatrix = transitionMatrix.selectRowsFromRowGroups(this->optimalChoices, true); storm::storage::SparseMatrix deterministicBackwardTransitions = deterministicMatrix.transpose(); @@ -263,17 +262,19 @@ namespace storm { std::vector b = storm::utility::vector::filterVector(deterministicStateRewards, maybeStates); // Now solve the resulting equation system. - std::unique_ptr> solver = linearEquationSolverFactory.create(env, std::move(submatrix)); + std::unique_ptr> solver = linearEquationSolverFactory.create(env, submatrix); auto req = solver->getRequirements(env); solver->clearBounds(); - if (obj.lowerResultBound) { + storm::storage::BitVector submatrixRowsWithSumLessOne = deterministicMatrix.getRowFilter(maybeStates, maybeStates) % maybeStates; + submatrixRowsWithSumLessOne.complement(); + this->setBoundsToSolver(*solver, req.requiresLowerBounds(), req.requiresUpperBounds(), objIndex, submatrix, submatrixRowsWithSumLessOne, b); + if (solver->hasLowerBound()) { req.clearLowerBounds(); - solver->setLowerBound(*obj.lowerResultBound); } - if (obj.upperResultBound) { - solver->setUpperBound(*obj.upperResultBound); + if (solver->hasUpperBound()) { req.clearUpperBounds(); } + STORM_LOG_THROW(req.empty(), storm::exceptions::UncheckedRequirementException, "At least one requirement of the LinearEquationSolver was not met."); solver->solveEquations(env, x, b); @@ -312,11 +313,26 @@ namespace storm { // Remove neutral end components, i.e., ECs in which no reward is earned. auto ecElimResult = storm::transformer::EndComponentEliminator::transform(transitionMatrix, subsystemStates, ecChoicesHint & newReward0Choices, reward0EStates); + storm::storage::BitVector rowsWithSumLessOne(ecElimResult.matrix.getRowCount(), false); + for (uint64_t row = 0; row < rowsWithSumLessOne.size(); ++row) { + if (ecElimResult.matrix.getRow(row).getNumberOfEntries() == 0) { + rowsWithSumLessOne.set(row, true); + } else { + for (auto const& entry : transitionMatrix.getRow(ecElimResult.newToOldRowMapping[row])) { + if (!subsystemStates.get(entry.getColumn())) { + rowsWithSumLessOne.set(row, true); + break; + } + } + } + } + ecQuotient = EcQuotient(); ecQuotient->matrix = std::move(ecElimResult.matrix); ecQuotient->ecqToOriginalChoiceMapping = std::move(ecElimResult.newToOldRowMapping); ecQuotient->originalToEcqStateMapping = std::move(ecElimResult.oldToNewStateMapping); ecQuotient->origReward0Choices = std::move(newReward0Choices); + ecQuotient->rowsWithSumLessOne = std::move(rowsWithSumLessOne); ecQuotient->auxStateValues.reserve(transitionMatrix.getRowGroupCount()); ecQuotient->auxStateValues.resize(ecQuotient->matrix.getRowGroupCount()); ecQuotient->auxChoiceValues.reserve(transitionMatrix.getRowCount()); @@ -324,6 +340,91 @@ namespace storm { } } + template + void StandardPcaaWeightVectorChecker::setBoundsToSolver(storm::solver::AbstractEquationSolver& solver, bool requiresLower, bool requiresUpper, uint64_t objIndex, storm::storage::SparseMatrix const& transitions, storm::storage::BitVector const& rowsWithSumLessOne, std::vector const& rewards) const { + + // Check whether bounds are already available + if (this->objectives[objIndex].lowerResultBound) { + solver.setLowerBound(this->objectives[objIndex].lowerResultBound.get()); + } + if (this->objectives[objIndex].upperResultBound) { + solver.setUpperBound(this->objectives[objIndex].upperResultBound.get()); + } + + if ((requiresLower && !solver.hasLowerBound()) || (requiresUpper && !solver.hasUpperBound())) { + computeAndSetBoundsToSolver(solver, requiresLower, requiresUpper, transitions, rowsWithSumLessOne, rewards); + } + } + + template + void StandardPcaaWeightVectorChecker::setBoundsToSolver(storm::solver::AbstractEquationSolver& solver, bool requiresLower, bool requiresUpper, std::vector const& weightVector, storm::storage::BitVector const& objectiveFilter, storm::storage::SparseMatrix const& transitions, storm::storage::BitVector const& rowsWithSumLessOne, std::vector const& rewards) const { + + // Check whether bounds are already available + boost::optional lowerBound = this->computeWeightedResultBound(true, weightVector, objectiveFilter); + if (lowerBound) { + solver.setLowerBound(lowerBound.get()); + } + boost::optional upperBound = this->computeWeightedResultBound(false, weightVector, objectiveFilter); + if (upperBound) { + solver.setUpperBound(upperBound.get()); + } + + if ((requiresLower && !solver.hasLowerBound()) || (requiresUpper && !solver.hasUpperBound())) { + computeAndSetBoundsToSolver(solver, requiresLower, requiresUpper, transitions, rowsWithSumLessOne, rewards); + } + } + + template + void StandardPcaaWeightVectorChecker::computeAndSetBoundsToSolver(storm::solver::AbstractEquationSolver& solver, bool requiresLower, bool requiresUpper, storm::storage::SparseMatrix const& transitions, storm::storage::BitVector const& rowsWithSumLessOne, std::vector const& rewards) const { + + // Compute the one step target probs + std::vector oneStepTargetProbs(transitions.getRowCount(), storm::utility::zero()); + for (auto const& row : rowsWithSumLessOne) { + oneStepTargetProbs[row] = storm::utility::one() - transitions.getRowSum(row); + } + + if (requiresLower && !solver.hasLowerBound()) { + // Compute lower bounds + std::vector negativeRewards; + negativeRewards.reserve(transitions.getRowCount()); + uint64_t row = 0; + for (auto const& rew : rewards) { + if (rew < storm::utility::zero()) { + negativeRewards.resize(row, storm::utility::zero()); + negativeRewards.push_back(-rew); + } + ++row; + } + if (!negativeRewards.empty()) { + negativeRewards.resize(row, storm::utility::zero()); + std::vector lowerBounds = storm::modelchecker::helper::DsMpiMdpUpperRewardBoundsComputer(transitions, negativeRewards, oneStepTargetProbs).computeUpperBounds(); + storm::utility::vector::scaleVectorInPlace(lowerBounds, -storm::utility::one()); + solver.setLowerBounds(std::move(lowerBounds)); + } else { + solver.setLowerBound(storm::utility::zero()); + } + } + + // Compute upper bounds + if (requiresUpper && !solver.hasUpperBound()) { + std::vector positiveRewards; + positiveRewards.reserve(transitions.getRowCount()); + uint64_t row = 0; + for (auto const& rew : rewards) { + if (rew > storm::utility::zero()) { + positiveRewards.resize(row, storm::utility::zero()); + positiveRewards.push_back(rew); + } + ++row; + } + if (!positiveRewards.empty()) { + positiveRewards.resize(row, storm::utility::zero()); + solver.setUpperBound(storm::modelchecker::helper::BaierUpperRewardBoundsComputer(transitions, positiveRewards, oneStepTargetProbs).computeUpperBound()); + } else { + solver.setUpperBound(storm::utility::zero()); + } + } + } template void StandardPcaaWeightVectorChecker::transformReducedSolutionToOriginalModel(storm::storage::SparseMatrix const& reducedMatrix, diff --git a/src/storm/modelchecker/multiobjective/pcaa/StandardPcaaWeightVectorChecker.h b/src/storm/modelchecker/multiobjective/pcaa/StandardPcaaWeightVectorChecker.h index 86cd79c7f..2f0a27f33 100644 --- a/src/storm/modelchecker/multiobjective/pcaa/StandardPcaaWeightVectorChecker.h +++ b/src/storm/modelchecker/multiobjective/pcaa/StandardPcaaWeightVectorChecker.h @@ -1,5 +1,6 @@ #pragma once +#include "storm/solver/AbstractEquationSolver.h" #include "storm/storage/BitVector.h" #include "storm/storage/SparseMatrix.h" #include "storm/storage/Scheduler.h" @@ -93,6 +94,11 @@ namespace storm { void updateEcQuotient(std::vector const& weightedRewardVector); + + void setBoundsToSolver(storm::solver::AbstractEquationSolver& solver, bool requiresLower, bool requiresUpper, uint64_t objIndex, storm::storage::SparseMatrix const& transitions, storm::storage::BitVector const& rowsWithSumLessOne, std::vector const& rewards) const; + void setBoundsToSolver(storm::solver::AbstractEquationSolver& solver, bool requiresLower, bool requiresUpper, std::vector const& weightVector, storm::storage::BitVector const& objectiveFilter, storm::storage::SparseMatrix const& transitions, storm::storage::BitVector const& rowsWithSumLessOne, std::vector const& rewards) const; + void computeAndSetBoundsToSolver(storm::solver::AbstractEquationSolver& solver, bool requiresLower, bool requiresUpper, storm::storage::SparseMatrix const& transitions, storm::storage::BitVector const& rowsWithSumLessOne, std::vector const& rewards) const; + /*! * Transforms the results of a min-max-solver that considers a reduced model (without end components) to a result for the original (unreduced) model */ @@ -141,6 +147,7 @@ namespace storm { std::vector ecqToOriginalChoiceMapping; std::vector originalToEcqStateMapping; storm::storage::BitVector origReward0Choices; + storm::storage::BitVector rowsWithSumLessOne; std::vector auxStateValues; std::vector auxChoiceValues;