You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
 
 
 
 

395 lines
29 KiB

#include "storm/modelchecker/multiobjective/pcaa/SparsePcaaWeightVectorChecker.h"
#include <map>
#include "storm/adapters/RationalFunctionAdapter.h"
#include "storm/models/sparse/Mdp.h"
#include "storm/models/sparse/MarkovAutomaton.h"
#include "storm/models/sparse/StandardRewardModel.h"
#include "storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.h"
#include "storm/solver/MinMaxLinearEquationSolver.h"
#include "storm/transformer/EndComponentEliminator.h"
#include "storm/utility/graph.h"
#include "storm/utility/macros.h"
#include "storm/utility/vector.h"
#include "storm/logic/Formulas.h"
#include "storm/exceptions/IllegalFunctionCallException.h"
#include "storm/exceptions/UnexpectedException.h"
#include "storm/exceptions/NotImplementedException.h"
namespace storm {
namespace modelchecker {
namespace multiobjective {
template <class SparseModelType>
SparsePcaaWeightVectorChecker<SparseModelType>::SparsePcaaWeightVectorChecker(SparseModelType const& model,
std::vector<Objective<typename SparseModelType::ValueType>> const& objectives,
storm::storage::BitVector const& possibleECActions,
storm::storage::BitVector const& possibleBottomStates) :
model(model),
objectives(objectives),
possibleECActions(possibleECActions),
actionsWithoutRewardInUnboundedPhase(model.getNumberOfChoices(), true),
possibleBottomStates(possibleBottomStates),
objectivesWithNoUpperTimeBound(objectives.size(), false),
discreteActionRewards(objectives.size()),
checkHasBeenCalled(false),
objectiveResults(objectives.size()),
offsetsToUnderApproximation(objectives.size()),
offsetsToOverApproximation(objectives.size()),
optimalChoices(model.getNumberOfStates(), 0) {
// set data for unbounded objectives
for(uint_fast64_t objIndex = 0; objIndex < objectives.size(); ++objIndex) {
auto const& formula = *objectives[objIndex].formula;
if (formula.getSubformula().isTotalRewardFormula()) {
objectivesWithNoUpperTimeBound.set(objIndex, true);
STORM_LOG_ASSERT(formula.isRewardOperatorFormula() && formula.asRewardOperatorFormula().hasRewardModelName(), "Unexpected type of operator formula.");
actionsWithoutRewardInUnboundedPhase &= model.getRewardModel(formula.asRewardOperatorFormula().getRewardModelName()).getChoicesWithZeroReward(model.getTransitionMatrix());
}
}
}
template <class SparseModelType>
void SparsePcaaWeightVectorChecker<SparseModelType>::check(std::vector<ValueType> const& weightVector) {
checkHasBeenCalled=true;
STORM_LOG_INFO("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>());
boost::optional<ValueType> weightedLowerResultBound = storm::utility::zero<ValueType>();
boost::optional<ValueType> weightedUpperResultBound = storm::utility::zero<ValueType>();
for (auto objIndex : objectivesWithNoUpperTimeBound) {
auto const& obj = objectives[objIndex];
if (storm::solver::minimize(objectives[objIndex].formula->getOptimalityType())) {
if (obj.lowerResultBound && weightedUpperResultBound) {
weightedUpperResultBound.get() -= weightVector[objIndex] * obj.lowerResultBound.get();
} else {
weightedUpperResultBound = boost::none;
}
if (obj.upperResultBound && weightedLowerResultBound) {
weightedLowerResultBound.get() -= weightVector[objIndex] * obj.upperResultBound.get();
} else {
weightedLowerResultBound = boost::none;
}
storm::utility::vector::addScaledVector(weightedRewardVector, discreteActionRewards[objIndex], -weightVector[objIndex]);
} else {
if (obj.lowerResultBound && weightedLowerResultBound) {
weightedLowerResultBound.get() += weightVector[objIndex] * obj.lowerResultBound.get();
} else {
weightedLowerResultBound = boost::none;
}
if (obj.upperResultBound && weightedUpperResultBound) {
weightedUpperResultBound.get() += weightVector[objIndex] * obj.upperResultBound.get();
} else {
weightedUpperResultBound = boost::none;
}
storm::utility::vector::addScaledVector(weightedRewardVector, discreteActionRewards[objIndex], weightVector[objIndex]);
}
}
unboundedWeightedPhase(weightedRewardVector, weightedLowerResultBound, weightedUpperResultBound);
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.formula->getSubformula().isTotalRewardFormula()) {
boundedPhase(weightVector, weightedRewardVector);
break;
}
}
STORM_LOG_INFO("Weight vector check done. Lower bounds for results in initial state: " << storm::utility::vector::toString(storm::utility::vector::convertNumericVector<double>(getUnderApproximationOfInitialStateResults())));
// Validate that the results are sufficiently precise
ValueType resultingWeightedPrecision = storm::utility::abs<ValueType>(storm::utility::vector::dotProduct(getOverApproximationOfInitialStateResults(), weightVector) - storm::utility::vector::dotProduct(getUnderApproximationOfInitialStateResults(), weightVector));
resultingWeightedPrecision /= storm::utility::sqrt(storm::utility::vector::dotProduct(weightVector, weightVector));
STORM_LOG_THROW(resultingWeightedPrecision <= weightedPrecision, storm::exceptions::UnexpectedException, "The desired precision was not reached");
}
template <class SparseModelType>
void SparsePcaaWeightVectorChecker<SparseModelType>::setWeightedPrecision(ValueType const& weightedPrecision) {
this->weightedPrecision = weightedPrecision;
}
template <class SparseModelType>
typename SparsePcaaWeightVectorChecker<SparseModelType>::ValueType const& SparsePcaaWeightVectorChecker<SparseModelType>::getWeightedPrecision() const {
return this->weightedPrecision;
}
template <class SparseModelType>
std::vector<typename SparsePcaaWeightVectorChecker<SparseModelType>::ValueType> SparsePcaaWeightVectorChecker<SparseModelType>::getUnderApproximationOfInitialStateResults() 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->offsetsToUnderApproximation[objIndex]);
}
return res;
}
template <class SparseModelType>
std::vector<typename SparsePcaaWeightVectorChecker<SparseModelType>::ValueType> SparsePcaaWeightVectorChecker<SparseModelType>::getOverApproximationOfInitialStateResults() 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->offsetsToOverApproximation[objIndex]);
}
return res;
}
template <class SparseModelType>
storm::storage::Scheduler<typename SparsePcaaWeightVectorChecker<SparseModelType>::ValueType> SparsePcaaWeightVectorChecker<SparseModelType>::computeScheduler() 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.formula->getSubformula().isTotalRewardFormula(), storm::exceptions::NotImplementedException, "Scheduler retrival is only implemented for objectives without time-bound.");
}
storm::storage::Scheduler<ValueType> result(this->optimalChoices.size());
uint_fast64_t state = 0;
for (auto const& choice : optimalChoices) {
result.setChoice(choice, state);
++state;
}
return result;
}
template <class SparseModelType>
void SparsePcaaWeightVectorChecker<SparseModelType>::unboundedWeightedPhase(std::vector<ValueType> const& weightedRewardVector, boost::optional<ValueType> const& lowerResultBound, boost::optional<ValueType> const& upperResultBound) {
if (this->objectivesWithNoUpperTimeBound.empty() || !storm::utility::vector::hasNonZeroEntry(weightedRewardVector)) {
this->weightedResult = std::vector<ValueType>(model.getNumberOfStates(), storm::utility::zero<ValueType>());
this->optimalChoices = std::vector<uint_fast64_t>(model.getNumberOfStates(), 0);
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().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, possibleECActions & zeroRewardActions, possibleBottomStates);
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);
if (lowerResultBound) {
solver->setLowerBound(*lowerResultBound);
}
if (upperResultBound) {
solver->setUpperBound(*upperResultBound);
}
solver->solveEquations(subResult, subRewardVector);
this->weightedResult = std::vector<ValueType>(model.getNumberOfStates());
transformReducedSolutionToOriginalModel(ecEliminatorResult.matrix, subResult, solver->getSchedulerChoices(), ecEliminatorResult.newToOldRowMapping, ecEliminatorResult.oldToNewStateMapping, this->weightedResult, this->optimalChoices);
}
template <class SparseModelType>
void SparsePcaaWeightVectorChecker<SparseModelType>::unboundedIndividualPhase(std::vector<ValueType> const& weightVector) {
if (objectivesWithNoUpperTimeBound.getNumberOfSetBits() == 1 && storm::utility::isOne(weightVector[*objectivesWithNoUpperTimeBound.begin()])) {
uint_fast64_t objIndex = *objectivesWithNoUpperTimeBound.begin();
objectiveResults[objIndex] = weightedResult;
if (storm::solver::minimize(objectives[objIndex].formula->getOptimalityType())) {
storm::utility::vector::scaleVectorInPlace(objectiveResults[objIndex], -storm::utility::one<ValueType>());
}
for (uint_fast64_t objIndex2 = 0; objIndex2 < objectives.size(); ++objIndex2) {
if (objIndex != objIndex2) {
objectiveResults[objIndex2] = std::vector<ValueType>(model.getNumberOfStates(), storm::utility::zero<ValueType>());
}
}
} else {
storm::storage::SparseMatrix<ValueType> deterministicMatrix = model.getTransitionMatrix().selectRowsFromRowGroups(this->optimalChoices, 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)) {
auto const& obj = objectives[objIndex];
if (objectivesWithNoUpperTimeBound.get(objIndex)) {
offsetsToUnderApproximation[objIndex] = storm::utility::zero<ValueType>();
offsetsToOverApproximation[objIndex] = storm::utility::zero<ValueType>();
storm::utility::vector::selectVectorValues(deterministicStateRewards, this->optimalChoices, model.getTransitionMatrix().getRowGroupIndices(), discreteActionRewards[objIndex]);
storm::storage::BitVector statesWithRewards = ~storm::utility::vector::filterZero(deterministicStateRewards);
// As maybestates we pick the states from which a state with reward is reachable
storm::storage::BitVector maybeStates = 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;
ValueType scalingFactor = storm::utility::one<ValueType>() / sumOfWeightsOfUncheckedObjectives;
if (storm::solver::minimize(obj.formula->getOptimalityType())) {
scalingFactor *= -storm::utility::one<ValueType>();
}
storm::utility::vector::scaleVectorInPlace(objectiveResults[objIndex], scalingFactor);
storm::utility::vector::clip(objectiveResults[objIndex], obj.lowerResultBound, obj.upperResultBound);
}
// Make sure that the objectiveResult is initialized correctly
objectiveResults[objIndex].resize(model.getNumberOfStates(), storm::utility::zero<ValueType>());
if (!maybeStates.empty()) {
storm::storage::SparseMatrix<ValueType> submatrix = deterministicMatrix.getSubmatrix(
true, maybeStates, maybeStates, true);
// Converting the matrix from the fixpoint notation to the form needed for the equation
// system. That is, we go from x = A*x + b to (I-A)x = b.
submatrix.convertToEquationSystem();
// Prepare solution vector and rhs of the equation system.
std::vector<ValueType> x = storm::utility::vector::filterVector(objectiveResults[objIndex], maybeStates);
std::vector<ValueType> b = storm::utility::vector::filterVector(deterministicStateRewards, maybeStates);
// Now solve the resulting equation system.
std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(std::move(submatrix));
if (obj.lowerResultBound) {
solver->setLowerBound(*obj.lowerResultBound);
}
if (obj.upperResultBound) {
solver->setUpperBound(*obj.upperResultBound);
}
solver->solveEquations(x, b);
// Set the result for this objective accordingly
storm::utility::vector::setVectorValues<ValueType>(objectiveResults[objIndex], maybeStates, x);
}
storm::utility::vector::setVectorValues<ValueType>(objectiveResults[objIndex], ~maybeStates, storm::utility::zero<ValueType>());
// 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 bottomStates(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 a bottom state, i.e., the chosen row stays inside its EC.
bool stateIsBottom = possibleBottomStates.get(state);
for(auto const& entry : model.getTransitionMatrix().getRow(chosenRowInOriginalModel)) {
stateIsBottom &= originalToReducedStateMapping[entry.getColumn()] == stateInReducedModel;
}
if(stateIsBottom) {
bottomStates.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 possible bottom states are made bottom and the remaining states have to reach a bottom state with prob. one
if(possibleBottomStates.get(state)) {
bottomStates.set(state);
} else {
statesWithUndefSched.set(state);
}
}
}
// Handle bottom states
for(auto state : bottomStates) {
bool foundRowForState = false;
// Find a row with zero rewards that only leads to bottom 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 rowOnlyLeadsToBottomStates = true;
bool rowStaysInEC = true;
for( auto const& entry : model.getTransitionMatrix().getRow(row)) {
rowOnlyLeadsToBottomStates &= bottomStates.get(entry.getColumn());
rowStaysInEC &= originalToReducedStateMapping[entry.getColumn()] == stateInReducedModel;
}
if(rowOnlyLeadsToBottomStates && (rowStaysInEC || !statesThatShouldStayInTheirEC.get(state)) && actionsWithoutRewardInUnboundedPhase.get(row)) {
foundRowForState = true;
originalOptimalChoices[state] = row - model.getTransitionMatrix().getRowGroupIndices()[state];
break;
}
}
STORM_LOG_ASSERT(foundRowForState, "Could not find a suitable choice for a bottom state.");
}
// Handle remaining states with still undef. scheduler (either EC states or non-subsystem states)
while(!statesWithUndefSched.empty()) {
for(auto state : statesWithUndefSched) {
// Iteratively Try to find a choice such that at least one successor has a defined scheduler.
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
}
}
}