Browse Source

Merge branch 'quantiles_refactor'

main
Tim Quatmann 6 years ago
parent
commit
ee6942c563
  1. 20
      resources/examples/testfiles/dtmc/quantiles_simple_dtmc.pm
  2. 80
      resources/examples/testfiles/mdp/quantiles_firewire.nm
  3. 72
      resources/examples/testfiles/mdp/quantiles_resources.nm
  4. 31
      resources/examples/testfiles/mdp/quantiles_simple_mdp.nm
  5. 25
      src/storm-parsers/parser/FormulaParserGrammar.cpp
  6. 7
      src/storm-parsers/parser/FormulaParserGrammar.h
  7. 20
      src/storm/logic/QuantileFormula.cpp
  8. 9
      src/storm/logic/QuantileFormula.h
  9. 32
      src/storm/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp
  10. 1
      src/storm/modelchecker/prctl/SparseDtmcPrctlModelChecker.h
  11. 2
      src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp
  12. 186
      src/storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.cpp
  13. 60
      src/storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.h
  14. 11
      src/storm/modelchecker/prctl/helper/rewardbounded/EpochManager.cpp
  15. 1
      src/storm/modelchecker/prctl/helper/rewardbounded/EpochManager.h
  16. 96
      src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp
  17. 29
      src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h
  18. 17
      src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.cpp
  19. 8
      src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.h
  20. 666
      src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp
  21. 39
      src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h
  22. 302
      src/test/storm/modelchecker/QuantileQueryTest.cpp

20
resources/examples/testfiles/dtmc/quantiles_simple_dtmc.pm

@ -0,0 +1,20 @@
dtmc
module main
s : [0..3] init 0;
[] s=0 -> 0.4 : (s'=1) + 0.4 : (s'=2) + 0.2 : (s'=3);
[a] s=1 -> 1 : (s'=0);
[b] s=2 -> 1 : (s'=0);
[] s=3 -> 1 : (s'=3);
endmodule
rewards "first"
[a] true : 1;
endrewards
rewards "second"
[b] true : 2;
endrewards

80
resources/examples/testfiles/mdp/quantiles_firewire.nm

@ -0,0 +1,80 @@
// integer semantics version of abstract firewire protocol
// gxn 23/05/2001
mdp
// wire delay
const int delay;
// probability of choosing fast and slow
const double fast = 0.5;
const double slow = 1-fast;
// largest constant the clock of the system is compared to
const int kx = 167;
module abstract_firewire
// clock
x : [0..kx+1];
// local state
s : [0..9];
// 0 -start_start
// 1 -fast_start
// 2 -start_fast
// 3 -start_slow
// 4 -slow_start
// 5 -fast_fast
// 6 -fast_slow
// 7 -slow_fast
// 8 -slow_slow
// 9 -done
// initial state
[time] s=0 & x<delay -> (x'=min(x+1,kx+1));
[round] s=0 -> fast : (s'=1) + slow : (s'=4);
[round] s=0 -> fast : (s'=2) + slow : (s'=3);
// fast_start
[time] s=1 & x<delay -> (x'=min(x+1,kx+1));
[] s=1 -> fast : (s'=5) & (x'=0) + slow : (s'=6) & (x'=0);
// start_fast
[time] s=2 & x<delay -> (x'=min(x+1,kx+1));
[] s=2 -> fast : (s'=5) & (x'=0) + slow : (s'=7) & (x'=0);
// start_slow
[time] s=3 & x<delay -> (x'=min(x+1,kx+1));
[] s=3 -> fast : (s'=6) & (x'=0) + slow : (s'=8) & (x'=0);
// slow_start
[time] s=4 & x<delay -> (x'=min(x+1,kx+1));
[] s=4 -> fast : (s'=7) & (x'=0) + slow : (s'=8) & (x'=0);
// fast_fast
[time] s=5 & (x<85) -> (x'=min(x+1,kx+1));
[] s=5 & (x>=76) -> (s'=0) & (x'=0);
[] s=5 & (x>=76-delay) -> (s'=9) & (x'=0);
// fast_slow
[time] s=6 & x<167 -> (x'=min(x+1,kx+1));
[] s=6 & x>=159-delay -> (s'=9) & (x'=0);
// slow_fast
[time] s=7 & x<167 -> (x'=min(x+1,kx+1));
[] s=7 & x>=159-delay -> (s'=9) & (x'=0);
// slow_slow
[time] s=8 & x<167 -> (x'=min(x+1,kx+1));
[] s=8 & x>=159 -> (s'=0) & (x'=0);
[] s=8 & x>=159-delay -> (s'=9) & (x'=0);
// done
[] s=9 -> (s'=s);
endmodule
// labels
label "done" = s=9;
//reward structures
// time
rewards "time"
[time] true : 1;
endrewards
// number of rounds
rewards "rounds"
[round] true : 1;
endrewards

72
resources/examples/testfiles/mdp/quantiles_resources.nm

@ -0,0 +1,72 @@
mdp
const int WIDTH = 5;
const int HEIGHT = 5;
const int XINIT = 3;
const int YINIT = 1;
const double pAttack = 1/10;
formula left_of_gold = x=2 & y=5;
formula right_of_gold = x=4 & y=5;
formula below_of_gold = (x=3 & y=4);
formula above_of_gold = false;
formula left_of_gem = (x=4 & y=4);
formula right_of_gem = false;
formula below_of_gem = (x=5 & y=3);
formula above_of_gem = (x=5 & y=5);
formula left_of_home = x=2 & y=1;
formula right_of_home = x=4 & y=1;
formula above_of_home = x=3 & y=2;
formula below_of_home = false;
formula left_of_enemy = (x=3 & y=5) | (x=2 & y=4);
formula right_of_enemy = (x=5 & y=5) | (x=4 & y=4);
formula above_of_enemy = x=3 & y=5;
formula below_of_enemy = (x=3 & y=3) | (x=4 & y=4);
module robot
gold : bool init false;
gem : bool init false;
attacked : bool init false;
x : [1..WIDTH] init XINIT;
y : [1..HEIGHT] init YINIT;
[right] !left_of_enemy & x<WIDTH -> (attacked'=false) & (x'=x+1) & (gold' = (gold & !left_of_home) | left_of_gold) & (gem' = (gem & !left_of_home) | left_of_gem);
[left] !right_of_enemy & x>1 -> (attacked'=false) & (x'=x-1) & (gold' = (gold & !right_of_home) | right_of_gold) & (gem' = (gem & !right_of_home) | right_of_gem);
[top] !below_of_enemy & y<HEIGHT -> (attacked'=false) & (y'=y+1) & (gold' = (gold & !below_of_home) | below_of_gold) & (gem' = (gem & !below_of_home) | below_of_gem);
[down] !above_of_enemy & y>1 -> (attacked'=false) & (y'=y-1) & (gold' = (gold & !above_of_home) | above_of_gold) & (gem' = (gem & !above_of_home) | above_of_gem);
[right] left_of_enemy & x<WIDTH -> pAttack : (attacked'=true) & (x'=XINIT) & (y'=YINIT) & (gold'=false) & (gem'=false) + (1-pAttack) : (attacked'=false) & (x'=x+1) & (gold' = (gold & !left_of_home) | left_of_gold) & (gem' = (gem & !left_of_home) | left_of_gem);
[left] right_of_enemy & x>1 -> pAttack : (attacked'=true) & (x'=XINIT) & (y'=YINIT) & (gold'=false) & (gem'=false) + (1-pAttack) : (attacked'=false) & (x'=x-1) & (gold' = (gold & !right_of_home) | right_of_gold) & (gem' = (gem & !right_of_home) | right_of_gem);
[top] below_of_enemy & y<HEIGHT -> pAttack : (attacked'=true) & (x'=XINIT) & (y'=YINIT) & (gold'=false) & (gem'=false) + (1-pAttack) : (attacked'=false) & (y'=y+1) & (gold' = (gold & !below_of_home) | below_of_gold) & (gem' = (gem & !below_of_home) | below_of_gem);
[down] above_of_enemy & y>1 -> pAttack : (attacked'=true) & (x'=XINIT) & (y'=YINIT) & (gold'=false) & (gem'=false) + (1-pAttack) : (attacked'=false) & (y'=y-1) & (gold' = (gold & !above_of_home) | above_of_gold) & (gem' = (gem & !above_of_home) | above_of_gem);
endmodule
rewards "attacks"
attacked : 1;
endrewards
rewards "steps"
[left] true : 1;
[right] true : 1;
[top] true : 1;
[down] true : 1;
endrewards
rewards "gold"
[right] left_of_home & gold : 1;
[left] right_of_home & gold : 1;
[top] below_of_home & gold : 1;
[down] above_of_home & gold : 1;
endrewards
rewards "gem"
[right] left_of_home & gem : 1;
[left] right_of_home & gem : 1;
[top] below_of_home & gem : 1;
[down] above_of_home & gem: 1;
endrewards

31
resources/examples/testfiles/mdp/quantiles_simple_mdp.nm

@ -0,0 +1,31 @@
mdp
module main
s : [0..6] init 6;
[] s=6 -> 1 : (s'=0);
[] s=6 -> 1 : (s'=5);
[a] s=5 -> 1 : (s'=5);
[b] s=5 -> 1 : (s'=1);
[] s=0 -> 1/10: (s'=1) + 9/10: (s'=3);
[] s=0 -> 1/10: (s'=1) + 9/10: (s'=4);
[a] s=3 -> 1: (s'=0);
[b] s=3 -> 1: (s'=0);
[] s=1 -> 1: (s'=2);
[] s=2 -> 1: (s'=2);
endmodule
rewards "first"
[a] true : 1;
endrewards
rewards "second"
[b] s=3 : 1;
[b] s=5 : 2;
endrewards
rewards "third"
true : 0.7;
endrewards

25
src/storm-parsers/parser/FormulaParserGrammar.cpp

@ -127,9 +127,9 @@ namespace storm {
identifier %= qi::as_string[qi::raw[qi::lexeme[((qi::alpha | qi::char_('_') | qi::char_('.')) >> *(qi::alnum | qi::char_('_')))]]];
identifier.name("identifier");
quantileBoundVariable = ((qi::lit("min")[qi::_a = storm::solver::OptimizationDirection::Minimize] | qi::lit("max")[qi::_a = storm::solver::OptimizationDirection::Maximize]) >> identifier)[qi::_val = phoenix::bind(&FormulaParserGrammar::createQuantileBoundVariables, phoenix::ref(*this), qi::_a, qi::_1)];
quantileBoundVariable.name("Quantile bound variable");
quantileFormula = (qi::lit("quantile") > qi::lit("(") >> (quantileBoundVariable % qi::lit(",")) >> qi::lit(",") >> stateFormula > qi::lit(")"))[qi::_val = phoenix::bind(&FormulaParserGrammar::createQuantileFormula, phoenix::ref(*this), qi::_1, qi::_2)];
quantileBoundVariable = (-(qi::lit("min")[qi::_a = storm::solver::OptimizationDirection::Minimize] | qi::lit("max")[qi::_a = storm::solver::OptimizationDirection::Maximize]) >> identifier >> qi::lit(","))[qi::_val = phoenix::bind(&FormulaParserGrammar::createQuantileBoundVariables, phoenix::ref(*this), qi::_a, qi::_1)];
quantileBoundVariable.name("quantile bound variable");
quantileFormula = (qi::lit("quantile") > qi::lit("(") >> *(quantileBoundVariable) >> stateFormula > qi::lit(")"))[qi::_val = phoenix::bind(&FormulaParserGrammar::createQuantileFormula, phoenix::ref(*this), qi::_1, qi::_2)];
quantileFormula.name("Quantile formula");
stateFormula = (orStateFormula | multiFormula | quantileFormula);
@ -424,17 +424,22 @@ namespace storm {
}
}
std::pair<storm::solver::OptimizationDirection, storm::expressions::Variable> FormulaParserGrammar::createQuantileBoundVariables(storm::solver::OptimizationDirection const& dir, std::string const& variableName) {
storm::expressions::Variable FormulaParserGrammar::createQuantileBoundVariables(boost::optional<storm::solver::OptimizationDirection> const& dir, std::string const& variableName) {
STORM_LOG_ASSERT(manager, "Mutable expression manager required to define quantile bound variable.");
STORM_LOG_THROW(!manager->hasVariable(variableName), storm::exceptions::WrongFormatException, "Invalid quantile variable name '" << variableName << "' in property: variable already exists.");
storm::expressions::Variable var = manager->declareRationalVariable(variableName);
storm::expressions::Variable var;
if (manager->hasVariable(variableName)) {
var = manager->getVariable(variableName);
STORM_LOG_THROW(quantileFormulaVariables.count(var) > 0, storm::exceptions::WrongFormatException, "Invalid quantile variable name '" << variableName << "' in quantile formula: variable already exists.");
} else {
var = manager->declareRationalVariable(variableName);
quantileFormulaVariables.insert(var);
}
STORM_LOG_WARN_COND(!dir.is_initialized(), "Optimization direction '" << dir.get() << "' for quantile variable " << variableName << " is ignored. This information will be derived from the subformula of the quantile.");
addIdentifierExpression(variableName, var);
std::pair<storm::solver::OptimizationDirection, storm::expressions::Variable> result(dir, var);
return result;
return var;
}
std::shared_ptr<storm::logic::Formula const> FormulaParserGrammar::createQuantileFormula(std::vector<std::pair<storm::solver::OptimizationDirection, storm::expressions::Variable>> const& boundVariables, std::shared_ptr<storm::logic::Formula const> const& subformula) {
std::shared_ptr<storm::logic::Formula const> FormulaParserGrammar::createQuantileFormula(std::vector<storm::expressions::Variable> const& boundVariables, std::shared_ptr<storm::logic::Formula const> const& subformula) {
return std::shared_ptr<storm::logic::Formula const>(new storm::logic::QuantileFormula(boundVariables, subformula));
}

7
src/storm-parsers/parser/FormulaParserGrammar.h

@ -197,7 +197,7 @@ namespace storm {
qi::rule<Iterator, std::shared_ptr<storm::logic::Formula const>(), Skipper> longRunAverageRewardFormula;
qi::rule<Iterator, std::shared_ptr<storm::logic::Formula const>(), Skipper> multiFormula;
qi::rule<Iterator, std::pair<storm::solver::OptimizationDirection, storm::expressions::Variable>(), qi::locals<storm::solver::OptimizationDirection>, Skipper> quantileBoundVariable;
qi::rule<Iterator, storm::expressions::Variable(), qi::locals<boost::optional<storm::solver::OptimizationDirection>>, Skipper> quantileBoundVariable;
qi::rule<Iterator, std::shared_ptr<storm::logic::Formula const>(), Skipper> quantileFormula;
// Parser that is used to recognize doubles only (as opposed to Spirit's double_ parser).
@ -232,8 +232,8 @@ namespace storm {
std::shared_ptr<storm::logic::Formula const> createBinaryBooleanStateFormula(std::shared_ptr<storm::logic::Formula const> const& leftSubformula, std::shared_ptr<storm::logic::Formula const> const& rightSubformula, storm::logic::BinaryBooleanStateFormula::OperatorType operatorType);
std::shared_ptr<storm::logic::Formula const> createUnaryBooleanStateFormula(std::shared_ptr<storm::logic::Formula const> const& subformula, boost::optional<storm::logic::UnaryBooleanStateFormula::OperatorType> const& operatorType);
std::shared_ptr<storm::logic::Formula const> createMultiFormula(std::vector<std::shared_ptr<storm::logic::Formula const>> const& subformulas);
std::pair<storm::solver::OptimizationDirection, storm::expressions::Variable> createQuantileBoundVariables(storm::solver::OptimizationDirection const& dir, std::string const& variableName);
std::shared_ptr<storm::logic::Formula const> createQuantileFormula(std::vector<std::pair<storm::solver::OptimizationDirection, storm::expressions::Variable>> const& boundVariables, std::shared_ptr<storm::logic::Formula const> const& subformula);
storm::expressions::Variable createQuantileBoundVariables(boost::optional<storm::solver::OptimizationDirection> const& dir, std::string const& variableName);
std::shared_ptr<storm::logic::Formula const> createQuantileFormula(std::vector<storm::expressions::Variable> const& boundVariables, std::shared_ptr<storm::logic::Formula const> const& subformula);
std::set<storm::expressions::Variable> getUndefinedConstants(std::shared_ptr<storm::logic::Formula const> const& formula) const;
storm::jani::Property createProperty(boost::optional<std::string> const& propertyName, storm::modelchecker::FilterType const& filterType, std::shared_ptr<storm::logic::Formula const> const& formula, std::shared_ptr<storm::logic::Formula const> const& states);
@ -245,6 +245,7 @@ namespace storm {
uint64_t propertyCount;
std::set<storm::expressions::Variable> undefinedConstants;
std::set<storm::expressions::Variable> quantileFormulaVariables;
};
}

20
src/storm/logic/QuantileFormula.cpp

@ -7,7 +7,7 @@
namespace storm {
namespace logic {
QuantileFormula::QuantileFormula(std::vector<std::pair<storm::solver::OptimizationDirection, storm::expressions::Variable>> const& boundVariables, std::shared_ptr<Formula const> subformula) : boundVariables(boundVariables), subformula(subformula) {
QuantileFormula::QuantileFormula(std::vector<storm::expressions::Variable> const& boundVariables, std::shared_ptr<Formula const> subformula) : boundVariables(boundVariables), subformula(subformula) {
STORM_LOG_THROW(!boundVariables.empty(), storm::exceptions::InvalidArgumentException, "Quantile formula without bound variables are invalid.");
}
@ -45,28 +45,18 @@ namespace storm {
storm::expressions::Variable const& QuantileFormula::getBoundVariable() const {
STORM_LOG_THROW(boundVariables.size() == 1, storm::exceptions::InvalidArgumentException, "Requested unique bound variables. However, there are multiple bound variables defined.");
return boundVariables.front().second;
return boundVariables.front();
}
storm::expressions::Variable const& QuantileFormula::getBoundVariable(uint64_t index) const {
STORM_LOG_THROW(index < boundVariables.size(), storm::exceptions::InvalidArgumentException, "Requested bound variable with index" << index << ". However, there are only " << boundVariables.size() << " bound variables.");
return boundVariables[index].second;
return boundVariables[index];
}
std::vector<std::pair<storm::solver::OptimizationDirection, storm::expressions::Variable>> const& QuantileFormula::getBoundVariables() const {
std::vector<storm::expressions::Variable> const& QuantileFormula::getBoundVariables() const {
return boundVariables;
}
storm::solver::OptimizationDirection const& QuantileFormula::getOptimizationDirection() const {
STORM_LOG_THROW(boundVariables.size() == 1, storm::exceptions::InvalidArgumentException, "Requested unique optimization direction of the bound variables. However, there are multiple bound variables defined.");
return boundVariables.front().first;
}
storm::solver::OptimizationDirection const& QuantileFormula::getOptimizationDirection(uint64_t index) const {
STORM_LOG_THROW(index < boundVariables.size(), storm::exceptions::InvalidArgumentException, "Requested optimization direction with index" << index << ". However, there are only " << boundVariables.size() << " bound variables.");
return boundVariables[index].first;
}
boost::any QuantileFormula::accept(FormulaVisitor const& visitor, boost::any const& data) const {
return visitor.visit(*this, data);
}
@ -86,7 +76,7 @@ namespace storm {
std::ostream& QuantileFormula::writeToStream(std::ostream& out) const {
out << "quantile(";
for (auto const& bv : boundVariables) {
out << (storm::solver::minimize(bv.first) ? "min" : "max") << " " << bv.second.getName() << ", ";
out << bv.getName() << ", ";
}
subformula->writeToStream(out);
out << ")";

9
src/storm/logic/QuantileFormula.h

@ -1,13 +1,12 @@
#pragma once
#include "storm/logic/Formula.h"
#include "storm/solver/OptimizationDirection.h"
namespace storm {
namespace logic {
class QuantileFormula : public Formula {
public:
QuantileFormula(std::vector<std::pair<storm::solver::OptimizationDirection, storm::expressions::Variable>> const& boundVariables, std::shared_ptr<Formula const> subformula);
QuantileFormula(std::vector<storm::expressions::Variable> const& boundVariables, std::shared_ptr<Formula const> subformula);
virtual ~QuantileFormula();
@ -23,9 +22,7 @@ namespace storm {
storm::expressions::Variable const& getBoundVariable() const;
storm::expressions::Variable const& getBoundVariable(uint64_t index) const;
std::vector<std::pair<storm::solver::OptimizationDirection, storm::expressions::Variable>> const& getBoundVariables() const;
storm::solver::OptimizationDirection const& getOptimizationDirection() const;
storm::solver::OptimizationDirection const& getOptimizationDirection(uint64_t index) const;
std::vector<storm::expressions::Variable> const& getBoundVariables() const;
virtual boost::any accept(FormulaVisitor const& visitor, boost::any const& data) const override;
virtual void gatherAtomicExpressionFormulas(std::vector<std::shared_ptr<AtomicExpressionFormula const>>& atomicExpressionFormulas) const override;
@ -34,7 +31,7 @@ namespace storm {
virtual std::ostream& writeToStream(std::ostream& out) const override;
private:
std::vector<std::pair<storm::solver::OptimizationDirection, storm::expressions::Variable>> boundVariables;
std::vector<storm::expressions::Variable> boundVariables;
std::shared_ptr<Formula const> subformula;
};
}

32
src/storm/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp

@ -7,9 +7,11 @@
#include "storm/modelchecker/results/ExplicitQualitativeCheckResult.h"
#include "storm/modelchecker/results/ExplicitQuantitativeCheckResult.h"
#include "storm/modelchecker/results/ExplicitParetoCurveCheckResult.h"
#include "storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.h"
#include "storm/modelchecker/csl/helper/SparseCtmcCslHelper.h"
#include "storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h"
#include "storm/logic/FragmentSpecification.h"
@ -31,7 +33,13 @@ namespace storm {
template<typename SparseDtmcModelType>
bool SparseDtmcPrctlModelChecker<SparseDtmcModelType>::canHandle(CheckTask<storm::logic::Formula, ValueType> const& checkTask) const {
storm::logic::Formula const& formula = checkTask.getFormula();
return formula.isInFragment(storm::logic::prctl().setLongRunAverageRewardFormulasAllowed(true).setLongRunAverageProbabilitiesAllowed(true).setConditionalProbabilityFormulasAllowed(true).setConditionalRewardFormulasAllowed(true).setTotalRewardFormulasAllowed(true).setOnlyEventuallyFormuluasInConditionalFormulasAllowed(true).setRewardBoundedUntilFormulasAllowed(true).setRewardBoundedCumulativeRewardFormulasAllowed(true).setMultiDimensionalBoundedUntilFormulasAllowed(true).setMultiDimensionalCumulativeRewardFormulasAllowed(true).setTimeOperatorsAllowed(true).setReachbilityTimeFormulasAllowed(true));
if (formula.isInFragment(storm::logic::prctl().setLongRunAverageRewardFormulasAllowed(true).setLongRunAverageProbabilitiesAllowed(true).setConditionalProbabilityFormulasAllowed(true).setConditionalRewardFormulasAllowed(true).setTotalRewardFormulasAllowed(true).setOnlyEventuallyFormuluasInConditionalFormulasAllowed(true).setRewardBoundedUntilFormulasAllowed(true).setRewardBoundedCumulativeRewardFormulasAllowed(true).setMultiDimensionalBoundedUntilFormulasAllowed(true).setMultiDimensionalCumulativeRewardFormulasAllowed(true).setTimeOperatorsAllowed(true).setReachbilityTimeFormulasAllowed(true))) {
return true;
} else if (formula.isInFragment(storm::logic::quantiles())) {
if (this->getModel().getInitialStates().getNumberOfSetBits() > 1) return false;
return true;
}
return false;
}
template<typename SparseDtmcModelType>
@ -184,6 +192,28 @@ namespace storm {
return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(std::move(numericResult)));
}
template<>
std::unique_ptr<CheckResult> SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<storm::RationalFunction>>::checkQuantileFormula(Environment const& env, CheckTask<storm::logic::QuantileFormula, ValueType> const& checkTask) {
STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "Quantiles for parametric models are not implemented.");
}
template<typename SparseDtmcModelType>
std::unique_ptr<CheckResult> SparseDtmcPrctlModelChecker<SparseDtmcModelType>::checkQuantileFormula(Environment const& env, CheckTask<storm::logic::QuantileFormula, ValueType> const& checkTask) {
STORM_LOG_THROW(checkTask.isOnlyInitialStatesRelevantSet(), storm::exceptions::InvalidOperationException, "Computing quantiles is only supported for the initial states of a model.");
STORM_LOG_THROW(this->getModel().getInitialStates().getNumberOfSetBits() == 1, storm::exceptions::InvalidOperationException, "Quantiles not supported on models with multiple initial states.");
uint64_t initialState = *this->getModel().getInitialStates().begin();
helper::rewardbounded::QuantileHelper<SparseDtmcModelType> qHelper(this->getModel(), checkTask.getFormula());
auto res = qHelper.computeQuantile(env);
if (res.size() == 1 && res.front().size() == 1) {
return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(initialState, std::move(res.front().front())));
} else {
return std::unique_ptr<CheckResult>(new ExplicitParetoCurveCheckResult<ValueType>(initialState, std::move(res)));
}
}
template class SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>>;
#ifdef STORM_HAVE_CARL

1
src/storm/modelchecker/prctl/SparseDtmcPrctlModelChecker.h

@ -34,6 +34,7 @@ namespace storm {
virtual std::unique_ptr<CheckResult> computeConditionalRewards(Environment const& env, storm::logic::RewardMeasureType rewardMeasureType, CheckTask<storm::logic::ConditionalFormula, ValueType> const& checkTask) override;
virtual std::unique_ptr<CheckResult> computeLongRunAverageRewards(Environment const& env, storm::logic::RewardMeasureType rewardMeasureType, CheckTask<storm::logic::LongRunAverageRewardFormula, ValueType> const& checkTask) override;
virtual std::unique_ptr<CheckResult> computeReachabilityTimes(Environment const& env, storm::logic::RewardMeasureType rewardMeasureType, CheckTask<storm::logic::EventuallyFormula, ValueType> const& checkTask) override;
virtual std::unique_ptr<CheckResult> checkQuantileFormula(Environment const& env, CheckTask<storm::logic::QuantileFormula, ValueType> const& checkTask) override;
};
} // namespace modelchecker

2
src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp

@ -228,7 +228,7 @@ namespace storm {
template<typename SparseMdpModelType>
std::unique_ptr<CheckResult> SparseMdpPrctlModelChecker<SparseMdpModelType>::checkQuantileFormula(Environment const& env, CheckTask<storm::logic::QuantileFormula, ValueType> const& checkTask) {
STORM_LOG_THROW(checkTask.isOnlyInitialStatesRelevantSet(), storm::exceptions::InvalidOperationException, "Computing quantiles is only supported for models with a single initial states.");
STORM_LOG_THROW(checkTask.isOnlyInitialStatesRelevantSet(), storm::exceptions::InvalidOperationException, "Computing quantiles is only supported for the initial states of a model.");
STORM_LOG_THROW(this->getModel().getInitialStates().getNumberOfSetBits() == 1, storm::exceptions::InvalidOperationException, "Quantiles not supported on models with multiple initial states.");
uint64_t initialState = *this->getModel().getInitialStates().begin();

186
src/storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.cpp

@ -0,0 +1,186 @@
#include <storm/utility/ExpressionHelper.h>
#include "storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.h"
#include "storm/utility/macros.h"
#include "storm/exceptions/IllegalArgumentException.h"
#include "storm/utility/solver.h"
#include "storm/solver/SmtSolver.h"
namespace storm {
namespace modelchecker {
namespace helper {
namespace rewardbounded {
CostLimit::CostLimit(uint64_t const &costLimit) : value(costLimit) {
// Intentionally left empty
}
bool CostLimit::isInfinity() const {
return value == std::numeric_limits<uint64_t>::max();
}
uint64_t const& CostLimit::get() const {
STORM_LOG_ASSERT(!isInfinity(), "Tried to get an infinite cost limit as int.");
return value;
}
uint64_t& CostLimit::get() {
STORM_LOG_ASSERT(!isInfinity(), "Tried to get an infinite cost limit as int.");
return value;
}
bool CostLimit::operator<(CostLimit const& other) const {
// Since infinity is represented by the max integer, we can compare this way.
return value < other.value;
}
bool CostLimit::operator==(CostLimit const& other) const {
return value == other.value;
}
CostLimit CostLimit::infinity() {
return CostLimit(std::numeric_limits<uint64_t>::max());
}
bool CostLimitClosure::CostLimitsCompare::operator()(storm::modelchecker::helper::rewardbounded::CostLimits const& lhs, storm::modelchecker::helper::rewardbounded::CostLimits const& rhs) const {
for (uint64_t i = 0; i < lhs.size(); ++i) {
if (lhs[i] < rhs[i]) {
return true;
} else if (rhs[i] < lhs[i]) {
return false;
}
}
return false;
}
CostLimitClosure::CostLimitClosure(storm::storage::BitVector const &downwardDimensions)
: downwardDimensions(downwardDimensions) {
// Intentionally left empty
}
bool CostLimitClosure::insert(CostLimits const& costLimits) {
// Iterate over all points in the generator and check whether they dominate the given point or vice versa
// TODO: make this more efficient by exploiting the order of the generator set.
std::vector<CostLimits> pointsToErase;
for (auto const& b : generator) {
if (dominates(b, costLimits)) {
// The given point is already contained in this closure.
// Since domination is transitive, this should not happen:
STORM_LOG_ASSERT(pointsToErase.empty(), "Inconsistent generator of CostLimitClosure.");
return false;
}
if (dominates(costLimits, b)) {
// b will be dominated by the new point so we erase it later.
// Note that b != newPoint holds if we reach this point
pointsToErase.push_back(b);
}
}
for (auto const& b : pointsToErase) {
generator.erase(b);
}
generator.insert(std::move(costLimits));
return true;
}
bool CostLimitClosure::contains(CostLimits const& costLimits) const {
// Iterate over all points in the generator and check whether they dominate the given point.
// TODO: make this more efficient by exploiting the order of the generator set.
for (auto const&b : generator) {
if (dominates(b,costLimits)) {
return true;
}
}
return false;
}
bool CostLimitClosure::containsUpwardClosure(CostLimits const& costLimits) const {
CostLimits infinityProjection(costLimits);
for (auto const& dim : downwardDimensions) {
infinityProjection[dim] = CostLimit::infinity();
}
return contains(infinityProjection);
}
bool CostLimitClosure::empty() const {
return generator.empty();
}
bool CostLimitClosure::full() const {
CostLimits p(dimension(), CostLimit(0));
for (auto const& dim : downwardDimensions) {
p[dim] = CostLimit::infinity();
}
return contains(p);
}
bool CostLimitClosure::dominates(CostLimits const& lhs, CostLimits const& rhs) const {
for (uint64_t i = 0; i < lhs.size(); ++i) {
if (downwardDimensions.get(i)) {
if (lhs[i] < rhs[i]) {
return false;
}
} else {
if (rhs[i] < lhs[i]) {
return false;
}
}
}
return true;
}
std::vector<CostLimits> CostLimitClosure::getDominatingCostLimits(CostLimits const& costLimits) const {
std::vector<CostLimits> result;
for (auto const &b : generator) {
if (dominates(b, costLimits)) {
result.push_back(b);
}
}
return result;
}
typename CostLimitClosure::GeneratorType const &CostLimitClosure::getGenerator() const {
return generator;
}
uint64_t CostLimitClosure::dimension() const {
return downwardDimensions.size();
}
bool CostLimitClosure::unionFull(CostLimitClosure const& first, CostLimitClosure const& second) {
assert(first.dimension() == second.dimension());
uint64_t dimension = first.dimension();
auto manager = std::make_shared<storm::expressions::ExpressionManager>();
auto solver = storm::utility::solver::getSmtSolver(*manager);
std::vector<storm::expressions::Expression> point;
storm::expressions::Expression zero = manager->integer(0);
for (uint64_t i = 0; i < dimension; ++i) {
point.push_back(manager->declareIntegerVariable("x" + std::to_string(i)).getExpression());
solver->add(point.back() >= zero);
}
for (auto const& cl : {first, second}) {
for (auto const& q : cl.getGenerator()) {
storm::expressions::Expression pointNotDominated;
for (uint64_t i = 0; i < point.size(); ++i) {
if (!cl.downwardDimensions.get(i) || !q[i].isInfinity()) {
assert(!q[i].isInfinity());
storm::expressions::Expression qi = manager->integer(q[i].get());
storm::expressions::Expression piNotDominated = cl.downwardDimensions.get(i) ? point[i] > qi : point[i] < qi;
if (piNotDominated.isInitialized()) {
pointNotDominated = pointNotDominated || piNotDominated;
} else {
pointNotDominated = piNotDominated;
}
}
}
if (pointNotDominated.isInitialized()) {
solver->add(pointNotDominated);
} else {
solver->add(manager->boolean(false));
}
}
}
return solver->check() == storm::solver::SmtSolver::CheckResult::Unsat;;
}
}
}
}
}

60
src/storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.h

@ -0,0 +1,60 @@
#pragma once
#include <vector>
#include <set>
#include <string>
#include "storm/storage/BitVector.h"
namespace storm {
namespace modelchecker {
namespace helper {
namespace rewardbounded {
class CostLimit {
public:
CostLimit(uint64_t const& costLimit);
bool isInfinity() const;
uint64_t const& get() const;
uint64_t& get();
bool operator<(CostLimit const& other) const;
bool operator==(CostLimit const& other) const;
static CostLimit infinity();
private:
uint64_t value;
};
typedef std::vector<CostLimit> CostLimits;
class CostLimitClosure {
public:
struct CostLimitsCompare {
bool operator() (CostLimits const& lhs, CostLimits const& rhs) const;
};
typedef std::set<CostLimits, CostLimitsCompare> GeneratorType;
explicit CostLimitClosure(storm::storage::BitVector const& downwardDimensions);
bool insert(CostLimits const& costLimits);
bool contains(CostLimits const& costLimits) const;
bool containsUpwardClosure(CostLimits const& costLimits) const;
bool dominates(CostLimits const& lhs, CostLimits const& rhs) const;
bool empty() const; // True if there is no point in this
bool full() const; // True if all points lie in this.
std::vector<CostLimits> getDominatingCostLimits(CostLimits const& costLimits) const;
GeneratorType const& getGenerator() const;
uint64_t dimension() const;
/*!
* Returns true if the union of the two closures is full, i.e., contains every point.
*/
static bool unionFull(CostLimitClosure const& first, CostLimitClosure const& second);
private:
/// The dimensions that are downwards closed, i.e., if x is in the closure, then also all y <= x
storm::storage::BitVector downwardDimensions;
GeneratorType generator;
};
}
}
}
}

11
src/storm/modelchecker/prctl/helper/rewardbounded/EpochManager.cpp

@ -210,6 +210,17 @@ namespace storm {
return (epoch >> (dimension * bitsPerDimension)) & dimensionBitMask;
}
uint64_t EpochManager::getSumOfDimensions(Epoch const& epoch) const {
STORM_LOG_ASSERT(dimensionCount > 0, "Invoked EpochManager with zero dimension count.");
uint64_t sumOfDimensions = 0;
for (uint64_t dim = 0; dim < getDimensionCount(); ++dim) {
if (!isBottomDimension(epoch, dim)) {
sumOfDimensions += getDimensionOfEpoch(epoch, dim) + 1;
}
}
return sumOfDimensions;
}
std::string EpochManager::toString(Epoch const& epoch) const {
STORM_LOG_ASSERT(dimensionCount > 0, "Invoked EpochManager with zero dimension count.");
std::string res = "<" + (isBottomDimension(epoch, 0) ? "_" : std::to_string(getDimensionOfEpoch(epoch, 0)));

1
src/storm/modelchecker/prctl/helper/rewardbounded/EpochManager.h

@ -45,6 +45,7 @@ namespace storm {
bool isBottomDimension(Epoch const& epoch, uint64_t const& dimension) const;
bool isBottomDimensionEpochClass(EpochClass const& epochClass, uint64_t const& dimension) const;
uint64_t getDimensionOfEpoch(Epoch const& epoch, uint64_t const& dimension) const; // assumes that the dimension is not bottom
uint64_t getSumOfDimensions(Epoch const& epoch) const; // assumes that the dimension is not bottom
std::string toString(Epoch const& epoch) const;

96
src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp

@ -64,8 +64,6 @@ namespace storm {
template<typename ValueType, bool SingleObjectiveMode>
void MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::initialize(std::set<storm::expressions::Variable> const& infinityBoundVariables) {
maxSolutionsStored = 0;
STORM_LOG_ASSERT(!SingleObjectiveMode || (this->objectives.size() == 1), "Enabled single objective mode but there are multiple objectives.");
std::vector<Epoch> epochSteps;
initializeObjectives(epochSteps, infinityBoundVariables);
@ -278,7 +276,8 @@ namespace storm {
STORM_LOG_THROW(SingleObjectiveMode, storm::exceptions::NotSupportedException, "Letting lower bounds approach infinity is only supported in single objective mode."); // It most likely also works with multiple objectives with the same shape. However, we haven't checked this yet.
STORM_LOG_THROW(objectives.front().formula->isProbabilityOperatorFormula(), storm::exceptions::NotSupportedException, "Letting lower bounds approach infinity is only supported for probability operator formulas");
auto const& probabilityOperatorFormula = objectives.front().formula->asProbabilityOperatorFormula();
STORM_LOG_THROW(probabilityOperatorFormula.getSubformula().isBoundedUntilFormula() && probabilityOperatorFormula.hasOptimalityType() && storm::solver::maximize(probabilityOperatorFormula.getOptimalityType()), storm::exceptions::NotSupportedException, "Letting lower bounds approach infinity is only supported for maximizing bounded until probabilities.");
STORM_LOG_THROW(probabilityOperatorFormula.getSubformula().isBoundedUntilFormula(), storm::exceptions::NotSupportedException, "Letting lower bounds approach infinity is only supported for bounded until probabilities.");
STORM_LOG_THROW(!model.isNondeterministicModel() || (probabilityOperatorFormula.hasOptimalityType() && storm::solver::maximize(probabilityOperatorFormula.getOptimalityType())), storm::exceptions::NotSupportedException, "Letting lower bounds approach infinity is only supported for maximizing bounded until probabilities.");
STORM_LOG_THROW(upperBoundedDimensions.empty() || !probabilityOperatorFormula.getSubformula().asBoundedUntilFormula().hasMultiDimensionalSubformulas(), storm::exceptions::NotSupportedException, "Letting lower bounds approach infinity is only supported if the formula has either only lower bounds or if it has a single goal state."); // This would fail because the upper bounded dimension(s) might be satisfied already. One should erase epoch steps in the epoch model (after applying the goal-unfolding).
storm::storage::BitVector choicesWithoutUpperBoundedStep(model.getNumberOfChoices(), true);
@ -332,48 +331,33 @@ namespace storm {
}
template<typename ValueType, bool SingleObjectiveMode>
std::vector<typename MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::Epoch> MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::getEpochComputationOrder(Epoch const& startEpoch) {
std::vector<typename MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::Epoch> MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::getEpochComputationOrder(Epoch const& startEpoch, bool stopAtComputedEpochs) {
// Perform a DFS to find all the reachable epochs
std::vector<Epoch> dfsStack;
std::set<Epoch, std::function<bool(Epoch const&, Epoch const&)>> collectedEpochs(std::bind(&EpochManager::epochClassZigZagOrder, &epochManager, std::placeholders::_1, std::placeholders::_2));
collectedEpochs.insert(startEpoch);
dfsStack.push_back(startEpoch);
if (!stopAtComputedEpochs || epochSolutions.count(startEpoch) == 0) {
collectedEpochs.insert(startEpoch);
dfsStack.push_back(startEpoch);
}
while (!dfsStack.empty()) {
Epoch currentEpoch = dfsStack.back();
dfsStack.pop_back();
for (auto const& step : possibleEpochSteps) {
Epoch successorEpoch = epochManager.getSuccessorEpoch(currentEpoch, step);
/*
for (auto const& e : collectedEpochs) {
std::cout << "Comparing " << epochManager.toString(e) << " and " << epochManager.toString(successorEpoch) << std::endl;
if (epochManager.epochClassZigZagOrder(e, successorEpoch)) {
std::cout << " " << epochManager.toString(e) << " < " << epochManager.toString(successorEpoch) << std::endl;
if (!stopAtComputedEpochs || epochSolutions.count(successorEpoch) == 0) {
if (collectedEpochs.insert(successorEpoch).second) {
dfsStack.push_back(std::move(successorEpoch));
}
if (epochManager.epochClassZigZagOrder(successorEpoch, e)) {
std::cout << " " << epochManager.toString(e) << " > " << epochManager.toString(successorEpoch) << std::endl;
}
}
*/
if (collectedEpochs.insert(successorEpoch).second) {
dfsStack.push_back(std::move(successorEpoch));
}
}
}
/*
std::cout << "Resulting order: ";
for (auto const& e : collectedEpochs) {
std::cout << epochManager.toString(e) << ", ";
}
std::cout << std::endl;
*/
return std::vector<Epoch>(collectedEpochs.begin(), collectedEpochs.end());
}
template<typename ValueType, bool SingleObjectiveMode>
EpochModel<ValueType, SingleObjectiveMode>& MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::setCurrentEpoch(Epoch const& epoch) {
STORM_LOG_DEBUG("Setting model for epoch " << epochManager.toString(epoch));
STORM_LOG_THROW(getProb1Objectives().empty(), storm::exceptions::InvalidOperationException, "The objective " << objectives[*getProb1Objectives().begin()].formula << " is satisfied already in the initial state(s). This special case can not be handled by the reward unfolding."); // This case should have been handled 'from the outside'
// Check if we need to update the current epoch class
if (!currentEpoch || !epochManager.compareEpochClass(epoch, currentEpoch.get())) {
@ -629,8 +613,6 @@ namespace storm {
epochModel.objectiveRewardFilter.push_back(storm::utility::vector::filterZero(objRewards));
epochModel.objectiveRewardFilter.back().complement();
}
epochModelSizes.push_back(epochModel.epochMatrix.getRowGroupCount());
}
@ -669,6 +651,20 @@ namespace storm {
storm::utility::vector::addScaledVector(solution, solutionToAdd, scalingFactor);
}
template<typename ValueType, bool SingleObjectiveMode>
template<bool SO, typename std::enable_if<SO, int>::type>
void MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::setSolutionEntry(SolutionType& solution, uint64_t objIndex, ValueType const& value) const {
STORM_LOG_ASSERT(objIndex == 0, "Invalid objective index in single objective mode");
solution = value;
}
template<typename ValueType, bool SingleObjectiveMode>
template<bool SO, typename std::enable_if<!SO, int>::type>
void MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::setSolutionEntry(SolutionType& solution, uint64_t objIndex, ValueType const& value) const {
STORM_LOG_ASSERT(objIndex < solution.size(), "Invalid objective index " << objIndex << ".");
solution[objIndex] = value;
}
template<typename ValueType, bool SingleObjectiveMode>
template<bool SO, typename std::enable_if<SO, int>::type>
std::string MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::solutionToString(SolutionType const& solution) const {
@ -697,11 +693,7 @@ namespace storm {
template<typename ValueType, bool SingleObjectiveMode>
ValueType MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::getRequiredEpochModelPrecision(Epoch const& startEpoch, ValueType const& precision) {
uint64_t sumOfDimensions = 0;
for (uint64_t dim = 0; dim < epochManager.getDimensionCount(); ++dim) {
sumOfDimensions += epochManager.getDimensionOfEpoch(startEpoch, dim) + 1;
}
return precision / storm::utility::convertNumber<ValueType>(sumOfDimensions);
return precision / storm::utility::convertNumber<ValueType>(epochManager.getSumOfDimensions(startEpoch) + 1);
}
template<typename ValueType, bool SingleObjectiveMode>
@ -836,9 +828,6 @@ namespace storm {
solution.productStateToSolutionVectorMap = productStateToEpochModelInStateMap;
solution.solutions = std::move(inStateSolutions);
epochSolutions[currentEpoch.get()] = std::move(solution);
maxSolutionsStored = std::max((uint64_t) epochSolutions.size(), maxSolutionsStored);
}
template<typename ValueType, bool SingleObjectiveMode>
@ -858,21 +847,39 @@ namespace storm {
template<typename ValueType, bool SingleObjectiveMode>
typename MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::SolutionType const& MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::getStateSolution(EpochSolution const& epochSolution, uint64_t const& productState) {
STORM_LOG_ASSERT(productState < epochSolution.productStateToSolutionVectorMap->size(), "Requested solution at an unexisting product state.");
STORM_LOG_ASSERT((*epochSolution.productStateToSolutionVectorMap)[productState] < epochSolution.solutions.size(), "Requested solution for epoch at a state for which no solution was stored.");
STORM_LOG_ASSERT(getProb1Objectives().empty(), "One of the objectives is already satisfied in the initial state. This special case is not handled."); // This case should have been handled 'from the outside'
STORM_LOG_ASSERT((*epochSolution.productStateToSolutionVectorMap)[productState] < epochSolution.solutions.size(), "Requested solution for epoch at product state " << productState << " for which no solution was stored.");
return epochSolution.solutions[(*epochSolution.productStateToSolutionVectorMap)[productState]];
}
template<typename ValueType, bool SingleObjectiveMode>
typename MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::SolutionType const& MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::getInitialStateResult(Epoch const& epoch) {
typename MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::SolutionType MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::getInitialStateResult(Epoch const& epoch) {
STORM_LOG_ASSERT(model.getInitialStates().getNumberOfSetBits() == 1, "The model has multiple initial states.");
return getStateSolution(epoch, productModel->getInitialProductState(*model.getInitialStates().begin(), model.getInitialStates()));
return getInitialStateResult(epoch, *model.getInitialStates().begin());
}
template<typename ValueType, bool SingleObjectiveMode>
typename MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::SolutionType const& MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::getInitialStateResult(Epoch const& epoch, uint64_t initialStateIndex) {
typename MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::SolutionType MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::getInitialStateResult(Epoch const& epoch, uint64_t initialStateIndex) {
STORM_LOG_ASSERT(model.getInitialStates().get(initialStateIndex), "The given model state is not an initial state.");
return getStateSolution(epoch, productModel->getInitialProductState(initialStateIndex, model.getInitialStates()));
auto result = getStateSolution(epoch, productModel->getInitialProductState(initialStateIndex, model.getInitialStates(), epochManager.getEpochClass(epoch)));
for (uint64_t objIndex = 0; objIndex < objectives.size(); ++objIndex) {
if (productModel->getProb1InitialStates(objIndex) && productModel->getProb1InitialStates(objIndex)->get(initialStateIndex)) {
// Check whether the objective can actually hold in this epoch
bool objectiveHolds = true;
for (auto const& dim : objectiveDimensions[objIndex]) {
if (dimensions[dim].boundType == DimensionBoundType::LowerBound && !epochManager.isBottomDimension(epoch, dim)) {
objectiveHolds = false;
} else if (dimensions[dim].boundType == DimensionBoundType::UpperBound && epochManager.isBottomDimension(epoch, dim)) {
objectiveHolds = false;
}
STORM_LOG_ASSERT(dimensions[dim].boundType != DimensionBoundType::LowerBoundInfinity, "Unexpected bound type at this point.");
}
if (objectiveHolds) {
setSolutionEntry(result, objIndex, storm::utility::one<ValueType>());
}
}
}
return result;
}
template<typename ValueType, bool SingleObjectiveMode>
@ -885,11 +892,6 @@ namespace storm {
return dimensions.at(dim);
}
template<typename ValueType, bool SingleObjectiveMode>
storm::storage::BitVector const& MultiDimensionalRewardUnfolding<ValueType, SingleObjectiveMode>::getProb1Objectives() const {
return productModel->getProb1Objectives();
}
template class MultiDimensionalRewardUnfolding<double, true>;
template class MultiDimensionalRewardUnfolding<double, false>;
template class MultiDimensionalRewardUnfolding<storm::RationalNumber, true>;

29
src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h

@ -54,7 +54,11 @@ namespace storm {
*/
Epoch getStartEpoch(bool setUnknownDimsToBottom = false);
std::vector<Epoch> getEpochComputationOrder(Epoch const& startEpoch);
/*!
* Computes a sequence of epochs that need to be analyzed to get a result at the start epoch.
* @param stopAtComputedEpochs if set, the search for epochs that need to be computed is stopped at epochs that already have been computed earlier.
*/
std::vector<Epoch> getEpochComputationOrder(Epoch const& startEpoch, bool stopAtComputedEpochs = false);
EpochModel<ValueType, SingleObjectiveMode>& setCurrentEpoch(Epoch const& epoch);
@ -72,18 +76,12 @@ namespace storm {
boost::optional<ValueType> getLowerObjectiveBound(uint64_t objectiveIndex = 0);
void setSolutionForCurrentEpoch(std::vector<SolutionType>&& inStateSolutions);
SolutionType const& getInitialStateResult(Epoch const& epoch); // Assumes that the initial state is unique
SolutionType const& getInitialStateResult(Epoch const& epoch, uint64_t initialStateIndex);
SolutionType getInitialStateResult(Epoch const& epoch); // Assumes that the initial state is unique
SolutionType getInitialStateResult(Epoch const& epoch, uint64_t initialStateIndex);
EpochManager const& getEpochManager() const;
Dimension<ValueType> const& getDimension(uint64_t dim) const;
/*!
* Returns objectives that are always satisfied (i.e., have probability one) in all initial states.
* These objectives can not be handled by this as they can not be translated into expected rewards.
*/
storm::storage::BitVector const& getProb1Objectives() const;
private:
void setCurrentEpochClass(Epoch const& epoch);
@ -105,6 +103,14 @@ namespace storm {
template<bool SO = SingleObjectiveMode, typename std::enable_if<!SO, int>::type = 0>
void addScaledSolution(SolutionType& solution, SolutionType const& solutionToAdd, ValueType const& scalingFactor) const;
template<bool SO = SingleObjectiveMode, typename std::enable_if<SO, int>::type = 0>
void setSolutionEntry(SolutionType& solution, uint64_t objIndex, ValueType const& value) const;
template<bool SO = SingleObjectiveMode, typename std::enable_if<!SO, int>::type = 0>
void setSolutionEntry(SolutionType& solution, uint64_t objIndex, ValueType const& value) const;
template<bool SO = SingleObjectiveMode, typename std::enable_if<SO, int>::type = 0>
std::string solutionToString(SolutionType const& solution) const;
template<bool SO = SingleObjectiveMode, typename std::enable_if<!SO, int>::type = 0>
@ -136,11 +142,6 @@ namespace storm {
std::vector<Dimension<ValueType>> dimensions;
std::vector<storm::storage::BitVector> objectiveDimensions;
storm::utility::Stopwatch swInit, swFindSol, swInsertSol, swSetEpoch, swSetEpochClass, swAux1, swAux2, swAux3, swAux4;
std::vector<uint64_t> epochModelSizes;
uint64_t maxSolutionsStored;
};
}
}

17
src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.cpp

@ -21,7 +21,7 @@ namespace storm {
namespace rewardbounded {
template<typename ValueType>
ProductModel<ValueType>::ProductModel(storm::models::sparse::Model<ValueType> const& model, std::vector<storm::modelchecker::multiobjective::Objective<ValueType>> const& objectives, std::vector<Dimension<ValueType>> const& dimensions, std::vector<storm::storage::BitVector> const& objectiveDimensions, EpochManager const& epochManager, std::vector<Epoch> const& originalModelSteps) : dimensions(dimensions), objectiveDimensions(objectiveDimensions), epochManager(epochManager), memoryStateManager(dimensions.size()), prob1Objectives(objectives.size(), false) {
ProductModel<ValueType>::ProductModel(storm::models::sparse::Model<ValueType> const& model, std::vector<storm::modelchecker::multiobjective::Objective<ValueType>> const& objectives, std::vector<Dimension<ValueType>> const& dimensions, std::vector<storm::storage::BitVector> const& objectiveDimensions, EpochManager const& epochManager, std::vector<Epoch> const& originalModelSteps) : dimensions(dimensions), objectiveDimensions(objectiveDimensions), epochManager(epochManager), memoryStateManager(dimensions.size()), prob1InitialStates(objectives.size(), boost::none) {
for (uint64_t dim = 0; dim < dimensions.size(); ++dim) {
if (!dimensions[dim].memoryLabel) {
@ -167,10 +167,9 @@ namespace storm {
if (memStateBV.full()) {
storm::storage::BitVector initialTransitionStates = model.getInitialStates() & transitionStates;
// At this point we can check whether there is an initial state that already satisfies all subObjectives.
// Such a situation can not be reduced (easily) to an expected reward computation.
if (memStatePrimeBV.empty() && !initialTransitionStates.empty() && !objectiveContainsLowerBound && !initialTransitionStates.isDisjointFrom(constraintStates)) {
STORM_LOG_THROW(model.getInitialStates() == initialTransitionStates, storm::exceptions::NotSupportedException, "The objective " << *objectives[objIndex].formula << " is already satisfied in an initial state. This special case is not supported.");
prob1Objectives.set(objIndex, true);
// Such a situation can not be reduced (easily) to an expected reward computation and thus requires special treatment
if (memStatePrimeBV.empty() && !initialTransitionStates.empty()) {
prob1InitialStates[objIndex] = initialTransitionStates;
}
for (auto const& initState : initialTransitionStates) {
@ -331,11 +330,11 @@ namespace storm {
}
template<typename ValueType>
uint64_t ProductModel<ValueType>::getInitialProductState(uint64_t const& initialModelState, storm::storage::BitVector const& initialModelStates) const {
uint64_t ProductModel<ValueType>::getInitialProductState(uint64_t const& initialModelState, storm::storage::BitVector const& initialModelStates, EpochClass const& epochClass) const {
auto productInitStateIt = getProduct().getInitialStates().begin();
productInitStateIt += initialModelStates.getNumberOfSetBitsBeforeIndex(initialModelState);
STORM_LOG_ASSERT(getModelState(*productInitStateIt) == initialModelState, "Could not find the corresponding initial state in the product model.");
return transformProductState(*productInitStateIt, epochManager.getEpochClass(epochManager.getZeroEpoch()), memoryStateManager.getInitialMemoryState());
return transformProductState(*productInitStateIt, epochClass, memoryStateManager.getInitialMemoryState());
}
template<typename ValueType>
@ -649,8 +648,8 @@ namespace storm {
}
template<typename ValueType>
storm::storage::BitVector const& ProductModel<ValueType>::getProb1Objectives() {
return prob1Objectives;
boost::optional<storm::storage::BitVector> const& ProductModel<ValueType>::getProb1InitialStates(uint64_t objectiveIndex) const {
return prob1InitialStates[objectiveIndex];
}
template class ProductModel<double>;

8
src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.h

@ -33,7 +33,7 @@ namespace storm {
bool productStateExists(uint64_t const& modelState, uint64_t const& memoryState) const;
uint64_t getProductState(uint64_t const& modelState, uint64_t const& memoryState) const;
uint64_t getInitialProductState(uint64_t const& initialModelState, storm::storage::BitVector const& initialModelStates) const;
uint64_t getInitialProductState(uint64_t const& initialModelState, storm::storage::BitVector const& initialModelStates, EpochClass const& epochClass) const;
uint64_t getModelState(uint64_t const& productState) const;
MemoryState getMemoryState(uint64_t const& productState) const;
MemoryStateManager const& getMemoryStateManager() const;
@ -47,8 +47,8 @@ namespace storm {
MemoryState transformMemoryState(MemoryState const& memoryState, EpochClass const& epochClass, MemoryState const& predecessorMemoryState) const;
uint64_t transformProductState(uint64_t const& productState, EpochClass const& epochClass, MemoryState const& predecessorMemoryState) const;
// returns objectives that have probability one, already in the initial state.
storm::storage::BitVector const& getProb1Objectives();
/// returns the initial states (with respect to the original model) that already satisfy the given objective with probability one, assuming that the cost bounds at the current epoch allow for the objective to be satisfied.
boost::optional<storm::storage::BitVector> const& getProb1InitialStates(uint64_t objectiveIndex) const;
private:
@ -77,7 +77,7 @@ namespace storm {
std::vector<uint64_t> productToModelStateMap;
std::vector<MemoryState> productToMemoryStateMap;
std::vector<uint64_t> choiceToStateMap;
storm::storage::BitVector prob1Objectives; /// Objectives that are already satisfied in the initial state
std::vector<boost::optional<storm::storage::BitVector>> prob1InitialStates; /// For each objective the set of initial states that already satisfy the objective
};
}
}

666
src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.cpp

@ -8,6 +8,7 @@
#include "storm/environment/solver/MinMaxSolverEnvironment.h"
#include "storm/models/sparse/Mdp.h"
#include "storm/models/sparse/Dtmc.h"
#include "storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h"
#include "storm/storage/expressions/ExpressionManager.h"
#include "storm/storage/expressions/Expressions.h"
@ -21,6 +22,7 @@
#include "storm/logic/BoundedUntilFormula.h"
#include "storm/exceptions/NotSupportedException.h"
#include "storm/exceptions/UnexpectedException.h"
namespace storm {
namespace modelchecker {
@ -29,25 +31,38 @@ namespace storm {
template<typename ModelType>
QuantileHelper<ModelType>::QuantileHelper(ModelType const& model, storm::logic::QuantileFormula const& quantileFormula) : model(model), quantileFormula(quantileFormula) {
// Intentionally left empty
// Do all kinds of sanity check.
std::set<storm::expressions::Variable> quantileVariables;
for (auto const& quantileVariable : quantileFormula.getBoundVariables()) {
STORM_LOG_THROW(quantileVariables.count(quantileVariable) == 0, storm::exceptions::NotSupportedException, "Quantile formula considers the same bound variable twice.");
quantileVariables.insert(quantileVariable);
}
STORM_LOG_THROW(quantileFormula.getSubformula().isProbabilityOperatorFormula(), storm::exceptions::NotSupportedException, "Quantile formula needs probability operator inside. The formula " << quantileFormula << " is not supported.");
auto const& probOpFormula = quantileFormula.getSubformula().asProbabilityOperatorFormula();
STORM_LOG_THROW(probOpFormula.getBound().comparisonType == storm::logic::ComparisonType::Greater || probOpFormula.getBound().comparisonType == storm::logic::ComparisonType::LessEqual, storm::exceptions::NotSupportedException, "Probability operator inside quantile formula needs to have bound > or <=.");
STORM_LOG_THROW(probOpFormula.hasBound(), storm::exceptions::InvalidOperationException, "Probability operator inside quantile formula needs to have a bound.");
STORM_LOG_THROW(!model.isNondeterministicModel() || probOpFormula.hasOptimalityType(), storm::exceptions::InvalidOperationException, "Probability operator inside quantile formula needs to have an optimality type.");
STORM_LOG_WARN_COND(probOpFormula.getBound().comparisonType == storm::logic::ComparisonType::Greater || probOpFormula.getBound().comparisonType == storm::logic::ComparisonType::LessEqual, "Probability operator inside quantile formula needs to have bound > or <=. The specified comparison type might lead to non-termination."); // This has to do with letting bound variables approach infinity, e.g., Pr>0.7 [F "goal"] holds iff Pr>0.7 [F<=B "goal"] holds for some B.
STORM_LOG_THROW(probOpFormula.getSubformula().isBoundedUntilFormula(), storm::exceptions::NotSupportedException, "Quantile formula needs bounded until probability operator formula as subformula. The formula " << quantileFormula << " is not supported.");
auto const& boundedUntilFormula = quantileFormula.getSubformula().asProbabilityOperatorFormula().getSubformula().asBoundedUntilFormula();
// Only > and <= are supported for upper bounds. This is to make sure that Pr>0.7 [F "goal"] holds iff Pr>0.7 [F<=B "goal"] holds for some B.
// Only >= and < are supported for lower bounds. (EC construction..)
// TODO
// Bounds are either constants or variables that are declared in the quantile formula.
// Prop op has optimality type
// No Prmin with lower cost bounds: Ec construction fails. In single obj we would do 1-Prmax[F "nonRewOrNonGoalEC"] but this invalidates other lower/upper cost bounds.
/* Todo: Current assumptions:
* Subformula is always prob operator with bounded until
* Each bound variable occurs at most once (change this?)
* cost bounds are assumed to be always non-strict (the epochs returned by the reward unfolding assume strict lower and non-strict upper cost bounds, though...)
* 'reasonable' quantile (e.g. not quantile(max B, Pmax>0.5 [F <=B G]) (we could also filter out these things early on)
*/
auto const& boundedUntilFormula = probOpFormula.getSubformula().asBoundedUntilFormula();
std::set<storm::expressions::Variable> boundVariables;
for (uint64_t dim = 0; dim < boundedUntilFormula.getDimension(); ++dim) {
storm::expressions::Expression boundExpression;
if (boundedUntilFormula.hasUpperBound(dim)) {
STORM_LOG_THROW(!boundedUntilFormula.hasLowerBound(dim), storm::exceptions::NotSupportedException, "Interval bounds are not supported within quantile formulas.");
STORM_LOG_THROW(!boundedUntilFormula.isUpperBoundStrict(dim), storm::exceptions::NotSupportedException, "Only non-strict upper reward bounds are supported for quantiles.");
boundExpression = boundedUntilFormula.getUpperBound(dim);
} else if (boundedUntilFormula.hasLowerBound(dim)) {
STORM_LOG_THROW(!boundedUntilFormula.isLowerBoundStrict(dim), storm::exceptions::NotSupportedException, "Only non-strict lower reward bounds are supported for quantiles.");
boundExpression = boundedUntilFormula.getLowerBound(dim);
}
if (boundExpression.isInitialized() && boundExpression.containsVariables()) {
STORM_LOG_THROW(boundExpression.isVariable(), storm::exceptions::NotSupportedException, "Non-trivial bound expressions such as '" << boundExpression << "' are not supported. Either specify a constant or a quantile variable.");
storm::expressions::Variable const& boundVariable = boundExpression.getBaseExpression().asVariableExpression().getVariable();
STORM_LOG_THROW(boundVariables.count(boundVariable) == 0, storm::exceptions::NotSupportedException, "Variable " << boundExpression << " occurs at multiple reward bounds.");
boundVariables.insert(boundVariable);
STORM_LOG_THROW(quantileVariables.count(boundVariable) == 1, storm::exceptions::NotSupportedException, "The formula contains undefined constant '" << boundExpression << "'.");
}
}
}
enum class BoundTransformation {
@ -56,7 +71,7 @@ namespace storm {
GreaterEqualZero,
LessEqualZero
};
std::shared_ptr<storm::logic::ProbabilityOperatorFormula> transformBoundedUntilOperator(storm::logic::ProbabilityOperatorFormula const& boundedUntilOperator, std::vector<BoundTransformation> const& transformations) {
std::shared_ptr<storm::logic::ProbabilityOperatorFormula> transformBoundedUntilOperator(storm::logic::ProbabilityOperatorFormula const& boundedUntilOperator, std::vector<BoundTransformation> const& transformations, bool complementQuery = false) {
auto const& origBoundedUntil = boundedUntilOperator.getSubformula().asBoundedUntilFormula();
STORM_LOG_ASSERT(transformations.size() == origBoundedUntil.getDimension(), "Tried to replace the bound of a dimension that is higher than the number of dimensions of the formula.");
std::vector<std::shared_ptr<storm::logic::Formula const>> leftSubformulas, rightSubformulas;
@ -105,7 +120,11 @@ namespace storm {
} else {
newBoundedUntil = std::make_shared<storm::logic::BoundedUntilFormula>(origBoundedUntil.getLeftSubformula().asSharedPointer(), origBoundedUntil.getRightSubformula().asSharedPointer(), lowerBounds, upperBounds, timeBoundReferences);
}
return std::make_shared<storm::logic::ProbabilityOperatorFormula>(newBoundedUntil, boundedUntilOperator.getOperatorInformation());
storm::logic::OperatorInformation newOpInfo(boundedUntilOperator.getOperatorInformation().optimalityType, boundedUntilOperator.getBound());
if (complementQuery) {
newOpInfo.bound->comparisonType = storm::logic::invert(newOpInfo.bound->comparisonType);
}
return std::make_shared<storm::logic::ProbabilityOperatorFormula>(newBoundedUntil, newOpInfo);
}
/// Increases the precision of solver results
@ -116,9 +135,12 @@ namespace storm {
env.solver().minMax().setPrecision(env.solver().minMax().getPrecision() * factor);
}
/// Computes a lower/ upper bound on the actual result of a minmax or linear equation solver
/*!
* Computes a lower / upper bound on the actual result of a sound minmax or linear equation solver
*
*/
template<typename ValueType>
std::pair<ValueType, ValueType> getLowerUpperBound(storm::Environment const& env, ValueType const& value, bool minMax = true) {
std::pair<ValueType, ValueType> getLowerUpperBound(storm::Environment const& env, ValueType const& factor, ValueType const& value, bool minMax = true) {
ValueType prec;
bool relative;
if (minMax) {
@ -128,6 +150,7 @@ namespace storm {
prec = storm::utility::convertNumber<ValueType>(env.solver().getPrecisionOfLinearEquationSolver(env.solver().getLinearEquationSolverType()).first.get());
relative = env.solver().getPrecisionOfLinearEquationSolver(env.solver().getLinearEquationSolverType()).second.get();
}
prec *= factor;
if (relative) {
return std::make_pair<ValueType, ValueType>(value * (1/(prec + 1)), value * (1 + prec/(prec +1)));
} else {
@ -152,19 +175,7 @@ namespace storm {
}
return res;
}
template<typename ModelType>
storm::solver::OptimizationDirection const& QuantileHelper<ModelType>::getOptimizationDirForDimension(uint64_t const& dim) const {
storm::expressions::Variable const& dimVar = getVariableForDimension(dim);
for (auto const& boundVar : quantileFormula.getBoundVariables()) {
if (boundVar.second == dimVar) {
return boundVar.first;
}
}
STORM_LOG_THROW(false, storm::exceptions::InvalidOperationException, "The bound variable '" << dimVar.getName() << "' is not specified within the quantile formula '" << quantileFormula << "'.");
return quantileFormula.getOptimizationDirection();
}
template<typename ModelType>
storm::expressions::Variable const& QuantileHelper<ModelType>::getVariableForDimension(uint64_t const& dim) const {
auto const& boundedUntil = quantileFormula.getSubformula().asProbabilityOperatorFormula().getSubformula().asBoundedUntilFormula();
@ -172,432 +183,269 @@ namespace storm {
}
template<typename ModelType>
storm::storage::BitVector QuantileHelper<ModelType>::getDimensionsForVariable(storm::expressions::Variable const& var) const {
auto const& boundedUntil = quantileFormula.getSubformula().asProbabilityOperatorFormula().getSubformula().asBoundedUntilFormula();
storm::storage::BitVector result(boundedUntil.getDimension(), false);
for (uint64_t dim = 0; dim < boundedUntil.getDimension(); ++dim) {
if (boundedUntil.hasLowerBound(dim) && boundedUntil.getLowerBound(dim).isVariable() && boundedUntil.getLowerBound(dim).getBaseExpression().asVariableExpression().getVariable() == var) {
result.set(dim, true);
}
if (boundedUntil.hasUpperBound(dim) && boundedUntil.getUpperBound(dim).isVariable() && boundedUntil.getUpperBound(dim).getBaseExpression().asVariableExpression().getVariable() == var) {
result.set(dim, true);
}
}
return result;
}
template<typename ModelType>
std::vector<std::vector<typename ModelType::ValueType>> QuantileHelper<ModelType>::computeQuantile(Environment const& env) const {
std::vector<std::vector<typename ModelType::ValueType>> QuantileHelper<ModelType>::computeQuantile(Environment const& env) {
numCheckedEpochs = 0;
numPrecisionRefinements = 0;
swEpochAnalysis.reset();
swExploration.reset();
cachedSubQueryResults.clear();
std::vector<std::vector<ValueType>> result;
Environment envCpy = env; // It might be necessary to increase the precision during the computation
numCheckedEpochs = 0; numPrecisionRefinements = 0;
if (getOpenDimensions().getNumberOfSetBits() == 1) {
uint64_t dimension = *getOpenDimensions().begin();
bool zeroSatisfiesFormula = checkLimitValue(envCpy, storm::storage::BitVector(getDimension(), false));
bool infSatisfiesFormula = checkLimitValue(envCpy, getOpenDimensions());
if (zeroSatisfiesFormula != infSatisfiesFormula) {
while (true) {
auto quantileRes = computeQuantileForDimension(envCpy, dimension);
if (quantileRes) {
result = {{storm::utility::convertNumber<ValueType>(quantileRes->first) * quantileRes->second}};
break;
}
increasePrecision(envCpy);
++numPrecisionRefinements;
}
} else if (zeroSatisfiesFormula) {
// thus also infSatisfiesFormula is true, i.e., all bound values satisfy the formula
if (storm::solver::minimize(getOptimizationDirForDimension(dimension))) {
result = {{ storm::utility::zero<ValueType>() }};
} else {
result = {{ storm::utility::infinity<ValueType>()}};
// Call the internal recursive function
auto internalResult = computeQuantile(envCpy, getOpenDimensions(), false);
// Translate the result by applying the scaling factors and permutation.
std::vector<uint64_t> permutation;
for (auto const& v : quantileFormula.getBoundVariables()) {
uint64_t openDim = 0;
for (auto const& dim : getOpenDimensions()) {
if (getVariableForDimension(dim) == v) {
permutation.push_back(openDim);
break;
}
} else {
// i.e., no bound value satisfies the formula
result = {{}};
++openDim;
}
} else if (getOpenDimensions().getNumberOfSetBits() == 2) {
result = computeTwoDimensionalQuantile(envCpy);
} else {
STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "The quantile formula considers " << getOpenDimensions().getNumberOfSetBits() << " open dimensions. Only one- or two-dimensional quantiles are supported.");
}
assert(permutation.size() == getOpenDimensions().getNumberOfSetBits());
for (auto const& costLimits : internalResult.first.getGenerator()) {
std::vector<ValueType> resultPoint;
for (auto const& dim : permutation) {
CostLimit const& cl = costLimits[dim];
resultPoint.push_back(cl.isInfinity() ? storm::utility::infinity<ValueType>() : storm::utility::convertNumber<ValueType>(cl.get()) * internalResult.second[dim]);
}
result.push_back(resultPoint);
}
if (storm::settings::getModule<storm::settings::modules::CoreSettings>().isShowStatisticsSet()) {
std::cout << "Number of checked epochs: " << numCheckedEpochs << std::endl;
std::cout << "Number of required precision refinements: " << numPrecisionRefinements << std::endl;
std::cout << "Time for epoch exploration: " << swExploration << " seconds." << std::endl;
std::cout << "\tTime for epoch model analysis: " << swEpochAnalysis << " seconds." << std::endl;
}
return result;
}
template<typename ValueType>
void filterDominatedPoints(std::vector<std::vector<ValueType>>& points, std::vector<storm::solver::OptimizationDirection> const& dirs) {
std::vector<std::vector<ValueType>> result;
// Note: this is slow and not inplace but also most likely not performance critical
for (auto const& p1 : points) {
bool p1Dominated = false;
for (auto const& p2 : points) {
assert(p1.size() == p2.size());
bool p2DominatesP1 = false;
for (uint64_t i = 0; i < dirs.size(); ++i) {
if (storm::solver::minimize(dirs[i]) ? p2[i] <= p1[i] : p2[i] >= p1[i]) {
if (p2[i] != p1[i]) {
p2DominatesP1 = true;
}
} else {
p2DominatesP1 = false;
break;
}
}
if (p2DominatesP1) {
p1Dominated = true;
break;
}
}
if (!p1Dominated) {
result.push_back(p1);
}
}
points = std::move(result);
}
template<typename ModelType>
std::vector<std::vector<typename ModelType::ValueType>> QuantileHelper<ModelType>::computeTwoDimensionalQuantile(Environment& env) const {
std::vector<std::vector<ValueType>> result;
std::pair<CostLimitClosure, std::vector<typename QuantileHelper<ModelType>::ValueType>> QuantileHelper<ModelType>::computeQuantile(Environment& env, storm::storage::BitVector const& consideredDimensions, bool complementaryQuery) {
STORM_LOG_ASSERT(consideredDimensions.isSubsetOf(getOpenDimensions()), "Considered dimensions for a quantile query should be a subset of the set of dimensions without a fixed bound.");
auto const& probOpFormula = quantileFormula.getSubformula().asProbabilityOperatorFormula();
auto const& boundedUntilFormula = probOpFormula.getSubformula().asBoundedUntilFormula();
std::vector<storm::storage::BitVector> dimensionsAsBitVector;
std::vector<uint64_t> dimensions;
std::vector<storm::solver::OptimizationDirection> optimizationDirections;
std::vector<bool> lowerCostBounds;
for (auto const& dirVar : quantileFormula.getBoundVariables()) {
dimensionsAsBitVector.push_back(getDimensionsForVariable(dirVar.second));
STORM_LOG_THROW(dimensionsAsBitVector.back().getNumberOfSetBits() == 1, storm::exceptions::NotSupportedException, "There is not exactly one reward bound referring to quantile variable '" << dirVar.second.getName() << "'.");
dimensions.push_back(*dimensionsAsBitVector.back().begin());
lowerCostBounds.push_back(boundedUntilFormula.hasLowerBound(dimensions.back()));
optimizationDirections.push_back(dirVar.first);
storm::storage::BitVector cacheKey = consideredDimensions;
cacheKey.resize(cacheKey.size() + 1, complementaryQuery);
auto cacheIt = cachedSubQueryResults.find(cacheKey);
if (cacheIt != cachedSubQueryResults.end()) {
return cacheIt->second;
}
STORM_LOG_ASSERT(dimensions.size() == 2, "Expected to have exactly two open dimensions.");
if (optimizationDirections[0] == optimizationDirections[1]) {
if (lowerCostBounds[0] == lowerCostBounds[1]) {
// TODO: Assert that we have a reasonable probability bound. STORM_LOG_THROW(storm::solver::minimize(optimizationDirections[0])
bool infInfProbSatisfiesFormula = checkLimitValue(env, storm::storage::BitVector(getDimension(), false));
bool zeroZeroProbSatisfiesFormula = checkLimitValue(env, dimensionsAsBitVector[0] | dimensionsAsBitVector[1]);
if (infInfProbSatisfiesFormula != zeroZeroProbSatisfiesFormula) {
std::vector<std::pair<int64_t, typename ModelType::ValueType>> singleQuantileValues;
for (uint64_t i = 0; i < 2; ++i) {
// find out whether the bounds B_i = 0 and B_1-i = infinity satisfy the formula
bool zeroInfSatisfiesFormula = checkLimitValue(env, dimensionsAsBitVector[1-i]);
if (zeroInfSatisfiesFormula == storm::solver::minimize(optimizationDirections[0])) {
// There is bound value b such that the point B_i=0 and B_1-i = b is part of the result
singleQuantileValues.emplace_back(0, storm::utility::zero<ValueType>());
auto boundedUntilOp = transformBoundedUntilOperator(quantileFormula.getSubformula().asProbabilityOperatorFormula(), std::vector<BoundTransformation>(getDimension(), BoundTransformation::None), complementaryQuery);
std::set<storm::expressions::Variable> infinityVariables;
storm::storage::BitVector lowerBoundedDimensions(getDimension());
storm::storage::BitVector downwardClosedDimensions(getDimension());
bool hasLowerValueBound = storm::logic::isLowerBound(boundedUntilOp->getComparisonType());
for (auto const& d : getOpenDimensions()) {
if (consideredDimensions.get(d)) {
bool hasLowerCostBound = boundedUntilOp->getSubformula().asBoundedUntilFormula().hasLowerBound(d);
lowerBoundedDimensions.set(d, hasLowerCostBound);
downwardClosedDimensions.set(d, hasLowerCostBound == hasLowerValueBound);
} else {
infinityVariables.insert(getVariableForDimension(d));
}
}
downwardClosedDimensions = downwardClosedDimensions % consideredDimensions;
CostLimitClosure satCostLimits(downwardClosedDimensions), unsatCostLimits(~downwardClosedDimensions);
// Initialize the (un)sat cost limits to guarantee termination
bool onlyUpperCostBounds = lowerBoundedDimensions.empty();
bool onlyLowerCostBounds = lowerBoundedDimensions == consideredDimensions;
if (onlyUpperCostBounds || onlyLowerCostBounds) {
for (auto const& k : consideredDimensions) {
storm::storage::BitVector subQueryDimensions = consideredDimensions;
subQueryDimensions.set(k, false);
bool subQueryComplement = complementaryQuery != ((onlyUpperCostBounds && hasLowerValueBound) || (onlyLowerCostBounds && !hasLowerValueBound));
auto subQueryResult = computeQuantile(env, subQueryDimensions, subQueryComplement);
for (auto const& subQueryCostLimit : subQueryResult.first.getGenerator()) {
CostLimits initPoint;
uint64_t i = 0;
for (auto const& dim : consideredDimensions) {
if (dim == k) {
initPoint.push_back(CostLimit::infinity());
} else {
// Compute quantile where 1-i is set to infinity
singleQuantileValues.push_back(computeQuantileForDimension(env, dimensions[i]).get());
if (!storm::solver::minimize(optimizationDirections[i])) {
// When maximizing bounds, the computed single dimensional quantile value is sat for all values of the other bound.
// Increase the computed value so that there is at least one unsat assignment of bounds with B[i] = singleQuantileValues[i]
++singleQuantileValues.back().first;
}
initPoint.push_back(subQueryCostLimit[i]);
++i;
}
// Decrease the value for lower cost bounds to convert >= to >
if (lowerCostBounds[i]) {
--singleQuantileValues[i].first;
}
}
std::vector<int64_t> currentEpochValues = {(int64_t) singleQuantileValues[0].first, (int64_t) singleQuantileValues[1].first}; // signed int is needed to allow lower bounds >-1 (aka >=0).
while (!exploreTwoDimensionalQuantile(env, singleQuantileValues, currentEpochValues, result)) {
increasePrecision(env);
++numPrecisionRefinements;
}
} else if (infInfProbSatisfiesFormula) {
// then also zeroZeroProb satisfies the formula, i.e., all bound values satisfy the formula
if (storm::solver::minimize(optimizationDirections[0])) {
result = {{storm::utility::zero<ValueType>(), storm::utility::zero<ValueType>()}};
if (subQueryComplement == complementaryQuery) {
satCostLimits.insert(initPoint);
} else {
result = {{ storm::utility::infinity<ValueType>(), storm::utility::infinity<ValueType>()}};
unsatCostLimits.insert(initPoint);
}
} else {
// i.e., no bound value satisfies the formula
result = {{}};
}
} else {
// TODO: this is an "unreasonable" case
}
} else {
// TODO: find reasonable min/max cases
STORM_LOG_WARN("Quantile formula considers mixtures of upper and lower reward-bounds. Termination is not guaranteed.");
}
return result;
}
template<typename ModelType>
bool QuantileHelper<ModelType>::exploreTwoDimensionalQuantile(Environment const& env, std::vector<std::pair<int64_t, typename ModelType::ValueType>> const& startEpochValues, std::vector<int64_t>& currentEpochValues, std::vector<std::vector<ValueType>>& resultPoints) const {
// Initialize some data for easy access
auto const& probOpFormula = quantileFormula.getSubformula().asProbabilityOperatorFormula();
auto const& boundedUntilFormula = probOpFormula.getSubformula().asBoundedUntilFormula();
std::vector<storm::storage::BitVector> dimensionsAsBitVector;
std::vector<uint64_t> dimensions;
std::vector<storm::solver::OptimizationDirection> optimizationDirections;
std::vector<bool> lowerCostBounds;
for (auto const& dirVar : quantileFormula.getBoundVariables()) {
dimensionsAsBitVector.push_back(getDimensionsForVariable(dirVar.second));
STORM_LOG_THROW(dimensionsAsBitVector.back().getNumberOfSetBits() == 1, storm::exceptions::NotSupportedException, "There is not exactly one reward bound referring to quantile variable '" << dirVar.second.getName() << "'.");
dimensions.push_back(*dimensionsAsBitVector.back().begin());
lowerCostBounds.push_back(boundedUntilFormula.hasLowerBound(dimensions.back()));
optimizationDirections.push_back(dirVar.first);
}
STORM_LOG_ASSERT(dimensions.size() == 2, "Expected to have exactly two open dimensions.");
MultiDimensionalRewardUnfolding<ValueType, true> rewardUnfolding(model, transformBoundedUntilOperator(quantileFormula.getSubformula().asProbabilityOperatorFormula(), std::vector<BoundTransformation>(getDimension(), BoundTransformation::None)));
// initialize data that will be needed for each epoch
auto lowerBound = rewardUnfolding.getLowerObjectiveBound();
auto upperBound = rewardUnfolding.getUpperObjectiveBound();
std::vector<ValueType> x, b;
std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> minMaxSolver;
if (currentEpochValues[0] < 0 && currentEpochValues[1] < 0) {
// This case can only happen in these cases:
assert(lowerCostBounds[0]);
assert(currentEpochValues[0] == -1);
assert(currentEpochValues[1] == -1);
// This case has been checked already, so we can skip it.
// Skipping this is actually necessary, since the rewardUnfolding will handle formulas like [F{"a"}>A,{"b"}>B "init"] incorrectly if A=B=-1.
++currentEpochValues[1];
}
std::set<typename EpochManager::Epoch> exploredEpochs;
// Loop until the goal precision is reached.
STORM_LOG_DEBUG("Computing quantile for dimensions: " << consideredDimensions);
while (true) {
auto currEpoch = rewardUnfolding.getStartEpoch(true);
for (uint64_t i = 0; i < 2; ++i) {
if (currentEpochValues[i] >= 0) {
rewardUnfolding.getEpochManager().setDimensionOfEpoch(currEpoch, dimensions[i], (uint64_t) currentEpochValues[i]);
} else {
rewardUnfolding.getEpochManager().setBottomDimension(currEpoch, dimensions[i]);
}
}
auto epochOrder = rewardUnfolding.getEpochComputationOrder(currEpoch);
for (auto const& epoch : epochOrder) {
if (exploredEpochs.count(epoch) == 0) {
auto& epochModel = rewardUnfolding.setCurrentEpoch(epoch);
rewardUnfolding.setSolutionForCurrentEpoch(epochModel.analyzeSingleObjective(env, quantileFormula.getSubformula().asProbabilityOperatorFormula().getOptimalityType(), x, b, minMaxSolver, lowerBound, upperBound));
++numCheckedEpochs;
exploredEpochs.insert(epoch);
// initialize reward unfolding and data that will be needed for each epoch
MultiDimensionalRewardUnfolding<ValueType, true> rewardUnfolding(model, boundedUntilOp, infinityVariables);
if (computeQuantile(env, consideredDimensions, *boundedUntilOp, lowerBoundedDimensions, satCostLimits, unsatCostLimits, rewardUnfolding)) {
std::vector<ValueType> scalingFactors;
for (auto const& dim : consideredDimensions) {
scalingFactors.push_back(rewardUnfolding.getDimension(dim).scalingFactor);
}
std::pair<CostLimitClosure, std::vector<ValueType>> result(satCostLimits, scalingFactors);
cachedSubQueryResults.emplace(cacheKey, result);
return result;
}
ValueType currValue = rewardUnfolding.getInitialStateResult(currEpoch);
bool propertySatisfied;
if (env.solver().isForceSoundness()) {
auto lowerUpperValue = getLowerUpperBound(env, currValue);
propertySatisfied = probOpFormula.getBound().isSatisfied(lowerUpperValue.first);
if (propertySatisfied != probOpFormula.getBound().isSatisfied(lowerUpperValue.second)) {
// unclear result due to insufficient precision.
return false;
}
} else {
propertySatisfied = probOpFormula.getBound().isSatisfied(currValue);
}
// If the reward bounds should be as small as possible, we should stop as soon as the property is satisfied.
// If the reward bounds should be as large as possible, we should stop as soon as the property is violated and then go a step backwards
if (storm::solver::minimize(optimizationDirections[0]) == propertySatisfied) {
// We found another point for the result! Translate it to the original domain
auto point = currentEpochValues;
std::vector<ValueType> convertedPoint;
for (uint64_t i = 0; i < 2; ++i) {
if (lowerCostBounds[i]) {
// Translate > to >=
++point[i];
}
if (i == 1 && storm::solver::maximize(optimizationDirections[i])) {
// When maximizing, we actually searched for each x-value the smallest y-value that leads to a property violation. Hence, decreasing y by one means property satisfaction
--point[i];
}
if (point[i] < 0) {
// Skip this point
convertedPoint.clear();
continue;
}
convertedPoint.push_back(storm::utility::convertNumber<ValueType>(point[i]) * rewardUnfolding.getDimension(dimensions[i]).scalingFactor);
}
if (!convertedPoint.empty()) {
resultPoints.push_back(std::move(convertedPoint));
}
if (currentEpochValues[1] == startEpochValues[1].first) {
break;
} else {
++currentEpochValues[0];
currentEpochValues[1] = startEpochValues[1].first;
}
} else {
++currentEpochValues[1];
}
STORM_LOG_WARN("Restarting quantile computation after " << swExploration << " seconds due to insufficient precision.");
++numPrecisionRefinements;
increasePrecision(env);
}
// When maximizing, there are border cases where one dimension can be arbitrarily large
for (uint64_t i = 0; i < 2; ++i) {
if (storm::solver::maximize(optimizationDirections[i]) && startEpochValues[i].first > 0) {
std::vector<ValueType> newResultPoint(2);
newResultPoint[i] = storm::utility::convertNumber<ValueType>(startEpochValues[i].first - 1) * startEpochValues[i].second;
newResultPoint[1-i] = storm::utility::infinity<ValueType>();
resultPoints.push_back(newResultPoint);
}
}
filterDominatedPoints(resultPoints, optimizationDirections);
return true;
}
template<typename ModelType>
boost::optional<std::pair<uint64_t, typename ModelType::ValueType>> QuantileHelper<ModelType>::computeQuantileForDimension(Environment const& env, uint64_t dimension) const {
// We assume that there is one bound value that violates the quantile and one bound value that satisfies it.
// Let all other open bounds approach infinity
std::set<storm::expressions::Variable> infinityVariables;
for (auto const& d : getOpenDimensions()) {
if (d != dimension) {
infinityVariables.insert(getVariableForDimension(d));
}
bool getNextCandidateCostLimit(CostLimit const& candidateCostLimitSum, CostLimits& current) {
if (current.size() == 0) {
return false;
}
auto transformedFormula = transformBoundedUntilOperator(quantileFormula.getSubformula().asProbabilityOperatorFormula(), std::vector<BoundTransformation>(getDimension(), BoundTransformation::None));
MultiDimensionalRewardUnfolding<ValueType, true> rewardUnfolding(model, transformedFormula, infinityVariables);
bool upperCostBound = quantileFormula.getSubformula().asProbabilityOperatorFormula().getSubformula().asBoundedUntilFormula().hasUpperBound(dimension);
bool minimizingRewardBound = storm::solver::minimize(getOptimizationDirForDimension(dimension));
// initialize data that will be needed for each epoch
auto lowerBound = rewardUnfolding.getLowerObjectiveBound();
auto upperBound = rewardUnfolding.getUpperObjectiveBound();
std::vector<ValueType> x, b;
std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> minMaxSolver;
uint64_t epochValue = 0;
std::set<typename EpochManager::Epoch> exploredEpochs;
while (true) {
auto currEpoch = rewardUnfolding.getStartEpoch(true);
rewardUnfolding.getEpochManager().setDimensionOfEpoch(currEpoch, dimension, (uint64_t) epochValue);
auto epochOrder = rewardUnfolding.getEpochComputationOrder(currEpoch);
for (auto const& epoch : epochOrder) {
if (exploredEpochs.count(epoch) == 0) {
auto& epochModel = rewardUnfolding.setCurrentEpoch(epoch);
rewardUnfolding.setSolutionForCurrentEpoch(epochModel.analyzeSingleObjective(env, quantileFormula.getSubformula().asProbabilityOperatorFormula().getOptimalityType(), x, b, minMaxSolver, lowerBound, upperBound));
++numCheckedEpochs;
exploredEpochs.insert(epoch);
}
}
ValueType currValue = rewardUnfolding.getInitialStateResult(currEpoch);
bool propertySatisfied;
if (env.solver().isForceSoundness()) {
auto lowerUpperValue = getLowerUpperBound(env, currValue);
propertySatisfied = transformedFormula->getBound().isSatisfied(lowerUpperValue.first);
if (propertySatisfied != transformedFormula->getBound().isSatisfied(lowerUpperValue.second)) {
// unclear result due to insufficient precision.
return boost::none;
}
} else {
propertySatisfied = transformedFormula->getBound().isSatisfied(currValue);
}
// If the reward bound should be as small as possible, we should stop as soon as the property is satisfied.
// If the reward bound should be as large as possible, we should stop as soon as sthe property is violated and then go a step backwards
if (minimizingRewardBound && propertySatisfied) {
// We found a satisfying value!
if (!upperCostBound) {
// The rewardunfolding assumes strict lower cost bounds while we assume non-strict ones. Hence, >B becomes >=(B+1).
++epochValue;
}
break;
} else if (!minimizingRewardBound && !propertySatisfied) {
// We found a non-satisfying value. Go one step back to get the largest satisfying value.
// ... however, lower cost bounds need to be converted from strict bounds >B to non strict bounds >=(B+1).
// Hence, no operation is necessary in this case.
if (upperCostBound) {
STORM_LOG_ASSERT(epochValue > 0, "The property does not seem to be satisfiable. This case should have been excluded earlier.");
--epochValue;
}
break;
uint64_t iSum = current.front().get();
if (iSum == candidateCostLimitSum.get()) {
return false;
}
for (uint64_t i = 1; i < current.size(); ++i) {
iSum += current[i].get();
if (iSum == candidateCostLimitSum.get()) {
++current[i-1].get();
uint64_t newVal = current[i].get() - 1;
current[i].get() = 0;
current.back().get() = newVal;
return true;
}
++epochValue;
}
return std::make_pair(epochValue, rewardUnfolding.getDimension(dimension).scalingFactor);
STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "The entries of the current cost limit candidate do not sum up to the current candidate sum.");
return false;
}
template<typename ModelType>
bool QuantileHelper<ModelType>::checkLimitValue(Environment& env, storm::storage::BitVector const& infDimensions) const {
auto const& probabilityBound = quantileFormula.getSubformula().asProbabilityOperatorFormula().getBound();
// Increase the precision until we get a conclusive result
while (true) {
ValueType numericResult = computeLimitValue(env, infDimensions);
if (env.solver().isForceSoundness()) {
auto lowerUpper = getLowerUpperBound(env, numericResult);
bool lowerSat = probabilityBound.isSatisfied(lowerUpper.first);
bool upperSat = probabilityBound.isSatisfied(lowerUpper.second);
if (lowerSat == upperSat) {
return lowerSat;
bool translateEpochToCostLimits(EpochManager::Epoch const& epoch, EpochManager::Epoch const& startEpoch,storm::storage::BitVector const& consideredDimensions, storm::storage::BitVector const& lowerBoundedDimensions, EpochManager const& epochManager, CostLimits& epochAsCostLimits) {
for (uint64_t dim = 0; dim < consideredDimensions.size(); ++dim) {
if (consideredDimensions.get(dim)) {
if (lowerBoundedDimensions.get(dim)) {
if (epochManager.isBottomDimension(epoch, dim)) {
epochAsCostLimits.push_back(CostLimit(0));
} else {
epochAsCostLimits.push_back(CostLimit(epochManager.getDimensionOfEpoch(epoch, dim) + 1));
}
} else {
increasePrecision(env);
++numPrecisionRefinements;
if (epochManager.isBottomDimension(epoch, dim)) {
return false;
} else {
epochAsCostLimits.push_back(CostLimit(epochManager.getDimensionOfEpoch(epoch, dim)));
}
}
} else {
return probabilityBound.isSatisfied(numericResult);
if (epochManager.isBottomDimension(epoch, dim)) {
if (!epochManager.isBottomDimension(startEpoch, dim)) {
return false;
}
} else if (epochManager.getDimensionOfEpoch(epoch, dim) != epochManager.getDimensionOfEpoch(startEpoch, dim)) {
return false;
}
}
}
return true;
}
template<typename ModelType>
typename ModelType::ValueType QuantileHelper<ModelType>::computeLimitValue(Environment const& env, storm::storage::BitVector const& infDimensions) const {
// To compute the limit for an upper bounded dimension, we can simply drop the bound
// To compute the limit for a lower bounded dimension, we only consider epoch steps that lie in an end component and check for the bound >0 instead.
// Notice, however, that this approach becomes problematic if, at the same time, we consider an upper bounded dimension with bound value zero.
std::vector<BoundTransformation> bts(getDimension(), BoundTransformation::None);
std::set<storm::expressions::Variable> infinityVariables;
for (auto const& d : getOpenDimensions()) {
bool upperCostBound = quantileFormula.getSubformula().asProbabilityOperatorFormula().getSubformula().asBoundedUntilFormula().hasUpperBound(d);
if (infDimensions.get(d)) {
infinityVariables.insert(getVariableForDimension(d));
} else {
if (upperCostBound) {
bts[d] = BoundTransformation::LessEqualZero;
} else {
bts[d] = BoundTransformation::GreaterEqualZero;
}
}
}
auto transformedFormula = transformBoundedUntilOperator(quantileFormula.getSubformula().asProbabilityOperatorFormula(), bts);
MultiDimensionalRewardUnfolding<ValueType, true> rewardUnfolding(model, transformedFormula, infinityVariables);
if (!rewardUnfolding.getProb1Objectives().empty()) {
assert(rewardUnfolding.getProb1Objectives().size() == 1);
// The probability is one.
return storm::utility::one<ValueType>();
}
// Get lower and upper bounds for the solution.
bool QuantileHelper<ModelType>::computeQuantile(Environment& env, storm::storage::BitVector const& consideredDimensions, storm::logic::ProbabilityOperatorFormula const& boundedUntilOperator, storm::storage::BitVector const& lowerBoundedDimensions, CostLimitClosure& satCostLimits, CostLimitClosure& unsatCostLimits, MultiDimensionalRewardUnfolding<ValueType, true>& rewardUnfolding) {
auto lowerBound = rewardUnfolding.getLowerObjectiveBound();
auto upperBound = rewardUnfolding.getUpperObjectiveBound();
// Initialize epoch models
auto initEpoch = rewardUnfolding.getStartEpoch();
auto epochOrder = rewardUnfolding.getEpochComputationOrder(initEpoch);
// initialize data that will be needed for each epoch
std::vector<ValueType> x, b;
std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> minMaxSolver;
for (auto const& epoch : epochOrder) {
auto& epochModel = rewardUnfolding.setCurrentEpoch(epoch);
rewardUnfolding.setSolutionForCurrentEpoch(epochModel.analyzeSingleObjective(env, quantileFormula.getSubformula().asProbabilityOperatorFormula().getOptimalityType(), x, b, minMaxSolver, lowerBound, upperBound));
++numCheckedEpochs;
std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> minMaxSolver; // Needed for MDP
std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> linEqSolver; // Needed for DTMC
if (!model.isNondeterministicModel()) {
rewardUnfolding.setEquationSystemFormatForEpochModel(storm::solver::GeneralLinearEquationSolverFactory<ValueType>().getEquationProblemFormat(env));
}
ValueType numericResult = rewardUnfolding.getInitialStateResult(initEpoch);
STORM_LOG_TRACE("Limit probability for infinity dimensions " << infDimensions << " is " << numericResult << ".");
return numericResult;
swExploration.start();
bool progress = true;
for (CostLimit candidateCostLimitSum(0); progress; ++candidateCostLimitSum.get()) {
CostLimits currentCandidate(satCostLimits.dimension(), CostLimit(0));
if (!currentCandidate.empty()) {
currentCandidate.back() = candidateCostLimitSum;
}
// We can still have progress if one of the closures is empty and the other is not full.
// This ensures that we do not terminate too early in case that the (un)satCostLimits are initially non-empty.
progress = (satCostLimits.empty() && !unsatCostLimits.full()) || (unsatCostLimits.empty() && !satCostLimits.full());
do {
if (!satCostLimits.contains(currentCandidate) && !unsatCostLimits.contains(currentCandidate)) {
progress = true;
// Transform candidate cost limits to an appropriate start epoch
auto startEpoch = rewardUnfolding.getStartEpoch(true);
auto costLimitIt = currentCandidate.begin();
for (auto const& dim : consideredDimensions) {
if (lowerBoundedDimensions.get(dim)) {
if (costLimitIt->get() > 0) {
rewardUnfolding.getEpochManager().setDimensionOfEpoch(startEpoch, dim, costLimitIt->get() - 1);
} else {
rewardUnfolding.getEpochManager().setBottomDimension(startEpoch, dim);
}
} else {
rewardUnfolding.getEpochManager().setDimensionOfEpoch(startEpoch, dim, costLimitIt->get());
}
++costLimitIt;
}
STORM_LOG_DEBUG("Checking start epoch " << rewardUnfolding.getEpochManager().toString(startEpoch) << ".");
auto epochSequence = rewardUnfolding.getEpochComputationOrder(startEpoch, true);
for (auto const& epoch : epochSequence) {
++numCheckedEpochs;
swEpochAnalysis.start();
auto& epochModel = rewardUnfolding.setCurrentEpoch(epoch);
if (model.isNondeterministicModel()) {
rewardUnfolding.setSolutionForCurrentEpoch(epochModel.analyzeSingleObjective(env, boundedUntilOperator.getOptimalityType(), x, b, minMaxSolver, lowerBound, upperBound));
} else {
rewardUnfolding.setSolutionForCurrentEpoch(epochModel.analyzeSingleObjective(env, x, b, linEqSolver, lowerBound, upperBound));
}
swEpochAnalysis.stop();
CostLimits epochAsCostLimits;
if (translateEpochToCostLimits(epoch, startEpoch, consideredDimensions, lowerBoundedDimensions, rewardUnfolding.getEpochManager(), epochAsCostLimits)) {
ValueType currValue = rewardUnfolding.getInitialStateResult(epoch);
bool propertySatisfied;
if (env.solver().isForceSoundness()) {
ValueType sumOfEpochDimensions = storm::utility::convertNumber<ValueType>(rewardUnfolding.getEpochManager().getSumOfDimensions(epoch) + 1);
auto lowerUpperValue = getLowerUpperBound(env, sumOfEpochDimensions, currValue);
propertySatisfied = boundedUntilOperator.getBound().isSatisfied(lowerUpperValue.first);
if (propertySatisfied != boundedUntilOperator.getBound().isSatisfied(lowerUpperValue.second)) {
// unclear result due to insufficient precision.
swExploration.stop();
return false;
}
} else {
propertySatisfied = boundedUntilOperator.getBound().isSatisfied(currValue);
}
if (propertySatisfied) {
satCostLimits.insert(epochAsCostLimits);
} else {
unsatCostLimits.insert(epochAsCostLimits);
}
}
}
}
} while (getNextCandidateCostLimit(candidateCostLimitSum, currentCandidate));
if (!progress) {
progress = !CostLimitClosure::unionFull(satCostLimits, unsatCostLimits);
}
}
swExploration.stop();
return true;
}
template class QuantileHelper<storm::models::sparse::Mdp<double>>;
template class QuantileHelper<storm::models::sparse::Mdp<storm::RationalNumber>>;
template class QuantileHelper<storm::models::sparse::Dtmc<double>>;
template class QuantileHelper<storm::models::sparse::Dtmc<storm::RationalNumber>>;
}
}

39
src/storm/modelchecker/prctl/helper/rewardbounded/QuantileHelper.h

@ -1,15 +1,18 @@
#pragma once
#include <map>
#include <boost/optional.hpp>
#include "storm/logic/QuantileFormula.h"
#include "boost/optional.hpp"
#include "storm/logic/ProbabilityOperatorFormula.h"
#include "storm/modelchecker/prctl/helper/rewardbounded/CostLimitClosure.h"
#include "storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h"
#include "storm/storage/BitVector.h"
#include "storm/utility/Stopwatch.h"
namespace storm {
class Environment;
namespace storage {
class BitVector;
}
namespace modelchecker {
namespace helper {
namespace rewardbounded {
@ -20,28 +23,14 @@ namespace storm {
public:
QuantileHelper(ModelType const& model, storm::logic::QuantileFormula const& quantileFormula);
std::vector<std::vector<ValueType>> computeQuantile(Environment const& env) const;
std::vector<std::vector<ValueType>> computeQuantile(Environment const& env);
private:
std::vector<std::vector<ValueType>> computeTwoDimensionalQuantile(Environment& env) const;
bool exploreTwoDimensionalQuantile(Environment const& env, std::vector<std::pair<int64_t, typename ModelType::ValueType>> const& startEpochValues, std::vector<int64_t>& currentEpochValues, std::vector<std::vector<ValueType>>& resultPoints) const;
/*!
* Computes the limit probability, where the given dimensions approach infinity and the remaining dimensions are set to zero.
*/
ValueType computeLimitValue(Environment const& env, storm::storage::BitVector const& infDimensions) const;
std::pair<CostLimitClosure, std::vector<ValueType>> computeQuantile(Environment& env, storm::storage::BitVector const& consideredDimensions, bool complementaryQuery);
bool computeQuantile(Environment& env, storm::storage::BitVector const& consideredDimensions, storm::logic::ProbabilityOperatorFormula const& boundedUntilOperator, storm::storage::BitVector const& lowerBoundedDimensions, CostLimitClosure& satCostLimits, CostLimitClosure& unsatCostLimits, MultiDimensionalRewardUnfolding<ValueType, true>& rewardUnfolding);
/*!
* Computes the limit probability, where the given dimensions approach infinity and the remaining dimensions are set to zero.
* The computed value is compared to the probability threshold.
* In sound mode, precision is iteratively increased in case of 'inconsistent' results.
*/
bool checkLimitValue(Environment& env, storm::storage::BitVector const& infDimensions) const;
/// Computes the quantile with respect to the given dimension.
/// boost::none is returned in case of insufficient precision.
boost::optional<std::pair<uint64_t, typename ModelType::ValueType>> computeQuantileForDimension(Environment const& env, uint64_t dim) const;
/*!
* Gets the number of dimensions of the underlying boudned until formula
*/
@ -54,15 +43,17 @@ namespace storm {
storm::storage::BitVector getOpenDimensions() const;
storm::storage::BitVector getDimensionsForVariable(storm::expressions::Variable const& var) const;
storm::solver::OptimizationDirection const& getOptimizationDirForDimension(uint64_t const& dim) const;
storm::expressions::Variable const& getVariableForDimension(uint64_t const& dim) const;
ModelType const& model;
storm::logic::QuantileFormula const& quantileFormula;
std::map<storm::storage::BitVector, std::pair<CostLimitClosure, std::vector<ValueType>>> cachedSubQueryResults;
/// Statistics
mutable uint64_t numCheckedEpochs;
mutable uint64_t numPrecisionRefinements;
mutable storm::utility::Stopwatch swEpochAnalysis;
mutable storm::utility::Stopwatch swExploration;
};
}
}

302
src/test/storm/modelchecker/QuantileQueryTest.cpp

@ -0,0 +1,302 @@
#include "gtest/gtest.h"
#include "storm-config.h"
#include "test/storm_gtest.h"
#include "storm/api/builder.h"
#include "storm-parsers/api/model_descriptions.h"
#include "storm/api/properties.h"
#include "storm-parsers/api/properties.h"
#include "storm/parser/CSVParser.h"
#include "storm/models/sparse/Mdp.h"
#include "storm/models/sparse/StandardRewardModel.h"
#include "storm/modelchecker/prctl/SparseMdpPrctlModelChecker.h"
#include "storm/modelchecker/prctl/SparseDtmcPrctlModelChecker.h"
#include "storm/modelchecker/results/ExplicitParetoCurveCheckResult.h"
#include "storm/modelchecker/results/ExplicitQualitativeCheckResult.h"
#include "storm/modelchecker/results/ExplicitQuantitativeCheckResult.h"
#include "storm/environment/solver/MinMaxSolverEnvironment.h"
#include "storm/logic/Formulas.h"
#include "storm/storage/jani/Property.h"
namespace {
enum class MdpEngine {PrismSparse, JaniSparse, JitSparse, Hybrid, PrismDd, JaniDd};
class UnsoundEnvironment {
public:
typedef double ValueType;
static storm::Environment createEnvironment() {
storm::Environment env;
env.solver().minMax().setPrecision(storm::utility::convertNumber<storm::RationalNumber>(1e-10));
return env;
}
};
class SoundEnvironment {
public:
typedef double ValueType;
static storm::Environment createEnvironment() {
storm::Environment env;
env.solver().setForceSoundness(true);
return env;
}
};
class ExactEnvironment {
public:
typedef storm::RationalNumber ValueType;
static storm::Environment createEnvironment() {
storm::Environment env;
return env;
}
};
template<typename TestType>
class QuantileQueryTest : public ::testing::Test {
public:
typedef typename TestType::ValueType ValueType;
QuantileQueryTest() : _environment(TestType::createEnvironment()) {}
storm::Environment const& env() const { return _environment; }
ValueType parseNumber(std::string const& input) const {
if (input.find("inf") != std::string::npos) {
return storm::utility::infinity<ValueType>();
}
return storm::utility::convertNumber<ValueType>(input);
}
template <typename MT>
std::pair<std::shared_ptr<MT>, std::vector<std::shared_ptr<storm::logic::Formula const>>> buildModelFormulas(std::string const& pathToPrismFile, std::string const& formulasAsString, std::string const& constantDefinitionString = "") const {
std::pair<std::shared_ptr<MT>, std::vector<std::shared_ptr<storm::logic::Formula const>>> result;
storm::prism::Program program = storm::api::parseProgram(pathToPrismFile);
program = storm::utility::prism::preprocess(program, constantDefinitionString);
result.second = storm::api::extractFormulasFromProperties(storm::api::parsePropertiesForPrismProgram(formulasAsString, program));
result.first = storm::api::buildSparseModel<ValueType>(program, result.second)->template as<MT>();
return result;
}
std::vector<storm::modelchecker::CheckTask<storm::logic::Formula, ValueType>> getTasks(std::vector<std::shared_ptr<storm::logic::Formula const>> const& formulas) const {
std::vector<storm::modelchecker::CheckTask<storm::logic::Formula, ValueType>> result;
for (auto const& f : formulas) {
result.emplace_back(*f, true); // only initial states are relevant
}
return result;
}
template <typename MT>
typename std::enable_if<std::is_same<MT, storm::models::sparse::Mdp<ValueType>>::value, std::shared_ptr<storm::modelchecker::AbstractModelChecker<MT>>>::type
createModelChecker(std::shared_ptr<MT> const& model) const {
return std::make_shared<storm::modelchecker::SparseMdpPrctlModelChecker<MT>>(*model);
}
template <typename MT>
typename std::enable_if<std::is_same<MT, storm::models::sparse::Dtmc<ValueType>>::value, std::shared_ptr<storm::modelchecker::AbstractModelChecker<MT>>>::type
createModelChecker(std::shared_ptr<MT> const& model) const {
return std::make_shared<storm::modelchecker::SparseDtmcPrctlModelChecker<MT>>(*model);
}
std::pair<bool, std::string> compareResult(std::shared_ptr<storm::models::sparse::Model<ValueType>> const& model, std::unique_ptr<storm::modelchecker::CheckResult>& result, std::vector<std::string> const& expected) {
bool equal = true;
std::string errorMessage = "";
auto filter = getInitialStateFilter(model);
result->filter(*filter);
std::vector<std::vector<ValueType>> resultPoints;
if (result->isExplicitParetoCurveCheckResult()) {
resultPoints = result->asExplicitParetoCurveCheckResult<ValueType>().getPoints();
} else {
if (!result->isExplicitQuantitativeCheckResult()) {
return { false, "The CheckResult has unexpected type."};
}
resultPoints = {{ result->asExplicitQuantitativeCheckResult<ValueType>().getMax() }};
}
std::vector<std::vector<ValueType>> expectedPoints;
for (auto const& pointAsString : expected) {
std::vector<ValueType> point;
for (auto const& entry : storm::parser::parseCommaSeperatedValues(pointAsString)) {
point.push_back(parseNumber(entry));
}
expectedPoints.push_back(std::move(point));
}
for (auto const& resPoint : resultPoints) {
bool contained = false;
for (auto const& expPoint : expectedPoints) {
if (resPoint == expPoint) {
contained = true;
break;
}
}
if (!contained) {
equal = false;
errorMessage += "The result " + storm::utility::vector::toString(resPoint) + " is not contained in the expected solution. ";
}
}
for (auto const& expPoint : expectedPoints) {
bool contained = false;
for (auto const& resPoint : resultPoints) {
if (resPoint == expPoint) {
contained = true;
break;
}
}
if (!contained) {
equal = false;
errorMessage += "The expected " + storm::utility::vector::toString(expPoint) + " is not contained in the obtained solution. ";
}
}
return {equal, errorMessage};
}
private:
storm::Environment _environment;
std::unique_ptr<storm::modelchecker::QualitativeCheckResult> getInitialStateFilter(std::shared_ptr<storm::models::sparse::Model<ValueType>> const& model) const {
return std::make_unique<storm::modelchecker::ExplicitQualitativeCheckResult>(model->getInitialStates());
}
};
typedef ::testing::Types<
UnsoundEnvironment,
SoundEnvironment,
ExactEnvironment
> TestingTypes;
TYPED_TEST_CASE(QuantileQueryTest, TestingTypes);
TYPED_TEST(QuantileQueryTest, simple_Dtmc) {
typedef storm::models::sparse::Dtmc<typename TestFixture::ValueType> ModelType;
std::string formulasString = "quantile(max A, max B, P>0.95 [F{\"first\"}<=A,{\"second\"}<=B s=3]);\n";
auto modelFormulas = this->template buildModelFormulas<ModelType>(STORM_TEST_RESOURCES_DIR "/dtmc/quantiles_simple_dtmc.pm", formulasString);
auto model = std::move(modelFormulas.first);
auto tasks = this->getTasks(modelFormulas.second);
EXPECT_EQ(4ul, model->getNumberOfStates());
EXPECT_EQ(6ul, model->getNumberOfTransitions());
auto checker = this->template createModelChecker<ModelType>(model);
std::unique_ptr<storm::modelchecker::CheckResult> result;
std::vector<std::string> expectedResult;
std::pair<bool, std::string> compare;
uint64_t taskId = 0;
expectedResult.clear();
expectedResult.push_back("7, 18");
expectedResult.push_back("8, 16");
expectedResult.push_back("9, 14");
result = checker->check(this->env(), tasks[taskId++]);
compare = this->compareResult(model, result, expectedResult);
EXPECT_TRUE(compare.first) << compare.second;
}
TYPED_TEST(QuantileQueryTest, simple_Mdp) {
typedef storm::models::sparse::Mdp<typename TestFixture::ValueType> ModelType;
std::string formulasString = "quantile(B1, B2, Pmax>0.5 [F{\"first\"}>=B1,{\"second\"}>=B2 s=1]);\n";
formulasString += "quantile(B1, B2, Pmax>0.5 [F{\"first\"}<=B1,{\"second\"}<=B2 s=1]);\n";
formulasString += "quantile(B1, B2, Pmin>0.5 [F{\"first\"}<=B1,{\"second\"}<=B2 s=1]);\n";
formulasString += "quantile(B1, Pmax>0.5 [F{\"second\"}>=B1 s=1]);\n";
formulasString += "quantile(B1, B2, B3, Pmax>0.5 [F{\"first\"}<=B1,{\"second\"}<=B2,{\"third\"}<=B3 s=1]);\n";
auto modelFormulas = this->template buildModelFormulas<ModelType>(STORM_TEST_RESOURCES_DIR "/mdp/quantiles_simple_mdp.nm", formulasString);
auto model = std::move(modelFormulas.first);
auto tasks = this->getTasks(modelFormulas.second);
EXPECT_EQ(7ul, model->getNumberOfStates());
EXPECT_EQ(13ul, model->getNumberOfTransitions());
auto checker = this->template createModelChecker<ModelType>(model);
std::unique_ptr<storm::modelchecker::CheckResult> result;
std::vector<std::string> expectedResult;
std::pair<bool, std::string> compare;
uint64_t taskId = 0;
expectedResult.clear();
expectedResult.push_back(" 0, 6");
expectedResult.push_back(" 1, 5");
expectedResult.push_back(" 2, 4");
expectedResult.push_back(" 3, 3");
expectedResult.push_back("inf, 2");
result = checker->check(this->env(), tasks[taskId++]);
compare = this->compareResult(model, result, expectedResult);
EXPECT_TRUE(compare.first) << compare.second;
expectedResult.clear();
expectedResult.push_back(" 0, 2");
expectedResult.push_back(" 5, 1");
expectedResult.push_back(" 6, 0");
result = checker->check(this->env(), tasks[taskId++]);
compare = this->compareResult(model, result, expectedResult);
EXPECT_TRUE(compare.first) << compare.second;
expectedResult.clear();
result = checker->check(this->env(), tasks[taskId++]);
compare = this->compareResult(model, result, expectedResult);
EXPECT_TRUE(compare.first) << compare.second;
expectedResult.clear();
expectedResult.push_back(" 6");
result = checker->check(this->env(), tasks[taskId++]);
compare = this->compareResult(model, result, expectedResult);
EXPECT_TRUE(compare.first) << compare.second;
expectedResult.clear();
expectedResult.push_back(" 0, 2, 1.4");
expectedResult.push_back(" 5, 1, 9.8");
expectedResult.push_back(" 6, 0, 9.8");
result = checker->check(this->env(), tasks[taskId++]);
compare = this->compareResult(model, result, expectedResult);
EXPECT_TRUE(compare.first) << compare.second;
}
TYPED_TEST(QuantileQueryTest, firewire) {
typedef storm::models::sparse::Mdp<typename TestFixture::ValueType> ModelType;
std::string formulasString = "quantile(min TIME, min ROUNDS, Pmax>0.95 [F{\"time\"}<=TIME,{\"rounds\"}<=ROUNDS \"done\"]);\n";
auto modelFormulas = this->template buildModelFormulas<ModelType>(STORM_TEST_RESOURCES_DIR "/mdp/quantiles_firewire.nm", formulasString, "delay=36");
auto model = std::move(modelFormulas.first);
auto tasks = this->getTasks(modelFormulas.second);
EXPECT_EQ(776ul, model->getNumberOfStates());
EXPECT_EQ(1411ul, model->getNumberOfTransitions());
auto checker = this->template createModelChecker<ModelType>(model);
std::unique_ptr<storm::modelchecker::CheckResult> result;
std::vector<std::string> expectedResult;
std::pair<bool, std::string> compare;
uint64_t taskId = 0;
expectedResult.clear();
expectedResult.push_back("123, 1");
result = checker->check(this->env(), tasks[taskId++]);
compare = this->compareResult(model, result, expectedResult);
EXPECT_TRUE(compare.first) << compare.second;
}
TYPED_TEST(QuantileQueryTest, resources) {
typedef storm::models::sparse::Mdp<typename TestFixture::ValueType> ModelType;
std::string formulasString = "quantile(max GOLD, max GEM, Pmax>0.95 [F{\"gold\"}>=GOLD,{\"gem\"}>=GEM,{\"steps\"}<=100 true]);\n";
auto modelFormulas = this->template buildModelFormulas<ModelType>(STORM_TEST_RESOURCES_DIR "/mdp/quantiles_resources.nm", formulasString);
auto model = std::move(modelFormulas.first);
auto tasks = this->getTasks(modelFormulas.second);
EXPECT_EQ(94ul, model->getNumberOfStates());
EXPECT_EQ(326ul, model->getNumberOfTransitions());
auto checker = this->template createModelChecker<ModelType>(model);
std::unique_ptr<storm::modelchecker::CheckResult> result;
std::vector<std::string> expectedResult;
std::pair<bool, std::string> compare;
uint64_t taskId = 0;
expectedResult.clear();
expectedResult.push_back("0, 10");
expectedResult.push_back("1, 9");
expectedResult.push_back("4, 8");
expectedResult.push_back("7, 7");
expectedResult.push_back("8, 4");
expectedResult.push_back("9, 2");
expectedResult.push_back("10, 0");
result = checker->check(this->env(), tasks[taskId++]);
compare = this->compareResult(model, result, expectedResult);
EXPECT_TRUE(compare.first) << compare.second;
}
}
|||||||
100:0
Loading…
Cancel
Save