diff --git a/src/modelchecker/multiobjective/Pcaa.cpp b/src/modelchecker/multiobjective/Pcaa.cpp new file mode 100644 index 000000000..1937835d4 --- /dev/null +++ b/src/modelchecker/multiobjective/Pcaa.cpp @@ -0,0 +1,60 @@ +#include "src/modelchecker/multiobjective/pcaa.h" + +#include "src/utility/macros.h" + +#include "src/models/sparse/Mdp.h" +#include "src/models/sparse/MarkovAutomaton.h" +#include "src/models/sparse/StandardRewardModel.h" +#include "src/modelchecker/multiobjective/pcaa/SparsePcaaPreprocessor.h" +#include "src/modelchecker/multiobjective/pcaa/SparsePcaaAchievabilityQuery.h" +#include "src/modelchecker/multiobjective/pcaa/SparsePcaaQuantitativeQuery.h" +#include "src/modelchecker/multiobjective/pcaa/SparsePcaaParetoQuery.h" + +#include "src/exceptions/InvalidArgumentException.h" + + +namespace storm { + namespace modelchecker { + namespace multiobjective { + + template + std::unique_ptr performPcaa(SparseModelType const& model, storm::logic::MultiObjectiveFormula const& formula) { + STORM_LOG_ASSERT(model.getInitialStates().getNumberOfSetBits() == 1, "Multi-objective Model checking on model with multiple initial states is not supported."); + +#ifdef STORM_HAVE_CARL + auto preprocessorResult = SparsePcaaPreprocessor::preprocess(model, formula); + STORM_LOG_DEBUG("Preprocessing done. Result: " << preprocessorResult); + + std::unique_ptr> query; + switch (preprocessorResult.queryType) { + case SparsePcaaPreprocessorReturnType::QueryType::Achievability: + query = std::unique_ptr> (new SparsePcaaAchievabilityQuery(preprocessorResult)); + break; + case SparsePcaaPreprocessorReturnType::QueryType::Quantitative: + query = std::unique_ptr> (new SparsePcaaQuantitativeQuery(preprocessorResult)); + break; + case SparsePcaaPreprocessorReturnType::QueryType::Pareto: + query = std::unique_ptr> (new SparsePcaaParetoQuery(preprocessorResult)); + break; + default: + STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "Unsupported multi-objective Query Type."); + break; + } + + return query->check(); +#else + STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "Multi-objective model checking requires carl."); + return nullptr; +#endif + } + + template std::unique_ptr performPcaa>(storm::models::sparse::Mdp const& model, storm::logic::MultiObjectiveFormula const& formula); + template std::unique_ptr performPcaa>(storm::models::sparse::MarkovAutomaton const& model, storm::logic::MultiObjectiveFormula const& formula); +#ifdef STORM_HAVE_CARL + template std::unique_ptr performPcaa>(storm::models::sparse::Mdp const& model, storm::logic::MultiObjectiveFormula const& formula); + // template std::unique_ptr performPcaa>(storm::models::sparse::MarkovAutomaton const& model, storm::logic::MultiObjectiveFormula const& formula); +#endif + + } + } +} diff --git a/src/modelchecker/multiobjective/Pcaa.h b/src/modelchecker/multiobjective/Pcaa.h new file mode 100644 index 000000000..c1e4a4566 --- /dev/null +++ b/src/modelchecker/multiobjective/Pcaa.h @@ -0,0 +1,20 @@ +#ifndef STORM_MODELCHECKER_MULTIOBJECTIVE_PCAA_H_ +#define STORM_MODELCHECKER_MULTIOBJECTIVE_PCAA_H_ + +#include + +#include "src/modelchecker/results/CheckResult.h" +#include "src/logic/Formulas.h" + +namespace storm { + namespace modelchecker { + namespace multiobjective { + + template + std::unique_ptr performPcaa(SparseModelType const& model, storm::logic::MultiObjectiveFormula const& formula); + + } + } +} + +#endif /* STORM_MODELCHECKER_MULTIOBJECTIVE_PCAA_H_ */ diff --git a/src/modelchecker/multiobjective/SparseMaMultiObjectiveModelChecker.cpp b/src/modelchecker/multiobjective/SparseMaMultiObjectiveModelChecker.cpp index e00cb79b5..f16b208b5 100644 --- a/src/modelchecker/multiobjective/SparseMaMultiObjectiveModelChecker.cpp +++ b/src/modelchecker/multiobjective/SparseMaMultiObjectiveModelChecker.cpp @@ -9,12 +9,7 @@ #include "src/modelchecker/results/ExplicitQualitativeCheckResult.h" #include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h" -#include "src/modelchecker/multiobjective/helper/SparseMultiObjectivePreprocessor.h" -#include "src/modelchecker/multiobjective/helper/SparseMaMultiObjectiveWeightVectorChecker.h" -#include "src/modelchecker/multiobjective/helper/SparseMultiObjectiveHelper.h" -#include "src/modelchecker/multiobjective/helper/SparseMultiObjectivePostprocessor.h" - -#include "src/utility/Stopwatch.h" +#include "src/modelchecker/multiobjective/pcaa.h" #include "src/exceptions/InvalidArgumentException.h" @@ -39,27 +34,8 @@ namespace storm { std::unique_ptr SparseMaMultiObjectiveModelChecker::checkMultiObjectiveFormula(CheckTask const& checkTask) { STORM_LOG_ASSERT(this->getModel().getInitialStates().getNumberOfSetBits() == 1, "Multi-objective Model checking on model with multiple initial states is not supported."); STORM_LOG_THROW(this->getModel().isClosed(), storm::exceptions::InvalidArgumentException, "Unable to check multi-objective formula in non-closed Markov automaton."); - std::unique_ptr result; - -#ifdef STORM_HAVE_CARL - storm::utility::Stopwatch swPreprocessing; - auto preprocessorData = helper::SparseMultiObjectivePreprocessor::preprocess(checkTask.getFormula(), this->getModel()); - swPreprocessing.pause(); - STORM_LOG_DEBUG("Preprocessing done. Data: " << preprocessorData); - - storm::utility::Stopwatch swHelper; - std::shared_ptr> weightVectorChecker( new helper::SparseMaMultiObjectiveWeightVectorChecker(preprocessorData)); - auto resultData = helper::SparseMultiObjectiveHelper::check(preprocessorData, weightVectorChecker); - swHelper.pause(); - STORM_LOG_DEBUG("Modelchecking done."); - - storm::utility::Stopwatch swPostprocessing; - result = helper::SparseMultiObjectivePostprocessor::postprocess(preprocessorData, resultData, swPreprocessing, swHelper, swPostprocessing); - STORM_LOG_DEBUG("Postprocessing done."); -#else - STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "Multi-objective model checking requires carl."); -#endif - return result; + + return multiobjective::performPcaa(this->getModel(), checkTask.getFormula()); } diff --git a/src/modelchecker/multiobjective/SparseMdpMultiObjectiveModelChecker.cpp b/src/modelchecker/multiobjective/SparseMdpMultiObjectiveModelChecker.cpp index d61596f1a..9f561b06e 100644 --- a/src/modelchecker/multiobjective/SparseMdpMultiObjectiveModelChecker.cpp +++ b/src/modelchecker/multiobjective/SparseMdpMultiObjectiveModelChecker.cpp @@ -9,14 +9,7 @@ #include "src/modelchecker/results/ExplicitQualitativeCheckResult.h" #include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h" -#include "src/modelchecker/multiobjective/helper/SparseMultiObjectivePreprocessor.h" -#include "src/modelchecker/multiobjective/helper/SparseMdpMultiObjectiveWeightVectorChecker.h" -#include "src/modelchecker/multiobjective/helper/SparseMultiObjectiveHelper.h" -#include "src/modelchecker/multiobjective/helper/SparseMultiObjectivePostprocessor.h" - -#include "src/utility/Stopwatch.h" - -#include "src/exceptions/NotImplementedException.h" +#include "src/modelchecker/multiobjective/pcaa.h" namespace storm { namespace modelchecker { @@ -38,27 +31,7 @@ namespace storm { template std::unique_ptr SparseMdpMultiObjectiveModelChecker::checkMultiObjectiveFormula(CheckTask const& checkTask) { STORM_LOG_ASSERT(this->getModel().getInitialStates().getNumberOfSetBits() == 1, "Multi Objective Model checking on model with multiple initial states is not supported."); - std::unique_ptr result; - -#ifdef STORM_HAVE_CARL - storm::utility::Stopwatch swPreprocessing; - auto preprocessorData = helper::SparseMultiObjectivePreprocessor::preprocess(checkTask.getFormula(), this->getModel()); - swPreprocessing.pause(); - STORM_LOG_DEBUG("Preprocessing done. Data: " << preprocessorData); - - storm::utility::Stopwatch swHelper; - std::shared_ptr> weightVectorChecker( new helper::SparseMdpMultiObjectiveWeightVectorChecker(preprocessorData)); - auto resultData = helper::SparseMultiObjectiveHelper::check(preprocessorData, weightVectorChecker); - swHelper.pause(); - STORM_LOG_DEBUG("Modelchecking done."); - - storm::utility::Stopwatch swPostprocessing; - result = helper::SparseMultiObjectivePostprocessor::postprocess(preprocessorData, resultData, swPreprocessing, swHelper, swPostprocessing); - STORM_LOG_DEBUG("Postprocessing done."); -#else - STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "Multi-objective model checking requires carl."); -#endif - return result; + return multiobjective::performPcaa(this->getModel(), checkTask.getFormula()); } diff --git a/src/modelchecker/multiobjective/helper/SparseMaMultiObjectiveWeightVectorChecker.cpp b/src/modelchecker/multiobjective/helper/SparseMaMultiObjectiveWeightVectorChecker.cpp index 6840c178f..b81fbba1c 100644 --- a/src/modelchecker/multiobjective/helper/SparseMaMultiObjectiveWeightVectorChecker.cpp +++ b/src/modelchecker/multiobjective/helper/SparseMaMultiObjectiveWeightVectorChecker.cpp @@ -303,17 +303,6 @@ namespace storm { std::unique_ptr::MinMaxSolverData> SparseMaMultiObjectiveWeightVectorChecker::initMinMaxSolverData(SubModel const& PS) const { std::unique_ptr result(new MinMaxSolverData()); - storm::storage::BitVector choicesStayingInPS(PS.getNumberOfChoices(), false); - for(uint_fast64_t choice = 0; choice < PS.toPS.getRowCount(); ++choice) { - if(storm::utility::isOne(PS.toPS.getRowSum(choice))) { - choicesStayingInPS.set(choice); - } - } - auto ecEliminatorResult = storm::transformer::EndComponentEliminator::transform(PS.toPS, choicesStayingInPS & storm::utility::vector::filterZero(PS.weightedRewardVector), storm::storage::BitVector(PS.getNumberOfStates(), true)); - result->matrix = std::move(ecEliminatorResult.matrix); - result->toPSChoiceMapping = std::move(ecEliminatorResult.newToOldRowMapping); - result->fromPSStateMapping = std::move(ecEliminatorResult.oldToNewStateMapping); - result->b.resize(result->matrix.getRowCount()); result->x.resize(result->matrix.getRowGroupCount()); for(uint_fast64_t state=0; state < result->fromPSStateMapping.size(); ++state) { diff --git a/src/modelchecker/multiobjective/helper/SparseMultiObjectiveWeightVectorChecker.cpp b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveWeightVectorChecker.cpp index 3c748b16e..e7d3fddf3 100644 --- a/src/modelchecker/multiobjective/helper/SparseMultiObjectiveWeightVectorChecker.cpp +++ b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveWeightVectorChecker.cpp @@ -105,7 +105,7 @@ namespace storm { } // Remove end components in which no reward is earned - auto ecEliminatorResult = storm::transformer::EndComponentEliminator::transform(data.preprocessedModel.getTransitionMatrix(), storm::utility::vector::filterZero(weightedRewardVector), storm::storage::BitVector(data.preprocessedModel.getTransitionMatrix().getRowGroupCount(), true)); + auto ecEliminatorResult = storm::transformer::EndComponentEliminator::transform(data.preprocessedModel.getTransitionMatrix(), storm::utility::vector::filterZero(weightedRewardVector), storm::storage::BitVector(data.preprocessedModel.getTransitionMatrix().getRowGroupCount(), true), storm::storage::BitVector(data.preprocessedModel.getTransitionMatrix().getRowGroupCount(), true)); std::vector subRewardVector(ecEliminatorResult.newToOldRowMapping.size()); storm::utility::vector::selectVectorValues(subRewardVector, ecEliminatorResult.newToOldRowMapping, weightedRewardVector); diff --git a/src/modelchecker/multiobjective/pcaa/PCAAObjective.h b/src/modelchecker/multiobjective/pcaa/PCAAObjective.h index c953d7062..5ff73230c 100644 --- a/src/modelchecker/multiobjective/pcaa/PCAAObjective.h +++ b/src/modelchecker/multiobjective/pcaa/PCAAObjective.h @@ -10,7 +10,7 @@ namespace storm { namespace modelchecker { namespace multiobjective { template - struct PCAAObjective { + struct PcaaObjective { // the original input formula std::shared_ptr originalFormula; @@ -35,8 +35,6 @@ namespace storm { boost::optional lowerTimeBound; boost::optional upperTimeBound; - bool rewardFinitenessChecked; - void printToStream(std::ostream& out) const { out << std::setw(30) << originalFormula->toString(); out << " \t(toOrigVal:" << std::setw(3) << toOriginalValueTransformationFactor << "*x +" << std::setw(3) << toOriginalValueTransformationOffset << ", \t"; diff --git a/src/modelchecker/multiobjective/pcaa/SparseMaPcaaWeightVectorChecker.cpp b/src/modelchecker/multiobjective/pcaa/SparseMaPcaaWeightVectorChecker.cpp new file mode 100644 index 000000000..324cd1396 --- /dev/null +++ b/src/modelchecker/multiobjective/pcaa/SparseMaPcaaWeightVectorChecker.cpp @@ -0,0 +1,423 @@ +#include "src/modelchecker/multiobjective/pcaa/SparseMaPcaaWeightVectorChecker.h" + +#include + +#include "src/adapters/CarlAdapter.h" +#include "src/models/sparse/MarkovAutomaton.h" +#include "src/models/sparse/StandardRewardModel.h" +#include "src/transformer/EndComponentEliminator.h" +#include "src/utility/macros.h" +#include "src/utility/vector.h" + +#include "src/exceptions/InvalidOperationException.h" + +namespace storm { + namespace modelchecker { + namespace multiobjective { + + template + SparseMaPcaaWeightVectorChecker::SparseMaPcaaWeightVectorChecker(SparseMaModelType const& model, + std::vector> const& objectives, + storm::storage::BitVector const& actionsWithNegativeReward, + storm::storage::BitVector const& ecActions, + storm::storage::BitVector const& possiblyRecurrentStates) : + SparsePcaaWeightVectorChecker(model, objectives, actionsWithNegativeReward, ecActions, possiblyRecurrentStates) { + // Set the (discretized) state action rewards. + this->discreteActionRewards.resize(objectives.size()); + for(auto objIndex : this->objectivesWithNoUpperTimeBound) { + typename SparseMaModelType::RewardModelType const& rewModel = this->model.getRewardModel(this->objectives[objIndex].rewardModelName); + STORM_LOG_ASSERT(!rewModel.hasTransitionRewards(), "Preprocessed Reward model has transition rewards which is not expected."); + this->discreteActionRewards[objIndex] = rewModel.hasStateActionRewards() ? rewModel.getStateActionRewardVector() : std::vector(this->model.getTransitionMatrix().getRowCount(), storm::utility::zero()); + if(rewModel.hasStateRewards()) { + // Note that state rewards are earned over time and thus play no role for probabilistic states + for(auto markovianState : this->model.getMarkovianStates()) { + this->discreteActionRewards[objIndex][this->model.getTransitionMatrix().getRowGroupIndices()[markovianState]] += rewModel.getStateReward(markovianState) / this->model.getExitRate(markovianState); + } + } + } + } + + template + void SparseMaPcaaWeightVectorChecker::boundedPhase(std::vector const& weightVector, std::vector& weightedRewardVector) { + + // Split the preprocessed model into transitions from/to probabilistic/Markovian states. + SubModel MS = createSubModel(true, weightedRewardVector); + SubModel PS = createSubModel(false, weightedRewardVector); + + // Apply digitization to Markovian transitions + ValueType digitizationConstant = getDigitizationConstant(); + digitize(MS, digitizationConstant); + + // Get for each occurring (digitized) timeBound the indices of the objectives with that bound. + TimeBoundMap lowerTimeBounds; + TimeBoundMap upperTimeBounds; + digitizeTimeBounds(lowerTimeBounds, upperTimeBounds, digitizationConstant); + + // Initialize a minMaxSolver to compute an optimal scheduler (w.r.t. PS) for each epoch + // No EC elimination is necessary as we assume non-zenoness + std::unique_ptr minMax = initMinMaxSolver(PS); + + // create a linear equation solver for the model induced by the optimal choice vector. + // the solver will be updated whenever the optimal choice vector has changed. + std::unique_ptr linEq = initLinEqSolver(PS); + + // Store the optimal choices of PS as computed by the minMax solver. + std::vector optimalChoicesAtCurrentEpoch(PS.getNumberOfStates(), std::numeric_limits::max()); + + // Stores the objectives for which we need to compute values in the current time epoch. + storm::storage::BitVector consideredObjectives = this->objectivesWithNoUpperTimeBound; + + auto lowerTimeBoundIt = lowerTimeBounds.begin(); + auto upperTimeBoundIt = upperTimeBounds.begin(); + uint_fast64_t currentEpoch = std::max(lowerTimeBounds.empty() ? 0 : lowerTimeBoundIt->first, upperTimeBounds.empty() ? 0 : upperTimeBoundIt->first); + while(true) { + // Update the objectives that are considered at the current time epoch as well as the (weighted) reward vectors. + updateDataToCurrentEpoch(MS, PS, *minMax, consideredObjectives, currentEpoch, weightVector, lowerTimeBoundIt, lowerTimeBounds, upperTimeBoundIt, upperTimeBounds); + + // Compute the values that can be obtained at probabilistic states in the current time epoch + performPSStep(PS, MS, *minMax, *linEq, optimalChoicesAtCurrentEpoch, consideredObjectives); + + // Compute values that can be obtained at Markovian states after letting one (digitized) time unit pass. + // Only perform such a step if there is time left. + if(currentEpoch>0) { + performMSStep(MS, PS, consideredObjectives); + --currentEpoch; + } else { + break; + } + } + + // compose the results from MS and PS + storm::utility::vector::setVectorValues(this->weightedResult, MS.states, MS.weightedSolutionVector); + storm::utility::vector::setVectorValues(this->weightedResult, PS.states, PS.weightedSolutionVector); + for(uint_fast64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) { + storm::utility::vector::setVectorValues(this->objectiveResults[objIndex], MS.states, MS.objectiveSolutionVectors[objIndex]); + storm::utility::vector::setVectorValues(this->objectiveResults[objIndex], PS.states, PS.objectiveSolutionVectors[objIndex]); + } + } + + + template + typename SparseMaPcaaWeightVectorChecker::SubModel SparseMaPcaaWeightVectorChecker::createSubModel(bool createMS, std::vector const& weightedRewardVector) const { + SubModel result; + + storm::storage::BitVector probabilisticStates = ~this->model.getMarkovianStates(); + result.states = createMS ? this->model.getMarkovianStates() : probabilisticStates; + result.choices = this->model.getTransitionMatrix().getRowIndicesOfRowGroups(result.states); + STORM_LOG_ASSERT(!createMS || result.states.getNumberOfSetBits() == result.choices.getNumberOfSetBits(), "row groups for Markovian states should consist of exactly one row"); + + //We need to add diagonal entries for selfloops on Markovian states. + result.toMS = this->model.getTransitionMatrix().getSubmatrix(true, result.states, this->model.getMarkovianStates(), createMS); + result.toPS = this->model.getTransitionMatrix().getSubmatrix(true, result.states, probabilisticStates, false); + STORM_LOG_ASSERT(result.getNumberOfStates() == result.states.getNumberOfSetBits() && result.getNumberOfStates() == result.toMS.getRowGroupCount() && result.getNumberOfStates() == result.toPS.getRowGroupCount(), "Invalid state count for subsystem"); + STORM_LOG_ASSERT(result.getNumberOfChoices() == result.choices.getNumberOfSetBits() && result.getNumberOfChoices() == result.toMS.getRowCount() && result.getNumberOfChoices() == result.toPS.getRowCount(), "Invalid state count for subsystem"); + + result.weightedRewardVector.resize(result.getNumberOfChoices()); + storm::utility::vector::selectVectorValues(result.weightedRewardVector, result.choices, weightedRewardVector); + result.objectiveRewardVectors.resize(this->objectives.size()); + for(uint_fast64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) { + std::vector& objVector = result.objectiveRewardVectors[objIndex]; + objVector = std::vector(result.weightedRewardVector.size(), storm::utility::zero()); + if(this->objectivesWithNoUpperTimeBound.get(objIndex)) { + storm::utility::vector::selectVectorValues(objVector, result.choices, this->discreteActionRewards[objIndex]); + } else { + typename SparseMaModelType::RewardModelType const& rewModel = this->model.getRewardModel(this->objectives[objIndex].rewardModelName); + STORM_LOG_ASSERT(!rewModel.hasTransitionRewards(), "Preprocessed Reward model has transition rewards which is not expected."); + STORM_LOG_ASSERT(!rewModel.hasStateRewards(), "State rewards for bounded objectives for MAs are not expected (bounded rewards are not supported)."); + if(rewModel.hasStateActionRewards()) { + storm::utility::vector::selectVectorValues(objVector, result.choices, rewModel.getStateActionRewardVector()); + } + } + } + + result.weightedSolutionVector.resize(result.getNumberOfStates()); + storm::utility::vector::selectVectorValues(result.weightedSolutionVector, result.states, this->weightedResult); + result.objectiveSolutionVectors.resize(this->objectives.size()); + for(uint_fast64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) { + result.objectiveSolutionVectors[objIndex].resize(result.weightedSolutionVector.size()); + storm::utility::vector::selectVectorValues(result.objectiveSolutionVectors[objIndex], result.states, this->objectiveResults[objIndex]); + } + + result.auxChoiceValues.resize(result.getNumberOfChoices()); + + return result; + } + + template + template ::SupportsExponential, int>::type> + VT SparseMaPcaaWeightVectorChecker::getDigitizationConstant() const { + STORM_LOG_DEBUG("Retrieving digitization constant"); + // We need to find a delta such that for each objective it holds that lowerbound/delta , upperbound/delta are natural numbers and + // If there is a lower and an upper bound: + // 1 - e^(-maxRate lowerbound) * (1 + maxRate delta) ^ (lowerbound / delta) + 1-e^(-maxRate upperbound) * (1 + maxRate delta) ^ (upperbound / delta) + (1-e^(-maxRate delta) <= maximumLowerUpperBoundGap + // If there is only an upper bound: + // 1-e^(-maxRate upperbound) * (1 + maxRate delta) ^ (upperbound / delta) <= maximumLowerUpperBoundGap + + // Initialize some data for fast and easy access + VT const maxRate = this->model.getMaximalExitRate(); + std::vector> eToPowerOfMinusMaxRateTimesBound; + VT smallestNonZeroBound = storm::utility::zero(); + for(auto const& obj : this->objectives) { + eToPowerOfMinusMaxRateTimesBound.emplace_back(); + if(obj.lowerTimeBound){ + STORM_LOG_ASSERT(!storm::utility::isZero(*obj.lowerTimeBound), "Got zero-valued lower bound."); // should have been handled in preprocessing + STORM_LOG_ASSERT(!obj.upperTimeBound || *obj.lowerTimeBound < *obj.upperTimeBound, "Got point intervall or empty intervall on time bounded property which is not supported"); // should also have been handled in preprocessing + eToPowerOfMinusMaxRateTimesBound.back().first = std::exp(-maxRate * (*obj.lowerTimeBound)); + smallestNonZeroBound = storm::utility::isZero(smallestNonZeroBound) ? *obj.lowerTimeBound : std::min(smallestNonZeroBound, *obj.lowerTimeBound); + } + if(obj.upperTimeBound){ + STORM_LOG_ASSERT(!storm::utility::isZero(*obj.upperTimeBound), "Got zero-valued upper bound."); // should have been handled in preprocessing + eToPowerOfMinusMaxRateTimesBound.back().second = std::exp(-maxRate * (*obj.upperTimeBound)); + smallestNonZeroBound = storm::utility::isZero(smallestNonZeroBound) ? *obj.upperTimeBound : std::min(smallestNonZeroBound, *obj.upperTimeBound); + } + } + if(storm::utility::isZero(smallestNonZeroBound)) { + // There are no time bounds. In this case, one is a valid digitization constant. + return storm::utility::one(); + } + + // We brute-force a delta, since a direct computation is apparently not easy. + // Also note that the number of times this loop runs is a lower bound for the number of minMaxSolver invocations. + // Hence, this brute-force approach will most likely not be a bottleneck. + uint_fast64_t smallestStepBound = 1; + VT delta = smallestNonZeroBound / smallestStepBound; + while(true) { + bool deltaValid = true; + for(auto const& obj : this->objectives) { + if((obj.lowerTimeBound && *obj.lowerTimeBound/delta != std::floor(*obj.lowerTimeBound/delta)) || + (obj.upperTimeBound && *obj.upperTimeBound/delta != std::floor(*obj.upperTimeBound/delta))) { + deltaValid = false; + break; + } + } + if(deltaValid) { + for(uint_fast64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) { + auto const& obj = this->objectives[objIndex]; + VT precisionOfObj = storm::utility::zero(); + if(obj.lowerTimeBound) { + precisionOfObj += storm::utility::one() - (eToPowerOfMinusMaxRateTimesBound[objIndex].first * storm::utility::pow(storm::utility::one() + maxRate * delta, *obj.lowerTimeBound / delta) ) + + storm::utility::one() - std::exp(-maxRate * delta); + } + if(obj.upperTimeBound) { + precisionOfObj += storm::utility::one() - (eToPowerOfMinusMaxRateTimesBound[objIndex].second * storm::utility::pow(storm::utility::one() + maxRate * delta, *obj.upperTimeBound / delta) ); + } + if(precisionOfObj > this->maximumLowerUpperBoundGap) { + deltaValid = false; + break; + } + } + } + if(deltaValid) { + break; + } + ++smallestStepBound; + STORM_LOG_ASSERT(delta>smallestNonZeroBound / smallestStepBound, "Digitization constant is expected to become smaller in every iteration"); + delta = smallestNonZeroBound / smallestStepBound; + } + STORM_LOG_DEBUG("Found digitization constant: " << delta << ". At least " << smallestStepBound << " digitization steps will be necessarry"); + return delta; + } + + template + template ::SupportsExponential, int>::type> + VT SparseMaPcaaWeightVectorChecker::getDigitizationConstant() const { + STORM_LOG_THROW(false, storm::exceptions::InvalidOperationException, "Computing bounded probabilities of MAs is unsupported for this value type."); + } + + template + template ::SupportsExponential, int>::type> + void SparseMaPcaaWeightVectorChecker::digitize(SubModel& MS, VT const& digitizationConstant) const { + std::vector rateVector(MS.getNumberOfChoices()); + storm::utility::vector::selectVectorValues(rateVector, MS.states, this->model.getExitRates()); + for(uint_fast64_t row = 0; row < rateVector.size(); ++row) { + VT const eToMinusRateTimesDelta = std::exp(-rateVector[row] * digitizationConstant); + for(auto& entry : MS.toMS.getRow(row)) { + entry.setValue((storm::utility::one() - eToMinusRateTimesDelta) * entry.getValue()); + if(entry.getColumn() == row) { + entry.setValue(entry.getValue() + eToMinusRateTimesDelta); + } + } + for(auto& entry : MS.toPS.getRow(row)) { + entry.setValue((storm::utility::one() - eToMinusRateTimesDelta) * entry.getValue()); + } + MS.weightedRewardVector[row] *= storm::utility::one() - eToMinusRateTimesDelta; + for(auto& objVector : MS.objectiveRewardVectors) { + objVector[row] *= storm::utility::one() - eToMinusRateTimesDelta; + } + } + } + + template + template ::SupportsExponential, int>::type> + void SparseMaPcaaWeightVectorChecker::digitize(SubModel& subModel, VT const& digitizationConstant) const { + STORM_LOG_THROW(false, storm::exceptions::InvalidOperationException, "Computing bounded probabilities of MAs is unsupported for this value type."); + } + + template + template ::SupportsExponential, int>::type> + void SparseMaPcaaWeightVectorChecker::digitizeTimeBounds(TimeBoundMap& lowerTimeBounds, TimeBoundMap& upperTimeBounds, VT const& digitizationConstant) { + + VT const maxRate = this->model.getMaximalExitRate(); + for(uint_fast64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) { + auto const& obj = this->objectives[objIndex]; + VT errorTowardsZero; + VT errorAwayFromZero; + if(obj.lowerTimeBound) { + uint_fast64_t digitizedBound = storm::utility::convertNumber((*obj.lowerTimeBound)/digitizationConstant); + auto timeBoundIt = lowerTimeBounds.insert(std::make_pair(digitizedBound, storm::storage::BitVector(this->objectives.size(), false))).first; + timeBoundIt->second.set(objIndex); + VT digitizationError = storm::utility::one(); + digitizationError -= std::exp(-maxRate * (*obj.lowerTimeBound)) * storm::utility::pow(storm::utility::one() + maxRate * digitizationConstant, digitizedBound); + errorTowardsZero = -digitizationError; + errorAwayFromZero = storm::utility::one() - std::exp(-maxRate * digitizationConstant);; + } else { + errorTowardsZero = storm::utility::zero(); + errorAwayFromZero = storm::utility::zero(); + } + if(obj.upperTimeBound) { + uint_fast64_t digitizedBound = storm::utility::convertNumber((*obj.upperTimeBound)/digitizationConstant); + auto timeBoundIt = upperTimeBounds.insert(std::make_pair(digitizedBound, storm::storage::BitVector(this->objectives.size(), false))).first; + timeBoundIt->second.set(objIndex); + VT digitizationError = storm::utility::one(); + digitizationError -= std::exp(-maxRate * (*obj.upperTimeBound)) * storm::utility::pow(storm::utility::one() + maxRate * digitizationConstant, digitizedBound); + errorAwayFromZero += digitizationError; + } + STORM_LOG_ASSERT(errorTowardsZero + errorAwayFromZero <= this->maximumLowerUpperBoundGap, "Precision not sufficient."); + if (obj.rewardsArePositive) { + this->offsetsToLowerBound[objIndex] = -errorTowardsZero; + this->offsetsToUpperBound[objIndex] = errorAwayFromZero; + } else { + this->offsetsToLowerBound[objIndex] = -errorAwayFromZero; + this->offsetsToUpperBound[objIndex] = errorTowardsZero; + } + } + } + + template + template ::SupportsExponential, int>::type> + void SparseMaPcaaWeightVectorChecker::digitizeTimeBounds(TimeBoundMap& lowerTimeBounds, TimeBoundMap& upperTimeBounds, VT const& digitizationConstant) { + STORM_LOG_THROW(false, storm::exceptions::InvalidOperationException, "Computing bounded probabilities of MAs is unsupported for this value type."); + } + + template + std::unique_ptr::MinMaxSolverData> SparseMaPcaaWeightVectorChecker::initMinMaxSolver(SubModel const& PS) const { + std::unique_ptr result(new MinMaxSolverData()); + storm::solver::GeneralMinMaxLinearEquationSolverFactory minMaxSolverFactory; + result->solver = minMaxSolverFactory.create(PS.toPS); + result->solver->setOptimizationDirection(storm::solver::OptimizationDirection::Maximize); + result->solver->setTrackScheduler(true); + result->solver->allocateAuxMemory(storm::solver::MinMaxLinearEquationSolverOperation::SolveEquations); + + result->b.resize(PS.getNumberOfChoices()); + + return result; + } + + template + std::unique_ptr::LinEqSolverData> SparseMaPcaaWeightVectorChecker::initLinEqSolver(SubModel const& PS) const { + std::unique_ptr result(new LinEqSolverData()); + // We choose Jacobi since we call the solver very frequently on 'easy' inputs (note that jacobi without preconditioning has very little overhead). + result->factory.getSettings().setSolutionMethod(storm::solver::GmmxxLinearEquationSolverSettings::SolutionMethod::Jacobi); + result->factory.getSettings().setPreconditioner(storm::solver::GmmxxLinearEquationSolverSettings::Preconditioner::None); + result->b.resize(PS.getNumberOfStates()); + return result; + } + + template + void SparseMaPcaaWeightVectorChecker::updateDataToCurrentEpoch(SubModel& MS, SubModel& PS, MinMaxSolverData& minMax, storm::storage::BitVector& consideredObjectives, uint_fast64_t const& currentEpoch, std::vector const& weightVector, TimeBoundMap::iterator& lowerTimeBoundIt, TimeBoundMap const& lowerTimeBounds, TimeBoundMap::iterator& upperTimeBoundIt, TimeBoundMap const& upperTimeBounds) { + + //Note that lower time bounds are always strict. Hence, we need to react when the current epoch equals the stored bound. + if(lowerTimeBoundIt != lowerTimeBounds.end() && currentEpoch == lowerTimeBoundIt->first) { + for(auto objIndex : lowerTimeBoundIt->second) { + // No more reward is earned for this objective. + storm::utility::vector::addScaledVector(MS.weightedRewardVector, MS.objectiveRewardVectors[objIndex], -weightVector[objIndex]); + storm::utility::vector::addScaledVector(PS.weightedRewardVector, PS.objectiveRewardVectors[objIndex], -weightVector[objIndex]); + MS.objectiveRewardVectors[objIndex] = std::vector(MS.objectiveRewardVectors[objIndex].size(), storm::utility::zero()); + PS.objectiveRewardVectors[objIndex] = std::vector(PS.objectiveRewardVectors[objIndex].size(), storm::utility::zero()); + } + ++lowerTimeBoundIt; + } + + if(upperTimeBoundIt != upperTimeBounds.end() && currentEpoch == upperTimeBoundIt->first) { + consideredObjectives |= upperTimeBoundIt->second; + for(auto objIndex : upperTimeBoundIt->second) { + // This objective now plays a role in the weighted sum + storm::utility::vector::addScaledVector(MS.weightedRewardVector, MS.objectiveRewardVectors[objIndex], weightVector[objIndex]); + storm::utility::vector::addScaledVector(PS.weightedRewardVector, PS.objectiveRewardVectors[objIndex], weightVector[objIndex]); + } + ++upperTimeBoundIt; + } + + // Update the solver data + PS.toMS.multiplyWithVector(MS.weightedSolutionVector, minMax.b); + storm::utility::vector::addVectors(minMax.b, PS.weightedRewardVector, minMax.b); + } + + template + void SparseMaPcaaWeightVectorChecker::performPSStep(SubModel& PS, SubModel const& MS, MinMaxSolverData& minMax, LinEqSolverData& linEq, std::vector& optimalChoicesAtCurrentEpoch, storm::storage::BitVector const& consideredObjectives) const { + // compute a choice vector for the probabilistic states that is optimal w.r.t. the weighted reward vector + minMax.solver->solveEquations(PS.weightedSolutionVector, minMax.b); + auto newScheduler = minMax.solver->getScheduler(); + // check whether the linEqSolver needs to be updated, i.e., whether the scheduler has changed + if(linEq.solver == nullptr || newScheduler->getChoices() != optimalChoicesAtCurrentEpoch) { + optimalChoicesAtCurrentEpoch = newScheduler->getChoices(); + linEq.solver = nullptr; + storm::storage::SparseMatrix linEqMatrix = PS.toPS.selectRowsFromRowGroups(optimalChoicesAtCurrentEpoch, true); + linEqMatrix.convertToEquationSystem(); + linEq.solver = linEq.factory.create(std::move(linEqMatrix)); + linEq.solver->allocateAuxMemory(storm::solver::LinearEquationSolverOperation::SolveEquations); + } + + // Get the results for the individual objectives. + // Note that we do not consider an estimate for each objective (as done in the unbounded phase) since the results from the previous epoch are already pretty close + for(auto objIndex : consideredObjectives) { + auto const& objectiveRewardVectorPS = PS.objectiveRewardVectors[objIndex]; + auto const& objectiveSolutionVectorMS = MS.objectiveSolutionVectors[objIndex]; + // compute rhs of equation system, i.e., PS.toMS * x + Rewards + // To safe some time, only do this for the obtained optimal choices + auto itGroupIndex = PS.toPS.getRowGroupIndices().begin(); + auto itChoiceOffset = optimalChoicesAtCurrentEpoch.begin(); + for(auto& bValue : linEq.b) { + uint_fast64_t row = (*itGroupIndex) + (*itChoiceOffset); + bValue = objectiveRewardVectorPS[row]; + for(auto const& entry : PS.toMS.getRow(row)){ + bValue += entry.getValue() * objectiveSolutionVectorMS[entry.getColumn()]; + } + ++itGroupIndex; + ++itChoiceOffset; + } + linEq.solver->solveEquations(PS.objectiveSolutionVectors[objIndex], linEq.b); + } + } + + template + void SparseMaPcaaWeightVectorChecker::performMSStep(SubModel& MS, SubModel const& PS, storm::storage::BitVector const& consideredObjectives) const { + + MS.toMS.multiplyWithVector(MS.weightedSolutionVector, MS.auxChoiceValues); + storm::utility::vector::addVectors(MS.weightedRewardVector, MS.auxChoiceValues, MS.weightedSolutionVector); + MS.toPS.multiplyWithVector(PS.weightedSolutionVector, MS.auxChoiceValues); + storm::utility::vector::addVectors(MS.weightedSolutionVector, MS.auxChoiceValues, MS.weightedSolutionVector); + + for(auto objIndex : consideredObjectives) { + MS.toMS.multiplyWithVector(MS.objectiveSolutionVectors[objIndex], MS.auxChoiceValues); + storm::utility::vector::addVectors(MS.objectiveRewardVectors[objIndex], MS.auxChoiceValues, MS.objectiveSolutionVectors[objIndex]); + MS.toPS.multiplyWithVector(PS.objectiveSolutionVectors[objIndex], MS.auxChoiceValues); + storm::utility::vector::addVectors(MS.objectiveSolutionVectors[objIndex], MS.auxChoiceValues, MS.objectiveSolutionVectors[objIndex]); + } + } + + + template class SparseMaPcaaWeightVectorChecker>; + template double SparseMaPcaaWeightVectorChecker>::getDigitizationConstant() const; + template void SparseMaPcaaWeightVectorChecker>::digitize(SparseMaPcaaWeightVectorChecker>::SubModel& subModel, double const& digitizationConstant) const; + template void SparseMaPcaaWeightVectorChecker>::digitizeTimeBounds(SparseMaPcaaWeightVectorChecker>::TimeBoundMap& lowerTimeBounds, SparseMaPcaaWeightVectorChecker>::TimeBoundMap& upperTimeBounds, double const& digitizationConstant); +#ifdef STORM_HAVE_CARL +// template class SparseMaPcaaWeightVectorChecker>; + // template storm::RationalNumber SparseMaPcaaWeightVectorChecker>::getDigitizationConstant() const; + // template void SparseMaPcaaWeightVectorChecker>::digitize(SparseMaPcaaWeightVectorChecker>::SubModel& subModel, storm::RationalNumber const& digitizationConstant) const; +// template void SparseMaPcaaWeightVectorChecker>::digitizeTimeBounds(SparseMaPcaaWeightVectorChecker>::TimeBoundMap& lowerTimeBounds, SparseMaPcaaWeightVectorChecker>::TimeBoundMap& upperTimeBounds, storm::RationalNumber const& digitizationConstant); +#endif + + } + } +} diff --git a/src/modelchecker/multiobjective/pcaa/SparseMaPcaaWeightVectorChecker.h b/src/modelchecker/multiobjective/pcaa/SparseMaPcaaWeightVectorChecker.h new file mode 100644 index 000000000..24166e023 --- /dev/null +++ b/src/modelchecker/multiobjective/pcaa/SparseMaPcaaWeightVectorChecker.h @@ -0,0 +1,157 @@ +#ifndef STORM_MODELCHECKER_MULTIOBJECTIVE_PCAA_SPARSEMAPCAAWEIGHTVECTORCHECKER_H_ +#define STORM_MODELCHECKER_MULTIOBJECTIVE_PCAA_SPARSEMAPCAAWEIGHTVECTORCHECKER_H_ + +#include +#include + +#include "src/modelchecker/multiobjective/pcaa/SparsePcaaWeightVectorChecker.h" +#include "src/solver/LinearEquationSolver.h" +#include "src/solver/GmmxxLinearEquationSolver.h" +#include "src/solver/MinMaxLinearEquationSolver.h" +#include "src/utility/NumberTraits.h" + +namespace storm { + namespace modelchecker { + namespace multiobjective { + + /*! + * Helper Class that takes preprocessed Pcaa data and a weight vector and ... + * - computes the maximal expected reward w.r.t. the weighted sum of the rewards of the individual objectives + * - extracts the scheduler that induces this maximum + * - computes for each objective the value induced by this scheduler + */ + template + class SparseMaPcaaWeightVectorChecker : public SparsePcaaWeightVectorChecker { + public: + typedef typename SparseMaModelType::ValueType ValueType; + + SparseMaPcaaWeightVectorChecker(SparseMaModelType const& model, + std::vector> const& objectives, + storm::storage::BitVector const& actionsWithNegativeReward, + storm::storage::BitVector const& ecActions, + storm::storage::BitVector const& possiblyRecurrentStates); + + private: + + /* + * Stores (digitized) time bounds in descending order + */ + typedef std::map> TimeBoundMap; + + /* + * Stores the ingredients of a sub model + */ + struct SubModel { + storm::storage::BitVector states; // The states that are part of this sub model + storm::storage::BitVector choices; // The choices that are part of this sub model + + storm::storage::SparseMatrix toMS; // Transitions to Markovian states + storm::storage::SparseMatrix toPS; // Transitions to probabilistic states + + std::vector weightedRewardVector; + std::vector> objectiveRewardVectors; + + std::vector weightedSolutionVector; + std::vector> objectiveSolutionVectors; + + std::vector auxChoiceValues; //stores auxiliary values for every choice + + uint_fast64_t getNumberOfStates() const { return toMS.getRowGroupCount(); }; + uint_fast64_t getNumberOfChoices() const { return toMS.getRowCount(); }; + }; + + /* + * Stores the data that is relevant to invoke the minMaxSolver and retrieve the result. + */ + struct MinMaxSolverData { + std::unique_ptr> solver; + std::vector b; + }; + + struct LinEqSolverData { + storm::solver::GmmxxLinearEquationSolverFactory factory; + std::unique_ptr> solver; + std::vector b; + }; + + /*! + * + * @param weightVector the weight vector of the current check + * @param weightedRewardVector the weighted rewards considering the unbounded objectives. Will be invalidated after calling this. + */ + virtual void boundedPhase(std::vector const& weightVector, std::vector& weightedRewardVector) override; + + /*! + * Retrieves the data for a submodel of the data->preprocessedModel + * @param createMS if true, the submodel containing the Markovian states is created. + * if false, the submodel containing the probabilistic states is created. + */ + SubModel createSubModel(bool createMS, std::vector const& weightedRewardVector) const; + + /*! + * Retrieves the delta used for digitization + */ + template ::SupportsExponential, int>::type = 0> + VT getDigitizationConstant() const; + template ::SupportsExponential, int>::type = 0> + VT getDigitizationConstant() const; + + /*! + * Digitizes the given matrix and vectors w.r.t. the given digitization constant and the given rate vector. + */ + template ::SupportsExponential, int>::type = 0> + void digitize(SubModel& subModel, VT const& digitizationConstant) const; + template ::SupportsExponential, int>::type = 0> + void digitize(SubModel& subModel, VT const& digitizationConstant) const; + + /* + * Fills the given maps with the digitized time bounds. Also sets the offsetsToLowerBound / offsetsToUpperBound values + * according to the digitization error + */ + template ::SupportsExponential, int>::type = 0> + void digitizeTimeBounds(TimeBoundMap& lowerTimeBounds, TimeBoundMap& upperTimeBounds, VT const& digitizationConstant); + template ::SupportsExponential, int>::type = 0> + void digitizeTimeBounds(TimeBoundMap& lowerTimeBounds, TimeBoundMap& upperTimeBounds, VT const& digitizationConstant); + + + /*! + * Initializes the data for the MinMax solver + */ + std::unique_ptr initMinMaxSolver(SubModel const& PS) const; + + /*! + * Initializes the data for the LinEq solver + */ + std::unique_ptr initLinEqSolver(SubModel const& PS) const; + + /* + * Updates the reward vectors within the split model, + * the reward vector of the reduced PStoPS model, and + * objectives that are considered at the current time epoch. + */ + void updateDataToCurrentEpoch(SubModel& MS, SubModel& PS, MinMaxSolverData& minMax, storm::storage::BitVector& consideredObjectives, uint_fast64_t const& currentEpoch, std::vector const& weightVector, TimeBoundMap::iterator& lowerTimeBoundIt, TimeBoundMap const& lowerTimeBounds, TimeBoundMap::iterator& upperTimeBoundIt, TimeBoundMap const& upperTimeBounds); + + /* + * Performs a step for the probabilistic states, that is + * * Compute an optimal scheduler for the weighted reward sum + * * Compute the values for the individual objectives w.r.t. that scheduler + * + * The resulting values represent the rewards at probabilistic states that are obtained at the current time epoch. + */ + void performPSStep(SubModel& PS, SubModel const& MS, MinMaxSolverData& minMax, LinEqSolverData& linEq, std::vector& optimalChoicesAtCurrentEpoch, storm::storage::BitVector const& consideredObjectives) const; + + /* + * Performs a step for the Markovian states, that is + * * Compute values for the weighted reward sum as well as for the individual objectives + * + * The resulting values represent the rewards at Markovian states that are obtained after one (digitized) time unit has passed. + */ + void performMSStep(SubModel& MS, SubModel const& PS, storm::storage::BitVector const& consideredObjectives) const; + + }; + + } + } +} + +#endif /* STORM_MODELCHECKER_MULTIOBJECTIVE_PCAA_SPARSEMAPCAAWEIGHTVECTORCHECKER_H_ */ diff --git a/src/modelchecker/multiobjective/pcaa/SparseMdpPcaaWeightVectorChecker.cpp b/src/modelchecker/multiobjective/pcaa/SparseMdpPcaaWeightVectorChecker.cpp new file mode 100644 index 000000000..9529bf890 --- /dev/null +++ b/src/modelchecker/multiobjective/pcaa/SparseMdpPcaaWeightVectorChecker.cpp @@ -0,0 +1,122 @@ +#include "src/modelchecker/multiobjective/pcaa/SparseMdpPcaaWeightVectorChecker.h" + +#include "src/adapters/CarlAdapter.h" +#include "src/models/sparse/Mdp.h" +#include "src/models/sparse/StandardRewardModel.h" +#include "src/utility/macros.h" +#include "src/utility/vector.h" + + +namespace storm { + namespace modelchecker { + namespace multiobjective { + + template + SparseMdpPcaaWeightVectorChecker::SparseMdpPcaaWeightVectorChecker(SparseMdpModelType const& model, + std::vector> const& objectives, + storm::storage::BitVector const& actionsWithNegativeReward, + storm::storage::BitVector const& ecActions, + storm::storage::BitVector const& possiblyRecurrentStates) : + SparsePcaaWeightVectorChecker(model, objectives, actionsWithNegativeReward, ecActions, possiblyRecurrentStates) { + // set the state action rewards + for(uint_fast64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) { + typename SparseMdpModelType::RewardModelType const& rewModel = this->model.getRewardModel(this->objectives[objIndex].rewardModelName); + STORM_LOG_ASSERT(!rewModel.hasTransitionRewards(), "Reward model has transition rewards which is not expected."); + this->discreteActionRewards[objIndex] = rewModel.getTotalRewardVector(this->model.getTransitionMatrix()); + } + } + + template + void SparseMdpPcaaWeightVectorChecker::boundedPhase(std::vector const& weightVector, std::vector& weightedRewardVector) { + // Allocate some memory so this does not need to happen for each time epoch + std::vector optimalChoicesInCurrentEpoch(this->model.getNumberOfStates()); + std::vector choiceValues(weightedRewardVector.size()); + std::vector temporaryResult(this->model.getNumberOfStates()); + std::vector zeroReward(weightedRewardVector.size(), storm::utility::zero()); + // Get for each occurring timeBound the indices of the objectives with that bound. + std::map> lowerTimeBounds; + std::map> upperTimeBounds; + for(uint_fast64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) { + auto const& obj = this->objectives[objIndex]; + if(obj.lowerTimeBound) { + auto timeBoundIt = lowerTimeBounds.insert(std::make_pair(storm::utility::convertNumber(*obj.lowerTimeBound), storm::storage::BitVector(this->objectives.size(), false))).first; + STORM_LOG_WARN_COND(storm::utility::convertNumber(timeBoundIt->first) == (*obj.lowerTimeBound), "Rounded non-integral bound " << *obj.lowerTimeBound << " to " << timeBoundIt->first << "."); + timeBoundIt->second.set(objIndex); + } + if(obj.upperTimeBound) { + auto timeBoundIt = upperTimeBounds.insert(std::make_pair(storm::utility::convertNumber(*obj.upperTimeBound), storm::storage::BitVector(this->objectives.size(), false))).first; + STORM_LOG_WARN_COND(storm::utility::convertNumber(timeBoundIt->first) == (*obj.upperTimeBound), "Rounded non-integral bound " << *obj.upperTimeBound << " to " << timeBoundIt->first << "."); + timeBoundIt->second.set(objIndex); + + // There is no error for the values of these objectives. + this->offsetsToLowerBound[objIndex] = storm::utility::zero(); + this->offsetsToUpperBound[objIndex] = storm::utility::zero(); + } + } + + // Stores the objectives for which we need to compute values in the current time epoch. + storm::storage::BitVector consideredObjectives = this->objectivesWithNoUpperTimeBound; + + // Stores objectives for which the current epoch passed their lower bound + storm::storage::BitVector lowerBoundViolatedObjectives(consideredObjectives.size(), false); + + auto lowerTimeBoundIt = lowerTimeBounds.begin(); + auto upperTimeBoundIt = upperTimeBounds.begin(); + uint_fast64_t currentEpoch = std::max(lowerTimeBounds.empty() ? 0 : lowerTimeBoundIt->first - 1, upperTimeBounds.empty() ? 0 : upperTimeBoundIt->first); // consider lowerBound - 1 since we are interested in the first epoch that passes the bound + + while(currentEpoch > 0) { + //For lower time bounds we need to react when the currentEpoch passed the bound + // Hence, we substract 1 from the lower time bounds. + if(lowerTimeBoundIt != lowerTimeBounds.end() && currentEpoch == lowerTimeBoundIt->first - 1) { + lowerBoundViolatedObjectives |= lowerTimeBoundIt->second; + for(auto objIndex : lowerTimeBoundIt->second) { + // No more reward is earned for this objective. + storm::utility::vector::addScaledVector(weightedRewardVector, this->discreteActionRewards[objIndex], -weightVector[objIndex]); + } + ++lowerTimeBoundIt; + } + + if(upperTimeBoundIt != upperTimeBounds.end() && currentEpoch == upperTimeBoundIt->first) { + consideredObjectives |= upperTimeBoundIt->second; + for(auto objIndex : upperTimeBoundIt->second) { + // This objective now plays a role in the weighted sum + storm::utility::vector::addScaledVector(weightedRewardVector, this->discreteActionRewards[objIndex], weightVector[objIndex]); + } + ++upperTimeBoundIt; + } + + // Get values and scheduler for weighted sum of objectives + this->model.getTransitionMatrix().multiplyWithVector(this->weightedResult, choiceValues); + storm::utility::vector::addVectors(choiceValues, weightedRewardVector, choiceValues); + storm::utility::vector::reduceVectorMax(choiceValues, this->weightedResult, this->model.getTransitionMatrix().getRowGroupIndices(), &optimalChoicesInCurrentEpoch); + + // get values for individual objectives + // TODO we could compute the result for one of the objectives from the weighted result, the given weight vector, and the remaining objective results. + for(auto objIndex : consideredObjectives) { + std::vector& objectiveResult = this->objectiveResults[objIndex]; + std::vector const& objectiveRewards = lowerBoundViolatedObjectives.get(objIndex) ? zeroReward : this->discreteActionRewards[objIndex]; + auto rowGroupIndexIt = this->model.getTransitionMatrix().getRowGroupIndices().begin(); + auto optimalChoiceIt = optimalChoicesInCurrentEpoch.begin(); + for(ValueType& stateValue : temporaryResult){ + uint_fast64_t row = (*rowGroupIndexIt) + (*optimalChoiceIt); + ++rowGroupIndexIt; + ++optimalChoiceIt; + stateValue = objectiveRewards[row]; + for(auto const& entry : this->model.getTransitionMatrix().getRow(row)) { + stateValue += entry.getValue() * objectiveResult[entry.getColumn()]; + } + } + objectiveResult.swap(temporaryResult); + } + --currentEpoch; + } + } + + template class SparseMdpPcaaWeightVectorChecker>; +#ifdef STORM_HAVE_CARL + template class SparseMdpPcaaWeightVectorChecker>; +#endif + + } + } +} diff --git a/src/modelchecker/multiobjective/pcaa/SparseMdpPcaaWeightVectorChecker.h b/src/modelchecker/multiobjective/pcaa/SparseMdpPcaaWeightVectorChecker.h new file mode 100644 index 000000000..e6565f413 --- /dev/null +++ b/src/modelchecker/multiobjective/pcaa/SparseMdpPcaaWeightVectorChecker.h @@ -0,0 +1,48 @@ +#ifndef STORM_MODELCHECKER_MULTIOBJECTIVE_PCAA_SPARSEMDPPCAAWEIGHTVECTORCHECKER_H_ +#define STORM_MODELCHECKER_MULTIOBJECTIVE_PCAA_SPARSEMDPPCAAWEIGHTVECTORCHECKER_H_ + +#include + +#include "src/modelchecker/multiobjective/pcaa/SparsePcaaWeightVectorChecker.h" + +namespace storm { + namespace modelchecker { + namespace multiobjective { + + /*! + * Helper Class that takes preprocessed Pcaa data and a weight vector and ... + * - computes the maximal expected reward w.r.t. the weighted sum of the rewards of the individual objectives + * - extracts the scheduler that induces this maximum + * - computes for each objective the value induced by this scheduler + */ + template + class SparseMdpPcaaWeightVectorChecker : public SparsePcaaWeightVectorChecker { + public: + typedef typename SparseMdpModelType::ValueType ValueType; + + SparseMdpPcaaWeightVectorChecker(SparseMdpModelType const& model, + std::vector> const& objectives, + storm::storage::BitVector const& actionsWithNegativeReward, + storm::storage::BitVector const& ecActions, + storm::storage::BitVector const& possiblyRecurrentStates); + + private: + + /*! + * For each time epoch (starting with the maximal stepBound occurring in the objectives), this method + * - determines the objectives that are relevant in the current time epoch + * - determines the maximizing scheduler for the weighted reward vector of these objectives + * - computes the values of these objectives w.r.t. this scheduler + * + * @param weightVector the weight vector of the current check + * @param weightedRewardVector the weighted rewards considering the unbounded objectives. Will be invalidated after calling this. + */ + virtual void boundedPhase(std::vector const& weightVector, std::vector& weightedRewardVector) override; + + }; + + } + } +} + +#endif /* STORM_MODELCHECKER_MULTIOBJECTIVE_PCAA_SPARSEMDPPCAAWEIGHTVECTORCHECKER_H_ */ diff --git a/src/modelchecker/multiobjective/pcaa/SparsePCAAPreprocessor.cpp b/src/modelchecker/multiobjective/pcaa/SparsePCAAPreprocessor.cpp index 4f33fc39a..3201de777 100644 --- a/src/modelchecker/multiobjective/pcaa/SparsePCAAPreprocessor.cpp +++ b/src/modelchecker/multiobjective/pcaa/SparsePCAAPreprocessor.cpp @@ -1,4 +1,4 @@ - #include "src/modelchecker/multiobjective/pcaa/SparsePCAAPreprocessor.h" + #include "src/modelchecker/multiobjective/pcaa/SparsePcaaPreprocessor.h" #include "src/models/sparse/Mdp.h" #include "src/models/sparse/MarkovAutomaton.h" @@ -20,7 +20,7 @@ namespace storm { namespace multiobjective { template - typename SparsePCAAPreprocessor::ReturnType SparsePCAAPreprocessor::preprocess(storm::logic::MultiObjectiveFormula const& originalFormula, SparseModelType const& originalModel) { + typename SparsePcaaPreprocessor::ReturnType SparsePcaaPreprocessor::preprocess(SparseModelType const& originalModel, storm::logic::MultiObjectiveFormula const& originalFormula) { ReturnType result(originalFormula, originalModel, SparseModelType(originalModel)); result.newToOldStateIndexMapping = storm::utility::vector::buildVectorForRange(0, originalModel.getNumberOfStates()); @@ -29,7 +29,7 @@ namespace storm { for(auto const& subFormula : originalFormula.getSubformulas()){ STORM_LOG_DEBUG("Preprocessing objective " << *subFormula<< "."); result.objectives.emplace_back(); - PCAAObjective& currentObjective = result.objectives.back(); + PcaaObjective& currentObjective = result.objectives.back(); currentObjective.originalFormula = subFormula; if(currentObjective.originalFormula->isOperatorFormula()) { preprocessFormula(currentObjective.originalFormula->asOperatorFormula(), result, currentObjective); @@ -62,13 +62,26 @@ namespace storm { result.preprocessedModel.removeRewardModel(rewModel); } - ensureRewardFiniteness(result); + //Get actions to which a positive or negative reward is assigned for at least one objective + result.actionsWithPositiveReward = storm::storage::BitVector(result.preprocessedModel.getNumberOfChoices(), false); + result.actionsWithNegativeReward = storm::storage::BitVector(result.preprocessedModel.getNumberOfChoices(), false); + for(uint_fast64_t objIndex = 0; objIndex < result.objectives.size(); ++objIndex) { + if(result.objectives[objIndex].rewardsArePositive) { + result.actionsWithPositiveReward |= ~storm::utility::vector::filterZero(result.preprocessedModel.getRewardModel(result.objectives[objIndex].rewardModelName).getTotalRewardVector(result.preprocessedModel.getTransitionMatrix())); + } else { + result.actionsWithNegativeReward |= ~storm::utility::vector::filterZero(result.preprocessedModel.getRewardModel(result.objectives[objIndex].rewardModelName).getTotalRewardVector(result.preprocessedModel.getTransitionMatrix())); + } + } + + auto backwardTransitions = result.preprocessedModel.getBackwardTransitions(); + analyzeEndComponents(result, backwardTransitions); + ensureRewardFiniteness(result, backwardTransitions); return result; } template - void SparsePCAAPreprocessor::updatePreprocessedModel(ReturnType& result, SparseModelType& newPreprocessedModel, std::vector& newToOldStateIndexMapping) { + void SparsePcaaPreprocessor::updatePreprocessedModel(ReturnType& result, SparseModelType& newPreprocessedModel, std::vector& newToOldStateIndexMapping) { result.preprocessedModel = std::move(newPreprocessedModel); // the given newToOldStateIndexMapping reffers to the indices of the former preprocessedModel as 'old indices' for(auto & preprocessedModelStateIndex : newToOldStateIndexMapping){ @@ -78,7 +91,7 @@ namespace storm { } template - void SparsePCAAPreprocessor::preprocessFormula(storm::logic::OperatorFormula const& formula, ReturnType& result, PCAAObjective& currentObjective) { + void SparsePcaaPreprocessor::preprocessFormula(storm::logic::OperatorFormula const& formula, ReturnType& result, PcaaObjective& currentObjective) { // Get a unique name for the new reward model. currentObjective.rewardModelName = "objective" + std::to_string(result.objectives.size()); @@ -125,8 +138,7 @@ namespace storm { } template - void SparsePCAAPreprocessor::preprocessFormula(storm::logic::ProbabilityOperatorFormula const& formula, ReturnType& result, PCAAObjective& currentObjective) { - currentObjective.rewardFinitenessChecked = true; + void SparsePcaaPreprocessor::preprocessFormula(storm::logic::ProbabilityOperatorFormula const& formula, ReturnType& result, PcaaObjective& currentObjective) { if(formula.getSubformula().isUntilFormula()){ preprocessFormula(formula.getSubformula().asUntilFormula(), result, currentObjective); @@ -142,16 +154,13 @@ namespace storm { } template - void SparsePCAAPreprocessor::preprocessFormula(storm::logic::RewardOperatorFormula const& formula, ReturnType& result, PCAAObjective& currentObjective) { + void SparsePcaaPreprocessor::preprocessFormula(storm::logic::RewardOperatorFormula const& formula, ReturnType& result, PcaaObjective& currentObjective) { // Check if the reward model is uniquely specified STORM_LOG_THROW((formula.hasRewardModelName() && result.preprocessedModel.hasRewardModel(formula.getRewardModelName())) || result.preprocessedModel.hasUniqueRewardModel(), storm::exceptions::InvalidPropertyException, "The reward model is not unique and the formula " << formula << " does not specify a reward model."); - // reward finiteness has to be checked later iff infinite reward is possible for the subformula - currentObjective.rewardFinitenessChecked = formula.getSubformula().isCumulativeRewardFormula(); - if(formula.getSubformula().isEventuallyFormula()){ - preprocessFormula(formula.getSubformula().asEventuallyFormula(), result, currentObjective, false, false, formula.getOptionalRewardModelName()); + preprocessFormula(formula.getSubformula().asEventuallyFormula(), result, currentObjective, formula.getOptionalRewardModelName()); } else if(formula.getSubformula().isCumulativeRewardFormula()) { preprocessFormula(formula.getSubformula().asCumulativeRewardFormula(), result, currentObjective, formula.getOptionalRewardModelName()); } else if(formula.getSubformula().isTotalRewardFormula()) { @@ -162,22 +171,19 @@ namespace storm { } template - void SparsePCAAPreprocessor::preprocessFormula(storm::logic::TimeOperatorFormula const& formula, ReturnType& result, PCAAObjective& currentObjective) { + void SparsePcaaPreprocessor::preprocessFormula(storm::logic::TimeOperatorFormula const& formula, ReturnType& result, PcaaObjective& currentObjective) { // Time formulas are only supported for Markov automata STORM_LOG_THROW(result.originalModel.isOfType(storm::models::ModelType::MarkovAutomaton), storm::exceptions::InvalidPropertyException, "Time operator formulas are only supported for Markov automata."); - // reward finiteness does not need to be checked if we want to minimize time - currentObjective.rewardFinitenessChecked = !currentObjective.rewardsArePositive; - if(formula.getSubformula().isEventuallyFormula()){ - preprocessFormula(formula.getSubformula().asEventuallyFormula(), result, currentObjective, false, false); + preprocessFormula(formula.getSubformula().asEventuallyFormula(), result, currentObjective); } else { STORM_LOG_THROW(false, storm::exceptions::InvalidPropertyException, "The subformula of " << formula << " is not supported."); } } template - void SparsePCAAPreprocessor::preprocessFormula(storm::logic::UntilFormula const& formula, ReturnType& result, PCAAObjective& currentObjective) { + void SparsePcaaPreprocessor::preprocessFormula(storm::logic::UntilFormula const& formula, ReturnType& result, PcaaObjective& currentObjective) { CheckTask phiTask(formula.getLeftSubformula()); CheckTask psiTask(formula.getRightSubformula()); storm::modelchecker::SparsePropositionalModelChecker mc(result.preprocessedModel); @@ -215,7 +221,7 @@ namespace storm { } template - void SparsePCAAPreprocessor::preprocessFormula(storm::logic::BoundedUntilFormula const& formula, ReturnType& result, PCAAObjective& currentObjective) { + void SparsePcaaPreprocessor::preprocessFormula(storm::logic::BoundedUntilFormula const& formula, ReturnType& result, PcaaObjective& currentObjective) { if(formula.hasDiscreteTimeBound()) { currentObjective.upperTimeBound = storm::utility::convertNumber(formula.getDiscreteTimeBound()); @@ -232,11 +238,11 @@ namespace storm { } currentObjective.upperTimeBound = storm::utility::convertNumber(formula.getIntervalBounds().second); } - preprocessFormula(storm::logic::UntilFormula(formula.getLeftSubformula().asSharedPointer(), formula.getRightSubformula().asSharedPointer()), result, currentObjective, false, false); + preprocessFormula(storm::logic::UntilFormula(formula.getLeftSubformula().asSharedPointer(), formula.getRightSubformula().asSharedPointer()), result, currentObjective); } template - void SparsePCAAPreprocessor::preprocessFormula(storm::logic::GloballyFormula const& formula, ReturnType& result, PCAAObjective& currentObjective) { + void SparsePcaaPreprocessor::preprocessFormula(storm::logic::GloballyFormula const& formula, ReturnType& result, PcaaObjective& currentObjective) { // The formula will be transformed to an until formula for the complementary event. // If the original formula minimizes, the complementary one will maximize and vice versa. // Hence, the decision whether to consider positive or negative rewards flips. @@ -247,14 +253,13 @@ namespace storm { auto negatedSubformula = std::make_shared(storm::logic::UnaryBooleanStateFormula::OperatorType::Not, formula.getSubformula().asSharedPointer()); - // We need to swap the two flags isProb0Formula and isProb1Formula - preprocessFormula(storm::logic::UntilFormula(storm::logic::Formula::getTrueFormula(), negatedSubformula), result, currentObjective, isProb1Formula, isProb0Formula); + preprocessFormula(storm::logic::UntilFormula(storm::logic::Formula::getTrueFormula(), negatedSubformula), result, currentObjective); } template - void SparsePCAAPreprocessor::preprocessFormula(storm::logic::EventuallyFormula const& formula, ReturnType& result, PCAAObjective& currentObjective, boost::optional const& optionalRewardModelName) { + void SparsePcaaPreprocessor::preprocessFormula(storm::logic::EventuallyFormula const& formula, ReturnType& result, PcaaObjective& currentObjective, boost::optional const& optionalRewardModelName) { if(formula.isReachabilityProbabilityFormula()){ - preprocessFormula(storm::logic::UntilFormula(storm::logic::Formula::getTrueFormula(), formula.getSubformula().asSharedPointer()), result, currentObjective, isProb0Formula, isProb1Formula); + preprocessFormula(storm::logic::UntilFormula(storm::logic::Formula::getTrueFormula(), formula.getSubformula().asSharedPointer()), result, currentObjective); return; } @@ -297,7 +302,7 @@ namespace storm { } template - void SparsePCAAPreprocessor::preprocessFormula(storm::logic::CumulativeRewardFormula const& formula, ReturnType& result, PCAAObjective& currentObjective, boost::optional const& optionalRewardModelName) { + void SparsePcaaPreprocessor::preprocessFormula(storm::logic::CumulativeRewardFormula const& formula, ReturnType& result, PcaaObjective& currentObjective, boost::optional const& optionalRewardModelName) { STORM_LOG_THROW(result.originalModel.isOfType(storm::models::ModelType::Mdp), storm::exceptions::InvalidPropertyException, "Cumulative reward formulas are not supported for the given model type."); STORM_LOG_THROW(formula.hasDiscreteTimeBound(), storm::exceptions::InvalidPropertyException, "Expected a cumulativeRewardFormula with a discrete time bound but got " << formula << "."); STORM_LOG_THROW(formula.getDiscreteTimeBound()>0, storm::exceptions::InvalidPropertyException, "Expected a cumulativeRewardFormula with a positive discrete time bound but got " << formula << "."); @@ -317,7 +322,7 @@ namespace storm { } template - void SparsePCAAPreprocessor::preprocessFormula(storm::logic::TotalRewardFormula const& formula, ReturnType& result, PCAAObjective& currentObjective, boost::optional const& optionalRewardModelName) { + void SparsePcaaPreprocessor::preprocessFormula(storm::logic::TotalRewardFormula const& formula, ReturnType& result, PcaaObjective& currentObjective, boost::optional const& optionalRewardModelName) { RewardModelType objectiveRewards = result.preprocessedModel.getRewardModel(optionalRewardModelName ? optionalRewardModelName.get() : ""); objectiveRewards.reduceToStateBasedRewards(result.preprocessedModel.getTransitionMatrix(), false); if(!currentObjective.rewardsArePositive){ @@ -333,87 +338,76 @@ namespace storm { template - void SparsePCAAPreprocessor::ensureRewardFiniteness(ReturnType& result) { - storm::storage::BitVector actionsWithNegativeReward(result.preprocessedModel.getTransitionMatrix().getRowCount(), false); - storm::storage::BitVector actionsWithPositiveReward(result.preprocessedModel.getTransitionMatrix().getRowCount(), false); - for(uint_fast64_t objIndex = 0; objIndex < result.objectives.size(); ++objIndex) { - if (!result.objectives[objIndex].rewardFinitenessChecked) { - result.objectives[objIndex].rewardFinitenessChecked = true; - if(result.objectives[objIndex].rewardsArePositive) { - actionsWithPositiveReward |= storm::utility::vector::filter(result.preprocessedModel.getRewardModel(result.objectives[objIndex].rewardModelName).getTotalRewardVector(result.preprocessedModel.getTransitionMatrix()), [&] (ValueType const& value) -> bool { return !storm::utility::isZero(value);}); - } else { - actionsWithNegativeReward |= storm::utility::vector::filter(result.preprocessedModel.getRewardModel(result.objectives[objIndex].rewardModelName).getTotalRewardVector(result.preprocessedModel.getTransitionMatrix()), [&] (ValueType const& value) -> bool { return !storm::utility::isZero(value);}); - } - } - } - if(actionsWithPositiveReward.empty() && actionsWithNegativeReward.empty()) { - // No rewards for which we need to ensure finiteness - result.possibleEcActions = storm::storage::BitVector(result.preprocessedModel.getNumberOfChoices(), true); - result.statesWhichCanBeVisitedInfinitelyOften = storm::storage::BitVector(result.preprocessedModel.getNumberOfStates(), true); - return; - } + void SparsePcaaPreprocessor::analyzeEndComponents(ReturnType& result, storm::storage::SparseMatrix const& backwardTransitions) { - result.possibleEcActions = storm::storage::BitVector(result.preprocessedModel.getNumberOfChoices(), false); - result.statesWhichCanBeVisitedInfinitelyOften = storm::storage::BitVector(result.preprocessedModel.getNumberOfStates(), false); - auto backwardTransitions = result.preprocessedModel.getBackwardTransitions(); + result.ecActions = storm::storage::BitVector(result.preprocessedModel.getNumberOfChoices(), false); + std::vector ecs; auto mecDecomposition = storm::storage::MaximalEndComponentDecomposition(result.preprocessedModel.getTransitionMatrix(), backwardTransitions); STORM_LOG_ASSERT(!mecDecomposition.empty(), "Empty maximal end component decomposition."); - std::vector ecs; ecs.reserve(mecDecomposition.size()); for(auto& mec : mecDecomposition) { for(auto const& stateActionsPair : mec) { for(auto const& action : stateActionsPair.second) { - result.possibleEcActions.set(action); - STORM_LOG_THROW(!actionsWithPositiveReward.get(action), storm::exceptions::InvalidPropertyException, "Infinite reward: Found an end componet that induces infinite reward for at least one objective"); + result.ecActions.set(action); } } ecs.push_back(std::move(mec)); } - storm::storage::BitVector currentECStates(result.preprocessedModel.getNumberOfStates(), false); + result.possiblyRecurrentStates = storm::storage::BitVector(result.preprocessedModel.getNumberOfStates(), false); + storm::storage::BitVector neutralActions = ~(result.actionsWithNegativeReward | result.actionsWithPositiveReward); + storm::storage::BitVector statesInCurrentECWithNeutralAction(result.preprocessedModel.getNumberOfStates(), false); for(uint_fast64_t ecIndex = 0; ecIndex < ecs.size(); ++ecIndex) { //we will insert new ecs in the vector (thus no iterators for the loop) bool currentECIsNeutral = true; for(auto const& stateActionsPair : ecs[ecIndex]) { bool stateHasNeutralActionWithinEC = false; for(auto const& action : stateActionsPair.second) { - stateHasNeutralActionWithinEC |= !actionsWithNegativeReward.get(action); + stateHasNeutralActionWithinEC |= neutralActions.get(action); } - currentECStates.set(stateActionsPair.first, stateHasNeutralActionWithinEC); + statesInCurrentECWithNeutralAction.set(stateActionsPair.first, stateHasNeutralActionWithinEC); currentECIsNeutral &= stateHasNeutralActionWithinEC; } if(currentECIsNeutral) { - result.statesWhichCanBeVisitedInfinitelyOften |= currentECStates; + result.possiblyRecurrentStates |= statesInCurrentECWithNeutralAction; }else{ // Check if the ec contains neutral sub ecs. This is done by adding the subECs to our list of ECs - auto subECs = storm::storage::MaximalEndComponentDecomposition(result.preprocessedModel.getTransitionMatrix(), backwardTransitions, currentECStates); + auto subECs = storm::storage::MaximalEndComponentDecomposition(result.preprocessedModel.getTransitionMatrix(), backwardTransitions, statesInCurrentECWithNeutralAction); ecs.reserve(ecs.size() + subECs.size()); for(auto& ec : subECs){ ecs.push_back(std::move(ec)); } } - currentECStates.clear(); + statesInCurrentECWithNeutralAction.clear(); } + } + + template + void SparsePcaaPreprocessor::ensureRewardFiniteness(ReturnType& result, storm::storage::SparseMatrix const& backwardTransitions) { + STORM_LOG_THROW((result.actionsWithPositiveReward & result.ecActions).empty(), storm::exceptions::InvalidPropertyException, "Infinite reward: There is one maximizing objective for which infinite reward is possible. This is not supported."); //Check whether the states that can be visited inf. often are reachable with prob. 1 under some scheduler - storm::storage::BitVector statesReachingNegativeRewardsFinitelyOftenForSomeScheduler = storm::utility::graph::performProb1E(result.preprocessedModel.getTransitionMatrix(), result.preprocessedModel.getTransitionMatrix().getRowGroupIndices(), backwardTransitions, storm::storage::BitVector(result.preprocessedModel.getNumberOfStates()), result.statesWhichCanBeVisitedInfinitelyOften); - STORM_LOG_Throw(!(statesReachingNegativeRewardsFinitelyOftenForSomeScheduler & result.preprocessedModel.getInitialStates()).empty(), storm::exceptions::InvalidPropertyException, "Infinite Rewards: For every scheduler, the induced reward for one or more of the objectives that minimize rewards is infinity."); + storm::storage::BitVector statesReachingNegativeRewardsFinitelyOftenForSomeScheduler = storm::utility::graph::performProb1E(result.preprocessedModel.getTransitionMatrix(), result.preprocessedModel.getTransitionMatrix().getRowGroupIndices(), backwardTransitions, storm::storage::BitVector(result.preprocessedModel.getNumberOfStates(), true), result.possiblyRecurrentStates); + STORM_LOG_THROW(!(statesReachingNegativeRewardsFinitelyOftenForSomeScheduler & result.preprocessedModel.getInitialStates()).empty(), storm::exceptions::InvalidPropertyException, "Infinite Rewards: For every scheduler, the induced reward for one or more of the objectives that minimize rewards is infinity."); if(!statesReachingNegativeRewardsFinitelyOftenForSomeScheduler.full()) { - auto subsystemBuilderResult = storm::transformer::SubsystemBuilder::transform(result.preprocessedModel, statesReachingNegativeRewardsFinitelyOftenForSomeScheduler, storm::storage::BitVector(result.preprocessedModel.getTransitionMatrix().getRowCount(), true)); + // Remove the states that for any scheduler have one objective with infinite expected reward. + auto subsystemBuilderResult = storm::transformer::SubsystemBuilder::transform(result.preprocessedModel, statesReachingNegativeRewardsFinitelyOftenForSomeScheduler, storm::storage::BitVector(result.preprocessedModel.getNumberOfChoices(), true)); updatePreprocessedModel(result, *subsystemBuilderResult.model, subsystemBuilderResult.newToOldStateIndexMapping); - result.possibleEcActions = result.possibleEcActions % subsystemBuilderResult.subsystemActions; - result.statesWhichCanBeVisitedInfinitelyOften = result.statesWhichCanBeVisitedInfinitelyOften % statesReachingNegativeRewardsFinitelyOftenForSomeScheduler; + result.ecActions = result.ecActions % subsystemBuilderResult.subsystemActions; + result.actionsWithPositiveReward = result.actionsWithPositiveReward % subsystemBuilderResult.subsystemActions; + result.actionsWithNegativeReward = result.actionsWithNegativeReward % subsystemBuilderResult.subsystemActions; + result.possiblyRecurrentStates = result.possiblyRecurrentStates % statesReachingNegativeRewardsFinitelyOftenForSomeScheduler; } } - template class SparsePCAAPreprocessor>; - template class SparsePCAAPreprocessor>; + template class SparsePcaaPreprocessor>; + template class SparsePcaaPreprocessor>; #ifdef STORM_HAVE_CARL - template class SparsePCAAPreprocessor>; - template class SparsePCAAPreprocessor>; + template class SparsePcaaPreprocessor>; + template class SparsePcaaPreprocessor>; #endif } } diff --git a/src/modelchecker/multiobjective/pcaa/SparsePCAAPreprocessor.h b/src/modelchecker/multiobjective/pcaa/SparsePCAAPreprocessor.h index 336c129d0..aa8832cc5 100644 --- a/src/modelchecker/multiobjective/pcaa/SparsePCAAPreprocessor.h +++ b/src/modelchecker/multiobjective/pcaa/SparsePCAAPreprocessor.h @@ -5,7 +5,7 @@ #include "src/logic/Formulas.h" #include "src/storage/BitVector.h" -#include "src/modelchecker/multiobjective/pcaa/SparsePCAAPreprocessorReturnType.h" +#include "src/modelchecker/multiobjective/pcaa/SparsePcaaPreprocessorReturnType.h" namespace storm { namespace modelchecker { @@ -15,18 +15,18 @@ namespace storm { * This class invokes the necessary preprocessing for the Pareto Curve Approximation Algorithm (PCAA) */ template - class SparsePCAAPreprocessor { + class SparsePcaaPreprocessor { public: typedef typename SparseModelType::ValueType ValueType; typedef typename SparseModelType::RewardModelType RewardModelType; - typedef SparsePCAAPreprocessorReturnType ReturnType; + typedef SparsePcaaPreprocessorReturnType ReturnType; /*! * Preprocesses the given model w.r.t. the given formulas. * @param originalModel The considered model * @param originalFormula the considered formula. The subformulas should only contain one OperatorFormula at top level, i.e., the formula is simple. */ - static ReturnType preprocess(storm::logic::MultiObjectiveFormula const& originalFormula, SparseModelType const& originalModel); + static ReturnType preprocess(SparseModelType const& originalModel, storm::logic::MultiObjectiveFormula const& originalFormula); private: /*! @@ -50,23 +50,29 @@ namespace storm { * @param currentObjective the currently considered objective. The given formula should be a a (sub)formula of this objective * @param optionalRewardModelName the reward model name that is considered for the formula (if available) */ - static void preprocessFormula(storm::logic::OperatorFormula const& formula, ReturnType& result, PCAAObjective& currentObjective); - static void preprocessFormula(storm::logic::ProbabilityOperatorFormula const& formula, ReturnType& result, PCAAObjective& currentObjective); - static void preprocessFormula(storm::logic::RewardOperatorFormula const& formula, ReturnType& result, PCAAObjective& currentObjective); - static void preprocessFormula(storm::logic::TimeOperatorFormula const& formula, ReturnType& result, PCAAObjective& currentObjective); - static void preprocessFormula(storm::logic::UntilFormula const& formula, ReturnType& result, PCAAObjective& currentObjective); - static void preprocessFormula(storm::logic::BoundedUntilFormula const& formula, ReturnType& result, PCAAObjective& currentObjective); - static void preprocessFormula(storm::logic::GloballyFormula const& formula, ReturnType& result, PCAAObjective& currentObjective); - static void preprocessFormula(storm::logic::EventuallyFormula const& formula, ReturnType& result, PCAAObjective& currentObjective, boost::optional const& optionalRewardModelName = boost::none); - static void preprocessFormula(storm::logic::CumulativeRewardFormula const& formula, ReturnType& result, PCAAObjective& currentObjective, boost::optional const& optionalRewardModelName = boost::none); - static void preprocessFormula(storm::logic::TotalRewardFormula const& formula, ReturnType& result, PCAAObjective& currentObjective, boost::optional const& optionalRewardModelName = boost::none); + static void preprocessFormula(storm::logic::OperatorFormula const& formula, ReturnType& result, PcaaObjective& currentObjective); + static void preprocessFormula(storm::logic::ProbabilityOperatorFormula const& formula, ReturnType& result, PcaaObjective& currentObjective); + static void preprocessFormula(storm::logic::RewardOperatorFormula const& formula, ReturnType& result, PcaaObjective& currentObjective); + static void preprocessFormula(storm::logic::TimeOperatorFormula const& formula, ReturnType& result, PcaaObjective& currentObjective); + static void preprocessFormula(storm::logic::UntilFormula const& formula, ReturnType& result, PcaaObjective& currentObjective); + static void preprocessFormula(storm::logic::BoundedUntilFormula const& formula, ReturnType& result, PcaaObjective& currentObjective); + static void preprocessFormula(storm::logic::GloballyFormula const& formula, ReturnType& result, PcaaObjective& currentObjective); + static void preprocessFormula(storm::logic::EventuallyFormula const& formula, ReturnType& result, PcaaObjective& currentObjective, boost::optional const& optionalRewardModelName = boost::none); + static void preprocessFormula(storm::logic::CumulativeRewardFormula const& formula, ReturnType& result, PcaaObjective& currentObjective, boost::optional const& optionalRewardModelName = boost::none); + static void preprocessFormula(storm::logic::TotalRewardFormula const& formula, ReturnType& result, PcaaObjective& currentObjective, boost::optional const& optionalRewardModelName = boost::none); /*! - * Checks whether the occurring reward properties are guaranteed to be finite for all states. - * if not, the input is rejected. - * Also applies further preprocessing steps regarding End Component Elimination + * Analyzes the end components of the preprocessed model. That is: + * -get the set of actions that are part of an end component + * -Find the states that can be visited infinitely often without inducing infinite reward */ - static void ensureRewardFiniteness(ReturnType& result); + static void analyzeEndComponents(ReturnType& result, storm::storage::SparseMatrix const& backwardTransitions); + + /*! + * Checks whether the occurring expected rewards are finite. If not, the input is rejected. + * Also removes all states for which no finite reward wrt. all objectives is possible + */ + static void ensureRewardFiniteness(ReturnType& result, storm::storage::SparseMatrix const& backwardTransitions); }; } diff --git a/src/modelchecker/multiobjective/pcaa/SparsePCAAPreprocessorReturnType.h b/src/modelchecker/multiobjective/pcaa/SparsePCAAPreprocessorReturnType.h index 02f3ed780..271c85ebc 100644 --- a/src/modelchecker/multiobjective/pcaa/SparsePCAAPreprocessorReturnType.h +++ b/src/modelchecker/multiobjective/pcaa/SparsePCAAPreprocessorReturnType.h @@ -8,7 +8,7 @@ #include #include "src/logic/Formulas.h" -#include "src/modelchecker/multiobjective/pcaa/PCAAObjective.h" +#include "src/modelchecker/multiobjective/pcaa/PcaaObjective.h" #include "src/models/sparse/MarkovAutomaton.h" #include "src/storage/BitVector.h" #include "src/utility/macros.h" @@ -21,7 +21,7 @@ namespace storm { namespace multiobjective { template - struct SparsePCAAPreprocessorReturnType { + struct SparsePcaaPreprocessorReturnType { enum class QueryType { Achievability, Quantitative, Pareto }; @@ -29,25 +29,28 @@ namespace storm { SparseModelType const& originalModel; SparseModelType preprocessedModel; QueryType queryType; - - // Maps any state of the preprocessed model to the corresponding state of the original Model - std::vector newToOldStateIndexMapping; - - // The actions of the preprocessed model that can be part of an EC - storm::storage::BitVector possibleEcActions; - - // The set of states of the preprocessed model for which there is a scheduler such that - // the state is visited infinitely often and the induced reward is finite for any objective - storm::storage::BitVector statesWhichCanBeVisitedInfinitelyOften; // The (preprocessed) objectives - std::vector> objectives; + std::vector> objectives; // The index of the objective that is to be maximized (or minimized) in case of a quantitative Query boost::optional indexOfOptimizingObjective; + + // Maps any state of the preprocessed model to the corresponding state of the original Model + std::vector newToOldStateIndexMapping; + // The actions of the preprocessed model that have positive reward assigned for at least one objective + storm::storage::BitVector actionsWithPositiveReward; + // The actions of the preprocessed model that have negative reward assigned for at least one objective + storm::storage::BitVector actionsWithNegativeReward; + // The actions of the preprocessed model that are part of an EC + storm::storage::BitVector ecActions; + + // The set of states of the preprocessed model for which there is a scheduler such that + // the state is visited infinitely often and the induced reward is finite for all objectives + storm::storage::BitVector possiblyRecurrentStates; - SparsePCAAPreprocessorReturnType(storm::logic::MultiObjectiveFormula const& originalFormula, SparseModelType const& originalModel, SparseModelType&& preprocessedModel) : originalFormula(originalFormula), originalModel(originalModel), preprocessedModel(preprocessedModel) { + SparsePcaaPreprocessorReturnType(storm::logic::MultiObjectiveFormula const& originalFormula, SparseModelType const& originalModel, SparseModelType&& preprocessedModel) : originalFormula(originalFormula), originalModel(originalModel), preprocessedModel(preprocessedModel) { // Intentionally left empty } @@ -77,7 +80,7 @@ namespace storm { out << "---------------------------------------------------------------------------------------------------------------------------------------" << std::endl; } - friend std::ostream& operator<<(std::ostream& out, SparsePCAAPreprocessorReturnType const& ret) { + friend std::ostream& operator<<(std::ostream& out, SparsePcaaPreprocessorReturnType const& ret) { ret.printToStream(out); return out; } diff --git a/src/modelchecker/multiobjective/pcaa/SparsePcaaQuantitativeQuery.cpp b/src/modelchecker/multiobjective/pcaa/SparsePcaaQuantitativeQuery.cpp new file mode 100644 index 000000000..0d4d879f3 --- /dev/null +++ b/src/modelchecker/multiobjective/pcaa/SparsePcaaQuantitativeQuery.cpp @@ -0,0 +1,169 @@ +#include "src/modelchecker/multiobjective/pcaa/SparsePcaaQuantitativeQuery.h" + +#include "src/adapters/CarlAdapter.h" +#include "src/models/sparse/Mdp.h" +#include "src/models/sparse/MarkovAutomaton.h" +#include "src/models/sparse/StandardRewardModel.h" +#include "src/modelchecker/results/ExplicitQualitativeCheckResult.h" +#include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h" +#include "src/utility/constants.h" +#include "src/utility/vector.h" +#include "src/settings//SettingsManager.h" +#include "src/settings/modules/MultiObjectiveSettings.h" +#include "src/settings/modules/GeneralSettings.h" + +namespace storm { + namespace modelchecker { + namespace multiobjective { + + + template + SparsePcaaQuantitativeQuery::SparsePcaaQuantitativeQuery(SparsePcaaPreprocessorReturnType& preprocessorResult) : SparsePcaaQuery(preprocessorResult) { + STORM_LOG_ASSERT(preprocessorResult.queryType==SparsePcaaPreprocessorReturnType::QueryType::Quantitative, "Invalid query Type"); + STORM_LOG_ASSERT(preprocessorResult.indexOfOptimizingObjective, "Detected quantitative query but index of optimizing objective is not set."); + indexOfOptimizingObjective = *preprocessorResult.indexOfOptimizingObjective; + initializeThresholdData(); + + // Set the maximum gap between lower and upper bound of the weightVectorChecker result. + // This is the maximal edge length of the box we have to consider around each computed point + // We pick the gap such that the maximal distance between two points within this box is less than the given precision divided by two. + typename SparseModelType::ValueType gap = storm::utility::convertNumber(storm::settings::getModule().getPrecision()); + gap /= (storm::utility::one() + storm::utility::one()); + gap /= storm::utility::sqrt(static_cast(this->objectives.size())); + this->weightVectorChecker->setMaximumLowerUpperBoundGap(gap); + + } + + template + void SparsePcaaQuantitativeQuery::initializeThresholdData() { + thresholds.reserve(this->objectives.size()); + strictThresholds = storm::storage::BitVector(this->objectives.size(), false); + std::vector> thresholdConstraints; + thresholdConstraints.reserve(this->objectives.size()-1); + for(uint_fast64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) { + if(this->objectives[objIndex].threshold) { + thresholds.push_back(storm::utility::convertNumber(*this->objectives[objIndex].threshold)); + WeightVector normalVector(this->objectives.size(), storm::utility::zero()); + normalVector[objIndex] = -storm::utility::one(); + thresholdConstraints.emplace_back(std::move(normalVector), -thresholds.back()); + strictThresholds.set(objIndex, this->objectives[objIndex].thresholdIsStrict); + } else { + thresholds.push_back(storm::utility::zero()); + } + } + // Note: If we have a single objective (i.e., no objectives with thresholds), thresholdsAsPolytope gets no constraints + thresholdsAsPolytope = storm::storage::geometry::Polytope::create(thresholdConstraints); + } + + template + std::unique_ptr SparsePcaaQuantitativeQuery::check() { + + + // First find one solution that achieves the given thresholds ... + if(this->checkAchievability()) { + // ... then improve it + GeometryValueType result = this->improveSolution(); + + // transform the obtained result for the preprocessed model to a result w.r.t. the original model and return the checkresult + typename SparseModelType::ValueType resultForOriginalModel = + storm::utility::convertNumber(result) * + this->objectives[indexOfOptimizingObjective].toOriginalValueTransformationFactor + + this->objectives[indexOfOptimizingObjective].toOriginalValueTransformationOffset; + return std::unique_ptr(new ExplicitQuantitativeCheckResult(this->originalModel.getInitialStates().getNextSetIndex(0), resultForOriginalModel)); + } else { + return std::unique_ptr(new ExplicitQualitativeCheckResult(this->originalModel.getInitialStates().getNextSetIndex(0), false)); + } + + + } + + template + bool SparsePcaaQuantitativeQuery::checkAchievability() { + // We don't care for the optimizing objective at this point + this->diracWeightVectorsToBeChecked.set(indexOfOptimizingObjective, false); + + while(!this->maxStepsPerformed()){ + WeightVector separatingVector = this->findSeparatingVector(thresholds); + this->performRefinementStep(std::move(separatingVector)); + //Pick the threshold for the optimizing objective low enough so valid solutions are not excluded + thresholds[indexOfOptimizingObjective] = std::min(thresholds[indexOfOptimizingObjective], this->refinementSteps.back().lowerBoundPoint[indexOfOptimizingObjective]); + if(!checkIfThresholdsAreSatisfied(this->overApproximation)){ + return false; + } + if(checkIfThresholdsAreSatisfied(this->underApproximation)){ + return true; + } + } + STORM_LOG_ERROR("Could not check whether thresholds are achievable: Exceeded maximum number of refinement steps"); + return false; + } + + template + GeometryValueType SparsePcaaQuantitativeQuery::improveSolution() { + this->diracWeightVectorsToBeChecked.clear(); // Only check weight vectors that can actually improve the solution + + WeightVector directionOfOptimizingObjective(this->objectives.size(), storm::utility::zero()); + directionOfOptimizingObjective[indexOfOptimizingObjective] = storm::utility::one(); + + // Improve the found solution. + // Note that we do not care whether a threshold is strict anymore, because the resulting optimum should be + // the supremum over all strategies. Hence, one could combine a scheduler inducing the optimum value (but possibly violating strict + // thresholds) and (with very low probability) a scheduler that satisfies all (possibly strict) thresholds. + GeometryValueType result = storm::utility::zero(); + while(!this->maxStepsPerformed()) { + std::pair optimizationRes = this->underApproximation->intersection(thresholdsAsPolytope)->optimize(directionOfOptimizingObjective); + STORM_LOG_THROW(optimizationRes.second, storm::exceptions::UnexpectedException, "The underapproximation is either unbounded or empty."); + result = optimizationRes.first[indexOfOptimizingObjective]; + STORM_LOG_DEBUG("Best solution found so far is ~" << storm::utility::convertNumber(result) << "."); + thresholds[indexOfOptimizingObjective] = result; + //Compute an upper bound for the optimum and check for convergence + optimizationRes = this->overApproximation->intersection(thresholdsAsPolytope)->optimize(directionOfOptimizingObjective); + if(optimizationRes.second) { + GeometryValueType precisionOfResult = optimizationRes.first[indexOfOptimizingObjective] - result; + if(precisionOfResult < storm::utility::convertNumber(storm::settings::getModule().getPrecision())) { + // Goal precision reached! + return result; + } else { + STORM_LOG_DEBUG("Solution can be improved by at most " << storm::utility::convertNumber(precisionOfResult)); + } + } + WeightVector separatingVector = this->findSeparatingVector(thresholds); + this->performRefinementStep(std::move(separatingVector)); + } + STORM_LOG_ERROR("Could not reach the desired precision: Exceeded maximum number of refinement steps"); + return result; + } + + template + bool SparsePcaaQuantitativeQuery::checkIfThresholdsAreSatisfied(std::shared_ptr> const& polytope) { + std::vector> halfspaces = polytope->getHalfspaces(); + for(auto const& h : halfspaces) { + GeometryValueType distance = h.distance(thresholds); + if(distance < storm::utility::zero()) { + return false; + } + if(distance == storm::utility::zero()) { + // In this case, the thresholds point is on the boundary of the polytope. + // Check if this is problematic for the strict thresholds + for(auto strictThreshold : strictThresholds) { + if(h.normalVector()[strictThreshold] > storm::utility::zero()) { + return false; + } + } + } + } + return true; + } + + + +#ifdef STORM_HAVE_CARL + template class SparsePcaaQuantitativeQuery, storm::RationalNumber>; + template class SparsePcaaQuantitativeQuery, storm::RationalNumber>; + + template class SparsePcaaQuantitativeQuery, storm::RationalNumber>; + // template class SparsePcaaQuantitativeQuery, storm::RationalNumber>; +#endif + } + } +} diff --git a/src/modelchecker/multiobjective/pcaa/SparsePcaaQuantitativeQuery.h b/src/modelchecker/multiobjective/pcaa/SparsePcaaQuantitativeQuery.h new file mode 100644 index 000000000..b2313e529 --- /dev/null +++ b/src/modelchecker/multiobjective/pcaa/SparsePcaaQuantitativeQuery.h @@ -0,0 +1,65 @@ +#ifndef STORM_MODELCHECKER_MULTIOBJECTIVE_PCAA_SPARSEPCAAQUANTITATIVEQUERY_H_ +#define STORM_MODELCHECKER_MULTIOBJECTIVE_PCAA_SPARSEPCAAQUANTITATIVEQUERY_H_ + +#include "src/modelchecker/multiobjective/pcaa/SparsePcaaQuery.h" + +namespace storm { + namespace modelchecker { + namespace multiobjective { + + /* + * This class represents a query for the Pareto curve approximation algorithm (Pcaa). + * It implements the necessary computations for the different query types. + */ + template + class SparsePcaaQuantitativeQuery : public SparsePcaaQuery { + public: + + // Typedefs for simple geometric objects + typedef std::vector Point; + typedef std::vector WeightVector; + + /* + * Creates a new query for the Pareto curve approximation algorithm (Pcaa) + * @param preprocessorResult the result from preprocessing + */ + SparsePcaaQuantitativeQuery(SparsePcaaPreprocessorReturnType& preprocessorResult); + + + /* + * Invokes the computation and retrieves the result + */ + virtual std::unique_ptr check() override; + + private: + + void initializeThresholdData(); + + /* + * Returns whether the given thresholds are achievable. + */ + bool checkAchievability(); + + /* + * Given that the thresholds are achievable, this function further refines the approximations and returns the optimized value + */ + GeometryValueType improveSolution(); + + /* + * Returns true iff there is one point in the given polytope that satisfies the given thresholds. + * It is assumed that the given polytope contains the downward closure of its vertices. + */ + bool checkIfThresholdsAreSatisfied(std::shared_ptr> const& polytope); + + uint_fast64_t indexOfOptimizingObjective; + + Point thresholds; + storm::storage::BitVector strictThresholds; + std::shared_ptr> thresholdsAsPolytope; + }; + + } + } +} + +#endif /* STORM_MODELCHECKER_MULTIOBJECTIVE_PCAA_SPARSEPCAAQUANTITATIVEQUERY_H_ */ diff --git a/src/modelchecker/multiobjective/pcaa/SparsePcaaQuery.cpp b/src/modelchecker/multiobjective/pcaa/SparsePcaaQuery.cpp new file mode 100644 index 000000000..45390e711 --- /dev/null +++ b/src/modelchecker/multiobjective/pcaa/SparsePcaaQuery.cpp @@ -0,0 +1,205 @@ +#include "src/modelchecker/multiobjective/pcaa/SparsePcaaQuery.h" + +#include "src/adapters/CarlAdapter.h" +#include "src/models/sparse/Mdp.h" +#include "src/models/sparse/MarkovAutomaton.h" +#include "src/models/sparse/StandardRewardModel.h" +#include "src/modelchecker/multiobjective/pcaa/PcaaObjective.h" +#include "src/modelchecker/multiobjective/pcaa/SparseMdpPcaaWeightVectorChecker.h" +#include "src/modelchecker/multiobjective/pcaa/SparseMaPcaaWeightVectorChecker.h" +#include "src/utility/constants.h" +#include "src/utility/vector.h" +#include "src/settings//SettingsManager.h" +#include "src/settings/modules/MultiObjectiveSettings.h" + +#include "src/exceptions/UnexpectedException.h" + +namespace storm { + namespace modelchecker { + namespace multiobjective { + + + template + SparsePcaaQuery::SparsePcaaQuery(SparsePcaaPreprocessorReturnType& preprocessorResult) : + originalModel(preprocessorResult.originalModel), originalFormula(preprocessorResult.originalFormula), + preprocessedModel(std::move(preprocessorResult.preprocessedModel)), objectives(std::move(preprocessorResult.objectives)) { + initializeWeightVectorChecker(preprocessedModel, objectives, preprocessorResult.actionsWithNegativeReward, preprocessorResult.ecActions, preprocessorResult.possiblyRecurrentStates); + this->diracWeightVectorsToBeChecked = storm::storage::BitVector(this->objectives.size(), true); + this->overApproximation = storm::storage::geometry::Polytope::createUniversalPolytope(); + this->underApproximation = storm::storage::geometry::Polytope::createEmptyPolytope(); + } + + template<> + void SparsePcaaQuery, storm::RationalNumber>::initializeWeightVectorChecker(storm::models::sparse::Mdp const& model, std::vector> const& objectives, storm::storage::BitVector const& actionsWithNegativeReward, storm::storage::BitVector const& ecActions, storm::storage::BitVector const& possiblyRecurrentStates) { + this->weightVectorChecker = std::unique_ptr>>(new SparseMdpPcaaWeightVectorChecker>(model, objectives, actionsWithNegativeReward, ecActions, possiblyRecurrentStates)); + } + + template<> + void SparsePcaaQuery, storm::RationalNumber>::initializeWeightVectorChecker(storm::models::sparse::Mdp const& model, std::vector> const& objectives, storm::storage::BitVector const& actionsWithNegativeReward, storm::storage::BitVector const& ecActions, storm::storage::BitVector const& possiblyRecurrentStates) { + this->weightVectorChecker = std::unique_ptr>>(new SparseMdpPcaaWeightVectorChecker>(model, objectives, actionsWithNegativeReward, ecActions, possiblyRecurrentStates)); + } + + template<> + void SparsePcaaQuery, storm::RationalNumber>::initializeWeightVectorChecker(storm::models::sparse::MarkovAutomaton const& model, std::vector> const& objectives, storm::storage::BitVector const& actionsWithNegativeReward, storm::storage::BitVector const& ecActions, storm::storage::BitVector const& possiblyRecurrentStates) { + this->weightVectorChecker = std::unique_ptr>>(new SparseMaPcaaWeightVectorChecker>(model, objectives, actionsWithNegativeReward, ecActions, possiblyRecurrentStates)); + } + + // template<> + // void SparsePcaaQuery, storm::RationalNumber>::initializeWeightVectorChecker(storm::models::sparse::MarkovAutomaton const& model, std::vector> const& objectives, storm::storage::BitVector const& actionsWithNegativeReward, storm::storage::BitVector const& ecActions, storm::storage::BitVector const& possiblyRecurrentStates) { + // this->weightVectorChecker = std::unique_ptr>>(new SparseMaPcaaWeightVectorChecker>(model, objectives, actionsWithNegativeReward, ecActions, possiblyRecurrentStates)); + // } + + + template + typename SparsePcaaQuery::WeightVector SparsePcaaQuery::findSeparatingVector(Point const& pointToBeSeparated) { + STORM_LOG_DEBUG("Searching a weight vector to seperate the point given by " << storm::utility::vector::toString(storm::utility::vector::convertNumericVector(pointToBeSeparated)) << "."); + + if(underApproximation->isEmpty()) { + // In this case, every weight vector is separating + uint_fast64_t objIndex = diracWeightVectorsToBeChecked.getNextSetIndex(0) % pointToBeSeparated.size(); + WeightVector result(pointToBeSeparated.size(), storm::utility::zero()); + result[objIndex] = storm::utility::one(); + diracWeightVectorsToBeChecked.set(objIndex, false); + return result; + } + + // Reaching this point means that the underApproximation contains halfspaces. The seperating vector has to be the normal vector of one of these halfspaces. + // We pick one with maximal distance to the given point. However, Dirac weight vectors that only assign a non-zero weight to a single objective take precedence. + std::vector> halfspaces = underApproximation->getHalfspaces(); + uint_fast64_t farestHalfspaceIndex = halfspaces.size(); + GeometryValueType farestDistance = -storm::utility::one(); + bool foundSeparatingDiracVector = false; + for(uint_fast64_t halfspaceIndex = 0; halfspaceIndex < halfspaces.size(); ++halfspaceIndex) { + GeometryValueType distance = -halfspaces[halfspaceIndex].euclideanDistance(pointToBeSeparated); + if(distance >= storm::utility::zero()) { + storm::storage::BitVector nonZeroVectorEntries = ~storm::utility::vector::filterZero(halfspaces[halfspaceIndex].normalVector()); + bool isSingleObjectiveVector = nonZeroVectorEntries.getNumberOfSetBits() == 1 && diracWeightVectorsToBeChecked.get(nonZeroVectorEntries.getNextSetIndex(0)); + // Check if this halfspace is a better candidate than the current one + if((!foundSeparatingDiracVector && isSingleObjectiveVector ) || (foundSeparatingDiracVector==isSingleObjectiveVector && distance>farestDistance)) { + foundSeparatingDiracVector = foundSeparatingDiracVector || isSingleObjectiveVector; + farestHalfspaceIndex = halfspaceIndex; + farestDistance = distance; + } + } + } + if(foundSeparatingDiracVector) { + diracWeightVectorsToBeChecked &= storm::utility::vector::filterZero(halfspaces[farestHalfspaceIndex].normalVector()); + } + + STORM_LOG_THROW(farestHalfspaceIndex(halfspaces[farestHalfspaceIndex].normalVector())) << "."); + return halfspaces[farestHalfspaceIndex].normalVector(); + } + + template + void SparsePcaaQuery::performRefinementStep(WeightVector&& direction) { + // Normalize the direction vector so that the entries sum up to one + storm::utility::vector::scaleVectorInPlace(direction, storm::utility::one() / std::accumulate(direction.begin(), direction.end(), storm::utility::zero())); + // Check if we already did a refinement step with that direction vector. If this is the case, we increase the precision. + // We start with the most recent steps to consider the most recent result for this direction vector + boost::optional oldMaximumLowerUpperBoundGap; + for(auto stepIt = refinementSteps.rbegin(); stepIt != refinementSteps.rend(); ++stepIt) { + if(stepIt->weightVector == direction) { + STORM_LOG_WARN("Performing multiple refinement steps with the same direction vector."); + oldMaximumLowerUpperBoundGap = weightVectorChecker->getMaximumLowerUpperBoundGap(); + std::vector lowerUpperDistances = stepIt->upperBoundPoint; + storm::utility::vector::subtractVectors(lowerUpperDistances, stepIt->lowerBoundPoint, lowerUpperDistances); + // shorten the distance between lower and upper bound for the new result by multiplying the current distance with 0.5 + // TODO: try other values/strategies? + GeometryValueType distance = storm::utility::sqrt(storm::utility::vector::dotProduct(lowerUpperDistances, lowerUpperDistances)); + weightVectorChecker->setMaximumLowerUpperBoundGap(std::min(*oldMaximumLowerUpperBoundGap, storm::utility::convertNumber(distance) * storm::utility::convertNumber(0.5))); + break; + } + } + weightVectorChecker->check(storm::utility::vector::convertNumericVector(direction)); + if(oldMaximumLowerUpperBoundGap) { + // Reset the precision back to the previous values + weightVectorChecker->setMaximumLowerUpperBoundGap(*oldMaximumLowerUpperBoundGap); + } + STORM_LOG_DEBUG("weighted objectives checker result (lower bounds) is " << storm::utility::vector::toString(storm::utility::vector::convertNumericVector(weightVectorChecker->getLowerBoundsOfInitialStateResults()))); + RefinementStep step; + step.weightVector = direction; + step.lowerBoundPoint = storm::utility::vector::convertNumericVector(weightVectorChecker->getLowerBoundsOfInitialStateResults()); + step.upperBoundPoint = storm::utility::vector::convertNumericVector(weightVectorChecker->getUpperBoundsOfInitialStateResults()); + refinementSteps.push_back(std::move(step)); + + updateOverApproximation(); + updateUnderApproximation(); + } + + template + void SparsePcaaQuery::updateOverApproximation() { + storm::storage::geometry::Halfspace h(refinementSteps.back().weightVector, storm::utility::vector::dotProduct(refinementSteps.back().weightVector, refinementSteps.back().upperBoundPoint)); + + // Due to numerical issues, it might be the case that the updated overapproximation does not contain the underapproximation, + // e.g., when the new point is strictly contained in the underapproximation. Check if this is the case. + GeometryValueType maximumOffset = h.offset(); + for(auto const& step : refinementSteps){ + maximumOffset = std::max(maximumOffset, storm::utility::vector::dotProduct(h.normalVector(), step.lowerBoundPoint)); + } + if(maximumOffset > h.offset()){ + // We correct the issue by shifting the halfspace such that it contains the underapproximation + h.offset() = maximumOffset; + STORM_LOG_WARN("Numerical issues: The overapproximation would not contain the underapproximation. Hence, a halfspace is shifted by " << storm::utility::convertNumber(h.euclideanDistance(refinementSteps.back().upperBoundPoint)) << "."); + } + overApproximation = overApproximation->intersection(h); + STORM_LOG_DEBUG("Updated OverApproximation to " << overApproximation->toString(true)); + } + + template + void SparsePcaaQuery::updateUnderApproximation() { + std::vector paretoPoints; + paretoPoints.reserve(refinementSteps.size()); + for(auto const& step : refinementSteps) { + paretoPoints.push_back(step.lowerBoundPoint); + } + underApproximation = storm::storage::geometry::Polytope::createDownwardClosure(paretoPoints); + STORM_LOG_DEBUG("Updated UnderApproximation to " << underApproximation->toString(true)); + } + + template + bool SparsePcaaQuery::maxStepsPerformed() const { + return storm::settings::getModule().isMaxStepsSet() && + this->refinementSteps.size() >= storm::settings::getModule().getMaxSteps(); + } + + + template + typename SparsePcaaQuery::Point SparsePcaaQuery::transformPointToOriginalModel(Point const& point) const { + Point result; + result.reserve(point.size()); + for(uint_fast64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) { + result.push_back(point[objIndex] * storm::utility::convertNumber(this->objectives[objIndex].toOriginalValueTransformationFactor) + storm::utility::convertNumber(this->objectives[objIndex].toOriginalValueTransformationOffset)); + } + return result; + } + + template + std::shared_ptr> SparsePcaaQuery::transformPolytopeToOriginalModel(std::shared_ptr> const& polytope) const { + if(polytope->isEmpty()) { + return storm::storage::geometry::Polytope::createEmptyPolytope(); + } + if(polytope->isUniversal()) { + return storm::storage::geometry::Polytope::createUniversalPolytope(); + } + uint_fast64_t numObjectives = this->objectives.size(); + std::vector> transformationMatrix(numObjectives, std::vector(numObjectives, storm::utility::zero())); + std::vector transformationVector; + transformationVector.reserve(numObjectives); + for(uint_fast64_t objIndex = 0; objIndex < numObjectives; ++objIndex) { + transformationMatrix[objIndex][objIndex] = storm::utility::convertNumber(this->objectives[objIndex].toOriginalValueTransformationFactor); + transformationVector.push_back(storm::utility::convertNumber(this->objectives[objIndex].toOriginalValueTransformationOffset)); + } + return polytope->affineTransformation(transformationMatrix, transformationVector); + } + +#ifdef STORM_HAVE_CARL + template class SparsePcaaQuery, storm::RationalNumber>; + template class SparsePcaaQuery, storm::RationalNumber>; + + template class SparsePcaaQuery, storm::RationalNumber>; + // template class SparsePcaaQuery, storm::RationalNumber>; +#endif + } + } +} diff --git a/src/modelchecker/multiobjective/pcaa/SparsePcaaQuery.h b/src/modelchecker/multiobjective/pcaa/SparsePcaaQuery.h new file mode 100644 index 000000000..4aba85b3b --- /dev/null +++ b/src/modelchecker/multiobjective/pcaa/SparsePcaaQuery.h @@ -0,0 +1,122 @@ +#ifndef STORM_MODELCHECKER_MULTIOBJECTIVE_PCAA_SPARSEPCAAQUERY_H_ +#define STORM_MODELCHECKER_MULTIOBJECTIVE_PCAA_SPARSEPCAAQUERY_H_ + +#include "src/modelchecker/results/CheckResult.h" +#include "src/modelchecker/multiobjective/pcaa/SparsePcaaPreprocessorReturnType.h" +#include "src/modelchecker/multiobjective/pcaa/SparsePcaaWeightVectorChecker.h" +#include "src/storage/geometry/Polytope.h" + +namespace storm { + namespace modelchecker { + namespace multiobjective { + + /* + * This class represents a query for the Pareto curve approximation algorithm (Pcaa). + * It implements the necessary computations for the different query types. + */ + template + class SparsePcaaQuery { + public: + + // Typedefs for simple geometric objects + typedef std::vector Point; + typedef std::vector WeightVector; + + /* + * Invokes the computation and retrieves the result + */ + virtual std::unique_ptr check() = 0; + + protected: + + /* + * Initializes the weight vector checker with the provided data from preprocessing + */ + void initializeWeightVectorChecker(SparseModelType const& model, + std::vector> const& objectives, + storm::storage::BitVector const& actionsWithNegativeReward, + storm::storage::BitVector const& ecActions, + storm::storage::BitVector const& possiblyRecurrentStates); + + /* + * Represents the information obtained in a single iteration of the algorithm + */ + struct RefinementStep { + WeightVector weightVector; + Point lowerBoundPoint; + Point upperBoundPoint; + }; + + /* + * Creates a new query for the Pareto curve approximation algorithm (Pcaa) + * @param preprocessorResult the result from preprocessing + */ + SparsePcaaQuery(SparsePcaaPreprocessorReturnType& preprocessorResult); + + /* + * Returns a weight vector w that separates the under approximation from the given point p, i.e., + * For each x in the under approximation, it holds that w*x <= w*p + * + * @param pointToBeSeparated the point that is to be seperated + */ + WeightVector findSeparatingVector(Point const& pointToBeSeparated); + + /* + * Refines the current result w.r.t. the given direction vector. + */ + void performRefinementStep(WeightVector&& direction); + + /* + * Updates the overapproximation after a refinement step has been performed + * + * @note The last entry of this->refinementSteps should be the newest step whose information is not yet included in the approximation. + */ + void updateOverApproximation(); + + /* + * Updates the underapproximation after a refinement step has been performed + * + * @note The last entry of this->refinementSteps should be the newest step whose information is not yet included in the approximation. + */ + void updateUnderApproximation(); + + /* + * Returns true iff the maximum number of refinement steps (as possibly specified in the settings) has been reached + */ + bool maxStepsPerformed() const; + + /* + * Transforms the given point (or polytope) to values w.r.t. the original model (e.g. negates negative rewards for minimizing objectives). + */ + Point transformPointToOriginalModel(Point const& polytope) const; + std::shared_ptr> transformPolytopeToOriginalModel(std::shared_ptr> const& polytope) const; + + + SparseModelType const& originalModel; + storm::logic::MultiObjectiveFormula const& originalFormula; + + SparseModelType preprocessedModel; + std::vector> objectives; + + + // The corresponding weight vector checker + std::unique_ptr> weightVectorChecker; + + //The results in each iteration of the algorithm + std::vector refinementSteps; + //Overapproximation of the set of achievable values + std::shared_ptr> overApproximation; + //Underapproximation of the set of achievable values + std::shared_ptr> underApproximation; + + // stores for each objective whether it still makes sense to check for this objective individually (i.e., with weight vector given by w_{i}>0 iff i=objIndex ) + storm::storage::BitVector diracWeightVectorsToBeChecked; + + + }; + + } + } +} + +#endif /* STORM_MODELCHECKER_MULTIOBJECTIVE_PCAA_SPARSEPCAAQUERY_H_ */ diff --git a/src/modelchecker/multiobjective/pcaa/SparsePcaaWeightVectorChecker.cpp b/src/modelchecker/multiobjective/pcaa/SparsePcaaWeightVectorChecker.cpp new file mode 100644 index 000000000..2ee474dd9 --- /dev/null +++ b/src/modelchecker/multiobjective/pcaa/SparsePcaaWeightVectorChecker.cpp @@ -0,0 +1,311 @@ +#include "src/modelchecker/multiobjective/pcaa/SparsePcaaWeightVectorChecker.h" + +#include + +#include "src/adapters/CarlAdapter.h" +#include "src/models/sparse/Mdp.h" +#include "src/models/sparse/MarkovAutomaton.h" +#include "src/models/sparse/StandardRewardModel.h" +#include "src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.h" +#include "src/solver/MinMaxLinearEquationSolver.h" +#include "src/transformer/EndComponentEliminator.h" +#include "src/utility/graph.h" +#include "src/utility/macros.h" +#include "src/utility/vector.h" + +#include "src/exceptions/IllegalFunctionCallException.h" +#include "src/exceptions/NotImplementedException.h" + +namespace storm { + namespace modelchecker { + namespace multiobjective { + + + template + SparsePcaaWeightVectorChecker::SparsePcaaWeightVectorChecker(SparseModelType const& model, + std::vector> const& objectives, + storm::storage::BitVector const& actionsWithNegativeReward, + storm::storage::BitVector const& ecActions, + storm::storage::BitVector const& possiblyRecurrentStates) : + model(model), + objectives(objectives), + actionsWithNegativeReward(actionsWithNegativeReward), + ecActions(ecActions), + possiblyRecurrentStates(possiblyRecurrentStates), + objectivesWithNoUpperTimeBound(objectives.size()), + discreteActionRewards(objectives.size()), + checkHasBeenCalled(false), + objectiveResults(objectives.size()), + offsetsToLowerBound(objectives.size()), + offsetsToUpperBound(objectives.size()) { + + // set the unbounded objectives + for(uint_fast64_t objIndex = 0; objIndex < objectives.size(); ++objIndex) { + objectivesWithNoUpperTimeBound.set(objIndex, !objectives[objIndex].upperTimeBound); + } + } + + + template + void SparsePcaaWeightVectorChecker::check(std::vector const& weightVector) { + checkHasBeenCalled=true; + STORM_LOG_DEBUG("Invoked WeightVectorChecker with weights " << std::endl << "\t" << storm::utility::vector::toString(storm::utility::vector::convertNumericVector(weightVector))); + std::vector weightedRewardVector(model.getTransitionMatrix().getRowCount(), storm::utility::zero()); + for(auto objIndex : objectivesWithNoUpperTimeBound) { + storm::utility::vector::addScaledVector(weightedRewardVector, discreteActionRewards[objIndex], weightVector[objIndex]); + } + unboundedWeightedPhase(weightedRewardVector); + unboundedIndividualPhase(weightVector); + // Only invoke boundedPhase if necessarry, i.e., if there is at least one objective with a time bound + for(auto const& obj : this->objectives) { + if(obj.lowerTimeBound || obj.upperTimeBound) { + boundedPhase(weightVector, weightedRewardVector); + break; + } + } + STORM_LOG_DEBUG("Weight vector check done. Lower bounds for results in initial state: " << storm::utility::vector::toString(storm::utility::vector::convertNumericVector(getLowerBoundsOfInitialStateResults()))); + } + + template + void SparsePcaaWeightVectorChecker::setMaximumLowerUpperBoundGap(ValueType const& value) { + this->maximumLowerUpperBoundGap = value; + } + + template + typename SparsePcaaWeightVectorChecker::ValueType const& SparsePcaaWeightVectorChecker::getMaximumLowerUpperBoundGap() const { + return this->maximumLowerUpperBoundGap; + } + + template + std::vector::ValueType> SparsePcaaWeightVectorChecker::getLowerBoundsOfInitialStateResults() const { + STORM_LOG_THROW(checkHasBeenCalled, storm::exceptions::IllegalFunctionCallException, "Tried to retrieve results but check(..) has not been called before."); + uint_fast64_t initstate = *this->model.getInitialStates().begin(); + std::vector res; + res.reserve(this->objectives.size()); + for(uint_fast64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) { + res.push_back(this->objectiveResults[objIndex][initstate] + this->offsetsToLowerBound[objIndex]); + } + return res; + } + + template + std::vector::ValueType> SparsePcaaWeightVectorChecker::getUpperBoundsOfInitialStateResults() const { + STORM_LOG_THROW(checkHasBeenCalled, storm::exceptions::IllegalFunctionCallException, "Tried to retrieve results but check(..) has not been called before."); + uint_fast64_t initstate = *this->model.getInitialStates().begin(); + std::vector res; + res.reserve(this->objectives.size()); + for(uint_fast64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) { + res.push_back(this->objectiveResults[objIndex][initstate] + this->offsetsToUpperBound[objIndex]); + } + return res; + } + + template + storm::storage::TotalScheduler const& SparsePcaaWeightVectorChecker::getScheduler() const { + STORM_LOG_THROW(this->checkHasBeenCalled, storm::exceptions::IllegalFunctionCallException, "Tried to retrieve results but check(..) has not been called before."); + for(auto const& obj : this->objectives) { + STORM_LOG_THROW(!obj.lowerTimeBound && !obj.upperTimeBound, storm::exceptions::NotImplementedException, "Scheduler retrival is not implemented for timeBounded objectives."); + } + return scheduler; + } + + template + void SparsePcaaWeightVectorChecker::unboundedWeightedPhase(std::vector const& weightedRewardVector) { + if(this->objectivesWithNoUpperTimeBound.empty()) { + this->weightedResult = std::vector(model.getNumberOfStates(), storm::utility::zero()); + this->scheduler = storm::storage::TotalScheduler(model.getNumberOfStates()); + return; + } + + + // Only consider the states from which a transition with non-zero reward is reachable. (The remaining states always have reward zero). + storm::storage::BitVector zeroRewardActions = storm::utility::vector::filterZero(weightedRewardVector); + storm::storage::BitVector nonZeroRewardActions = ~zeroRewardActions; + storm::storage::BitVector nonZeroRewardStates(model.getNumberOfStates(), false); + for(uint_fast64_t state = 0; state < model.getNumberOfStates(); ++state){ + if(nonZeroRewardActions.getNextSetIndex(model.getTransitionMatrix().getRowGroupIndices()[state]) < model.getTransitionMatrix().getRowGroupIndices()[state+1]) { + nonZeroRewardStates.set(state); + } + } + storm::storage::BitVector subsystemStates = storm::utility::graph::performProbGreater0E(model.getTransitionMatrix(), model.getTransitionMatrix().getRowGroupIndices(), model.getTransitionMatrix().transpose(true), storm::storage::BitVector(model.getNumberOfStates(), true), nonZeroRewardStates); + + // Remove neutral end components, i.e., ECs in which no reward is earned. + auto ecEliminatorResult = storm::transformer::EndComponentEliminator::transform(model.getTransitionMatrix(), subsystemStates, ecActions & zeroRewardActions, possiblyRecurrentStates); + + std::vector subRewardVector(ecEliminatorResult.newToOldRowMapping.size()); + storm::utility::vector::selectVectorValues(subRewardVector, ecEliminatorResult.newToOldRowMapping, weightedRewardVector); + std::vector subResult(ecEliminatorResult.matrix.getRowGroupCount()); + + storm::solver::GeneralMinMaxLinearEquationSolverFactory solverFactory; + std::unique_ptr> solver = solverFactory.create(ecEliminatorResult.matrix); + solver->setOptimizationDirection(storm::solver::OptimizationDirection::Maximize); + solver->setTrackScheduler(true); + solver->solveEquations(subResult, subRewardVector); + + this->weightedResult = std::vector(model.getNumberOfStates()); + std::vector optimalChoices(model.getNumberOfStates()); + + transformReducedSolutionToOriginalModel(ecEliminatorResult.matrix, subResult, solver->getScheduler()->getChoices(), ecEliminatorResult.newToOldRowMapping, ecEliminatorResult.oldToNewStateMapping, this->weightedResult, optimalChoices); + + this->scheduler = storm::storage::TotalScheduler(std::move(optimalChoices)); + } + + template + void SparsePcaaWeightVectorChecker::unboundedIndividualPhase(std::vector const& weightVector) { + + storm::storage::SparseMatrix deterministicMatrix = model.getTransitionMatrix().selectRowsFromRowGroups(this->scheduler.getChoices(), true); + storm::storage::SparseMatrix deterministicBackwardTransitions = deterministicMatrix.transpose(); + std::vector deterministicStateRewards(deterministicMatrix.getRowCount()); + storm::solver::GeneralLinearEquationSolverFactory linearEquationSolverFactory; + + // We compute an estimate for the results of the individual objectives which is obtained from the weighted result and the result of the objectives computed so far. + // Note that weightedResult = Sum_{i=1}^{n} w_i * objectiveResult_i. + std::vector weightedSumOfUncheckedObjectives = weightedResult; + ValueType sumOfWeightsOfUncheckedObjectives = storm::utility::vector::sum_if(weightVector, objectivesWithNoUpperTimeBound); + + for(uint_fast64_t const& objIndex : storm::utility::vector::getSortedIndices(weightVector)) { + if(objectivesWithNoUpperTimeBound.get(objIndex)){ + offsetsToLowerBound[objIndex] = storm::utility::zero(); + offsetsToUpperBound[objIndex] = storm::utility::zero(); + storm::utility::vector::selectVectorValues(deterministicStateRewards, this->scheduler.getChoices(), model.getTransitionMatrix().getRowGroupIndices(), discreteActionRewards[objIndex]); + storm::storage::BitVector statesWithRewards = ~storm::utility::vector::filterZero(deterministicStateRewards); + // As target states, we pick the states from which no reward is reachable. + storm::storage::BitVector targetStates = ~storm::utility::graph::performProbGreater0(deterministicBackwardTransitions, storm::storage::BitVector(deterministicMatrix.getRowCount(), true), statesWithRewards); + + // Compute the estimate for this objective + if(!storm::utility::isZero(weightVector[objIndex])) { + objectiveResults[objIndex] = weightedSumOfUncheckedObjectives; + storm::utility::vector::scaleVectorInPlace(objectiveResults[objIndex], storm::utility::one() / sumOfWeightsOfUncheckedObjectives); + } + + // Make sure that the objectiveResult is initialized in some way + objectiveResults[objIndex].resize(model.getNumberOfStates(), storm::utility::zero()); + + // Invoke the linear equation solver + objectiveResults[objIndex] = storm::modelchecker::helper::SparseDtmcPrctlHelper::computeReachabilityRewards(deterministicMatrix, + deterministicBackwardTransitions, + deterministicStateRewards, + targetStates, + false, //no qualitative checking, + linearEquationSolverFactory, + objectiveResults[objIndex]); + // Update the estimate for the next objectives. + if(!storm::utility::isZero(weightVector[objIndex])) { + storm::utility::vector::addScaledVector(weightedSumOfUncheckedObjectives, objectiveResults[objIndex], -weightVector[objIndex]); + sumOfWeightsOfUncheckedObjectives -= weightVector[objIndex]; + } + } else { + objectiveResults[objIndex] = std::vector(model.getNumberOfStates(), storm::utility::zero()); + } + } + } + + template + void SparsePcaaWeightVectorChecker::transformReducedSolutionToOriginalModel(storm::storage::SparseMatrix const& reducedMatrix, + std::vector const& reducedSolution, + std::vector const& reducedOptimalChoices, + std::vector const& reducedToOriginalChoiceMapping, + std::vector const& originalToReducedStateMapping, + std::vector& originalSolution, + std::vector& originalOptimalChoices) const { + + storm::storage::BitVector recurrentStates(model.getTransitionMatrix().getRowGroupCount(), false); + storm::storage::BitVector statesThatShouldStayInTheirEC(model.getTransitionMatrix().getRowGroupCount(), false); + storm::storage::BitVector statesWithUndefSched(model.getTransitionMatrix().getRowGroupCount(), false); + + // Handle all the states for which the choice in the original model is uniquely given by the choice in the reduced model + // Also store some information regarding the remaining states + for(uint_fast64_t state = 0; state < model.getTransitionMatrix().getRowGroupCount(); ++state) { + // Check if the state exists in the reduced model, i.e., the mapping retrieves a valid index + uint_fast64_t stateInReducedModel = originalToReducedStateMapping[state]; + if(stateInReducedModel < reducedMatrix.getRowGroupCount()) { + originalSolution[state] = reducedSolution[stateInReducedModel]; + uint_fast64_t chosenRowInReducedModel = reducedMatrix.getRowGroupIndices()[stateInReducedModel] + reducedOptimalChoices[stateInReducedModel]; + uint_fast64_t chosenRowInOriginalModel = reducedToOriginalChoiceMapping[chosenRowInReducedModel]; + // Check if the state is recurrent, i.e., the chosen row stays inside this EC. + bool stateIsRecurrent = possiblyRecurrentStates.get(state); + for(auto const& entry : model.getTransitionMatrix().getRow(chosenRowInOriginalModel)) { + stateIsRecurrent &= originalToReducedStateMapping[entry.getColumn()] == stateInReducedModel; + } + if(stateIsRecurrent) { + recurrentStates.set(state); + statesThatShouldStayInTheirEC.set(state); + } else { + // Check if the chosen row originaly belonged to the current state (and not to another state of the EC) + if(chosenRowInOriginalModel >= model.getTransitionMatrix().getRowGroupIndices()[state] && + chosenRowInOriginalModel < model.getTransitionMatrix().getRowGroupIndices()[state+1]) { + originalOptimalChoices[state] = chosenRowInOriginalModel - model.getTransitionMatrix().getRowGroupIndices()[state]; + } else { + statesWithUndefSched.set(state); + statesThatShouldStayInTheirEC.set(state); + } + } + } else { + // if the state does not exist in the reduced model, it means that the (weighted) result is always zero, independent of the scheduler. + originalSolution[state] = storm::utility::zero(); + // However, it might be the case that infinite reward is induced for an objective with weight 0. + // To avoid this, all possibly recurrent states are made recurrent and the remaining states have to reach a recurrent state with prob. one + if(possiblyRecurrentStates.get(state)) { + recurrentStates.set(state); + } else { + statesWithUndefSched.set(state); + } + } + } + + // Handle recurrent states + for(auto state : recurrentStates) { + bool foundRowForState = false; + // Find a row with zero rewards that only leads to recurrent states. + // If the state should stay in its EC, we also need to make sure that all successors map to the same state in the reduced model + uint_fast64_t stateInReducedModel = originalToReducedStateMapping[state]; + for(uint_fast64_t row = model.getTransitionMatrix().getRowGroupIndices()[state]; row < model.getTransitionMatrix().getRowGroupIndices()[state+1]; ++row) { + bool rowOnlyLeadsToRecurrentStates = true; + bool rowStaysInEC = true; + for( auto const& entry : model.getTransitionMatrix().getRow(row)) { + rowOnlyLeadsToRecurrentStates &= recurrentStates.get(entry.getColumn()); + rowStaysInEC &= originalToReducedStateMapping[entry.getColumn()] == stateInReducedModel; + } + if(rowOnlyLeadsToRecurrentStates && (rowStaysInEC || !statesThatShouldStayInTheirEC.get(state)) && !actionsWithNegativeReward.get(row)) { + foundRowForState = true; + originalOptimalChoices[state] = row - model.getTransitionMatrix().getRowGroupIndices()[state]; + break; + } + } + STORM_LOG_ASSERT(foundRowForState, "Could not find a suitable choice for a recurrent state."); + } + + // Handle remaining states with still undef. scheduler (either EC states or non-subsystem states) + while(!statesWithUndefSched.empty()) { + for(auto state : statesWithUndefSched) { + // Try to find a choice such that at least one successor has a defined scheduler. + // This way, a non-recurrent state will never become recurrent + uint_fast64_t stateInReducedModel = originalToReducedStateMapping[state]; + for(uint_fast64_t row = model.getTransitionMatrix().getRowGroupIndices()[state]; row < model.getTransitionMatrix().getRowGroupIndices()[state+1]; ++row) { + bool rowStaysInEC = true; + bool rowLeadsToDefinedScheduler = false; + for(auto const& entry : model.getTransitionMatrix().getRow(row)) { + rowStaysInEC &= ( stateInReducedModel == originalToReducedStateMapping[entry.getColumn()]); + rowLeadsToDefinedScheduler |= !statesWithUndefSched.get(entry.getColumn()); + } + if(rowLeadsToDefinedScheduler && (rowStaysInEC || !statesThatShouldStayInTheirEC.get(state))) { + originalOptimalChoices[state] = row - model.getTransitionMatrix().getRowGroupIndices()[state]; + statesWithUndefSched.set(state, false); + } + } + } + } + } + + + + template class SparsePcaaWeightVectorChecker>; + template class SparsePcaaWeightVectorChecker>; +#ifdef STORM_HAVE_CARL + template class SparsePcaaWeightVectorChecker>; + template class SparsePcaaWeightVectorChecker>; +#endif + + } + } +} diff --git a/src/modelchecker/multiobjective/pcaa/SparsePcaaWeightVectorChecker.h b/src/modelchecker/multiobjective/pcaa/SparsePcaaWeightVectorChecker.h new file mode 100644 index 000000000..df69ec059 --- /dev/null +++ b/src/modelchecker/multiobjective/pcaa/SparsePcaaWeightVectorChecker.h @@ -0,0 +1,155 @@ +#ifndef STORM_MODELCHECKER_MULTIOBJECTIVE_PCAA_SPARSEPCAAWEIGHTVECTORCHECKER_H_ +#define STORM_MODELCHECKER_MULTIOBJECTIVE_PCAA_SPARSEPCAAWEIGHTVECTORCHECKER_H_ + + +#include "src/storage/BitVector.h" +#include "src/storage/SparseMatrix.h" +#include "src/storage/TotalScheduler.h" +#include "src/modelchecker/multiobjective/pcaa/PcaaObjective.h" +#include "src/utility/vector.h" + +namespace storm { + namespace modelchecker { + namespace multiobjective { + + /*! + * Helper Class that takes preprocessed Pcaa data and a weight vector and ... + * - computes the maximal expected reward w.r.t. the weighted sum of the rewards of the individual objectives + * - extracts the scheduler that induces this maximum + * - computes for each objective the value induced by this scheduler + */ + template + class SparsePcaaWeightVectorChecker { + public: + typedef typename SparseModelType::ValueType ValueType; + + /* + * Creates a weight vextor checker. + * + * @param model The (preprocessed) model + * @param objectives The (preprocessed) objectives + * @param actionsWithNegativeReward The actions that have negative reward assigned for at least one objective + * @param ecActions The actions that are part of an EC + * @param possiblyRecurrentStates The states for which it is posible to visit them infinitely often (without inducing inf. reward) + * + */ + + SparsePcaaWeightVectorChecker(SparseModelType const& model, + std::vector> const& objectives, + storm::storage::BitVector const& actionsWithNegativeReward, + storm::storage::BitVector const& ecActions, + storm::storage::BitVector const& possiblyRecurrentStates); + + /*! + * - computes the maximal expected reward w.r.t. the weighted sum of the rewards of the individual objectives + * - extracts the scheduler that induces this maximum + * - computes for each objective the value induced by this scheduler + */ + void check(std::vector const& weightVector); + + /*! + * Sets the maximum gap that is allowed between the lower and upper bound of the result of some objective. + */ + void setMaximumLowerUpperBoundGap(ValueType const& value); + + /*! + * Retrieves the maximum gap that is allowed between the lower and upper bound of the result of some objective. + */ + ValueType const& getMaximumLowerUpperBoundGap() const; + + /*! + * Retrieves the results of the individual objectives at the initial state of the given model. + * Note that check(..) has to be called before retrieving results. Otherwise, an exception is thrown. + * Also note that there is no guarantee that the lower/upper bounds are sound + * as long as the underlying solution methods are unsound (e.g., standard value iteration). + */ + std::vector getLowerBoundsOfInitialStateResults() const; + std::vector getUpperBoundsOfInitialStateResults() const; + + /*! + * Retrieves a scheduler that induces the current values + * Note that check(..) has to be called before retrieving the scheduler. Otherwise, an exception is thrown. + */ + storm::storage::TotalScheduler const& getScheduler() const; + + + protected: + + /*! + * Determines the scheduler that maximizes the weighted reward vector of the unbounded objectives + * + * @param weightedRewardVector the weighted rewards (only considering the unbounded objectives) + */ + void unboundedWeightedPhase(std::vector const& weightedRewardVector); + + /*! + * Computes the values of the objectives that do not have a stepBound w.r.t. the scheduler computed in the unboundedWeightedPhase + * + * @param weightVector the weight vector of the current check + */ + void unboundedIndividualPhase(std::vector const& weightVector); + + /*! + * For each time epoch (starting with the maximal stepBound occurring in the objectives), this method + * - determines the objectives that are relevant in the current time epoch + * - determines the maximizing scheduler for the weighted reward vector of these objectives + * - computes the values of these objectives w.r.t. this scheduler + * + * @param weightVector the weight vector of the current check + * @param weightedRewardVector the weighted rewards considering the unbounded objectives. Will be invalidated after calling this. + */ + virtual void boundedPhase(std::vector const& weightVector, std::vector& weightedRewardVector) = 0; + + /*! + * Transforms the results of a min-max-solver that considers a reduced model (without end components) to a result for the original (unreduced) model + */ + void transformReducedSolutionToOriginalModel(storm::storage::SparseMatrix const& reducedMatrix, + std::vector const& reducedSolution, + std::vector const& reducedOptimalChoices, + std::vector const& reducedToOriginalChoiceMapping, + std::vector const& originalToReducedStateMapping, + std::vector& originalSolution, + std::vector& originalOptimalChoices) const; + + // The (preprocessed) model + SparseModelType const& model; + // The (preprocessed) objectives + std::vector> const& objectives; + + // The actions that have negative reward assigned for at least one objective + storm::storage::BitVector actionsWithNegativeReward; + // The actions that are part of an EC + storm::storage::BitVector ecActions; + // The states for which it is allowed to visit them infinitely often + // Put differently, if one of the states is part of a neutral EC, it is possible to + // stay in this EC forever (withoud inducing infinite reward for some objective). + storm::storage::BitVector possiblyRecurrentStates; + // stores the indices of the objectives for which there is no upper time bound + storm::storage::BitVector objectivesWithNoUpperTimeBound; + // stores the (discretized) state action rewards for each objective. + std::vector> discreteActionRewards; + // stores the maximum gap that is allowed between the lower and upper bound of the result of some objective. + ValueType maximumLowerUpperBoundGap; + + // Memory for the solution of the most recent call of check(..) + // becomes true after the first call of check(..) + bool checkHasBeenCalled; + // The result for the weighted reward vector (for all states of the model) + std::vector weightedResult; + // The lower bounds of the results for the individual objectives (w.r.t. all states of the model) + std::vector> objectiveResults; + // Stores for each objective the distance between the computed result (w.r.t. the initial state) and a lower/upper bound for the actual result. + // The distances are stored as a (possibly negative) offset that has to be added to to the objectiveResults. + // Note that there is no guarantee that the lower/upper bounds are sound as long as the underlying solution method is not sound (e.g. standard value iteration). + std::vector offsetsToLowerBound; + std::vector offsetsToUpperBound; + // The scheduler that maximizes the weighted rewards + storm::storage::TotalScheduler scheduler; + + }; + + } + } +} + +#endif /* STORM_MODELCHECKER_MULTIOBJECTIVE_PCAA_SPARSEPCAAWEIGHTVECTORCHECKER_H_ */ diff --git a/src/transformer/EndComponentEliminator.h b/src/transformer/EndComponentEliminator.h index d60dd9c2c..a73c5cf0e 100644 --- a/src/transformer/EndComponentEliminator.h +++ b/src/transformer/EndComponentEliminator.h @@ -18,30 +18,35 @@ namespace storm { // The resulting matrix storm::storage::SparseMatrix matrix; // Index mapping that gives for each row of the resulting matrix the corresponding row in the original matrix. - // For the empty rows added to EC states, some row of the original matrix that stays inside the EC is given. + // For the empty rows added to EC states, an arbitrary row of the original matrix that stays inside the EC is given. std::vector newToOldRowMapping; // Gives for each state (=rowGroup) of the original matrix the corresponding state in the resulting matrix. - // States of a removed end component are mapped to the state that substitutes the EC. - // If the given state does not exist in the result (e.g. it belonged to a bottom EC), the returned value will be std::numeric_limits::max(), i.e., an invalid index. + // States of a removed ECs are mapped to the state that substitutes the EC. + // If the given state does not exist in the result (i.e., it is not in the provided subsystem), the returned value will be std::numeric_limits::max(), i.e., an invalid index. std::vector oldToNewStateMapping; }; /* * Identifies end components and substitutes them by a single state. * - * Only ECs for which the given bit vector is true for all choices are considered. - * For each such EC (that is not contained by another EC), we add a new state and redirect all incoming and outgoing + * Only states in the given subsystem are kept. Transitions leading to a state outside of the subsystem will be + * removed (but the corresponding row is kept, possibly yielding empty rows). + * The ECs are then identified on the subsystem. + * + * Only ECs for which possibleECRows is true for all choices are considered. + * Furthermore, the rows that contain a transition leading outside of the subsystem are not considered for an EC. + * + * For each such EC (that is not contained in another EC), we add a new state and redirect all incoming and outgoing * transitions of the EC to (and from) this state. - * If allowEmptyRow is true for at least one state of an EC, an empty row is added to the new state (representing the choice to stay at the EC forever). - * States for which all reachable choices are selected by the given bit vector will be removed! + * If addEmptyRowStates is true for at least one state of an eliminated EC, an empty row is added to the new state (representing the choice to stay at the EC forever). */ - static EndComponentEliminatorReturnType transform(storm::storage::SparseMatrix const& originalMatrix, storm::storage::BitVector const& consideredRows, storm::storage::BitVector const& allowEmptyRow) { + static EndComponentEliminatorReturnType transform(storm::storage::SparseMatrix const& originalMatrix, storm::storage::BitVector const& subsystemStates, storm::storage::BitVector const& possibleECRows, storm::storage::BitVector const& addEmptyRowStates) { STORM_LOG_DEBUG("Invoked EndComponentEliminator on matrix with " << originalMatrix.getRowGroupCount() << " row groups."); - storm::storage::BitVector keptStates = computeStatesFromWhichNonConsideredRowIsReachable(originalMatrix, consideredRows); - storm::storage::MaximalEndComponentDecomposition ecs = computeECs(originalMatrix, consideredRows, keptStates); + storm::storage::MaximalEndComponentDecomposition ecs = computeECs(originalMatrix, possibleECRows, subsystemStates); //further shrink the set of kept states by removing all states that are part of an EC + storm::storage::BitVector keptStates = subsystemStates; for (auto const& ec : ecs) { for (auto const& stateActionsPair : ec) { keptStates.set(stateActionsPair.first, false); @@ -71,7 +76,7 @@ namespace storm { result.newToOldRowMapping.push_back(row); } } - ecGetsEmptyRow |= allowEmptyRow.get(stateActionsPair.first); + ecGetsEmptyRow |= addEmptyRowStates.get(stateActionsPair.first); } if(ecGetsEmptyRow) { STORM_LOG_ASSERT(result.newToOldRowMapping.size() < originalMatrix.getRowCount(), "Didn't expect to see more rows in the reduced matrix than in the original one."); @@ -89,31 +94,20 @@ namespace storm { private: - static storm::storage::BitVector computeStatesFromWhichNonConsideredRowIsReachable(storm::storage::SparseMatrix const& originalMatrix, storm::storage::BitVector const& consideredRows) { - storm::storage::BitVector nonConsideredRows = ~consideredRows; - storm::storage::BitVector groupsWithNonConsideredRow(originalMatrix.getRowGroupCount(), false); - for(uint_fast64_t rowGroup = 0; rowGroup < originalMatrix.getRowGroupCount(); ++rowGroup){ - if(nonConsideredRows.getNextSetIndex(originalMatrix.getRowGroupIndices()[rowGroup]) < originalMatrix.getRowGroupIndices()[rowGroup+1]) { - groupsWithNonConsideredRow.set(rowGroup); - } - } - return storm::utility::graph::performProbGreater0E(originalMatrix, originalMatrix.getRowGroupIndices(), originalMatrix.transpose(true), storm::storage::BitVector(originalMatrix.getRowGroupCount(), true), groupsWithNonConsideredRow); - } - - static storm::storage::MaximalEndComponentDecomposition computeECs(storm::storage::SparseMatrix const& originalMatrix, storm::storage::BitVector const& consideredRows, storm::storage::BitVector const& keptStates) { - // Get an auxiliary matrix to identify the correct end components w.r.t. the considered choices and the kept states. + static storm::storage::MaximalEndComponentDecomposition computeECs(storm::storage::SparseMatrix const& originalMatrix, storm::storage::BitVector const& possibleECRows, storm::storage::BitVector const& subsystemStates) { + // Get an auxiliary matrix to identify the correct end components w.r.t. the possible EC rows and the subsystem. // This is done by redirecting choices that can never be part of an EC to a sink state. - // Such choices are either non-considered choices or choices that lead to a state that is not in keptStates. + // Such choices are either non-possible EC rows or rows that lead to a state that is not in the subsystem. uint_fast64_t sinkState = originalMatrix.getRowGroupCount(); storm::storage::SparseMatrixBuilder builder(originalMatrix.getRowCount() + 1, originalMatrix.getColumnCount() + 1, originalMatrix.getEntryCount() + 1, false, true, originalMatrix.getRowGroupCount()+1); uint_fast64_t row = 0; for(uint_fast64_t rowGroup = 0; rowGroup < originalMatrix.getRowGroupCount(); ++rowGroup) { builder.newRowGroup(row); for (; row < originalMatrix.getRowGroupIndices()[rowGroup+1]; ++row ){ - bool keepRow = consideredRows.get(row); - if(keepRow) { //Also check whether all successors are kept + bool keepRow = possibleECRows.get(row); + if(keepRow) { //Also check whether all successors are in the subsystem for(auto const& entry : originalMatrix.getRow(row)){ - keepRow &= keptStates.get(entry.getColumn()); + keepRow &= subsystemStates.get(entry.getColumn()); } } if(keepRow) { @@ -131,11 +125,11 @@ namespace storm { storm::storage::SparseMatrix backwardsTransitions = auxiliaryMatrix.transpose(true); storm::storage::BitVector sinkStateAsBitVector(auxiliaryMatrix.getRowGroupCount(), false); sinkStateAsBitVector.set(sinkState); - storm::storage::BitVector subsystem = keptStates; - subsystem.resize(keptStates.size() + 1, true); - // The states for which sinkState is reachable under any scheduler can not be part of an EC - subsystem &= ~(storm::utility::graph::performProbGreater0A(auxiliaryMatrix, auxiliaryMatrix.getRowGroupIndices(), backwardsTransitions, subsystem, sinkStateAsBitVector)); - return storm::storage::MaximalEndComponentDecomposition(auxiliaryMatrix, backwardsTransitions, subsystem); + storm::storage::BitVector auxSubsystemStates = subsystemStates; + auxSubsystemStates.resize(subsystemStates.size() + 1, true); + // The states for which sinkState is reachable under every scheduler can not be part of an EC + auxSubsystemStates &= ~(storm::utility::graph::performProbGreater0A(auxiliaryMatrix, auxiliaryMatrix.getRowGroupIndices(), backwardsTransitions, auxSubsystemStates, sinkStateAsBitVector)); + return storm::storage::MaximalEndComponentDecomposition(auxiliaryMatrix, backwardsTransitions, auxSubsystemStates); } static storm::storage::SparseMatrix buildTransformedMatrix(storm::storage::SparseMatrix const& originalMatrix, diff --git a/test/functional/transformer/EndComponentEliminatorTest.cpp b/test/functional/transformer/EndComponentEliminatorTest.cpp index a4ce99a39..8785df336 100644 --- a/test/functional/transformer/EndComponentEliminatorTest.cpp +++ b/test/functional/transformer/EndComponentEliminatorTest.cpp @@ -34,7 +34,7 @@ TEST(NeutralECRemover, SimpleModelTest) { storm::storage::SparseMatrix matrix; ASSERT_NO_THROW(matrix = builder.build()); - storm::storage::BitVector consideredRows(12, true); + storm::storage::BitVector possibleEcRows(12, true); consideredRows.set(3, false); consideredRows.set(6, false); consideredRows.set(7, false); @@ -45,7 +45,7 @@ TEST(NeutralECRemover, SimpleModelTest) { allowEmptyRows.set(1, false); allowEmptyRows.set(4, false); - auto res = storm::transformer::EndComponentEliminator::transform(matrix, consideredRows, allowEmptyRows); + auto res = storm::transformer::EndComponentEliminator::transform(matrix, possibleEcRows, allowEmptyRows); // Expected data storm::storage::SparseMatrixBuilder expectedBuilder(8, 3, 8, true, true, 3);