TimQu
9 years ago
22 changed files with 1997 additions and 207 deletions
-
60src/modelchecker/multiobjective/Pcaa.cpp
-
20src/modelchecker/multiobjective/Pcaa.h
-
30src/modelchecker/multiobjective/SparseMaMultiObjectiveModelChecker.cpp
-
31src/modelchecker/multiobjective/SparseMdpMultiObjectiveModelChecker.cpp
-
11src/modelchecker/multiobjective/helper/SparseMaMultiObjectiveWeightVectorChecker.cpp
-
2src/modelchecker/multiobjective/helper/SparseMultiObjectiveWeightVectorChecker.cpp
-
4src/modelchecker/multiobjective/pcaa/PCAAObjective.h
-
423src/modelchecker/multiobjective/pcaa/SparseMaPcaaWeightVectorChecker.cpp
-
157src/modelchecker/multiobjective/pcaa/SparseMaPcaaWeightVectorChecker.h
-
122src/modelchecker/multiobjective/pcaa/SparseMdpPcaaWeightVectorChecker.cpp
-
48src/modelchecker/multiobjective/pcaa/SparseMdpPcaaWeightVectorChecker.h
-
130src/modelchecker/multiobjective/pcaa/SparsePCAAPreprocessor.cpp
-
42src/modelchecker/multiobjective/pcaa/SparsePCAAPreprocessor.h
-
33src/modelchecker/multiobjective/pcaa/SparsePCAAPreprocessorReturnType.h
-
169src/modelchecker/multiobjective/pcaa/SparsePcaaQuantitativeQuery.cpp
-
65src/modelchecker/multiobjective/pcaa/SparsePcaaQuantitativeQuery.h
-
205src/modelchecker/multiobjective/pcaa/SparsePcaaQuery.cpp
-
122src/modelchecker/multiobjective/pcaa/SparsePcaaQuery.h
-
311src/modelchecker/multiobjective/pcaa/SparsePcaaWeightVectorChecker.cpp
-
155src/modelchecker/multiobjective/pcaa/SparsePcaaWeightVectorChecker.h
-
60src/transformer/EndComponentEliminator.h
-
4test/functional/transformer/EndComponentEliminatorTest.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<typename SparseModelType> |
||||
|
std::unique_ptr<CheckResult> 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<SparseModelType>::preprocess(model, formula); |
||||
|
STORM_LOG_DEBUG("Preprocessing done. Result: " << preprocessorResult); |
||||
|
|
||||
|
std::unique_ptr<SparsePcaaQuery<SparseModelType, storm::RationalNumber>> query; |
||||
|
switch (preprocessorResult.queryType) { |
||||
|
case SparsePcaaPreprocessorReturnType<SparseModelType>::QueryType::Achievability: |
||||
|
query = std::unique_ptr<SparsePcaaQuery<SparseModelType, storm::RationalNumber>> (new SparsePcaaAchievabilityQuery<SparseModelType, storm::RationalNumber>(preprocessorResult)); |
||||
|
break; |
||||
|
case SparsePcaaPreprocessorReturnType<SparseModelType>::QueryType::Quantitative: |
||||
|
query = std::unique_ptr<SparsePcaaQuery<SparseModelType, storm::RationalNumber>> (new SparsePcaaQuantitativeQuery<SparseModelType, storm::RationalNumber>(preprocessorResult)); |
||||
|
break; |
||||
|
case SparsePcaaPreprocessorReturnType<SparseModelType>::QueryType::Pareto: |
||||
|
query = std::unique_ptr<SparsePcaaQuery<SparseModelType, storm::RationalNumber>> (new SparsePcaaParetoQuery<SparseModelType, storm::RationalNumber>(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<CheckResult> performPcaa<storm::models::sparse::Mdp<double>>(storm::models::sparse::Mdp<double> const& model, storm::logic::MultiObjectiveFormula const& formula); |
||||
|
template std::unique_ptr<CheckResult> performPcaa<storm::models::sparse::MarkovAutomaton<double>>(storm::models::sparse::MarkovAutomaton<double> const& model, storm::logic::MultiObjectiveFormula const& formula); |
||||
|
#ifdef STORM_HAVE_CARL
|
||||
|
template std::unique_ptr<CheckResult> performPcaa<storm::models::sparse::Mdp<storm::RationalNumber>>(storm::models::sparse::Mdp<storm::RationalNumber> const& model, storm::logic::MultiObjectiveFormula const& formula); |
||||
|
// template std::unique_ptr<CheckResult> performPcaa<storm::models::sparse::MarkovAutomaton<storm::RationalNumber>>(storm::models::sparse::MarkovAutomaton<storm::RationalNumber> const& model, storm::logic::MultiObjectiveFormula const& formula);
|
||||
|
#endif
|
||||
|
|
||||
|
} |
||||
|
} |
||||
|
} |
||||
@ -0,0 +1,20 @@ |
|||||
|
#ifndef STORM_MODELCHECKER_MULTIOBJECTIVE_PCAA_H_ |
||||
|
#define STORM_MODELCHECKER_MULTIOBJECTIVE_PCAA_H_ |
||||
|
|
||||
|
#include <memory> |
||||
|
|
||||
|
#include "src/modelchecker/results/CheckResult.h" |
||||
|
#include "src/logic/Formulas.h" |
||||
|
|
||||
|
namespace storm { |
||||
|
namespace modelchecker { |
||||
|
namespace multiobjective { |
||||
|
|
||||
|
template<typename SparseModelType> |
||||
|
std::unique_ptr<CheckResult> performPcaa(SparseModelType const& model, storm::logic::MultiObjectiveFormula const& formula); |
||||
|
|
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
#endif /* STORM_MODELCHECKER_MULTIOBJECTIVE_PCAA_H_ */ |
||||
@ -0,0 +1,423 @@ |
|||||
|
#include "src/modelchecker/multiobjective/pcaa/SparseMaPcaaWeightVectorChecker.h"
|
||||
|
|
||||
|
#include <cmath>
|
||||
|
|
||||
|
#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 <class SparseMaModelType> |
||||
|
SparseMaPcaaWeightVectorChecker<SparseMaModelType>::SparseMaPcaaWeightVectorChecker(SparseMaModelType const& model, |
||||
|
std::vector<PcaaObjective<ValueType>> const& objectives, |
||||
|
storm::storage::BitVector const& actionsWithNegativeReward, |
||||
|
storm::storage::BitVector const& ecActions, |
||||
|
storm::storage::BitVector const& possiblyRecurrentStates) : |
||||
|
SparsePcaaWeightVectorChecker<SparseMaModelType>(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<ValueType>(this->model.getTransitionMatrix().getRowCount(), storm::utility::zero<ValueType>()); |
||||
|
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 <class SparseMaModelType> |
||||
|
void SparseMaPcaaWeightVectorChecker<SparseMaModelType>::boundedPhase(std::vector<ValueType> const& weightVector, std::vector<ValueType>& 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<MinMaxSolverData> 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<LinEqSolverData> linEq = initLinEqSolver(PS); |
||||
|
|
||||
|
// Store the optimal choices of PS as computed by the minMax solver.
|
||||
|
std::vector<uint_fast64_t> optimalChoicesAtCurrentEpoch(PS.getNumberOfStates(), std::numeric_limits<uint_fast64_t>::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 <class SparseMaModelType> |
||||
|
typename SparseMaPcaaWeightVectorChecker<SparseMaModelType>::SubModel SparseMaPcaaWeightVectorChecker<SparseMaModelType>::createSubModel(bool createMS, std::vector<ValueType> 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<ValueType>& objVector = result.objectiveRewardVectors[objIndex]; |
||||
|
objVector = std::vector<ValueType>(result.weightedRewardVector.size(), storm::utility::zero<ValueType>()); |
||||
|
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 <class SparseMaModelType> |
||||
|
template <typename VT, typename std::enable_if<storm::NumberTraits<VT>::SupportsExponential, int>::type> |
||||
|
VT SparseMaPcaaWeightVectorChecker<SparseMaModelType>::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<std::pair<VT, VT>> eToPowerOfMinusMaxRateTimesBound; |
||||
|
VT smallestNonZeroBound = storm::utility::zero<VT>(); |
||||
|
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<VT>(); |
||||
|
} |
||||
|
|
||||
|
// 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<VT>(); |
||||
|
if(obj.lowerTimeBound) { |
||||
|
precisionOfObj += storm::utility::one<VT>() - (eToPowerOfMinusMaxRateTimesBound[objIndex].first * storm::utility::pow(storm::utility::one<VT>() + maxRate * delta, *obj.lowerTimeBound / delta) ) |
||||
|
+ storm::utility::one<VT>() - std::exp(-maxRate * delta); |
||||
|
} |
||||
|
if(obj.upperTimeBound) { |
||||
|
precisionOfObj += storm::utility::one<VT>() - (eToPowerOfMinusMaxRateTimesBound[objIndex].second * storm::utility::pow(storm::utility::one<VT>() + 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 <class SparseMaModelType> |
||||
|
template <typename VT, typename std::enable_if<!storm::NumberTraits<VT>::SupportsExponential, int>::type> |
||||
|
VT SparseMaPcaaWeightVectorChecker<SparseMaModelType>::getDigitizationConstant() const { |
||||
|
STORM_LOG_THROW(false, storm::exceptions::InvalidOperationException, "Computing bounded probabilities of MAs is unsupported for this value type."); |
||||
|
} |
||||
|
|
||||
|
template <class SparseMaModelType> |
||||
|
template <typename VT, typename std::enable_if<storm::NumberTraits<VT>::SupportsExponential, int>::type> |
||||
|
void SparseMaPcaaWeightVectorChecker<SparseMaModelType>::digitize(SubModel& MS, VT const& digitizationConstant) const { |
||||
|
std::vector<VT> 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<VT>() - eToMinusRateTimesDelta) * entry.getValue()); |
||||
|
if(entry.getColumn() == row) { |
||||
|
entry.setValue(entry.getValue() + eToMinusRateTimesDelta); |
||||
|
} |
||||
|
} |
||||
|
for(auto& entry : MS.toPS.getRow(row)) { |
||||
|
entry.setValue((storm::utility::one<VT>() - eToMinusRateTimesDelta) * entry.getValue()); |
||||
|
} |
||||
|
MS.weightedRewardVector[row] *= storm::utility::one<VT>() - eToMinusRateTimesDelta; |
||||
|
for(auto& objVector : MS.objectiveRewardVectors) { |
||||
|
objVector[row] *= storm::utility::one<VT>() - eToMinusRateTimesDelta; |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
template <class SparseMaModelType> |
||||
|
template <typename VT, typename std::enable_if<!storm::NumberTraits<VT>::SupportsExponential, int>::type> |
||||
|
void SparseMaPcaaWeightVectorChecker<SparseMaModelType>::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 <class SparseMaModelType> |
||||
|
template <typename VT, typename std::enable_if<storm::NumberTraits<VT>::SupportsExponential, int>::type> |
||||
|
void SparseMaPcaaWeightVectorChecker<SparseMaModelType>::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<uint_fast64_t>((*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<VT>(); |
||||
|
digitizationError -= std::exp(-maxRate * (*obj.lowerTimeBound)) * storm::utility::pow(storm::utility::one<VT>() + maxRate * digitizationConstant, digitizedBound); |
||||
|
errorTowardsZero = -digitizationError; |
||||
|
errorAwayFromZero = storm::utility::one<VT>() - std::exp(-maxRate * digitizationConstant);; |
||||
|
} else { |
||||
|
errorTowardsZero = storm::utility::zero<VT>(); |
||||
|
errorAwayFromZero = storm::utility::zero<VT>(); |
||||
|
} |
||||
|
if(obj.upperTimeBound) { |
||||
|
uint_fast64_t digitizedBound = storm::utility::convertNumber<uint_fast64_t>((*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<VT>(); |
||||
|
digitizationError -= std::exp(-maxRate * (*obj.upperTimeBound)) * storm::utility::pow(storm::utility::one<VT>() + 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 <class SparseMaModelType> |
||||
|
template <typename VT, typename std::enable_if<!storm::NumberTraits<VT>::SupportsExponential, int>::type> |
||||
|
void SparseMaPcaaWeightVectorChecker<SparseMaModelType>::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 <class SparseMaModelType> |
||||
|
std::unique_ptr<typename SparseMaPcaaWeightVectorChecker<SparseMaModelType>::MinMaxSolverData> SparseMaPcaaWeightVectorChecker<SparseMaModelType>::initMinMaxSolver(SubModel const& PS) const { |
||||
|
std::unique_ptr<MinMaxSolverData> result(new MinMaxSolverData()); |
||||
|
storm::solver::GeneralMinMaxLinearEquationSolverFactory<ValueType> 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 <class SparseMaModelType> |
||||
|
std::unique_ptr<typename SparseMaPcaaWeightVectorChecker<SparseMaModelType>::LinEqSolverData> SparseMaPcaaWeightVectorChecker<SparseMaModelType>::initLinEqSolver(SubModel const& PS) const { |
||||
|
std::unique_ptr<LinEqSolverData> 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<ValueType>::SolutionMethod::Jacobi); |
||||
|
result->factory.getSettings().setPreconditioner(storm::solver::GmmxxLinearEquationSolverSettings<ValueType>::Preconditioner::None); |
||||
|
result->b.resize(PS.getNumberOfStates()); |
||||
|
return result; |
||||
|
} |
||||
|
|
||||
|
template <class SparseMaModelType> |
||||
|
void SparseMaPcaaWeightVectorChecker<SparseMaModelType>::updateDataToCurrentEpoch(SubModel& MS, SubModel& PS, MinMaxSolverData& minMax, storm::storage::BitVector& consideredObjectives, uint_fast64_t const& currentEpoch, std::vector<ValueType> 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<ValueType>(MS.objectiveRewardVectors[objIndex].size(), storm::utility::zero<ValueType>()); |
||||
|
PS.objectiveRewardVectors[objIndex] = std::vector<ValueType>(PS.objectiveRewardVectors[objIndex].size(), storm::utility::zero<ValueType>()); |
||||
|
} |
||||
|
++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 <class SparseMaModelType> |
||||
|
void SparseMaPcaaWeightVectorChecker<SparseMaModelType>::performPSStep(SubModel& PS, SubModel const& MS, MinMaxSolverData& minMax, LinEqSolverData& linEq, std::vector<uint_fast64_t>& 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<ValueType> 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 <class SparseMaModelType> |
||||
|
void SparseMaPcaaWeightVectorChecker<SparseMaModelType>::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<storm::models::sparse::MarkovAutomaton<double>>; |
||||
|
template double SparseMaPcaaWeightVectorChecker<storm::models::sparse::MarkovAutomaton<double>>::getDigitizationConstant<double>() const; |
||||
|
template void SparseMaPcaaWeightVectorChecker<storm::models::sparse::MarkovAutomaton<double>>::digitize<double>(SparseMaPcaaWeightVectorChecker<storm::models::sparse::MarkovAutomaton<double>>::SubModel& subModel, double const& digitizationConstant) const; |
||||
|
template void SparseMaPcaaWeightVectorChecker<storm::models::sparse::MarkovAutomaton<double>>::digitizeTimeBounds<double>(SparseMaPcaaWeightVectorChecker<storm::models::sparse::MarkovAutomaton<double>>::TimeBoundMap& lowerTimeBounds, SparseMaPcaaWeightVectorChecker<storm::models::sparse::MarkovAutomaton<double>>::TimeBoundMap& upperTimeBounds, double const& digitizationConstant); |
||||
|
#ifdef STORM_HAVE_CARL
|
||||
|
// template class SparseMaPcaaWeightVectorChecker<storm::models::sparse::MarkovAutomaton<storm::RationalNumber>>;
|
||||
|
// template storm::RationalNumber SparseMaPcaaWeightVectorChecker<storm::models::sparse::MarkovAutomaton<storm::RationalNumber>>::getDigitizationConstant<storm::RationalNumber>() const;
|
||||
|
// template void SparseMaPcaaWeightVectorChecker<storm::models::sparse::MarkovAutomaton<storm::RationalNumber>>::digitize<storm::RationalNumber>(SparseMaPcaaWeightVectorChecker<storm::models::sparse::MarkovAutomaton<storm::RationalNumber>>::SubModel& subModel, storm::RationalNumber const& digitizationConstant) const;
|
||||
|
// template void SparseMaPcaaWeightVectorChecker<storm::models::sparse::MarkovAutomaton<storm::RationalNumber>>::digitizeTimeBounds<storm::RationalNumber>(SparseMaPcaaWeightVectorChecker<storm::models::sparse::MarkovAutomaton<double>>::TimeBoundMap& lowerTimeBounds, SparseMaPcaaWeightVectorChecker<storm::models::sparse::MarkovAutomaton<double>>::TimeBoundMap& upperTimeBounds, storm::RationalNumber const& digitizationConstant);
|
||||
|
#endif
|
||||
|
|
||||
|
} |
||||
|
} |
||||
|
} |
||||
@ -0,0 +1,157 @@ |
|||||
|
#ifndef STORM_MODELCHECKER_MULTIOBJECTIVE_PCAA_SPARSEMAPCAAWEIGHTVECTORCHECKER_H_ |
||||
|
#define STORM_MODELCHECKER_MULTIOBJECTIVE_PCAA_SPARSEMAPCAAWEIGHTVECTORCHECKER_H_ |
||||
|
|
||||
|
#include <vector> |
||||
|
#include <type_traits> |
||||
|
|
||||
|
#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 SparseMaModelType> |
||||
|
class SparseMaPcaaWeightVectorChecker : public SparsePcaaWeightVectorChecker<SparseMaModelType> { |
||||
|
public: |
||||
|
typedef typename SparseMaModelType::ValueType ValueType; |
||||
|
|
||||
|
SparseMaPcaaWeightVectorChecker(SparseMaModelType const& model, |
||||
|
std::vector<PcaaObjective<ValueType>> 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<uint_fast64_t, storm::storage::BitVector, std::greater<uint_fast64_t>> 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<ValueType> toMS; // Transitions to Markovian states |
||||
|
storm::storage::SparseMatrix<ValueType> toPS; // Transitions to probabilistic states |
||||
|
|
||||
|
std::vector<ValueType> weightedRewardVector; |
||||
|
std::vector<std::vector<ValueType>> objectiveRewardVectors; |
||||
|
|
||||
|
std::vector<ValueType> weightedSolutionVector; |
||||
|
std::vector<std::vector<ValueType>> objectiveSolutionVectors; |
||||
|
|
||||
|
std::vector<ValueType> 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<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver; |
||||
|
std::vector<ValueType> b; |
||||
|
}; |
||||
|
|
||||
|
struct LinEqSolverData { |
||||
|
storm::solver::GmmxxLinearEquationSolverFactory<ValueType> factory; |
||||
|
std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver; |
||||
|
std::vector<ValueType> 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<ValueType> const& weightVector, std::vector<ValueType>& 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<ValueType> const& weightedRewardVector) const; |
||||
|
|
||||
|
/*! |
||||
|
* Retrieves the delta used for digitization |
||||
|
*/ |
||||
|
template <typename VT = ValueType, typename std::enable_if<storm::NumberTraits<VT>::SupportsExponential, int>::type = 0> |
||||
|
VT getDigitizationConstant() const; |
||||
|
template <typename VT = ValueType, typename std::enable_if<!storm::NumberTraits<VT>::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 <typename VT = ValueType, typename std::enable_if<storm::NumberTraits<VT>::SupportsExponential, int>::type = 0> |
||||
|
void digitize(SubModel& subModel, VT const& digitizationConstant) const; |
||||
|
template <typename VT = ValueType, typename std::enable_if<!storm::NumberTraits<VT>::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 <typename VT = ValueType, typename std::enable_if<storm::NumberTraits<VT>::SupportsExponential, int>::type = 0> |
||||
|
void digitizeTimeBounds(TimeBoundMap& lowerTimeBounds, TimeBoundMap& upperTimeBounds, VT const& digitizationConstant); |
||||
|
template <typename VT = ValueType, typename std::enable_if<!storm::NumberTraits<VT>::SupportsExponential, int>::type = 0> |
||||
|
void digitizeTimeBounds(TimeBoundMap& lowerTimeBounds, TimeBoundMap& upperTimeBounds, VT const& digitizationConstant); |
||||
|
|
||||
|
|
||||
|
/*! |
||||
|
* Initializes the data for the MinMax solver |
||||
|
*/ |
||||
|
std::unique_ptr<MinMaxSolverData> initMinMaxSolver(SubModel const& PS) const; |
||||
|
|
||||
|
/*! |
||||
|
* Initializes the data for the LinEq solver |
||||
|
*/ |
||||
|
std::unique_ptr<LinEqSolverData> 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<ValueType> 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<uint_fast64_t>& 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_ */ |
||||
@ -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 <class SparseMdpModelType> |
||||
|
SparseMdpPcaaWeightVectorChecker<SparseMdpModelType>::SparseMdpPcaaWeightVectorChecker(SparseMdpModelType const& model, |
||||
|
std::vector<PcaaObjective<ValueType>> const& objectives, |
||||
|
storm::storage::BitVector const& actionsWithNegativeReward, |
||||
|
storm::storage::BitVector const& ecActions, |
||||
|
storm::storage::BitVector const& possiblyRecurrentStates) : |
||||
|
SparsePcaaWeightVectorChecker<SparseMdpModelType>(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 <class SparseMdpModelType> |
||||
|
void SparseMdpPcaaWeightVectorChecker<SparseMdpModelType>::boundedPhase(std::vector<ValueType> const& weightVector, std::vector<ValueType>& weightedRewardVector) { |
||||
|
// Allocate some memory so this does not need to happen for each time epoch
|
||||
|
std::vector<uint_fast64_t> optimalChoicesInCurrentEpoch(this->model.getNumberOfStates()); |
||||
|
std::vector<ValueType> choiceValues(weightedRewardVector.size()); |
||||
|
std::vector<ValueType> temporaryResult(this->model.getNumberOfStates()); |
||||
|
std::vector<ValueType> zeroReward(weightedRewardVector.size(), storm::utility::zero<ValueType>()); |
||||
|
// Get for each occurring timeBound the indices of the objectives with that bound.
|
||||
|
std::map<uint_fast64_t, storm::storage::BitVector, std::greater<uint_fast64_t>> lowerTimeBounds; |
||||
|
std::map<uint_fast64_t, storm::storage::BitVector, std::greater<uint_fast64_t>> 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<uint_fast64_t>(*obj.lowerTimeBound), storm::storage::BitVector(this->objectives.size(), false))).first; |
||||
|
STORM_LOG_WARN_COND(storm::utility::convertNumber<ValueType>(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<uint_fast64_t>(*obj.upperTimeBound), storm::storage::BitVector(this->objectives.size(), false))).first; |
||||
|
STORM_LOG_WARN_COND(storm::utility::convertNumber<ValueType>(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<ValueType>(); |
||||
|
this->offsetsToUpperBound[objIndex] = storm::utility::zero<ValueType>(); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// 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<ValueType>& objectiveResult = this->objectiveResults[objIndex]; |
||||
|
std::vector<ValueType> 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<storm::models::sparse::Mdp<double>>; |
||||
|
#ifdef STORM_HAVE_CARL
|
||||
|
template class SparseMdpPcaaWeightVectorChecker<storm::models::sparse::Mdp<storm::RationalNumber>>; |
||||
|
#endif
|
||||
|
|
||||
|
} |
||||
|
} |
||||
|
} |
||||
@ -0,0 +1,48 @@ |
|||||
|
#ifndef STORM_MODELCHECKER_MULTIOBJECTIVE_PCAA_SPARSEMDPPCAAWEIGHTVECTORCHECKER_H_ |
||||
|
#define STORM_MODELCHECKER_MULTIOBJECTIVE_PCAA_SPARSEMDPPCAAWEIGHTVECTORCHECKER_H_ |
||||
|
|
||||
|
#include <vector> |
||||
|
|
||||
|
#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 SparseMdpModelType> |
||||
|
class SparseMdpPcaaWeightVectorChecker : public SparsePcaaWeightVectorChecker<SparseMdpModelType> { |
||||
|
public: |
||||
|
typedef typename SparseMdpModelType::ValueType ValueType; |
||||
|
|
||||
|
SparseMdpPcaaWeightVectorChecker(SparseMdpModelType const& model, |
||||
|
std::vector<PcaaObjective<ValueType>> 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<ValueType> const& weightVector, std::vector<ValueType>& weightedRewardVector) override; |
||||
|
|
||||
|
}; |
||||
|
|
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
#endif /* STORM_MODELCHECKER_MULTIOBJECTIVE_PCAA_SPARSEMDPPCAAWEIGHTVECTORCHECKER_H_ */ |
||||
@ -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 <class SparseModelType, typename GeometryValueType> |
||||
|
SparsePcaaQuantitativeQuery<SparseModelType, GeometryValueType>::SparsePcaaQuantitativeQuery(SparsePcaaPreprocessorReturnType<SparseModelType>& preprocessorResult) : SparsePcaaQuery<SparseModelType, GeometryValueType>(preprocessorResult) { |
||||
|
STORM_LOG_ASSERT(preprocessorResult.queryType==SparsePcaaPreprocessorReturnType<SparseModelType>::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<typename SparseModelType::ValueType>(storm::settings::getModule<storm::settings::modules::MultiObjectiveSettings>().getPrecision()); |
||||
|
gap /= (storm::utility::one<typename SparseModelType::ValueType>() + storm::utility::one<typename SparseModelType::ValueType>()); |
||||
|
gap /= storm::utility::sqrt(static_cast<typename SparseModelType::ValueType>(this->objectives.size())); |
||||
|
this->weightVectorChecker->setMaximumLowerUpperBoundGap(gap); |
||||
|
|
||||
|
} |
||||
|
|
||||
|
template <class SparseModelType, typename GeometryValueType> |
||||
|
void SparsePcaaQuantitativeQuery<SparseModelType, GeometryValueType>::initializeThresholdData() { |
||||
|
thresholds.reserve(this->objectives.size()); |
||||
|
strictThresholds = storm::storage::BitVector(this->objectives.size(), false); |
||||
|
std::vector<storm::storage::geometry::Halfspace<GeometryValueType>> 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<GeometryValueType>(*this->objectives[objIndex].threshold)); |
||||
|
WeightVector normalVector(this->objectives.size(), storm::utility::zero<GeometryValueType>()); |
||||
|
normalVector[objIndex] = -storm::utility::one<GeometryValueType>(); |
||||
|
thresholdConstraints.emplace_back(std::move(normalVector), -thresholds.back()); |
||||
|
strictThresholds.set(objIndex, this->objectives[objIndex].thresholdIsStrict); |
||||
|
} else { |
||||
|
thresholds.push_back(storm::utility::zero<GeometryValueType>()); |
||||
|
} |
||||
|
} |
||||
|
// Note: If we have a single objective (i.e., no objectives with thresholds), thresholdsAsPolytope gets no constraints
|
||||
|
thresholdsAsPolytope = storm::storage::geometry::Polytope<GeometryValueType>::create(thresholdConstraints); |
||||
|
} |
||||
|
|
||||
|
template <class SparseModelType, typename GeometryValueType> |
||||
|
std::unique_ptr<CheckResult> SparsePcaaQuantitativeQuery<SparseModelType, GeometryValueType>::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<typename SparseModelType::ValueType>(result) * |
||||
|
this->objectives[indexOfOptimizingObjective].toOriginalValueTransformationFactor + |
||||
|
this->objectives[indexOfOptimizingObjective].toOriginalValueTransformationOffset; |
||||
|
return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<typename SparseModelType::ValueType>(this->originalModel.getInitialStates().getNextSetIndex(0), resultForOriginalModel)); |
||||
|
} else { |
||||
|
return std::unique_ptr<CheckResult>(new ExplicitQualitativeCheckResult(this->originalModel.getInitialStates().getNextSetIndex(0), false)); |
||||
|
} |
||||
|
|
||||
|
|
||||
|
} |
||||
|
|
||||
|
template <class SparseModelType, typename GeometryValueType> |
||||
|
bool SparsePcaaQuantitativeQuery<SparseModelType, GeometryValueType>::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 <class SparseModelType, typename GeometryValueType> |
||||
|
GeometryValueType SparsePcaaQuantitativeQuery<SparseModelType, GeometryValueType>::improveSolution() { |
||||
|
this->diracWeightVectorsToBeChecked.clear(); // Only check weight vectors that can actually improve the solution
|
||||
|
|
||||
|
WeightVector directionOfOptimizingObjective(this->objectives.size(), storm::utility::zero<GeometryValueType>()); |
||||
|
directionOfOptimizingObjective[indexOfOptimizingObjective] = storm::utility::one<GeometryValueType>(); |
||||
|
|
||||
|
// 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<GeometryValueType>(); |
||||
|
while(!this->maxStepsPerformed()) { |
||||
|
std::pair<Point, bool> 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<double>(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<GeometryValueType>(storm::settings::getModule<storm::settings::modules::MultiObjectiveSettings>().getPrecision())) { |
||||
|
// Goal precision reached!
|
||||
|
return result; |
||||
|
} else { |
||||
|
STORM_LOG_DEBUG("Solution can be improved by at most " << storm::utility::convertNumber<double>(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 <class SparseModelType, typename GeometryValueType> |
||||
|
bool SparsePcaaQuantitativeQuery<SparseModelType, GeometryValueType>::checkIfThresholdsAreSatisfied(std::shared_ptr<storm::storage::geometry::Polytope<GeometryValueType>> const& polytope) { |
||||
|
std::vector<storm::storage::geometry::Halfspace<GeometryValueType>> halfspaces = polytope->getHalfspaces(); |
||||
|
for(auto const& h : halfspaces) { |
||||
|
GeometryValueType distance = h.distance(thresholds); |
||||
|
if(distance < storm::utility::zero<GeometryValueType>()) { |
||||
|
return false; |
||||
|
} |
||||
|
if(distance == storm::utility::zero<GeometryValueType>()) { |
||||
|
// 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<GeometryValueType>()) { |
||||
|
return false; |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
return true; |
||||
|
} |
||||
|
|
||||
|
|
||||
|
|
||||
|
#ifdef STORM_HAVE_CARL
|
||||
|
template class SparsePcaaQuantitativeQuery<storm::models::sparse::Mdp<double>, storm::RationalNumber>; |
||||
|
template class SparsePcaaQuantitativeQuery<storm::models::sparse::MarkovAutomaton<double>, storm::RationalNumber>; |
||||
|
|
||||
|
template class SparsePcaaQuantitativeQuery<storm::models::sparse::Mdp<storm::RationalNumber>, storm::RationalNumber>; |
||||
|
// template class SparsePcaaQuantitativeQuery<storm::models::sparse::MarkovAutomaton<storm::RationalNumber>, storm::RationalNumber>;
|
||||
|
#endif
|
||||
|
} |
||||
|
} |
||||
|
} |
||||
@ -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 SparseModelType, typename GeometryValueType> |
||||
|
class SparsePcaaQuantitativeQuery : public SparsePcaaQuery<SparseModelType, GeometryValueType> { |
||||
|
public: |
||||
|
|
||||
|
// Typedefs for simple geometric objects |
||||
|
typedef std::vector<GeometryValueType> Point; |
||||
|
typedef std::vector<GeometryValueType> WeightVector; |
||||
|
|
||||
|
/* |
||||
|
* Creates a new query for the Pareto curve approximation algorithm (Pcaa) |
||||
|
* @param preprocessorResult the result from preprocessing |
||||
|
*/ |
||||
|
SparsePcaaQuantitativeQuery(SparsePcaaPreprocessorReturnType<SparseModelType>& preprocessorResult); |
||||
|
|
||||
|
|
||||
|
/* |
||||
|
* Invokes the computation and retrieves the result |
||||
|
*/ |
||||
|
virtual std::unique_ptr<CheckResult> 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<storm::storage::geometry::Polytope<GeometryValueType>> const& polytope); |
||||
|
|
||||
|
uint_fast64_t indexOfOptimizingObjective; |
||||
|
|
||||
|
Point thresholds; |
||||
|
storm::storage::BitVector strictThresholds; |
||||
|
std::shared_ptr<storm::storage::geometry::Polytope<GeometryValueType>> thresholdsAsPolytope; |
||||
|
}; |
||||
|
|
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
#endif /* STORM_MODELCHECKER_MULTIOBJECTIVE_PCAA_SPARSEPCAAQUANTITATIVEQUERY_H_ */ |
||||
@ -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 <class SparseModelType, typename GeometryValueType> |
||||
|
SparsePcaaQuery<SparseModelType, GeometryValueType>::SparsePcaaQuery(SparsePcaaPreprocessorReturnType<SparseModelType>& 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<GeometryValueType>::createUniversalPolytope(); |
||||
|
this->underApproximation = storm::storage::geometry::Polytope<GeometryValueType>::createEmptyPolytope(); |
||||
|
} |
||||
|
|
||||
|
template<> |
||||
|
void SparsePcaaQuery<storm::models::sparse::Mdp<double>, storm::RationalNumber>::initializeWeightVectorChecker(storm::models::sparse::Mdp<double> const& model, std::vector<PcaaObjective<double>> const& objectives, storm::storage::BitVector const& actionsWithNegativeReward, storm::storage::BitVector const& ecActions, storm::storage::BitVector const& possiblyRecurrentStates) { |
||||
|
this->weightVectorChecker = std::unique_ptr<SparsePcaaWeightVectorChecker<storm::models::sparse::Mdp<double>>>(new SparseMdpPcaaWeightVectorChecker<storm::models::sparse::Mdp<double>>(model, objectives, actionsWithNegativeReward, ecActions, possiblyRecurrentStates)); |
||||
|
} |
||||
|
|
||||
|
template<> |
||||
|
void SparsePcaaQuery<storm::models::sparse::Mdp<storm::RationalNumber>, storm::RationalNumber>::initializeWeightVectorChecker(storm::models::sparse::Mdp<storm::RationalNumber> const& model, std::vector<PcaaObjective<storm::RationalNumber>> const& objectives, storm::storage::BitVector const& actionsWithNegativeReward, storm::storage::BitVector const& ecActions, storm::storage::BitVector const& possiblyRecurrentStates) { |
||||
|
this->weightVectorChecker = std::unique_ptr<SparsePcaaWeightVectorChecker<storm::models::sparse::Mdp<storm::RationalNumber>>>(new SparseMdpPcaaWeightVectorChecker<storm::models::sparse::Mdp<storm::RationalNumber>>(model, objectives, actionsWithNegativeReward, ecActions, possiblyRecurrentStates)); |
||||
|
} |
||||
|
|
||||
|
template<> |
||||
|
void SparsePcaaQuery<storm::models::sparse::MarkovAutomaton<double>, storm::RationalNumber>::initializeWeightVectorChecker(storm::models::sparse::MarkovAutomaton<double> const& model, std::vector<PcaaObjective<double>> const& objectives, storm::storage::BitVector const& actionsWithNegativeReward, storm::storage::BitVector const& ecActions, storm::storage::BitVector const& possiblyRecurrentStates) { |
||||
|
this->weightVectorChecker = std::unique_ptr<SparsePcaaWeightVectorChecker<storm::models::sparse::MarkovAutomaton<double>>>(new SparseMaPcaaWeightVectorChecker<storm::models::sparse::MarkovAutomaton<double>>(model, objectives, actionsWithNegativeReward, ecActions, possiblyRecurrentStates)); |
||||
|
} |
||||
|
|
||||
|
// template<>
|
||||
|
// void SparsePcaaQuery<storm::models::sparse::MarkovAutomaton<storm::RationalNumber>, storm::RationalNumber>::initializeWeightVectorChecker(storm::models::sparse::MarkovAutomaton<storm::RationalNumber> const& model, std::vector<PcaaObjective<storm::RationalNumber>> const& objectives, storm::storage::BitVector const& actionsWithNegativeReward, storm::storage::BitVector const& ecActions, storm::storage::BitVector const& possiblyRecurrentStates) {
|
||||
|
// this->weightVectorChecker = std::unique_ptr<SparsePcaaWeightVectorChecker<storm::models::sparse::MarkovAutomaton<storm::RationalNumber>>>(new SparseMaPcaaWeightVectorChecker<storm::models::sparse::MarkovAutomaton<storm::RationalNumber>>(model, objectives, actionsWithNegativeReward, ecActions, possiblyRecurrentStates));
|
||||
|
// }
|
||||
|
|
||||
|
|
||||
|
template <class SparseModelType, typename GeometryValueType> |
||||
|
typename SparsePcaaQuery<SparseModelType, GeometryValueType>::WeightVector SparsePcaaQuery<SparseModelType, GeometryValueType>::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<double>(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<GeometryValueType>()); |
||||
|
result[objIndex] = storm::utility::one<GeometryValueType>(); |
||||
|
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<storm::storage::geometry::Halfspace<GeometryValueType>> halfspaces = underApproximation->getHalfspaces(); |
||||
|
uint_fast64_t farestHalfspaceIndex = halfspaces.size(); |
||||
|
GeometryValueType farestDistance = -storm::utility::one<GeometryValueType>(); |
||||
|
bool foundSeparatingDiracVector = false; |
||||
|
for(uint_fast64_t halfspaceIndex = 0; halfspaceIndex < halfspaces.size(); ++halfspaceIndex) { |
||||
|
GeometryValueType distance = -halfspaces[halfspaceIndex].euclideanDistance(pointToBeSeparated); |
||||
|
if(distance >= storm::utility::zero<GeometryValueType>()) { |
||||
|
storm::storage::BitVector nonZeroVectorEntries = ~storm::utility::vector::filterZero<GeometryValueType>(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<GeometryValueType>(halfspaces[farestHalfspaceIndex].normalVector()); |
||||
|
} |
||||
|
|
||||
|
STORM_LOG_THROW(farestHalfspaceIndex<halfspaces.size(), storm::exceptions::UnexpectedException, "There is no seperating vector."); |
||||
|
STORM_LOG_DEBUG("Found separating weight vector: " << storm::utility::vector::toString(storm::utility::vector::convertNumericVector<double>(halfspaces[farestHalfspaceIndex].normalVector())) << "."); |
||||
|
return halfspaces[farestHalfspaceIndex].normalVector(); |
||||
|
} |
||||
|
|
||||
|
template <class SparseModelType, typename GeometryValueType> |
||||
|
void SparsePcaaQuery<SparseModelType, GeometryValueType>::performRefinementStep(WeightVector&& direction) { |
||||
|
// Normalize the direction vector so that the entries sum up to one
|
||||
|
storm::utility::vector::scaleVectorInPlace(direction, storm::utility::one<GeometryValueType>() / std::accumulate(direction.begin(), direction.end(), storm::utility::zero<GeometryValueType>())); |
||||
|
// 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<typename SparseModelType::ValueType> 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<GeometryValueType> 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<typename SparseModelType::ValueType>(distance) * storm::utility::convertNumber<typename SparseModelType::ValueType>(0.5))); |
||||
|
break; |
||||
|
} |
||||
|
} |
||||
|
weightVectorChecker->check(storm::utility::vector::convertNumericVector<typename SparseModelType::ValueType>(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<double>(weightVectorChecker->getLowerBoundsOfInitialStateResults()))); |
||||
|
RefinementStep step; |
||||
|
step.weightVector = direction; |
||||
|
step.lowerBoundPoint = storm::utility::vector::convertNumericVector<GeometryValueType>(weightVectorChecker->getLowerBoundsOfInitialStateResults()); |
||||
|
step.upperBoundPoint = storm::utility::vector::convertNumericVector<GeometryValueType>(weightVectorChecker->getUpperBoundsOfInitialStateResults()); |
||||
|
refinementSteps.push_back(std::move(step)); |
||||
|
|
||||
|
updateOverApproximation(); |
||||
|
updateUnderApproximation(); |
||||
|
} |
||||
|
|
||||
|
template <class SparseModelType, typename GeometryValueType> |
||||
|
void SparsePcaaQuery<SparseModelType, GeometryValueType>::updateOverApproximation() { |
||||
|
storm::storage::geometry::Halfspace<GeometryValueType> 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<double>(h.euclideanDistance(refinementSteps.back().upperBoundPoint)) << "."); |
||||
|
} |
||||
|
overApproximation = overApproximation->intersection(h); |
||||
|
STORM_LOG_DEBUG("Updated OverApproximation to " << overApproximation->toString(true)); |
||||
|
} |
||||
|
|
||||
|
template <class SparseModelType, typename GeometryValueType> |
||||
|
void SparsePcaaQuery<SparseModelType, GeometryValueType>::updateUnderApproximation() { |
||||
|
std::vector<Point> paretoPoints; |
||||
|
paretoPoints.reserve(refinementSteps.size()); |
||||
|
for(auto const& step : refinementSteps) { |
||||
|
paretoPoints.push_back(step.lowerBoundPoint); |
||||
|
} |
||||
|
underApproximation = storm::storage::geometry::Polytope<GeometryValueType>::createDownwardClosure(paretoPoints); |
||||
|
STORM_LOG_DEBUG("Updated UnderApproximation to " << underApproximation->toString(true)); |
||||
|
} |
||||
|
|
||||
|
template <class SparseModelType, typename GeometryValueType> |
||||
|
bool SparsePcaaQuery<SparseModelType, GeometryValueType>::maxStepsPerformed() const { |
||||
|
return storm::settings::getModule<storm::settings::modules::MultiObjectiveSettings>().isMaxStepsSet() && |
||||
|
this->refinementSteps.size() >= storm::settings::getModule<storm::settings::modules::MultiObjectiveSettings>().getMaxSteps(); |
||||
|
} |
||||
|
|
||||
|
|
||||
|
template<typename SparseModelType, typename GeometryValueType> |
||||
|
typename SparsePcaaQuery<SparseModelType, GeometryValueType>::Point SparsePcaaQuery<SparseModelType, GeometryValueType>::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<GeometryValueType>(this->objectives[objIndex].toOriginalValueTransformationFactor) + storm::utility::convertNumber<GeometryValueType>(this->objectives[objIndex].toOriginalValueTransformationOffset)); |
||||
|
} |
||||
|
return result; |
||||
|
} |
||||
|
|
||||
|
template<typename SparseModelType, typename GeometryValueType> |
||||
|
std::shared_ptr<storm::storage::geometry::Polytope<GeometryValueType>> SparsePcaaQuery<SparseModelType, GeometryValueType>::transformPolytopeToOriginalModel(std::shared_ptr<storm::storage::geometry::Polytope<GeometryValueType>> const& polytope) const { |
||||
|
if(polytope->isEmpty()) { |
||||
|
return storm::storage::geometry::Polytope<GeometryValueType>::createEmptyPolytope(); |
||||
|
} |
||||
|
if(polytope->isUniversal()) { |
||||
|
return storm::storage::geometry::Polytope<GeometryValueType>::createUniversalPolytope(); |
||||
|
} |
||||
|
uint_fast64_t numObjectives = this->objectives.size(); |
||||
|
std::vector<std::vector<GeometryValueType>> transformationMatrix(numObjectives, std::vector<GeometryValueType>(numObjectives, storm::utility::zero<GeometryValueType>())); |
||||
|
std::vector<GeometryValueType> transformationVector; |
||||
|
transformationVector.reserve(numObjectives); |
||||
|
for(uint_fast64_t objIndex = 0; objIndex < numObjectives; ++objIndex) { |
||||
|
transformationMatrix[objIndex][objIndex] = storm::utility::convertNumber<GeometryValueType>(this->objectives[objIndex].toOriginalValueTransformationFactor); |
||||
|
transformationVector.push_back(storm::utility::convertNumber<GeometryValueType>(this->objectives[objIndex].toOriginalValueTransformationOffset)); |
||||
|
} |
||||
|
return polytope->affineTransformation(transformationMatrix, transformationVector); |
||||
|
} |
||||
|
|
||||
|
#ifdef STORM_HAVE_CARL
|
||||
|
template class SparsePcaaQuery<storm::models::sparse::Mdp<double>, storm::RationalNumber>; |
||||
|
template class SparsePcaaQuery<storm::models::sparse::MarkovAutomaton<double>, storm::RationalNumber>; |
||||
|
|
||||
|
template class SparsePcaaQuery<storm::models::sparse::Mdp<storm::RationalNumber>, storm::RationalNumber>; |
||||
|
// template class SparsePcaaQuery<storm::models::sparse::MarkovAutomaton<storm::RationalNumber>, storm::RationalNumber>;
|
||||
|
#endif
|
||||
|
} |
||||
|
} |
||||
|
} |
||||
@ -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 SparseModelType, typename GeometryValueType> |
||||
|
class SparsePcaaQuery { |
||||
|
public: |
||||
|
|
||||
|
// Typedefs for simple geometric objects |
||||
|
typedef std::vector<GeometryValueType> Point; |
||||
|
typedef std::vector<GeometryValueType> WeightVector; |
||||
|
|
||||
|
/* |
||||
|
* Invokes the computation and retrieves the result |
||||
|
*/ |
||||
|
virtual std::unique_ptr<CheckResult> check() = 0; |
||||
|
|
||||
|
protected: |
||||
|
|
||||
|
/* |
||||
|
* Initializes the weight vector checker with the provided data from preprocessing |
||||
|
*/ |
||||
|
void initializeWeightVectorChecker(SparseModelType const& model, |
||||
|
std::vector<PcaaObjective<typename SparseModelType::ValueType>> 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<SparseModelType>& 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<storm::storage::geometry::Polytope<GeometryValueType>> transformPolytopeToOriginalModel(std::shared_ptr<storm::storage::geometry::Polytope<GeometryValueType>> const& polytope) const; |
||||
|
|
||||
|
|
||||
|
SparseModelType const& originalModel; |
||||
|
storm::logic::MultiObjectiveFormula const& originalFormula; |
||||
|
|
||||
|
SparseModelType preprocessedModel; |
||||
|
std::vector<PcaaObjective<typename SparseModelType::ValueType>> objectives; |
||||
|
|
||||
|
|
||||
|
// The corresponding weight vector checker |
||||
|
std::unique_ptr<SparsePcaaWeightVectorChecker<SparseModelType>> weightVectorChecker; |
||||
|
|
||||
|
//The results in each iteration of the algorithm |
||||
|
std::vector<RefinementStep> refinementSteps; |
||||
|
//Overapproximation of the set of achievable values |
||||
|
std::shared_ptr<storm::storage::geometry::Polytope<GeometryValueType>> overApproximation; |
||||
|
//Underapproximation of the set of achievable values |
||||
|
std::shared_ptr<storm::storage::geometry::Polytope<GeometryValueType>> 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_ */ |
||||
@ -0,0 +1,311 @@ |
|||||
|
#include "src/modelchecker/multiobjective/pcaa/SparsePcaaWeightVectorChecker.h"
|
||||
|
|
||||
|
#include <map>
|
||||
|
|
||||
|
#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 <class SparseModelType> |
||||
|
SparsePcaaWeightVectorChecker<SparseModelType>::SparsePcaaWeightVectorChecker(SparseModelType const& model, |
||||
|
std::vector<PcaaObjective<typename SparseModelType::ValueType>> 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 <class SparseModelType> |
||||
|
void SparsePcaaWeightVectorChecker<SparseModelType>::check(std::vector<ValueType> const& weightVector) { |
||||
|
checkHasBeenCalled=true; |
||||
|
STORM_LOG_DEBUG("Invoked WeightVectorChecker with weights " << std::endl << "\t" << storm::utility::vector::toString(storm::utility::vector::convertNumericVector<double>(weightVector))); |
||||
|
std::vector<ValueType> weightedRewardVector(model.getTransitionMatrix().getRowCount(), storm::utility::zero<ValueType>()); |
||||
|
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<double>(getLowerBoundsOfInitialStateResults()))); |
||||
|
} |
||||
|
|
||||
|
template <class SparseModelType> |
||||
|
void SparsePcaaWeightVectorChecker<SparseModelType>::setMaximumLowerUpperBoundGap(ValueType const& value) { |
||||
|
this->maximumLowerUpperBoundGap = value; |
||||
|
} |
||||
|
|
||||
|
template <class SparseModelType> |
||||
|
typename SparsePcaaWeightVectorChecker<SparseModelType>::ValueType const& SparsePcaaWeightVectorChecker<SparseModelType>::getMaximumLowerUpperBoundGap() const { |
||||
|
return this->maximumLowerUpperBoundGap; |
||||
|
} |
||||
|
|
||||
|
template <class SparseModelType> |
||||
|
std::vector<typename SparsePcaaWeightVectorChecker<SparseModelType>::ValueType> SparsePcaaWeightVectorChecker<SparseModelType>::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<ValueType> 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 <class SparseModelType> |
||||
|
std::vector<typename SparsePcaaWeightVectorChecker<SparseModelType>::ValueType> SparsePcaaWeightVectorChecker<SparseModelType>::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<ValueType> 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 <class SparseModelType> |
||||
|
storm::storage::TotalScheduler const& SparsePcaaWeightVectorChecker<SparseModelType>::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 <class SparseModelType> |
||||
|
void SparsePcaaWeightVectorChecker<SparseModelType>::unboundedWeightedPhase(std::vector<ValueType> const& weightedRewardVector) { |
||||
|
if(this->objectivesWithNoUpperTimeBound.empty()) { |
||||
|
this->weightedResult = std::vector<ValueType>(model.getNumberOfStates(), storm::utility::zero<ValueType>()); |
||||
|
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<ValueType>::transform(model.getTransitionMatrix(), subsystemStates, ecActions & zeroRewardActions, possiblyRecurrentStates); |
||||
|
|
||||
|
std::vector<ValueType> subRewardVector(ecEliminatorResult.newToOldRowMapping.size()); |
||||
|
storm::utility::vector::selectVectorValues(subRewardVector, ecEliminatorResult.newToOldRowMapping, weightedRewardVector); |
||||
|
std::vector<ValueType> subResult(ecEliminatorResult.matrix.getRowGroupCount()); |
||||
|
|
||||
|
storm::solver::GeneralMinMaxLinearEquationSolverFactory<ValueType> solverFactory; |
||||
|
std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = solverFactory.create(ecEliminatorResult.matrix); |
||||
|
solver->setOptimizationDirection(storm::solver::OptimizationDirection::Maximize); |
||||
|
solver->setTrackScheduler(true); |
||||
|
solver->solveEquations(subResult, subRewardVector); |
||||
|
|
||||
|
this->weightedResult = std::vector<ValueType>(model.getNumberOfStates()); |
||||
|
std::vector<uint_fast64_t> 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 <class SparseModelType> |
||||
|
void SparsePcaaWeightVectorChecker<SparseModelType>::unboundedIndividualPhase(std::vector<ValueType> const& weightVector) { |
||||
|
|
||||
|
storm::storage::SparseMatrix<ValueType> deterministicMatrix = model.getTransitionMatrix().selectRowsFromRowGroups(this->scheduler.getChoices(), true); |
||||
|
storm::storage::SparseMatrix<ValueType> deterministicBackwardTransitions = deterministicMatrix.transpose(); |
||||
|
std::vector<ValueType> deterministicStateRewards(deterministicMatrix.getRowCount()); |
||||
|
storm::solver::GeneralLinearEquationSolverFactory<ValueType> 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<ValueType> 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<ValueType>(); |
||||
|
offsetsToUpperBound[objIndex] = storm::utility::zero<ValueType>(); |
||||
|
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<ValueType>() / sumOfWeightsOfUncheckedObjectives); |
||||
|
} |
||||
|
|
||||
|
// Make sure that the objectiveResult is initialized in some way
|
||||
|
objectiveResults[objIndex].resize(model.getNumberOfStates(), storm::utility::zero<ValueType>()); |
||||
|
|
||||
|
// Invoke the linear equation solver
|
||||
|
objectiveResults[objIndex] = storm::modelchecker::helper::SparseDtmcPrctlHelper<ValueType>::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<ValueType>(model.getNumberOfStates(), storm::utility::zero<ValueType>()); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
template <class SparseModelType> |
||||
|
void SparsePcaaWeightVectorChecker<SparseModelType>::transformReducedSolutionToOriginalModel(storm::storage::SparseMatrix<ValueType> const& reducedMatrix, |
||||
|
std::vector<ValueType> const& reducedSolution, |
||||
|
std::vector<uint_fast64_t> const& reducedOptimalChoices, |
||||
|
std::vector<uint_fast64_t> const& reducedToOriginalChoiceMapping, |
||||
|
std::vector<uint_fast64_t> const& originalToReducedStateMapping, |
||||
|
std::vector<ValueType>& originalSolution, |
||||
|
std::vector<uint_fast64_t>& 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<ValueType>(); |
||||
|
// 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<storm::models::sparse::Mdp<double>>; |
||||
|
template class SparsePcaaWeightVectorChecker<storm::models::sparse::MarkovAutomaton<double>>; |
||||
|
#ifdef STORM_HAVE_CARL
|
||||
|
template class SparsePcaaWeightVectorChecker<storm::models::sparse::Mdp<storm::RationalNumber>>; |
||||
|
template class SparsePcaaWeightVectorChecker<storm::models::sparse::MarkovAutomaton<storm::RationalNumber>>; |
||||
|
#endif
|
||||
|
|
||||
|
} |
||||
|
} |
||||
|
} |
||||
@ -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 SparseModelType> |
||||
|
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<PcaaObjective<ValueType>> 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<ValueType> 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<ValueType> getLowerBoundsOfInitialStateResults() const; |
||||
|
std::vector<ValueType> 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<ValueType> 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<ValueType> 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<ValueType> const& weightVector, std::vector<ValueType>& 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<ValueType> const& reducedMatrix, |
||||
|
std::vector<ValueType> const& reducedSolution, |
||||
|
std::vector<uint_fast64_t> const& reducedOptimalChoices, |
||||
|
std::vector<uint_fast64_t> const& reducedToOriginalChoiceMapping, |
||||
|
std::vector<uint_fast64_t> const& originalToReducedStateMapping, |
||||
|
std::vector<ValueType>& originalSolution, |
||||
|
std::vector<uint_fast64_t>& originalOptimalChoices) const; |
||||
|
|
||||
|
// The (preprocessed) model |
||||
|
SparseModelType const& model; |
||||
|
// The (preprocessed) objectives |
||||
|
std::vector<PcaaObjective<ValueType>> 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<std::vector<ValueType>> 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<ValueType> weightedResult; |
||||
|
// The lower bounds of the results for the individual objectives (w.r.t. all states of the model) |
||||
|
std::vector<std::vector<ValueType>> 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<ValueType> offsetsToLowerBound; |
||||
|
std::vector<ValueType> offsetsToUpperBound; |
||||
|
// The scheduler that maximizes the weighted rewards |
||||
|
storm::storage::TotalScheduler scheduler; |
||||
|
|
||||
|
}; |
||||
|
|
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
#endif /* STORM_MODELCHECKER_MULTIOBJECTIVE_PCAA_SPARSEPCAAWEIGHTVECTORCHECKER_H_ */ |
||||
Write
Preview
Loading…
Cancel
Save
Reference in new issue