diff --git a/resources/examples/testfiles/mdp/cvar_bad_mec_mdp.nm b/resources/examples/testfiles/mdp/cvar_bad_mec_mdp.nm new file mode 100644 index 0000000000..529ef1306a --- /dev/null +++ b/resources/examples/testfiles/mdp/cvar_bad_mec_mdp.nm @@ -0,0 +1,16 @@ +mdp + +module main + s : [0..3] init 0; + + [] s=0 -> 1/2 : (s'=1) + 1/2 : (s'=2); + [] s=1 -> 1 : (s'=1); + [] s=2 -> 1 : (s'=3); + [] s=3 -> 1 : (s'=2); +endmodule + +label "target" = s=1; + +rewards "term" + s=1 : 4; +endrewards diff --git a/resources/examples/testfiles/mdp/cvar_branching_tradeoff_mdp.nm b/resources/examples/testfiles/mdp/cvar_branching_tradeoff_mdp.nm new file mode 100644 index 0000000000..2c2738887f --- /dev/null +++ b/resources/examples/testfiles/mdp/cvar_branching_tradeoff_mdp.nm @@ -0,0 +1,41 @@ +mdp + +module main + s : [0..10] init 0; + + // Start: direct safe/balanced/risky options, plus an adaptive branch. + [safe] s=0 -> 1 : (s'=3); + [balanced] s=0 -> 3/4 : (s'=1) + 1/4 : (s'=5); + [adaptive] s=0 -> 1/2 : (s'=8) + 1/2 : (s'=9); + [risky] s=0 -> 1/2 : (s'=10) + 1/2 : (s'=7); + + // Duplicate terminal rewards at 6 and 12 make threshold splitting observable. + [] s=1 -> 1 : (s'=1); + [] s=2 -> 1 : (s'=2); + [] s=3 -> 1 : (s'=3); + [] s=4 -> 1 : (s'=4); + [] s=5 -> 1 : (s'=5); + [] s=6 -> 1 : (s'=6); + [] s=7 -> 1 : (s'=7); + + // After the adaptive probabilistic split, the scheduler can react to good/bad news. + [cash] s=8 -> 1 : (s'=4); + [push] s=8 -> 1/2 : (s'=6) + 1/2 : (s'=7); + + [cash] s=9 -> 1 : (s'=2); + [push] s=9 -> 1/2 : (s'=10) + 1/2 : (s'=4); + + [] s=10 -> 1 : (s'=10); +endmodule + +label "target" = s=1 | s=2 | s=3 | s=4 | s=5 | s=6 | s=7 | s=10; + +rewards "term" + s=1 : 6; + s=2 : 6; + s=3 : 7; + s=4 : 10; + s=5 : 12; + s=6 : 12; + s=7 : 20; +endrewards diff --git a/resources/examples/testfiles/mdp/cvar_nonabsorbing_target_mdp.nm b/resources/examples/testfiles/mdp/cvar_nonabsorbing_target_mdp.nm new file mode 100644 index 0000000000..7dd60967b8 --- /dev/null +++ b/resources/examples/testfiles/mdp/cvar_nonabsorbing_target_mdp.nm @@ -0,0 +1,15 @@ +mdp + +module main + s : [0..2] init 0; + + [] s=0 -> 1 : (s'=1); + [] s=1 -> 1 : (s'=2); + [] s=2 -> 1 : (s'=2); +endmodule + +label "target" = s=1; + +rewards "term" + s=1 : 5; +endrewards diff --git a/resources/examples/testfiles/mdp/cvar_simple_mdp.nm b/resources/examples/testfiles/mdp/cvar_simple_mdp.nm new file mode 100644 index 0000000000..cbdd275071 --- /dev/null +++ b/resources/examples/testfiles/mdp/cvar_simple_mdp.nm @@ -0,0 +1,19 @@ +mdp + +module main + s : [0..3] init 0; + + [] s=0 -> 1/2 : (s'=1) + 1/2 : (s'=2); + [] s=0 -> 1 : (s'=3); + [] s=1 -> 1 : (s'=1); + [] s=2 -> 1 : (s'=2); + [] s=3 -> 1 : (s'=3); +endmodule + +label "target" = s=1 | s=2 | s=3; + +rewards "term" + s=1 : 1; + s=2 : 3; + s=3 : 2; +endrewards diff --git a/resources/examples/testfiles/mdp/cvar_ssp_branching_tradeoff_mdp.nm b/resources/examples/testfiles/mdp/cvar_ssp_branching_tradeoff_mdp.nm new file mode 100644 index 0000000000..4471601b7d --- /dev/null +++ b/resources/examples/testfiles/mdp/cvar_ssp_branching_tradeoff_mdp.nm @@ -0,0 +1,21 @@ + +mdp + +module main + s : [0..3] init 0; + + [safe] s=0 -> 1 : (s'=3); + [risky] s=0 -> 1/2 : (s'=1) + 1/2 : (s'=2); + [low] s=1 -> 1 : (s'=3); + [high] s=2 -> 1 : (s'=3); + [] s=3 -> 1 : (s'=3); +endmodule + +label "goal" = s=3; + +rewards "cost" + [safe] true : 6; + [risky] true : 1; + [low] true : 1; + [high] true : 7; +endrewards diff --git a/resources/examples/testfiles/mdp/cvar_ssp_deterministic_mdp.nm b/resources/examples/testfiles/mdp/cvar_ssp_deterministic_mdp.nm new file mode 100644 index 0000000000..4fbbc8dcfe --- /dev/null +++ b/resources/examples/testfiles/mdp/cvar_ssp_deterministic_mdp.nm @@ -0,0 +1,16 @@ +mdp + +module main + s : [0..2] init 0; + + [step0] s=0 -> 1 : (s'=1); + [step1] s=1 -> 1 : (s'=2); + [] s=2 -> 1 : (s'=2); +endmodule + +label "goal" = s=2; + +rewards "cost" + [step0] true : 2; + [step1] true : 3; +endrewards diff --git a/resources/examples/testfiles/mdp/cvar_ssp_reward_delayed_challenger_mdp.nm b/resources/examples/testfiles/mdp/cvar_ssp_reward_delayed_challenger_mdp.nm new file mode 100644 index 0000000000..9eaec9708d --- /dev/null +++ b/resources/examples/testfiles/mdp/cvar_ssp_reward_delayed_challenger_mdp.nm @@ -0,0 +1,20 @@ +mdp + +module main + s : [0..3] init 0; + + [det] s=0 -> 1 : (s'=3); + [lottery] s=0 -> 1/2 : (s'=1) + 1/2 : (s'=2); + [low] s=1 -> 1 : (s'=3); + [high] s=2 -> 1 : (s'=3); + [] s=3 -> 1 : (s'=3); +endmodule + +label "goal" = s=3; + +rewards "reward" + [det] true : 6; + [lottery] true : 1; + [low] true : 7; + [high] true : 99; +endrewards diff --git a/resources/examples/testfiles/mdp/cvar_ssp_reward_geometric_mdp.nm b/resources/examples/testfiles/mdp/cvar_ssp_reward_geometric_mdp.nm new file mode 100644 index 0000000000..fccc3fde34 --- /dev/null +++ b/resources/examples/testfiles/mdp/cvar_ssp_reward_geometric_mdp.nm @@ -0,0 +1,14 @@ +mdp + +module main + s : [0..1] init 0; + + [step] s=0 -> 1/2 : (s'=0) + 1/2 : (s'=1); + [] s=1 -> 1 : (s'=1); +endmodule + +label "goal" = s=1; + +rewards "reward" + [step] true : 1; +endrewards diff --git a/resources/examples/testfiles/mdp/cvar_ssp_reward_prob1a_choice_violation_mdp.nm b/resources/examples/testfiles/mdp/cvar_ssp_reward_prob1a_choice_violation_mdp.nm new file mode 100644 index 0000000000..dafa87e386 --- /dev/null +++ b/resources/examples/testfiles/mdp/cvar_ssp_reward_prob1a_choice_violation_mdp.nm @@ -0,0 +1,16 @@ +mdp + +module main + s : [0..1] init 0; + + [loop] s=0 -> 1 : (s'=0); + [exit] s=0 -> 1 : (s'=1); + [] s=1 -> 1 : (s'=1); +endmodule + +label "goal" = s=1; + +rewards "reward" + [loop] true : 1; + [exit] true : 1; +endrewards diff --git a/resources/examples/testfiles/mdp/cvar_ssp_reward_prob1a_probabilistic_violation_mdp.nm b/resources/examples/testfiles/mdp/cvar_ssp_reward_prob1a_probabilistic_violation_mdp.nm new file mode 100644 index 0000000000..293c61f115 --- /dev/null +++ b/resources/examples/testfiles/mdp/cvar_ssp_reward_prob1a_probabilistic_violation_mdp.nm @@ -0,0 +1,16 @@ +mdp + +module main + s : [0..2] init 0; + + [split] s=0 -> 1/2 : (s'=1) + 1/2 : (s'=2); + [loop] s=1 -> 1 : (s'=1); + [] s=2 -> 1 : (s'=2); +endmodule + +label "goal" = s=2; + +rewards "reward" + [split] true : 1; + [loop] true : 1; +endrewards diff --git a/resources/examples/testfiles/mdp/cvar_ssp_reward_safe_risky_mdp.nm b/resources/examples/testfiles/mdp/cvar_ssp_reward_safe_risky_mdp.nm new file mode 100644 index 0000000000..a4c957f93d --- /dev/null +++ b/resources/examples/testfiles/mdp/cvar_ssp_reward_safe_risky_mdp.nm @@ -0,0 +1,20 @@ +mdp + +module main + s : [0..3] init 0; + + [safe] s=0 -> 1 : (s'=3); + [risky] s=0 -> 1/5 : (s'=1) + 4/5 : (s'=2); + [low] s=1 -> 1 : (s'=3); + [high] s=2 -> 1 : (s'=3); + [] s=3 -> 1 : (s'=3); +endmodule + +label "goal" = s=3; + +rewards "reward" + [safe] true : 6; + [risky] true : 1; + [low] true : 1; + [high] true : 9; +endrewards diff --git a/resources/examples/testfiles/mdp/cvar_target_reaching_mec_mdp.nm b/resources/examples/testfiles/mdp/cvar_target_reaching_mec_mdp.nm new file mode 100644 index 0000000000..3563c87290 --- /dev/null +++ b/resources/examples/testfiles/mdp/cvar_target_reaching_mec_mdp.nm @@ -0,0 +1,23 @@ +mdp + +module main + s : [0..4] init 0; + + [] s=0 -> 1 : (s'=1); + + [cycle12] s=1 -> 1 : (s'=2); + [takeLow] s=1 -> 1 : (s'=3); + + [cycle21] s=2 -> 1 : (s'=1); + [takeRisk] s=2 -> 1/2 : (s'=3) + 1/2 : (s'=4); + + [] s=3 -> 1 : (s'=3); + [] s=4 -> 1 : (s'=4); +endmodule + +label "target" = s=3 | s=4; + +rewards "term" + s=3 : 4; + s=4 : 10; +endrewards diff --git a/src/storm-cli-utilities/model-handling.h b/src/storm-cli-utilities/model-handling.h index 64962fb2c1..7b8a1ff396 100644 --- a/src/storm-cli-utilities/model-handling.h +++ b/src/storm-cli-utilities/model-handling.h @@ -432,6 +432,9 @@ inline std::pair preprocessSymbolicIn SymbolicInput output = input; // Preprocess properties (if requested) + STORM_LOG_THROW(!(ioSettings.isPropertiesAsMultiSet() && ioSettings.isCvarSet()), storm::exceptions::InvalidArgumentException, + "Options '--propsasmulti' and '--cvar' can not be combined."); + if (ioSettings.isPropertiesAsMultiSet()) { STORM_LOG_THROW(!input.properties.empty(), storm::exceptions::InvalidArgumentException, "Can not translate properties to multi-objective formula because no properties were specified."); @@ -440,6 +443,14 @@ inline std::pair preprocessSymbolicIn output.properties = {storm::api::createMultiObjectiveProperty(output.properties, false)}; } + if (ioSettings.isCvarSet()) { + STORM_LOG_THROW(!input.properties.empty(), storm::exceptions::InvalidArgumentException, + "Can not translate properties to a CVaR formula because no properties were specified."); + STORM_LOG_THROW(output.properties.size() == 1, storm::exceptions::InvalidArgumentException, + "The '--cvar' option currently requires exactly one selected property."); + output.properties = {storm::api::createCvarProperty(output.properties.front(), ioSettings.getCvarAlpha())}; + } + // Substitute constant definitions in symbolic input. std::string constantDefinitionString = ioSettings.getConstantDefinitionString(); std::map constantDefinitions; diff --git a/src/storm/api/properties.cpp b/src/storm/api/properties.cpp index 70c23f05c2..9da0f9e686 100644 --- a/src/storm/api/properties.cpp +++ b/src/storm/api/properties.cpp @@ -2,14 +2,31 @@ #include -#include "storm/logic/Formula.h" +#include "storm/exceptions/InvalidArgumentException.h" +#include "storm/logic/Formulas.h" #include "storm/storage/SymbolicModelDescription.h" #include "storm/storage/jani/Model.h" #include "storm/storage/jani/Property.h" #include "storm/storage/prism/Program.h" +#include "storm/utility/constants.h" +#include "storm/utility/macros.h" namespace storm { namespace api { +namespace { + +storm::RationalNumber parseCvarAlpha(std::string const& input) { + std::string strippedInput = boost::algorithm::trim_copy(input); + STORM_LOG_THROW(!strippedInput.empty(), storm::exceptions::InvalidArgumentException, "Unable to parse CVaR alpha '" << input << "'."); + + storm::RationalNumber alpha = storm::utility::convertNumber(strippedInput); + + STORM_LOG_THROW(storm::utility::zero() < alpha && alpha < storm::utility::one(), + storm::exceptions::InvalidArgumentException, "The CVaR alpha must be in the open interval (0, 1)."); + return alpha; +} + +} // namespace std::vector substituteConstantsInProperties(std::vector const& properties, std::map const& substitution) { @@ -67,6 +84,17 @@ std::vector> extractFormulasFromPro return formulas; } +storm::jani::Property createCvarProperty(storm::jani::Property const& property, storm::RationalNumber const& alpha) { + STORM_LOG_WARN_COND(property.getFilter().isDefault(), + "Non-default property filter of property " << property.getName() << " will be dropped during conversion to CVaR property."); + auto cvarFormula = std::make_shared(alpha, property.getRawFormula()); + return storm::jani::Property(property.getName(), cvarFormula, property.getUndefinedConstants(), property.getComment()); +} + +storm::jani::Property createCvarProperty(storm::jani::Property const& property, std::string const& alpha) { + return createCvarProperty(property, parseCvarAlpha(alpha)); +} + storm::jani::Property createMultiObjectiveProperty(std::vector const& properties, bool lexicographic) { std::set undefConstants; std::string name = ""; diff --git a/src/storm/api/properties.h b/src/storm/api/properties.h index a852b523ed..c2dbd19394 100644 --- a/src/storm/api/properties.h +++ b/src/storm/api/properties.h @@ -7,6 +7,8 @@ #include #include +#include "storm/adapters/RationalNumberAdapter.h" + namespace storm { namespace jani { @@ -36,6 +38,8 @@ std::vector substituteTranscendentalNumbersInProperties(s std::vector filterProperties(std::vector const& properties, boost::optional> const& propertyFilter); std::vector> extractFormulasFromProperties(std::vector const& properties); +storm::jani::Property createCvarProperty(storm::jani::Property const& property, storm::RationalNumber const& alpha); +storm::jani::Property createCvarProperty(storm::jani::Property const& property, std::string const& alpha); storm::jani::Property createMultiObjectiveProperty(std::vector const& properties, bool lexicographic); } // namespace api diff --git a/src/storm/environment/SubEnvironment.cpp b/src/storm/environment/SubEnvironment.cpp index 35a23d35bb..eba49290ab 100644 --- a/src/storm/environment/SubEnvironment.cpp +++ b/src/storm/environment/SubEnvironment.cpp @@ -48,6 +48,7 @@ void SubEnvironment::assertInitialized() const { template class SubEnvironment; +template class SubEnvironment; template class SubEnvironment; template class SubEnvironment; template class SubEnvironment; diff --git a/src/storm/environment/modelchecker/AllModelCheckerEnvironments.h b/src/storm/environment/modelchecker/AllModelCheckerEnvironments.h index 30a6735b27..c1d068ae64 100644 --- a/src/storm/environment/modelchecker/AllModelCheckerEnvironments.h +++ b/src/storm/environment/modelchecker/AllModelCheckerEnvironments.h @@ -1,5 +1,6 @@ #pragma once #include "storm/environment/modelchecker/ConditionalModelCheckerEnvironment.h" +#include "storm/environment/modelchecker/CvarModelCheckerEnvironment.h" #include "storm/environment/modelchecker/ModelCheckerEnvironment.h" -#include "storm/environment/modelchecker/MultiObjectiveModelCheckerEnvironment.h" \ No newline at end of file +#include "storm/environment/modelchecker/MultiObjectiveModelCheckerEnvironment.h" diff --git a/src/storm/environment/modelchecker/CvarModelCheckerEnvironment.cpp b/src/storm/environment/modelchecker/CvarModelCheckerEnvironment.cpp new file mode 100644 index 0000000000..5046fea488 --- /dev/null +++ b/src/storm/environment/modelchecker/CvarModelCheckerEnvironment.cpp @@ -0,0 +1,34 @@ +#include "storm/environment/modelchecker/CvarModelCheckerEnvironment.h" + +#include "storm/settings/SettingsManager.h" +#include "storm/settings/modules/CvarSettings.h" + +namespace storm { + +CvarModelCheckerEnvironment::CvarModelCheckerEnvironment() { + auto const& cvarSettings = storm::settings::getModule(); + method = cvarSettings.getCvarMethod(); + interpretationSelection = cvarSettings.getInterpretationSelection(); +} + +CvarModelCheckerEnvironment::~CvarModelCheckerEnvironment() { + // Intentionally left empty +} + +storm::modelchecker::cvar::CvarMethod const& CvarModelCheckerEnvironment::getMethod() const { + return method; +} + +void CvarModelCheckerEnvironment::setMethod(storm::modelchecker::cvar::CvarMethod value) { + method = value; +} + +storm::modelchecker::cvar::CvarInterpretationSelection const& CvarModelCheckerEnvironment::getInterpretationSelection() const { + return interpretationSelection; +} + +void CvarModelCheckerEnvironment::setInterpretationSelection(storm::modelchecker::cvar::CvarInterpretationSelection value) { + interpretationSelection = value; +} + +} // namespace storm diff --git a/src/storm/environment/modelchecker/CvarModelCheckerEnvironment.h b/src/storm/environment/modelchecker/CvarModelCheckerEnvironment.h new file mode 100644 index 0000000000..ffea10b5f9 --- /dev/null +++ b/src/storm/environment/modelchecker/CvarModelCheckerEnvironment.h @@ -0,0 +1,23 @@ +#pragma once + +#include "storm/modelchecker/cvar/CvarInterpretation.h" +#include "storm/modelchecker/cvar/CvarMethod.h" + +namespace storm { + +class CvarModelCheckerEnvironment { + public: + CvarModelCheckerEnvironment(); + ~CvarModelCheckerEnvironment(); + + storm::modelchecker::cvar::CvarMethod const& getMethod() const; + void setMethod(storm::modelchecker::cvar::CvarMethod value); + storm::modelchecker::cvar::CvarInterpretationSelection const& getInterpretationSelection() const; + void setInterpretationSelection(storm::modelchecker::cvar::CvarInterpretationSelection value); + + private: + storm::modelchecker::cvar::CvarMethod method; + storm::modelchecker::cvar::CvarInterpretationSelection interpretationSelection; +}; + +} // namespace storm diff --git a/src/storm/environment/modelchecker/ModelCheckerEnvironment.cpp b/src/storm/environment/modelchecker/ModelCheckerEnvironment.cpp index bed60a93bc..16d4a22ae1 100644 --- a/src/storm/environment/modelchecker/ModelCheckerEnvironment.cpp +++ b/src/storm/environment/modelchecker/ModelCheckerEnvironment.cpp @@ -1,6 +1,7 @@ #include "storm/environment/modelchecker/ModelCheckerEnvironment.h" #include "storm/environment/modelchecker/ConditionalModelCheckerEnvironment.h" +#include "storm/environment/modelchecker/CvarModelCheckerEnvironment.h" #include "storm/environment/modelchecker/MultiObjectiveModelCheckerEnvironment.h" #include "storm/settings/SettingsManager.h" @@ -26,6 +27,14 @@ ModelCheckerEnvironment::~ModelCheckerEnvironment() { // Intentionally left empty } +CvarModelCheckerEnvironment& ModelCheckerEnvironment::cvar() { + return cvarModelCheckerEnvironment.get(); +} + +CvarModelCheckerEnvironment const& ModelCheckerEnvironment::cvar() const { + return cvarModelCheckerEnvironment.get(); +} + ConditionalModelCheckerEnvironment& ModelCheckerEnvironment::conditional() { return conditionalModelCheckerEnvironment.get(); } diff --git a/src/storm/environment/modelchecker/ModelCheckerEnvironment.h b/src/storm/environment/modelchecker/ModelCheckerEnvironment.h index 983d3d37a8..7ae5c3b518 100644 --- a/src/storm/environment/modelchecker/ModelCheckerEnvironment.h +++ b/src/storm/environment/modelchecker/ModelCheckerEnvironment.h @@ -5,6 +5,7 @@ #include "storm/environment/Environment.h" #include "storm/environment/SubEnvironment.h" +#include "storm/environment/modelchecker/CvarModelCheckerEnvironment.h" #include "storm/modelchecker/helper/infinitehorizon/SteadyStateDistributionAlgorithm.h" namespace storm { @@ -18,6 +19,9 @@ class ModelCheckerEnvironment { ModelCheckerEnvironment(); ~ModelCheckerEnvironment(); + CvarModelCheckerEnvironment& cvar(); + CvarModelCheckerEnvironment const& cvar() const; + ConditionalModelCheckerEnvironment& conditional(); ConditionalModelCheckerEnvironment const& conditional() const; @@ -33,6 +37,7 @@ class ModelCheckerEnvironment { void unsetLtl2daTool(); private: + SubEnvironment cvarModelCheckerEnvironment; SubEnvironment conditionalModelCheckerEnvironment; SubEnvironment multiObjectiveModelCheckerEnvironment; boost::optional ltl2daTool; diff --git a/src/storm/logic/CloneVisitor.cpp b/src/storm/logic/CloneVisitor.cpp index f3531aa60c..a0141c0add 100644 --- a/src/storm/logic/CloneVisitor.cpp +++ b/src/storm/logic/CloneVisitor.cpp @@ -76,6 +76,11 @@ boost::any CloneVisitor::visit(CumulativeRewardFormula const& f, boost::any cons return std::static_pointer_cast(std::make_shared(f)); } +boost::any CloneVisitor::visit(CvarFormula const& f, boost::any const& data) const { + std::shared_ptr subformula = boost::any_cast>(f.getSubformula().accept(*this, data)); + return std::static_pointer_cast(std::make_shared(f.getAlpha(), subformula)); +} + boost::any CloneVisitor::visit(EventuallyFormula const& f, boost::any const& data) const { std::shared_ptr subformula = boost::any_cast>(f.getSubformula().accept(*this, data)); if (f.hasRewardAccumulation()) { diff --git a/src/storm/logic/CloneVisitor.h b/src/storm/logic/CloneVisitor.h index c9eefef784..ad7b049cfd 100644 --- a/src/storm/logic/CloneVisitor.h +++ b/src/storm/logic/CloneVisitor.h @@ -19,6 +19,7 @@ class CloneVisitor : public FormulaVisitor { virtual boost::any visit(BoundedUntilFormula const& f, boost::any const& data) const override; virtual boost::any visit(ConditionalFormula const& f, boost::any const& data) const override; virtual boost::any visit(CumulativeRewardFormula const& f, boost::any const& data) const override; + virtual boost::any visit(CvarFormula const& f, boost::any const& data) const override; virtual boost::any visit(EventuallyFormula const& f, boost::any const& data) const override; virtual boost::any visit(TimeOperatorFormula const& f, boost::any const& data) const override; virtual boost::any visit(GloballyFormula const& f, boost::any const& data) const override; diff --git a/src/storm/logic/CvarFormula.cpp b/src/storm/logic/CvarFormula.cpp new file mode 100644 index 0000000000..5853aaf588 --- /dev/null +++ b/src/storm/logic/CvarFormula.cpp @@ -0,0 +1,77 @@ +#include "storm/logic/CvarFormula.h" + +#include +#include +#include + +#include "storm/exceptions/InvalidArgumentException.h" +#include "storm/logic/FormulaVisitor.h" +#include "storm/utility/constants.h" +#include "storm/utility/macros.h" + +namespace storm { +namespace logic { + +CvarFormula::CvarFormula(storm::RationalNumber const& alpha, std::shared_ptr subformula) : alpha(alpha), subformula(std::move(subformula)) { + STORM_LOG_THROW(this->subformula != nullptr, storm::exceptions::InvalidArgumentException, "A CVaR formula requires a subformula."); + STORM_LOG_THROW(storm::utility::zero() < this->alpha && this->alpha < storm::utility::one(), + storm::exceptions::InvalidArgumentException, "The CVaR alpha must be in the open interval (0, 1)."); +} + +CvarFormula::~CvarFormula() { + // Intentionally left empty. +} + +bool CvarFormula::isCvarFormula() const { + return true; +} + +bool CvarFormula::hasQuantitativeResult() const { + return true; +} + +bool CvarFormula::hasNumericalResult() const { + return true; +} + +bool CvarFormula::hasMultiDimensionalResult() const { + return false; +} + +storm::RationalNumber const& CvarFormula::getAlpha() const { + return alpha; +} + +Formula const& CvarFormula::getSubformula() const { + return *subformula; +} + +boost::any CvarFormula::accept(FormulaVisitor const& visitor, boost::any const& data) const { + return visitor.visit(*this, data); +} + +void CvarFormula::gatherAtomicExpressionFormulas(std::vector>& atomicExpressionFormulas) const { + subformula->gatherAtomicExpressionFormulas(atomicExpressionFormulas); +} + +void CvarFormula::gatherAtomicLabelFormulas(std::vector>& atomicLabelFormulas) const { + subformula->gatherAtomicLabelFormulas(atomicLabelFormulas); +} + +void CvarFormula::gatherReferencedRewardModels(std::set& referencedRewardModels) const { + subformula->gatherReferencedRewardModels(referencedRewardModels); +} + +void CvarFormula::gatherUsedVariables(std::set& usedVariables) const { + subformula->gatherUsedVariables(usedVariables); +} + +std::ostream& CvarFormula::writeToStream(std::ostream& out, bool /* allowParentheses */) const { + out << "cvar(" << alpha << ", "; + subformula->writeToStream(out); + out << ")"; + return out; +} + +} // namespace logic +} // namespace storm diff --git a/src/storm/logic/CvarFormula.h b/src/storm/logic/CvarFormula.h new file mode 100644 index 0000000000..e4e66c2f41 --- /dev/null +++ b/src/storm/logic/CvarFormula.h @@ -0,0 +1,38 @@ +#pragma once + +#include "storm/adapters/RationalNumberAdapter.h" +#include "storm/logic/StateFormula.h" + +namespace storm { +namespace logic { + +class CvarFormula : public StateFormula { + public: + CvarFormula(storm::RationalNumber const& alpha, std::shared_ptr subformula); + + virtual ~CvarFormula(); + + virtual bool isCvarFormula() const override; + + virtual bool hasQuantitativeResult() const override; + virtual bool hasNumericalResult() const; + virtual bool hasMultiDimensionalResult() const; + + storm::RationalNumber const& getAlpha() const; + Formula const& getSubformula() const; + + virtual boost::any accept(FormulaVisitor const& visitor, boost::any const& data) const override; + virtual void gatherAtomicExpressionFormulas(std::vector>& atomicExpressionFormulas) const override; + virtual void gatherAtomicLabelFormulas(std::vector>& atomicLabelFormulas) const override; + virtual void gatherReferencedRewardModels(std::set& referencedRewardModels) const override; + virtual void gatherUsedVariables(std::set& usedVariables) const override; + + virtual std::ostream& writeToStream(std::ostream& out, bool allowParentheses = false) const override; + + private: + storm::RationalNumber alpha; + std::shared_ptr subformula; +}; + +} // namespace logic +} // namespace storm diff --git a/src/storm/logic/ExtractMaximalStateFormulasVisitor.cpp b/src/storm/logic/ExtractMaximalStateFormulasVisitor.cpp index 571bb66f5d..fe6e2c2413 100644 --- a/src/storm/logic/ExtractMaximalStateFormulasVisitor.cpp +++ b/src/storm/logic/ExtractMaximalStateFormulasVisitor.cpp @@ -164,6 +164,13 @@ boost::any ExtractMaximalStateFormulasVisitor::visit(MultiObjectiveFormula const return result; } +boost::any ExtractMaximalStateFormulasVisitor::visit(CvarFormula const& f, boost::any const& data) const { + incrementNestingLevel(); + boost::any result = CloneVisitor::visit(f, data); + decrementNestingLevel(); + return result; +} + boost::any ExtractMaximalStateFormulasVisitor::visit(ProbabilityOperatorFormula const& f, boost::any const& data) const { incrementNestingLevel(); boost::any result = CloneVisitor::visit(f, data); diff --git a/src/storm/logic/ExtractMaximalStateFormulasVisitor.h b/src/storm/logic/ExtractMaximalStateFormulasVisitor.h index 82bffc792d..6723199a6c 100644 --- a/src/storm/logic/ExtractMaximalStateFormulasVisitor.h +++ b/src/storm/logic/ExtractMaximalStateFormulasVisitor.h @@ -30,6 +30,7 @@ class ExtractMaximalStateFormulasVisitor : public CloneVisitor { virtual boost::any visit(TimeOperatorFormula const& f, boost::any const& data) const override; virtual boost::any visit(LongRunAverageOperatorFormula const& f, boost::any const& data) const override; virtual boost::any visit(MultiObjectiveFormula const& f, boost::any const& data) const override; + virtual boost::any visit(CvarFormula const& f, boost::any const& data) const override; virtual boost::any visit(ProbabilityOperatorFormula const& f, boost::any const& data) const override; virtual boost::any visit(RewardOperatorFormula const& f, boost::any const& data) const override; diff --git a/src/storm/logic/Formula.cpp b/src/storm/logic/Formula.cpp index 440d9139d4..49de2b7300 100644 --- a/src/storm/logic/Formula.cpp +++ b/src/storm/logic/Formula.cpp @@ -33,6 +33,10 @@ bool Formula::isQuantileFormula() const { return false; } +bool Formula::isCvarFormula() const { + return false; +} + bool Formula::isBinaryStateFormula() const { return false; } @@ -250,6 +254,14 @@ QuantileFormula const& Formula::asQuantileFormula() const { return dynamic_cast(*this); } +CvarFormula& Formula::asCvarFormula() { + return dynamic_cast(*this); +} + +CvarFormula const& Formula::asCvarFormula() const { + return dynamic_cast(*this); +} + BinaryStateFormula& Formula::asBinaryStateFormula() { return dynamic_cast(*this); } diff --git a/src/storm/logic/Formula.h b/src/storm/logic/Formula.h index c6a4a7f771..15274c5286 100644 --- a/src/storm/logic/Formula.h +++ b/src/storm/logic/Formula.h @@ -55,6 +55,7 @@ class Formula : public std::enable_shared_from_this { virtual bool isMultiObjectiveFormula() const; virtual bool isQuantileFormula() const; + virtual bool isCvarFormula() const; // Operator formulas. virtual bool isOperatorFormula() const; @@ -127,6 +128,9 @@ class Formula : public std::enable_shared_from_this { QuantileFormula& asQuantileFormula(); QuantileFormula const& asQuantileFormula() const; + CvarFormula& asCvarFormula(); + CvarFormula const& asCvarFormula() const; + BinaryStateFormula& asBinaryStateFormula(); BinaryStateFormula const& asBinaryStateFormula() const; diff --git a/src/storm/logic/FormulaInformationVisitor.cpp b/src/storm/logic/FormulaInformationVisitor.cpp index 61eef9e4a4..1b8d4ef9ab 100644 --- a/src/storm/logic/FormulaInformationVisitor.cpp +++ b/src/storm/logic/FormulaInformationVisitor.cpp @@ -79,6 +79,14 @@ boost::any FormulaInformationVisitor::visit(CumulativeRewardFormula const& f, bo return result; } +boost::any FormulaInformationVisitor::visit(CvarFormula const& f, boost::any const& data) const { + if (recurseIntoOperators) { + return f.getSubformula().accept(*this, data); + } else { + return FormulaInformation(); + } +} + boost::any FormulaInformationVisitor::visit(EventuallyFormula const& f, boost::any const& data) const { FormulaInformation result = boost::any_cast(f.getSubformula().accept(*this, data)); if (f.getSubformula().isPathFormula()) { diff --git a/src/storm/logic/FormulaInformationVisitor.h b/src/storm/logic/FormulaInformationVisitor.h index 287c25d3d4..f3470a1ac4 100644 --- a/src/storm/logic/FormulaInformationVisitor.h +++ b/src/storm/logic/FormulaInformationVisitor.h @@ -24,6 +24,7 @@ class FormulaInformationVisitor : public FormulaVisitor { virtual boost::any visit(BoundedUntilFormula const& f, boost::any const& data) const override; virtual boost::any visit(ConditionalFormula const& f, boost::any const& data) const override; virtual boost::any visit(CumulativeRewardFormula const& f, boost::any const& data) const override; + virtual boost::any visit(CvarFormula const& f, boost::any const& data) const override; virtual boost::any visit(EventuallyFormula const& f, boost::any const& data) const override; virtual boost::any visit(TimeOperatorFormula const& f, boost::any const& data) const override; virtual boost::any visit(GloballyFormula const& f, boost::any const& data) const override; diff --git a/src/storm/logic/FormulaVisitor.h b/src/storm/logic/FormulaVisitor.h index cb8f71d0d4..c4eb4cb78e 100644 --- a/src/storm/logic/FormulaVisitor.h +++ b/src/storm/logic/FormulaVisitor.h @@ -20,6 +20,7 @@ class FormulaVisitor { virtual boost::any visit(BoundedUntilFormula const& f, boost::any const& data) const = 0; virtual boost::any visit(ConditionalFormula const& f, boost::any const& data) const = 0; virtual boost::any visit(CumulativeRewardFormula const& f, boost::any const& data) const = 0; + virtual boost::any visit(CvarFormula const& f, boost::any const& data) const = 0; virtual boost::any visit(EventuallyFormula const& f, boost::any const& data) const = 0; virtual boost::any visit(TimeOperatorFormula const& f, boost::any const& data) const = 0; virtual boost::any visit(GloballyFormula const& f, boost::any const& data) const = 0; diff --git a/src/storm/logic/Formulas.h b/src/storm/logic/Formulas.h index e2937bc7ce..eddd5fcaf5 100644 --- a/src/storm/logic/Formulas.h +++ b/src/storm/logic/Formulas.h @@ -11,6 +11,7 @@ #include "storm/logic/ComparisonType.h" #include "storm/logic/ConditionalFormula.h" #include "storm/logic/CumulativeRewardFormula.h" +#include "storm/logic/CvarFormula.h" #include "storm/logic/DiscountedCumulativeRewardFormula.h" #include "storm/logic/DiscountedTotalRewardFormula.h" #include "storm/logic/EventuallyFormula.h" diff --git a/src/storm/logic/FormulasForwardDeclarations.h b/src/storm/logic/FormulasForwardDeclarations.h index 4f1cec67dd..af4bc82235 100644 --- a/src/storm/logic/FormulasForwardDeclarations.h +++ b/src/storm/logic/FormulasForwardDeclarations.h @@ -14,6 +14,7 @@ class BooleanLiteralFormula; class BoundedUntilFormula; class ConditionalFormula; class CumulativeRewardFormula; +class CvarFormula; class EventuallyFormula; class TimeOperatorFormula; class GloballyFormula; diff --git a/src/storm/logic/FragmentChecker.cpp b/src/storm/logic/FragmentChecker.cpp index 8dc07f5f2d..a151cecea8 100644 --- a/src/storm/logic/FragmentChecker.cpp +++ b/src/storm/logic/FragmentChecker.cpp @@ -31,6 +31,9 @@ bool FragmentChecker::conformsToSpecification(Formula const& f, FragmentSpecific if (specification.isQuantileFormulaAtTopLevelRequired()) { result &= f.isQuantileFormula(); } + if (specification.isCvarFormulaAtTopLevelRequired()) { + result &= f.isCvarFormula(); + } return result; } @@ -151,6 +154,14 @@ boost::any FragmentChecker::visit(CumulativeRewardFormula const& f, boost::any c return result; } +boost::any FragmentChecker::visit(CvarFormula const& f, boost::any const& data) const { + InheritedInformation const& inherited = boost::any_cast(data); + if (!inherited.getSpecification().areCvarFormulasAllowed()) { + return false; + } + return f.getSubformula().accept(*this, data); +} + boost::any FragmentChecker::visit(EventuallyFormula const& f, boost::any const& data) const { InheritedInformation const& inherited = boost::any_cast(data); bool result = true; diff --git a/src/storm/logic/FragmentChecker.h b/src/storm/logic/FragmentChecker.h index eabf04cdb7..aadda36863 100644 --- a/src/storm/logic/FragmentChecker.h +++ b/src/storm/logic/FragmentChecker.h @@ -19,6 +19,7 @@ class FragmentChecker : public FormulaVisitor { virtual boost::any visit(BoundedUntilFormula const& f, boost::any const& data) const override; virtual boost::any visit(ConditionalFormula const& f, boost::any const& data) const override; virtual boost::any visit(CumulativeRewardFormula const& f, boost::any const& data) const override; + virtual boost::any visit(CvarFormula const& f, boost::any const& data) const override; virtual boost::any visit(EventuallyFormula const& f, boost::any const& data) const override; virtual boost::any visit(TimeOperatorFormula const& f, boost::any const& data) const override; virtual boost::any visit(GloballyFormula const& f, boost::any const& data) const override; diff --git a/src/storm/logic/FragmentSpecification.cpp b/src/storm/logic/FragmentSpecification.cpp index 63a546f6b9..3076518b61 100644 --- a/src/storm/logic/FragmentSpecification.cpp +++ b/src/storm/logic/FragmentSpecification.cpp @@ -239,6 +239,17 @@ FragmentSpecification quantiles() { return quantiles; } +FragmentSpecification cvars() { + FragmentSpecification cvars = propositional(); + + cvars.setCvarFormulasAllowed(true); + cvars.setCvarFormulaAtTopLevelRequired(true); + cvars.setRewardOperatorsAllowed(true); + cvars.setReachabilityRewardFormulasAllowed(true); + + return cvars; +} + FragmentSpecification::FragmentSpecification() { probabilityOperator = false; rewardOperator = false; @@ -247,6 +258,7 @@ FragmentSpecification::FragmentSpecification() { multiObjectiveFormula = false; quantileFormula = false; + cvarFormula = false; globallyFormula = false; reachabilityProbabilityFormula = false; @@ -295,8 +307,9 @@ FragmentSpecification::FragmentSpecification() { operatorAtTopLevelRequired = false; multiObjectiveFormulaAtTopLevelRequired = false; - operatorsAtTopLevelOfMultiObjectiveFormulasRequired = false; quantileFormulaAtTopLevelRequired = false; + cvarFormulaAtTopLevelRequired = false; + operatorsAtTopLevelOfMultiObjectiveFormulasRequired = false; rewardAccumulation = false; @@ -362,6 +375,15 @@ FragmentSpecification& FragmentSpecification::setQuantileFormulasAllowed(bool ne return *this; } +bool FragmentSpecification::areCvarFormulasAllowed() const { + return cvarFormula; +} + +FragmentSpecification& FragmentSpecification::setCvarFormulasAllowed(bool newValue) { + this->cvarFormula = newValue; + return *this; +} + bool FragmentSpecification::areGloballyFormulasAllowed() const { return globallyFormula; } @@ -741,6 +763,15 @@ FragmentSpecification& FragmentSpecification::setQuantileFormulaAtTopLevelRequir return *this; } +bool FragmentSpecification::isCvarFormulaAtTopLevelRequired() const { + return cvarFormulaAtTopLevelRequired; +} + +FragmentSpecification& FragmentSpecification::setCvarFormulaAtTopLevelRequired(bool newValue) { + cvarFormulaAtTopLevelRequired = newValue; + return *this; +} + bool FragmentSpecification::isRewardAccumulationAllowed() const { return rewardAccumulation; } diff --git a/src/storm/logic/FragmentSpecification.h b/src/storm/logic/FragmentSpecification.h index 9ec6d1708e..827ef9dfcb 100644 --- a/src/storm/logic/FragmentSpecification.h +++ b/src/storm/logic/FragmentSpecification.h @@ -36,6 +36,9 @@ class FragmentSpecification { bool areQuantileFormulasAllowed() const; FragmentSpecification& setQuantileFormulasAllowed(bool newValue); + bool areCvarFormulasAllowed() const; + FragmentSpecification& setCvarFormulasAllowed(bool newValue); + bool areGloballyFormulasAllowed() const; FragmentSpecification& setGloballyFormulasAllowed(bool newValue); @@ -156,6 +159,9 @@ class FragmentSpecification { bool isQuantileFormulaAtTopLevelRequired() const; FragmentSpecification& setQuantileFormulaAtTopLevelRequired(bool newValue); + bool isCvarFormulaAtTopLevelRequired() const; + FragmentSpecification& setCvarFormulaAtTopLevelRequired(bool newValue); + bool isRewardAccumulationAllowed() const; FragmentSpecification& setRewardAccumulationAllowed(bool newValue); @@ -181,6 +187,7 @@ class FragmentSpecification { bool multiObjectiveFormula; bool quantileFormula; + bool cvarFormula; bool globallyFormula; bool reachabilityProbabilityFormula; @@ -229,6 +236,7 @@ class FragmentSpecification { bool operatorAtTopLevelRequired; bool multiObjectiveFormulaAtTopLevelRequired; bool quantileFormulaAtTopLevelRequired; + bool cvarFormulaAtTopLevelRequired; bool operatorsAtTopLevelOfMultiObjectiveFormulasRequired; bool rewardAccumulation; @@ -281,5 +289,8 @@ FragmentSpecification lexObjective(); // Quantile formulas. FragmentSpecification quantiles(); +// CVaR formulas. +FragmentSpecification cvars(); + } // namespace logic } // namespace storm diff --git a/src/storm/logic/LiftableTransitionRewardsVisitor.cpp b/src/storm/logic/LiftableTransitionRewardsVisitor.cpp index 9bd6ec5d37..e35590b9d6 100644 --- a/src/storm/logic/LiftableTransitionRewardsVisitor.cpp +++ b/src/storm/logic/LiftableTransitionRewardsVisitor.cpp @@ -70,6 +70,10 @@ boost::any LiftableTransitionRewardsVisitor::visit(CumulativeRewardFormula const return true; } +boost::any LiftableTransitionRewardsVisitor::visit(CvarFormula const& f, boost::any const& data) const { + return f.getSubformula().accept(*this, data); +} + boost::any LiftableTransitionRewardsVisitor::visit(EventuallyFormula const& f, boost::any const& data) const { return f.getSubformula().accept(*this, data); } diff --git a/src/storm/logic/LiftableTransitionRewardsVisitor.h b/src/storm/logic/LiftableTransitionRewardsVisitor.h index 27cf39bb20..cd300aab0d 100644 --- a/src/storm/logic/LiftableTransitionRewardsVisitor.h +++ b/src/storm/logic/LiftableTransitionRewardsVisitor.h @@ -25,6 +25,7 @@ class LiftableTransitionRewardsVisitor : public FormulaVisitor { virtual boost::any visit(BoundedUntilFormula const& f, boost::any const& data) const override; virtual boost::any visit(ConditionalFormula const& f, boost::any const& data) const override; virtual boost::any visit(CumulativeRewardFormula const& f, boost::any const& data) const override; + virtual boost::any visit(CvarFormula const& f, boost::any const& data) const override; virtual boost::any visit(EventuallyFormula const& f, boost::any const& data) const override; virtual boost::any visit(TimeOperatorFormula const& f, boost::any const& data) const override; virtual boost::any visit(GloballyFormula const& f, boost::any const& data) const override; diff --git a/src/storm/logic/ToExpressionVisitor.cpp b/src/storm/logic/ToExpressionVisitor.cpp index 88da97df38..bb96e5fb21 100644 --- a/src/storm/logic/ToExpressionVisitor.cpp +++ b/src/storm/logic/ToExpressionVisitor.cpp @@ -63,6 +63,10 @@ boost::any ToExpressionVisitor::visit(CumulativeRewardFormula const&, boost::any STORM_LOG_THROW(false, storm::exceptions::InvalidOperationException, "Cannot assemble expression from formula that contains illegal elements."); } +boost::any ToExpressionVisitor::visit(CvarFormula const&, boost::any const&) const { + STORM_LOG_THROW(false, storm::exceptions::InvalidOperationException, "Cannot assemble expression from formula that contains illegal elements."); +} + boost::any ToExpressionVisitor::visit(EventuallyFormula const&, boost::any const&) const { STORM_LOG_THROW(false, storm::exceptions::InvalidOperationException, "Cannot assemble expression from formula that contains illegal elements."); } diff --git a/src/storm/logic/ToExpressionVisitor.h b/src/storm/logic/ToExpressionVisitor.h index 4ce6e1f70e..47c093fdc3 100644 --- a/src/storm/logic/ToExpressionVisitor.h +++ b/src/storm/logic/ToExpressionVisitor.h @@ -19,6 +19,7 @@ class ToExpressionVisitor : public FormulaVisitor { virtual boost::any visit(BoundedUntilFormula const& f, boost::any const& data) const override; virtual boost::any visit(ConditionalFormula const& f, boost::any const& data) const override; virtual boost::any visit(CumulativeRewardFormula const& f, boost::any const& data) const override; + virtual boost::any visit(CvarFormula const& f, boost::any const& data) const override; virtual boost::any visit(EventuallyFormula const& f, boost::any const& data) const override; virtual boost::any visit(TimeOperatorFormula const& f, boost::any const& data) const override; virtual boost::any visit(GloballyFormula const& f, boost::any const& data) const override; diff --git a/src/storm/logic/ToPrefixStringVisitor.cpp b/src/storm/logic/ToPrefixStringVisitor.cpp index ef393499bf..562eda13f1 100644 --- a/src/storm/logic/ToPrefixStringVisitor.cpp +++ b/src/storm/logic/ToPrefixStringVisitor.cpp @@ -104,6 +104,10 @@ boost::any ToPrefixStringVisitor::visit(CumulativeRewardFormula const&, boost::a STORM_LOG_THROW(false, storm::exceptions::InvalidOperationException, "Can not convert to prefix string."); } +boost::any ToPrefixStringVisitor::visit(CvarFormula const&, boost::any const&) const { + STORM_LOG_THROW(false, storm::exceptions::InvalidOperationException, "Can not convert to prefix string"); +} + boost::any ToPrefixStringVisitor::visit(EventuallyFormula const& f, boost::any const& data) const { std::string subexpression = boost::any_cast(f.getSubformula().accept(*this, data)); return std::string("F ") + subexpression; diff --git a/src/storm/logic/ToPrefixStringVisitor.h b/src/storm/logic/ToPrefixStringVisitor.h index 2fe860e01c..0f3a1bd498 100644 --- a/src/storm/logic/ToPrefixStringVisitor.h +++ b/src/storm/logic/ToPrefixStringVisitor.h @@ -19,6 +19,7 @@ class ToPrefixStringVisitor : public FormulaVisitor { virtual boost::any visit(BoundedUntilFormula const& f, boost::any const& data) const override; virtual boost::any visit(ConditionalFormula const& f, boost::any const& data) const override; virtual boost::any visit(CumulativeRewardFormula const& f, boost::any const& data) const override; + virtual boost::any visit(CvarFormula const& f, boost::any const& data) const override; virtual boost::any visit(EventuallyFormula const& f, boost::any const& data) const override; virtual boost::any visit(TimeOperatorFormula const& f, boost::any const& data) const override; virtual boost::any visit(GloballyFormula const& f, boost::any const& data) const override; diff --git a/src/storm/modelchecker/AbstractModelChecker.cpp b/src/storm/modelchecker/AbstractModelChecker.cpp index 165fe7ea3f..ccf47e33ad 100644 --- a/src/storm/modelchecker/AbstractModelChecker.cpp +++ b/src/storm/modelchecker/AbstractModelChecker.cpp @@ -296,6 +296,8 @@ std::unique_ptr AbstractModelChecker::checkStateFormula( STORM_LOG_ASSERT(mof.isTradeoff(), "Unexpected multi-objective formula type."); return this->checkMultiObjectiveFormula(env, checkTask.substituteFormula(mof)); } + } else if (stateFormula.isCvarFormula()) { + return this->checkCvarFormula(env, checkTask.substituteFormula(stateFormula.asCvarFormula())); } else if (stateFormula.isQuantileFormula()) { return this->checkQuantileFormula(env, checkTask.substituteFormula(stateFormula.asQuantileFormula())); } @@ -444,6 +446,13 @@ std::unique_ptr AbstractModelChecker::checkMultiObjectiv "This model checker (" << getClassName() << ") does not support the formula: " << checkTask.getFormula() << "."); } +template +std::unique_ptr AbstractModelChecker::checkCvarFormula(Environment const&, + CheckTask const& checkTask) { + STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, + "This model checker (" << getClassName() << ") does not support the formula: " << checkTask.getFormula() << "."); +} + template std::unique_ptr AbstractModelChecker::checkQuantileFormula(Environment const&, CheckTask const& checkTask) { diff --git a/src/storm/modelchecker/AbstractModelChecker.h b/src/storm/modelchecker/AbstractModelChecker.h index 51261bddce..e2e9ce1f3c 100644 --- a/src/storm/modelchecker/AbstractModelChecker.h +++ b/src/storm/modelchecker/AbstractModelChecker.h @@ -138,6 +138,9 @@ class AbstractModelChecker { virtual std::unique_ptr checkMultiObjectiveFormula(Environment const& env, CheckTask const& checkTask); + // The methods to check CVaR formulas. + virtual std::unique_ptr checkCvarFormula(Environment const& env, CheckTask const& checkTask); + // The methods to check quantile formulas. virtual std::unique_ptr checkQuantileFormula(Environment const& env, CheckTask const& checkTask); diff --git a/src/storm/modelchecker/CheckTask.h b/src/storm/modelchecker/CheckTask.h index ae1960341c..c7714b6b77 100644 --- a/src/storm/modelchecker/CheckTask.h +++ b/src/storm/modelchecker/CheckTask.h @@ -63,12 +63,13 @@ class CheckTask { * Calling this method has no effect if the provided formula is not an operator formula. */ void updateOperatorInformation() { - if (formula.get().isOperatorFormula()) { - storm::logic::OperatorFormula const& operatorFormula = formula.get().asOperatorFormula(); - if (operatorFormula.hasOptimalityType()) { - this->optimizationDirection = operatorFormula.getOptimalityType(); - } + storm::logic::Formula const* formulaForOperatorInformation = &formula.get(); + if (formulaForOperatorInformation->isCvarFormula()) { + formulaForOperatorInformation = &formulaForOperatorInformation->asCvarFormula().getSubformula(); + } + if (formulaForOperatorInformation->isOperatorFormula()) { + storm::logic::OperatorFormula const& operatorFormula = formulaForOperatorInformation->asOperatorFormula(); if (operatorFormula.hasBound()) { this->bound = operatorFormula.getBound(); } @@ -82,8 +83,8 @@ class CheckTask { : OptimizationDirection::Minimize; } - if (formula.get().isProbabilityOperatorFormula()) { - storm::logic::ProbabilityOperatorFormula const& probabilityOperatorFormula = formula.get().asProbabilityOperatorFormula(); + if (formulaForOperatorInformation->isProbabilityOperatorFormula()) { + storm::logic::ProbabilityOperatorFormula const& probabilityOperatorFormula = formulaForOperatorInformation->asProbabilityOperatorFormula(); if (probabilityOperatorFormula.hasBound()) { if (storm::utility::isZero(probabilityOperatorFormula.template getThresholdAs()) || @@ -91,8 +92,8 @@ class CheckTask { this->qualitative = true; } } - } else if (formula.get().isRewardOperatorFormula()) { - storm::logic::RewardOperatorFormula const& rewardOperatorFormula = formula.get().asRewardOperatorFormula(); + } else if (formulaForOperatorInformation->isRewardOperatorFormula()) { + storm::logic::RewardOperatorFormula const& rewardOperatorFormula = formulaForOperatorInformation->asRewardOperatorFormula(); this->rewardModel = rewardOperatorFormula.getOptionalRewardModelName(); if (rewardOperatorFormula.hasBound()) { diff --git a/src/storm/modelchecker/cvar/CvarClassification.h b/src/storm/modelchecker/cvar/CvarClassification.h new file mode 100644 index 0000000000..8de9d1203c --- /dev/null +++ b/src/storm/modelchecker/cvar/CvarClassification.h @@ -0,0 +1,60 @@ +#pragma once + +#include "storm/exceptions/InvalidPropertyException.h" +#include "storm/exceptions/NotImplementedException.h" +#include "storm/modelchecker/cvar/CvarMethod.h" +#include "storm/modelchecker/cvar/CvarQueryInformation.h" +#include "storm/storage/BitVector.h" +#include "storm/utility/macros.h" + +namespace storm { +namespace modelchecker { +namespace cvar { +/*! + * Selects the concrete CVaR backend induced by a query on a given + * model and reward structure. + * + * Weighted reachability is the currently implemented LP-based terminal-reward + * setting. SSP uses Pareto value iteration for accumulated costs until + * reaching the goal. + */ +enum class CvarBackendKind { WeightedReachability, Ssp }; + +/*! + * Selects the concrete CVaR backend to use. + * + * The selection can be overridden explicitly via the CVaR method setting. + * Otherwise, classification stays conservative: state-action reward models are + * routed to the SSP branch, while state-only reward models remain on the + * weighted-reachability path until SSP preprocessing is introduced. + */ +template +CvarBackendKind selectCvarBackend(SparseMdpModelType const& model, CvarQueryInformation const& queryInformation, storm::storage::BitVector const&, + CvarMethod method) { + std::string rewardModelName = queryInformation.rewardModelName.value_or(""); + auto const& rewardModel = model.getRewardModel(rewardModelName); + if (rewardModelName.empty()) { + rewardModelName = model.getUniqueRewardModelName(); + } + + if (method == CvarMethod::WeightedReachability) { + STORM_LOG_THROW(!rewardModel.hasStateActionRewards() && !rewardModel.hasTransitionRewards(), storm::exceptions::InvalidPropertyException, + "The weighted-reachability CVaR method requires state-based terminal rewards only."); + return CvarBackendKind::WeightedReachability; + } + + STORM_LOG_THROW(!rewardModel.hasTransitionRewards(), storm::exceptions::NotImplementedException, + "CVaR queries with transition rewards are not supported yet."); + + if (method == CvarMethod::SspParetoVi) { + return CvarBackendKind::Ssp; + } + + if (rewardModel.hasStateActionRewards()) { + return CvarBackendKind::Ssp; + } + return CvarBackendKind::WeightedReachability; +} +} // namespace cvar +} // namespace modelchecker +} // namespace storm diff --git a/src/storm/modelchecker/cvar/CvarComputationResult.h b/src/storm/modelchecker/cvar/CvarComputationResult.h new file mode 100644 index 0000000000..706e8d89dc --- /dev/null +++ b/src/storm/modelchecker/cvar/CvarComputationResult.h @@ -0,0 +1,19 @@ +#pragma once + +#include + +#include "storm/storage/Scheduler.h" + +namespace storm { +namespace modelchecker { +namespace cvar { + +template +struct CvarComputationResult { + ValueType value; + std::unique_ptr> scheduler; +}; + +} // namespace cvar +} // namespace modelchecker +} // namespace storm diff --git a/src/storm/modelchecker/cvar/CvarInterpretation.h b/src/storm/modelchecker/cvar/CvarInterpretation.h new file mode 100644 index 0000000000..eaed7778ee --- /dev/null +++ b/src/storm/modelchecker/cvar/CvarInterpretation.h @@ -0,0 +1,29 @@ +#pragma once + +#include "storm/solver/OptimizationDirection.h" +#include "storm/utility/macros.h" + +namespace storm { +namespace modelchecker { +namespace cvar { + +enum class CvarInterpretation { Cost, Reward }; + +enum class CvarInterpretationSelection { Auto, Cost, Reward }; + +inline CvarInterpretation resolveCvarInterpretation(CvarInterpretationSelection selection, storm::solver::OptimizationDirection optimizationDirection) { + switch (selection) { + case CvarInterpretationSelection::Auto: + return storm::solver::minimize(optimizationDirection) ? CvarInterpretation::Cost : CvarInterpretation::Reward; + case CvarInterpretationSelection::Cost: + return CvarInterpretation::Cost; + case CvarInterpretationSelection::Reward: + return CvarInterpretation::Reward; + } + STORM_LOG_ASSERT(false, "Encountered an unknown CVaR interpretation selection."); + return CvarInterpretation::Cost; +} + +} // namespace cvar +} // namespace modelchecker +} // namespace storm diff --git a/src/storm/modelchecker/cvar/CvarMethod.cpp b/src/storm/modelchecker/cvar/CvarMethod.cpp new file mode 100644 index 0000000000..7fd0f06a1d --- /dev/null +++ b/src/storm/modelchecker/cvar/CvarMethod.cpp @@ -0,0 +1,21 @@ +#include "storm/modelchecker/cvar/CvarMethod.h" + +namespace storm { +namespace modelchecker { +namespace cvar { + +std::string toString(CvarMethod method) { + switch (method) { + case CvarMethod::Auto: + return "auto"; + case CvarMethod::WeightedReachability: + return "weighted-reachability"; + case CvarMethod::SspParetoVi: + return "ssp-vi"; + } + return "unknown"; +} + +} // namespace cvar +} // namespace modelchecker +} // namespace storm diff --git a/src/storm/modelchecker/cvar/CvarMethod.h b/src/storm/modelchecker/cvar/CvarMethod.h new file mode 100644 index 0000000000..dc40e7a46b --- /dev/null +++ b/src/storm/modelchecker/cvar/CvarMethod.h @@ -0,0 +1,13 @@ +#pragma once + +#include "storm/utility/ExtendSettingEnumWithSelectionField.h" + +namespace storm { +namespace modelchecker { +namespace cvar { + +ExtendEnumsWithSelectionField(CvarMethod, Auto, WeightedReachability, SspParetoVi) std::string toString(CvarMethod method); + +} // namespace cvar +} // namespace modelchecker +} // namespace storm diff --git a/src/storm/modelchecker/cvar/CvarModelChecking.cpp b/src/storm/modelchecker/cvar/CvarModelChecking.cpp new file mode 100644 index 0000000000..641dea840d --- /dev/null +++ b/src/storm/modelchecker/cvar/CvarModelChecking.cpp @@ -0,0 +1,56 @@ +#include "storm/modelchecker/cvar/CvarModelChecking.h" + +#include "storm/adapters/RationalNumberAdapter.h" +#include "storm/environment/Environment.h" +#include "storm/exceptions/InvalidOperationException.h" +#include "storm/exceptions/InvalidPropertyException.h" +#include "storm/modelchecker/cvar/CvarInterpretation.h" +#include "storm/modelchecker/cvar/CvarQueryInformation.h" +#include "storm/modelchecker/cvar/helper/SparseCvarComputationHelper.h" +#include "storm/modelchecker/results/ExplicitQuantitativeCheckResult.h" +#include "storm/models/sparse/Mdp.h" +#include "storm/utility/macros.h" + +namespace storm { +namespace modelchecker { +namespace cvar { + +template +std::unique_ptr performCvarModelChecking(Environment const& env, SparseMdpModelType const& model, + CheckTask const& checkTask, + std::function const& formulaChecker) { + using ValueType = typename SparseMdpModelType::ValueType; + + STORM_LOG_THROW(checkTask.isOnlyInitialStatesRelevantSet(), storm::exceptions::InvalidOperationException, + "Computing CVaR is only supported for the initial states of a model."); + STORM_LOG_THROW(model.getInitialStates().getNumberOfSetBits() == 1, storm::exceptions::InvalidOperationException, + "CVaR is not supported on models with multiple initial states."); + STORM_LOG_THROW(checkTask.isOptimizationDirectionSet(), storm::exceptions::InvalidPropertyException, + "The embedded reward operator formula of a CVaR query must specify whether to minimize or maximize."); + + auto interpretation = resolveCvarInterpretation(env.modelchecker().cvar().getInterpretationSelection(), checkTask.getOptimizationDirection()); + auto cvarQueryInformation = extractCvarQueryInformation(checkTask.getFormula(), interpretation); + auto targetStates = formulaChecker(*cvarQueryInformation.targetFormula); + + SparseCvarComputationHelper cvarHelper(model, cvarQueryInformation, targetStates); + auto cvarResult = cvarHelper.computeCvar(env, checkTask.isProduceSchedulersSet()); + + std::unique_ptr result(new ExplicitQuantitativeCheckResult(*model.getInitialStates().begin(), std::move(cvarResult.value))); + if (checkTask.isProduceSchedulersSet() && cvarResult.scheduler) { + result->asExplicitQuantitativeCheckResult().setScheduler(std::move(cvarResult.scheduler)); + } + return result; +} + +template std::unique_ptr performCvarModelChecking>( + Environment const& env, storm::models::sparse::Mdp const& model, CheckTask const& checkTask, + std::function const& formulaChecker); + +template std::unique_ptr performCvarModelChecking>( + Environment const& env, storm::models::sparse::Mdp const& model, + CheckTask const& checkTask, + std::function const& formulaChecker); + +} // namespace cvar +} // namespace modelchecker +} // namespace storm diff --git a/src/storm/modelchecker/cvar/CvarModelChecking.h b/src/storm/modelchecker/cvar/CvarModelChecking.h new file mode 100644 index 0000000000..33795430b5 --- /dev/null +++ b/src/storm/modelchecker/cvar/CvarModelChecking.h @@ -0,0 +1,29 @@ +#pragma once + +#include +#include + +#include "storm/logic/CvarFormula.h" +#include "storm/modelchecker/CheckTask.h" +#include "storm/modelchecker/results/CheckResult.h" +#include "storm/storage/BitVector.h" + +namespace storm { + +class Environment; + +namespace logic { +class Formula; +} + +namespace modelchecker { +namespace cvar { + +template +std::unique_ptr performCvarModelChecking(Environment const& env, SparseMdpModelType const& model, + CheckTask const& checkTask, + std::function const& formulaChecker); + +} // namespace cvar +} // namespace modelchecker +} // namespace storm diff --git a/src/storm/modelchecker/cvar/CvarPreprocessingUtilities.h b/src/storm/modelchecker/cvar/CvarPreprocessingUtilities.h new file mode 100644 index 0000000000..9f66df31ad --- /dev/null +++ b/src/storm/modelchecker/cvar/CvarPreprocessingUtilities.h @@ -0,0 +1,61 @@ +#pragma once + +#include "storm/exceptions/InvalidPropertyException.h" +#include "storm/storage/BitVector.h" +#include "storm/storage/MaximalEndComponentDecomposition.h" +#include "storm/storage/SparseMatrix.h" +#include "storm/utility/graph.h" +#include "storm/utility/macros.h" + +namespace storm { +namespace modelchecker { +namespace cvar { + +template +void validateTargetStatesAreAbsorbing(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& targetStates) { + for (auto targetState : targetStates) { + for (uint64_t row = transitionMatrix.getRowGroupIndices()[targetState], endRow = transitionMatrix.getRowGroupIndices()[targetState + 1]; row < endRow; + ++row) { + for (auto const& entry : transitionMatrix.getRow(row)) { + STORM_LOG_THROW(entry.getColumn() == targetState, storm::exceptions::InvalidPropertyException, + "CVaR query currently requires all original target states to be absorbing."); + } + } + } +} + +template +storm::storage::MaximalEndComponentDecomposition computeReachableMecs(storm::storage::SparseMatrix const& transitionMatrix, + storm::storage::BitVector const& initialStates) { + storm::storage::BitVector allStates(transitionMatrix.getRowGroupCount(), true); + storm::storage::BitVector noStates(transitionMatrix.getRowGroupCount(), false); + auto reachableStates = storm::utility::graph::getReachableStates(transitionMatrix, initialStates, allStates, noStates); + auto backwardTransitions = transitionMatrix.transpose(true); + return storm::storage::MaximalEndComponentDecomposition(transitionMatrix, backwardTransitions, reachableStates); +} + +template +storm::storage::BitVector computeStatesThatCanReachTarget(storm::storage::SparseMatrix const& transitionMatrix, + storm::storage::BitVector const& targetStates) { + storm::storage::BitVector allStates(transitionMatrix.getRowGroupCount(), true); + auto backwardTransitions = transitionMatrix.transpose(true); + return storm::utility::graph::performProbGreater0(backwardTransitions, allStates, targetStates); +} + +template +storm::storage::BitVector computeBadMecStates(storm::storage::MaximalEndComponentDecomposition const& mecs, + storm::storage::BitVector const& statesThatCanReachTarget, uint64_t numberOfStates) { + storm::storage::BitVector badMecStates(numberOfStates, false); + for (auto const& mec : mecs) { + if (!mec.containsAnyState(statesThatCanReachTarget)) { + for (auto const& stateChoices : mec) { + badMecStates.set(stateChoices.first, true); + } + } + } + return badMecStates; +} + +} // namespace cvar +} // namespace modelchecker +} // namespace storm diff --git a/src/storm/modelchecker/cvar/CvarQueryInformation.cpp b/src/storm/modelchecker/cvar/CvarQueryInformation.cpp new file mode 100644 index 0000000000..63c47fed4b --- /dev/null +++ b/src/storm/modelchecker/cvar/CvarQueryInformation.cpp @@ -0,0 +1,42 @@ +#include "storm/modelchecker/cvar/CvarQueryInformation.h" + +#include "storm/exceptions/InvalidPropertyException.h" +#include "storm/logic/EventuallyFormula.h" +#include "storm/logic/RewardOperatorFormula.h" +#include "storm/utility/logging.h" +#include "storm/utility/macros.h" + +namespace storm { +namespace modelchecker { +namespace cvar { + +CvarQueryInformation extractCvarQueryInformation(storm::logic::CvarFormula const& formula, CvarInterpretation interpretation) { + storm::logic::Formula const& embeddedFormula = formula.getSubformula(); + STORM_LOG_THROW(embeddedFormula.isRewardOperatorFormula(), storm::exceptions::InvalidPropertyException, + "CVaR formulas currently require an embedded reward operator formula."); + + auto const& rewardOperator = embeddedFormula.asRewardOperatorFormula(); + STORM_LOG_THROW(rewardOperator.hasOptimalityType(), storm::exceptions::InvalidPropertyException, + "The embedded reward operator formula of a CVaR query must specify whether to minimize or maximize."); + STORM_LOG_THROW(!rewardOperator.hasBound(), storm::exceptions::InvalidPropertyException, + "The embedded reward operator formula of a CVaR query must not specify a threshold."); + + storm::logic::Formula const& rewardPathFormula = rewardOperator.getSubformula(); + STORM_LOG_THROW(rewardPathFormula.isReachabilityRewardFormula(), storm::exceptions::InvalidPropertyException, + "CVaR queries currently only support weighted reachability objectives."); + STORM_LOG_THROW(rewardPathFormula.isEventuallyFormula(), storm::exceptions::InvalidPropertyException, + "CVaR queries currently only support reachability reward formulas of the form Rmin/max=? [ F target ]."); + + auto const& eventuallyFormula = rewardPathFormula.asEventuallyFormula(); + STORM_LOG_THROW(eventuallyFormula.getSubformula().isStateFormula(), storm::exceptions::InvalidPropertyException, + "The target of the embedded reachability reward formula of a CVaR query must be a state formula."); + + auto const optimizationDirection = rewardOperator.getOptimalityType(); + auto const& optionalRewardModelName = rewardOperator.getOptionalRewardModelName(); + return {formula.getAlpha(), optimizationDirection, interpretation, + optionalRewardModelName ? std::optional(optionalRewardModelName.get()) : std::nullopt, + eventuallyFormula.getSubformula().asSharedPointer()}; +} +} // namespace cvar +} // namespace modelchecker +} // namespace storm diff --git a/src/storm/modelchecker/cvar/CvarQueryInformation.h b/src/storm/modelchecker/cvar/CvarQueryInformation.h new file mode 100644 index 0000000000..2b92a99865 --- /dev/null +++ b/src/storm/modelchecker/cvar/CvarQueryInformation.h @@ -0,0 +1,28 @@ +#pragma once + +#include +#include +#include + +#include "storm/adapters/RationalNumberAdapter.h" +#include "storm/logic/CvarFormula.h" +#include "storm/modelchecker/cvar/CvarInterpretation.h" +#include "storm/solver/OptimizationDirection.h" + +namespace storm { +namespace modelchecker { +namespace cvar { + +struct CvarQueryInformation { + storm::RationalNumber alpha; + storm::solver::OptimizationDirection optimizationDirection; + CvarInterpretation interpretation; + std::optional rewardModelName; + std::shared_ptr targetFormula; +}; + +CvarQueryInformation extractCvarQueryInformation(storm::logic::CvarFormula const& formula, CvarInterpretation interpretation); + +} // namespace cvar +} // namespace modelchecker +} // namespace storm diff --git a/src/storm/modelchecker/cvar/helper/SparseCvarComputationHelper.h b/src/storm/modelchecker/cvar/helper/SparseCvarComputationHelper.h new file mode 100644 index 0000000000..999c2b0cb2 --- /dev/null +++ b/src/storm/modelchecker/cvar/helper/SparseCvarComputationHelper.h @@ -0,0 +1,70 @@ +#pragma once + +#include "storm/environment/Environment.h" +#include "storm/environment/modelchecker/ModelCheckerEnvironment.h" +#include "storm/exceptions/InvalidPropertyException.h" +#include "storm/exceptions/NotImplementedException.h" +#include "storm/exceptions/UnexpectedException.h" +#include "storm/modelchecker/cvar/CvarClassification.h" +#include "storm/modelchecker/cvar/CvarComputationResult.h" +#include "storm/modelchecker/cvar/CvarQueryInformation.h" +#include "storm/modelchecker/cvar/helper/SparseSspCvarParetoViHelper.h" +#include "storm/modelchecker/cvar/helper/SparseSspRewardCvarParetoViHelper.h" +#include "storm/modelchecker/cvar/helper/SparseWeightedReachabilityCvarLpHelper.h" +#include "storm/modelchecker/cvar/preprocessing/SspCvarPreprocessor.h" +#include "storm/modelchecker/cvar/preprocessing/WeightedReachabilityCvarPreprocessor.h" +#include "storm/storage/BitVector.h" + +namespace storm { +namespace modelchecker { +namespace cvar { + +template +class SparseCvarComputationHelper { + public: + using ValueType = typename SparseMdpModelType::ValueType; + + SparseCvarComputationHelper(SparseMdpModelType const& model, CvarQueryInformation const& queryInformation, storm::storage::BitVector const& targetStates) + : model(model), queryInformation(queryInformation), targetStates(targetStates) { + // Intentionally left empty. + } + + CvarComputationResult computeCvar(Environment const& env, bool produceScheduler = false) const { + auto backendKind = selectCvarBackend(model, queryInformation, targetStates, env.modelchecker().cvar().getMethod()); + + switch (backendKind) { + case CvarBackendKind::WeightedReachability: { + auto weightedReachabilityPreprocessingResult = + preprocessing::preprocessWeightedReachabilityCvar(model, queryInformation, targetStates, produceScheduler); + SparseWeightedReachabilityCvarLpHelper cvarHelper(queryInformation, weightedReachabilityPreprocessingResult); + return cvarHelper.computeCvar(env, produceScheduler); + } + case CvarBackendKind::Ssp: { + if (queryInformation.optimizationDirection == storm::solver::OptimizationDirection::Minimize && + queryInformation.interpretation == CvarInterpretation::Cost) { + auto sspPreprocessingResult = preprocessing::preprocessSspCvar(env, model, queryInformation, targetStates); + SparseSspCvarParetoViHelper cvarHelper(queryInformation, sspPreprocessingResult); + return cvarHelper.computeCvar(env, produceScheduler); + } + if (queryInformation.optimizationDirection == storm::solver::OptimizationDirection::Maximize && + queryInformation.interpretation == CvarInterpretation::Reward) { + auto sspPreprocessingResult = preprocessing::preprocessSspRewardCvar(env, model, queryInformation, targetStates); + SparseSspRewardCvarParetoViHelper cvarHelper(queryInformation, sspPreprocessingResult); + return cvarHelper.computeCvar(env, produceScheduler); + } + STORM_LOG_THROW(false, storm::exceptions::InvalidPropertyException, + "CVaR SSP value iteration currently supports only minimizing costs and maximizing rewards."); + } + } + STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "Encountered an unknown CVaR backend."); + } + + private: + SparseMdpModelType const& model; + CvarQueryInformation const& queryInformation; + storm::storage::BitVector const& targetStates; +}; + +} // namespace cvar +} // namespace modelchecker +} // namespace storm diff --git a/src/storm/modelchecker/cvar/helper/SparseSspCvarParetoViHelper.h b/src/storm/modelchecker/cvar/helper/SparseSspCvarParetoViHelper.h new file mode 100644 index 0000000000..ffb0471005 --- /dev/null +++ b/src/storm/modelchecker/cvar/helper/SparseSspCvarParetoViHelper.h @@ -0,0 +1,133 @@ +#pragma once + +#include +#include +#include +#include + +#include "storm/environment/Environment.h" +#include "storm/exceptions/NotImplementedException.h" +#include "storm/exceptions/UnexpectedException.h" +#include "storm/modelchecker/cvar/CvarComputationResult.h" +#include "storm/modelchecker/cvar/CvarQueryInformation.h" +#include "storm/modelchecker/cvar/helper/SspParetoValueIterationOperator.h" +#include "storm/modelchecker/cvar/preprocessing/SspCvarPreprocessingResult.h" +#include "storm/utility/constants.h" +#include "storm/utility/macros.h" + +namespace storm { +namespace modelchecker { +namespace cvar { + +/*! + * Implements the Pareto-front value iteration for Conditional Value-at-Risk on an SSP-style total-cost objective. + * + * Supported CLI shape: + * storm --prism model.nm --prop 'R{"reward"}min/max=? [ F "target" ]' --cvar --cvar:method ssp + * + * The helper assumes that SSP-specific preprocessing has already normalized the sparse MDP to terminal-goal + * semantics, validated integer-valued positive choice costs outside goal states, and computed the classical + * expected cost-to-go vector used in the base layer. + * + * @see https://doi.org/10.1609/aaai.v36i9.21222 for the Pareto-front value iteration on which this implementation is based. + */ +template +class SparseSspCvarParetoViHelper { + public: + using ParetoFront = SspParetoFront; + using ParetoViOperator = SspParetoValueIterationOperator; + using FrontierLayer = std::vector; + using FrontierWindow = std::vector; + + SparseSspCvarParetoViHelper(CvarQueryInformation const& queryInformation, preprocessing::SspCvarPreprocessingResult const& preprocessingResult) + : queryInformation(queryInformation), preprocessingResult(preprocessingResult), paretoViOperator(preprocessingResult) { + // Intentionally left empty. + } + + CvarComputationResult computeCvar(Environment const&, bool produceScheduler = false) const { + STORM_LOG_THROW(!produceScheduler, storm::exceptions::NotImplementedException, + "Scheduler extraction for CVaR SSP value iteration is not implemented yet."); + + FrontierWindow frontierWindow = initializeFrontierWindow(); + std::optional bestCandidate = + extractCvarCandidateFromInitialFrontier(frontierWindow[0][preprocessingResult.initialState], 0, queryInformation.alpha); + FrontierLayer currentLayer(paretoViOperator.getStateCount()); + + for (uint64_t costBound = 1; !bestCandidate.has_value() || storm::utility::convertNumber(costBound) <= bestCandidate.value(); ++costBound) { + paretoViOperator.apply(costBound, frontierWindow, currentLayer); + auto currentCandidate = extractCvarCandidateFromInitialFrontier(currentLayer[preprocessingResult.initialState], costBound, queryInformation.alpha); + if (currentCandidate.has_value() && (!bestCandidate.has_value() || currentCandidate.value() < bestCandidate.value())) { + bestCandidate = currentCandidate; + } + swapFrontierLayerIntoWindow(costBound, currentLayer, frontierWindow); + } + + STORM_LOG_THROW(bestCandidate.has_value(), storm::exceptions::UnexpectedException, "CVaR SSP value iteration did not find a feasible candidate."); + return {bestCandidate.value(), nullptr}; + } + + private: + /*! + * Creates the initial frontier layer for a given bound n. + * + * For n < 0 we intentionally deviate from the paper's P_n = empty convention and instead use the + * mathematically direct excess-cost semantics induced by E[(X-n)^+]: because all costs are nonnegative, + * Pr[X <= n] = 0 and E[(X-n)^+] = E[X] - n. Hence the Pareto set collapses to the singleton + * (0, e*(s) - n), where e*(s) is the minimal expected cost-to-go from s. + */ + FrontierLayer createInitialFrontierLayer(int64_t costBound) const { + FrontierLayer baseLayer(paretoViOperator.getStateCount()); + ValueType const boundValue = storm::utility::convertNumber(costBound); + if (costBound >= 0) { + ParetoFront const targetFront = ParetoViOperator::createTargetFrontier(); + for (auto state : paretoViOperator.getReachableTargetStates()) { + baseLayer[state] = targetFront; + } + } else { + for (auto state : paretoViOperator.getReachableTargetStates()) { + baseLayer[state] = ParetoFront::singleton(storm::utility::zero(), preprocessingResult.expectedCostsToGoal[state] - boundValue); + } + } + for (auto state : paretoViOperator.getReachableNonTargetStates()) { + baseLayer[state] = ParetoFront::singleton(storm::utility::zero(), preprocessingResult.expectedCostsToGoal[state] - boundValue); + } + return baseLayer; + } + + FrontierWindow initializeFrontierWindow() const { + STORM_LOG_ASSERT(preprocessingResult.maximalChoiceCost > 0, "Expected a strictly positive maximal choice cost."); + FrontierWindow frontierWindow(preprocessingResult.maximalChoiceCost, FrontierLayer(paretoViOperator.getStateCount())); + for (int64_t costBound = 1 - static_cast(preprocessingResult.maximalChoiceCost); costBound <= 0; ++costBound) { + frontierWindow[ParetoViOperator::getWindowIndex(costBound, frontierWindow.size())] = createInitialFrontierLayer(costBound); + } + return frontierWindow; + } + + static void swapFrontierLayerIntoWindow(int64_t costBound, FrontierLayer& layer, FrontierWindow& frontierWindow) { + std::swap(frontierWindow[ParetoViOperator::getWindowIndex(costBound, frontierWindow.size())], layer); + } + + /*! + * Evaluates one initial-state frontier for a fixed cost bound n. + * + * Following Algorithm 1 and the total-cost extension in the paper, a frontier for bound n induces the CVaR + * candidate n + E / t, where E is the minimal continuation cost on the frontier at probability 1 - t. + */ + static std::optional extractCvarCandidateFromInitialFrontier(SspParetoFront const& initialFrontier, uint64_t costBound, + storm::RationalNumber const& alpha) { + ValueType const targetProbability = storm::utility::one() - storm::utility::convertNumber(alpha); + auto continuationCost = initialFrontier.getMinimalContinuationCostAtProbability(targetProbability); + if (!continuationCost.has_value()) { + return std::nullopt; + } + return storm::utility::convertNumber(costBound) + continuationCost.value() / storm::utility::convertNumber(alpha); + } + + CvarQueryInformation const& queryInformation; + preprocessing::SspCvarPreprocessingResult const& preprocessingResult; + ParetoViOperator paretoViOperator; +}; + +} // namespace cvar +} // namespace modelchecker +} // namespace storm diff --git a/src/storm/modelchecker/cvar/helper/SparseSspRewardCvarParetoViHelper.h b/src/storm/modelchecker/cvar/helper/SparseSspRewardCvarParetoViHelper.h new file mode 100644 index 0000000000..358d2da58c --- /dev/null +++ b/src/storm/modelchecker/cvar/helper/SparseSspRewardCvarParetoViHelper.h @@ -0,0 +1,147 @@ +#pragma once + +#include +#include +#include +#include + +#include "storm/environment/Environment.h" +#include "storm/exceptions/NotImplementedException.h" +#include "storm/exceptions/UnexpectedException.h" +#include "storm/modelchecker/cvar/CvarComputationResult.h" +#include "storm/modelchecker/cvar/CvarQueryInformation.h" +#include "storm/modelchecker/cvar/helper/SspParetoValueIterationOperator.h" +#include "storm/modelchecker/cvar/preprocessing/SspCvarPreprocessingResult.h" +#include "storm/utility/constants.h" +#include "storm/utility/macros.h" + +namespace storm { +namespace modelchecker { +namespace cvar { + +/*! + * Implements Pareto-front value iteration for lower-tail CVaR on an SSP-style total-reward objective. + * + * The helper assumes that preprocessing has normalized the model to terminal-goal semantics and validated + * strictly positive integer per-choice rewards outside goal states. The computed frontier stores tuples + * (p, e), where p is the probability of accumulating reward at most the current threshold and e is the + * expected shortfall below that threshold. + */ +template +class SparseSspRewardCvarParetoViHelper { + public: + using ParetoFront = SspRewardParetoFront; + using ParetoViOperator = SspParetoValueIterationOperator; + using FrontierLayer = std::vector; + using FrontierWindow = std::vector; + + SparseSspRewardCvarParetoViHelper(CvarQueryInformation const& queryInformation, + preprocessing::SspCvarPreprocessingResult const& preprocessingResult) + : queryInformation(queryInformation), preprocessingResult(preprocessingResult), paretoViOperator(preprocessingResult) { + // Intentionally left empty. + } + + CvarComputationResult computeCvar(Environment const&, bool produceScheduler = false) const { + STORM_LOG_THROW(!produceScheduler, storm::exceptions::NotImplementedException, + "Scheduler extraction for CVaR SSP value iteration is not implemented yet."); + + ValueType const alpha = storm::utility::convertNumber(queryInformation.alpha); + FrontierWindow frontierWindow = initializeFrontierWindow(); + auto bestCandidate = extractBestCvarCandidateFromInitialFrontier(frontierWindow[0][preprocessingResult.initialState], 0, alpha); + if (initialFrontierHasReachedAlpha(frontierWindow[0][preprocessingResult.initialState], alpha)) { + STORM_LOG_THROW(bestCandidate.has_value(), storm::exceptions::UnexpectedException, + "CVaR SSP reward value iteration reached the stopping criterion without a candidate."); + return {bestCandidate.value(), nullptr}; + } + + FrontierLayer currentLayer(paretoViOperator.getStateCount()); + for (uint64_t rewardThreshold = 1;; ++rewardThreshold) { + paretoViOperator.apply(rewardThreshold, frontierWindow, currentLayer); + auto currentCandidate = extractBestCvarCandidateFromInitialFrontier(currentLayer[preprocessingResult.initialState], rewardThreshold, alpha); + if (currentCandidate.has_value() && (!bestCandidate.has_value() || bestCandidate.value() < currentCandidate.value())) { + bestCandidate = currentCandidate; + } + + bool const reachedAlpha = initialFrontierHasReachedAlpha(currentLayer[preprocessingResult.initialState], alpha); + swapFrontierLayerIntoWindow(rewardThreshold, currentLayer, frontierWindow); + if (reachedAlpha) { + break; + } + } + + STORM_LOG_THROW(bestCandidate.has_value(), storm::exceptions::UnexpectedException, + "CVaR SSP reward value iteration did not find a feasible candidate."); + return {bestCandidate.value(), nullptr}; + } + + private: + FrontierLayer createInitialFrontierLayer(int64_t rewardThreshold) const { + FrontierLayer baseLayer(paretoViOperator.getStateCount()); + if (rewardThreshold >= 0) { + ParetoFront const targetFront = ParetoFront::singleton(storm::utility::one(), storm::utility::convertNumber(rewardThreshold)); + for (auto state : paretoViOperator.getReachableTargetStates()) { + baseLayer[state] = targetFront; + } + } else { + ParetoFront const targetFront = ParetoFront::singleton(storm::utility::zero(), storm::utility::zero()); + for (auto state : paretoViOperator.getReachableTargetStates()) { + baseLayer[state] = targetFront; + } + } + ParetoFront const nonTargetFront = ParetoFront::singleton(storm::utility::zero(), storm::utility::zero()); + for (auto state : paretoViOperator.getReachableNonTargetStates()) { + baseLayer[state] = nonTargetFront; + } + return baseLayer; + } + + FrontierWindow initializeFrontierWindow() const { + STORM_LOG_ASSERT(preprocessingResult.maximalChoiceCost > 0, "Expected a strictly positive maximal choice reward."); + FrontierWindow frontierWindow(preprocessingResult.maximalChoiceCost, FrontierLayer(paretoViOperator.getStateCount())); + for (int64_t rewardThreshold = 1 - static_cast(preprocessingResult.maximalChoiceCost); rewardThreshold <= 0; ++rewardThreshold) { + frontierWindow[ParetoViOperator::getWindowIndex(rewardThreshold, frontierWindow.size())] = createInitialFrontierLayer(rewardThreshold); + } + return frontierWindow; + } + + static void swapFrontierLayerIntoWindow(int64_t rewardThreshold, FrontierLayer& layer, FrontierWindow& frontierWindow) { + std::swap(frontierWindow[ParetoViOperator::getWindowIndex(rewardThreshold, frontierWindow.size())], layer); + } + + static std::optional extractBestCvarCandidateFromInitialFrontier(ParetoFront const& initialFrontier, uint64_t rewardThreshold, + ValueType const& alpha) { + if (initialFrontier.empty()) { + return std::nullopt; + } + ValueType const thresholdValue = storm::utility::convertNumber(rewardThreshold); + std::optional result; + for (auto const& point : initialFrontier) { + ValueType const candidate = thresholdValue - point.expectedCost / alpha; + if (!result.has_value() || result.value() < candidate) { + result = candidate; + } + } + return result; + } + + static bool initialFrontierHasReachedAlpha(ParetoFront const& initialFrontier, ValueType const& alpha) { + if (initialFrontier.empty()) { + return false; + } + storm::utility::ElementLess less; + for (auto const& point : initialFrontier) { + if (less(point.probability, alpha)) { + return false; + } + } + return true; + } + + CvarQueryInformation const& queryInformation; + preprocessing::SspCvarPreprocessingResult const& preprocessingResult; + ParetoViOperator paretoViOperator; +}; + +} // namespace cvar +} // namespace modelchecker +} // namespace storm diff --git a/src/storm/modelchecker/cvar/helper/SparseWeightedReachabilityCvarLpHelper.h b/src/storm/modelchecker/cvar/helper/SparseWeightedReachabilityCvarLpHelper.h new file mode 100644 index 0000000000..7d92f0817e --- /dev/null +++ b/src/storm/modelchecker/cvar/helper/SparseWeightedReachabilityCvarLpHelper.h @@ -0,0 +1,453 @@ +#pragma once + +#include +#include +#include +#include +#include + +#include "storm/adapters/RationalNumberAdapter.h" +#include "storm/environment/Environment.h" +#include "storm/environment/solver/MinMaxSolverEnvironment.h" +#include "storm/environment/solver/SolverEnvironment.h" +#include "storm/exceptions/NotImplementedException.h" +#include "storm/exceptions/UnexpectedException.h" +#include "storm/modelchecker/cvar/CvarComputationResult.h" +#include "storm/modelchecker/cvar/CvarQueryInformation.h" +#include "storm/modelchecker/cvar/preprocessing/WeightedReachabilityCvarPreprocessingResult.h" +#include "storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.h" +#include "storm/solver/LpSolver.h" +#include "storm/solver/SolveGoal.h" +#include "storm/storage/Scheduler.h" +#include "storm/storage/expressions/BinaryRelationType.h" +#include "storm/utility/ConstantsComparator.h" +#include "storm/utility/constants.h" +#include "storm/utility/graph.h" +#include "storm/utility/macros.h" +#include "storm/utility/solver.h" + +namespace storm { +namespace modelchecker { +namespace cvar { + +template +struct CvarThresholdData { + ValueType threshold; + storm::storage::BitVector targetStatesInStrictTail; + storm::storage::BitVector targetStatesAtThreshold; + storm::storage::BitVector targetStatesInTail; +}; + +template +struct CvarRewardBucket { + ValueType reward; + std::vector targetStates; +}; + +template +struct WeightedReachabilityCvarLpData { + storm::RationalNumber alpha; + storm::solver::OptimizationDirection optimizationDirection; + CvarInterpretation interpretation; + uint64_t initialState; + storm::storage::BitVector initialStates; + std::string rewardModelName; + storm::storage::BitVector targetStates; + std::vector terminalRewards; + std::vector> rewardBuckets; + storm::storage::SparseMatrix transitionMatrix; + storm::storage::SparseMatrix backwardChoices; + storm::storage::SparseMatrix backwardTransitions; +}; + +template +std::vector> collectRewardBuckets(storm::storage::BitVector const& targetStates, std::vector const& terminalRewards, + storm::storage::BitVector const& reachableStates) { + std::map> buckets; + for (auto state : targetStates) { + if (!reachableStates[state]) { + continue; + } + buckets[terminalRewards[state]].push_back(state); + } + + std::vector> result; + result.reserve(buckets.size()); + for (auto& bucket : buckets) { + result.push_back({bucket.first, std::move(bucket.second)}); + } + return result; +} + +/*! + * Solves an LP for Conditional Value-at-Risk on an MDP with a terminal reward objective. + * + * Supported CLI shape: + * storm --prism model.nm --prop 'R{"reward"}min/max=? [ F "target" ]' --cvar + * storm --prism model.nm --prop 'R{"reward"}min/max=? [ F "target" ]' --cvar --cvar:method wr + * + * The --cvar option requires exactly one selected property. That property must be an unbounded + * reward query with an optimization direction (min or max) and an eventually formula F phi whose target phi is a state + * formula. The reward model may be named explicitly (R{"reward"}...) or omitted if the model has a unique reward model. + * + * Requirements: sparse MDP, 0 < alpha < 1, exactly one initial state, state-based terminal rewards, + * reward 0 on non-target states, and absorbing original target states. + * + * @see https://doi.org/10.1145/3209108.3209176 Fig. 4 for a description of the algorithm as implemented (and slightly altered) here. + */ +template +class SparseWeightedReachabilityCvarLpHelper { + public: + SparseWeightedReachabilityCvarLpHelper(CvarQueryInformation const& queryInformation, + preprocessing::WeightedReachabilityCvarPreprocessingResult const& weightedReachabilityPreprocessingResult) + : lpData(createLpData(queryInformation, weightedReachabilityPreprocessingResult)) {} + + CvarComputationResult computeCvar(Environment const& env, bool produceScheduler = false) const { + STORM_LOG_THROW(!lpData.rewardBuckets.empty(), storm::exceptions::NotImplementedException, + "CVaR model checking requires at least one target reward threshold candidate."); + + std::optional bestValue; + std::optional bestThresholdIndex; + std::unique_ptr> bestScheduler; + auto candidateRange = computeCandidateRange(env); + for (uint64_t thresholdIndex = candidateRange.first; thresholdIndex < candidateRange.second; ++thresholdIndex) { + auto thresholdData = createThresholdData(thresholdIndex); + auto thresholdResult = buildLpForThreshold(env, thresholdData, false); + if (!thresholdResult.has_value()) { + continue; + } + if (!bestValue.has_value()) { + bestValue = thresholdResult->value; + bestThresholdIndex = thresholdIndex; + } else if ((storm::solver::minimize(lpData.optimizationDirection) && thresholdResult->value < bestValue.value()) || + (storm::solver::maximize(lpData.optimizationDirection) && thresholdResult->value > bestValue.value())) { + bestValue = thresholdResult->value; + bestThresholdIndex = thresholdIndex; + } + } + + STORM_LOG_THROW(bestValue.has_value(), storm::exceptions::UnexpectedException, + "CVaR model checking did not find a feasible LP for any threshold candidate."); + if (produceScheduler) { + STORM_LOG_ASSERT(bestThresholdIndex.has_value(), "Expected a threshold index for the best CVaR LP value."); + auto thresholdData = createThresholdData(bestThresholdIndex.value()); + auto thresholdResult = buildLpForThreshold(env, thresholdData, true); + STORM_LOG_THROW(thresholdResult.has_value(), storm::exceptions::UnexpectedException, + "The previously optimal CVaR LP threshold became infeasible when extracting a scheduler."); + bestScheduler = std::move(thresholdResult->scheduler); + } + return {bestValue.value() / storm::utility::convertNumber(lpData.alpha), std::move(bestScheduler)}; + } + + private: + static WeightedReachabilityCvarLpData createLpData( + CvarQueryInformation const& queryInformation, + preprocessing::WeightedReachabilityCvarPreprocessingResult const& weightedReachabilityPreprocessingResult) { + auto initialStates = createInitialStateBitVector(weightedReachabilityPreprocessingResult.transitionMatrix.getRowGroupCount(), + weightedReachabilityPreprocessingResult.initialState); + storm::storage::BitVector allStates(weightedReachabilityPreprocessingResult.transitionMatrix.getRowGroupCount(), true); + storm::storage::BitVector noStates(weightedReachabilityPreprocessingResult.transitionMatrix.getRowGroupCount(), false); + auto reachableStates = + storm::utility::graph::getReachableStates(weightedReachabilityPreprocessingResult.transitionMatrix, initialStates, allStates, noStates); + auto rewardBuckets = collectRewardBuckets(weightedReachabilityPreprocessingResult.effectiveTargetStates, + weightedReachabilityPreprocessingResult.terminalRewards, reachableStates); + return {queryInformation.alpha, + queryInformation.optimizationDirection, + queryInformation.interpretation, + weightedReachabilityPreprocessingResult.initialState, + std::move(initialStates), + weightedReachabilityPreprocessingResult.rewardModelName, + weightedReachabilityPreprocessingResult.effectiveTargetStates, + weightedReachabilityPreprocessingResult.terminalRewards, + std::move(rewardBuckets), + weightedReachabilityPreprocessingResult.transitionMatrix, + weightedReachabilityPreprocessingResult.transitionMatrix.transpose(), + weightedReachabilityPreprocessingResult.transitionMatrix.transpose(true)}; + } + + static storm::storage::BitVector createInitialStateBitVector(uint64_t stateCount, uint64_t initialState) { + storm::storage::BitVector initialStates(stateCount, false); + initialStates.set(initialState, true); + return initialStates; + } + + void addBucketStates(storm::storage::BitVector& states, uint64_t bucketIndex) const { + for (auto state : lpData.rewardBuckets[bucketIndex].targetStates) { + states.set(state, true); + } + } + + storm::storage::BitVector createBucketTargetStates(uint64_t bucketIndex) const { + storm::storage::BitVector states(lpData.transitionMatrix.getRowGroupCount(), false); + addBucketStates(states, bucketIndex); + return states; + } + + storm::storage::BitVector createPrefixTargetStates(uint64_t endBucketIndex) const { + storm::storage::BitVector states(lpData.transitionMatrix.getRowGroupCount(), false); + for (uint64_t bucketIndex = 0; bucketIndex < endBucketIndex; ++bucketIndex) { + addBucketStates(states, bucketIndex); + } + return states; + } + + storm::storage::BitVector createSuffixTargetStates(uint64_t startBucketIndex) const { + storm::storage::BitVector states(lpData.transitionMatrix.getRowGroupCount(), false); + for (uint64_t bucketIndex = startBucketIndex; bucketIndex < lpData.rewardBuckets.size(); ++bucketIndex) { + addBucketStates(states, bucketIndex); + } + return states; + } + + storm::storage::BitVector const& getCachedPrefixTargetStates(uint64_t endBucketIndex, std::map& cache) const { + auto cachedPrefix = cache.find(endBucketIndex); + if (cachedPrefix != cache.end()) { + return cachedPrefix->second; + } + return cache.emplace(endBucketIndex, createPrefixTargetStates(endBucketIndex)).first->second; + } + + storm::storage::BitVector const& getCachedSuffixTargetStates(uint64_t startBucketIndex, std::map& cache) const { + auto cachedSuffix = cache.find(startBucketIndex); + if (cachedSuffix != cache.end()) { + return cachedSuffix->second; + } + return cache.emplace(startBucketIndex, createSuffixTargetStates(startBucketIndex)).first->second; + } + + CvarThresholdData createThresholdData(uint64_t thresholdIndex) const { + auto targetStatesAtThreshold = createBucketTargetStates(thresholdIndex); + if (lpData.interpretation == CvarInterpretation::Reward) { + auto targetStatesInStrictTail = createPrefixTargetStates(thresholdIndex); + auto targetStatesInTail = targetStatesInStrictTail | targetStatesAtThreshold; + return {lpData.rewardBuckets[thresholdIndex].reward, std::move(targetStatesInStrictTail), std::move(targetStatesAtThreshold), + std::move(targetStatesInTail)}; + } + + auto targetStatesInStrictTail = createSuffixTargetStates(thresholdIndex + 1); + auto targetStatesInTail = targetStatesInStrictTail | targetStatesAtThreshold; + return {lpData.rewardBuckets[thresholdIndex].reward, std::move(targetStatesInStrictTail), std::move(targetStatesAtThreshold), + std::move(targetStatesInTail)}; + } + + ValueType computeReachabilityProbability(Environment const& env, storm::solver::OptimizationDirection direction, + storm::storage::BitVector const& targetStates) const { + if (targetStates.empty()) { + return storm::utility::zero(); + } + if (targetStates[lpData.initialState]) { + return storm::utility::one(); + } + + storm::storage::BitVector allStates(lpData.transitionMatrix.getRowGroupCount(), true); + auto result = storm::modelchecker::helper::SparseMdpPrctlHelper::computeUntilProbabilities( + env, storm::solver::SolveGoal(direction, lpData.initialStates), lpData.transitionMatrix, lpData.backwardTransitions, + allStates, targetStates, false, false); + return result.values[lpData.initialState]; + } + + std::pair computeCandidateRange(Environment const& env) const { + uint64_t const bucketCount = lpData.rewardBuckets.size(); + ValueType const alpha = storm::utility::convertNumber(lpData.alpha); + storm::utility::ConstantsComparator comparator(storm::utility::convertNumber(env.solver().minMax().getPrecision())); + + if (lpData.interpretation == CvarInterpretation::Reward) { + std::map prefixTargetStateCache; + + uint64_t lower = 0; + uint64_t upper = bucketCount; + while (lower < upper) { + uint64_t const mid = lower + (upper - lower) / 2; + auto const& targetStatesInTail = getCachedPrefixTargetStates(mid + 1, prefixTargetStateCache); + auto maxReachability = computeReachabilityProbability(env, storm::solver::OptimizationDirection::Maximize, targetStatesInTail); + if (comparator.isLess(maxReachability, alpha)) { + lower = mid + 1; + } else { + upper = mid; + } + } + uint64_t const firstNotTooLow = lower; + + lower = firstNotTooLow; + upper = bucketCount; + while (lower < upper) { + uint64_t const mid = lower + (upper - lower) / 2; + auto const& targetStatesInStrictTail = getCachedPrefixTargetStates(mid, prefixTargetStateCache); + auto minReachability = computeReachabilityProbability(env, storm::solver::OptimizationDirection::Minimize, targetStatesInStrictTail); + if (comparator.isLess(alpha, minReachability)) { + upper = mid; + } else { + lower = mid + 1; + } + } + + return {firstNotTooLow, lower}; + } + + std::map suffixTargetStateCache; + + uint64_t lower = 0; + uint64_t upper = bucketCount; + while (lower < upper) { + uint64_t const mid = lower + (upper - lower) / 2; + auto const& targetStatesInStrictTail = getCachedSuffixTargetStates(mid + 1, suffixTargetStateCache); + auto minReachability = computeReachabilityProbability(env, storm::solver::OptimizationDirection::Minimize, targetStatesInStrictTail); + if (comparator.isLess(alpha, minReachability)) { + lower = mid + 1; + } else { + upper = mid; + } + } + uint64_t const firstNotTooLow = lower; + + lower = firstNotTooLow; + upper = bucketCount; + while (lower < upper) { + uint64_t const mid = lower + (upper - lower) / 2; + auto const& targetStatesInTail = getCachedSuffixTargetStates(mid, suffixTargetStateCache); + auto maxReachability = computeReachabilityProbability(env, storm::solver::OptimizationDirection::Maximize, targetStatesInTail); + if (comparator.isLess(maxReachability, alpha)) { + upper = mid; + } else { + lower = mid + 1; + } + } + + return {firstNotTooLow, lower}; + } + + std::optional> buildLpForThreshold(Environment const& env, CvarThresholdData const& thresholdData, + bool produceScheduler) const { + using RawLpSolver = storm::solver::LpSolver; + using RawLpConstraint = storm::solver::RawLpConstraint; + + auto lpSolverFactory = storm::utility::solver::getLpSolverFactory(env); + auto solver = lpSolverFactory->createRaw(env, "cvar"); + solver->setOptimizationDirection(lpData.optimizationDirection); + + std::vector actionFlowVariables; + actionFlowVariables.reserve(lpData.transitionMatrix.getRowCount()); + for (uint64_t row = 0; row < lpData.transitionMatrix.getRowCount(); ++row) { + actionFlowVariables.push_back(solver->addLowerBoundedContinuousVariable("y_" + std::to_string(row), storm::utility::zero())); + } + + std::vector> recurrentFlowVariables(lpData.transitionMatrix.getRowGroupCount(), std::nullopt); + for (uint64_t state = 0; state < lpData.transitionMatrix.getRowGroupCount(); ++state) { + if (lpData.targetStates[state]) { + recurrentFlowVariables[state] = solver->addLowerBoundedContinuousVariable("x_" + std::to_string(state), storm::utility::zero()); + } + } + + std::vector> splitFlowVariables(lpData.transitionMatrix.getRowGroupCount(), std::nullopt); + for (uint64_t state = 0; state < lpData.transitionMatrix.getRowGroupCount(); ++state) { + if (thresholdData.targetStatesInTail[state]) { + splitFlowVariables[state] = + solver->addLowerBoundedContinuousVariable("xb_" + std::to_string(state), storm::utility::zero(), lpData.terminalRewards[state]); + } + } + + solver->update(); + + for (uint64_t state = 0; state < lpData.transitionMatrix.getRowGroupCount(); ++state) { + auto outgoingActions = lpData.transitionMatrix.getRowGroupIndices(state); + auto incomingActions = lpData.backwardChoices.getRow(state); + uint64_t reservedSize = outgoingActions.size() + incomingActions.getNumberOfEntries() + (recurrentFlowVariables[state].has_value() ? 1 : 0); + RawLpConstraint constraint(storm::expressions::RelationType::Equal, + state == lpData.initialState ? storm::utility::one() : storm::utility::zero(), reservedSize); + std::map actionCoefficients; + + for (auto const& incomingAction : incomingActions) { + actionCoefficients[actionFlowVariables[incomingAction.getColumn()]] -= incomingAction.getValue(); + } + for (auto const& action : outgoingActions) { + actionCoefficients[actionFlowVariables[action]] += storm::utility::one(); + } + for (auto const& actionCoefficient : actionCoefficients) { + if (!storm::utility::isZero(actionCoefficient.second)) { + constraint.addToLhs(actionCoefficient.first, actionCoefficient.second); + } + } + if (recurrentFlowVariables[state].has_value()) { + constraint.addToLhs(recurrentFlowVariables[state].value(), storm::utility::one()); + } + + solver->addConstraint("transient_flow_" + std::to_string(state), constraint); + } + + RawLpConstraint recurrentConstraint(storm::expressions::RelationType::Equal, storm::utility::one(), + lpData.targetStates.getNumberOfSetBits()); + + for (auto state : lpData.targetStates) { + recurrentConstraint.addToLhs(recurrentFlowVariables[state].value(), storm::utility::one()); + } + solver->addConstraint("recurrent_behaviour", recurrentConstraint); + + for (auto state : thresholdData.targetStatesInStrictTail) { + RawLpConstraint splitEqualityConstraint(storm::expressions::RelationType::Equal, storm::utility::zero(), 2); + splitEqualityConstraint.addToLhs(splitFlowVariables[state].value(), storm::utility::one()); + splitEqualityConstraint.addToLhs(recurrentFlowVariables[state].value(), -storm::utility::one()); + solver->addConstraint("split_eq_" + std::to_string(state), splitEqualityConstraint); + } + for (auto state : thresholdData.targetStatesAtThreshold) { + RawLpConstraint splitInequalityConstraint(storm::expressions::RelationType::LessOrEqual, storm::utility::zero(), 2); + splitInequalityConstraint.addToLhs(splitFlowVariables[state].value(), storm::utility::one()); + splitInequalityConstraint.addToLhs(recurrentFlowVariables[state].value(), -storm::utility::one()); + solver->addConstraint("split_le_" + std::to_string(state), splitInequalityConstraint); + } + + RawLpConstraint probabilityConsistentSplitConstraint(storm::expressions::RelationType::Equal, storm::utility::convertNumber(lpData.alpha), + thresholdData.targetStatesInTail.getNumberOfSetBits()); + for (auto state : thresholdData.targetStatesInTail) { + probabilityConsistentSplitConstraint.addToLhs(splitFlowVariables[state].value(), storm::utility::one()); + } + solver->addConstraint("probability_consistent_split", probabilityConsistentSplitConstraint); + + solver->optimize(); + if (solver->isInfeasible()) { + return std::nullopt; + } + STORM_LOG_THROW(!solver->isUnbounded(), storm::exceptions::UnexpectedException, + "The CVaR LP for threshold " << thresholdData.threshold << " is unbounded."); + STORM_LOG_THROW(solver->isOptimal(), storm::exceptions::UnexpectedException, + "The CVaR LP for threshold " << thresholdData.threshold << " did not reach an optimal solution."); + + std::unique_ptr> scheduler; + if (produceScheduler) { + scheduler = std::make_unique>(lpData.transitionMatrix.getRowGroupCount()); + for (uint64_t state = 0; state < lpData.transitionMatrix.getRowGroupCount(); ++state) { + if (lpData.targetStates[state]) { + scheduler->setDontCare(state); + continue; + } + + uint64_t firstRow = lpData.transitionMatrix.getRowGroupIndices()[state]; + uint64_t lastRow = lpData.transitionMatrix.getRowGroupIndices()[state + 1]; + storm::storage::Distribution actionDistribution; + actionDistribution.reserve(lastRow - firstRow); + + for (uint64_t row = firstRow; row < lastRow; ++row) { + auto flow = solver->getContinuousValue(actionFlowVariables[row]); + if (storm::utility::isAlmostZero(flow)) { + continue; + } + actionDistribution.addProbability(row - firstRow, flow); + } + + if (actionDistribution.size() == 0) { + scheduler->setDontCare(state); + continue; + } + + actionDistribution.normalize(); + scheduler->setChoice(storm::storage::SchedulerChoice(std::move(actionDistribution)), state); + } + } + + return CvarComputationResult{solver->getObjectiveValue(), std::move(scheduler)}; + } + + WeightedReachabilityCvarLpData lpData; +}; +} // namespace cvar +} // namespace modelchecker +} // namespace storm diff --git a/src/storm/modelchecker/cvar/helper/SspParetoFront.h b/src/storm/modelchecker/cvar/helper/SspParetoFront.h new file mode 100644 index 0000000000..2ec90c04e8 --- /dev/null +++ b/src/storm/modelchecker/cvar/helper/SspParetoFront.h @@ -0,0 +1,573 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "storm/utility/constants.h" +#include "storm/utility/macros.h" + +namespace storm { +namespace modelchecker { +namespace cvar { + +enum class SspParetoFrontKind { CostUpperTail, RewardLowerTail }; + +/*! + * Represents one SSP CVaR Pareto frontier. + * + * For the cost-side upper-tail algorithm, each point (p, E) represents an achievable tradeoff where: + * - p is the probability of reaching the goal within the current cost bound, and + * - E is the expected continuation cost beyond that bound. + * + * For the reward-side lower-tail algorithm, each point (p, e) represents an achievable tradeoff where: + * - p is the probability of accumulating reward at most the current threshold, and + * - e is the expected shortfall below that threshold. + * + * The class only stores extremal boundary points. The closed polygon itself is induced by these points together + * with convex closure. Cost frontiers use the original lower-right order; reward frontiers use lower-left order. + */ +template +class SspParetoFront { + public: + struct Point { + //! Probability of reaching the goal within the current cost bound. + ValueType probability; + //! Expected continuation cost beyond the current cost bound. + ValueType expectedCost; + + enum class DominanceResult { Incomparable, Dominates, Dominated, Equal }; + + DominanceResult getDominance(Point const& other) const { + if (probability == other.probability && expectedCost == other.expectedCost) { + return DominanceResult::Equal; + } + if constexpr (FrontKind == SspParetoFrontKind::CostUpperTail) { + if (probability >= other.probability && expectedCost <= other.expectedCost) { + return DominanceResult::Dominates; + } + if (probability <= other.probability && expectedCost >= other.expectedCost) { + return DominanceResult::Dominated; + } + } else { + if (probability <= other.probability && expectedCost <= other.expectedCost) { + return DominanceResult::Dominates; + } + if (probability >= other.probability && expectedCost >= other.expectedCost) { + return DominanceResult::Dominated; + } + } + return DominanceResult::Incomparable; + } + }; + + using container_type = std::vector; + using const_iterator = typename container_type::const_iterator; + + private: + struct AlreadyCanonicalTag {}; + struct AlreadySortedTag {}; + + explicit SspParetoFront(container_type points, AlreadyCanonicalTag) : points(std::move(points)) { + // Intentionally left empty. + } + + explicit SspParetoFront(container_type points, AlreadySortedTag) : points(std::move(points)) { + canonicalizeSorted(); + } + + public: + SspParetoFront() = default; + + explicit SspParetoFront(container_type points) : points(std::move(points)) { + canonicalize(); + } + + static SspParetoFront singleton(ValueType const& probability, ValueType const& expectedCost) { + return SspParetoFront(container_type{{probability, expectedCost}}, AlreadyCanonicalTag{}); + } + + bool empty() const { + return points.empty(); + } + + std::size_t size() const { + return points.size(); + } + + bool isSingleton() const { + return points.size() == 1; + } + + container_type const& getPoints() const { + return points; + } + + const_iterator begin() const { + return points.begin(); + } + + const_iterator end() const { + return points.end(); + } + + void clear() { + points.clear(); + } + + void addPoint(Point const& point) { + addPoint(Point{point}); + } + + void addPoint(Point&& point) { + auto it = points.begin(); + while (it != points.end()) { + switch (point.getDominance(*it)) { + case Point::DominanceResult::Equal: + case Point::DominanceResult::Dominated: + return; + case Point::DominanceResult::Dominates: + it = points.erase(it); + break; + case Point::DominanceResult::Incomparable: + ++it; + break; + } + } + points.push_back(std::move(point)); + canonicalize(); + } + + void addPoints(container_type const& additionalPoints) { + points.insert(points.end(), additionalPoints.begin(), additionalPoints.end()); + canonicalize(); + } + + SspParetoFront scaled(ValueType const& factor) const { + if (empty()) { + return SspParetoFront(); + } + if (storm::utility::isOne(factor)) { + return *this; + } + if (storm::utility::isZero(factor)) { + return singleton(storm::utility::zero(), storm::utility::zero()); + } + if (isSingleton()) { + auto const& point = points.front(); + return singleton(factor * point.probability, factor * point.expectedCost); + } + container_type scaledPoints; + scaledPoints.reserve(points.size()); + for (auto const& point : points) { + scaledPoints.push_back(Point{factor * point.probability, factor * point.expectedCost}); + } + if (storm::utility::isPositive(factor)) { + return SspParetoFront(std::move(scaledPoints), AlreadySortedTag{}); + } + return SspParetoFront(std::move(scaledPoints)); + } + + SspParetoFront minkowskiSum(SspParetoFront const& other) const { + if (empty() || other.empty()) { + return SspParetoFront(); + } + if (isSingleton()) { + return other.translated(points.front()); + } + if (other.isSingleton()) { + return translated(other.points.front()); + } + return minkowskiSumConvexChain(other, storm::utility::one()); + } + + SspParetoFront minkowskiSumScaled(SspParetoFront const& other, ValueType const& factor) const { + if (empty() || other.empty()) { + return SspParetoFront(); + } + if (storm::utility::isZero(factor)) { + return *this; + } + if (storm::utility::isOne(factor)) { + return minkowskiSum(other); + } + if (isSingleton()) { + return other.scaledTranslated(factor, points.front()); + } + if (other.isSingleton()) { + Point const& point = other.points.front(); + return translated(Point{factor * point.probability, factor * point.expectedCost}); + } + if (storm::utility::isPositive(factor)) { + return minkowskiSumConvexChain(other, factor); + } + + container_type sumPoints; + sumPoints.reserve(points.size() * other.points.size()); + for (auto const& left : points) { + for (auto const& right : other.points) { + sumPoints.push_back(Point{left.probability + factor * right.probability, left.expectedCost + factor * right.expectedCost}); + } + } + return SspParetoFront(std::move(sumPoints)); + } + + static SspParetoFront convexUnion(std::vector const& fronts) { + SspParetoFront const* onlyNonEmptyFront = nullptr; + std::size_t totalPointCount = 0; + for (auto const& front : fronts) { + if (front.empty()) { + continue; + } + onlyNonEmptyFront = onlyNonEmptyFront == nullptr ? &front : onlyNonEmptyFront; + totalPointCount += front.size(); + } + if (totalPointCount == 0) { + return SspParetoFront(); + } + if (onlyNonEmptyFront != nullptr && onlyNonEmptyFront->size() == totalPointCount) { + return *onlyNonEmptyFront; + } + return SspParetoFront(mergeSortedFrontPoints(fronts, totalPointCount), AlreadySortedTag{}); + } + + static SspParetoFront convexUnion(std::vector&& fronts) { + return convexUnionDestructive(fronts); + } + + static SspParetoFront convexUnionDestructive(std::vector& fronts) { + SspParetoFront* onlyNonEmptyFront = nullptr; + std::size_t totalPointCount = 0; + for (auto& front : fronts) { + if (front.empty()) { + continue; + } + onlyNonEmptyFront = onlyNonEmptyFront == nullptr ? &front : onlyNonEmptyFront; + totalPointCount += front.size(); + } + if (totalPointCount == 0) { + return SspParetoFront(); + } + if (onlyNonEmptyFront != nullptr && onlyNonEmptyFront->size() == totalPointCount) { + return std::move(*onlyNonEmptyFront); + } + return SspParetoFront(mergeSortedFrontPoints(fronts, totalPointCount), AlreadySortedTag{}); + } + + /*! + * Returns the minimal continuation cost E on the lower boundary at the given probability p. + * + * This is the paper-specific query used to evaluate a frontier for a fixed bound n: + * it yields the value min{E | (p, E) is contained in the convex Pareto set}. + */ + std::optional getMinimalContinuationCostAtProbability(ValueType const& probability) const { + if (empty()) { + return std::nullopt; + } + if (probability > points.back().probability) { + return std::nullopt; + } + if (probability <= points.front().probability) { + return points.front().expectedCost; + } + + auto rightIt = + std::lower_bound(points.begin(), points.end(), probability, [](Point const& point, ValueType const& value) { return point.probability < value; }); + STORM_LOG_ASSERT(rightIt != points.begin() && rightIt != points.end(), "Expected probability to lie inside the SSP Pareto-front range."); + + Point const& left = *(rightIt - 1); + Point const& right = *rightIt; + ValueType const probabilityDelta = right.probability - left.probability; + STORM_LOG_ASSERT(probabilityDelta > 0, "Expected SSP Pareto front points to be strictly sorted by probability."); + ValueType const interpolationFactor = (probability - left.probability) / probabilityDelta; + return left.expectedCost + interpolationFactor * (right.expectedCost - left.expectedCost); + } + + std::string toString() const { + std::stringstream stream; + stream << "{"; + bool first = true; + for (auto const& point : points) { + if (!first) { + stream << ", "; + } + first = false; + stream << "(" << point.probability << ", " << point.expectedCost << ")"; + } + stream << "}"; + return stream.str(); + } + + private: + SspParetoFront translated(Point const& offset) const { + if (empty()) { + return SspParetoFront(); + } + if (storm::utility::isZero(offset.probability) && storm::utility::isZero(offset.expectedCost)) { + return *this; + } + if (isSingleton()) { + Point const& point = points.front(); + return singleton(point.probability + offset.probability, point.expectedCost + offset.expectedCost); + } + + container_type translatedPoints; + translatedPoints.reserve(points.size()); + for (auto const& point : points) { + translatedPoints.push_back(Point{point.probability + offset.probability, point.expectedCost + offset.expectedCost}); + } + return SspParetoFront(std::move(translatedPoints), AlreadySortedTag{}); + } + + SspParetoFront scaledTranslated(ValueType const& factor, Point const& offset) const { + if (empty()) { + return SspParetoFront(); + } + if (storm::utility::isZero(factor)) { + return singleton(offset.probability, offset.expectedCost); + } + if (isSingleton()) { + Point const& point = points.front(); + return singleton(offset.probability + factor * point.probability, offset.expectedCost + factor * point.expectedCost); + } + + container_type scaledTranslatedPoints; + scaledTranslatedPoints.reserve(points.size()); + for (auto const& point : points) { + scaledTranslatedPoints.push_back(Point{offset.probability + factor * point.probability, offset.expectedCost + factor * point.expectedCost}); + } + if (storm::utility::isPositive(factor)) { + return SspParetoFront(std::move(scaledTranslatedPoints), AlreadySortedTag{}); + } + return SspParetoFront(std::move(scaledTranslatedPoints)); + } + + SspParetoFront minkowskiSumConvexChain(SspParetoFront const& other, ValueType const& factor) const { + STORM_LOG_ASSERT(storm::utility::isPositive(factor), "Expected a positive scaling factor for convex-chain Minkowski sum."); + + container_type sumPoints; + sumPoints.reserve(points.size() + other.points.size() - 1); + + Point current{points.front().probability + factor * other.points.front().probability, + points.front().expectedCost + factor * other.points.front().expectedCost}; + sumPoints.push_back(current); + + std::size_t leftEdge = 0; + std::size_t rightEdge = 0; + while (leftEdge + 1 < points.size() || rightEdge + 1 < other.points.size()) { + if (rightEdge + 1 == other.points.size()) { + addEdge(current, points[leftEdge], points[leftEdge + 1]); + ++leftEdge; + } else if (leftEdge + 1 == points.size()) { + addScaledEdge(current, other.points[rightEdge], other.points[rightEdge + 1], factor); + ++rightEdge; + } else { + int_fast8_t const slopeComparison = + compareEdgeSlopes(points[leftEdge], points[leftEdge + 1], other.points[rightEdge], other.points[rightEdge + 1]); + if (slopeComparison < 0) { + addEdge(current, points[leftEdge], points[leftEdge + 1]); + ++leftEdge; + } else if (slopeComparison > 0) { + addScaledEdge(current, other.points[rightEdge], other.points[rightEdge + 1], factor); + ++rightEdge; + } else { + addEdge(current, points[leftEdge], points[leftEdge + 1]); + addScaledEdge(current, other.points[rightEdge], other.points[rightEdge + 1], factor); + ++leftEdge; + ++rightEdge; + } + } + sumPoints.push_back(current); + } + + return SspParetoFront(std::move(sumPoints), AlreadySortedTag{}); + } + + void canonicalize() { + if (points.empty()) { + return; + } + sortPoints(); + canonicalizeSorted(); + } + + void canonicalizeSorted() { + if (points.empty()) { + return; + } + removeDuplicateAndDominatedPoints(); + removeNonExtremeConvexPoints(); + } + + void sortPoints() { + std::sort(points.begin(), points.end(), [](Point const& left, Point const& right) { + if (left.probability == right.probability) { + return left.expectedCost < right.expectedCost; + } + return left.probability < right.probability; + }); + } + + static bool pointLess(Point const& left, Point const& right) { + if (left.probability == right.probability) { + return left.expectedCost < right.expectedCost; + } + return left.probability < right.probability; + } + + static void addEdge(Point& point, Point const& edgeStart, Point const& edgeEnd) { + point.probability += edgeEnd.probability - edgeStart.probability; + point.expectedCost += edgeEnd.expectedCost - edgeStart.expectedCost; + } + + static void addScaledEdge(Point& point, Point const& edgeStart, Point const& edgeEnd, ValueType const& factor) { + point.probability += factor * (edgeEnd.probability - edgeStart.probability); + point.expectedCost += factor * (edgeEnd.expectedCost - edgeStart.expectedCost); + } + + static int_fast8_t compareEdgeSlopes(Point const& leftStart, Point const& leftEnd, Point const& rightStart, Point const& rightEnd) { + ValueType const leftProbabilityDelta = leftEnd.probability - leftStart.probability; + ValueType const rightProbabilityDelta = rightEnd.probability - rightStart.probability; + STORM_LOG_ASSERT(leftProbabilityDelta > 0 && rightProbabilityDelta > 0, "Expected strictly increasing probabilities in SSP Pareto front edges."); + + ValueType const leftScaledCostDelta = (leftEnd.expectedCost - leftStart.expectedCost) * rightProbabilityDelta; + ValueType const rightScaledCostDelta = (rightEnd.expectedCost - rightStart.expectedCost) * leftProbabilityDelta; + if (leftScaledCostDelta < rightScaledCostDelta) { + return -1; + } + if (rightScaledCostDelta < leftScaledCostDelta) { + return 1; + } + return 0; + } + + static container_type mergeSortedFrontPoints(std::vector const& fronts, std::size_t totalPointCount) { + struct MergeCursor { + const_iterator iterator; + const_iterator end; + }; + auto cursorGreater = [](MergeCursor const& left, MergeCursor const& right) { return pointLess(*right.iterator, *left.iterator); }; + std::vector initialCursors; + initialCursors.reserve(fronts.size()); + for (auto const& front : fronts) { + if (!front.empty()) { + initialCursors.push_back(MergeCursor{front.begin(), front.end()}); + } + } + std::priority_queue, decltype(cursorGreater)> cursors(cursorGreater, std::move(initialCursors)); + + container_type mergedPoints; + mergedPoints.reserve(totalPointCount); + while (!cursors.empty()) { + MergeCursor cursor = cursors.top(); + cursors.pop(); + mergedPoints.push_back(*cursor.iterator); + ++cursor.iterator; + if (cursor.iterator != cursor.end) { + cursors.push(cursor); + } + } + STORM_LOG_ASSERT(mergedPoints.size() == totalPointCount, "Unexpected number of points produced by sorted SSP Pareto-front merge."); + return mergedPoints; + } + + void removeDuplicateAndDominatedPoints() { + if (points.size() < 2) { + return; + } + + if constexpr (FrontKind == SspParetoFrontKind::CostUpperTail) { + std::size_t writeIndex = points.size(); + std::size_t index = points.size(); + bool hasBestExpectedCostSeenFromRight = false; + ValueType bestExpectedCostSeenFromRight{}; + while (index > 0) { + std::size_t const groupEnd = index; + ValueType const probability = points[groupEnd - 1].probability; + while (index > 0 && points[index - 1].probability == probability) { + --index; + } + + Point const& bestPointForProbability = points[index]; + if (!hasBestExpectedCostSeenFromRight || bestPointForProbability.expectedCost < bestExpectedCostSeenFromRight) { + --writeIndex; + points[writeIndex] = bestPointForProbability; + bestExpectedCostSeenFromRight = bestPointForProbability.expectedCost; + hasBestExpectedCostSeenFromRight = true; + } + } + points.erase(points.begin(), points.begin() + static_cast(writeIndex)); + } else { + std::size_t writeIndex = 0; + std::size_t index = 0; + bool hasBestExpectedCostSeenFromLeft = false; + ValueType bestExpectedCostSeenFromLeft{}; + while (index < points.size()) { + ValueType const probability = points[index].probability; + Point const& bestPointForProbability = points[index]; + while (index < points.size() && points[index].probability == probability) { + ++index; + } + + if (!hasBestExpectedCostSeenFromLeft || bestPointForProbability.expectedCost < bestExpectedCostSeenFromLeft) { + points[writeIndex] = bestPointForProbability; + ++writeIndex; + bestExpectedCostSeenFromLeft = bestPointForProbability.expectedCost; + hasBestExpectedCostSeenFromLeft = true; + } + } + points.resize(writeIndex); + } + } + + void removeNonExtremeConvexPoints() { + if (points.size() < 3) { + return; + } + std::size_t hullSize = 0; + for (std::size_t index = 0, endIndex = points.size(); index < endIndex; ++index) { + points[hullSize++] = points[index]; + while (hullSize >= 3 && liesOnOrAboveSegment(points[hullSize - 3], points[hullSize - 2], points[hullSize - 1])) { + points[hullSize - 2] = points[hullSize - 1]; + --hullSize; + } + } + points.resize(hullSize); + if constexpr (FrontKind == SspParetoFrontKind::CostUpperTail) { + STORM_LOG_ASSERT(std::adjacent_find(points.begin(), points.end(), + [](Point const& left, Point const& right) { + return left.probability >= right.probability || left.expectedCost >= right.expectedCost; + }) == points.end(), + "Expected cost SSP Pareto front points to be strictly ordered by increasing probability and increasing expected cost."); + } else { + STORM_LOG_ASSERT(std::adjacent_find(points.begin(), points.end(), + [](Point const& left, Point const& right) { + return left.probability >= right.probability || left.expectedCost <= right.expectedCost; + }) == points.end(), + "Expected reward SSP Pareto front points to be strictly ordered by increasing probability and decreasing expected shortfall."); + } + } + + static bool liesOnOrAboveSegment(Point const& left, Point const& middle, Point const& right) { + ValueType const leftToRightProbabilityDelta = right.probability - left.probability; + ValueType const leftToMiddleProbabilityDelta = middle.probability - left.probability; + STORM_LOG_ASSERT(leftToRightProbabilityDelta > 0 && leftToMiddleProbabilityDelta > 0, + "Expected SSP Pareto front probabilities to be strictly increasing."); + return (middle.expectedCost - left.expectedCost) * leftToRightProbabilityDelta >= + (right.expectedCost - left.expectedCost) * leftToMiddleProbabilityDelta; + } + + container_type points; +}; + +template +using SspRewardParetoFront = SspParetoFront; + +} // namespace cvar +} // namespace modelchecker +} // namespace storm diff --git a/src/storm/modelchecker/cvar/helper/SspParetoValueIterationOperator.h b/src/storm/modelchecker/cvar/helper/SspParetoValueIterationOperator.h new file mode 100644 index 0000000000..307f4583ed --- /dev/null +++ b/src/storm/modelchecker/cvar/helper/SspParetoValueIterationOperator.h @@ -0,0 +1,210 @@ +#pragma once + +#include +#include +#include +#include + +#include "storm/modelchecker/cvar/helper/SspParetoFront.h" +#include "storm/modelchecker/cvar/preprocessing/SspCvarPreprocessingResult.h" +#include "storm/utility/constants.h" +#include "storm/utility/macros.h" + +namespace storm { +namespace modelchecker { +namespace cvar { + +/*! + * Applies one cost-bound layer update for SSP CVaR Pareto-front value iteration. + * + * This class follows the role of Storm's generic value-iteration operators, but keeps the SSP-specific + * cost-window semantics explicit: each action row reads from the predecessor layer determined by the current + * cost bound minus that action's integer cost. + */ +template +class SspParetoValueIterationOperator { + public: + using ParetoFront = SspParetoFront; + using FrontierLayer = std::vector; + using FrontierWindow = std::vector; + + explicit SspParetoValueIterationOperator(preprocessing::SspCvarPreprocessingResult const& preprocessingResult) + : preprocessingResult(preprocessingResult), + choiceCostOffsets(createChoiceCostOffsets(preprocessingResult)), + reachableTargetStates(collectReachableStates(preprocessingResult, true)), + reachableNonTargetStates(collectReachableStates(preprocessingResult, false)) {} + + void apply(uint64_t costBound, FrontierWindow const& frontierWindow, FrontierLayer& outputLayer) const { + prepareOutputLayer(outputLayer); + + ParetoFront const targetFront = createTargetFrontier(costBound); + for (auto state : reachableTargetStates) { + outputLayer[state] = targetFront; + } + + std::vector actionFrontStorage; + ActionFrontReducer actionFrontReducer(actionFrontStorage); + for (auto state : reachableNonTargetStates) { + applyRowGroup(state, costBound, frontierWindow, outputLayer, actionFrontReducer); + } + } + + std::size_t getStateCount() const { + return preprocessingResult.transitionMatrix.getRowGroupCount(); + } + + std::vector const& getReachableTargetStates() const { + return reachableTargetStates; + } + + std::vector const& getReachableNonTargetStates() const { + return reachableNonTargetStates; + } + + static std::size_t getWindowIndex(int64_t costBound, std::size_t windowSize) { + STORM_LOG_ASSERT(windowSize > 0, "Expected a non-empty SSP frontier window."); + int64_t const signedWindowSize = static_cast(windowSize); + int64_t index = costBound % signedWindowSize; + if (index < 0) { + index += signedWindowSize; + } + return static_cast(index); + } + + static ParetoFront createTargetFrontier() { + return createTargetFrontier(0); + } + + static ParetoFront createTargetFrontier(uint64_t costBound) { + if constexpr (FrontKind == SspParetoFrontKind::CostUpperTail) { + return ParetoFront::singleton(storm::utility::one(), storm::utility::zero()); + } else { + return ParetoFront::singleton(storm::utility::one(), storm::utility::convertNumber(costBound)); + } + } + + private: + class ActionFrontReducer { + public: + explicit ActionFrontReducer(std::vector& actionFrontStorage) : actionFrontStorage(actionFrontStorage) {} + + void reset(std::size_t rowCount) { + actionFrontStorage.clear(); + rowCountHint = rowCount; + firstNonEmptyActionFront.clear(); + hasNonEmptyActionFront = false; + usesActionFrontStorage = false; + } + + void add(ParetoFront&& actionFront) { + addActionFront(std::move(actionFront)); + } + + void reduceInto(ParetoFront& outputFront) { + if (!hasNonEmptyActionFront) { + outputFront.clear(); + } else if (!usesActionFrontStorage) { + outputFront = std::move(firstNonEmptyActionFront); + } else { + outputFront = ParetoFront::convexUnionDestructive(actionFrontStorage); + } + } + + private: + void addActionFront(ParetoFront&& actionFront) { + if (actionFront.empty()) { + return; + } + if (!hasNonEmptyActionFront) { + firstNonEmptyActionFront = std::move(actionFront); + hasNonEmptyActionFront = true; + return; + } + if (!usesActionFrontStorage) { + actionFrontStorage.reserve(rowCountHint); + actionFrontStorage.push_back(std::move(firstNonEmptyActionFront)); + usesActionFrontStorage = true; + } + actionFrontStorage.push_back(std::move(actionFront)); + } + + std::vector& actionFrontStorage; + ParetoFront firstNonEmptyActionFront; + std::size_t rowCountHint{0}; + bool hasNonEmptyActionFront{false}; + bool usesActionFrontStorage{false}; + }; + + void prepareOutputLayer(FrontierLayer& layer) const { + if (layer.size() != getStateCount()) { + layer.resize(getStateCount()); + } + for (auto& front : layer) { + front.clear(); + } + } + + void applyRowGroup(uint64_t state, uint64_t costBound, FrontierWindow const& frontierWindow, FrontierLayer& outputLayer, + ActionFrontReducer& actionFrontReducer) const { + uint64_t const firstActionRow = preprocessingResult.transitionMatrix.getRowGroupIndices()[state]; + uint64_t const endActionRow = preprocessingResult.transitionMatrix.getRowGroupIndices()[state + 1]; + STORM_LOG_ASSERT(firstActionRow < endActionRow, "Expected at least one action row in each reachable non-target SSP state."); + + actionFrontReducer.reset(endActionRow - firstActionRow); + for (uint64_t actionRow = firstActionRow; actionRow < endActionRow; ++actionRow) { + actionFrontReducer.add(computeActionFront(actionRow, costBound, frontierWindow)); + } + actionFrontReducer.reduceInto(outputLayer[state]); + } + + ParetoFront computeActionFront(uint64_t actionRow, uint64_t costBound, FrontierWindow const& frontierWindow) const { + uint64_t const actionCost = choiceCostOffsets[actionRow]; + ParetoFront actionFront; + int64_t const predecessorBound = static_cast(costBound) - static_cast(actionCost); + FrontierLayer const& predecessorLayer = getFrontierLayerForBound(predecessorBound, frontierWindow); + + bool initialized = false; + for (auto const& transition : preprocessingResult.transitionMatrix.getRow(actionRow)) { + auto const& successorFront = predecessorLayer[transition.getColumn()]; + if (initialized) { + actionFront = actionFront.minkowskiSumScaled(successorFront, transition.getValue()); + } else { + actionFront = successorFront.scaled(transition.getValue()); + initialized = true; + } + } + return actionFront; + } + + FrontierLayer const& getFrontierLayerForBound(int64_t costBound, FrontierWindow const& frontierWindow) const { + return frontierWindow[getWindowIndex(costBound, frontierWindow.size())]; + } + + static std::vector createChoiceCostOffsets(preprocessing::SspCvarPreprocessingResult const& preprocessingResult) { + std::vector result; + result.reserve(preprocessingResult.choiceCosts.size()); + for (auto const& choiceCost : preprocessingResult.choiceCosts) { + result.push_back(storm::utility::convertNumber(choiceCost)); + } + return result; + } + + static std::vector collectReachableStates(preprocessing::SspCvarPreprocessingResult const& preprocessingResult, bool targetStates) { + std::vector result; + for (uint64_t state = 0; state < preprocessingResult.transitionMatrix.getRowGroupCount(); ++state) { + if (preprocessingResult.reachableStates[state] && preprocessingResult.targetStates[state] == targetStates) { + result.push_back(state); + } + } + return result; + } + + preprocessing::SspCvarPreprocessingResult const& preprocessingResult; + std::vector choiceCostOffsets; + std::vector reachableTargetStates; + std::vector reachableNonTargetStates; +}; + +} // namespace cvar +} // namespace modelchecker +} // namespace storm diff --git a/src/storm/modelchecker/cvar/preprocessing/SspCvarPreprocessingResult.h b/src/storm/modelchecker/cvar/preprocessing/SspCvarPreprocessingResult.h new file mode 100644 index 0000000000..e586798dc3 --- /dev/null +++ b/src/storm/modelchecker/cvar/preprocessing/SspCvarPreprocessingResult.h @@ -0,0 +1,31 @@ +#pragma once + +#include +#include + +#include "storm/storage/BitVector.h" +#include "storm/storage/SparseMatrix.h" + +namespace storm { +namespace modelchecker { +namespace cvar { +namespace preprocessing { + +template +struct SspCvarPreprocessingResult { + std::string rewardModelName; + uint64_t initialState; + storm::storage::BitVector targetStates; + storm::storage::BitVector reachableStates; + bool liftedStateRewardsToChoiceCosts; + bool normalizedTargetStatesToAbsorbing; + uint64_t maximalChoiceCost; + std::vector choiceCosts; + std::vector expectedCostsToGoal; + storm::storage::SparseMatrix transitionMatrix; +}; + +} // namespace preprocessing +} // namespace cvar +} // namespace modelchecker +} // namespace storm diff --git a/src/storm/modelchecker/cvar/preprocessing/SspCvarPreprocessor.h b/src/storm/modelchecker/cvar/preprocessing/SspCvarPreprocessor.h new file mode 100644 index 0000000000..75fc6f5d56 --- /dev/null +++ b/src/storm/modelchecker/cvar/preprocessing/SspCvarPreprocessor.h @@ -0,0 +1,249 @@ +#pragma once + +#include +#include +#include + +#include "storm/exceptions/InvalidPropertyException.h" +#include "storm/exceptions/NotImplementedException.h" +#include "storm/modelchecker/cvar/CvarQueryInformation.h" +#include "storm/modelchecker/cvar/preprocessing/SspCvarPreprocessingResult.h" +#include "storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.h" +#include "storm/models/sparse/StandardRewardModel.h" +#include "storm/solver/SolveGoal.h" +#include "storm/storage/BitVector.h" +#include "storm/storage/SparseMatrix.h" +#include "storm/utility/constants.h" +#include "storm/utility/graph.h" +#include "storm/utility/logging.h" +#include "storm/utility/macros.h" + +namespace storm { +namespace modelchecker { +namespace cvar { +namespace preprocessing { + +/*! + * Preprocesses a sparse MDP for the SSP CVaR backend. + * + * This normalizes the model to terminal-goal semantics and extracts the + * choice-based costs and graph information needed by the future Pareto-front VI. + */ +template +std::vector extractChoiceCostsForSsp(SparseMdpModelType const& model, + typename SparseMdpModelType::RewardModelType const& rewardModel, + storm::storage::BitVector const& targetStates) { + using ValueType = typename SparseMdpModelType::ValueType; + + std::vector choiceCosts(model.getNumberOfChoices(), storm::utility::zero()); + bool hasStateRewards = rewardModel.hasStateRewards(); + bool hasStateActionRewards = rewardModel.hasStateActionRewards(); + + for (uint64_t state = 0; state < model.getNumberOfStates(); ++state) { + if (targetStates[state]) { + continue; + } + + ValueType stateReward = hasStateRewards ? rewardModel.getStateReward(state) : storm::utility::zero(); + for (uint64_t row = model.getTransitionMatrix().getRowGroupIndices()[state], endRow = model.getTransitionMatrix().getRowGroupIndices()[state + 1]; + row < endRow; ++row) { + choiceCosts[row] = stateReward; + if (hasStateActionRewards) { + choiceCosts[row] += rewardModel.getStateActionReward(row); + } + } + } + + return choiceCosts; +} + +template +void validatePositiveChoiceCostsOutsideGoals(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& targetStates, + std::vector const& choiceCosts, std::string const& valueName = "choice costs") { + ValueType const zero = storm::utility::zero(); + for (uint64_t state = 0; state < transitionMatrix.getRowGroupCount(); ++state) { + if (targetStates[state]) { + continue; + } + for (uint64_t row = transitionMatrix.getRowGroupIndices()[state], endRow = transitionMatrix.getRowGroupIndices()[state + 1]; row < endRow; ++row) { + STORM_LOG_THROW(choiceCosts[row] > zero, storm::exceptions::InvalidPropertyException, + "CVaR SSP preprocessing currently requires strictly positive " << valueName << " outside goal states."); + } + } +} + +template +uint64_t validateAndComputeMaximalChoiceCostOutsideGoals(storm::storage::SparseMatrix const& transitionMatrix, + storm::storage::BitVector const& targetStates, std::vector const& choiceCosts, + std::string const& valueName = "choice costs") { + uint64_t maximalChoiceCost = 0; + for (uint64_t state = 0; state < transitionMatrix.getRowGroupCount(); ++state) { + if (targetStates[state]) { + continue; + } + for (uint64_t row = transitionMatrix.getRowGroupIndices()[state], endRow = transitionMatrix.getRowGroupIndices()[state + 1]; row < endRow; ++row) { + STORM_LOG_THROW(storm::utility::isInteger(choiceCosts[row]), storm::exceptions::InvalidPropertyException, + "CVaR SSP preprocessing currently requires integer-valued " << valueName << "."); + maximalChoiceCost = std::max(maximalChoiceCost, storm::utility::convertNumber(choiceCosts[row])); + } + } + return maximalChoiceCost; +} + +template +std::vector computeExpectedCostsToGoal(Environment const& env, storm::storage::SparseMatrix const& transitionMatrix, + storm::storage::SparseMatrix const& backwardTransitions, + storm::storage::BitVector const& targetStates, std::vector const& choiceCosts) { + storm::models::sparse::StandardRewardModel rewardModel(std::nullopt, std::vector(choiceCosts), std::nullopt); + auto result = storm::modelchecker::helper::SparseMdpPrctlHelper::computeReachabilityRewards( + env, storm::solver::SolveGoal(storm::OptimizationDirection::Minimize), transitionMatrix, backwardTransitions, rewardModel, + targetStates, false, false); + return std::move(result.values); +} + +template +bool normalizeTargetStatesToAbsorbing(storm::storage::SparseMatrix& transitionMatrix, storm::storage::BitVector const& targetStates) { + bool normalizedTargetStatesToAbsorbing = false; + for (auto targetState : targetStates) { + for (uint64_t row = transitionMatrix.getRowGroupIndices()[targetState], endRow = transitionMatrix.getRowGroupIndices()[targetState + 1]; row < endRow; + ++row) { + for (auto const& entry : transitionMatrix.getRow(row)) { + if (entry.getColumn() != targetState) { + normalizedTargetStatesToAbsorbing = true; + break; + } + } + if (normalizedTargetStatesToAbsorbing) { + break; + } + } + if (normalizedTargetStatesToAbsorbing) { + break; + } + } + if (normalizedTargetStatesToAbsorbing) { + STORM_LOG_INFO("CVaR SSP preprocessing makes target states absorbing to match terminal-goal semantics."); + transitionMatrix.makeRowGroupsAbsorbing(targetStates, true); + } + return normalizedTargetStatesToAbsorbing; +} + +template +void validateRewardAlmostSureReachability(storm::storage::SparseMatrix const& transitionMatrix, + storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& targetStates, + storm::storage::BitVector const& reachableStates) { + storm::storage::BitVector prob1AStates = + storm::utility::graph::performProb1A(transitionMatrix, transitionMatrix.getRowGroupIndices(), backwardTransitions, + storm::storage::BitVector(transitionMatrix.getRowGroupCount(), true), targetStates); + STORM_LOG_THROW(reachableStates.isSubsetOf(prob1AStates), storm::exceptions::InvalidPropertyException, + "CVaR SSP reward preprocessing currently requires all reachable states to reach a goal almost surely under all schedulers."); +} + +template +SspCvarPreprocessingResult preprocessSspCvar(Environment const& env, SparseMdpModelType const& model, + CvarQueryInformation const& queryInformation, + storm::storage::BitVector const& targetStates) { + using ValueType = typename SparseMdpModelType::ValueType; + + std::string rewardModelName = queryInformation.rewardModelName.value_or(""); + auto const& rewardModel = model.getRewardModel(rewardModelName); + if (rewardModelName.empty()) { + rewardModelName = model.getUniqueRewardModelName(); + } + + STORM_LOG_THROW(queryInformation.optimizationDirection == storm::solver::OptimizationDirection::Minimize, storm::exceptions::InvalidPropertyException, + "CVaR SSP preprocessing currently only supports minimizing total costs."); + STORM_LOG_THROW(queryInformation.interpretation == CvarInterpretation::Cost, storm::exceptions::InvalidPropertyException, + "CVaR SSP preprocessing currently only supports the cost interpretation."); + + STORM_LOG_THROW(!rewardModel.hasTransitionRewards(), storm::exceptions::NotImplementedException, + "CVaR SSP preprocessing does not support transition rewards."); + + bool liftedStateRewardsToChoiceCosts = rewardModel.hasStateRewards() && !rewardModel.hasStateActionRewards(); + if (liftedStateRewardsToChoiceCosts) { + STORM_LOG_INFO("CVaR SSP preprocessing lifts state rewards to equivalent per-choice costs."); + } + + auto transitionMatrix = model.getTransitionMatrix(); + bool normalizedTargetStatesToAbsorbing = normalizeTargetStatesToAbsorbing(transitionMatrix, targetStates); + + auto backwardTransitions = transitionMatrix.transpose(true); + auto reachableStates = storm::utility::graph::getReachableStates(transitionMatrix, model.getInitialStates(), + storm::storage::BitVector(transitionMatrix.getRowGroupCount(), true), + storm::storage::BitVector(transitionMatrix.getRowGroupCount(), false)); + auto properStates = storm::utility::graph::performProb1E(transitionMatrix, transitionMatrix.getRowGroupIndices(), backwardTransitions, + storm::storage::BitVector(transitionMatrix.getRowGroupCount(), true), targetStates); + STORM_LOG_THROW(reachableStates.isSubsetOf(properStates), storm::exceptions::InvalidPropertyException, + "CVaR SSP preprocessing currently requires a proper policy from every reachable state."); + auto choiceCosts = extractChoiceCostsForSsp(model, rewardModel, targetStates); + validatePositiveChoiceCostsOutsideGoals(transitionMatrix, targetStates, choiceCosts); + uint64_t maximalChoiceCost = validateAndComputeMaximalChoiceCostOutsideGoals(transitionMatrix, targetStates, choiceCosts); + auto expectedCostsToGoal = computeExpectedCostsToGoal(env, transitionMatrix, backwardTransitions, targetStates, choiceCosts); + + return {rewardModelName, + *model.getInitialStates().begin(), + targetStates, + std::move(reachableStates), + liftedStateRewardsToChoiceCosts, + normalizedTargetStatesToAbsorbing, + maximalChoiceCost, + std::move(choiceCosts), + std::move(expectedCostsToGoal), + std::move(transitionMatrix)}; +} + +template +SspCvarPreprocessingResult preprocessSspRewardCvar(Environment const&, SparseMdpModelType const& model, + CvarQueryInformation const& queryInformation, + storm::storage::BitVector const& targetStates) { + using ValueType = typename SparseMdpModelType::ValueType; + + std::string rewardModelName = queryInformation.rewardModelName.value_or(""); + auto const& rewardModel = model.getRewardModel(rewardModelName); + if (rewardModelName.empty()) { + rewardModelName = model.getUniqueRewardModelName(); + } + + STORM_LOG_THROW(queryInformation.optimizationDirection == storm::solver::OptimizationDirection::Maximize, storm::exceptions::InvalidPropertyException, + "CVaR SSP reward preprocessing currently only supports maximizing total rewards."); + STORM_LOG_THROW(queryInformation.interpretation == CvarInterpretation::Reward, storm::exceptions::InvalidPropertyException, + "CVaR SSP reward preprocessing requires the reward interpretation."); + + STORM_LOG_THROW(!rewardModel.hasTransitionRewards(), storm::exceptions::NotImplementedException, + "CVaR SSP preprocessing does not support transition rewards."); + + bool liftedStateRewardsToChoiceCosts = rewardModel.hasStateRewards() && !rewardModel.hasStateActionRewards(); + if (liftedStateRewardsToChoiceCosts) { + STORM_LOG_INFO("CVaR SSP preprocessing lifts state rewards to equivalent per-choice rewards."); + } + + auto transitionMatrix = model.getTransitionMatrix(); + bool normalizedTargetStatesToAbsorbing = normalizeTargetStatesToAbsorbing(transitionMatrix, targetStates); + auto backwardTransitions = transitionMatrix.transpose(true); + auto reachableStates = storm::utility::graph::getReachableStates(transitionMatrix, model.getInitialStates(), + storm::storage::BitVector(transitionMatrix.getRowGroupCount(), true), + storm::storage::BitVector(transitionMatrix.getRowGroupCount(), false)); + + uint64_t const initialState = *model.getInitialStates().begin(); + validateRewardAlmostSureReachability(transitionMatrix, backwardTransitions, targetStates, reachableStates); + + auto choiceRewards = extractChoiceCostsForSsp(model, rewardModel, targetStates); + validatePositiveChoiceCostsOutsideGoals(transitionMatrix, targetStates, choiceRewards, "choice rewards"); + uint64_t maximalChoiceReward = validateAndComputeMaximalChoiceCostOutsideGoals(transitionMatrix, targetStates, choiceRewards, "choice rewards"); + + return {rewardModelName, + initialState, + targetStates, + std::move(reachableStates), + liftedStateRewardsToChoiceCosts, + normalizedTargetStatesToAbsorbing, + maximalChoiceReward, + std::move(choiceRewards), + std::vector(), + std::move(transitionMatrix)}; +} + +} // namespace preprocessing +} // namespace cvar +} // namespace modelchecker +} // namespace storm diff --git a/src/storm/modelchecker/cvar/preprocessing/WeightedReachabilityCvarPreprocessingResult.h b/src/storm/modelchecker/cvar/preprocessing/WeightedReachabilityCvarPreprocessingResult.h new file mode 100644 index 0000000000..0a63d03762 --- /dev/null +++ b/src/storm/modelchecker/cvar/preprocessing/WeightedReachabilityCvarPreprocessingResult.h @@ -0,0 +1,29 @@ +#pragma once + +#include +#include + +#include "storm/storage/BitVector.h" +#include "storm/storage/SparseMatrix.h" + +namespace storm { +namespace modelchecker { +namespace cvar { +namespace preprocessing { + +template +struct WeightedReachabilityCvarPreprocessingResult { + std::string rewardModelName; + uint64_t initialState; + storm::storage::BitVector originalTargetStates; + storm::storage::BitVector effectiveTargetStates; + storm::storage::BitVector badMecStates; + uint64_t collapsedTargetReachingMecCount; + std::vector terminalRewards; + storm::storage::SparseMatrix transitionMatrix; +}; + +} // namespace preprocessing +} // namespace cvar +} // namespace modelchecker +} // namespace storm diff --git a/src/storm/modelchecker/cvar/preprocessing/WeightedReachabilityCvarPreprocessor.h b/src/storm/modelchecker/cvar/preprocessing/WeightedReachabilityCvarPreprocessor.h new file mode 100644 index 0000000000..fcb9181b61 --- /dev/null +++ b/src/storm/modelchecker/cvar/preprocessing/WeightedReachabilityCvarPreprocessor.h @@ -0,0 +1,157 @@ +#pragma once + +#include "storm/exceptions/InvalidOperationException.h" +#include "storm/exceptions/InvalidPropertyException.h" +#include "storm/modelchecker/cvar/CvarPreprocessingUtilities.h" +#include "storm/modelchecker/cvar/CvarQueryInformation.h" +#include "storm/modelchecker/cvar/preprocessing/WeightedReachabilityCvarPreprocessingResult.h" +#include "storm/storage/MaximalEndComponentDecomposition.h" +#include "storm/transformer/EndComponentEliminator.h" +#include "storm/utility/constants.h" +#include "storm/utility/graph.h" +#include "storm/utility/logging.h" +#include "storm/utility/macros.h" + +namespace storm { +namespace modelchecker { +namespace cvar { +namespace preprocessing { + +/*! + * Preprocesses a sparse MDP for the weighted-reachability CVaR LP backend. + * + * This enforces the weighted-reachability assumptions on the reward model and + * transition structure and returns the normalized terminal-reward instance that + * can be consumed by the LP helper. + */ +template +storm::storage::MaximalEndComponentDecomposition computeTargetReachingMecs(storm::storage::SparseMatrix const& transitionMatrix, + storm::storage::BitVector const& initialStates, + storm::storage::BitVector const& effectiveTargetStates) { + storm::storage::BitVector allStates(transitionMatrix.getRowGroupCount(), true); + storm::storage::BitVector nonTargetStates = ~effectiveTargetStates; + storm::storage::BitVector noStates(allStates.size(), false); + auto reachableStates = storm::utility::graph::getReachableStates(transitionMatrix, initialStates, allStates, noStates); + auto subsystemStates = reachableStates & nonTargetStates; + + storm::storage::BitVector possibleEcRows(transitionMatrix.getRowCount(), false); + for (auto state : subsystemStates) { + for (uint64_t row = transitionMatrix.getRowGroupIndices()[state], endRow = transitionMatrix.getRowGroupIndices()[state + 1]; row < endRow; ++row) { + possibleEcRows.set(row, true); + } + } + + auto backwardTransitions = transitionMatrix.transpose(true); + auto statesThatCanReachTarget = storm::utility::graph::performProbGreater0(backwardTransitions, allStates, effectiveTargetStates); + auto targetReachingStates = statesThatCanReachTarget & subsystemStates; + return storm::storage::MaximalEndComponentDecomposition(transitionMatrix, backwardTransitions, targetReachingStates, possibleEcRows); +} + +template +void applyTargetReachingMecCollapse(storm::storage::SparseMatrix& transitionMatrix, storm::storage::BitVector& effectiveTargetStates, + storm::storage::BitVector& badMecStates, std::vector& terminalRewards, uint64_t& initialState, + storm::storage::MaximalEndComponentDecomposition const& targetReachingMecs) { + storm::storage::BitVector allStates(transitionMatrix.getRowGroupCount(), true); + storm::storage::BitVector noSinkRows(transitionMatrix.getRowGroupCount(), false); + auto eliminationResult = storm::transformer::EndComponentEliminator::transform(transitionMatrix, targetReachingMecs, allStates, noSinkRows); + + storm::storage::BitVector newEffectiveTargetStates(eliminationResult.matrix.getRowGroupCount(), false); + storm::storage::BitVector newBadMecStates(eliminationResult.matrix.getRowGroupCount(), false); + std::vector newTerminalRewards(eliminationResult.matrix.getRowGroupCount(), storm::utility::zero()); + for (auto oldTargetState : effectiveTargetStates) { + auto newTargetState = eliminationResult.oldToNewStateMapping[oldTargetState]; + if (newTargetState < eliminationResult.matrix.getRowGroupCount()) { + newEffectiveTargetStates.set(newTargetState, true); + newTerminalRewards[newTargetState] = terminalRewards[oldTargetState]; + } + } + for (auto oldBadMecState : badMecStates) { + auto newBadMecState = eliminationResult.oldToNewStateMapping[oldBadMecState]; + if (newBadMecState < eliminationResult.matrix.getRowGroupCount()) { + newBadMecStates.set(newBadMecState, true); + } + } + + initialState = eliminationResult.oldToNewStateMapping[initialState]; + effectiveTargetStates = std::move(newEffectiveTargetStates); + badMecStates = std::move(newBadMecStates); + terminalRewards = std::move(newTerminalRewards); + transitionMatrix = std::move(eliminationResult.matrix); +} + +template +WeightedReachabilityCvarPreprocessingResult preprocessWeightedReachabilityCvar( + SparseMdpModelType const& model, CvarQueryInformation const& queryInformation, storm::storage::BitVector const& targetStates, + bool produceScheduler = false) { + using ValueType = typename SparseMdpModelType::ValueType; + + std::string rewardModelName = queryInformation.rewardModelName.value_or(""); + auto const& rewardModel = model.getRewardModel(rewardModelName); + if (rewardModelName.empty()) { + rewardModelName = model.getUniqueRewardModelName(); + } + + STORM_LOG_THROW(rewardModel.hasOnlyStateRewards(), storm::exceptions::InvalidPropertyException, + "CVaR queries currently only support weighted reachability with state-based terminal rewards."); + + std::vector const& stateRewards = rewardModel.getStateRewardVector(); + validateTargetStatesAreAbsorbing(model.getTransitionMatrix(), targetStates); + + bool hasNonZeroTargetReward = false; + for (uint64_t state = 0; state < stateRewards.size(); ++state) { + if (targetStates[state]) { + hasNonZeroTargetReward |= !storm::utility::isZero(stateRewards[state]); + } else { + STORM_LOG_THROW(storm::utility::isZero(stateRewards[state]), storm::exceptions::InvalidPropertyException, + "CVaR queries currently require terminal rewards, i.e. non-target states must have reward 0."); + } + } + + if (!hasNonZeroTargetReward) { + STORM_LOG_WARN("All target states have terminal reward 0 in reward model '" << rewardModelName << "'."); + } + + auto reachableMecs = computeReachableMecs(model.getTransitionMatrix(), model.getInitialStates()); + auto statesThatCanReachTarget = computeStatesThatCanReachTarget(model.getTransitionMatrix(), targetStates); + auto badMecStates = computeBadMecStates(reachableMecs, statesThatCanReachTarget, model.getNumberOfStates()); + auto effectiveTargetStates = targetStates | badMecStates; + auto terminalRewards = stateRewards; + auto transitionMatrix = model.getTransitionMatrix(); + uint64_t initialState = *model.getInitialStates().begin(); + + if (!badMecStates.empty()) { + STORM_LOG_INFO( + "CVaR preprocessing converted reachable end components that cannot reach the original target set into zero-reward absorbing terminal behaviour."); + transitionMatrix.makeRowGroupsAbsorbing(badMecStates, true); + for (auto state : badMecStates) { + terminalRewards[state] = storm::utility::zero(); + } + } + + uint64_t collapsedTargetReachingMecCount = 0; + auto targetReachingMecs = computeTargetReachingMecs(transitionMatrix, model.getInitialStates(), effectiveTargetStates); + if (!targetReachingMecs.empty()) { + STORM_LOG_THROW(!produceScheduler, storm::exceptions::InvalidOperationException, + "Cannot produce a CVaR scheduler because preprocessing has to collapse target-reaching end components."); + collapsedTargetReachingMecCount = targetReachingMecs.size(); + auto oldStateCount = transitionMatrix.getRowGroupCount(); + applyTargetReachingMecCollapse(transitionMatrix, effectiveTargetStates, badMecStates, terminalRewards, initialState, targetReachingMecs); + STORM_PRINT_AND_LOG("CVaR preprocessing collapsed " << collapsedTargetReachingMecCount + << " target-reaching end component(s), reducing the transition matrix from " << oldStateCount + << " to " << transitionMatrix.getRowGroupCount() << " states.\n"); + } + + return {rewardModelName, + initialState, + targetStates, + effectiveTargetStates, + badMecStates, + collapsedTargetReachingMecCount, + std::move(terminalRewards), + std::move(transitionMatrix)}; +} + +} // namespace preprocessing +} // namespace cvar +} // namespace modelchecker +} // namespace storm diff --git a/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp b/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp index a572bffbb6..f942b2cd0b 100644 --- a/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp +++ b/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp @@ -3,10 +3,12 @@ #include "storm/adapters/IntervalAdapter.h" #include "storm/adapters/RationalFunctionAdapter.h" #include "storm/adapters/RationalNumberAdapter.h" +#include "storm/environment/modelchecker/ModelCheckerEnvironment.h" #include "storm/exceptions/InvalidPropertyException.h" #include "storm/exceptions/NotImplementedException.h" #include "storm/exceptions/NotSupportedException.h" #include "storm/logic/FragmentSpecification.h" +#include "storm/modelchecker/cvar/CvarModelChecking.h" #include "storm/modelchecker/helper/conditional/ConditionalHelper.h" #include "storm/modelchecker/helper/finitehorizon/SparseStepBoundedHorizonHelper.h" #include "storm/modelchecker/helper/infinitehorizon/SparseNondeterministicInfiniteHorizonHelper.h" @@ -91,7 +93,8 @@ bool SparseMdpPrctlModelChecker::canHandleStatic(CheckTask SparseMdpPrctlModelChecker::che } } +template +std::unique_ptr SparseMdpPrctlModelChecker::checkCvarFormula( + Environment const& env, CheckTask const& checkTask) { + if constexpr (storm::IsIntervalType) { + STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "CVaR formulas are not supported for interval models."); + } else { + auto formulaChecker = [&](storm::logic::Formula const& formula) { + return this->check(env, formula)->template asExplicitQualitativeCheckResult().getTruthValuesVector(); + }; + return cvar::performCvarModelChecking(env, this->getModel(), checkTask, formulaChecker); + } +} + template std::unique_ptr SparseMdpPrctlModelChecker::checkQuantileFormula( Environment const& env, CheckTask const& checkTask) { diff --git a/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.h b/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.h index 5a7f10a1fe..987209b8e4 100644 --- a/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.h +++ b/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.h @@ -62,6 +62,7 @@ class SparseMdpPrctlModelChecker : public SparsePropositionalModelChecker const& checkTask) override; virtual std::unique_ptr checkLexObjectiveFormula(Environment const& env, CheckTask const& checkTask) override; + virtual std::unique_ptr checkCvarFormula(Environment const& env, CheckTask const& checkTask) override; virtual std::unique_ptr checkQuantileFormula(Environment const& env, CheckTask const& checkTask) override; }; diff --git a/src/storm/settings/SettingsManager.cpp b/src/storm/settings/SettingsManager.cpp index 366aac3a59..d944758fca 100644 --- a/src/storm/settings/SettingsManager.cpp +++ b/src/storm/settings/SettingsManager.cpp @@ -19,6 +19,7 @@ #include "storm/settings/modules/ConditionalSettings.h" #include "storm/settings/modules/CoreSettings.h" #include "storm/settings/modules/CuddSettings.h" +#include "storm/settings/modules/CvarSettings.h" #include "storm/settings/modules/DebugSettings.h" #include "storm/settings/modules/EigenEquationSolverSettings.h" #include "storm/settings/modules/EliminationSettings.h" @@ -699,6 +700,7 @@ void initializeAll(std::string const& name, std::string const& executableName) { // Register all known settings modules. storm::settings::addModule(); storm::settings::addModule(); + storm::settings::addModule(); storm::settings::addModule(); storm::settings::addModule(); storm::settings::addModule(); diff --git a/src/storm/settings/modules/CvarSettings.cpp b/src/storm/settings/modules/CvarSettings.cpp new file mode 100644 index 0000000000..3e15f02d92 --- /dev/null +++ b/src/storm/settings/modules/CvarSettings.cpp @@ -0,0 +1,63 @@ +#include "storm/settings/modules/CvarSettings.h" + +#include + +#include "storm/exceptions/IllegalArgumentValueException.h" +#include "storm/settings/ArgumentBuilder.h" +#include "storm/settings/OptionBuilder.h" +#include "storm/utility/macros.h" + +namespace storm { +namespace settings { +namespace modules { + +std::string const CvarSettings::moduleName = "cvar"; +std::string const CvarSettings::methodOptionName = "method"; +std::string const CvarSettings::interpretationOptionName = "interpretation"; + +CvarSettings::CvarSettings() : ModuleSettings(moduleName) { + std::vector methods = {"auto", "wr", "weighted-reachability", "ssp"}; + std::vector interpretations = {"auto", "cost", "reward"}; + this->addOption(storm::settings::OptionBuilder(moduleName, methodOptionName, true, "The method to be used for CVaR model checking.") + .setIsAdvanced() + .addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of the method to use.") + .addValidatorString(ArgumentValidatorFactory::createMultipleChoiceValidator(methods)) + .setDefaultValueString("auto") + .build()) + .build()); + this->addOption(storm::settings::OptionBuilder(moduleName, interpretationOptionName, true, "The interpretation to be used for CVaR model checking.") + .setIsAdvanced() + .addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of the interpretation to use.") + .addValidatorString(ArgumentValidatorFactory::createMultipleChoiceValidator(interpretations)) + .setDefaultValueString("auto") + .build()) + .build()); +} + +storm::modelchecker::cvar::CvarMethod CvarSettings::getCvarMethod() const { + std::string methodAsString = this->getOption(methodOptionName).getArgumentByName("name").getValueAsString(); + if (methodAsString == "auto") { + return storm::modelchecker::cvar::CvarMethod::Auto; + } else if (methodAsString == "wr" || methodAsString == "weighted-reachability") { + return storm::modelchecker::cvar::CvarMethod::WeightedReachability; + } else if (methodAsString == "ssp" || methodAsString == "ssp-vi" || methodAsString == "pareto-vi") { + return storm::modelchecker::cvar::CvarMethod::SspParetoVi; + } + STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown CVaR method '" << methodAsString << "'."); +} + +storm::modelchecker::cvar::CvarInterpretationSelection CvarSettings::getInterpretationSelection() const { + std::string interpretationAsString = this->getOption(interpretationOptionName).getArgumentByName("name").getValueAsString(); + if (interpretationAsString == "auto") { + return storm::modelchecker::cvar::CvarInterpretationSelection::Auto; + } else if (interpretationAsString == "cost") { + return storm::modelchecker::cvar::CvarInterpretationSelection::Cost; + } else if (interpretationAsString == "reward") { + return storm::modelchecker::cvar::CvarInterpretationSelection::Reward; + } + STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown CVaR interpretation '" << interpretationAsString << "'."); +} + +} // namespace modules +} // namespace settings +} // namespace storm diff --git a/src/storm/settings/modules/CvarSettings.h b/src/storm/settings/modules/CvarSettings.h new file mode 100644 index 0000000000..2e7f223004 --- /dev/null +++ b/src/storm/settings/modules/CvarSettings.h @@ -0,0 +1,47 @@ +#ifndef STORM_SETTINGS_MODULES_CVARSETTINGS_H_ +#define STORM_SETTINGS_MODULES_CVARSETTINGS_H_ + +#include "storm/modelchecker/cvar/CvarInterpretation.h" +#include "storm/modelchecker/cvar/CvarMethod.h" +#include "storm/settings/modules/ModuleSettings.h" + +namespace storm { +namespace settings { +namespace modules { + +/*! + * This class represents the settings for CVaR model checking. + */ +class CvarSettings : public ModuleSettings { + public: + /*! + * Creates a new set of CVaR model checking settings. + */ + CvarSettings(); + + /*! + * Retrieves the selected CVaR method. + * + * @return The selected CVaR method. + */ + storm::modelchecker::cvar::CvarMethod getCvarMethod() const; + + /*! + * Retrieves the selected CVaR interpretation policy. + * + * @return The selected CVaR interpretation policy. + */ + storm::modelchecker::cvar::CvarInterpretationSelection getInterpretationSelection() const; + + static std::string const moduleName; + + private: + static std::string const methodOptionName; + static std::string const interpretationOptionName; +}; + +} // namespace modules +} // namespace settings +} // namespace storm + +#endif /* STORM_SETTINGS_MODULES_CVARSETTINGS_H_ */ diff --git a/src/storm/settings/modules/IOSettings.cpp b/src/storm/settings/modules/IOSettings.cpp index b83b9884a0..74d1f5a564 100644 --- a/src/storm/settings/modules/IOSettings.cpp +++ b/src/storm/settings/modules/IOSettings.cpp @@ -59,6 +59,7 @@ const std::string IOSettings::qvbsInputOptionName = "qvbs"; const std::string IOSettings::qvbsInputOptionShortName = "qvbs"; const std::string IOSettings::qvbsRootOptionName = "qvbsroot"; const std::string IOSettings::propertiesAsMultiOptionName = "propsasmulti"; +const std::string IOSettings::cvarOptionName = "cvar"; const std::string IOSettings::uncertaintyResolutionModeName = "uncertainty-resolution"; @@ -281,6 +282,10 @@ IOSettings::IOSettings() : ModuleSettings(moduleName) { .setIsAdvanced() .build()); + this->addOption(storm::settings::OptionBuilder(moduleName, cvarOptionName, false, "Computes the conditional value-at-risk for the selected property.") + .addArgument(storm::settings::ArgumentBuilder::createStringArgument("alpha", "The size of the tail.").build()) + .build()); + std::vector uncertaintyResolutionModes = {"minimize", "maximize", "robust", "cooperative", "min", "max"}; this->addOption(storm::settings::OptionBuilder(moduleName, uncertaintyResolutionModeName, false, "Mode to resolve the uncertainty (intervals)") .addArgument(storm::settings::ArgumentBuilder::createStringArgument("mode", "Mode to resolve the uncertainty (intervals) by nature.") @@ -518,6 +523,14 @@ std::string IOSettings::getPropertyFilter() const { return this->getOption(propertyOptionName).getArgumentByName("filter").getValueAsString(); } +bool IOSettings::isCvarSet() const { + return this->getOption(cvarOptionName).getHasOptionBeenSet(); +} + +std::string IOSettings::getCvarAlpha() const { + return this->getOption(cvarOptionName).getArgumentByName("alpha").getValueAsString(); +} + bool IOSettings::isComputeSteadyStateDistributionSet() const { return this->getOption(steadyStateDistrOptionName).getHasOptionBeenSet(); } diff --git a/src/storm/settings/modules/IOSettings.h b/src/storm/settings/modules/IOSettings.h index 94d3f33c0e..b78cdc00c3 100644 --- a/src/storm/settings/modules/IOSettings.h +++ b/src/storm/settings/modules/IOSettings.h @@ -366,6 +366,20 @@ class IOSettings : public ModuleSettings { */ std::string getPropertyFilter() const; + /*! + * Retrieves whether the CVaR option was set. + * + * @return True if the CVaR option was set. + */ + bool isCvarSet() const; + + /*! + * Retrieves the alpha value specified with the CVaR option. + * + * @return The alpha value specified with the CVaR option. + */ + std::string getCvarAlpha() const; + /*! * Retrieves whether the steady-state distribution is to be computed. */ @@ -469,6 +483,7 @@ class IOSettings : public ModuleSettings { static const std::string qvbsInputOptionShortName; static const std::string qvbsRootOptionName; static const std::string propertiesAsMultiOptionName; + static const std::string cvarOptionName; static const std::string uncertaintyResolutionModeName; }; diff --git a/src/storm/storage/jani/visitor/JSONExporter.cpp b/src/storm/storage/jani/visitor/JSONExporter.cpp index 3ec72f6b95..6ad2855175 100644 --- a/src/storm/storage/jani/visitor/JSONExporter.cpp +++ b/src/storm/storage/jani/visitor/JSONExporter.cpp @@ -468,6 +468,10 @@ boost::any FormulaToJaniJson::visit(storm::logic::QuantileFormula const&, boost: STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "Jani currently does not support conversion of a Quantile formula."); } +boost::any FormulaToJaniJson::visit(storm::logic::CvarFormula const&, boost::any const&) const { + STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "JANI export does not support Storm's CLI-only CVaR property wrapper."); +} + boost::any FormulaToJaniJson::visit(storm::logic::NextFormula const& f, boost::any const& data) const { ExportJsonType opDecl; opDecl["op"] = "U"; diff --git a/src/storm/storage/jani/visitor/JSONExporter.h b/src/storm/storage/jani/visitor/JSONExporter.h index 6fe16eb5e3..cb5ad50f70 100644 --- a/src/storm/storage/jani/visitor/JSONExporter.h +++ b/src/storm/storage/jani/visitor/JSONExporter.h @@ -59,6 +59,7 @@ class FormulaToJaniJson : public storm::logic::FormulaVisitor { virtual boost::any visit(storm::logic::BoundedUntilFormula const& f, boost::any const& data) const; virtual boost::any visit(storm::logic::ConditionalFormula const& f, boost::any const& data) const; virtual boost::any visit(storm::logic::CumulativeRewardFormula const& f, boost::any const& data) const; + virtual boost::any visit(storm::logic::CvarFormula const& f, boost::any const& data) const; virtual boost::any visit(storm::logic::EventuallyFormula const& f, boost::any const& data) const; virtual boost::any visit(storm::logic::TimeOperatorFormula const& f, boost::any const& data) const; virtual boost::any visit(storm::logic::GameFormula const& f, boost::any const& data) const; diff --git a/src/test/storm/modelchecker/prctl/mdp/CvarQueryTest.cpp b/src/test/storm/modelchecker/prctl/mdp/CvarQueryTest.cpp new file mode 100644 index 0000000000..11f1cafa57 --- /dev/null +++ b/src/test/storm/modelchecker/prctl/mdp/CvarQueryTest.cpp @@ -0,0 +1,513 @@ +#include "storm-config.h" +#include "test/storm_gtest.h" + +#include + +#include "storm-parsers/api/model_descriptions.h" +#include "storm-parsers/api/properties.h" +#include "storm/adapters/RationalNumberAdapter.h" +#include "storm/api/builder.h" +#include "storm/api/properties.h" +#include "storm/environment/Environment.h" +#include "storm/environment/modelchecker/ModelCheckerEnvironment.h" +#include "storm/exceptions/InvalidArgumentException.h" +#include "storm/exceptions/InvalidPropertyException.h" +#include "storm/exceptions/NotImplementedException.h" +#include "storm/logic/CvarFormula.h" +#include "storm/modelchecker/CheckTask.h" +#include "storm/modelchecker/cvar/CvarInterpretation.h" +#include "storm/modelchecker/cvar/CvarMethod.h" +#include "storm/modelchecker/cvar/helper/SspParetoFront.h" +#include "storm/modelchecker/cvar/helper/SspParetoValueIterationOperator.h" +#include "storm/modelchecker/cvar/preprocessing/SspCvarPreprocessingResult.h" +#include "storm/modelchecker/prctl/SparseMdpPrctlModelChecker.h" +#include "storm/modelchecker/results/ExplicitQuantitativeCheckResult.h" +#include "storm/models/sparse/Mdp.h" +#include "storm/models/sparse/StandardRewardModel.h" +#include "storm/utility/constants.h" + +namespace { +bool hasLpSolver() { +#if !defined(STORM_HAVE_GLPK) && !defined(STORM_HAVE_GUROBI) && !defined(STORM_HAVE_Z3) && !defined(STORM_HAVE_SOPLEX) + return false; +#else + return true; +#endif +} + +bool hasExactLpSolver() { +#if !defined(STORM_HAVE_Z3) && !defined(STORM_HAVE_SOPLEX) + return false; +#else + return true; +#endif +} + +class CvarQueryTest : public ::testing::Test { + public: + void SetUp() override { +#ifndef STORM_HAVE_Z3 + GTEST_SKIP() << "Z3 not available."; +#endif + } +}; + +template +struct CvarTestInput { + std::shared_ptr> mdp; + std::shared_ptr formula; +}; + +template +CvarTestInput buildCvarInput(std::string const& modelPath, std::string const& propertyString, std::string const& alpha) { + storm::prism::Program program = storm::api::parseProgram(modelPath); + auto properties = storm::api::parsePropertiesForPrismProgram(propertyString, program); + std::vector cvarProperties = {storm::api::createCvarProperty(properties.front(), alpha)}; + auto formulas = storm::api::extractFormulasFromProperties(cvarProperties); + auto mdp = storm::api::buildSparseModel(program, formulas)->template as>(); + return {mdp, cvarProperties.front().getRawFormula()}; +} + +template +std::unique_ptr checkInitialStateResult(CvarTestInput const& input) { + storm::Environment env; + storm::modelchecker::SparseMdpPrctlModelChecker> checker(*input.mdp); + storm::modelchecker::CheckTask task(*input.formula, true); + return checker.check(env, task); +} + +template +ValueType checkInitialStateValue(CvarTestInput const& input) { + auto result = checkInitialStateResult(input); + return result->template asExplicitQuantitativeCheckResult().getMax(); +} + +template +ValueType checkInitialStateValueWithMethod(CvarTestInput const& input, storm::modelchecker::cvar::CvarMethod method) { + storm::Environment env; + env.modelchecker().cvar().setMethod(method); + env.modelchecker().cvar().setInterpretationSelection(storm::modelchecker::cvar::CvarInterpretationSelection::Auto); + storm::modelchecker::SparseMdpPrctlModelChecker> checker(*input.mdp); + storm::modelchecker::CheckTask task(*input.formula, true); + auto result = checker.check(env, task); + return result->template asExplicitQuantitativeCheckResult().getMax(); +} + +template +ValueType checkInitialStateValueWithMethodAndInterpretationSelection(CvarTestInput const& input, storm::modelchecker::cvar::CvarMethod method, + storm::modelchecker::cvar::CvarInterpretationSelection interpretationSelection) { + storm::Environment env; + env.modelchecker().cvar().setMethod(method); + env.modelchecker().cvar().setInterpretationSelection(interpretationSelection); + storm::modelchecker::SparseMdpPrctlModelChecker> checker(*input.mdp); + storm::modelchecker::CheckTask task(*input.formula, true); + auto result = checker.check(env, task); + return result->template asExplicitQuantitativeCheckResult().getMax(); +} + +template +void expectParetoFrontPoints(ParetoFront const& front, std::vector> const& expectedPoints) { + auto const& actualPoints = front.getPoints(); + ASSERT_EQ(expectedPoints.size(), actualPoints.size()); + for (std::size_t index = 0; index < expectedPoints.size(); ++index) { + EXPECT_NEAR(expectedPoints[index].first, actualPoints[index].probability, 1e-10); + EXPECT_NEAR(expectedPoints[index].second, actualPoints[index].expectedCost, 1e-10); + } +} + +storm::modelchecker::cvar::preprocessing::SspCvarPreprocessingResult buildTinySspPreprocessingResult() { + storm::storage::SparseMatrixBuilder builder(4, 3, 4, true, true, 3); + builder.newRowGroup(0); + builder.addNextValue(0, 1, 1.0); + builder.addNextValue(1, 1, 1.0); + builder.newRowGroup(2); + builder.addNextValue(2, 2, 1.0); + builder.newRowGroup(3); + builder.addNextValue(3, 2, 1.0); + + storm::modelchecker::cvar::preprocessing::SspCvarPreprocessingResult result; + result.rewardModelName = "cost"; + result.initialState = 0; + result.targetStates = storm::storage::BitVector(3, false); + result.targetStates.set(2, true); + result.reachableStates = storm::storage::BitVector(3, true); + result.liftedStateRewardsToChoiceCosts = true; + result.normalizedTargetStatesToAbsorbing = true; + result.maximalChoiceCost = 2; + result.choiceCosts = {1.0, 2.0, 1.0, 0.0}; + result.expectedCostsToGoal = {0.0, 0.0, 0.0}; + result.transitionMatrix = builder.build(); + return result; +} + +TEST(CvarSspParetoFrontTest, CanonicalizesDuplicateDominatedAndConvexRedundantPoints) { + using ParetoFront = storm::modelchecker::cvar::SspParetoFront; + + ParetoFront front({{0.4, 5.0}, {0.2, 4.0}, {0.2, 3.0}, {0.6, 6.0}, {0.8, 9.0}, {0.5, 7.0}}); + + expectParetoFrontPoints(front, {{0.2, 3.0}, {0.6, 6.0}, {0.8, 9.0}}); +} + +TEST(CvarSspParetoFrontTest, ScaledMinkowskiSumMergesConvexChains) { + using ParetoFront = storm::modelchecker::cvar::SspParetoFront; + + ParetoFront left({{0.0, 0.0}, {0.25, 1.0}, {0.5, 3.0}}); + ParetoFront right({{0.0, 0.0}, {0.25, 0.5}, {0.5, 2.0}}); + + auto result = left.minkowskiSumScaled(right, 0.5); + + expectParetoFrontPoints(result, {{0.0, 0.0}, {0.125, 0.25}, {0.375, 1.25}, {0.5, 2.0}, {0.75, 4.0}}); +} + +TEST(CvarSspParetoFrontTest, CanonicalizesAfterFloatingPointTranslationCollapse) { + using ParetoFront = storm::modelchecker::cvar::SspParetoFront; + + ParetoFront front({{0.0, 0.0}, {1e-17, 1.0}}); + + auto result = front.minkowskiSum(ParetoFront::singleton(1.0, 0.0)); + + expectParetoFrontPoints(result, {{1.0, 0.0}}); +} + +TEST(CvarSspParetoValueIterationOperatorTest, AppliesActionCostsAndUnionsActionFronts) { + using ParetoFront = storm::modelchecker::cvar::SspParetoFront; + using ParetoViOperator = storm::modelchecker::cvar::SspParetoValueIterationOperator; + using FrontierLayer = std::vector; + using FrontierWindow = std::vector; + + auto preprocessingResult = buildTinySspPreprocessingResult(); + ParetoViOperator paretoViOperator(preprocessingResult); + FrontierWindow frontierWindow(3, FrontierLayer(3)); + for (auto& layer : frontierWindow) { + layer[2] = ParetoViOperator::createTargetFrontier(); + } + frontierWindow[0][1] = ParetoFront::singleton(0.8, 4.0); + frontierWindow[1][1] = ParetoFront::singleton(0.5, 1.0); + + FrontierLayer outputLayer; + paretoViOperator.apply(2, frontierWindow, outputLayer); + + expectParetoFrontPoints(outputLayer[0], {{0.5, 1.0}, {0.8, 4.0}}); + expectParetoFrontPoints(outputLayer[1], {{1.0, 0.0}}); + expectParetoFrontPoints(outputLayer[2], {{1.0, 0.0}}); +} + +TEST(CvarSspRewardParetoFrontTest, KeepsLowerLeftIncomparablePoints) { + using ParetoFront = storm::modelchecker::cvar::SspRewardParetoFront; + + ParetoFront front({{0.1, 0.9}, {0.5, 0.5}, {0.7, 0.8}, {0.1, 1.1}, {0.9, 0.2}}); + + expectParetoFrontPoints(front, {{0.1, 0.9}, {0.5, 0.5}, {0.9, 0.2}}); +} + +TEST(CvarSspRewardParetoValueIterationOperatorTest, AppliesRewardShiftsAndUnionsActionFronts) { + using ParetoFront = storm::modelchecker::cvar::SspRewardParetoFront; + using ParetoViOperator = storm::modelchecker::cvar::SspParetoValueIterationOperator; + using FrontierLayer = std::vector; + using FrontierWindow = std::vector; + + auto preprocessingResult = buildTinySspPreprocessingResult(); + ParetoViOperator paretoViOperator(preprocessingResult); + FrontierWindow frontierWindow(3, FrontierLayer(3)); + for (uint64_t threshold = 0; threshold < frontierWindow.size(); ++threshold) { + frontierWindow[threshold][2] = ParetoViOperator::createTargetFrontier(threshold); + } + frontierWindow[2][1] = ParetoFront::singleton(0.5, 0.1); + frontierWindow[1][1] = ParetoFront::singleton(0.2, 0.5); + + FrontierLayer outputLayer; + paretoViOperator.apply(3, frontierWindow, outputLayer); + + expectParetoFrontPoints(outputLayer[0], {{0.2, 0.5}, {0.5, 0.1}}); + expectParetoFrontPoints(outputLayer[1], {{1.0, 2.0}}); + expectParetoFrontPoints(outputLayer[2], {{1.0, 3.0}}); +} + +TEST_F(CvarQueryTest, SimpleMdp) { + if (!hasLpSolver()) { + GTEST_SKIP() << "No LP solver available."; + } + + std::string alpha = "0.75"; + std::string modelPath = STORM_TEST_RESOURCES_DIR "/mdp/cvar_simple_mdp.nm"; + + auto maxInput = buildCvarInput(modelPath, "R{\"term\"}max=? [ F \"target\" ];", alpha); + double maxValue = checkInitialStateValue(maxInput); + EXPECT_NEAR(maxValue, 2.0, 1e-10); + + auto minInput = buildCvarInput(modelPath, "R{\"term\"}min=? [ F \"target\" ];", alpha); + double minValue = checkInitialStateValue(minInput); + EXPECT_NEAR(minValue, 2.0, 1e-10); +} + +TEST_F(CvarQueryTest, ReachableBadMecIsPreprocessedToZeroTerminalReward) { + if (!hasLpSolver()) { + GTEST_SKIP() << "No LP solver available."; + } + + std::string alpha = "0.5"; + std::string modelPath = STORM_TEST_RESOURCES_DIR "/mdp/cvar_bad_mec_mdp.nm"; + + auto maxInput = buildCvarInput(modelPath, "R{\"term\"}max=? [ F \"target\" ];", alpha); + double maxValue = checkInitialStateValue(maxInput); + EXPECT_NEAR(maxValue, 0.0, 1e-10); + + auto minInput = buildCvarInput(modelPath, "R{\"term\"}min=? [ F \"target\" ];", alpha); + double minValue = checkInitialStateValue(minInput); + EXPECT_NEAR(minValue, 4.0, 1e-10); +} + +TEST_F(CvarQueryTest, TargetReachingMecIsCollapsed) { + if (!hasLpSolver()) { + GTEST_SKIP() << "No LP solver available."; + } + + std::string modelPath = STORM_TEST_RESOURCES_DIR "/mdp/cvar_target_reaching_mec_mdp.nm"; + + auto maxInput = buildCvarInput(modelPath, "R{\"term\"}max=? [ F \"target\" ];", "0.75"); + EXPECT_NEAR(checkInitialStateValue(maxInput), 6.0, 1e-10); + + auto minInput = buildCvarInput(modelPath, "R{\"term\"}min=? [ F \"target\" ];", "0.75"); + EXPECT_NEAR(checkInitialStateValue(minInput), 4.0, 1e-10); +} + +TEST_F(CvarQueryTest, RejectsNonAbsorbingOriginalTargetStates) { + std::string alpha = "0.5"; + std::string modelPath = STORM_TEST_RESOURCES_DIR "/mdp/cvar_nonabsorbing_target_mdp.nm"; + + auto input = buildCvarInput(modelPath, "R{\"term\"}max=? [ F \"target\" ];", alpha); + + storm::Environment env; + storm::modelchecker::SparseMdpPrctlModelChecker> checker(*input.mdp); + storm::modelchecker::CheckTask task(*input.formula, true); + STORM_SILENT_EXPECT_THROW(checker.check(env, task), storm::exceptions::InvalidPropertyException); +} + +TEST_F(CvarQueryTest, BranchingTradeoffMdp) { + if (!hasLpSolver()) { + GTEST_SKIP() << "No LP solver available."; + } + + std::string modelPath = STORM_TEST_RESOURCES_DIR "/mdp/cvar_branching_tradeoff_mdp.nm"; + + auto maxHalfInput = buildCvarInput(modelPath, "R{\"term\"}max=? [ F \"target\" ];", "0.5"); + EXPECT_NEAR(checkInitialStateValue(maxHalfInput), 7.0, 1e-10); + + auto maxThreeQuarterInput = buildCvarInput(modelPath, "R{\"term\"}max=? [ F \"target\" ];", "0.75"); + EXPECT_NEAR(checkInitialStateValue(maxThreeQuarterInput), 8.0, 1e-10); + + auto minHalfInput = buildCvarInput(modelPath, "R{\"term\"}min=? [ F \"target\" ];", "0.5"); + EXPECT_NEAR(checkInitialStateValue(minHalfInput), 7.0, 1e-10); + + auto minThreeQuarterInput = buildCvarInput(modelPath, "R{\"term\"}min=? [ F \"target\" ];", "0.75"); + EXPECT_NEAR(checkInitialStateValue(minThreeQuarterInput), 7.0, 1e-10); +} + +TEST_F(CvarQueryTest, BranchingTradeoffMdpRationalNumbers) { + if (!hasExactLpSolver()) { + GTEST_SKIP() << "No exact LP solver available."; + } + + std::string modelPath = STORM_TEST_RESOURCES_DIR "/mdp/cvar_branching_tradeoff_mdp.nm"; + + auto maxInput = buildCvarInput(modelPath, "R{\"term\"}max=? [ F \"target\" ];", "0.75"); + EXPECT_EQ(storm::RationalNumber(8), checkInitialStateValue(maxInput)); + + auto minInput = buildCvarInput(modelPath, "R{\"term\"}min=? [ F \"target\" ];", "0.75"); + EXPECT_EQ(storm::RationalNumber(7), checkInitialStateValue(minInput)); +} + +TEST_F(CvarQueryTest, InterpretationOverridesOnSimpleMdp) { + if (!hasLpSolver()) { + GTEST_SKIP() << "No LP solver available."; + } + + std::string modelPath = STORM_TEST_RESOURCES_DIR "/mdp/cvar_simple_mdp.nm"; + + auto maxInput = buildCvarInput(modelPath, "R{\"term\"}max=? [ F \"target\" ];", "0.75"); + EXPECT_NEAR(checkInitialStateValueWithMethodAndInterpretationSelection(maxInput, storm::modelchecker::cvar::CvarMethod::WeightedReachability, + storm::modelchecker::cvar::CvarInterpretationSelection::Reward), + 2.0, 1e-10); + EXPECT_NEAR(checkInitialStateValueWithMethodAndInterpretationSelection(maxInput, storm::modelchecker::cvar::CvarMethod::WeightedReachability, + storm::modelchecker::cvar::CvarInterpretationSelection::Cost), + 7.0 / 3.0, 1e-10); + + auto minInput = buildCvarInput(modelPath, "R{\"term\"}min=? [ F \"target\" ];", "0.75"); + EXPECT_NEAR(checkInitialStateValueWithMethodAndInterpretationSelection(minInput, storm::modelchecker::cvar::CvarMethod::WeightedReachability, + storm::modelchecker::cvar::CvarInterpretationSelection::Cost), + 2.0, 1e-10); + EXPECT_NEAR(checkInitialStateValueWithMethodAndInterpretationSelection(minInput, storm::modelchecker::cvar::CvarMethod::WeightedReachability, + storm::modelchecker::cvar::CvarInterpretationSelection::Reward), + 5.0 / 3.0, 1e-10); +} + +TEST_F(CvarQueryTest, InterpretationOverridesOnSimpleMdpRationalNumbers) { + if (!hasExactLpSolver()) { + GTEST_SKIP() << "No exact LP solver available."; + } + + std::string modelPath = STORM_TEST_RESOURCES_DIR "/mdp/cvar_simple_mdp.nm"; + + auto maxInput = buildCvarInput(modelPath, "R{\"term\"}max=? [ F \"target\" ];", "0.75"); + EXPECT_EQ(storm::RationalNumber(2), + checkInitialStateValueWithMethodAndInterpretationSelection(maxInput, storm::modelchecker::cvar::CvarMethod::WeightedReachability, + storm::modelchecker::cvar::CvarInterpretationSelection::Reward)); + EXPECT_EQ(storm::RationalNumber("7/3"), + checkInitialStateValueWithMethodAndInterpretationSelection(maxInput, storm::modelchecker::cvar::CvarMethod::WeightedReachability, + storm::modelchecker::cvar::CvarInterpretationSelection::Cost)); + + auto minInput = buildCvarInput(modelPath, "R{\"term\"}min=? [ F \"target\" ];", "0.75"); + EXPECT_EQ(storm::RationalNumber(2), + checkInitialStateValueWithMethodAndInterpretationSelection(minInput, storm::modelchecker::cvar::CvarMethod::WeightedReachability, + storm::modelchecker::cvar::CvarInterpretationSelection::Cost)); + EXPECT_EQ(storm::RationalNumber("5/3"), + checkInitialStateValueWithMethodAndInterpretationSelection(minInput, storm::modelchecker::cvar::CvarMethod::WeightedReachability, + storm::modelchecker::cvar::CvarInterpretationSelection::Reward)); +} + +TEST_F(CvarQueryTest, EquivalentExactAlphaSyntaxes) { + if (!hasLpSolver()) { + GTEST_SKIP() << "No LP solver available."; + } + + std::string modelPath = STORM_TEST_RESOURCES_DIR "/mdp/cvar_simple_mdp.nm"; + + auto decimalInput = buildCvarInput(modelPath, "R{\"term\"}max=? [ F \"target\" ];", "0.75"); + auto fractionInput = buildCvarInput(modelPath, "R{\"term\"}max=? [ F \"target\" ];", "3/4"); + auto scientificInput = buildCvarInput(modelPath, "R{\"term\"}max=? [ F \"target\" ];", "7.5e-1"); + + EXPECT_NEAR(checkInitialStateValue(decimalInput), 2.0, 1e-10); + EXPECT_NEAR(checkInitialStateValue(fractionInput), 2.0, 1e-10); + EXPECT_NEAR(checkInitialStateValue(scientificInput), 2.0, 1e-10); +} + +TEST(CvarFormulaTest, RejectsInvalidAlphaSyntaxes) { + std::string modelPath = STORM_TEST_RESOURCES_DIR "/mdp/cvar_simple_mdp.nm"; + storm::prism::Program program = storm::api::parseProgram(modelPath); + auto properties = storm::api::parsePropertiesForPrismProgram("R{\"term\"}max=? [ F \"target\" ];", program); + + for (auto const& alpha : {"0", "1", "-0.1", "abc", "0/1", "1/1"}) { + STORM_SILENT_EXPECT_THROW(storm::api::createCvarProperty(properties.front(), std::string(alpha)), storm::exceptions::BaseException); + } +} + +TEST(CvarFormulaTest, CvarFormulaValidatesAlphaAndSubformula) { + std::string modelPath = STORM_TEST_RESOURCES_DIR "/mdp/cvar_simple_mdp.nm"; + storm::prism::Program program = storm::api::parseProgram(modelPath); + auto properties = storm::api::parsePropertiesForPrismProgram("R{\"term\"}max=? [ F \"target\" ];", program); + auto formula = properties.front().getRawFormula(); + + EXPECT_NO_THROW(storm::logic::CvarFormula(storm::RationalNumber("3/4"), formula)); + STORM_SILENT_EXPECT_THROW(storm::logic::CvarFormula(storm::utility::zero(), formula), storm::exceptions::InvalidArgumentException); + STORM_SILENT_EXPECT_THROW(storm::logic::CvarFormula(storm::utility::one(), formula), storm::exceptions::InvalidArgumentException); + STORM_SILENT_EXPECT_THROW(storm::logic::CvarFormula(storm::RationalNumber("1/2"), nullptr), storm::exceptions::InvalidArgumentException); +} + +TEST_F(CvarQueryTest, CheckTaskExtractsWrappedRewardMetadata) { + std::string modelPath = STORM_TEST_RESOURCES_DIR "/mdp/cvar_simple_mdp.nm"; + auto input = buildCvarInput(modelPath, "R{\"term\"}max=? [ F \"target\" ];", "3/4"); + + auto const& cvarFormula = input.formula->asCvarFormula(); + storm::modelchecker::CheckTask task(cvarFormula, true); + + ASSERT_TRUE(task.isOptimizationDirectionSet()); + EXPECT_EQ(storm::solver::OptimizationDirection::Maximize, task.getOptimizationDirection()); + ASSERT_TRUE(task.isRewardModelSet()); + EXPECT_EQ("term", task.getRewardModel()); + EXPECT_FALSE(task.isBoundSet()); +} + +TEST_F(CvarQueryTest, DeterministicSspPathMdp) { + std::string modelPath = STORM_TEST_RESOURCES_DIR "/mdp/cvar_ssp_deterministic_mdp.nm"; + auto input = buildCvarInput(modelPath, "R{\"cost\"}min=? [ F \"goal\" ];", "0.5"); + + double value = checkInitialStateValueWithMethod(input, storm::modelchecker::cvar::CvarMethod::SspParetoVi); + EXPECT_NEAR(value, 5.0, 1e-10); +} + +TEST_F(CvarQueryTest, DeterministicRewardSspPathMdp) { + std::string modelPath = STORM_TEST_RESOURCES_DIR "/mdp/cvar_ssp_deterministic_mdp.nm"; + auto input = buildCvarInput(modelPath, "R{\"cost\"}max=? [ F \"goal\" ];", "0.5"); + + EXPECT_NEAR(checkInitialStateValueWithMethod(input, storm::modelchecker::cvar::CvarMethod::SspParetoVi), 5.0, 1e-10); + EXPECT_NEAR(checkInitialStateValueWithMethodAndInterpretationSelection(input, storm::modelchecker::cvar::CvarMethod::SspParetoVi, + storm::modelchecker::cvar::CvarInterpretationSelection::Reward), + 5.0, 1e-10); +} + +TEST_F(CvarQueryTest, BranchingSspTradeoffMdp) { + std::string modelPath = STORM_TEST_RESOURCES_DIR "/mdp/cvar_ssp_branching_tradeoff_mdp.nm"; + + auto halfInput = buildCvarInput(modelPath, "R{\"cost\"}min=? [ F \"goal\" ];", "0.5"); + EXPECT_NEAR(checkInitialStateValueWithMethod(halfInput, storm::modelchecker::cvar::CvarMethod::SspParetoVi), 6.0, 1e-10); + + auto nineTenthsInput = buildCvarInput(modelPath, "R{\"cost\"}min=? [ F \"goal\" ];", "0.9"); + EXPECT_NEAR(checkInitialStateValueWithMethod(nineTenthsInput, storm::modelchecker::cvar::CvarMethod::SspParetoVi), 16.0 / 3.0, 1e-10); +} + +TEST_F(CvarQueryTest, SafeRiskyRewardSspMdp) { + std::string modelPath = STORM_TEST_RESOURCES_DIR "/mdp/cvar_ssp_reward_safe_risky_mdp.nm"; + + auto quarterInput = buildCvarInput(modelPath, "R{\"reward\"}max=? [ F \"goal\" ];", "0.25"); + EXPECT_NEAR(checkInitialStateValueWithMethod(quarterInput, storm::modelchecker::cvar::CvarMethod::SspParetoVi), 6.0, 1e-10); + + auto fourFifthsInput = buildCvarInput(modelPath, "R{\"reward\"}max=? [ F \"goal\" ];", "0.8"); + EXPECT_NEAR(checkInitialStateValueWithMethod(fourFifthsInput, storm::modelchecker::cvar::CvarMethod::SspParetoVi), 8.0, 1e-10); +} + +TEST_F(CvarQueryTest, DelayedChallengerRewardSspMdp) { + std::string modelPath = STORM_TEST_RESOURCES_DIR "/mdp/cvar_ssp_reward_delayed_challenger_mdp.nm"; + auto input = buildCvarInput(modelPath, "R{\"reward\"}max=? [ F \"goal\" ];", "0.5"); + + EXPECT_NEAR(checkInitialStateValueWithMethod(input, storm::modelchecker::cvar::CvarMethod::SspParetoVi), 8.0, 1e-10); +} + +TEST_F(CvarQueryTest, GeometricRewardSspMdp) { + std::string modelPath = STORM_TEST_RESOURCES_DIR "/mdp/cvar_ssp_reward_geometric_mdp.nm"; + auto input = buildCvarInput(modelPath, "R{\"reward\"}max=? [ F \"goal\" ];", "0.8"); + + EXPECT_NEAR(checkInitialStateValueWithMethod(input, storm::modelchecker::cvar::CvarMethod::SspParetoVi), 1.4375, 1e-10); +} + +TEST_F(CvarQueryTest, RejectsRewardSspProb1AChoiceViolation) { + std::string modelPath = STORM_TEST_RESOURCES_DIR "/mdp/cvar_ssp_reward_prob1a_choice_violation_mdp.nm"; + auto input = buildCvarInput(modelPath, "R{\"reward\"}max=? [ F \"goal\" ];", "0.5"); + + STORM_SILENT_EXPECT_THROW(checkInitialStateValueWithMethod(input, storm::modelchecker::cvar::CvarMethod::SspParetoVi), + storm::exceptions::InvalidPropertyException); +} + +TEST_F(CvarQueryTest, RejectsRewardSspProb1AProbabilisticViolation) { + std::string modelPath = STORM_TEST_RESOURCES_DIR "/mdp/cvar_ssp_reward_prob1a_probabilistic_violation_mdp.nm"; + auto input = buildCvarInput(modelPath, "R{\"reward\"}max=? [ F \"goal\" ];", "0.5"); + + STORM_SILENT_EXPECT_THROW(checkInitialStateValueWithMethod(input, storm::modelchecker::cvar::CvarMethod::SspParetoVi), + storm::exceptions::InvalidPropertyException); +} + +TEST_F(CvarQueryTest, RejectsUnsupportedSspInterpretationCombinations) { + std::string modelPath = STORM_TEST_RESOURCES_DIR "/mdp/cvar_ssp_deterministic_mdp.nm"; + + auto minRewardInput = buildCvarInput(modelPath, "R{\"cost\"}min=? [ F \"goal\" ];", "0.5"); + STORM_SILENT_EXPECT_THROW(checkInitialStateValueWithMethodAndInterpretationSelection(minRewardInput, storm::modelchecker::cvar::CvarMethod::SspParetoVi, + storm::modelchecker::cvar::CvarInterpretationSelection::Reward), + storm::exceptions::InvalidPropertyException); + + auto maxCostInput = buildCvarInput(modelPath, "R{\"cost\"}max=? [ F \"goal\" ];", "0.5"); + STORM_SILENT_EXPECT_THROW(checkInitialStateValueWithMethodAndInterpretationSelection(maxCostInput, storm::modelchecker::cvar::CvarMethod::SspParetoVi, + storm::modelchecker::cvar::CvarInterpretationSelection::Cost), + storm::exceptions::InvalidPropertyException); +} + +TEST_F(CvarQueryTest, RejectsSspTransitionRewards) { + std::string modelPath = STORM_TEST_RESOURCES_DIR "/mdp/cvar_ssp_deterministic_mdp.nm"; + auto input = buildCvarInput(modelPath, "R{\"cost\"}max=? [ F \"goal\" ];", "0.5"); + + storm::models::sparse::StandardRewardModel transitionRewardModel( + std::nullopt, std::make_optional(std::vector(input.mdp->getNumberOfChoices(), 1.0)), std::make_optional(input.mdp->getTransitionMatrix())); + input.mdp->getRewardModels()["cost"] = std::move(transitionRewardModel); + + STORM_SILENT_EXPECT_THROW(checkInitialStateValueWithMethod(input, storm::modelchecker::cvar::CvarMethod::SspParetoVi), + storm::exceptions::NotImplementedException); +} +} // namespace