diff --git a/src/modelchecker/multiobjective/SparseMdpMultiObjectiveModelChecker.cpp b/src/modelchecker/multiobjective/SparseMdpMultiObjectiveModelChecker.cpp index 5cd49308e..145b3b4d8 100644 --- a/src/modelchecker/multiobjective/SparseMdpMultiObjectiveModelChecker.cpp +++ b/src/modelchecker/multiobjective/SparseMdpMultiObjectiveModelChecker.cpp @@ -9,8 +9,9 @@ #include "src/modelchecker/results/ExplicitQualitativeCheckResult.h" #include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h" -#include "src/modelchecker/multiobjective/helper/SparseMdpMultiObjectivePreprocessingHelper.h" #include "src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerInformation.h" +#include "src/modelchecker/multiobjective/helper/SparseMdpMultiObjectivePreprocessingHelper.h" +#include "src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerHelper.h" #include "src/exceptions/NotImplementedException.h" @@ -35,11 +36,19 @@ namespace storm { std::unique_ptr SparseMdpMultiObjectiveModelChecker::checkMultiObjectiveFormula(CheckTask const& checkTask) { helper::SparseMultiObjectiveModelCheckerInformation info = helper::SparseMdpMultiObjectivePreprocessingHelper::preprocess(checkTask.getFormula(), this->getModel()); - std::cout << std::endl; + std::cout << "Preprocessed Information:" << std::endl; info.printInformationToStream(std::cout); +#ifdef STORM_HAVE_CARL + helper::SparseMultiObjectiveModelCheckerHelper helper(info); + helper.achievabilityQuery(); +#else + STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "Multi objective model checking currently requires carl."); +#endif + std::cout << "Information after helper call: " << std::endl; + info.printInformationToStream(std::cout); return std::unique_ptr(new ExplicitQualitativeCheckResult()); } diff --git a/src/modelchecker/multiobjective/helper/SparseMdpMultiObjectivePreprocessingHelper.cpp b/src/modelchecker/multiobjective/helper/SparseMdpMultiObjectivePreprocessingHelper.cpp index 7b438cba2..c5b44cf2b 100644 --- a/src/modelchecker/multiobjective/helper/SparseMdpMultiObjectivePreprocessingHelper.cpp +++ b/src/modelchecker/multiobjective/helper/SparseMdpMultiObjectivePreprocessingHelper.cpp @@ -22,10 +22,10 @@ namespace storm { namespace helper { template - SparseMultiObjectiveModelCheckerInformation SparseMdpMultiObjectivePreprocessingHelper::preprocess(storm::logic::MultiObjectiveFormula const& originalFormula, SparseMdpModelType originalModel) { + SparseMultiObjectiveModelCheckerInformation SparseMdpMultiObjectivePreprocessingHelper::preprocess(storm::logic::MultiObjectiveFormula const& originalFormula, SparseMdpModelType const& originalModel) { Information info(std::move(originalModel)); //Initialize the state mapping. - info.newToOldStateIndexMapping = storm::utility::vector::buildVectorForRange(0, info.model.getNumberOfStates()); + info.newToOldStateIndexMapping = storm::utility::vector::buildVectorForRange(0, info.preprocessedModel.getNumberOfStates()); //Gather information regarding the individual objectives for(auto const& subFormula : originalFormula.getSubFormulas()){ addObjective(subFormula, info); @@ -46,7 +46,7 @@ namespace storm { //We can now remove all original reward models to save some memory std::set origRewardModels = originalFormula.getReferencedRewardModels(); for (auto const& rewModel : origRewardModels){ - info.model.removeRewardModel(rewModel); + info.preprocessedModel.removeRewardModel(rewModel); } return info; } @@ -70,7 +70,7 @@ namespace storm { } objective.rewardModelName = "objective" + std::to_string(info.objectives.size()); // Make sure the new reward model gets a unique name - while(info.model.hasRewardModel(objective.rewardModelName)){ + while(info.preprocessedModel.hasRewardModel(objective.rewardModelName)){ objective.rewardModelName = "_" + objective.rewardModelName; } setStepBoundOfObjective(objective); @@ -145,8 +145,8 @@ namespace storm { STORM_LOG_THROW(info.negatedRewardsConsidered == currentObjective.originalFormulaMinimizes, storm::exceptions::InvalidPropertyException, "Decided to consider " << (info.negatedRewardsConsidered ? "negated" : "non-negated") << "rewards, but the formula " << formula << (currentObjective.originalFormulaMinimizes ? " minimizes" : "maximizes") << ". Reward objectives should either all minimize or all maximize."); // Check if the reward model is uniquely specified - STORM_LOG_THROW((formula.hasRewardModelName() && info.model.hasRewardModel(formula.getRewardModelName())) - || info.model.hasUniqueRewardModel(), storm::exceptions::InvalidPropertyException, "The reward model is not unique and the formula " << formula << " does not specify a reward model."); + STORM_LOG_THROW((formula.hasRewardModelName() && info.preprocessedModel.hasRewardModel(formula.getRewardModelName())) + || info.preprocessedModel.hasUniqueRewardModel(), storm::exceptions::InvalidPropertyException, "The reward model is not unique and the formula " << formula << " does not specify a reward model."); if(info.negatedRewardsConsidered && currentObjective.threshold){ *(currentObjective.threshold) *= -storm::utility::one(); @@ -166,13 +166,13 @@ namespace storm { void SparseMdpMultiObjectivePreprocessingHelper::preprocessFormula(storm::logic::UntilFormula const& formula, Information& info, typename Information::ObjectiveInformation& currentObjective) { CheckTask phiTask(formula.getLeftSubformula()); CheckTask psiTask(formula.getRightSubformula()); - storm::modelchecker::SparsePropositionalModelChecker mc(info.model); + storm::modelchecker::SparsePropositionalModelChecker mc(info.preprocessedModel); STORM_LOG_THROW(mc.canHandle(phiTask) && mc.canHandle(psiTask), storm::exceptions::InvalidPropertyException, "The subformulas of " << formula << " should be propositional."); storm::storage::BitVector phiStates = mc.check(phiTask)->asExplicitQualitativeCheckResult().getTruthValuesVector(); storm::storage::BitVector psiStates = mc.check(psiTask)->asExplicitQualitativeCheckResult().getTruthValuesVector(); - auto duplicatorResult = storm::transformer::StateDuplicator::transform(info.model, ~phiStates | psiStates); - info.model = std::move(*duplicatorResult.model); + auto duplicatorResult = storm::transformer::StateDuplicator::transform(info.preprocessedModel, ~phiStates | psiStates); + info.preprocessedModel = std::move(*duplicatorResult.model); // duplicatorResult.newToOldStateIndexMapping now reffers to the indices of the model we had before preprocessing this formula. // This might not be the actual original model. for(auto & originalModelStateIndex : duplicatorResult.newToOldStateIndexMapping){ @@ -181,21 +181,21 @@ namespace storm { info.newToOldStateIndexMapping = std::move(duplicatorResult.newToOldStateIndexMapping); // build stateAction reward vector that gives (one*transitionProbability) reward whenever a transition leads from the firstCopy to a psiState - storm::storage::BitVector newPsiStates(info.model.getNumberOfStates(), false); + storm::storage::BitVector newPsiStates(info.preprocessedModel.getNumberOfStates(), false); for(auto const& oldPsiState : psiStates){ //note that psiStates are always located in the second copy newPsiStates.set(duplicatorResult.secondCopyOldToNewStateIndexMapping[oldPsiState], true); } - std::vector objectiveRewards(info.model.getTransitionMatrix().getRowCount(), storm::utility::zero()); + std::vector objectiveRewards(info.preprocessedModel.getTransitionMatrix().getRowCount(), storm::utility::zero()); for (auto const& firstCopyState : duplicatorResult.firstCopy) { - for (uint_fast64_t row = info.model.getTransitionMatrix().getRowGroupIndices()[firstCopyState]; row < info.model.getTransitionMatrix().getRowGroupIndices()[firstCopyState + 1]; ++row) { - objectiveRewards[row] = info.model.getTransitionMatrix().getConstrainedRowSum(row, newPsiStates); + for (uint_fast64_t row = info.preprocessedModel.getTransitionMatrix().getRowGroupIndices()[firstCopyState]; row < info.preprocessedModel.getTransitionMatrix().getRowGroupIndices()[firstCopyState + 1]; ++row) { + objectiveRewards[row] = info.preprocessedModel.getTransitionMatrix().getConstrainedRowSum(row, newPsiStates); } } if(info.negatedRewardsConsidered){ storm::utility::vector::scaleVectorInPlace(objectiveRewards, -storm::utility::one()); } - info.model.addRewardModel(currentObjective.rewardModelName, RewardModelType(boost::none, objectiveRewards)); + info.preprocessedModel.addRewardModel(currentObjective.rewardModelName, RewardModelType(boost::none, objectiveRewards)); } template @@ -218,13 +218,13 @@ namespace storm { STORM_LOG_THROW(formula.isReachabilityRewardFormula(), storm::exceptions::InvalidPropertyException, "The formula " << formula << " neither considers reachability Probabilities nor reachability rewards"); CheckTask targetTask(formula.getSubformula()); - storm::modelchecker::SparsePropositionalModelChecker mc(info.model); + storm::modelchecker::SparsePropositionalModelChecker mc(info.preprocessedModel); STORM_LOG_THROW(mc.canHandle(targetTask), storm::exceptions::InvalidPropertyException, "The subformula of " << formula << " should be propositional."); storm::storage::BitVector targetStates = mc.check(targetTask)->asExplicitQualitativeCheckResult().getTruthValuesVector(); STORM_LOG_WARN("There is no check for reward finiteness yet"); //TODO - auto duplicatorResult = storm::transformer::StateDuplicator::transform(info.model, targetStates); - info.model = std::move(*duplicatorResult.model); + auto duplicatorResult = storm::transformer::StateDuplicator::transform(info.preprocessedModel, targetStates); + info.preprocessedModel = std::move(*duplicatorResult.model); // duplicatorResult.newToOldStateIndexMapping now reffers to the indices of the model we had before preprocessing this formula. // This might not be the actual original model. for(auto & originalModelStateIndex : duplicatorResult.newToOldStateIndexMapping){ @@ -233,21 +233,21 @@ namespace storm { info.newToOldStateIndexMapping = std::move(duplicatorResult.newToOldStateIndexMapping); // Add a reward model that gives zero reward to the states of the second copy. - std::vector objectiveRewards = info.model.getRewardModel(optionalRewardModelName ? optionalRewardModelName.get() : "").getTotalRewardVector(info.model.getTransitionMatrix()); + std::vector objectiveRewards = info.preprocessedModel.getRewardModel(optionalRewardModelName ? optionalRewardModelName.get() : "").getTotalRewardVector(info.preprocessedModel.getTransitionMatrix()); storm::utility::vector::setVectorValues(objectiveRewards, duplicatorResult.secondCopy, storm::utility::zero()); if(info.negatedRewardsConsidered){ storm::utility::vector::scaleVectorInPlace(objectiveRewards, -storm::utility::one()); } - info.model.addRewardModel(currentObjective.rewardModelName, RewardModelType(boost::none, objectiveRewards)); + info.preprocessedModel.addRewardModel(currentObjective.rewardModelName, RewardModelType(boost::none, objectiveRewards)); } template void SparseMdpMultiObjectivePreprocessingHelper::preprocessFormula(storm::logic::CumulativeRewardFormula const& formula, Information& info, typename Information::ObjectiveInformation& currentObjective, boost::optional const& optionalRewardModelName) { - std::vector objectiveRewards = info.model.getRewardModel(optionalRewardModelName ? optionalRewardModelName.get() : "").getTotalRewardVector(info.model.getTransitionMatrix()); + std::vector objectiveRewards = info.preprocessedModel.getRewardModel(optionalRewardModelName ? optionalRewardModelName.get() : "").getTotalRewardVector(info.preprocessedModel.getTransitionMatrix()); if(info.negatedRewardsConsidered){ storm::utility::vector::scaleVectorInPlace(objectiveRewards, -storm::utility::one()); } - info.model.addRewardModel(currentObjective.rewardModelName, RewardModelType(boost::none, objectiveRewards)); + info.preprocessedModel.addRewardModel(currentObjective.rewardModelName, RewardModelType(boost::none, objectiveRewards)); } template class SparseMdpMultiObjectivePreprocessingHelper>; diff --git a/src/modelchecker/multiobjective/helper/SparseMdpMultiObjectivePreprocessingHelper.h b/src/modelchecker/multiobjective/helper/SparseMdpMultiObjectivePreprocessingHelper.h index 5adf1bda9..cc1ea8925 100644 --- a/src/modelchecker/multiobjective/helper/SparseMdpMultiObjectivePreprocessingHelper.h +++ b/src/modelchecker/multiobjective/helper/SparseMdpMultiObjectivePreprocessingHelper.h @@ -25,7 +25,7 @@ namespace storm { * @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 Information preprocess(storm::logic::MultiObjectiveFormula const& originalFormula, SparseMdpModelType originalModel); + static Information preprocess(storm::logic::MultiObjectiveFormula const& originalFormula, SparseMdpModelType const& originalModel); private: diff --git a/src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerHelper.cpp b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerHelper.cpp new file mode 100644 index 000000000..28eb326bd --- /dev/null +++ b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerHelper.cpp @@ -0,0 +1,137 @@ +#include "src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerHelper.h" + +#include "src/adapters/CarlAdapter.h" +#include "src/models/sparse/Mdp.h" +#include "src/models/sparse/StandardRewardModel.h" +#include "src/utility/constants.h" + +#include "src/exceptions/UnexpectedException.h" + +namespace storm { + namespace modelchecker { + namespace helper { + + template + SparseMultiObjectiveModelCheckerHelper::SparseMultiObjectiveModelCheckerHelper(Information& info) : info(info), weightedObjectivesChecker(info) { + overApproximation = storm::storage::geometry::Polytope::createUniversalPolytope(); + underApproximation = storm::storage::geometry::Polytope::createEmptyPolytope(); + } + + template + SparseMultiObjectiveModelCheckerHelper::~SparseMultiObjectiveModelCheckerHelper() { + // Intentionally left empty + } + + template + void SparseMultiObjectiveModelCheckerHelper::achievabilityQuery() { + Point queryPoint; + queryPoint.reserve(info.objectives.size()); + for(auto const& obj : info.objectives) { + STORM_LOG_THROW(obj.threshold, storm::exceptions::UnexpectedException, "Can not perform achievabilityQuery: No threshold given for at least one objective."); + queryPoint.push_back(storm::utility::convertNumber(*obj.threshold)); + } + achievabilityQuery(queryPoint); + } + + template + void SparseMultiObjectiveModelCheckerHelper::achievabilityQuery(Point const& queryPoint) { + storm::storage::BitVector individualObjectivesToBeChecked(info.objectives.size(), true); + while(true) { //TODO introduce convergence criterion? (like max num of iterations) + storm::storage::geometry::Halfspace separatingHalfspace = findSeparatingHalfspace(queryPoint, individualObjectivesToBeChecked); + + weightedObjectivesChecker.check(storm::utility::vector::convertNumericVector(separatingHalfspace.normalVector())); + storm::storage::TotalScheduler scheduler = weightedObjectivesChecker.getScheduler(); + Point weightedObjectivesResult = storm::utility::vector::convertNumericVector(weightedObjectivesChecker.getInitialStateResultOfObjectives()); + + // Insert the computed scheduler and check whether we have already seen it before + auto schedulerInsertRes = schedulers.insert(std::make_pair(std::move(scheduler), paretoOptimalPoints.end())); + if(schedulerInsertRes.second){ + // The scheduler is new, so insert the newly computed pareto optimal point. + // To each scheduler, we assign the (unique) pareto optimal point it induces. + schedulerInsertRes.first->second = paretoOptimalPoints.insert(std::make_pair(std::move(weightedObjectivesResult), std::vector())).first; + } + // In the case where the scheduler is not new, we assume that the corresponding pareto optimal points for the old and new scheduler are equal + // Note that the values might not be exactly equal due to numerical issues. + + Point const& paretoPoint = schedulerInsertRes.first->second->first; + std::vector& weightVectors = schedulerInsertRes.first->second->second; + weightVectors.push_back(std::move(separatingHalfspace.normalVector())); + + updateOverApproximation(paretoPoint, weightVectors.back()); + if(!overApproximation->contains(queryPoint)){ + std::cout << "PROPERTY VIOLATED" << std::endl; + return; + } + updateUnderApproximation(); + if(underApproximation->contains(queryPoint)){ + std::cout << "PROPERTY SATISFIED" << std::endl; + return; + } + } + } + + template + storm::storage::geometry::Halfspace SparseMultiObjectiveModelCheckerHelper::findSeparatingHalfspace(Point const& pointToBeSeparated, storm::storage::BitVector& individualObjectivesToBeChecked) { + // First, we check 'simple' weight vectors that correspond to checking a single objective + while (!individualObjectivesToBeChecked.empty()){ + uint_fast64_t objectiveIndex = individualObjectivesToBeChecked.getNextSetIndex(0); + individualObjectivesToBeChecked.set(objectiveIndex, false); + + WeightVector normalVector; // = (0..0 1 0..0) + normalVector.reserve(info.objectives.size()); + for(uint_fast64_t i = 0; i() : storm::utility::zero()); + } + + storm::storage::geometry::Halfspace h(normalVector, storm::utility::vector::dotProduct(normalVector, pointToBeSeparated)); + bool hIsSeparating = true; + for(auto const& paretoPoint : paretoOptimalPoints){ + hIsSeparating &= h.contains(paretoPoint.first); + } + if(hIsSeparating) { + return h; + } + } + return underApproximation->findSeparatingHalfspace(pointToBeSeparated); + } + + template + void SparseMultiObjectiveModelCheckerHelper::updateOverApproximation(Point const& newPoint, WeightVector const& newWeightVector) { + storm::storage::geometry::Halfspace h(newWeightVector, storm::utility::vector::dotProduct(newWeightVector, newPoint)); + overApproximation = overApproximation->intersection(h); + STORM_LOG_DEBUG("Updated OverApproximation to polytope " << overApproximation->toString(true)); + } + + template + void SparseMultiObjectiveModelCheckerHelper::updateUnderApproximation() { + std::vector paretoPointsVec; + paretoPointsVec.reserve(paretoOptimalPoints.size()); + for(auto const& paretoPoint : paretoOptimalPoints) { + paretoPointsVec.push_back(paretoPoint.first); + } + + boost::optional upperBounds; + if(!paretoPointsVec.empty()){ + //Get the pointwise maximum of the pareto points + upperBounds = paretoPointsVec.front(); + for(auto paretoPointIt = paretoPointsVec.begin()+1; paretoPointIt != paretoPointsVec.end(); ++paretoPointIt){ + auto upperBoundIt = upperBounds->begin(); + for(auto const& paretoPointCoordinate : *paretoPointIt){ + if(paretoPointCoordinate>*upperBoundIt){ + *upperBoundIt = paretoPointCoordinate; + } + ++upperBoundIt; + } + } + } + + underApproximation = storm::storage::geometry::Polytope::create(paretoPointsVec)->downwardClosure(upperBounds); + STORM_LOG_DEBUG("Updated UnderApproximation to polytope " << overApproximation->toString(true)); + } + +#ifdef STORM_HAVE_CARL + template class SparseMultiObjectiveModelCheckerHelper, storm::RationalNumber>; +#endif + } + } +} diff --git a/src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerHelper.h b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerHelper.h new file mode 100644 index 000000000..24d7e43c8 --- /dev/null +++ b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerHelper.h @@ -0,0 +1,64 @@ +#ifndef STORM_MODELCHECKER_MULTIOBJECTIVE_HELPER_SPARSEMULTIOBJECTIVEMODELCHECKERHELPER_H_ +#define STORM_MODELCHECKER_MULTIOBJECTIVE_HELPER_SPARSEMULTIOBJECTIVEMODELCHECKERHELPER_H_ + +#include +#include +#include "src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerInformation.h" +#include "src/modelchecker/multiObjective/helper/SparseWeightedObjectivesModelCheckerHelper.h" +#include "src/storage/geometry/Polytope.h" +#include "src/storage/TotalScheduler.h" + + +namespace storm { + namespace modelchecker { + namespace helper { + + template + class SparseMultiObjectiveModelCheckerHelper { + public: + typedef typename SparseModelType::ValueType ModelValueType; + typedef typename SparseModelType::RewardModelType RewardModelType; + typedef SparseMultiObjectiveModelCheckerInformation Information; + + + typedef std::vector Point; + typedef std::vector WeightVector; + + SparseMultiObjectiveModelCheckerHelper(Information& info); + + ~SparseMultiObjectiveModelCheckerHelper(); + + void achievabilityQuery(); + + + private: + void achievabilityQuery(Point const& queryPoint); + + /* + * Returns a halfspace h that separates the underapproximation from the given point p, i.e., + * - p lies on the border of h and + * - for each x in the underApproximation, it holds that h contains x + * + * @param pointToBeSeparated the point that is to be seperated + * @param objectivesToBeCheckedIndividually stores for each objective whether it still makes sense to check for this objective individually (i.e., with weight vector given by w_{objIndex} = 1 ) + */ + storm::storage::geometry::Halfspace findSeparatingHalfspace(Point const& pointToBeSeparated, storm::storage::BitVector& individualObjectivesToBeChecked); + void updateOverApproximation(Point const& newPoint, WeightVector const& newWeightVector); + void updateUnderApproximation(); + + Information& info; + SparseWeightedObjectivesModelCheckerHelper weightedObjectivesChecker; + + //TODO: sort points as needed + std::map> paretoOptimalPoints; + std::unordered_map>::iterator> schedulers; + + std::shared_ptr> overApproximation; + std::shared_ptr> underApproximation; + }; + + } + } +} + +#endif /* STORM_MODELCHECKER_MULTIOBJECTIVE_HELPER_SPARSEMULTIOBJECTIVEMODELCHECKERHELPER_H_ */ diff --git a/src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerInformation.h b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerInformation.h index 3ef2a3a6d..c44f3c972 100644 --- a/src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerInformation.h +++ b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerInformation.h @@ -49,15 +49,12 @@ namespace storm { }; - SparseMultiObjectiveModelCheckerInformation(SparseModelType const& model) : model(model) { + SparseMultiObjectiveModelCheckerInformation(SparseModelType const& model) : preprocessedModel(model), originalModel(model) { //Intentionally left empty } - - SparseMultiObjectiveModelCheckerInformation(SparseModelType && model) : model(model) { - //Intentionally left empty - } - - SparseModelType model; + + SparseModelType preprocessedModel; + SparseModelType const& originalModel; std::vector newToOldStateIndexMapping; bool negatedRewardsConsidered; @@ -75,8 +72,11 @@ namespace storm { } out << "--------------------------------------------------------------" << std::endl; out << std::endl; + out << "Original Model Information:" << std::endl; + originalModel.printModelInformationToStream(out); + out << std::endl; out << "Preprocessed Model Information:" << std::endl; - model.printModelInformationToStream(out); + preprocessedModel.printModelInformationToStream(out); out << std::endl; if(negatedRewardsConsidered){ out << "The rewards in the preprocessed model are negated" << std::endl; diff --git a/src/modelchecker/multiobjective/helper/SparseWeightedObjectivesModelCheckerHelper.cpp b/src/modelchecker/multiobjective/helper/SparseWeightedObjectivesModelCheckerHelper.cpp new file mode 100644 index 000000000..7015a12e2 --- /dev/null +++ b/src/modelchecker/multiobjective/helper/SparseWeightedObjectivesModelCheckerHelper.cpp @@ -0,0 +1,163 @@ +#include "src/modelchecker/multiobjective/helper/SparseWeightedObjectivesModelCheckerHelper.h" + +#include "src/models/sparse/Mdp.h" +#include "src/models/sparse/StandardRewardModel.h" +#include "src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.h" +#include "src/solver/LinearEquationSolver.h" +#include "src/solver/MinMaxLinearEquationSolver.h" +#include "src/utility/graph.h" +#include "src/utility/macros.h" +#include "src/utility/solver.h" +#include "src/utility/vector.h" +#include "src/exceptions/IllegalFunctionCallException.h" + +namespace storm { + namespace modelchecker { + namespace helper { + + + template + SparseWeightedObjectivesModelCheckerHelper::SparseWeightedObjectivesModelCheckerHelper(Information const& info) : info(info), checkHasBeenCalled(false) , objectiveResults(info.objectives.size()){ + + + } + + template + void SparseWeightedObjectivesModelCheckerHelper::check(std::vector const& weightVector) { + + unboundedWeightedPhase(weightVector); + unboundedIndividualPhase(weightVector); + boundedPhase(weightVector); + checkHasBeenCalled=true; + } + + template + std::vector::ValueType> SparseWeightedObjectivesModelCheckerHelper::getInitialStateResultOfObjectives() const { + STORM_LOG_THROW(checkHasBeenCalled, storm::exceptions::IllegalFunctionCallException, "Tried to retrieve results but check(..) has not been called before."); + STORM_LOG_ASSERT(info.preprocessedModel.getInitialStates().getNumberOfSetBits()==1, "The considered model has multiple initial states"); + std::vector res; + res.reserve(objectiveResults.size()); + for(auto const& objResult : objectiveResults) { + res.push_back(objResult[*info.preprocessedModel.getInitialStates().begin()]); + } + return res; + } + + template + storm::storage::TotalScheduler const& SparseWeightedObjectivesModelCheckerHelper::getScheduler() const { + STORM_LOG_THROW(checkHasBeenCalled, storm::exceptions::IllegalFunctionCallException, "Tried to retrieve results but check(..) has not been called before."); + return scheduler; + } + + template + void SparseWeightedObjectivesModelCheckerHelper::unboundedWeightedPhase(std::vector const& weightVector) { + std::vector weightedRewardVector(info.preprocessedModel.getTransitionMatrix().getRowCount(), storm::utility::zero()); + for(uint_fast64_t objIndex = 0; objIndex < weightVector.size(); ++objIndex) { + if(!info.objectives[objIndex].stepBound){ + storm::utility::vector::addScaledVector(weightedRewardVector, info.preprocessedModel.getRewardModel(info.objectives[objIndex].rewardModelName).getStateActionRewardVector(), weightVector[objIndex]); + } + } + + // TODO check for +/- infty reward... + + + //Testing.. + storm::storage::BitVector test(64); + test.set(63); + std::cout << "Test set index with 0: " << test.getNextSetIndex(0) << std::endl; + std::cout << "Test set index with 63: " << test.getNextSetIndex(63) << std::endl; + std::cout << "Test set index with 64: " << test.getNextSetIndex(64) << std::endl; + + storm::storage::BitVector actionsWithRewards = storm::utility::vector::filter(weightedRewardVector, [&] (ValueType const& value) -> bool {return !storm::utility::isZero(value);}); + storm::storage::BitVector statesWithRewards(info.preprocessedModel.getNumberOfStates(), false); + uint_fast64_t currActionIndex = actionsWithRewards.getNextSetIndex(0); + auto endOfRowGroup = info.preprocessedModel.getTransitionMatrix().getRowGroupIndices().begin() + 1; + for(uint_fast64_t state = 0; stateweightedResult = std::vector(info.preprocessedModel.getNumberOfStates()); + this->scheduler = storm::storage::TotalScheduler(info.preprocessedModel.getNumberOfStates()); + + storm::storage::SparseMatrix submatrix = info.preprocessedModel.getTransitionMatrix().getSubmatrix(true, maybeStates, maybeStates, false); + std::vector b(submatrix.getRowCount()); + storm::utility::vector::selectVectorValues(b, maybeStates, weightedRewardVector); + std::vector x(submatrix.getRowGroupCount(), storm::utility::zero()); + + storm::utility::solver::MinMaxLinearEquationSolverFactory solverFactory; + std::unique_ptr> solver = solverFactory.create(submatrix, true); + solver->solveEquationSystem(x, b); + + storm::utility::vector::setVectorValues(this->weightedResult, maybeStates, x); + storm::utility::vector::setVectorValues(this->weightedResult, ~maybeStates, storm::utility::zero()); + + uint_fast64_t currentSubState = 0; + for (auto maybeState : maybeStates) { + this->scheduler.setChoice(maybeState, solver->getScheduler().getChoice(currentSubState)); + ++currentSubState; + } + // Note that the choices for the ~maybeStates are arbitrary as no states with rewards are reachable any way. + + } + + template + void SparseWeightedObjectivesModelCheckerHelper::unboundedIndividualPhase(std::vector const& weightVector) { + + storm::storage::SparseMatrix deterministicMatrix = info.preprocessedModel.getTransitionMatrix().selectRowsFromRowGroups(this->scheduler.getChoices(), true); + storm::storage::SparseMatrix deterministicBackwardTransitions = deterministicMatrix.transpose(); + std::vector deterministicStateRewards(deterministicMatrix.getRowCount()); + storm::utility::solver::LinearEquationSolverFactory linearEquationSolverFactory; + //TODO check if all but one entry of weightVector is zero + for(uint_fast64_t objIndex = 0; objIndex < weightVector.size(); ++objIndex) { + if(!info.objectives[objIndex].stepBound){ + + storm::utility::vector::selectVectorValues(deterministicStateRewards, this->scheduler.getChoices(), info.preprocessedModel.getTransitionMatrix().getRowGroupIndices(), info.preprocessedModel.getRewardModel(info.objectives[objIndex].rewardModelName).getStateActionRewardVector()); + + storm::storage::BitVector statesWithRewards = storm::utility::vector::filter(deterministicStateRewards, [&] (ValueType const& value) -> bool {return !storm::utility::isZero(value);}); + // As target states, we take the states from which no reward is reachable. + STORM_LOG_WARN("TODO: target state selection is currently only valid for reachability properties..."); + //TODO: we should be able to give some hint to the solver.. + storm::storage::BitVector targetStates = storm::utility::graph::performProbGreater0(deterministicBackwardTransitions, storm::storage::BitVector(deterministicMatrix.getRowCount(), true), statesWithRewards); + objectiveResults[objIndex] = storm::modelchecker::helper::SparseDtmcPrctlHelper::computeReachabilityRewards(deterministicMatrix, + deterministicBackwardTransitions, + deterministicStateRewards, + targetStates, + false, //no qualitative checking, + linearEquationSolverFactory); + } + } + } + + template + void SparseWeightedObjectivesModelCheckerHelper::boundedPhase(std::vector const& weightVector) { + STORM_LOG_WARN("bounded properties not yet implemented"); + } + + + + template class SparseWeightedObjectivesModelCheckerHelper>; + + + } + } +} diff --git a/src/modelchecker/multiobjective/helper/SparseWeightedObjectivesModelCheckerHelper.h b/src/modelchecker/multiobjective/helper/SparseWeightedObjectivesModelCheckerHelper.h new file mode 100644 index 000000000..94f7b9e76 --- /dev/null +++ b/src/modelchecker/multiobjective/helper/SparseWeightedObjectivesModelCheckerHelper.h @@ -0,0 +1,90 @@ +#ifndef STORM_MODELCHECKER_MULTIOBJECTIVE_HELPER_SPARSEWEIGHTEDOBJECTIVESMODELCHECKERHELPER_H_ +#define STORM_MODELCHECKER_MULTIOBJECTIVE_HELPER_SPARSEWEIGHTEDOBJECTIVESMODELCHECKERHELPER_H_ + +#include + +#include "src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerInformation.h" +#include "src/storage/TotalScheduler.h" + +namespace storm { + namespace modelchecker { + namespace helper { + + /*! + * Helper Class that takes a MultiObjectiveInformation 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 SparseWeightedObjectivesModelCheckerHelper { + public: + typedef typename SparseModelType::ValueType ValueType; + typedef typename SparseModelType::RewardModelType RewardModelType; + typedef SparseMultiObjectiveModelCheckerInformation Information; + + SparseWeightedObjectivesModelCheckerHelper(Information const& info); + + /*! + * - 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); + + /*! + * Getter methods for the results of the most recent call of check(..) + * Note that check(..) has to be called before retrieving results. Otherwise, an exception is thrown. + */ + // The results of the individual objectives at the initial state of the given model + std::vector getInitialStateResultOfObjectives() const; + // A scheduler that induces the optimal values + storm::storage::TotalScheduler const& getScheduler() const; + + + private: + + /*! + * Determines the scheduler that maximizes the weighted reward vector of the unbounded objectives + * + * @param weightVector the weight vector of the current check + */ + void unboundedWeightedPhase(std::vector const& weightVector); + + /*! + * 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 + */ + void boundedPhase(std::vector const& weightVector); + + // stores the considered information of the multi-objective model checking problem + Information const& info; + + // 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 results for the individual objectives (for all states of the model) + std::vector> objectiveResults; + // The scheduler that maximizes the weighted rewards + storm::storage::TotalScheduler scheduler; + + }; + + } + } +} + +#endif /* STORM_MODELCHECKER_MULTIOBJECTIVE_HELPER_SPARSEWEIGHTEDOBJECTIVESMODELCHECKERHELPER_H_ */ diff --git a/src/storage/TotalScheduler.cpp b/src/storage/TotalScheduler.cpp index 6ae99d212..d689ae050 100644 --- a/src/storage/TotalScheduler.cpp +++ b/src/storage/TotalScheduler.cpp @@ -1,5 +1,6 @@ #include "src/storage/TotalScheduler.h" #include "src/exceptions/InvalidArgumentException.h" +#include "src/utility/Hash.h" namespace storm { namespace storage { @@ -15,6 +16,10 @@ namespace storm { // Intentionally left empty. } + bool TotalScheduler::operator==(storm::storage::TotalScheduler const& other) const { + return this->choices == other.choices; + } + void TotalScheduler::setChoice(uint_fast64_t state, uint_fast64_t choice) { if (state > choices.size()) { throw storm::exceptions::InvalidArgumentException() << "Invalid call to TotalScheduler::setChoice: scheduler cannot not define a choice for state " << state << "."; @@ -34,6 +39,10 @@ namespace storm { return choices[state]; } + std::vector const& TotalScheduler::getChoices() const { + return choices; + } + std::ostream& operator<<(std::ostream& out, TotalScheduler const& scheduler) { out << "total scheduler (defined on " << scheduler.choices.size() << " states) [ "; for (uint_fast64_t state = 0; state < scheduler.choices.size() - 1; ++state) { @@ -44,6 +53,11 @@ namespace storm { } return out; } - } -} \ No newline at end of file +} + +namespace std { + std::size_t hash::operator()(storm::storage::TotalScheduler const& totalScheduler) const { + return storm::utility::Hash::getHash(totalScheduler.choices); + } +} diff --git a/src/storage/TotalScheduler.h b/src/storage/TotalScheduler.h index 6f95edf71..f44c91205 100644 --- a/src/storage/TotalScheduler.h +++ b/src/storage/TotalScheduler.h @@ -11,6 +11,12 @@ namespace storm { class TotalScheduler : public Scheduler { public: + + TotalScheduler(TotalScheduler const& other) = default; + TotalScheduler(TotalScheduler&& other) = default; + TotalScheduler& operator=(TotalScheduler const& other) = default; + TotalScheduler& operator=(TotalScheduler&& other) = default; + /*! * Creates a total scheduler that defines a choice for the given number of states. By default, all choices * are initialized to zero. @@ -32,6 +38,8 @@ namespace storm { * @param choices A vector whose i-th entry defines the choice of state i. */ TotalScheduler(std::vector&& choices); + + bool operator==(TotalScheduler const& other) const; void setChoice(uint_fast64_t state, uint_fast64_t choice) override; @@ -39,7 +47,10 @@ namespace storm { uint_fast64_t getChoice(uint_fast64_t state) const override; + std::vector const& getChoices() const; + friend std::ostream& operator<<(std::ostream& out, TotalScheduler const& scheduler); + friend struct std::hash; private: // A vector that stores the choice for each state. @@ -48,4 +59,12 @@ namespace storm { } // namespace storage } // namespace storm + +namespace std { + template <> + struct hash { + std::size_t operator()(storm::storage::TotalScheduler const& totalScheduler) const; + }; +} + #endif /* STORM_STORAGE_TOTALSCHEDULER_H_ */ diff --git a/src/storage/geometry/Polytope.cpp b/src/storage/geometry/Polytope.cpp index 791eb6d86..0b7c75d69 100644 --- a/src/storage/geometry/Polytope.cpp +++ b/src/storage/geometry/Polytope.cpp @@ -22,6 +22,16 @@ namespace storm { return create(boost::none, points); } + template + std::shared_ptr> Polytope::createUniversalPolytope() { + return create(std::vector>(), boost::none); + } + + template + std::shared_ptr> Polytope::createEmptyPolytope() { + return create(boost::none, std::vector()); + } + #ifdef STORM_HAVE_CARL template <> std::shared_ptr> Polytope::create(boost::optional>> const& halfspaces, @@ -118,6 +128,12 @@ namespace storm { return nullptr; } + template + Halfspace Polytope::findSeparatingHalfspace(Point const& pointToBeSeparated) const { + STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "Functionality not implemented."); + return Halfspace(Point(), storm::utility::zero()); + } + template std::string Polytope::toString(bool numbersAsDouble) const { auto halfspaces = this->getHalfspaces(); diff --git a/src/storage/geometry/Polytope.h b/src/storage/geometry/Polytope.h index be113c934..e95c894c2 100644 --- a/src/storage/geometry/Polytope.h +++ b/src/storage/geometry/Polytope.h @@ -31,6 +31,16 @@ namespace storm { */ static std::shared_ptr> create(std::vector const& points); + /*! + * Creates the universal polytope (i.e., the set R^n) + */ + static std::shared_ptr> createUniversalPolytope(); + + /*! + * Creates the empty polytope (i.e., emptyset) + */ + static std::shared_ptr> createEmptyPolytope(); + /*! * Returns the vertices of this polytope. */ @@ -88,6 +98,17 @@ namespace storm { */ virtual std::shared_ptr> downwardClosure(boost::optional const& upperBounds = boost::none) const; + /*! + * Returns a halfspace h that separates this polytope from the given point p, i.e., + * - p lies on the border of h and + * - for each x in this polytope, it holds that h contains x + * + * A halfspace with maximal distance to the polytope is preferred. + * + * @param pointToBeSeparated the point that is to be seperated + */ + virtual Halfspace findSeparatingHalfspace(Point const& pointToBeSeparated) const; + /* * Returns a string representation of this polytope. * If the given flag is true, the occurring numbers are converted to double before printing to increase readability diff --git a/src/utility/vector.h b/src/utility/vector.h index 8fc12959c..37500feb7 100644 --- a/src/utility/vector.h +++ b/src/utility/vector.h @@ -339,14 +339,26 @@ namespace storm { /*! * Multiplies each element of the given vector with the given factor and writes the result into the vector. * - * @param target The first summand and target vector. - * @param summand The second summand. + * @param target The operand and target vector. + * @param factor The scaling factor */ template void scaleVectorInPlace(std::vector& target, ValueType2 const& factor) { applyPointwise(target, target, [&] (ValueType1 const& argument) -> ValueType1 { return argument * factor; }); } + /*! + * Adds each element of the first vector and (the corresponding element of the second vector times the given factor) and writes the result into the first vector. + * + * @param firstOperand The first operand. + * @param secondOperand The second operand + * @param factor The factor for the elements of the second operand + */ + template + void addScaledVector(std::vector& firstOperand, std::vector const& secondOperand, InValueType3 const& factor) { + applyPointwise(firstOperand, secondOperand, firstOperand, [&] (InValueType1 const& val1, InValueType2 const& val2) -> InValueType1 { return val1 + (factor * val2); }); + } + /*! * Computes the dot product (aka scalar product) and returns the result * @@ -692,17 +704,30 @@ namespace storm { return subVector; } + /*! + * Converts the given vector to the given ValueType + * Assumes that both, TargetType and SourceType are numeric + */ + template + std::vector convertNumericVector(std::vector const& oldVector) { + std::vector resultVector; + resultVector.reserve(oldVector.size()); + for(auto const& oldValue : oldVector){ + resultVector.push_back(storm::utility::convertNumber(oldValue)); + } + return resultVector; + } + /*! * Converts the given vector to the given ValueType */ template - std::vector toValueType(std::vector const& oldVector) { + typename std::vector toValueType(std::vector const& oldVector) { std::vector resultVector; - resultVector.resize(oldVector.size()); - for (size_t i = 0, size = oldVector.size(); i < size; ++i) { - resultVector.at(i) = static_cast(oldVector.at(i)); - } - + resultVector.reserve(oldVector.size()); + for(auto const& oldValue : oldVector){ + resultVector.push_back(static_cast(oldValue)); + } return resultVector; }