Browse Source

Merge branch 'master' into sound-vi

tempestpy_adaptions
TimQu 7 years ago
parent
commit
ce4d75b142
  1. 24
      CHANGELOG.md
  2. 36
      doc/checklist_new_release.md
  3. 7
      resources/3rdparty/CMakeLists.txt
  4. 2
      resources/3rdparty/carl/CMakeLists.txt
  5. 2
      resources/3rdparty/sylvan/CMakeLists.txt
  6. 2
      src/storm/adapters/Smt2ExpressionAdapter.h
  7. 159
      src/storm/counterexamples/SMTMinimalLabelSetGenerator.h
  8. 57
      src/storm/logic/ComparisonType.h
  9. 24
      src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.cpp
  10. 30
      src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp
  11. 5
      src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.h
  12. 5
      src/storm/solver/AbstractEquationSolver.cpp
  13. 5
      src/storm/solver/AbstractEquationSolver.h
  14. 4
      src/storm/solver/NativeLinearEquationSolver.cpp
  15. 2
      src/storm/solver/SmtlibSmtSolver.cpp
  16. 8
      src/storm/solver/SymbolicMinMaxLinearEquationSolver.cpp
  17. 10
      src/storm/solver/SymbolicNativeLinearEquationSolver.cpp
  18. 6
      src/storm/storage/dd/Add.cpp
  19. 7
      src/storm/storage/dd/cudd/InternalCuddAdd.cpp
  20. 5
      src/storm/storage/dd/cudd/InternalCuddAdd.h
  21. 7
      src/storm/storage/dd/sylvan/InternalSylvanAdd.cpp
  22. 2
      src/storm/storage/dd/sylvan/InternalSylvanAdd.h
  23. 2
      src/storm/storage/sparse/PrismChoiceOrigins.cpp
  24. 2
      src/storm/storage/sparse/PrismChoiceOrigins.h
  25. 275
      src/storm/utility/numerical.cpp
  26. 134
      src/storm/utility/numerical.h
  27. 2
      version.cmake

24
CHANGELOG.md

@ -4,11 +4,12 @@ Changelog
This changelog lists only the most important changes. Smaller (bug)fixes as well as non-mature features are not part of the changelog.
The releases of major and minor versions contain an overview of changes since the last major/minor update.
Version 1.2.x
-------------
### Version 1.2
- C++ api changes: Building model takes BuilderOptions instead of extended list of Booleans, does not depend on settings anymore.
- storm-cli-utilities now contains cli related stuff, instead of storm-lib
### Version 1.2.0 (2017/12)
- C++ api changes: Building model takes `BuilderOptions` instead of extended list of Booleans, does not depend on settings anymore.
- `storm-cli-utilities` now contains cli related stuff, instead of `storm-lib`
- Symbolic (MT/BDD) bisimulation
- Fixed issue related to variable names that can not be used in Exprtk.
- DRN parser improved
@ -20,9 +21,10 @@ The releases of major and minor versions contain an overview of changes since th
- Performance improvements for sparse model building
- Performance improvements for conditional properties on MDPs
- Automatically convert MA without probabilistic states into CTMC
- storm-pars: support for welldefinedness constraints in mdps.
- storm-dft: split DFT settings into IO settings and fault tree settings
- storm-dft: removed obsolete explicit model builder for DFTs
- Fixed implemention of Fox and Glynn' algorithm
- `storm-pars`: support for welldefinedness constraints in mdps.
- `storm-dft`: split DFT settings into IO settings and fault tree settings
- `storm-dft`: removed obsolete explicit model builder for DFTs
- Features for developers:
* Solvers can now expose requirements
* unbounded reachability and reachability rewards now correctly respect solver requirements
@ -34,10 +36,9 @@ Version 1.1.x
-------------
### Version 1.1.0 (2017/8)
- Support for long-run average rewards on MDPs and Markov automata using a value-iteration based approach.
- Storm can now check MDPs and Markov Automata (i.e. MinMax equation systems) via Linear Programming.
- Parametric model checking is now handled in a separated library/executable called storm-pars.
- Parametric model checking is now handled in a separated library/executable called `storm-pars`.
- Wellformedness constraints on PMCs:
* include constraints from rewards
* are in smtlib2
@ -48,15 +49,14 @@ Version 1.1.x
- Support for parsing/building models given in the explicit input format of IMCA.
- Storm now overwrites files if asked to write files to a specific location.
- Changes in build process to accommodate for changes in carl. Also, more robust against issues with carl.
- USE_POPCNT removed in favor of FORCE_POPCNT. The popcnt instruction is used if available due to march=native, unless portable is set.
Then, using FORCE_POPCNT enables the use of the SSE 4.2 instruction
- `USE_POPCNT` removed in favor of `FORCE_POPCNT`. The popcnt instruction is used if available due to `march=native`, unless portable is set.
Then, using `FORCE_POPCNT` enables the use of the SSE 4.2 instruction
Version 1.0.x
-------------
### Version 1.0.1 (2017/4)
- Multi-objective model checking support now fully included
- Several improvements in parameter lifting
- Several improvements in JANI parsing

36
doc/checklist_new_release.md

@ -0,0 +1,36 @@
The following steps should be performed before releasing a new storm version.
Note that in most case a simultaneous release of [carl](https://github.com/smtrat/carl), [storm](https://github.com/moves-rwth/storm), [pycarl](https://github.com/moves-rwth/pycarl/) and [stormpy](https://github.com/moves-rwth/stormpy/) is preferred.
1. Update `CHANGELOG.md`
To get all the commits from an author since the last tag execute:
```console
git log last_tag..HEAD --author "author_name"
```
2. Update used carl version:
* Update `GIT_TAG` in `resources/3rdparty/carl/CMakeLists.txt`
* Maybe update `CARL_MINVERSION` in `resources/3rdparty/CMakeLists.txt`
3. Check that storm builds without errors and all tests are successful
* [Travis](https://travis-ci.org/moves-rwth/storm) should run successfully
4. Set new storm version:
* Set new storm version in `version.cmake`
5. Set new tag in git
```console
git tag -a new_version
git push origin new_version
```
Next we push the tag to GitHub. This step requires the GitHub repo to to be configured as a remote.
```console
git remote add github https://github.com/moves-rwth/storm.git
git push github new_version
```
The new tag should now be visible on [GitHub](https://github.com/moves-rwth/storm/tags)
6. [Add new release](https://github.com/moves-rwth/storm/releases/new) in GitHub
7. Update [Homebrew formula](https://github.com/moves-rwth/homebrew-storm)
8. Announce new storm version on [website](http://www.stormchecker.org/news.html)

7
resources/3rdparty/CMakeLists.txt

@ -208,7 +208,7 @@ include(${STORM_3RDPARTY_SOURCE_DIR}/include_cudd.cmake)
#############################################################
set(STORM_HAVE_CARL OFF)
set(CARL_MINVERSION "17.08")
set(CARL_MINVERSION "17.10")
if (NOT STORM_FORCE_SHIPPED_CARL)
if (NOT "${STORM_CARL_DIR_HINT}" STREQUAL "")
find_package(carl QUIET PATHS ${STORM_CARL_DIR_HINT} NO_DEFAULT_PATH)
@ -280,7 +280,8 @@ else()
add_dependencies(resources carl)
set(carl_INCLUDE_DIR "${STORM_3RDPARTY_BINARY_DIR}/carl/include/")
set(carl_LIBRARIES ${STORM_3RDPARTY_BINARY_DIR}/carl/lib/libcarl${DYNAMIC_EXT})
set(carl_DIR "${STORM_3RDPARTY_BINARY_DIR}/carl/")
set(carl_LIBRARIES ${STORM_3RDPARTY_BINARY_DIR}/carl/lib/libcarl${DYNAMIC_EXT})
set(STORM_HAVE_CARL ON)
message(STATUS "Storm - Linking with shipped carl ${carl_VERSION} (include: ${carl_INCLUDE_DIR}, library ${carl_LIBRARIES}, CARL_USE_CLN_NUMBERS: ${CARL_USE_CLN_NUMBERS}, CARL_USE_GINAC: ${CARL_USE_GINAC}).")
@ -418,7 +419,7 @@ ExternalProject_Add(
DOWNLOAD_COMMAND ""
PREFIX "sylvan"
SOURCE_DIR ${STORM_3RDPARTY_SOURCE_DIR}/sylvan
CMAKE_ARGS -DCMAKE_C_COMPILER=${CMAKE_C_COMPILER} -DGMP_LOCATION=${GMP_LIB_LOCATION} -DGMP_INCLUDE=${GMP_INCLUDE_DIR} -DCMAKE_CXX_COMPILER=${CMAKE_CXX_COMPILER} -DSYLVAN_BUILD_DOCS=OFF -DSYLVAN_BUILD_EXAMPLES=OFF -DCMAKE_BUILD_TYPE=${SYLVAN_BUILD_TYPE} -DCMAKE_POSITION_INDEPENDENT_CODE=ON -DUSE_CARL=ON -Dcarl_INCLUDE_DIR=${carl_INCLUDE_DIR} -DSYLVAN_PORTABLE=${STORM_PORTABLE} -Dcarl_LIBRARIES=${carl_LIBRARIES} -DBoost_INCLUDE_DIRS=${Boost_INCLUDE_DIRS} -DBUILD_SHARED_LIBS=OFF -DSYLVAN_BUILD_TESTS=OFF
CMAKE_ARGS -DCMAKE_C_COMPILER=${CMAKE_C_COMPILER} -DGMP_LOCATION=${GMP_LIB_LOCATION} -DGMP_INCLUDE=${GMP_INCLUDE_DIR} -DCMAKE_CXX_COMPILER=${CMAKE_CXX_COMPILER} -DSYLVAN_BUILD_DOCS=OFF -DSYLVAN_BUILD_EXAMPLES=OFF -DCMAKE_BUILD_TYPE=${SYLVAN_BUILD_TYPE} -DCMAKE_POSITION_INDEPENDENT_CODE=ON -DUSE_CARL=ON -Dcarl_DIR=${carl_DIR} -DSYLVAN_PORTABLE=${STORM_PORTABLE} -DBoost_INCLUDE_DIRS=${Boost_INCLUDE_DIRS} -DBUILD_SHARED_LIBS=OFF -DSYLVAN_BUILD_TESTS=OFF
BINARY_DIR ${STORM_3RDPARTY_BINARY_DIR}/sylvan
BUILD_IN_SOURCE 0
INSTALL_COMMAND ""

2
resources/3rdparty/carl/CMakeLists.txt

@ -8,7 +8,7 @@ message(STORM_3RDPARTY_BINARY_DIR: ${STORM_3RDPARTY_BINARY_DIR})
ExternalProject_Add(carl-config
GIT_REPOSITORY https://github.com/smtrat/carl
GIT_TAG 17.10
GIT_TAG 17.12
PREFIX here
SOURCE_DIR source_dir
BINARY_DIR ${STORM_3RDPARTY_BINARY_DIR}/carl

2
resources/3rdparty/sylvan/CMakeLists.txt

@ -61,8 +61,8 @@ include_directories("${Boost_INCLUDE_DIRS}")
if(USE_CARL)
add_definitions(-DSYLVAN_HAVE_CARL)
include_directories("${carl_INCLUDE_DIR}")
message(STATUS "Sylvan - Using CArL.")
find_package(carl REQUIRED PATHS ${carl_DIR} NO_DEFAULT_PATH)
else()
message(STATUS "Sylvan - Not using CArL.")
endif()

2
src/storm/adapters/Smt2ExpressionAdapter.h

@ -133,7 +133,7 @@ namespace storm {
STORM_LOG_DEBUG("Declaring the variable " + variableString);
declaredVariables.back().insert(variableString);
std::string varDeclaration = "( declare-fun " + variableString + " () ";
switch (variableToCheck.getType()){
switch (variableToCheck.type()){
case carl::VariableType::VT_BOOL:
varDeclaration += "Bool";
break;

159
src/storm/counterexamples/SMTMinimalLabelSetGenerator.h

@ -17,6 +17,9 @@
#include "storm/settings/SettingsManager.h"
#include "storm/settings/modules/CoreSettings.h"
#include "storm/utility/macros.h"
#include "storm/exceptions/NotSupportedException.h"
#include "storm/utility/counterexamples.h"
#include "storm/utility/cli.h"
@ -596,28 +599,6 @@ namespace storm {
// cuts.
std::unique_ptr<storm::solver::SmtSolver> localSolver(new storm::solver::Z3SmtSolver(program.getManager()));
storm::expressions::ExpressionManager const& localManager = program.getManager();
//
// // Create a context and register all variables of the program with their correct type.
// z3::context localContext;
// z3::solver localSolver(localContext);
// std::map<std::string, z3::expr> solverVariables;
// for (auto const& booleanVariable : program.getGlobalBooleanVariables()) {
// solverVariables.emplace(booleanVariable.getName(), localContext.bool_const(booleanVariable.getName().c_str()));
// }
// for (auto const& integerVariable : program.getGlobalIntegerVariables()) {
// solverVariables.emplace(integerVariable.getName(), localContext.int_const(integerVariable.getName().c_str()));
// }
//
// for (auto const& module : program.getModules()) {
// for (auto const& booleanVariable : module.getBooleanVariables()) {
// solverVariables.emplace(booleanVariable.getName(), localContext.bool_const(booleanVariable.getName().c_str()));
// }
// for (auto const& integerVariable : module.getIntegerVariables()) {
// solverVariables.emplace(integerVariable.getName(), localContext.int_const(integerVariable.getName().c_str()));
// }
// }
//
// storm::adapters::Z3ExpressionAdapter expressionAdapter(localContext, false, solverVariables);
// Then add the constraints for bounds of the integer variables..
for (auto const& integerVariable : program.getGlobalIntegerVariables()) {
@ -1592,9 +1573,7 @@ namespace storm {
* Returns the submdp obtained from removing all choices that do not originate from the specified filterLabelSet.
* Also returns the Labelsets of the submdp
*/
static std::pair<storm::models::sparse::Mdp<T>, std::vector<boost::container::flat_set<uint_fast64_t>>> restrictMdpToLabelSet(storm::models::sparse::Mdp<T> const& mdp, std::vector<boost::container::flat_set<uint_fast64_t>> const& labelSets, boost::container::flat_set<uint_fast64_t> const& filterLabelSet) {
STORM_LOG_THROW(mdp.getNumberOfChoices() == labelSets.size(), storm::exceptions::InvalidArgumentException, "The given number of labels does not match the number of choices.");
static std::pair<storm::models::sparse::Mdp<T>, std::vector<boost::container::flat_set<uint_fast64_t>>> restrictMdpToLabelSet(storm::models::sparse::Mdp<T> const& mdp, boost::container::flat_set<uint_fast64_t> const& filterLabelSet) {
std::vector<boost::container::flat_set<uint_fast64_t>> resultLabelSet;
storm::storage::SparseMatrixBuilder<T> transitionMatrixBuilder(0, mdp.getTransitionMatrix().getColumnCount(), 0, true, true, mdp.getTransitionMatrix().getRowGroupCount());
@ -1604,7 +1583,8 @@ namespace storm {
for(uint_fast64_t state = 0; state < mdp.getNumberOfStates(); ++state) {
bool stateHasValidChoice = false;
for (uint_fast64_t choice = mdp.getTransitionMatrix().getRowGroupIndices()[state]; choice < mdp.getTransitionMatrix().getRowGroupIndices()[state + 1]; ++choice) {
bool choiceValid = std::includes(filterLabelSet.begin(), filterLabelSet.end(), labelSets[choice].begin(), labelSets[choice].end());
auto const& choiceLabelSet = mdp.getChoiceOrigins()->asPrismChoiceOrigins().getCommandSet(choice);
bool choiceValid = std::includes(filterLabelSet.begin(), filterLabelSet.end(), choiceLabelSet.begin(), choiceLabelSet.end());
// If the choice is valid, copy over all its elements.
if (choiceValid) {
@ -1615,7 +1595,7 @@ namespace storm {
for (auto const& entry : mdp.getTransitionMatrix().getRow(choice)) {
transitionMatrixBuilder.addNextValue(currentRow, entry.getColumn(), entry.getValue());
}
resultLabelSet.push_back(labelSets[choice]);
resultLabelSet.push_back(choiceLabelSet);
++currentRow;
}
}
@ -1636,8 +1616,6 @@ namespace storm {
}
public:
/*!
* Computes the minimal command set that is needed in the given MDP to exceed the given probability threshold for satisfying phi until psi.
*
@ -1645,9 +1623,8 @@ namespace storm {
* @param mdp The MDP in which to find the minimal command set.
* @param phiStates A bit vector characterizing all phi states in the model.
* @param psiStates A bit vector characterizing all psi states in the model.
* @param probabilityThreshold The probability value that must be achieved or exceeded.
* @param strictBound A flag indicating whether the probability must be achieved (in which case the flag must be set) or strictly exceeded
* (if the flag is set to false).
* @param probabilityThreshold The threshold that is to be achieved or exceeded.
* @param strictBound Indicates whether the threshold needs to be achieved (true) or exceeded (false).
* @param checkThresholdFeasible If set, it is verified that the model can actually achieve/exceed the given probability value. If this check
* is made and fails, an exception is thrown.
*/
@ -1672,7 +1649,7 @@ namespace storm {
// (0) Obtain the label sets for each choice.
// The label set of a choice corresponds to the set of prism commands that induce the choice.
STORM_LOG_THROW(mdp.hasChoiceOrigins(), storm::exceptions::InvalidArgumentException, "Restriction to minimal command set is impossible for model without choice origns.");
STORM_LOG_THROW(mdp.hasChoiceOrigins(), storm::exceptions::InvalidArgumentException, "Restriction to minimal command set is impossible for model without choice origins.");
STORM_LOG_THROW(mdp.getChoiceOrigins()->isPrismChoiceOrigins(), storm::exceptions::InvalidArgumentException, "Restriction to command set is impossible for model without prism choice origins.");
storm::storage::sparse::PrismChoiceOrigins const& choiceOrigins = mdp.getChoiceOrigins()->asPrismChoiceOrigins();
std::vector<boost::container::flat_set<uint_fast64_t>> labelSets;
@ -1698,8 +1675,8 @@ namespace storm {
RelevancyInformation relevancyInformation = determineRelevantStatesAndLabels(mdp, labelSets, phiStates, psiStates);
// (3) Create a solver.
std::shared_ptr<storm::expressions::ExpressionManager> manager(new storm::expressions::ExpressionManager());
std::unique_ptr<storm::solver::SmtSolver> solver(new storm::solver::Z3SmtSolver(*manager));
std::shared_ptr<storm::expressions::ExpressionManager> manager = std::make_shared<storm::expressions::ExpressionManager>();
std::unique_ptr<storm::solver::SmtSolver> solver = std::make_unique<storm::solver::Z3SmtSolver>(*manager);
// (4) Create the variables for the relevant commands.
VariableInformation variableInformation = createVariables(manager, mdp, psiStates, relevancyInformation, includeReachabilityEncoding);
@ -1747,7 +1724,7 @@ namespace storm {
// Restrict the given MDP to the current set of labels and compute the reachability probability.
modelCheckingClock = std::chrono::high_resolution_clock::now();
commandSet.insert(relevancyInformation.knownLabels.begin(), relevancyInformation.knownLabels.end());
auto subMdpChoiceOrigins = restrictMdpToLabelSet(mdp, labelSets, commandSet);
auto subMdpChoiceOrigins = restrictMdpToLabelSet(mdp, commandSet);
storm::models::sparse::Mdp<T> const& subMdp = subMdpChoiceOrigins.first;
std::vector<boost::container::flat_set<uint_fast64_t>> const& subLabelSets = subMdpChoiceOrigins.second;
@ -1811,25 +1788,94 @@ namespace storm {
#endif
}
static std::shared_ptr<PrismHighLevelCounterexample> computeCounterexample(Environment const& env,storm::prism::Program program, storm::models::sparse::Mdp<T> const& mdp, std::shared_ptr<storm::logic::Formula const> const& formula) {
static void extendCommandSetLowerBound(storm::models::sparse::Mdp<T> const& mdp, boost::container::flat_set<uint_fast64_t>& commandSet, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates) {
auto startTime = std::chrono::high_resolution_clock::now();
// Create sub-MDP that only contains the choices allowed by the given command set.
storm::models::sparse::Mdp<T> subMdp = restrictMdpToLabelSet(mdp, commandSet).first;
// Then determine all prob0E(psi) states that are reachable in the sub-MDP.
storm::storage::BitVector reachableProb0EStates = storm::utility::graph::getReachableStates(subMdp.getTransitionMatrix(), subMdp.getInitialStates(), phiStates, psiStates);
// Create a queue of reachable prob0E(psi) states so we can check which commands need to be added
// to give them a strategy that avoids psi states.
std::queue<uint_fast64_t> prob0EWorklist;
for (auto const& e : reachableProb0EStates) {
prob0EWorklist.push(e);
}
// As long as there are reachable prob0E(psi) states, we add commands so they can stay within
// prob0E(states).
while (!prob0EWorklist.empty()) {
uint_fast64_t state = prob0EWorklist.front();
prob0EWorklist.pop();
// Now iterate over the original choices of the prob0E(psi) state and add at least one.
bool hasLabeledChoice = false;
uint64_t smallestCommandSetSize = 0;
uint64_t smallestCommandChoice = mdp.getTransitionMatrix().getRowGroupIndices()[state];
// Determine the choice with the least amount of commands (bad heuristic).
for (uint64_t choice = smallestCommandChoice; choice < mdp.getTransitionMatrix().getRowGroupIndices()[state + 1]; ++choice) {
bool onlyProb0ESuccessors = true;
for (auto const& successorEntry : mdp.getTransitionMatrix().getRow(choice)) {
if (!psiStates.get(successorEntry.getColumn())) {
onlyProb0ESuccessors = false;
break;
}
}
if (onlyProb0ESuccessors) {
auto const& labelSet = mdp.getChoiceOrigins()->asPrismChoiceOrigins().getCommandSet(choice);
hasLabeledChoice |= !labelSet.empty();
if (smallestCommandChoice == 0 || labelSet.size() < smallestCommandSetSize) {
smallestCommandSetSize = labelSet.size();
smallestCommandChoice = choice;
}
}
}
if (hasLabeledChoice) {
// Take all labels of the selected choice.
auto const& labelSet = mdp.getChoiceOrigins()->asPrismChoiceOrigins().getCommandSet(smallestCommandChoice);
commandSet.insert(labelSet.begin(), labelSet.end());
// Check for which successor states choices need to be added
for (auto const& successorEntry : mdp.getTransitionMatrix().getRow(smallestCommandChoice)) {
if (!storm::utility::isZero(successorEntry.getValue())) {
if (!reachableProb0EStates.get(successorEntry.getColumn())) {
reachableProb0EStates.set(successorEntry.getColumn());
prob0EWorklist.push(successorEntry.getColumn());
}
}
}
}
}
auto endTime = std::chrono::high_resolution_clock::now();
std::cout << std::endl << "Extended command for lower bounded property to size " << commandSet.size() << " in " << std::chrono::duration_cast<std::chrono::milliseconds>(endTime - startTime).count() << "ms." << std::endl;
}
static std::shared_ptr<PrismHighLevelCounterexample> computeCounterexample(Environment const& env, storm::prism::Program program, storm::models::sparse::Mdp<T> const& mdp, std::shared_ptr<storm::logic::Formula const> const& formula) {
#ifdef STORM_HAVE_Z3
std::cout << std::endl << "Generating minimal label counterexample for formula " << *formula << std::endl;
STORM_LOG_THROW(formula->isProbabilityOperatorFormula(), storm::exceptions::InvalidPropertyException, "Counterexample generation does not support this kind of formula. Expecting a probability operator as the outermost formula element.");
storm::logic::ProbabilityOperatorFormula const& probabilityOperator = formula->asProbabilityOperatorFormula();
STORM_LOG_THROW(probabilityOperator.hasBound(), storm::exceptions::InvalidPropertyException, "Counterexample generation only supports bounded formulas.");
storm::logic::ComparisonType comparisonType = probabilityOperator.getComparisonType();
STORM_LOG_THROW(comparisonType == storm::logic::ComparisonType::Less || comparisonType == storm::logic::ComparisonType::LessEqual, storm::exceptions::InvalidPropertyException, "Counterexample generation only supports formulas with an upper probability bound.");
STORM_LOG_THROW(probabilityOperator.getSubformula().isUntilFormula() || probabilityOperator.getSubformula().isEventuallyFormula(), storm::exceptions::InvalidPropertyException, "Path formula is required to be of the form 'phi U psi' for counterexample generation.");
storm::logic::ComparisonType comparisonType = probabilityOperator.getComparisonType();
bool strictBound = comparisonType == storm::logic::ComparisonType::Less;
double threshold = probabilityOperator.getThresholdAs<double>();
double threshold = probabilityOperator.getThresholdAs<T>();
storm::storage::BitVector phiStates;
storm::storage::BitVector psiStates;
storm::modelchecker::SparseMdpPrctlModelChecker<storm::models::sparse::Mdp<T>> modelchecker(mdp);
if (probabilityOperator.getSubformula().isUntilFormula()) {
STORM_LOG_THROW(!storm::logic::isLowerBound(comparisonType), storm::exceptions::NotSupportedException, "Lower bounds in counterexamples are only supported for eventually formulas.");
storm::logic::UntilFormula const& untilFormula = probabilityOperator.getSubformula().asUntilFormula();
std::unique_ptr<storm::modelchecker::CheckResult> leftResult = modelchecker.check(env, untilFormula.getLeftSubformula());
@ -1851,12 +1897,43 @@ namespace storm {
psiStates = subQualitativeResult.getTruthValuesVector();
}
bool lowerBoundedFormula = false;
if (storm::logic::isLowerBound(comparisonType)) {
// If the formula specifies a lower bound, we need to modify the phi and psi states.
// More concretely, we convert P(min)>lambda(F psi) to P(max)<(1-lambda)(G !psi) = P(max)<(1-lambda)(!psi U prob0E(psi))
// where prob0E(psi) is the set of states for which there exists a strategy \sigma_0 that avoids
// reaching psi states completely.
// This means that from all states in prob0E(psi) we need to include labels such that \sigma_0
// is actually included in the resulting MDP. This prevents us from guaranteeing the minimality of
// the returned counterexample, so we warn about that.
STORM_LOG_WARN("Generating counterexample for lower-bounded property. The resulting command set need not be minimal.");
// Modify bound appropriately.
comparisonType = storm::logic::invertPreserveStrictness(comparisonType);
threshold = storm::utility::one<T>() - threshold;
// Modify the phi and psi states appropriately.
storm::storage::BitVector statesWithProbability0E = storm::utility::graph::performProb0E(mdp.getTransitionMatrix(), mdp.getTransitionMatrix().getRowGroupIndices(), mdp.getBackwardTransitions(), phiStates, psiStates);
phiStates = ~psiStates;
psiStates = std::move(statesWithProbability0E);
// Remember our transformation so we can add commands to guarantee that the prob0E(a) states actually
// have a strategy that voids a states.
lowerBoundedFormula = true;
}
// Delegate the actual computation work to the function of equal name.
auto startTime = std::chrono::high_resolution_clock::now();
auto commandSet = getMinimalCommandSet(env, program, mdp, phiStates, psiStates, threshold, strictBound, true, storm::settings::getModule<storm::settings::modules::CounterexampleGeneratorSettings>().isEncodeReachabilitySet());
auto endTime = std::chrono::high_resolution_clock::now();
std::cout << std::endl << "Computed minimal command set of size " << commandSet.size() << " in " << std::chrono::duration_cast<std::chrono::milliseconds>(endTime - startTime).count() << "ms." << std::endl;
// Extend the command set properly.
if (lowerBoundedFormula) {
extendCommandSetLowerBound(mdp, commandSet, phiStates, psiStates);
}
return std::make_shared<PrismHighLevelCounterexample>(program.restrictCommands(commandSet));
#else

57
src/storm/logic/ComparisonType.h

@ -4,32 +4,45 @@
#include <iostream>
namespace storm {
namespace logic {
enum class ComparisonType { Less, LessEqual, Greater, GreaterEqual };
inline bool isStrict(ComparisonType t) {
return (t == ComparisonType::Less || t == ComparisonType::Greater);
}
inline bool isLowerBound(ComparisonType t) {
return (t == ComparisonType::Greater || t == ComparisonType::GreaterEqual);
namespace logic {
enum class ComparisonType { Less, LessEqual, Greater, GreaterEqual };
inline bool isStrict(ComparisonType t) {
return (t == ComparisonType::Less || t == ComparisonType::Greater);
}
inline bool isLowerBound(ComparisonType t) {
return (t == ComparisonType::Greater || t == ComparisonType::GreaterEqual);
}
inline ComparisonType invert(ComparisonType t) {
switch(t) {
case ComparisonType::Less:
return ComparisonType::GreaterEqual;
case ComparisonType::LessEqual:
return ComparisonType::Greater;
case ComparisonType::Greater:
return ComparisonType::LessEqual;
case ComparisonType::GreaterEqual:
return ComparisonType::Less;
}
}
inline ComparisonType invert(ComparisonType t) {
switch(t) {
case ComparisonType::Less:
return ComparisonType::GreaterEqual;
case ComparisonType::LessEqual:
return ComparisonType::Greater;
case ComparisonType::Greater:
return ComparisonType::LessEqual;
case ComparisonType::GreaterEqual:
return ComparisonType::Less;
}
inline ComparisonType invertPreserveStrictness(ComparisonType t) {
switch(t) {
case ComparisonType::Less:
return ComparisonType::Greater;
case ComparisonType::LessEqual:
return ComparisonType::GreaterEqual;
case ComparisonType::Greater:
return ComparisonType::Less;
case ComparisonType::GreaterEqual:
return ComparisonType::LessEqual;
}
}
std::ostream& operator<<(std::ostream& out, ComparisonType const& comparisonType);
}
std::ostream& operator<<(std::ostream& out, ComparisonType const& comparisonType);
}
}
#endif /* STORM_LOGIC_COMPARISONTYPE_H_ */

24
src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.cpp

@ -632,19 +632,21 @@ namespace storm {
}
// Use Fox-Glynn to get the truncation points and the weights.
std::tuple<uint_fast64_t, uint_fast64_t, ValueType, std::vector<ValueType>> foxGlynnResult = storm::utility::numerical::getFoxGlynnCutoff(lambda, 1e+300, storm::settings::getModule<storm::settings::modules::GeneralSettings>().getPrecision() / 8.0);
STORM_LOG_DEBUG("Fox-Glynn cutoff points: left=" << std::get<0>(foxGlynnResult) << ", right=" << std::get<1>(foxGlynnResult));
// std::tuple<uint_fast64_t, uint_fast64_t, ValueType, std::vector<ValueType>> foxGlynnResult = storm::utility::numerical::getFoxGlynnCutoff(lambda, 1e+300, storm::settings::getModule<storm::settings::modules::GeneralSettings>().getPrecision() / 8.0);
storm::utility::numerical::FoxGlynnResult<ValueType> foxGlynnResult = storm::utility::numerical::foxGlynn(lambda, storm::settings::getModule<storm::settings::modules::GeneralSettings>().getPrecision() / 8.0);
STORM_LOG_DEBUG("Fox-Glynn cutoff points: left=" << foxGlynnResult.left << ", right=" << foxGlynnResult.right);
// Scale the weights so they add up to one.
for (auto& element : std::get<3>(foxGlynnResult)) {
element /= std::get<2>(foxGlynnResult);
for (auto& element : foxGlynnResult.weights) {
element /= foxGlynnResult.totalWeight;
}
// If the cumulative reward is to be computed, we need to adjust the weights.
if (useMixedPoissonProbabilities) {
ValueType sum = storm::utility::zero<ValueType>();
for (auto& element : std::get<3>(foxGlynnResult)) {
for (auto& element : foxGlynnResult.weights) {
sum += element;
element = (1 - sum) / uniformizationRate;
}
@ -654,10 +656,10 @@ namespace storm {
// Initialize result.
std::vector<ValueType> result;
uint_fast64_t startingIteration = std::get<0>(foxGlynnResult);
uint_fast64_t startingIteration = foxGlynnResult.left;
if (startingIteration == 0) {
result = values;
storm::utility::vector::scaleVectorInPlace(result, std::get<3>(foxGlynnResult)[0]);
storm::utility::vector::scaleVectorInPlace(result, foxGlynnResult.weights.front());
++startingIteration;
} else {
if (useMixedPoissonProbabilities) {
@ -672,9 +674,9 @@ namespace storm {
std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(env, std::move(uniformizedMatrix), storm::solver::LinearEquationSolverTask::Multiply);
solver->setCachingEnabled(true);
if (!useMixedPoissonProbabilities && std::get<0>(foxGlynnResult) > 1) {
if (!useMixedPoissonProbabilities && foxGlynnResult.left > 1) {
// Perform the matrix-vector multiplications (without adding).
solver->repeatedMultiply(values, addVector, std::get<0>(foxGlynnResult) - 1);
solver->repeatedMultiply(values, addVector, foxGlynnResult.left - 1);
} else if (useMixedPoissonProbabilities) {
std::function<ValueType(ValueType const&, ValueType const&)> addAndScale = [&uniformizationRate] (ValueType const& a, ValueType const& b) { return a + b / uniformizationRate; };
@ -689,10 +691,10 @@ namespace storm {
// multiplication, scale and add the result.
ValueType weight = 0;
std::function<ValueType(ValueType const&, ValueType const&)> addAndScale = [&weight] (ValueType const& a, ValueType const& b) { return a + weight * b; };
for (uint_fast64_t index = startingIteration; index <= std::get<1>(foxGlynnResult); ++index) {
for (uint_fast64_t index = startingIteration; index <= foxGlynnResult.right; ++index) {
solver->repeatedMultiply(values, addVector, 1);
weight = std::get<3>(foxGlynnResult)[index - std::get<0>(foxGlynnResult)];
weight = foxGlynnResult.weights[index - foxGlynnResult.left];
storm::utility::vector::applyPointwise(result, values, result, addAndScale);
}

30
src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp

@ -463,6 +463,22 @@ namespace storm {
return result;
}
template<typename ValueType, typename RewardModelType>
storm::storage::BitVector SparseDtmcPrctlHelper<ValueType, RewardModelType>::BaierTransformedModel::getNewRelevantStates() const {
storm::storage::BitVector newRelevantStates(transitionMatrix.get().getRowCount());
for (uint64_t i = 0; i < this->beforeStates.getNumberOfSetBits(); ++i) {
newRelevantStates.set(i);
}
return newRelevantStates;
}
template<typename ValueType, typename RewardModelType>
storm::storage::BitVector SparseDtmcPrctlHelper<ValueType, RewardModelType>::BaierTransformedModel::getNewRelevantStates(storm::storage::BitVector const& oldRelevantStates) const {
storm::storage::BitVector result = oldRelevantStates % this->beforeStates;
result.resize(transitionMatrix.get().getRowCount());
return result;
}
template<typename ValueType, typename RewardModelType>
std::vector<ValueType> SparseDtmcPrctlHelper<ValueType, RewardModelType>::computeConditionalProbabilities(Environment const& env, storm::solver::SolveGoal<ValueType>&& goal, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& targetStates, storm::storage::BitVector const& conditionStates, bool qualitative, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
@ -482,6 +498,13 @@ namespace storm {
// Now compute reachability probabilities in the transformed model.
storm::storage::SparseMatrix<ValueType> const& newTransitionMatrix = transformedModel.transitionMatrix.get();
storm::storage::BitVector newRelevantValues;
if (goal.hasRelevantValues()) {
newRelevantValues = transformedModel.getNewRelevantStates(goal.relevantValues());
} else {
newRelevantValues = transformedModel.getNewRelevantStates();
}
goal.setRelevantValues(std::move(newRelevantValues));
std::vector<ValueType> conditionalProbabilities = computeUntilProbabilities(env, std::move(goal), newTransitionMatrix, newTransitionMatrix.transpose(), storm::storage::BitVector(newTransitionMatrix.getRowCount(), true), transformedModel.targetStates.get(), qualitative, linearEquationSolverFactory);
storm::utility::vector::setVectorValues(result, transformedModel.beforeStates, conditionalProbabilities);
@ -508,6 +531,13 @@ namespace storm {
// Now compute reachability probabilities in the transformed model.
storm::storage::SparseMatrix<ValueType> const& newTransitionMatrix = transformedModel.transitionMatrix.get();
storm::storage::BitVector newRelevantValues;
if (goal.hasRelevantValues()) {
newRelevantValues = transformedModel.getNewRelevantStates(goal.relevantValues());
} else {
newRelevantValues = transformedModel.getNewRelevantStates();
}
goal.setRelevantValues(std::move(newRelevantValues));
std::vector<ValueType> conditionalRewards = computeReachabilityRewards(env, std::move(goal), newTransitionMatrix, newTransitionMatrix.transpose(), transformedModel.stateRewards.get(), transformedModel.targetStates.get(), qualitative, linearEquationSolverFactory);
storm::utility::vector::setVectorValues(result, transformedModel.beforeStates, conditionalRewards);
}

5
src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.h

@ -59,7 +59,10 @@ namespace storm {
BaierTransformedModel() : noTargetStates(false) {
// Intentionally left empty.
}
storm::storage::BitVector getNewRelevantStates(storm::storage::BitVector const& oldRelevantStates) const;
storm::storage::BitVector getNewRelevantStates() const;
storm::storage::BitVector beforeStates;
boost::optional<storm::storage::SparseMatrix<ValueType>> transitionMatrix;
boost::optional<storm::storage::BitVector> targetStates;

5
src/storm/solver/AbstractEquationSolver.cpp

@ -139,6 +139,11 @@ namespace storm {
lowerBounds = values;
}
template<typename ValueType>
void AbstractEquationSolver<ValueType>::setLowerBounds(std::vector<ValueType>&& values) {
lowerBounds = std::move(values);
}
template<typename ValueType>
void AbstractEquationSolver<ValueType>::setUpperBounds(std::vector<ValueType> const& values) {
upperBounds = values;

5
src/storm/solver/AbstractEquationSolver.h

@ -117,6 +117,11 @@ namespace storm {
*/
void setLowerBounds(std::vector<ValueType> const& values);
/*!
* Sets lower bounds for the solution that can potentially be used by the solver.
*/
void setLowerBounds(std::vector<ValueType>&& values);
/*!
* Sets upper bounds for the solution that can potentially be used by the solver.
*/

4
src/storm/solver/NativeLinearEquationSolver.cpp

@ -937,9 +937,9 @@ namespace storm {
template<typename ValueType>
template<typename RationalType, typename ImpreciseType>
bool NativeLinearEquationSolver<ValueType>::sharpen(uint64_t precision, storm::storage::SparseMatrix<RationalType> const& A, std::vector<ImpreciseType> const& x, std::vector<RationalType> const& b, std::vector<RationalType>& tmp) {
for (uint64_t p = 0; p <= precision; ++p) {
for (uint64_t p = 1; p <= precision; ++p) {
storm::utility::kwek_mehlhorn::sharpen(p, x, tmp);
if (NativeLinearEquationSolver<RationalType>::isSolution(A, tmp, b)) {
return true;
}

2
src/storm/solver/SmtlibSmtSolver.cpp

@ -103,7 +103,7 @@ namespace storm {
}
void SmtlibSmtSolver::add(const storm::RationalFunctionVariable& variable, bool value){
STORM_LOG_THROW((variable.getType()==carl::VariableType::VT_BOOL), storm::exceptions::IllegalArgumentException, "Tried to add a constraint that consists of a non-boolean variable.");
STORM_LOG_THROW((variable.type()==carl::VariableType::VT_BOOL), storm::exceptions::IllegalArgumentException, "Tried to add a constraint that consists of a non-boolean variable.");
std::set<storm::RationalFunctionVariable> variableSet;
variableSet.insert(variable);
std::vector<std::string> const varDeclarations = expressionAdapter->checkForUndeclaredVariables(variableSet);

8
src/storm/solver/SymbolicMinMaxLinearEquationSolver.cpp

@ -142,9 +142,10 @@ namespace storm {
template<typename RationalType, typename ImpreciseType>
storm::dd::Add<DdType, RationalType> SymbolicMinMaxLinearEquationSolver<DdType, ValueType>::solveEquationsRationalSearchHelper(Environment const& env, OptimizationDirection dir, SymbolicMinMaxLinearEquationSolver<DdType, RationalType> const& rationalSolver, SymbolicMinMaxLinearEquationSolver<DdType, ImpreciseType> const& impreciseSolver, storm::dd::Add<DdType, RationalType> const& rationalB, storm::dd::Add<DdType, ImpreciseType> const& x, storm::dd::Add<DdType, ImpreciseType> const& b) const {
// Storage for the rational sharpened vector.
// Storage for the rational sharpened vector and the power iteration intermediate vector.
storm::dd::Add<DdType, RationalType> sharpenedX;
storm::dd::Add<DdType, ImpreciseType> currentX = x;
// The actual rational search.
uint64_t overallIterations = 0;
uint64_t valueIterationInvocations = 0;
@ -153,7 +154,7 @@ namespace storm {
bool relative = env.solver().minMax().getRelativeTerminationCriterion();
SolverStatus status = SolverStatus::InProgress;
while (status == SolverStatus::InProgress && overallIterations < maxIter) {
typename SymbolicMinMaxLinearEquationSolver<DdType, ImpreciseType>::ValueIterationResult viResult = impreciseSolver.performValueIteration(dir, x, b, storm::utility::convertNumber<ImpreciseType, ValueType>(precision), relative, maxIter);
typename SymbolicMinMaxLinearEquationSolver<DdType, ImpreciseType>::ValueIterationResult viResult = impreciseSolver.performValueIteration(dir, currentX, b, storm::utility::convertNumber<ImpreciseType, ValueType>(precision), relative, maxIter);
++valueIterationInvocations;
STORM_LOG_TRACE("Completed " << valueIterationInvocations << " value iteration invocations, the last one with precision " << precision << " completed in " << viResult.iterations << " iterations.");
@ -170,6 +171,7 @@ namespace storm {
if (isSolution) {
status = SolverStatus::Converged;
} else {
currentX = viResult.values;
precision /= storm::utility::convertNumber<ValueType, uint64_t>(10);
}
}

10
src/storm/solver/SymbolicNativeLinearEquationSolver.cpp

@ -174,10 +174,10 @@ namespace storm {
storm::dd::Add<DdType, RationalType> sharpenedX;
for (uint64_t p = 1; p < precision; ++p) {
for (uint64_t p = 1; p <= precision; ++p) {
sharpenedX = x.sharpenKwekMehlhorn(p);
isSolution = rationalSolver.isSolutionFixedPoint(sharpenedX, rationalB);
if (isSolution) {
break;
}
@ -190,8 +190,9 @@ namespace storm {
template<typename RationalType, typename ImpreciseType>
storm::dd::Add<DdType, RationalType> SymbolicNativeLinearEquationSolver<DdType, ValueType>::solveEquationsRationalSearchHelper(Environment const& env, SymbolicNativeLinearEquationSolver<DdType, RationalType> const& rationalSolver, SymbolicNativeLinearEquationSolver<DdType, ImpreciseType> const& impreciseSolver, storm::dd::Add<DdType, RationalType> const& rationalB, storm::dd::Add<DdType, ImpreciseType> const& x, storm::dd::Add<DdType, ImpreciseType> const& b) const {
// Storage for the rational sharpened vector.
// Storage for the rational sharpened vector and the power iteration intermediate vector.
storm::dd::Add<DdType, RationalType> sharpenedX;
storm::dd::Add<DdType, ImpreciseType> currentX = x;
// The actual rational search.
uint64_t overallIterations = 0;
@ -201,7 +202,7 @@ namespace storm {
bool relative = env.solver().native().getRelativeTerminationCriterion();
SolverStatus status = SolverStatus::InProgress;
while (status == SolverStatus::InProgress && overallIterations < maxIter) {
typename SymbolicNativeLinearEquationSolver<DdType, ImpreciseType>::PowerIterationResult result = impreciseSolver.performPowerIteration(x, b, storm::utility::convertNumber<ImpreciseType, ValueType>(precision), relative, maxIter - overallIterations);
typename SymbolicNativeLinearEquationSolver<DdType, ImpreciseType>::PowerIterationResult result = impreciseSolver.performPowerIteration(currentX, b, storm::utility::convertNumber<ImpreciseType, ValueType>(precision), relative, maxIter - overallIterations);
++powerIterationInvocations;
STORM_LOG_TRACE("Completed " << powerIterationInvocations << " power iteration invocations, the last one with precision " << precision << " completed in " << result.iterations << " iterations.");
@ -218,6 +219,7 @@ namespace storm {
if (isSolution) {
status = SolverStatus::Converged;
} else {
currentX = result.values;
precision /= storm::utility::convertNumber<ValueType, uint64_t>(10);
}
}

6
src/storm/storage/dd/Add.cpp

@ -941,7 +941,7 @@ namespace storm {
template<DdType LibraryType, typename ValueType>
std::ostream& operator<<(std::ostream& out, Add<LibraryType, ValueType> const& add) {
out << "ADD with " << add.getNonZeroCount() << " nnz, " << add.getNodeCount() << " nodes, " << add.getLeafCount() << " leaves" << std::endl;
out << "ADD [" << add.getInternalAdd().getStringId() << "] with " << add.getNonZeroCount() << " nnz, " << add.getNodeCount() << " nodes, " << add.getLeafCount() << " leaves" << std::endl;
std::vector<std::string> variableNames;
for (auto const& variable : add.getContainedMetaVariables()) {
variableNames.push_back(variable.getName());
@ -993,15 +993,19 @@ namespace storm {
}
template class Add<storm::dd::DdType::CUDD, double>;
template std::ostream& operator<<(std::ostream& out, Add<storm::dd::DdType::CUDD, double> const& add);
template class Add<storm::dd::DdType::CUDD, uint_fast64_t>;
template class Add<storm::dd::DdType::Sylvan, double>;
template std::ostream& operator<<(std::ostream& out, Add<storm::dd::DdType::Sylvan, double> const& add);
template class Add<storm::dd::DdType::Sylvan, uint_fast64_t>;
#ifdef STORM_HAVE_CARL
template class Add<storm::dd::DdType::CUDD, storm::RationalNumber>;
template std::ostream& operator<<(std::ostream& out, Add<storm::dd::DdType::CUDD, storm::RationalNumber> const& add);
template class Add<storm::dd::DdType::Sylvan, storm::RationalNumber>;
template std::ostream& operator<<(std::ostream& out, Add<storm::dd::DdType::Sylvan, storm::RationalNumber> const& add);
template class Add<storm::dd::DdType::Sylvan, storm::RationalFunction>;
template Add<storm::dd::DdType::CUDD, storm::RationalNumber> Add<storm::dd::DdType::CUDD, storm::RationalNumber>::toValueType<storm::RationalNumber>() const;

7
src/storm/storage/dd/cudd/InternalCuddAdd.cpp

@ -420,6 +420,13 @@ namespace storm {
DdNode* InternalAdd<DdType::CUDD, ValueType>::getCuddDdNode() const {
return this->getCuddAdd().getNode();
}
template<typename ValueType>
std::string InternalAdd<DdType::CUDD, ValueType>::getStringId() const {
std::stringstream ss;
ss << this->getCuddDdNode();
return ss.str();
}
template<typename ValueType>
Odd InternalAdd<DdType::CUDD, ValueType>::createOdd(std::vector<uint_fast64_t> const& ddVariableIndices) const {

5
src/storm/storage/dd/cudd/InternalCuddAdd.h

@ -629,6 +629,11 @@ namespace storm {
*/
DdNode* getCuddDdNode() const;
/*!
* Retrieves a string representation of an ID for thid ADD.
*/
std::string getStringId() const;
private:
/*!
* Performs a recursive step to perform the given function between the given DD-based vector and the given

7
src/storm/storage/dd/sylvan/InternalSylvanAdd.cpp

@ -1151,6 +1151,13 @@ namespace storm {
return sylvanMtbdd;
}
template<typename ValueType>
std::string InternalAdd<DdType::Sylvan, ValueType>::getStringId() const {
std::stringstream ss;
ss << this->getSylvanMtbdd().GetMTBDD();
return ss.str();
}
template class InternalAdd<DdType::Sylvan, double>;
template class InternalAdd<DdType::Sylvan, uint_fast64_t>;

2
src/storm/storage/dd/sylvan/InternalSylvanAdd.h

@ -620,6 +620,8 @@ namespace storm {
*/
static ValueType getValue(MTBDD const& node);
std::string getStringId() const;
private:
/*!
* Recursively builds the ODD from an ADD.

2
src/storm/storage/sparse/PrismChoiceOrigins.cpp

@ -111,4 +111,4 @@ namespace storm {
}
}
}
}
}

2
src/storm/storage/sparse/PrismChoiceOrigins.h

@ -67,4 +67,4 @@ namespace storm {
};
}
}
}
}

275
src/storm/utility/numerical.cpp

@ -0,0 +1,275 @@
#include "storm/utility/numerical.h"
#include <cmath>
#include <boost/math/constants/constants.hpp>
#include "storm/utility/macros.h"
#include "storm/utility/constants.h"
#include "storm/exceptions/InvalidArgumentException.h"
#include "storm/exceptions/PrecisionExceededException.h"
namespace storm {
namespace utility {
namespace numerical {
template<typename ValueType>
FoxGlynnResult<ValueType>::FoxGlynnResult() : left(0), right(0), totalWeight(storm::utility::zero<ValueType>()) {
// Intentionally left empty.
}
/*!
* The following implementation of Fox and Glynn's algorithm is taken from David Jansen's patched version
* in MRMC, which is based on his paper:
*
* https://pms.cs.ru.nl/iris-diglib/src/getContent.php?id=2011-Jansen-UnderstandingFoxGlynn
*
* We have only adapted the code to match more of C++'s and our coding guidelines.
*/
template<typename ValueType>
FoxGlynnResult<ValueType> foxGlynnFinder(ValueType lambda, ValueType epsilon) {
ValueType tau = std::numeric_limits<ValueType>::min();
ValueType omega = std::numeric_limits<ValueType>::max();
ValueType const sqrt_2_pi = boost::math::constants::root_two_pi<ValueType>();
ValueType const log10_e = std::log10(boost::math::constants::e<ValueType>());
uint64_t m = static_cast<uint64_t>(lambda);
int64_t left = 0;
int64_t right = 0;
// tau is only used in underflow checks, which we are going to do in the logarithm domain.
tau = log(tau);
// In error bound comparisons, we always compare with epsilon*sqrt_2_pi.
epsilon *= sqrt_2_pi;
// Compute left truncation point.
if (m < 25) {
// For lambda below 25 the exponential can be smaller than tau. If that is the case we expect
// underflows and warn the user.
if (-lambda <= tau) {
STORM_LOG_WARN("Fox-Glynn: 0 < lambda < 25, underflow near Poi(" << lambda << ", 0) = " << std::exp(-lambda) << ". The results are unreliable.");
}
// Zero is used as left truncation point for lambda <= 25.
left = 0;
} else {
// Compute the left truncation point for lambda >= 25 (for lambda < 25 we use zero as left truncation point).
ValueType const bl = (1 + 1 / lambda) * std::exp((1/lambda) * 0.125);
ValueType const sqrt_lambda = std::sqrt(lambda);
int64_t k;
// Start looking for the left truncation point:
// * start search at k=4 (taken from original Fox-Glynn paper)
// * increase the left truncation point until we fulfil the error condition
for (k = 4;; ++k) {
ValueType max_err;
left = m - static_cast<int64_t>(std::ceil(k*sqrt_lambda + 0.5));
// For small lambda the above calculation can yield negative truncation points, crop them here.
if (left <= 0) {
left = 0;
break;
}
// Note that Propositions 2-4 in Fox--Glynn mix up notation: they write Phi where they mean
// 1 - Phi. (In Corollaries 1 and 2, phi is used correctly again.)
max_err = bl * exp(-0.5 * (k*k)) / k;
if (max_err * 2 <= epsilon) {
// If the error on the left hand side is smaller, we can be more lenient on the right hand
// side. To this end, we now set epsilon to the part of the error that has not yet been eaten
// up by the left-hand truncation.
epsilon -= max_err;
break;
}
}
// Finally the left truncation point is found.
}
// Compute right truncation point.
{
ValueType lambda_max;
int64_t m_max, k;
// According to Fox-Glynn, if lambda < 400 we should take lambda = 400, otherwise use the original
// value. This is for computing the right truncation point.
if (m < 400) {
lambda_max = 400;
m_max = 400;
epsilon *= 0.662608824988162441697980;
/* i.e. al = (1+1/400) * exp(1/16) * sqrt_2; epsilon /= al; */
} else {
lambda_max = lambda;
m_max = m;
epsilon *= (1 - 1 / (lambda + 1)) * 0.664265347050632847802225;
/* i.e. al = (1+1/lambda) * exp(1/16) * sqrt_2; epsilon /= al; */
}
// Find right truncation point.
// This loop is a modification to the original Fox-Glynn paper.
// The search for the right truncation point is only terminated by the error condition and not by
// the stop index from the FG paper. This can yield more accurate results if necessary.
for (k = 4;; ++k) {
// dkl_inv is between 1 - 1e-33 and 1 if lambda_max >= 400 and k >= 4; this will always be
// rounded to 1.0. We therefore leave the factor out.
// double dkl_inv=1 - exp(-266/401.0 * (k*sqrt(2*lambda_max) + 1.5));
// actually: "k * (dkl_inv*epsilon/al) >= exp(-0.5 * k^2)", but epsilon has been changed appropriately.
if (k * epsilon >= exp(-0.5*(k*k))) {
break;
}
}
right = m_max + static_cast<int64_t>(std::ceil(k * std::sqrt(2 * lambda_max) + 0.5));
if (right > m_max + static_cast<int64_t>(std::ceil((lambda_max + 1) * 0.5))) {
STORM_LOG_WARN("Fox-Glynn: right = " << right << " >> lambda = " << lambda_max << ", cannot bound the right tail. The results are unreliable.");
}
}
// Time to set the initial value for weights.
FoxGlynnResult<ValueType> fgresult;
fgresult.left = static_cast<uint64_t>(left);
fgresult.right = static_cast<uint64_t>(right);
fgresult.weights.resize(fgresult.right - fgresult.left + 1);
fgresult.weights[m - left] = omega / (1.0e+10 * (right - left));
if (m >= 25) {
// Perform underflow check.
ValueType result, log_c_m_inf;
int64_t i;
// we are going to compare with tau - log(w[m]).
tau -= std::log(fgresult.weights[m - left]);
// We take the c_m_inf = 0.14627 / sqrt( m ), as for lambda >= 25
// c_m = 1 / ( sqrt( 2.0 * pi * m ) ) * exp( m - lambda - 1 / ( 12.0 * m ) ) => c_m_inf.
// Note that m-lambda is in the interval (-1,0], and -1/(12*m) is in [-1/(12*25),0).
// So, exp(m-lambda - 1/(12*m)) is in (exp(-1-1/(12*25)),exp(0)).
// Therefore, we can improve the lower bound on c_m to exp(-1-1/(12*25)) / sqrt(2*pi) = ~0.14627.
// Its logarithm is -1 - 1/(12*25) - log(2*pi) * 0.5 = ~ -1.922272 (rounded towards -infinity).
log_c_m_inf = -1.922272 - log((double) m) * 0.5;
// We use FG's Proposition 6 directly (and not Corollary 4 i and ii), as k_prime may be too large
// if pFG->left == 0.
i = m - left;
// Equivalent to 2*i <= m, equivalent to i <= lambda/2.
if (i <= left) {
// Use Proposition 6 (i). Note that Fox--Glynn are off by one in the proof of this proposition;
// they sum up to i-1, but should have summed up to i. */
result = log_c_m_inf
- i * (i+1) * (0.5 + (2*i+1)/(6*lambda)) / lambda;
} else {
// Use Corollary 4 (iii). Note that k_prime <= sqrt(m+1)/m is a misprint for k_prime <= m/sqrt(m+1),
// which is equivalent to left >= 0, which holds trivially.
result = -lambda;
if (left != 0) {
// Also use Proposition 6 (ii).
double result_1 = log_c_m_inf + i * log(1 - i/(double) (m+1));
// Take the maximum.
if (result_1 > result) {
result = result_1;
}
}
}
if (result <= tau) {
int64_t const log10_result = static_cast<int64_t>(std::floor(result * log10_e));
STORM_LOG_WARN("Fox-Glynn: lambda >= 25, underflow near Poi(" << lambda << "," << left << ") <= " << std::exp(result - log10_result/log10_e) << log10_result << ". The results are unreliable.");
}
// We still have to perform an underflow check for the right truncation point when lambda >= 400.
if (m >= 400) {
// Use Proposition 5 of Fox--Glynn.
i = right - m;
result = log_c_m_inf - i * (i + 1) / (2 * lambda);
if (result <= tau) {
int64_t const log10_result = static_cast<int64_t>(std::floor(result * log10_e));
STORM_LOG_WARN("Fox-Glynn: lambda >= 25, underflow near Poi(" << lambda << "," << right << ") <= " << std::exp(result - log10_result/log10_e) << log10_result << ". The results are unreliable.");
}
}
}
return fgresult;
}
template<typename ValueType>
FoxGlynnResult<ValueType> foxGlynnWeighter(ValueType lambda, ValueType epsilon) {
ValueType tau = std::numeric_limits<ValueType>::min();
// The magic m point.
uint64_t m = static_cast<uint64_t>(lambda);
int64_t j, t;
FoxGlynnResult<ValueType> result = foxGlynnFinder(lambda, epsilon);
// Fill the left side of the array.
for (j = m - result.left; j > 0; --j) {
result.weights[j - 1] = (j + result.left) / lambda * result.weights[j];
}
t = result.right - result.left;
// Fill the right side of the array, have two cases lambda < 400 & lambda >= 400.
if (m < 400) {
// Perform the underflow check, according to Fox-Glynn.
STORM_LOG_THROW(result.right <= 600, storm::exceptions::PrecisionExceededException, "Fox-Glynn: " << result.right << " > 600, underflow is possible.");
// Compute weights.
for (j = m - result.left; j < t; ++j) {
ValueType q = lambda / (j + 1 + result.left);
if (result.weights[j] > tau / q) {
result.weights[j + 1] = q * result.weights[j];
} else {
t = j;
result.right = j + result.left;
// It's time to compute W.
break;
}
}
} else {
// Compute weights.
for (j = m - result.left; j < t; ++j) {
result.weights[j + 1] = lambda / (j + 1 + result.left) * result.weights[j];
}
}
// It is time to compute the normalization weight W.
result.totalWeight = storm::utility::zero<ValueType>();
j = 0;
// t was set above.
while(j < t) {
if (result.weights[j] <= result.weights[t]) {
result.totalWeight += result.weights[j];
j++;
} else {
result.totalWeight += result.weights[t];
t--;
}
}
result.totalWeight += result.weights[j];
STORM_LOG_TRACE("Fox-Glynn: ltp = " << result.left << ", rtp = " << result.right << ", w = " << result.totalWeight << ".");
return result;
}
template<typename ValueType>
FoxGlynnResult<ValueType> foxGlynn(ValueType lambda, ValueType epsilon) {
STORM_LOG_THROW(lambda > 0, storm::exceptions::InvalidArgumentException, "Fox-Glynn requires positive lambda.");
return foxGlynnWeighter(lambda, epsilon);
}
template FoxGlynnResult<double> foxGlynn(double lambda, double epsilon);
}
}
}

134
src/storm/utility/numerical.h

@ -1,131 +1,23 @@
#include <vector>
#include <tuple>
#include <cmath>
#include "storm/utility/macros.h"
#include "storm/utility/constants.h"
#include "storm/exceptions/InvalidArgumentException.h"
#include "storm/exceptions/OutOfRangeException.h"
#include <cstdint>
namespace storm {
namespace utility {
namespace numerical {
template<typename ValueType>
struct FoxGlynnResult {
FoxGlynnResult();
uint64_t left;
uint64_t right;
ValueType totalWeight;
std::vector<ValueType> weights;
};
template<typename ValueType>
std::tuple<uint_fast64_t, uint_fast64_t, ValueType, std::vector<ValueType>> getFoxGlynnCutoff(ValueType lambda, ValueType overflow, ValueType accuracy) {
STORM_LOG_THROW(lambda != storm::utility::zero<ValueType>(), storm::exceptions::InvalidArgumentException, "Error in Fox-Glynn algorithm: lambda must not be zero.");
FoxGlynnResult<ValueType> foxGlynn(ValueType lambda, ValueType epsilon);
// This code is a modified version of the one in PRISM. According to their implementation, for lambda
// smaller than 400, we compute the result using the naive method.
if (lambda < 400) {
ValueType eToPowerMinusLambda = std::exp(-lambda);
ValueType targetValue = (1 - accuracy) / eToPowerMinusLambda;
std::vector<ValueType> weights;
ValueType exactlyKEvents = 1;
ValueType atMostKEvents = exactlyKEvents;
weights.push_back(exactlyKEvents * eToPowerMinusLambda);
uint_fast64_t k = 1;
do {
exactlyKEvents *= lambda / k;
atMostKEvents += exactlyKEvents;
weights.push_back(exactlyKEvents * eToPowerMinusLambda);
++k;
} while (atMostKEvents < targetValue);
return std::make_tuple(0, k - 1, 1.0, weights);
} else {
STORM_LOG_THROW(accuracy >= 1e-10, storm::exceptions::InvalidArgumentException, "Error in Fox-Glynn algorithm: the accuracy must not be below 1e-10.");
// Factor from Fox&Glynn's paper. The paper does not explain where it comes from.
ValueType factor = 1e+10;
// Now start the Finder algorithm to find the truncation points.
ValueType m = std::floor(lambda);
uint_fast64_t leftTruncationPoint = 0, rightTruncationPoint = 0;
{
// Factors used by the corollaries explained in Fox & Glynns paper.
// Square root of pi.
ValueType sqrtpi = 1.77245385090551602729;
// Square root of 2.
ValueType sqrt2 = 1.41421356237309504880;
// Set up a_\lambda, b_\lambda, and the square root of lambda.
ValueType aLambda = 0, bLambda = 0, sqrtLambda = 0;
if (m < 400) {
sqrtLambda = std::sqrt(400.0);
aLambda = (1.0 + 1.0 / 400.0) * std::exp(0.0625) * sqrt2;
bLambda = (1.0 + 1.0 / 400.0) * std::exp(0.125 / 400.0);
} else {
sqrtLambda = std::sqrt(lambda);
aLambda = (1.0 + 1.0 / lambda) * std::exp(0.0625) * sqrt2;
bLambda = (1.0 + 1.0 / lambda) * std::exp(0.125 / lambda);
}
// Use Corollary 1 from the paper to find the right truncation point.
uint_fast64_t k = 4;
ValueType dkl = 1.0 / (1 - std::exp(-(2.0 / 9.0) * (k * sqrt2 * sqrtLambda + 1.5)));
// According to David Jansen the upper bound can be ignored to achieve more accurate results.
// Right hand side of the equation in Corollary 1.
while ((accuracy / 2.0) < (aLambda * dkl * std::exp(-(k*k / 2.0)) / (k * sqrt2 * sqrtpi))) {
++k;
// d(k,Lambda) from the paper.
dkl = 1.0 / (1 - std::exp(-(2.0 / 9.0)*(k * sqrt2 * sqrtLambda + 1.5)));
}
// Left hand side of the equation in Corollary 1.
rightTruncationPoint = static_cast<uint_fast64_t>(std::ceil((m + k * sqrt2 * sqrtLambda + 1.5)));
// Use Corollary 2 to find left truncation point.
k = 4;
// Right hand side of the equation in Corollary 2.
while ((accuracy / 2.0) < ((bLambda * std::exp(-(k*k / 2.0))) / (k * sqrt2 * sqrtpi))) {
++k;
}
// Left hand side of the equation in Corollary 2.
leftTruncationPoint = static_cast<uint_fast64_t>(std::max(std::trunc(m - k * sqrtLambda - 1.5), storm::utility::zero<ValueType>()));
// Check for underflow.
STORM_LOG_THROW(std::trunc((m - k * sqrtLambda - 1.5)) >= 0, storm::exceptions::OutOfRangeException, "Error in Fox-Glynn algorithm: Underflow of left truncation point.");
}
std::vector<ValueType> weights(rightTruncationPoint - leftTruncationPoint + 1);
weights[m - leftTruncationPoint] = overflow / (factor * (rightTruncationPoint - leftTruncationPoint));
for (uint_fast64_t j = m; j > leftTruncationPoint; --j) {
weights[j - 1 - leftTruncationPoint] = (j / lambda) * weights[j - leftTruncationPoint];
}
for (uint_fast64_t j = m; j < rightTruncationPoint; ++j) {
weights[j + 1 - leftTruncationPoint] = (lambda / (j + 1)) * weights[j - leftTruncationPoint];
}
// Compute the total weight and start with smallest to avoid roundoff errors.
ValueType totalWeight = storm::utility::zero<ValueType>();
uint_fast64_t s = leftTruncationPoint;
uint_fast64_t t = rightTruncationPoint;
while (s < t) {
if (weights[s - leftTruncationPoint] <= weights[t - leftTruncationPoint]) {
totalWeight += weights[s - leftTruncationPoint];
++s;
} else {
totalWeight += weights[t - leftTruncationPoint];
--t;
}
}
totalWeight += weights[s - leftTruncationPoint];
return std::make_tuple(leftTruncationPoint, rightTruncationPoint, totalWeight, weights);
}
}
}
}
}

2
version.cmake

@ -1,4 +1,4 @@
set(STORM_VERSION_MAJOR 1)
set(STORM_VERSION_MINOR 1)
set(STORM_VERSION_MINOR 2)
set(STORM_VERSION_PATCH 0)
Loading…
Cancel
Save