Merge remote-tracking branch 'origin/prism-pomdp' into prism-pomdp

CHANGELOG.md

@@ -7,8 +7,12 @@ The releases of major and minor versions contain an overview of changes since th
 Version 1.4.x
 -------------
 
-### Version 1.4.1 (under development)
+## Version 1.4.2 (under development)
+- n.a.
 
+### Version 1.4.1 (2019/12)
+- Implemented long run average (LRA) computation for DTMCs/CTMCs via value iteration and via gain/bias equations.
+- Added several LRA related settings in a new settings module. Note that `--minmax:lramethod` has been replaced by `--lra:nondetmethod`.
 
 ### Version 1.4.0 (2019/11)
 - Added support for multi-dimensional quantile queries.
@@ -44,6 +48,7 @@ Version 1.4.x
 - Fixed compilation for macOS Mojave and higher.
 - Several bug fixes.
 
+
 Version 1.3.x
 -------------
 
@@ -79,6 +84,7 @@ Version 1.3.x
 - `storm-gspn` improved
 - Sound value iteration
 
+
 Version 1.2.x
 -------------
 

README.md

@@ -1,6 +1,9 @@
 Storm - A Modern Probabilistic Model Checker
 ============================================
 [![Build Status](https://travis-ci.org/moves-rwth/storm.svg?branch=master)](https://travis-ci.org/moves-rwth/storm)
+[![GitHub release](https://img.shields.io/github/release/moves-rwth/storm.svg)](https://github.com/moves-rwth/storm/releases/)
+[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.3574502.svg)](https://doi.org/10.5281/zenodo.3574502)
+
 
 For more instructions, check out the documentation found in [Getting Started](http://www.stormchecker.org/getting-started.html).
 

doc/checklist_new_release.md

@@ -1,11 +1,12 @@
 The following steps should be performed before releasing a new storm version.
 Note that in most cases 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:
+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"
    ```
+   * Set release month
 
 2. Update used carl version:
    * Update `GIT_TAG` in `resources/3rdparty/carl/CMakeLists.txt`
@@ -32,6 +33,15 @@ Note that in most cases a simultaneous release of [carl](https://github.com/smtr
 
 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).
+7. Update `stable` branch:
+   ```console
+   git checkout stable
+   git merge master
+   git push origin stable
+   ```
+
+8. Update [Homebrew formula](https://github.com/moves-rwth/homebrew-storm).
+
+9. Announce new storm version on [website](http://www.stormchecker.org/news.html).
 
-8. Announce new storm version on [website](http://www.stormchecker.org/news.html).
+10. Create [Docker containers](https://hub.docker.com/r/movesrwth/storm) for new version.

resources/examples/testfiles/ctmc/kanban.prism

@@ -0,0 +1,111 @@
+// kanban manufacturing example [CT96]
+// dxp/gxn 3/2/00
+
+ctmc
+
+// number of tokens 
+const int t = 1; 
+
+// rates
+const double in1      = 1.0;
+const double out4     = 0.9;
+const double synch123 = 0.4;
+const double synch234 = 0.5;
+const double back     = 0.3;
+const double redo1    = 0.36;
+const double redo2    = 0.42;
+const double redo3    = 0.39;
+const double redo4    = 0.33;
+const double ok1      = 0.84;
+const double ok2      = 0.98;
+const double ok3      = 0.91;
+const double ok4      = 0.77;
+
+module k1
+
+	w1 : [0..t];
+	x1 : [0..t];
+	y1 : [0..t];
+	z1 : [0..t];
+	
+	[in]   (w1<t) & (x1<t) -> in1 	 : (w1'=w1+1) & (x1'=x1+1);
+	[]   (x1>0) & (y1<t) -> redo1	 : (x1'=x1-1) & (y1'=y1+1);
+	[]   (x1>0) & (z1<t) -> ok1 	 : (x1'=x1-1) & (z1'=z1+1);
+	[]   (y1>0) & (x1<t) -> back	 : (y1'=y1-1) & (x1'=x1+1);
+	[s1] (z1>0) & (w1>0) -> synch123 : (z1'=z1-1) & (w1'=w1-1);
+
+endmodule
+
+module k2
+
+	w2 : [0..t];
+	x2 : [0..t];
+	y2 : [0..t];
+	z2 : [0..t];
+	
+	[s1] (w2<t) & (x2<t) -> 1	 : (w2'=w2+1) & (x2'=x2+1);
+	[]   (x2>0) & (y2<t) -> redo2 : (x2'=x2-1) & (y2'=y2+1);
+	[]   (x2>0) & (z2<t) -> ok2   : (x2'=x2-1) & (z2'=z2+1);
+	[]   (y2>0) & (x2<t) -> back  : (y2'=y2-1) & (x2'=x2+1);
+	[s2] (z2>0) & (w2>0) -> 1	 : (z2'=z2-1) & (w2'=w2-1);
+
+endmodule
+
+module k3
+
+	w3 : [0..t];
+	x3 : [0..t];
+	y3 : [0..t];
+	z3 : [0..t];
+	
+	[s1] (w3<t) & (x3<t) -> 1	 : (w3'=w3+1) & (x3'=x3+1);
+	[]   (x3>0) & (y3<t) -> redo3 : (x3'=x3-1) & (y3'=y3+1);
+	[]   (x3>0) & (z3<t) -> ok3   : (x3'=x3-1) & (z3'=z3+1);
+	[]   (y3>0) & (x3<t) -> back  : (y3'=y3-1) & (x3'=x3+1);
+	[s2] (z3>0) & (w3>0) -> 1	 : (z3'=z3-1) & (w3'=w3-1);
+
+
+endmodule
+
+module k4
+
+	w4 : [0..t];
+	x4 : [0..t];
+	y4 : [0..t];
+	z4 : [0..t];
+	
+	[s2] (w4<t) & (x4<t) -> synch234 : (w4'=w4+1) & (x4'=x4+1);
+	[]   (x4>0) & (y4<t) -> redo4	: (x4'=x4-1) & (y4'=y4+1);
+	[]   (x4>0) & (z4<t) -> ok4  	: (x4'=x4-1) & (z4'=z4+1);
+	[]   (y4>0) & (x4<t) -> back 	: (y4'=y4-1) & (x4'=x4+1);
+	[]   (z4>0) & (w4>0) -> out4 	: (z4'=z4-1) & (w4'=w4-1);
+
+endmodule
+
+// reward structures
+
+// tokens in cell1
+rewards "tokens_cell1"
+	true : x1+y1+z1;
+endrewards
+
+// tokens in cell2
+rewards "tokens_cell2"
+	true : x2+y2+z2;
+endrewards
+
+// tokens in cell3
+rewards "tokens_cell3"
+	true : x3+y3+z3;
+endrewards
+
+// tokens in cell4
+rewards "tokens_cell4"
+	true : x4+y4+z4;
+endrewards
+
+// throughput of the system
+rewards "throughput"
+	[in] true : 1;
+endrewards
+

resources/examples/testfiles/ctmc/lrarewards.sm

@@ -0,0 +1,11 @@
+ctmc
+
+module main
+	s : [0..3];
+	[] s < 3 -> 3/2: (s'=s+1);
+	[] s > 0 -> 3: (s'=s-1);
+endmodule
+
+rewards
+	true: s;
+endrewards

src/storm-conv/settings/modules/ConversionGeneralSettings.cpp

@@ -23,7 +23,7 @@ namespace storm {
 
             ConversionGeneralSettings::ConversionGeneralSettings() : ModuleSettings(moduleName) {
                 this->addOption(storm::settings::OptionBuilder(moduleName, helpOptionName, false, "Shows available options, arguments and descriptions.").setShortName(helpOptionShortName)
-                                .addArgument(storm::settings::ArgumentBuilder::createStringArgument("filter", "'frequent' for frequently used options, 'all' for the complete help, or a regular expression to show help for all matching entities.").setDefaultValueString("frequent").build()).build());
+                                .addArgument(storm::settings::ArgumentBuilder::createStringArgument("filter", "'frequent' for frequently used options, 'all' for the complete help, or a regular expression to show help for all matching entities.").setDefaultValueString("frequent").makeOptional().build()).build());
                 this->addOption(storm::settings::OptionBuilder(moduleName, versionOptionName, false, "Prints the version information.").build());
                 this->addOption(storm::settings::OptionBuilder(moduleName, verboseOptionName, false, "Enables more verbose output.").setShortName(verboseOptionShortName).build());
                 this->addOption(storm::settings::OptionBuilder(moduleName, debugOptionName, false, "Enables verbose and debug output.").build());

src/storm-conv/settings/modules/ConversionInputSettings.cpp

@@ -30,7 +30,7 @@ namespace storm {
                 // General
                 this->addOption(storm::settings::OptionBuilder(moduleName, propertyOptionName, false, "Specifies the properties to be checked on the model.").setShortName(propertyOptionShortName)
                                         .addArgument(storm::settings::ArgumentBuilder::createStringArgument("property or filename", "The formula or the file containing the formulas.").build())
-                                        .addArgument(storm::settings::ArgumentBuilder::createStringArgument("filter", "The names of the properties to check.").setDefaultValueString("all").build())
+                                        .addArgument(storm::settings::ArgumentBuilder::createStringArgument("filter", "The names of the properties to check.").setDefaultValueString("all").makeOptional().build())
                                         .build());
                 this->addOption(storm::settings::OptionBuilder(moduleName, constantsOptionName, false, "Specifies the constant replacements to use in symbolic models. Note that this requires the model to be given as an symbolic model (e.g., via --" + prismInputOptionName + ").").setShortName(constantsOptionShortName)
                                 .addArgument(storm::settings::ArgumentBuilder::createStringArgument("values", "A comma separated list of constants and their value, e.g. a=1,b=2,c=3.").setDefaultValueString("").build()).build());
@@ -44,7 +44,7 @@ namespace storm {
                 this->addOption(storm::settings::OptionBuilder(moduleName, janiInputOptionName, false, "Parses the model given in the JANI format.")
                                 .addArgument(storm::settings::ArgumentBuilder::createStringArgument("filename", "The name of the file from which to read the JANI input.").addValidatorString(ArgumentValidatorFactory::createExistingFileValidator()).build()).build());
                 this->addOption(storm::settings::OptionBuilder(moduleName, janiPropertyOptionName, false, "Specifies the properties from the jani model (given by --" + janiInputOptionName + ") to be checked.").setShortName(janiPropertyOptionShortName)
-                                .addArgument(storm::settings::ArgumentBuilder::createStringArgument("values", "A comma separated list of properties to be checked").setDefaultValueString("").build()).build());
+                                .addArgument(storm::settings::ArgumentBuilder::createStringArgument("values", "A comma separated list of properties to be checked").setDefaultValueString("").makeOptional().build()).build());
 
             }
             

src/storm-dft/modelchecker/dft/DFTModelChecker.cpp

@@ -466,13 +466,17 @@ namespace storm {
                 singleModelCheckingTimer.start();
                 //STORM_PRINT_AND_LOG("Model checking property " << *property << " ..." << std::endl);
                 std::unique_ptr<storm::modelchecker::CheckResult> result(
-                        storm::api::verifyWithSparseEngine<ValueType>(model, storm::api::createTask<ValueType>(property,
-                                                                                                               true)));
-                STORM_LOG_ASSERT(result, "Result does not exist.");
-                result->filter(storm::modelchecker::ExplicitQualitativeCheckResult(model->getInitialStates()));
-                ValueType resultValue = result->asExplicitQuantitativeCheckResult<ValueType>().getValueMap().begin()->second;
+                        storm::api::verifyWithSparseEngine<ValueType>(model, storm::api::createTask<ValueType>(property,true)));
+
+                if (result) {
+                    result->filter(storm::modelchecker::ExplicitQualitativeCheckResult(model->getInitialStates()));
+                    ValueType resultValue = result->asExplicitQuantitativeCheckResult<ValueType>().getValueMap().begin()->second;
+                    results.push_back(resultValue);
+                } else {
+                    STORM_LOG_WARN("The property '" << *property << "' could not be checked with the current settings.");
+                    results.push_back(-storm::utility::one<ValueType>());
+                }
                 //STORM_PRINT_AND_LOG("Result (initial states): " << resultValue << std::endl);
-                results.push_back(resultValue);
                 singleModelCheckingTimer.stop();
                 //STORM_PRINT_AND_LOG("Time for model checking: " << singleModelCheckingTimer << "." << std::endl);
             }

src/storm-dft/settings/modules/DftIOSettings.cpp

@@ -50,7 +50,7 @@ namespace storm {
                                         .addArgument(storm::settings::ArgumentBuilder::createDoubleArgument("endtime", "The timebound to end with.")
                                                              .addValidatorDouble(ArgumentValidatorFactory::createDoubleGreaterEqualValidator(0.0)).build())
                                         .addArgument(storm::settings::ArgumentBuilder::createDoubleArgument("inc", "The value to increment with to get the next timepoint.")
-                                                             .addValidatorDouble(ArgumentValidatorFactory::createDoubleGreaterEqualValidator(0.0)).build())
+                                                             .addValidatorDouble(ArgumentValidatorFactory::createDoubleGreaterEqualValidator(0.0)).setDefaultValueDouble(1.0).makeOptional().build())
                                         .build());
                 this->addOption(storm::settings::OptionBuilder(moduleName, minValueOptionName, false, "Compute minimal value in case of non-determinism.").build());
                 this->addOption(storm::settings::OptionBuilder(moduleName, maxValueOptionName, false, "Compute maximal value in case of non-determinism.").build());

src/storm-pars-cli/storm-pars.cpp

@@ -744,7 +744,14 @@ namespace storm {
 
                 STORM_LOG_THROW(regions.size() <= 1, storm::exceptions::InvalidArgumentException, "Monotonicity analysis only allowed on single region");
                 storm::analysis::MonotonicityChecker<ValueType> monotonicityChecker = storm::analysis::MonotonicityChecker<ValueType>(model, formulas, regions, monSettings.isValidateAssumptionsSet(), monSettings.getNumberOfSamples(), monSettings.getMonotonicityAnalysisPrecision());
-                monotonicityChecker.checkMonotonicity();
+                if (ioSettings.isExportMonotonicitySet()) {
+                    std::ofstream outfile;
+                    utility::openFile(ioSettings.getExportMonotonicityFilename(), outfile);
+                    monotonicityChecker.checkMonotonicity(outfile);
+                    utility::closeFile(outfile);
+                } else {
+                    monotonicityChecker.checkMonotonicity(std::cout);
+                }
                 monotonicityWatch.stop();
                 STORM_PRINT(std::endl << "Total time for monotonicity checking: " << monotonicityWatch << "." << std::endl
                                     << std::endl);

src/storm-pars/analysis/MonotonicityChecker.cpp

@@ -67,19 +67,18 @@ namespace storm {
         }
 
         template <typename ValueType>
-        std::map<storm::analysis::Order*, std::map<typename utility::parametric::VariableType<ValueType>::type, std::pair<bool, bool>>> MonotonicityChecker<ValueType>::checkMonotonicity() {
+        std::map<storm::analysis::Order*, std::map<typename utility::parametric::VariableType<ValueType>::type, std::pair<bool, bool>>> MonotonicityChecker<ValueType>::checkMonotonicity(std::ostream& outfile) {
             auto map = createOrder();
             std::shared_ptr<storm::models::sparse::Model<ValueType>> sparseModel = model->as<storm::models::sparse::Model<ValueType>>();
             auto matrix = sparseModel->getTransitionMatrix();
-            return checkMonotonicity(map, matrix);
+            return checkMonotonicity(outfile, map, matrix);
         }
 
         template <typename ValueType>
-        std::map<storm::analysis::Order*, std::map<typename utility::parametric::VariableType<ValueType>::type, std::pair<bool, bool>>> MonotonicityChecker<ValueType>::checkMonotonicity(std::map<storm::analysis::Order*, std::vector<std::shared_ptr<storm::expressions::BinaryRelationExpression>>> map, storm::storage::SparseMatrix<ValueType> matrix) {
+        std::map<storm::analysis::Order*, std::map<typename utility::parametric::VariableType<ValueType>::type, std::pair<bool, bool>>> MonotonicityChecker<ValueType>::checkMonotonicity(std::ostream& outfile, std::map<storm::analysis::Order*, std::vector<std::shared_ptr<storm::expressions::BinaryRelationExpression>>> map, storm::storage::SparseMatrix<ValueType> matrix) {
             storm::utility::Stopwatch monotonicityCheckWatch(true);
             std::map<storm::analysis::Order *, std::map<typename utility::parametric::VariableType<ValueType>::type, std::pair<bool, bool>>> result;
 
-            outfile.open(filename, std::ios_base::app);
 
             if (map.size() == 0) {
                 // Nothing is known
@@ -166,7 +165,6 @@ namespace storm {
             outfile << ", ";
 
             monotonicityCheckWatch.stop();
-            outfile.close();
             return result;
         }
 

src/storm-pars/analysis/MonotonicityChecker.h

@@ -41,7 +41,7 @@ namespace storm {
             /*!
              * Checks for model and formula as provided in constructor for monotonicity
              */
-            std::map<storm::analysis::Order*, std::map<typename utility::parametric::VariableType<ValueType>::type, std::pair<bool, bool>>> checkMonotonicity();
+            std::map<storm::analysis::Order*, std::map<typename utility::parametric::VariableType<ValueType>::type, std::pair<bool, bool>>> checkMonotonicity(std::ostream& outfile);
 
             /*!
              * Checks if monotonicity can be found in this order. Unordered states are not checked
@@ -113,7 +113,7 @@ namespace storm {
             }
 
         private:
-            std::map<storm::analysis::Order*, std::map<typename utility::parametric::VariableType<ValueType>::type, std::pair<bool, bool>>> checkMonotonicity(std::map<storm::analysis::Order*, std::vector<std::shared_ptr<storm::expressions::BinaryRelationExpression>>> map, storm::storage::SparseMatrix<ValueType> matrix);
+            std::map<storm::analysis::Order*, std::map<typename utility::parametric::VariableType<ValueType>::type, std::pair<bool, bool>>> checkMonotonicity(std::ostream& outfile, std::map<storm::analysis::Order*, std::vector<std::shared_ptr<storm::expressions::BinaryRelationExpression>>> map, storm::storage::SparseMatrix<ValueType> matrix);
 
             std::map<typename utility::parametric::VariableType<ValueType>::type, std::pair<bool, bool>> analyseMonotonicity(uint_fast64_t i, Order* order, storm::storage::SparseMatrix<ValueType> matrix) ;
 
@@ -141,10 +141,6 @@ namespace storm {
 
             OrderExtender<ValueType> *extender;
 
-            std::ofstream outfile;
-
-            std::string filename = "monotonicity.txt";
-
             double precision;
 
             storm::storage::ParameterRegion<ValueType> region;

src/storm-pars/settings/ParsSettings.cpp

@@ -18,6 +18,7 @@
 #include "storm/settings/modules/NativeEquationSolverSettings.h"
 #include "storm/settings/modules/TopologicalEquationSolverSettings.h"
 #include "storm/settings/modules/EliminationSettings.h"
+#include "storm/settings/modules/LongRunAverageSolverSettings.h"
 #include "storm/settings/modules/MinMaxEquationSolverSettings.h"
 #include "storm/settings/modules/GameSolverSettings.h"
 #include "storm/settings/modules/BisimulationSettings.h"
@@ -48,6 +49,7 @@ namespace storm {
             storm::settings::addModule<storm::settings::modules::NativeEquationSolverSettings>();
             storm::settings::addModule<storm::settings::modules::TopologicalEquationSolverSettings>();
             storm::settings::addModule<storm::settings::modules::EliminationSettings>();
+            storm::settings::addModule<storm::settings::modules::LongRunAverageSolverSettings>();
             storm::settings::addModule<storm::settings::modules::MinMaxEquationSolverSettings>();
             storm::settings::addModule<storm::settings::modules::GameSolverSettings>();
             storm::settings::addModule<storm::settings::modules::BisimulationSettings>();

src/storm-pars/settings/modules/RegionSettings.cpp

@@ -35,12 +35,12 @@ namespace storm {
                 
                 this->addOption(storm::settings::OptionBuilder(moduleName, refineOptionName, false, "Enables region refinement.")
                                 .addArgument(storm::settings::ArgumentBuilder::createDoubleArgument("coverage-threshold", "Refinement converges if the fraction of unknown area falls below this threshold.").setDefaultValueDouble(0.05).addValidatorDouble(storm::settings::ArgumentValidatorFactory::createDoubleRangeValidatorIncluding(0.0,1.0)).build())
-                                .addArgument(storm::settings::ArgumentBuilder::createIntegerArgument("depth-limit", "If given, limits the number of times a region is refined.").setDefaultValueInteger(-1).build()).build());
+                                .addArgument(storm::settings::ArgumentBuilder::createIntegerArgument("depth-limit", "If given, limits the number of times a region is refined.").setDefaultValueInteger(-1).makeOptional().build()).build());
                 
                 std::vector<std::string> directions = {"min", "max"};
                 this->addOption(storm::settings::OptionBuilder(moduleName, extremumOptionName, false, "Computes the extremum within the region.")
                                 .addArgument(storm::settings::ArgumentBuilder::createStringArgument("direction", "The optimization direction").addValidatorString(storm::settings::ArgumentValidatorFactory::createMultipleChoiceValidator(directions)).build())
-                                .addArgument(storm::settings::ArgumentBuilder::createDoubleArgument("precision", "The desired precision").setDefaultValueDouble(0.05).addValidatorDouble(storm::settings::ArgumentValidatorFactory::createDoubleRangeValidatorIncluding(0.0,1.0)).build()).build());
+                                .addArgument(storm::settings::ArgumentBuilder::createDoubleArgument("precision", "The desired precision").setDefaultValueDouble(0.05).makeOptional().addValidatorDouble(storm::settings::ArgumentValidatorFactory::createDoubleRangeValidatorIncluding(0.0,1.0)).build()).build());
                 
                 std::vector<std::string> engines = {"pl", "exactpl", "validatingpl"};
                 this->addOption(storm::settings::OptionBuilder(moduleName, checkEngineOptionName, true, "Sets which engine is used for analyzing regions.")

src/storm-parsers/parser/ExpressionParser.cpp

@@ -81,7 +81,7 @@ namespace storm {
             literalExpression = qi::lit("true")[qi::_val = phoenix::bind(&ExpressionCreator::createBooleanLiteralExpression, phoenix::ref(*expressionCreator), true, qi::_pass)]
                 | qi::lit("false")[qi::_val = phoenix::bind(&ExpressionCreator::createBooleanLiteralExpression, phoenix::ref(*expressionCreator), false, qi::_pass)]
                 | rationalLiteral_[qi::_val = phoenix::bind(&ExpressionCreator::createRationalLiteralExpression, phoenix::ref(*expressionCreator), qi::_1, qi::_pass)]
-                | qi::int_[qi::_val = phoenix::bind(&ExpressionCreator::createIntegerLiteralExpression, phoenix::ref(*expressionCreator), qi::_1, qi::_pass)];
+                | qi::long_long[qi::_val = phoenix::bind(&ExpressionCreator::createIntegerLiteralExpression, phoenix::ref(*expressionCreator), qi::_1, qi::_pass)];
             literalExpression.name("literal expression");
             
             atomicExpression = floorCeilExpression | roundExpression | prefixPowerModuloExpression | minMaxExpression | (qi::lit("(") >> expression >> qi::lit(")")) | identifierExpression | literalExpression;

src/storm-parsers/parser/ValueParser.h

@@ -19,6 +19,8 @@ namespace storm {
              * Constructor.
              */
             ValueParser() : manager(new storm::expressions::ExpressionManager()), parser(*manager), evaluator(*manager) {
+                // Set empty mapping to enable expression creation even without parameters
+                parser.setIdentifierMapping(identifierMapping);
             }
 
             /*!

src/storm-pgcl/settings/modules/PGCLSettings.cpp

@@ -32,7 +32,7 @@ namespace storm {
                 this->addOption(storm::settings::OptionBuilder(moduleName, programVariableRestrictionsOptionName, false, "Restrictions of program variables").setShortName(programVariableRestrictionShortOptionName).addArgument(storm::settings::ArgumentBuilder::createStringArgument("description", "description of the variable restrictions").build()).build());
                 this->addOption(storm::settings::OptionBuilder(moduleName, propertyOptionName, false, "Specifies the properties to be checked on the model.").setShortName(propertyOptionShortName)
                         .addArgument(storm::settings::ArgumentBuilder::createStringArgument("property or filename", "The formula or the file containing the formulas.").build())
-                        .addArgument(storm::settings::ArgumentBuilder::createStringArgument("filter", "The names of the properties to check.").setDefaultValueString("all").build())
+                        .addArgument(storm::settings::ArgumentBuilder::createStringArgument("filter", "The names of the properties to check.").setDefaultValueString("all").makeOptional().build())
                         .build());
             }
             

src/storm/environment/SubEnvironment.cpp

@@ -8,6 +8,7 @@
 #include "storm/environment/solver/EigenSolverEnvironment.h"
 #include "storm/environment/solver/GmmxxSolverEnvironment.h"
 #include "storm/environment/solver/NativeSolverEnvironment.h"
+#include "storm/environment/solver/LongRunAverageSolverEnvironment.h"
 #include "storm/environment/solver/MinMaxSolverEnvironment.h"
 #include "storm/environment/solver/MultiplierEnvironment.h"
 #include "storm/environment/solver/GameSolverEnvironment.h"
@@ -63,6 +64,7 @@ namespace storm {
     template class SubEnvironment<EigenSolverEnvironment>;
     template class SubEnvironment<GmmxxSolverEnvironment>;
     template class SubEnvironment<NativeSolverEnvironment>;
+    template class SubEnvironment<LongRunAverageSolverEnvironment>;
     template class SubEnvironment<MinMaxSolverEnvironment>;
     template class SubEnvironment<MultiplierEnvironment>;
     template class SubEnvironment<GameSolverEnvironment>;

src/storm/environment/solver/LongRunAverageSolverEnvironment.cpp

@@ -0,0 +1,94 @@
+#include "storm/environment/solver/LongRunAverageSolverEnvironment.h"
+
+#include "storm/settings/SettingsManager.h"
+#include "storm/settings/modules/LongRunAverageSolverSettings.h"
+#include "storm/utility/constants.h"
+#include "storm/utility/macros.h"
+
+namespace storm {
+    
+    LongRunAverageSolverEnvironment::LongRunAverageSolverEnvironment() {
+        auto const& lraSettings = storm::settings::getModule<storm::settings::modules::LongRunAverageSolverSettings>();
+        detMethod = lraSettings.getDetLraMethod();
+        detMethodSetFromDefault = lraSettings.isDetLraMethodSetFromDefaultValue();
+        nondetMethod = lraSettings.getNondetLraMethod();
+        nondetMethodSetFromDefault = lraSettings.isNondetLraMethodSetFromDefaultValue();
+        precision = storm::utility::convertNumber<storm::RationalNumber>(lraSettings.getPrecision());
+        relative = lraSettings.isRelativePrecision();
+        if (lraSettings.isMaximalIterationCountSet()) {
+            maxIters = lraSettings.getMaximalIterationCount();
+        }
+        aperiodicFactor = lraSettings.getAperiodicFactor();
+    }
+    
+    LongRunAverageSolverEnvironment::~LongRunAverageSolverEnvironment() {
+        // Intentionally left empty
+    }
+    
+    storm::solver::LraMethod const& LongRunAverageSolverEnvironment::getDetLraMethod() const {
+        return detMethod;
+    }
+    
+    bool const& LongRunAverageSolverEnvironment::isDetLraMethodSetFromDefault() const {
+        return detMethodSetFromDefault;
+    }
+    
+    void LongRunAverageSolverEnvironment::setDetLraMethod(storm::solver::LraMethod value, bool isSetFromDefault) {
+        detMethod = value;
+        detMethodSetFromDefault = isSetFromDefault;
+    }
+    
+    storm::solver::LraMethod const& LongRunAverageSolverEnvironment::getNondetLraMethod() const {
+        return nondetMethod;
+    }
+    
+    bool const& LongRunAverageSolverEnvironment::isNondetLraMethodSetFromDefault() const {
+        return nondetMethodSetFromDefault;
+    }
+    
+    void LongRunAverageSolverEnvironment::setNondetLraMethod(storm::solver::LraMethod value, bool isSetFromDefault) {
+        nondetMethod = value;
+        nondetMethodSetFromDefault = isSetFromDefault;
+    }
+    
+    storm::RationalNumber const& LongRunAverageSolverEnvironment::getPrecision() const {
+        return precision;
+    }
+    
+    void LongRunAverageSolverEnvironment::setPrecision(storm::RationalNumber value) {
+        precision  = value;
+    }
+    
+    bool const& LongRunAverageSolverEnvironment::getRelativeTerminationCriterion() const {
+        return relative;
+    }
+    
+    void LongRunAverageSolverEnvironment::setRelativeTerminationCriterion(bool value) {
+        relative = value;
+    }
+    
+    bool LongRunAverageSolverEnvironment::isMaximalIterationCountSet() const {
+        return maxIters.is_initialized();
+    }
+    
+    uint64_t LongRunAverageSolverEnvironment::getMaximalIterationCount() const {
+        return maxIters.get();
+    }
+    
+    void LongRunAverageSolverEnvironment::setMaximalIterationCount(uint64_t value) {
+        maxIters = value;
+    }
+    
+    void LongRunAverageSolverEnvironment::unsetMaximalIterationCount() {
+        maxIters = boost::none;
+    }
+    
+    storm::RationalNumber const& LongRunAverageSolverEnvironment::getAperiodicFactor() const {
+        return aperiodicFactor;
+    }
+    
+    void LongRunAverageSolverEnvironment::setAperiodicFactor(storm::RationalNumber value) {
+        aperiodicFactor  = value;
+    }
+
+}

src/storm/environment/solver/LongRunAverageSolverEnvironment.h

@@ -0,0 +1,49 @@
+#pragma once
+
+#include "storm/environment/solver/SolverEnvironment.h"
+#include "storm/solver/SolverSelectionOptions.h"
+
+namespace storm {
+    
+    class LongRunAverageSolverEnvironment {
+    public:
+        
+        LongRunAverageSolverEnvironment();
+        ~LongRunAverageSolverEnvironment();
+        
+        storm::solver::LraMethod const& getDetLraMethod() const;
+        bool const& isDetLraMethodSetFromDefault() const;
+        void setDetLraMethod(storm::solver::LraMethod value, bool isSetFromDefault = false);
+
+        storm::solver::LraMethod const& getNondetLraMethod() const;
+        bool const& isNondetLraMethodSetFromDefault() const;
+        void setNondetLraMethod(storm::solver::LraMethod value, bool isSetFromDefault = false);
+        
+        storm::RationalNumber const& getPrecision() const;
+        void setPrecision(storm::RationalNumber value);
+        bool const& getRelativeTerminationCriterion() const;
+        void setRelativeTerminationCriterion(bool value);
+
+        bool isMaximalIterationCountSet() const;
+        uint64_t getMaximalIterationCount() const;
+        void setMaximalIterationCount(uint64_t value);
+        void unsetMaximalIterationCount();
+        
+        storm::RationalNumber const& getAperiodicFactor() const;
+        void setAperiodicFactor(storm::RationalNumber value);
+
+    private:
+        storm::solver::LraMethod detMethod;
+        bool detMethodSetFromDefault;
+        
+        storm::solver::LraMethod nondetMethod;
+        bool nondetMethodSetFromDefault;
+        
+        storm::RationalNumber precision;
+        bool relative;
+        boost::optional<uint64_t> maxIters;
+        
+        storm::RationalNumber aperiodicFactor;
+    };
+}
+

src/storm/environment/solver/MinMaxSolverEnvironment.cpp

@@ -22,8 +22,6 @@ namespace storm {
         STORM_LOG_ASSERT(considerRelativeTerminationCriterion || minMaxSettings.getConvergenceCriterion() == storm::settings::modules::MinMaxEquationSolverSettings::ConvergenceCriterion::Absolute, "Unknown convergence criterion");
         multiplicationStyle = minMaxSettings.getValueIterationMultiplicationStyle();
         symmetricUpdates = minMaxSettings.isForceIntervalIterationSymmetricUpdatesSet();
-        lraMethod = minMaxSettings.getLraMethod();
-        lraMethodSetFromDefault = minMaxSettings.isLraMethodSetFromDefaultValue();
     }
 
     MinMaxSolverEnvironment::~MinMaxSolverEnvironment() {
@@ -83,19 +81,4 @@ namespace storm {
         symmetricUpdates = value;
     }
     
-    storm::solver::LraMethod const& MinMaxSolverEnvironment::getLraMethod() const {
-        return lraMethod;
-    }
-    
-    bool const& MinMaxSolverEnvironment::isLraMethodSetFromDefault() const {
-        return lraMethodSetFromDefault;
-    }
-    
-    void MinMaxSolverEnvironment::setLraMethod(storm::solver::LraMethod value, bool isSetFromDefault) {
-        lraMethod = value;
-        lraMethodSetFromDefault = isSetFromDefault;
-    }
-    
-
-    
 }

src/storm/environment/solver/MinMaxSolverEnvironment.h

@@ -27,9 +27,6 @@ namespace storm {
         void setMultiplicationStyle(storm::solver::MultiplicationStyle value);
         bool isSymmetricUpdatesSet() const;
         void setSymmetricUpdates(bool value);
-        storm::solver::LraMethod const& getLraMethod() const;
-        bool const& isLraMethodSetFromDefault() const;
-        void setLraMethod(storm::solver::LraMethod value, bool isSetFromDefault = false);
         
     private:
         storm::solver::MinMaxMethod minMaxMethod;
@@ -39,8 +36,6 @@ namespace storm {
         bool considerRelativeTerminationCriterion;
         storm::solver::MultiplicationStyle multiplicationStyle;
         bool symmetricUpdates;
-        storm::solver::LraMethod lraMethod;
-        bool lraMethodSetFromDefault;
     };
 }
 

src/storm/environment/solver/SolverEnvironment.cpp

@@ -1,5 +1,6 @@
 #include "storm/environment/solver/SolverEnvironment.h"
 
+#include "storm/environment/solver/LongRunAverageSolverEnvironment.h"
 #include "storm/environment/solver/MinMaxSolverEnvironment.h"
 #include "storm/environment/solver/MultiplierEnvironment.h"
 #include "storm/environment/solver/EigenSolverEnvironment.h"
@@ -29,6 +30,14 @@ namespace storm {
         // Intentionally left empty
     }
     
+    LongRunAverageSolverEnvironment& SolverEnvironment::lra() {
+        return longRunAverageSolverEnvironment.get();
+    }
+    
+    LongRunAverageSolverEnvironment const& SolverEnvironment::lra() const {
+        return longRunAverageSolverEnvironment.get();
+    }
+    
     MinMaxSolverEnvironment& SolverEnvironment::minMax() {
         return minMaxSolverEnvironment.get();
     }

src/storm/environment/solver/SolverEnvironment.h

@@ -14,6 +14,7 @@ namespace storm {
     class EigenSolverEnvironment;
     class GmmxxSolverEnvironment;
     class NativeSolverEnvironment;
+    class LongRunAverageSolverEnvironment;
     class MinMaxSolverEnvironment;
     class MultiplierEnvironment;
     class GameSolverEnvironment;
@@ -31,6 +32,8 @@ namespace storm {
         GmmxxSolverEnvironment const& gmmxx() const;
         NativeSolverEnvironment& native();
         NativeSolverEnvironment const& native() const;
+        LongRunAverageSolverEnvironment& lra();
+        LongRunAverageSolverEnvironment const& lra() const;
         MinMaxSolverEnvironment& minMax();
         MinMaxSolverEnvironment const& minMax() const;
         MultiplierEnvironment& multiplier();
@@ -56,6 +59,7 @@ namespace storm {
         SubEnvironment<NativeSolverEnvironment> nativeSolverEnvironment;
         SubEnvironment<GameSolverEnvironment> gameSolverEnvironment;
         SubEnvironment<TopologicalSolverEnvironment> topologicalSolverEnvironment;
+        SubEnvironment<LongRunAverageSolverEnvironment> longRunAverageSolverEnvironment;
         SubEnvironment<MinMaxSolverEnvironment> minMaxSolverEnvironment;
         SubEnvironment<MultiplierEnvironment> multiplierEnvironment;
       

src/storm/modelchecker/csl/SparseCtmcCslModelChecker.cpp

@@ -132,16 +132,14 @@ namespace storm {
             std::unique_ptr<CheckResult> subResultPointer = this->check(env, stateFormula);
             ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult();
             
-            storm::storage::SparseMatrix<ValueType> probabilityMatrix = storm::modelchecker::helper::SparseCtmcCslHelper::computeProbabilityMatrix(this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector());
-            std::vector<ValueType> numericResult = storm::modelchecker::helper::SparseCtmcCslHelper::computeLongRunAverageProbabilities(env, storm::solver::SolveGoal<ValueType>(this->getModel(), checkTask), probabilityMatrix, subResult.getTruthValuesVector(), &this->getModel().getExitRateVector());
+            std::vector<ValueType> numericResult = storm::modelchecker::helper::SparseCtmcCslHelper::computeLongRunAverageProbabilities(env, storm::solver::SolveGoal<ValueType>(this->getModel(), checkTask), this->getModel().getTransitionMatrix(), subResult.getTruthValuesVector(), &this->getModel().getExitRateVector());
             return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(std::move(numericResult)));
         }
 
         template <typename SparseCtmcModelType>
         std::unique_ptr<CheckResult> SparseCtmcCslModelChecker<SparseCtmcModelType>::computeLongRunAverageRewards(Environment const& env, storm::logic::RewardMeasureType rewardMeasureType, CheckTask<storm::logic::LongRunAverageRewardFormula, ValueType> const& checkTask) {
-            storm::storage::SparseMatrix<ValueType> probabilityMatrix = storm::modelchecker::helper::SparseCtmcCslHelper::computeProbabilityMatrix(this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector());
             auto rewardModel = storm::utility::createFilteredRewardModel(this->getModel(), checkTask);
-            std::vector<ValueType> numericResult = storm::modelchecker::helper::SparseCtmcCslHelper::computeLongRunAverageRewards(env, storm::solver::SolveGoal<ValueType>(this->getModel(), checkTask), probabilityMatrix, rewardModel.get(), &this->getModel().getExitRateVector());
+            std::vector<ValueType> numericResult = storm::modelchecker::helper::SparseCtmcCslHelper::computeLongRunAverageRewards(env, storm::solver::SolveGoal<ValueType>(this->getModel(), checkTask), this->getModel().getTransitionMatrix(), rewardModel.get(), &this->getModel().getExitRateVector());
             return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(std::move(numericResult)));
         }
         

src/storm/modelchecker/csl/helper/HybridCtmcCslHelper.cpp

@@ -339,19 +339,18 @@ namespace storm {
             
             template<storm::dd::DdType DdType, class ValueType>
             std::unique_ptr<CheckResult> HybridCtmcCslHelper::computeLongRunAverageProbabilities(Environment const& env, storm::models::symbolic::Ctmc<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& rateMatrix, storm::dd::Add<DdType, ValueType> const& exitRateVector, storm::dd::Bdd<DdType> const& psiStates) {
-                storm::dd::Add<DdType, ValueType> probabilityMatrix = computeProbabilityMatrix(rateMatrix, exitRateVector);
                 
                 storm::utility::Stopwatch conversionWatch(true);
                 
                 // Create ODD for the translation.
                 storm::dd::Odd odd = model.getReachableStates().createOdd();
                 
-                storm::storage::SparseMatrix<ValueType> explicitProbabilityMatrix = probabilityMatrix.toMatrix(odd, odd);
+                storm::storage::SparseMatrix<ValueType> explicitRateMatrix = rateMatrix.toMatrix(odd, odd);
                 std::vector<ValueType> explicitExitRateVector = exitRateVector.toVector(odd);
                 conversionWatch.stop();
                 STORM_LOG_INFO("Converting symbolic matrix/vector to explicit representation done in " << conversionWatch.getTimeInMilliseconds() << "ms.");
 
-                std::vector<ValueType> result = storm::modelchecker::helper::SparseCtmcCslHelper::computeLongRunAverageProbabilities(env, storm::solver::SolveGoal<ValueType>(), explicitProbabilityMatrix, psiStates.toVector(odd), &explicitExitRateVector);
+                std::vector<ValueType> result = storm::modelchecker::helper::SparseCtmcCslHelper::computeLongRunAverageProbabilities(env, storm::solver::SolveGoal<ValueType>(), explicitRateMatrix, psiStates.toVector(odd), &explicitExitRateVector);
 
                 return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType, ValueType>(model.getReachableStates(), model.getManager().getBddZero(), model.getManager().template getAddZero<ValueType>(), model.getReachableStates(), std::move(odd), std::move(result)));
             }
@@ -367,12 +366,12 @@ namespace storm {
                 // Create ODD for the translation.
                 storm::dd::Odd odd = model.getReachableStates().createOdd();
                 
-                storm::storage::SparseMatrix<ValueType> explicitProbabilityMatrix = probabilityMatrix.toMatrix(odd, odd);
+                storm::storage::SparseMatrix<ValueType> explicitRateMatrix = rateMatrix.toMatrix(odd, odd);
                 std::vector<ValueType> explicitExitRateVector = exitRateVector.toVector(odd);
                 conversionWatch.stop();
                 STORM_LOG_INFO("Converting symbolic matrix/vector to explicit representation done in " << conversionWatch.getTimeInMilliseconds() << "ms.");
 
-                std::vector<ValueType> result = storm::modelchecker::helper::SparseCtmcCslHelper::computeLongRunAverageRewards(env, storm::solver::SolveGoal<ValueType>(), explicitProbabilityMatrix, rewardModel.getTotalRewardVector(probabilityMatrix, model.getColumnVariables(), exitRateVector, true).toVector(odd), &explicitExitRateVector);
+                std::vector<ValueType> result = storm::modelchecker::helper::SparseCtmcCslHelper::computeLongRunAverageRewards(env, storm::solver::SolveGoal<ValueType>(), explicitRateMatrix, rewardModel.getTotalRewardVector(probabilityMatrix, model.getColumnVariables(), exitRateVector, true).toVector(odd), &explicitExitRateVector);
                 
                 return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType, ValueType>(model.getReachableStates(), model.getManager().getBddZero(), model.getManager().template getAddZero<ValueType>(), model.getReachableStates(), std::move(odd), std::move(result)));
             }

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

@@ -14,6 +14,8 @@
 #include "storm/storage/StronglyConnectedComponentDecomposition.h"
 
 #include "storm/adapters/RationalFunctionAdapter.h"
+#include "storm/environment/solver/LongRunAverageSolverEnvironment.h"
+#include "storm/environment/solver/TopologicalSolverEnvironment.h"
 
 #include "storm/utility/macros.h"
 #include "storm/utility/vector.h"
@@ -25,6 +27,7 @@
 #include "storm/exceptions/InvalidPropertyException.h"
 #include "storm/exceptions/FormatUnsupportedBySolverException.h"
 #include "storm/exceptions/UncheckedRequirementException.h"
+#include "storm/exceptions/NotSupportedException.h"
 
 namespace storm {
     namespace modelchecker {
@@ -370,10 +373,10 @@ namespace storm {
             }
             
             template <typename ValueType>
-            std::vector<ValueType> SparseCtmcCslHelper::computeLongRunAverageProbabilities(Environment const& env, storm::solver::SolveGoal<ValueType>&& goal, storm::storage::SparseMatrix<ValueType> const& probabilityMatrix, storm::storage::BitVector const& psiStates, std::vector<ValueType> const* exitRateVector) {
+            std::vector<ValueType> SparseCtmcCslHelper::computeLongRunAverageProbabilities(Environment const& env, storm::solver::SolveGoal<ValueType>&& goal, storm::storage::SparseMatrix<ValueType> const& rateMatrix, storm::storage::BitVector const& psiStates, std::vector<ValueType> const* exitRateVector) {
                 
                 // If there are no goal states, we avoid the computation and directly return zero.
-                uint_fast64_t numberOfStates = probabilityMatrix.getRowCount();
+                uint_fast64_t numberOfStates = rateMatrix.getRowCount();
                 if (psiStates.empty()) {
                     return std::vector<ValueType>(numberOfStates, storm::utility::zero<ValueType>());
                 }
@@ -386,7 +389,7 @@ namespace storm {
                 ValueType zero = storm::utility::zero<ValueType>();
                 ValueType one = storm::utility::one<ValueType>();
                 
-                return computeLongRunAverages<ValueType>(env, std::move(goal), probabilityMatrix,
+                return computeLongRunAverages<ValueType>(env, std::move(goal), rateMatrix,
                                               [&zero, &one, &psiStates] (storm::storage::sparse::state_type const& state) -> ValueType {
                                                   if (psiStates.get(state)) {
                                                       return one;
@@ -397,16 +400,33 @@ namespace storm {
             }
             
             template <typename ValueType, typename RewardModelType>
-            std::vector<ValueType> SparseCtmcCslHelper::computeLongRunAverageRewards(Environment const& env, storm::solver::SolveGoal<ValueType>&& goal, storm::storage::SparseMatrix<ValueType> const& probabilityMatrix, RewardModelType const& rewardModel, std::vector<ValueType> const* exitRateVector) {
+            std::vector<ValueType> SparseCtmcCslHelper::computeLongRunAverageRewards(Environment const& env, storm::solver::SolveGoal<ValueType>&& goal, storm::storage::SparseMatrix<ValueType> const& rateMatrix, RewardModelType const& rewardModel, std::vector<ValueType> const* exitRateVector) {
                 // Only compute the result if the model has a state-based reward model.
                 STORM_LOG_THROW(!rewardModel.empty(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
 
-                return computeLongRunAverageRewards(env, std::move(goal), probabilityMatrix, rewardModel.getTotalRewardVector(probabilityMatrix, *exitRateVector), exitRateVector);
+                return computeLongRunAverages<ValueType>(env, std::move(goal), rateMatrix,
+                        [&] (storm::storage::sparse::state_type const& state) -> ValueType {
+                            ValueType result = rewardModel.hasStateRewards() ? rewardModel.getStateReward(state) : storm::utility::zero<ValueType>();
+                            if (rewardModel.hasStateActionRewards()) {
+                                // State action rewards are multiplied with the exit rate r(s). Then, multiplying the reward with the expected time we stay at s (i.e. 1/r(s)) yields the original state reward
+                                if (exitRateVector) {
+                                    result += rewardModel.getStateActionReward(state) * (*exitRateVector)[state];
+                                } else {
+                                    result += rewardModel.getStateActionReward(state);
+                                }
+                            }
+                            if (rewardModel.hasTransitionRewards()) {
+                                // Transition rewards are already multiplied with the rates
+                                result += rateMatrix.getPointwiseProductRowSum(rewardModel.getTransitionRewardMatrix(), state);
+                            }
+                            return result;
+                        },
+                        exitRateVector);
             }
             
             template <typename ValueType>
-            std::vector<ValueType> SparseCtmcCslHelper::computeLongRunAverageRewards(Environment const& env, storm::solver::SolveGoal<ValueType>&& goal, storm::storage::SparseMatrix<ValueType> const& probabilityMatrix, std::vector<ValueType> const& stateRewardVector, std::vector<ValueType> const* exitRateVector) {
-                return computeLongRunAverages<ValueType>(env, std::move(goal), probabilityMatrix,
+            std::vector<ValueType> SparseCtmcCslHelper::computeLongRunAverageRewards(Environment const& env, storm::solver::SolveGoal<ValueType>&& goal, storm::storage::SparseMatrix<ValueType> const& rateMatrix, std::vector<ValueType> const& stateRewardVector, std::vector<ValueType> const* exitRateVector) {
+                return computeLongRunAverages<ValueType>(env, std::move(goal), rateMatrix,
                                                          [&stateRewardVector] (storm::storage::sparse::state_type const& state) -> ValueType {
                                                              return stateRewardVector[state];
                                                          },
@@ -414,183 +434,57 @@ namespace storm {
             }
             
             template <typename ValueType>
-            std::vector<ValueType> SparseCtmcCslHelper::computeLongRunAverages(Environment const& env, storm::solver::SolveGoal<ValueType>&& goal, storm::storage::SparseMatrix<ValueType> const& probabilityMatrix, std::function<ValueType (storm::storage::sparse::state_type const& state)> const& valueGetter, std::vector<ValueType> const* exitRateVector){
-                uint_fast64_t numberOfStates = probabilityMatrix.getRowCount();
-
+            std::vector<ValueType> SparseCtmcCslHelper::computeLongRunAverages(Environment const& env, storm::solver::SolveGoal<ValueType>&& goal, storm::storage::SparseMatrix<ValueType> const& rateMatrix, std::function<ValueType (storm::storage::sparse::state_type const& state)> const& valueGetter, std::vector<ValueType> const* exitRateVector){
+                storm::storage::SparseMatrix<ValueType> probabilityMatrix;
+                if (exitRateVector) {
+                    probabilityMatrix = computeProbabilityMatrix(rateMatrix, *exitRateVector);
+                } else {
+                    probabilityMatrix = rateMatrix;
+                }
+                uint_fast64_t numberOfStates = rateMatrix.getRowCount();
+            
                 // Start by decomposing the CTMC into its BSCCs.
-                storm::storage::StronglyConnectedComponentDecomposition<ValueType> bsccDecomposition(probabilityMatrix, storm::storage::StronglyConnectedComponentDecompositionOptions().onlyBottomSccs());
+                storm::storage::StronglyConnectedComponentDecomposition<ValueType> bsccDecomposition(rateMatrix, storm::storage::StronglyConnectedComponentDecompositionOptions().onlyBottomSccs());
                 
                 STORM_LOG_DEBUG("Found " << bsccDecomposition.size() << " BSCCs.");
-                
-                // Get some data members for convenience.
-                ValueType one = storm::utility::one<ValueType>();
-                ValueType zero = storm::utility::zero<ValueType>();
-                
+
                 // Prepare the vector holding the LRA values for each of the BSCCs.
-                std::vector<ValueType> bsccLra(bsccDecomposition.size(), zero);
+                std::vector<ValueType> bsccLra;
+                bsccLra.reserve(bsccDecomposition.size());
                 
-                // First we check which states are in BSCCs.
-                storm::storage::BitVector statesInBsccs(numberOfStates);
-                storm::storage::BitVector firstStatesInBsccs(numberOfStates);
+                auto underlyingSolverEnvironment = env;
+                auto precision = env.solver().lra().getPrecision();
+                if (env.solver().isForceSoundness()) {
+                    // For sound computations, the error in the MECS plus the error in the remaining system should be less then the user defined precision.
+                    precision /= storm::utility::convertNumber<storm::RationalNumber>(2);
+                    underlyingSolverEnvironment.solver().lra().setPrecision(precision);
+                }
+                underlyingSolverEnvironment.solver().setLinearEquationSolverPrecision(precision, env.solver().lra().getRelativeTerminationCriterion());
                 
-                for (uint_fast64_t currentBsccIndex = 0; currentBsccIndex < bsccDecomposition.size(); ++currentBsccIndex) {
-                    storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[currentBsccIndex];
-                    
-                    // Gather information for later use.
-                    bool first = true;
+                // Keep track of the maximal and minimal value occuring in one of the BSCCs
+                ValueType maxValue, minValue;
+                storm::storage::BitVector statesInBsccs(numberOfStates);
+                for (auto const& bscc : bsccDecomposition) {
                     for (auto const& state : bscc) {
                         statesInBsccs.set(state);
-                        if (first) {
-                            firstStatesInBsccs.set(state);
-                        }
-                        first = false;
+                    }
+                    bsccLra.push_back(computeLongRunAveragesForBscc<ValueType>(underlyingSolverEnvironment, bscc, rateMatrix, valueGetter, exitRateVector));
+                    if (bsccLra.size() == 1) {
+                        maxValue = bsccLra.back();
+                        minValue = bsccLra.back();
+                    } else {
+                        maxValue = std::max(bsccLra.back(), maxValue);
+                        minValue = std::min(bsccLra.back(), minValue);
                     }
                 }
-                storm::storage::BitVector statesNotInBsccs = ~statesInBsccs;
                 
+                storm::storage::BitVector statesNotInBsccs = ~statesInBsccs;
                 STORM_LOG_DEBUG("Found " << statesInBsccs.getNumberOfSetBits() << " states in BSCCs.");
                 
-                // Prepare a vector holding the index within all states that are in BSCCs for every state.
-                std::vector<uint_fast64_t> indexInStatesInBsccs;
-                
-                // Prepare a vector that maps the index within the set of all states in BSCCs to the index of the containing BSCC.
-                std::vector<uint_fast64_t> stateToBsccIndexMap;
-                
-                if (!statesInBsccs.empty()) {
-                    firstStatesInBsccs = firstStatesInBsccs % statesInBsccs;
-                    
-                    // Then we construct an equation system that yields the steady state probabilities for all states in BSCCs.
-                    storm::storage::SparseMatrix<ValueType> bsccEquationSystem = probabilityMatrix.getSubmatrix(false, statesInBsccs, statesInBsccs, true);
-                    
-                    // Since in the fix point equation, we need to multiply the vector from the left, we convert this to a
-                    // multiplication from the right by transposing the system.
-                    bsccEquationSystem = bsccEquationSystem.transpose(false, true);
-                    
-                    // Create an auxiliary structure that makes it easy to look up the indices within the set of BSCC states.
-                    uint_fast64_t lastIndex = 0;
-                    uint_fast64_t currentNumberOfSetBits = 0;
-                    indexInStatesInBsccs.reserve(probabilityMatrix.getRowCount());
-                    for (auto index : statesInBsccs) {
-                        while (lastIndex <= index) {
-                            indexInStatesInBsccs.push_back(currentNumberOfSetBits);
-                            ++lastIndex;
-                        }
-                        ++currentNumberOfSetBits;
-                    }
-                    
-                    stateToBsccIndexMap.resize(statesInBsccs.getNumberOfSetBits());
-                    for (uint_fast64_t currentBsccIndex = 0; currentBsccIndex < bsccDecomposition.size(); ++currentBsccIndex) {
-                        storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[currentBsccIndex];
-                        for (auto const& state : bscc) {
-                            stateToBsccIndexMap[indexInStatesInBsccs[state]] = currentBsccIndex;
-                        }
-                    }
-                    
-                    // Build a different system depending on the problem format of the equation solver.
-                    // Check solver requirements.
-                    storm::solver::GeneralLinearEquationSolverFactory<ValueType> linearEquationSolverFactory;
-                    auto requirements = linearEquationSolverFactory.getRequirements(env);
-                    requirements.clearLowerBounds();
-                    requirements.clearUpperBounds();
-                    STORM_LOG_THROW(!requirements.hasEnabledCriticalRequirement(), storm::exceptions::UncheckedRequirementException, "Solver requirements " + requirements.getEnabledRequirementsAsString() + " not checked.");
-                    
-                    bool fixedPointSystem = false;
-                    if (linearEquationSolverFactory.getEquationProblemFormat(env) == storm::solver::LinearEquationSolverProblemFormat::FixedPointSystem) {
-                        fixedPointSystem = true;
-                    }
-
-                    // Now build the final equation system matrix, the initial guess and the right-hand side in one go.
-                    std::vector<ValueType> bsccEquationSystemRightSide(bsccEquationSystem.getColumnCount(), zero);
-                    storm::storage::SparseMatrixBuilder<ValueType> builder;
-                    for (uint_fast64_t row = 0; row < bsccEquationSystem.getRowCount(); ++row) {
-                    
-                        // If the current row is the first one belonging to a BSCC, we substitute it by the constraint that the
-                        // values for states of this BSCC must sum to one. However, in order to have a non-zero value on the
-                        // diagonal, we add the constraint of the BSCC that produces a 1 on the diagonal.
-                        if (firstStatesInBsccs.get(row)) {
-                            uint_fast64_t requiredBscc = stateToBsccIndexMap[row];
-                            storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[requiredBscc];
-
-                            if (fixedPointSystem) {
-                                for (auto const& state : bscc) {
-                                    if (row == indexInStatesInBsccs[state]) {
-                                        builder.addNextValue(row, indexInStatesInBsccs[state], zero);
-                                    } else {
-                                        builder.addNextValue(row, indexInStatesInBsccs[state], -one);
-                                    }
-                                }
-                            } else {
-                                for (auto const& state : bscc) {
-                                    builder.addNextValue(row, indexInStatesInBsccs[state], one);
-                                }
-                            }
-                            
-                            bsccEquationSystemRightSide[row] = one;
-                        } else {
-                            // Otherwise, we copy the row, and subtract 1 from the diagonal (only for the equation solver format).
-                            for (auto& entry : bsccEquationSystem.getRow(row)) {
-                                if (fixedPointSystem || entry.getColumn() != row) {
-                                    builder.addNextValue(row, entry.getColumn(), entry.getValue());
-                                } else {
-                                    builder.addNextValue(row, entry.getColumn(), entry.getValue() - one);
-                                }
-                            }
-                        }
-                    }
-                    
-                    // Create the initial guess for the LRAs. We take a uniform distribution over all states in a BSCC.
-                    std::vector<ValueType> bsccEquationSystemSolution(bsccEquationSystem.getColumnCount(), zero);
-                    for (uint_fast64_t bsccIndex = 0; bsccIndex < bsccDecomposition.size(); ++bsccIndex) {
-                        storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[bsccIndex];
-                        
-                        for (auto const& state : bscc) {
-                            bsccEquationSystemSolution[indexInStatesInBsccs[state]] = one / bscc.size();
-                        }
-                    }
-                    
-                    bsccEquationSystem = builder.build();
-                    {
-                        std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(env, std::move(bsccEquationSystem));
-                        solver->setLowerBound(storm::utility::zero<ValueType>());
-                        solver->setUpperBound(storm::utility::one<ValueType>());
-                        solver->solveEquations(env, bsccEquationSystemSolution, bsccEquationSystemRightSide);
-                    }
-                    
-                    // If exit rates were given, we need to 'fix' the results to also account for the timing behaviour.
-                    if (exitRateVector != nullptr) {
-                        std::vector<ValueType> bsccTotalValue(bsccDecomposition.size(), zero);
-                        for (auto stateIter = statesInBsccs.begin(); stateIter != statesInBsccs.end(); ++stateIter) {
-                            bsccTotalValue[stateToBsccIndexMap[indexInStatesInBsccs[*stateIter]]] += bsccEquationSystemSolution[indexInStatesInBsccs[*stateIter]] * (one / (*exitRateVector)[*stateIter]);
-                        }
-                        
-                        for (auto stateIter = statesInBsccs.begin(); stateIter != statesInBsccs.end(); ++stateIter) {
-                            bsccEquationSystemSolution[indexInStatesInBsccs[*stateIter]] = (bsccEquationSystemSolution[indexInStatesInBsccs[*stateIter]] * (one / (*exitRateVector)[*stateIter])) / bsccTotalValue[stateToBsccIndexMap[indexInStatesInBsccs[*stateIter]]];
-                        }
-                    }
-                                        
-                    // Calculate LRA Value for each BSCC from steady state distribution in BSCCs.
-                    for (uint_fast64_t bsccIndex = 0; bsccIndex < bsccDecomposition.size(); ++bsccIndex) {
-                        storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[bsccIndex];
-                        
-                        for (auto const& state : bscc) {
-                            bsccLra[stateToBsccIndexMap[indexInStatesInBsccs[state]]] += valueGetter(state) * bsccEquationSystemSolution[indexInStatesInBsccs[state]];
-                        }
-                    }
-                    
-                    for (uint_fast64_t bsccIndex = 0; bsccIndex < bsccDecomposition.size(); ++bsccIndex) {
-                        STORM_LOG_DEBUG("Found LRA " << bsccLra[bsccIndex] << " for BSCC " << bsccIndex << ".");
-                    }
-                } else {
-                    for (uint_fast64_t bsccIndex = 0; bsccIndex < bsccDecomposition.size(); ++bsccIndex) {
-                        storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[bsccIndex];
-                        
-                        // At this point, all BSCCs are known to contain exactly one state, which is why we can set all values
-                        // directly (based on whether or not the contained state is a psi state).
-                        bsccLra[bsccIndex] = valueGetter(*bscc.begin());
-                    }
-                    
-                    for (uint_fast64_t bsccIndex = 0; bsccIndex < bsccDecomposition.size(); ++bsccIndex) {
-                        STORM_LOG_DEBUG("Found LRA " << bsccLra[bsccIndex] << " for BSCC " << bsccIndex << ".");
+                std::vector<uint64_t> stateToBsccMap(statesInBsccs.size(), -1);
+                for (uint64_t bsccIndex = 0; bsccIndex < bsccDecomposition.size(); ++bsccIndex) {
+                    for (auto const& state : bsccDecomposition[bsccIndex]) {
+                        stateToBsccMap[state] = bsccIndex;
                     }
                 }
                 
@@ -605,30 +499,42 @@ namespace storm {
                     rewardRightSide.reserve(statesNotInBsccs.getNumberOfSetBits());
                     
                     for (auto state : statesNotInBsccs) {
-                        ValueType reward = zero;
-                        for (auto entry : probabilityMatrix.getRow(state)) {
+                        ValueType reward = storm::utility::zero<ValueType>();
+                        for (auto entry : rateMatrix.getRow(state)) {
                             if (statesInBsccs.get(entry.getColumn())) {
-                                reward += entry.getValue() * bsccLra[stateToBsccIndexMap[indexInStatesInBsccs[entry.getColumn()]]];
+                                if (exitRateVector) {
+                                    reward += (entry.getValue() / (*exitRateVector)[state]) * bsccLra[stateToBsccMap[entry.getColumn()]];
+                                } else {
+                                    reward += entry.getValue() * bsccLra[stateToBsccMap[entry.getColumn()]];
+                                }
                             }
                         }
                         rewardRightSide.push_back(reward);
                     }
                     
-                    storm::storage::SparseMatrix<ValueType> rewardEquationSystemMatrix = probabilityMatrix.getSubmatrix(false, statesNotInBsccs, statesNotInBsccs, true);
-                    rewardEquationSystemMatrix.convertToEquationSystem();
-                    
-                    rewardSolution = std::vector<ValueType>(rewardEquationSystemMatrix.getColumnCount(), one);
-                    
-                    {
-                        storm::solver::GeneralLinearEquationSolverFactory<ValueType> linearEquationSolverFactory;
-                        // Check solver requirements
-                        auto requirements = linearEquationSolverFactory.getRequirements(env);
-                        STORM_LOG_THROW(!requirements.hasEnabledCriticalRequirement(), storm::exceptions::UncheckedRequirementException, "Solver requirements " + requirements.getEnabledRequirementsAsString() + " not checked.");
-                        // Check whether we have the right input format for the solver.
-                        STORM_LOG_THROW(linearEquationSolverFactory.getEquationProblemFormat(env) == storm::solver::LinearEquationSolverProblemFormat::EquationSystem, storm::exceptions::FormatUnsupportedBySolverException, "The selected solver does not support the required format.");
-                        std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(env, std::move(rewardEquationSystemMatrix));
-                        solver->solveEquations(env, rewardSolution, rewardRightSide);
+                    // Compute reachability rewards
+                    storm::solver::GeneralLinearEquationSolverFactory<ValueType> linearEquationSolverFactory;
+                    bool isEqSysFormat = linearEquationSolverFactory.getEquationProblemFormat(underlyingSolverEnvironment) == storm::solver::LinearEquationSolverProblemFormat::EquationSystem;
+                    storm::storage::SparseMatrix<ValueType> rewardEquationSystemMatrix = rateMatrix.getSubmatrix(false, statesNotInBsccs, statesNotInBsccs, isEqSysFormat);
+                    if (exitRateVector) {
+                        uint64_t localRow = 0;
+                        for (auto const& globalRow : statesNotInBsccs) {
+                            for (auto& entry : rewardEquationSystemMatrix.getRow(localRow)) {
+                                entry.setValue(entry.getValue() / (*exitRateVector)[globalRow]);
+                            }
+                            ++localRow;
+                        }
+                    }
+                    if (isEqSysFormat) {
+                        rewardEquationSystemMatrix.convertToEquationSystem();
                     }
+                    rewardSolution = std::vector<ValueType>(rewardEquationSystemMatrix.getColumnCount(), (maxValue + minValue) / storm::utility::convertNumber<ValueType,uint64_t>(2));
+                    std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(underlyingSolverEnvironment, std::move(rewardEquationSystemMatrix));
+                    solver->setBounds(minValue, maxValue);
+                    // Check solver requirements
+                    auto requirements = solver->getRequirements(underlyingSolverEnvironment);
+                    STORM_LOG_THROW(!requirements.hasEnabledCriticalRequirement(), storm::exceptions::UncheckedRequirementException, "Solver requirements " + requirements.getEnabledRequirementsAsString() + " not checked.");
+                    solver->solveEquations(underlyingSolverEnvironment, rewardSolution, rewardRightSide);
                 }
                 
                 // Fill the result vector.
@@ -653,6 +559,365 @@ namespace storm {
                 return result;
             }
 
+            template <typename ValueType>
+            ValueType SparseCtmcCslHelper::computeLongRunAveragesForBscc(Environment const& env, storm::storage::StronglyConnectedComponent const& bscc, storm::storage::SparseMatrix<ValueType> const& rateMatrix, std::function<ValueType (storm::storage::sparse::state_type const& state)> const& valueGetter, std::vector<ValueType> const* exitRateVector) {
+
+                // Catch the case where all values are the same (this includes the special case where the bscc is of size 1)
+                auto it = bscc.begin();
+                ValueType val = valueGetter(*it);
+                for (++it; it != bscc.end(); ++it) {
+                    if (valueGetter(*it) != val) {
+                        break;
+                    }
+                }
+                if (it == bscc.end()) {
+                    // All entries have the same LRA
+                    return val;
+                }
+                
+                storm::solver::LraMethod method = env.solver().lra().getDetLraMethod();
+                if (storm::NumberTraits<ValueType>::IsExact && env.solver().lra().isDetLraMethodSetFromDefault() && method == storm::solver::LraMethod::ValueIteration) {
+                    method = storm::solver::LraMethod::GainBiasEquations;
+                    STORM_LOG_INFO("Selecting " << storm::solver::toString(method) << " as the solution technique for long-run properties to guarantee exact results. If you want to override this, please explicitly specify a different LRA method.");
+                } else if (env.solver().isForceSoundness() && env.solver().lra().isDetLraMethodSetFromDefault() && method != storm::solver::LraMethod::ValueIteration) {
+                    method = storm::solver::LraMethod::ValueIteration;
+                    STORM_LOG_INFO("Selecting " << storm::solver::toString(method) << " as the solution technique for long-run properties to guarantee sound results. If you want to override this, please explicitly specify a different LRA method.");
+                }
+                STORM_LOG_TRACE("Computing LRA for BSCC of size " << bscc.size() << " using '" << storm::solver::toString(method) << "'.");
+                if (method == storm::solver::LraMethod::ValueIteration) {
+                    return computeLongRunAveragesForBsccVi<ValueType>(env, bscc, rateMatrix, valueGetter, exitRateVector);
+                } else if (method == storm::solver::LraMethod::LraDistributionEquations) {
+                    // We only need the first element of the pair as the lra distribution is not relevant at this point.
+                    return computeLongRunAveragesForBsccLraDistr<ValueType>(env, bscc, rateMatrix, valueGetter, exitRateVector).first;
+                }
+                STORM_LOG_WARN_COND(method == storm::solver::LraMethod::GainBiasEquations, "Unsupported lra method selected. Defaulting to " << storm::solver::toString(storm::solver::LraMethod::GainBiasEquations) << ".");
+                // We don't need the bias values
+                return computeLongRunAveragesForBsccGainBias<ValueType>(env, bscc, rateMatrix, valueGetter, exitRateVector).first;
+            }
+            
+            template <>
+            storm::RationalFunction SparseCtmcCslHelper::computeLongRunAveragesForBsccVi<storm::RationalFunction>(Environment const&, storm::storage::StronglyConnectedComponent const&, storm::storage::SparseMatrix<storm::RationalFunction> const&, std::function<storm::RationalFunction (storm::storage::sparse::state_type const& state)> const&, std::vector<storm::RationalFunction> const*) {
+                STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "The requested Method for LRA computation is not supported for parametric models.");
+            }
+                
+            template <typename ValueType>
+            ValueType SparseCtmcCslHelper::computeLongRunAveragesForBsccVi(Environment const& env, storm::storage::StronglyConnectedComponent const& bscc, storm::storage::SparseMatrix<ValueType> const& rateMatrix, std::function<ValueType (storm::storage::sparse::state_type const& state)> const& valueGetter, std::vector<ValueType> const* exitRateVector) {
+                
+                // Initialize data about the bscc
+                storm::storage::BitVector bsccStates(rateMatrix.getRowGroupCount(), false);
+                for (auto const& state : bscc) {
+                    bsccStates.set(state);
+                }
+                
+                // Get the uniformization rate
+                ValueType uniformizationRate = storm::utility::one<ValueType>();
+                if (exitRateVector) {
+                    uniformizationRate = storm::utility::vector::max_if(*exitRateVector, bsccStates);
+                }
+                // To ensure that the model is aperiodic, we need to make sure that every Markovian state gets a self loop.
+                // Hence, we increase the uniformization rate a little.
+                uniformizationRate *= (storm::utility::one<ValueType>() + storm::utility::convertNumber<ValueType>(env.solver().lra().getAperiodicFactor()));
+
+                // Get the transitions of the submodel
+                typename storm::storage::SparseMatrix<ValueType> bsccMatrix = rateMatrix.getSubmatrix(true, bsccStates, bsccStates, true);
+                
+                // Uniformize the transitions
+                uint64_t subState = 0;
+                for (auto state : bsccStates) {
+                    for (auto& entry : bsccMatrix.getRow(subState)) {
+                        if (entry.getColumn() == subState) {
+                            if (exitRateVector) {
+                                entry.setValue(storm::utility::one<ValueType>() + (entry.getValue() - (*exitRateVector)[state]) / uniformizationRate);
+                            } else {
+                                entry.setValue(storm::utility::one<ValueType>() + (entry.getValue() - storm::utility::one<ValueType>()) / uniformizationRate);
+                            }
+                        } else {
+                            entry.setValue(entry.getValue() / uniformizationRate);
+                        }
+                    }
+                    ++subState;
+                }
+
+                // Compute the rewards obtained in a single uniformization step
+                std::vector<ValueType> markovianRewards;
+                markovianRewards.reserve(bsccMatrix.getRowCount());
+                for (auto const& state : bsccStates) {
+                    markovianRewards.push_back(valueGetter(state) / uniformizationRate);
+                }
+                
+                // start the iterations
+                ValueType precision = storm::utility::convertNumber<ValueType>(env.solver().lra().getPrecision()) / uniformizationRate;
+                bool relative = env.solver().lra().getRelativeTerminationCriterion();
+                if (!relative) {
+                    precision /= uniformizationRate;
+                }
+                std::vector<ValueType> x(bsccMatrix.getRowCount(), storm::utility::zero<ValueType>());
+                std::vector<ValueType> xPrime(x.size());
+                
+                auto multiplier = storm::solver::MultiplierFactory<ValueType>().create(env, bsccMatrix);
+                ValueType maxDiff, minDiff;
+                uint64_t iter = 0;
+                boost::optional<uint64_t> maxIter;
+                if (env.solver().lra().isMaximalIterationCountSet()) {
+                    maxIter = env.solver().lra().getMaximalIterationCount();
+                }
+                while (!maxIter.is_initialized() || iter < maxIter.get()) {
+                    ++iter;
+                    // Compute the values for the markovian states. We also keep track of the maximal and minimal difference between two values (for convergence checking)
+                    multiplier->multiply(env, x, &markovianRewards, xPrime);
+                    
+                    // update xPrime and check for convergence
+                    // to avoid large (and numerically unstable) x-values, we substract a reference value.
+                    auto xIt = x.begin();
+                    auto xPrimeIt = xPrime.begin();
+                    ValueType refVal = *xPrimeIt;
+                    maxDiff = *xPrimeIt - *xIt;
+                    minDiff = maxDiff;
+                    *xPrimeIt -= refVal;
+                    *xIt = *xPrimeIt;
+                    for (++xIt, ++xPrimeIt; xIt != x.end(); ++xIt, ++xPrimeIt) {
+                        ValueType diff = *xPrimeIt - *xIt;
+                        maxDiff = std::max(maxDiff, diff);
+                        minDiff = std::min(minDiff, diff);
+                        *xPrimeIt -= refVal;
+                        *xIt = *xPrimeIt;
+                    }
+                    
+                    // Check for convergence. The uniformization rate is already incorporated into the precision parameter
+                    if ((maxDiff - minDiff) <= (relative ? (precision * minDiff) : precision)) {
+                        break;
+                    }
+                }
+                if (maxIter.is_initialized() && iter == maxIter.get()) {
+                    STORM_LOG_WARN("LRA computation did not converge within " << iter << " iterations.");
+                } else {
+                    STORM_LOG_TRACE("LRA computation converged after " << iter << " iterations.");
+                }
+                return (maxDiff + minDiff) * uniformizationRate / (storm::utility::convertNumber<ValueType>(2.0));
+            }
+            
+            template <typename ValueType>
+            std::pair<ValueType, std::vector<ValueType>> SparseCtmcCslHelper::computeLongRunAveragesForBsccGainBias(Environment const& env, storm::storage::StronglyConnectedComponent const& bscc, storm::storage::SparseMatrix<ValueType> const& rateMatrix, std::function<ValueType (storm::storage::sparse::state_type const& state)> const& valueGetter, std::vector<ValueType> const* exitRateVector) {
+                // We build the equation system as in Line 3 of Algorithm 3 from
+                // Kretinsky, Meggendorfer: Efficient Strategy Iteration for Mean Payoff in Markov Decision Processes (ATVA 2017)
+                // The first variable corresponds to the gain of the bscc whereas the subsequent variables yield the bias for each state s_1, s_2, ....
+                // No bias variable for s_0 is needed since it is always set to zero, yielding an nxn equation system matrix
+                // To make this work for CTMC, we could uniformize the model. This preserves LRA and ensures that we can compute the
+                // LRA as for a DTMC (the soujourn time in each state is the same). If we then multiply the equations with the uniformization rate,
+                // the uniformization rate cancels out. Hence, we obtain the equation system below.
+                
+                // Get a mapping from global state indices to local ones.
+                std::unordered_map<uint64_t, uint64_t> toLocalIndexMap;
+                uint64_t localIndex = 0;
+                for (auto const& globalIndex : bscc) {
+                    toLocalIndexMap[globalIndex] = localIndex;
+                    ++localIndex;
+                }
+                
+                // Prepare an environment for the underlying equation solver
+                auto subEnv = env;
+                if (subEnv.solver().getLinearEquationSolverType() == storm::solver::EquationSolverType::Topological) {
+                    // Topological solver does not make any sense since the BSCC is connected.
+                    subEnv.solver().setLinearEquationSolverType(subEnv.solver().topological().getUnderlyingEquationSolverType());
+                }
+                subEnv.solver().setLinearEquationSolverPrecision(env.solver().lra().getPrecision(), env.solver().lra().getRelativeTerminationCriterion());
+                
+                // Build the equation system matrix and vector.
+                storm::solver::GeneralLinearEquationSolverFactory<ValueType> linearEquationSolverFactory;
+                bool isEquationSystemFormat = linearEquationSolverFactory.getEquationProblemFormat(subEnv) == storm::solver::LinearEquationSolverProblemFormat::EquationSystem;
+                storm::storage::SparseMatrixBuilder<ValueType> builder(bscc.size(), bscc.size());
+                std::vector<ValueType> eqSysVector;
+                eqSysVector.reserve(bscc.size());
+                // The first row asserts that the weighted bias variables and the reward at s_0 sum up to the gain
+                uint64_t row = 0;
+                ValueType entryValue;
+                for (auto const& globalState : bscc) {
+                    // Coefficient for the gain variable
+                    if (isEquationSystemFormat) {
+                        // '1-0' in row 0 and -(-1) in other rows
+                        builder.addNextValue(row, 0, storm::utility::one<ValueType>());
+                    } else if (row > 0) {
+                        // No coeficient in row 0, othwerise substract the gain
+                        builder.addNextValue(row, 0, -storm::utility::one<ValueType>());
+                    }
+                    // Compute weighted sum over successor state. As this is a BSCC, each successor state will again be in the BSCC.
+                    auto diagonalValue = storm::utility::zero<ValueType>();
+                    if (row > 0) {
+                        if (isEquationSystemFormat) {
+                            diagonalValue = exitRateVector ? (*exitRateVector)[globalState] : storm::utility::one<ValueType>();
+                        } else {
+                            diagonalValue = storm::utility::one<ValueType>() - (exitRateVector ? (*exitRateVector)[globalState] : storm::utility::one<ValueType>());
+                        }
+                    }
+                    bool needDiagonalEntry = !storm::utility::isZero(diagonalValue);
+                    for (auto const& entry : rateMatrix.getRow(globalState)) {
+                        uint64_t col = toLocalIndexMap[entry.getColumn()];
+                        if (col == 0) {
+                            //Skip transition to state_0. This corresponds to setting the bias of state_0 to zero
+                            continue;
+                        }
+                        entryValue = entry.getValue();
+                        if (isEquationSystemFormat) {
+                            entryValue = -entryValue;
+                        }
+                        if (needDiagonalEntry && col >= row) {
+                            if (col == row) {
+                                entryValue += diagonalValue;
+                            } else { // col > row
+                                builder.addNextValue(row, row, diagonalValue);
+                            }
+                            needDiagonalEntry = false;
+                        }
+                        builder.addNextValue(row, col, entryValue);
+                    }
+                    if (needDiagonalEntry) {
+                        builder.addNextValue(row, row, diagonalValue);
+                    }
+                    eqSysVector.push_back(valueGetter(globalState));
+                    ++row;
+                }
+
+                // Create a linear equation solver
+                auto solver = linearEquationSolverFactory.create(subEnv, builder.build());
+                // Check solver requirements.
+                auto requirements = solver->getRequirements(subEnv);
+                STORM_LOG_THROW(!requirements.hasEnabledCriticalRequirement(), storm::exceptions::UncheckedRequirementException, "Solver requirements " + requirements.getEnabledRequirementsAsString() + " not checked.");
+                // Todo: Find bounds on the bias variables. Just inserting the maximal value from the vector probably does not work.
+                
+                std::vector<ValueType> eqSysSol(bscc.size(), storm::utility::zero<ValueType>());
+                // Take the mean of the rewards as an initial guess for the gain
+                //eqSysSol.front() = std::accumulate(eqSysVector.begin(), eqSysVector.end(), storm::utility::zero<ValueType>()) / storm::utility::convertNumber<ValueType, uint64_t>(bscc.size());
+                solver->solveEquations(subEnv, eqSysSol, eqSysVector);
+                
+                ValueType gain = eqSysSol.front();
+                // insert bias value for state 0
+                eqSysSol.front() = storm::utility::zero<ValueType>();
+                // Return the gain and the bias values
+                return std::pair<ValueType, std::vector<ValueType>>(std::move(gain), std::move(eqSysSol));
+            }
+            
+            template <typename ValueType>
+            std::pair<ValueType, std::vector<ValueType>> SparseCtmcCslHelper::computeLongRunAveragesForBsccLraDistr(Environment const& env, storm::storage::StronglyConnectedComponent const& bscc, storm::storage::SparseMatrix<ValueType> const& rateMatrix, std::function<ValueType (storm::storage::sparse::state_type const& state)> const& valueGetter, std::vector<ValueType> const* exitRateVector) {
+
+                // Let A be ab auxiliary Matrix with A[s,s] =  R(s,s) - r(s) & A[s,s'] = R(s,s') for s,s' in BSCC and s!=s'.
+                // We build and solve the equation system for
+                // x*A=0 &  x_0+...+x_n=1  <=>  A^t*x=0=x-x & x_0+...+x_n=1  <=> (1+A^t)*x = x & 1-x_0-...-x_n-1=x_n
+                // Then, x[i] will be the fraction of the time we are in state i.
+                
+                // This method assumes that this BSCC consist of more than one state
+                if (bscc.size() == 1) {
+                    return { valueGetter(*bscc.begin()), {storm::utility::one<ValueType>()} };
+                }
+                
+                // Prepare an environment for the underlying linear equation solver
+                auto subEnv = env;
+                if (subEnv.solver().getLinearEquationSolverType() == storm::solver::EquationSolverType::Topological) {
+                    // Topological solver does not make any sense since the BSCC is connected.
+                    subEnv.solver().setLinearEquationSolverType(subEnv.solver().topological().getUnderlyingEquationSolverType());
+                }
+                subEnv.solver().setLinearEquationSolverPrecision(env.solver().lra().getPrecision(), env.solver().lra().getRelativeTerminationCriterion());
+                
+                // Get a mapping from global state indices to local ones as well as a bitvector containing states within the BSCC.
+                std::unordered_map<uint64_t, uint64_t> toLocalIndexMap;
+                storm::storage::BitVector bsccStates(rateMatrix.getRowCount(), false);
+                uint64_t localIndex = 0;
+                for (auto const& globalIndex : bscc) {
+                    bsccStates.set(globalIndex, true);
+                    toLocalIndexMap[globalIndex] = localIndex;
+                    ++localIndex;
+                }
+                
+                // Build the auxiliary Matrix A.
+                auto auxMatrix = rateMatrix.getSubmatrix(false, bsccStates, bsccStates, true); // add diagonal entries!
+                uint64_t row = 0;
+                for (auto const& globalIndex : bscc) {
+                    for (auto& entry : auxMatrix.getRow(row)) {
+                        if (entry.getColumn() == row) {
+                            // This value is non-zero since we have a BSCC with more than one state
+                            if (exitRateVector) {
+                                entry.setValue(entry.getValue() - (*exitRateVector)[globalIndex]);
+                            } else {
+                                entry.setValue(entry.getValue() - storm::utility::one<ValueType>());
+                            }
+                        }
+                    }
+                    ++row;
+                }
+                assert(row == auxMatrix.getRowCount());
+                
+                // We need to consider A^t. This will not delete diagonal entries since they are non-zero.
+                auxMatrix = auxMatrix.transpose();
+                
+                // Check whether we need the fixpoint characterization
+                storm::solver::GeneralLinearEquationSolverFactory<ValueType> linearEquationSolverFactory;
+                bool isFixpointFormat = linearEquationSolverFactory.getEquationProblemFormat(subEnv) == storm::solver::LinearEquationSolverProblemFormat::FixedPointSystem;
+                if (isFixpointFormat) {
+                    // Add a 1 on the diagonal
+                    for (row = 0; row < auxMatrix.getRowCount(); ++row) {
+                        for (auto& entry : auxMatrix.getRow(row)) {
+                            if (entry.getColumn() == row) {
+                                entry.setValue(storm::utility::one<ValueType>() + entry.getValue());
+                            }
+                        }
+                    }
+                }
+                
+                // We now build the equation system matrix.
+                // We can drop the last row of A and add ones in this row instead to assert that the variables sum up to one
+                // Phase 1: replace the existing entries of the last row with ones
+                uint64_t col = 0;
+                uint64_t lastRow = auxMatrix.getRowCount() - 1;
+                for (auto& entry : auxMatrix.getRow(lastRow)) {
+                    entry.setColumn(col);
+                    if (isFixpointFormat) {
+                        if (col == lastRow) {
+                            entry.setValue(storm::utility::zero<ValueType>());
+                        } else {
+                            entry.setValue(-storm::utility::one<ValueType>());
+                        }
+                    } else {
+                        entry.setValue(storm::utility::one<ValueType>());
+                    }
+                    ++col;
+                }
+                storm::storage::SparseMatrixBuilder<ValueType> builder(std::move(auxMatrix));
+                for (; col <= lastRow; ++col) {
+                    if (isFixpointFormat) {
+                        if (col != lastRow) {
+                            builder.addNextValue(lastRow, col, -storm::utility::one<ValueType>());
+                        }
+                    } else {
+                        builder.addNextValue(lastRow, col, storm::utility::one<ValueType>());
+                    }
+                }
+                
+                std::vector<ValueType> bsccEquationSystemRightSide(bscc.size(), storm::utility::zero<ValueType>());
+                bsccEquationSystemRightSide.back() = storm::utility::one<ValueType>();
+                
+                // Create a linear equation solver
+                auto solver = linearEquationSolverFactory.create(subEnv,  builder.build());
+                solver->setBounds(storm::utility::zero<ValueType>(), storm::utility::one<ValueType>());
+                // Check solver requirements.
+                auto requirements = solver->getRequirements(subEnv);
+                requirements.clearLowerBounds();
+                requirements.clearUpperBounds();
+                STORM_LOG_THROW(!requirements.hasEnabledCriticalRequirement(), storm::exceptions::UncheckedRequirementException, "Solver requirements " + requirements.getEnabledRequirementsAsString() + " not checked.");
+                
+                std::vector<ValueType> lraDistr(bscc.size(), storm::utility::one<ValueType>() / storm::utility::convertNumber<ValueType, uint64_t>(bscc.size()));
+                solver->solveEquations(subEnv, lraDistr, bsccEquationSystemRightSide);
+                
+                // Calculate final LRA Value
+                ValueType result = storm::utility::zero<ValueType>();
+                auto solIt = lraDistr.begin();
+                for (auto const& globalState : bscc) {
+                    result += valueGetter(globalState) * (*solIt);
+                    ++solIt;
+                }
+                assert(solIt == lraDistr.end());
+
+                return std::pair<ValueType, std::vector<ValueType>>(std::move(result), std::move(lraDistr));
+            }
+            
             template <typename ValueType, typename std::enable_if<storm::NumberTraits<ValueType>::SupportsExponential, int>::type>
             std::vector<ValueType> SparseCtmcCslHelper::computeAllTransientProbabilities(Environment const& env, storm::storage::SparseMatrix<ValueType> const& rateMatrix, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, std::vector<ValueType> const& exitRates, double timeBound) {
 
@@ -873,9 +1138,9 @@ namespace storm {
             
             template std::vector<double> SparseCtmcCslHelper::computeTotalRewards(Environment const& env, storm::solver::SolveGoal<double>&& goal, storm::storage::SparseMatrix<double> const& rateMatrix, storm::storage::SparseMatrix<double> const& backwardTransitions, std::vector<double> const& exitRateVector, storm::models::sparse::StandardRewardModel<double> const& rewardModel, bool qualitative);
             
-            template std::vector<double> SparseCtmcCslHelper::computeLongRunAverageProbabilities(Environment const& env, storm::solver::SolveGoal<double>&& goal, storm::storage::SparseMatrix<double> const& probabilityMatrix, storm::storage::BitVector const& psiStates, std::vector<double> const* exitRateVector);
-            template std::vector<double> SparseCtmcCslHelper::computeLongRunAverageRewards(Environment const& env, storm::solver::SolveGoal<double>&& goal, storm::storage::SparseMatrix<double> const& probabilityMatrix, storm::models::sparse::StandardRewardModel<double> const& rewardModel, std::vector<double> const* exitRateVector);
-            template std::vector<double> SparseCtmcCslHelper::computeLongRunAverageRewards(Environment const& env, storm::solver::SolveGoal<double>&& goal, storm::storage::SparseMatrix<double> const& probabilityMatrix, std::vector<double> const& stateRewardVector, std::vector<double> const* exitRateVector);
+            template std::vector<double> SparseCtmcCslHelper::computeLongRunAverageProbabilities(Environment const& env, storm::solver::SolveGoal<double>&& goal, storm::storage::SparseMatrix<double> const& rateMatrix, storm::storage::BitVector const& psiStates, std::vector<double> const* exitRateVector);
+            template std::vector<double> SparseCtmcCslHelper::computeLongRunAverageRewards(Environment const& env, storm::solver::SolveGoal<double>&& goal, storm::storage::SparseMatrix<double> const& rateMatrix, storm::models::sparse::StandardRewardModel<double> const& rewardModel, std::vector<double> const* exitRateVector);
+            template std::vector<double> SparseCtmcCslHelper::computeLongRunAverageRewards(Environment const& env, storm::solver::SolveGoal<double>&& goal, storm::storage::SparseMatrix<double> const& rateMatrix, std::vector<double> const& stateRewardVector, std::vector<double> const* exitRateVector);
             
             template std::vector<double> SparseCtmcCslHelper::computeCumulativeRewards(Environment const& env, storm::solver::SolveGoal<double>&& goal, storm::storage::SparseMatrix<double> const& rateMatrix, std::vector<double> const& exitRateVector, storm::models::sparse::StandardRewardModel<double> const& rewardModel, double timeBound);
 
@@ -910,14 +1175,14 @@ namespace storm {
             template std::vector<storm::RationalNumber> SparseCtmcCslHelper::computeTotalRewards(Environment const& env, storm::solver::SolveGoal<storm::RationalNumber>&& goal, storm::storage::SparseMatrix<storm::RationalNumber> const& rateMatrix, storm::storage::SparseMatrix<storm::RationalNumber> const& backwardTransitions, std::vector<storm::RationalNumber> const& exitRateVector, storm::models::sparse::StandardRewardModel<storm::RationalNumber> const& rewardModel, bool qualitative);
             template std::vector<storm::RationalFunction> SparseCtmcCslHelper::computeTotalRewards(Environment const& env, storm::solver::SolveGoal<storm::RationalFunction>&& goal, storm::storage::SparseMatrix<storm::RationalFunction> const& rateMatrix, storm::storage::SparseMatrix<storm::RationalFunction> const& backwardTransitions, std::vector<storm::RationalFunction> const& exitRateVector, storm::models::sparse::StandardRewardModel<storm::RationalFunction> const& rewardModel, bool qualitative);
 
-            template std::vector<storm::RationalNumber> SparseCtmcCslHelper::computeLongRunAverageProbabilities(Environment const& env, storm::solver::SolveGoal<storm::RationalNumber>&& goal, storm::storage::SparseMatrix<storm::RationalNumber> const& probabilityMatrix, storm::storage::BitVector const& psiStates, std::vector<storm::RationalNumber> const* exitRateVector);
-            template std::vector<storm::RationalFunction> SparseCtmcCslHelper::computeLongRunAverageProbabilities(Environment const& env, storm::solver::SolveGoal<storm::RationalFunction>&& goal, storm::storage::SparseMatrix<storm::RationalFunction> const& probabilityMatrix, storm::storage::BitVector const& psiStates, std::vector<storm::RationalFunction> const* exitRateVector);
+            template std::vector<storm::RationalNumber> SparseCtmcCslHelper::computeLongRunAverageProbabilities(Environment const& env, storm::solver::SolveGoal<storm::RationalNumber>&& goal, storm::storage::SparseMatrix<storm::RationalNumber> const& rateMatrix, storm::storage::BitVector const& psiStates, std::vector<storm::RationalNumber> const* exitRateVector);
+            template std::vector<storm::RationalFunction> SparseCtmcCslHelper::computeLongRunAverageProbabilities(Environment const& env, storm::solver::SolveGoal<storm::RationalFunction>&& goal, storm::storage::SparseMatrix<storm::RationalFunction> const& rateMatrix, storm::storage::BitVector const& psiStates, std::vector<storm::RationalFunction> const* exitRateVector);
             
-            template std::vector<storm::RationalNumber> SparseCtmcCslHelper::computeLongRunAverageRewards(Environment const& env, storm::solver::SolveGoal<storm::RationalNumber>&& goal, storm::storage::SparseMatrix<storm::RationalNumber> const& probabilityMatrix, storm::models::sparse::StandardRewardModel<RationalNumber> const& rewardModel, std::vector<storm::RationalNumber> const* exitRateVector);
-            template std::vector<storm::RationalFunction> SparseCtmcCslHelper::computeLongRunAverageRewards(Environment const& env, storm::solver::SolveGoal<storm::RationalFunction>&& goal, storm::storage::SparseMatrix<storm::RationalFunction> const& probabilityMatrix, storm::models::sparse::StandardRewardModel<RationalFunction> const& rewardModel, std::vector<storm::RationalFunction> const* exitRateVector);
+            template std::vector<storm::RationalNumber> SparseCtmcCslHelper::computeLongRunAverageRewards(Environment const& env, storm::solver::SolveGoal<storm::RationalNumber>&& goal, storm::storage::SparseMatrix<storm::RationalNumber> const& rateMatrix, storm::models::sparse::StandardRewardModel<RationalNumber> const& rewardModel, std::vector<storm::RationalNumber> const* exitRateVector);
+            template std::vector<storm::RationalFunction> SparseCtmcCslHelper::computeLongRunAverageRewards(Environment const& env, storm::solver::SolveGoal<storm::RationalFunction>&& goal, storm::storage::SparseMatrix<storm::RationalFunction> const& rateMatrix, storm::models::sparse::StandardRewardModel<RationalFunction> const& rewardModel, std::vector<storm::RationalFunction> const* exitRateVector);
 
-            template std::vector<storm::RationalNumber> SparseCtmcCslHelper::computeLongRunAverageRewards(Environment const& env, storm::solver::SolveGoal<storm::RationalNumber>&& goal, storm::storage::SparseMatrix<storm::RationalNumber> const& probabilityMatrix, std::vector<storm::RationalNumber> const& stateRewardVector, std::vector<storm::RationalNumber> const* exitRateVector);
-            template std::vector<storm::RationalFunction> SparseCtmcCslHelper::computeLongRunAverageRewards(Environment const& env, storm::solver::SolveGoal<storm::RationalFunction>&& goal, storm::storage::SparseMatrix<storm::RationalFunction> const& probabilityMatrix, std::vector<storm::RationalFunction> const& stateRewardVector, std::vector<storm::RationalFunction> const* exitRateVector);
+            template std::vector<storm::RationalNumber> SparseCtmcCslHelper::computeLongRunAverageRewards(Environment const& env, storm::solver::SolveGoal<storm::RationalNumber>&& goal, storm::storage::SparseMatrix<storm::RationalNumber> const& rateMatrix, std::vector<storm::RationalNumber> const& stateRewardVector, std::vector<storm::RationalNumber> const* exitRateVector);
+            template std::vector<storm::RationalFunction> SparseCtmcCslHelper::computeLongRunAverageRewards(Environment const& env, storm::solver::SolveGoal<storm::RationalFunction>&& goal, storm::storage::SparseMatrix<storm::RationalFunction> const& rateMatrix, std::vector<storm::RationalFunction> const& stateRewardVector, std::vector<storm::RationalFunction> const* exitRateVector);
 
             template std::vector<storm::RationalNumber> SparseCtmcCslHelper::computeCumulativeRewards(Environment const& env, storm::solver::SolveGoal<storm::RationalNumber>&& goal, storm::storage::SparseMatrix<storm::RationalNumber> const& rateMatrix, std::vector<storm::RationalNumber> const& exitRateVector, storm::models::sparse::StandardRewardModel<storm::RationalNumber> const& rewardModel, double timeBound);
             template std::vector<storm::RationalFunction> SparseCtmcCslHelper::computeCumulativeRewards(Environment const& env, storm::solver::SolveGoal<storm::RationalFunction>&& goal, storm::storage::SparseMatrix<storm::RationalFunction> const& rateMatrix, std::vector<storm::RationalFunction> const& exitRateVector, storm::models::sparse::StandardRewardModel<storm::RationalFunction> const& rewardModel, double timeBound);

src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.h

@@ -14,6 +14,10 @@ namespace storm {
     
     class Environment;
     
+    namespace storage {
+        class StronglyConnectedComponent;
+    }
+    
     namespace modelchecker {
         namespace helper {
             class SparseCtmcCslHelper {
@@ -52,13 +56,13 @@ namespace storm {
                 static std::vector<ValueType> computeTotalRewards(Environment const& env, storm::solver::SolveGoal<ValueType>&& goal, storm::storage::SparseMatrix<ValueType> const& rateMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::vector<ValueType> const& exitRateVector, RewardModelType const& rewardModel, bool qualitative);
 
                 template <typename ValueType>
-                static std::vector<ValueType> computeLongRunAverageProbabilities(Environment const& env, storm::solver::SolveGoal<ValueType>&& goal, storm::storage::SparseMatrix<ValueType> const& probabilityMatrix, storm::storage::BitVector const& psiStates, std::vector<ValueType> const* exitRateVector);
+                static std::vector<ValueType> computeLongRunAverageProbabilities(Environment const& env, storm::solver::SolveGoal<ValueType>&& goal, storm::storage::SparseMatrix<ValueType> const& rateMatrix, storm::storage::BitVector const& psiStates, std::vector<ValueType> const* exitRateVector);
 
                 template <typename ValueType, typename RewardModelType>
-                static std::vector<ValueType> computeLongRunAverageRewards(Environment const& env, storm::solver::SolveGoal<ValueType>&& goal, storm::storage::SparseMatrix<ValueType> const& probabilityMatrix, RewardModelType const& rewardModel, std::vector<ValueType> const* exitRateVector);
+                static std::vector<ValueType> computeLongRunAverageRewards(Environment const& env, storm::solver::SolveGoal<ValueType>&& goal, storm::storage::SparseMatrix<ValueType> const& rateMatrix, RewardModelType const& rewardModel, std::vector<ValueType> const* exitRateVector);
 
                 template <typename ValueType>
-                static std::vector<ValueType> computeLongRunAverageRewards(Environment const& env, storm::solver::SolveGoal<ValueType>&& goal, storm::storage::SparseMatrix<ValueType> const& probabilityMatrix, std::vector<ValueType> const& stateRewardVector, std::vector<ValueType> const* exitRateVector);
+                static std::vector<ValueType> computeLongRunAverageRewards(Environment const& env, storm::solver::SolveGoal<ValueType>&& goal, storm::storage::SparseMatrix<ValueType> const& rateMatrix, std::vector<ValueType> const& stateRewardVector, std::vector<ValueType> const* exitRateVector);
 
                 template <typename ValueType>
                 static std::vector<ValueType> computeReachabilityTimes(Environment const& env, storm::solver::SolveGoal<ValueType>&& goal, storm::storage::SparseMatrix<ValueType> const& rateMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::vector<ValueType> const& exitRateVector, storm::storage::BitVector const& targetStates, bool qualitative);
@@ -119,7 +123,16 @@ namespace storm {
                 
             private:
                 template <typename ValueType>
-                static std::vector<ValueType> computeLongRunAverages(Environment const& env, storm::solver::SolveGoal<ValueType>&& goal, storm::storage::SparseMatrix<ValueType> const& probabilityMatrix, std::function<ValueType (storm::storage::sparse::state_type const& state)> const& valueGetter, std::vector<ValueType> const* exitRateVector);
+                static std::vector<ValueType> computeLongRunAverages(Environment const& env, storm::solver::SolveGoal<ValueType>&& goal, storm::storage::SparseMatrix<ValueType> const& rateMatrix, std::function<ValueType (storm::storage::sparse::state_type const& state)> const& valueGetter, std::vector<ValueType> const* exitRateVector);
+                template <typename ValueType>
+                static ValueType computeLongRunAveragesForBscc(Environment const& env, storm::storage::StronglyConnectedComponent const& bscc, storm::storage::SparseMatrix<ValueType> const& rateMatrix, std::function<ValueType (storm::storage::sparse::state_type const& state)> const& valueGetter, std::vector<ValueType> const* exitRateVector);
+                template <typename ValueType>
+                static ValueType computeLongRunAveragesForBsccVi(Environment const& env, storm::storage::StronglyConnectedComponent const& bscc, storm::storage::SparseMatrix<ValueType> const& rateMatrix, std::function<ValueType (storm::storage::sparse::state_type const& state)> const& valueGetter, std::vector<ValueType> const* exitRateVector);
+                template <typename ValueType>
+                static std::pair<ValueType, std::vector<ValueType>> computeLongRunAveragesForBsccGainBias(Environment const& env, storm::storage::StronglyConnectedComponent const& bscc, storm::storage::SparseMatrix<ValueType> const& rateMatrix, std::function<ValueType (storm::storage::sparse::state_type const& state)> const& valueGetter, std::vector<ValueType> const* exitRateVector);
+                template <typename ValueType>
+                static std::pair<ValueType, std::vector<ValueType>> computeLongRunAveragesForBsccLraDistr(Environment const& env, storm::storage::StronglyConnectedComponent const& bscc, storm::storage::SparseMatrix<ValueType> const& rateMatrix, std::function<ValueType (storm::storage::sparse::state_type const& state)> const& valueGetter, std::vector<ValueType> const* exitRateVector);
+                
             };
         }
     }

src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp

@@ -13,6 +13,9 @@
 
 #include "storm/environment/Environment.h"
 #include "storm/environment/solver/MinMaxSolverEnvironment.h"
+#include "storm/environment/solver/TopologicalSolverEnvironment.h"
+#include "storm/environment/solver/LongRunAverageSolverEnvironment.h"
+#include "storm/environment/solver/EigenSolverEnvironment.h"
 
 #include "storm/utility/macros.h"
 #include "storm/utility/vector.h"
@@ -703,7 +706,9 @@ namespace storm {
                 auto underlyingSolverEnvironment = env;
                 if (env.solver().isForceSoundness()) {
                     // For sound computations, the error in the MECS plus the error in the remaining system should be less then the user defined precsion.
-                    underlyingSolverEnvironment.solver().minMax().setPrecision(env.solver().minMax().getPrecision() / storm::utility::convertNumber<storm::RationalNumber>(2));
+                    underlyingSolverEnvironment.solver().minMax().setPrecision(env.solver().lra().getPrecision() / storm::utility::convertNumber<storm::RationalNumber>(2));
+                    underlyingSolverEnvironment.solver().minMax().setRelativeTerminationCriterion(env.solver().lra().getRelativeTerminationCriterion());
+                    underlyingSolverEnvironment.solver().lra().setPrecision(env.solver().lra().getPrecision() / storm::utility::convertNumber<storm::RationalNumber>(2));
                 }
                 
                 for (uint64_t currentMecIndex = 0; currentMecIndex < mecDecomposition.size(); ++currentMecIndex) {
@@ -877,11 +882,11 @@ namespace storm {
                 }
                 
                 // Solve MEC with the method specified in the settings
-                storm::solver::LraMethod method = env.solver().minMax().getLraMethod();
-                if (storm::NumberTraits<ValueType>::IsExact && env.solver().minMax().isLraMethodSetFromDefault() && method != storm::solver::LraMethod::LinearProgramming) {
+                storm::solver::LraMethod method = env.solver().lra().getNondetLraMethod();
+                if (storm::NumberTraits<ValueType>::IsExact && env.solver().lra().isNondetLraMethodSetFromDefault() && method != storm::solver::LraMethod::LinearProgramming) {
                     STORM_LOG_INFO("Selecting 'LP' as the solution technique for long-run properties to guarantee exact results. If you want to override this, please explicitly specify a different LRA method.");
                     method = storm::solver::LraMethod::LinearProgramming;
-                } else if (env.solver().isForceSoundness() && env.solver().minMax().isLraMethodSetFromDefault() && method != storm::solver::LraMethod::ValueIteration) {
+                } else if (env.solver().isForceSoundness() && env.solver().lra().isNondetLraMethodSetFromDefault() && method != storm::solver::LraMethod::ValueIteration) {
                     STORM_LOG_INFO("Selecting 'VI' as the solution technique for long-run properties to guarantee sound results. If you want to override this, please explicitly specify a different LRA method.");
                     method = storm::solver::LraMethod::ValueIteration;
                 }
@@ -982,7 +987,7 @@ namespace storm {
                 ValueType uniformizationRate = storm::utility::vector::max_if(exitRateVector, markovianMecStates);
                 // To ensure that the model is aperiodic, we need to make sure that every Markovian state gets a self loop.
                 // Hence, we increase the uniformization rate a little.
-                uniformizationRate += storm::utility::one<ValueType>(); // Todo: try other values such as *=1.01
+                uniformizationRate *= (storm::utility::one<ValueType>() + storm::utility::convertNumber<ValueType>(env.solver().lra().getAperiodicFactor()));
 
                 // Get the transitions of the submodel, that is
                 // * a matrix aMarkovian with all (uniformized) transitions from Markovian mec states to all Markovian mec states.
@@ -1042,24 +1047,35 @@ namespace storm {
                 
                 // start the iterations
                 
-                ValueType precision = storm::utility::convertNumber<ValueType>(env.solver().minMax().getPrecision()) / uniformizationRate;
-                bool relative = env.solver().minMax().getRelativeTerminationCriterion();
+                ValueType precision = storm::utility::convertNumber<ValueType>(env.solver().lra().getPrecision()) / uniformizationRate;
+                bool relative = env.solver().lra().getRelativeTerminationCriterion();
                 std::vector<ValueType> v(aMarkovian.getRowCount(), storm::utility::zero<ValueType>());
                 std::vector<ValueType> w = v;
                 std::vector<ValueType> x, b;
+                auto solverEnv = env;
                 std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver;
                 if (hasProbabilisticStates) {
+                    if (env.solver().isForceSoundness()) {
+                        // To get correct results, the inner equation systems are solved exactly.
+                        // TODO investigate how an error would propagate
+                        solverEnv.solver().minMax().setMethod(storm::solver::MinMaxMethod::Topological);
+                        solverEnv.solver().topological().setUnderlyingMinMaxMethod(storm::solver::MinMaxMethod::PolicyIteration);
+                        solverEnv.solver().setLinearEquationSolverType(storm::solver::EquationSolverType::Topological);
+                        solverEnv.solver().topological().setUnderlyingEquationSolverType(storm::solver::EquationSolverType::Eigen);
+                        solverEnv.solver().eigen().setMethod(storm::solver::EigenLinearEquationSolverMethod::SparseLU);
+                    }
+                    
                     x.resize(aProbabilistic.getRowGroupCount(), storm::utility::zero<ValueType>());
                     b = probabilisticChoiceRewards;
                     
                     // Check for requirements of the solver.
                     // The solution is unique as we assume non-zeno MAs.
                     storm::solver::GeneralMinMaxLinearEquationSolverFactory<ValueType> minMaxLinearEquationSolverFactory;
-                    storm::solver::MinMaxLinearEquationSolverRequirements requirements = minMaxLinearEquationSolverFactory.getRequirements(env, true, true, dir);
+                    storm::solver::MinMaxLinearEquationSolverRequirements requirements = minMaxLinearEquationSolverFactory.getRequirements(solverEnv, true, true, dir);
                     requirements.clearLowerBounds();
                     STORM_LOG_THROW(!requirements.hasEnabledCriticalRequirement(), storm::exceptions::UncheckedRequirementException, "Solver requirements " + requirements.getEnabledRequirementsAsString() + " not checked.");
     
-                    solver = minMaxLinearEquationSolverFactory.create(env, std::move(aProbabilistic));
+                    solver = minMaxLinearEquationSolverFactory.create(solverEnv, std::move(aProbabilistic));
                     solver->setLowerBound(storm::utility::zero<ValueType>());
                     solver->setHasUniqueSolution(true);
                     solver->setHasNoEndComponents(true);
@@ -1067,10 +1083,16 @@ namespace storm {
                     solver->setCachingEnabled(true);
                 }
                 
-                while (true) {
+                uint64_t iter = 0;
+                boost::optional<uint64_t> maxIter;
+                if (env.solver().lra().isMaximalIterationCountSet()) {
+                    maxIter = env.solver().lra().getMaximalIterationCount();
+                }
+                while (!maxIter.is_initialized() || iter < maxIter.get()) {
+                    ++iter;
                     // Compute the expected total rewards for the probabilistic states
                     if (hasProbabilisticStates) {
-                        solver->solveEquations(env, dir, x, b);
+                        solver->solveEquations(solverEnv, dir, x, b);
                     }
                     // now compute the values for the markovian states. We also keep track of the maximal and minimal difference between two values (for convergence checking)
                     auto vIt = v.begin();
@@ -1106,8 +1128,12 @@ namespace storm {
                         storm::utility::vector::addVectors(b, probabilisticChoiceRewards, b);
                     }
                 }
+                if (maxIter.is_initialized() && iter == maxIter.get()) {
+                    STORM_LOG_WARN("LRA computation did not converge within " << iter << " iterations.");
+                } else {
+                    STORM_LOG_TRACE("LRA computation converged after " << iter << " iterations.");
+                }
                 return v.front() * uniformizationRate;
-            
             }
             
             template std::vector<double> SparseMarkovAutomatonCslHelper::computeBoundedUntilProbabilities(Environment const& env, OptimizationDirection dir, storm::storage::SparseMatrix<double> const& transitionMatrix, std::vector<double> const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, std::pair<double, double> const& boundsPair);

src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp

@@ -38,6 +38,7 @@
 #include "storm/transformer/EndComponentEliminator.h"
 
 #include "storm/environment/solver/MinMaxSolverEnvironment.h"
+#include "storm/environment/solver/LongRunAverageSolverEnvironment.h"
 
 #include "storm/exceptions/InvalidStateException.h"
 #include "storm/exceptions/InvalidPropertyException.h"
@@ -1251,8 +1252,10 @@ namespace storm {
                 
                 auto underlyingSolverEnvironment = env;
                 if (env.solver().isForceSoundness()) {
-                    // For sound computations, the error in the MECS plus the error in the remaining system should be less then the user defined precsion.
-                    underlyingSolverEnvironment.solver().minMax().setPrecision(env.solver().minMax().getPrecision() / storm::utility::convertNumber<storm::RationalNumber>(2));
+                    // For sound computations, the error in the MECS plus the error in the remaining system should be less than the user defined precsion.
+                    underlyingSolverEnvironment.solver().lra().setPrecision(env.solver().lra().getPrecision() / storm::utility::convertNumber<storm::RationalNumber>(2));
+                    underlyingSolverEnvironment.solver().minMax().setPrecision(env.solver().lra().getPrecision() / storm::utility::convertNumber<storm::RationalNumber>(2));
+                    underlyingSolverEnvironment.solver().minMax().setRelativeTerminationCriterion(env.solver().lra().getRelativeTerminationCriterion());
                 }
                 
                 for (uint_fast64_t currentMecIndex = 0; currentMecIndex < mecDecomposition.size(); ++currentMecIndex) {
@@ -1473,11 +1476,11 @@ namespace storm {
                 }
                 
                 // Solve MEC with the method specified in the settings
-                storm::solver::LraMethod method = env.solver().minMax().getLraMethod();
-                if (storm::NumberTraits<ValueType>::IsExact && env.solver().minMax().isLraMethodSetFromDefault() && method != storm::solver::LraMethod::LinearProgramming) {
+                storm::solver::LraMethod method = env.solver().lra().getNondetLraMethod();
+                if (storm::NumberTraits<ValueType>::IsExact && env.solver().lra().isNondetLraMethodSetFromDefault() && method != storm::solver::LraMethod::LinearProgramming) {
                     STORM_LOG_INFO("Selecting 'LP' as the solution technique for long-run properties to guarantee exact results. If you want to override this, please explicitly specify a different LRA method.");
                     method = storm::solver::LraMethod::LinearProgramming;
-                } else if (env.solver().isForceSoundness() && env.solver().minMax().isLraMethodSetFromDefault() && method != storm::solver::LraMethod::ValueIteration) {
+                } else if (env.solver().isForceSoundness() && env.solver().lra().isNondetLraMethodSetFromDefault() && method != storm::solver::LraMethod::ValueIteration) {
                     STORM_LOG_INFO("Selecting 'VI' as the solution technique for long-run properties to guarantee sound results. If you want to override this, please explicitly specify a different LRA method.");
                     method = storm::solver::LraMethod::ValueIteration;
                 }
@@ -1518,7 +1521,7 @@ namespace storm {
                 std::vector<ValueType> choiceRewards;
                 choiceRewards.reserve(mecChoices.getNumberOfSetBits());
                 uint64_t currRow = 0;
-                ValueType selfLoopProb = storm::utility::convertNumber<ValueType>(0.1); // todo try other values
+                ValueType selfLoopProb = storm::utility::convertNumber<ValueType>(env.solver().lra().getAperiodicFactor());
                 ValueType scalingFactor = storm::utility::one<ValueType>() - selfLoopProb;
                 for (auto const& mecState : mecStates) {
                     mecTransitionBuilder.newRowGroup(currRow);
@@ -1553,14 +1556,21 @@ namespace storm {
                 STORM_LOG_ASSERT(mecTransitions.isProbabilistic(), "The MEC-Matrix is not probabilistic.");
                 
                 // start the iterations
-                ValueType precision = storm::utility::convertNumber<ValueType>(env.solver().minMax().getPrecision()) / scalingFactor;
-                bool relative = env.solver().minMax().getRelativeTerminationCriterion();
+                ValueType precision = storm::utility::convertNumber<ValueType>(env.solver().lra().getPrecision()) / scalingFactor;
+                bool relative = env.solver().lra().getRelativeTerminationCriterion();
                 std::vector<ValueType> x(mecTransitions.getRowGroupCount(), storm::utility::zero<ValueType>());
                 std::vector<ValueType> xPrime = x;
                 
                 auto multiplier = storm::solver::MultiplierFactory<ValueType>().create(env, mecTransitions);
                 ValueType maxDiff, minDiff;
-                while (true) {
+                
+                uint64_t iter = 0;
+                boost::optional<uint64_t> maxIter;
+                if (env.solver().lra().isMaximalIterationCountSet()) {
+                    maxIter = env.solver().lra().getMaximalIterationCount();
+                }
+                while (!maxIter.is_initialized() || iter < maxIter.get()) {
+                    ++iter;
                     // Compute the obtained rewards for the next step
                     multiplier->multiplyAndReduce(env, dir, x, &choiceRewards, x);
                     
@@ -1585,6 +1595,11 @@ namespace storm {
                         break;
                     }
                 }
+                if (maxIter.is_initialized() && iter == maxIter.get()) {
+                    STORM_LOG_WARN("LRA computation did not converge within " << iter << " iterations.");
+                } else {
+                    STORM_LOG_TRACE("LRA computation converged after " << iter << " iterations.");
+                }
                 
                 if (scheduler) {
                     std::vector<uint_fast64_t> localMecChoices(mecTransitions.getRowGroupCount(), 0);

src/storm/models/symbolic/Ctmc.cpp

@@ -84,6 +84,17 @@ namespace storm {
                 return exitRates.get();
             }
             
+            template<storm::dd::DdType Type, typename ValueType>
+            void Ctmc<Type, ValueType>::reduceToStateBasedRewards() {
+                for (auto& rewardModel : this->getRewardModels()) {
+                    if (rewardModel.second.hasStateActionRewards()) {
+                        rewardModel.second.getStateActionRewardVector() *= getExitRateVector();
+                    }
+                    rewardModel.second.reduceToStateBasedRewards(this->getTransitionMatrix(), this->getRowVariables(), this->getColumnVariables(), true);
+                }
+            }
+            
+            
             template<storm::dd::DdType Type, typename ValueType>
             template<typename NewValueType>
             std::shared_ptr<Ctmc<Type, NewValueType>> Ctmc<Type, ValueType>::toValueType() const {

src/storm/models/symbolic/Ctmc.h

@@ -141,6 +141,8 @@ namespace storm {
                  */
                 storm::dd::Add<Type, ValueType> const& getExitRateVector() const;
                 
+                virtual void reduceToStateBasedRewards() override;
+
                 template<typename NewValueType>
                 std::shared_ptr<Ctmc<Type, NewValueType>> toValueType() const;
                 

src/storm/models/symbolic/DeterministicModel.cpp

@@ -45,13 +45,6 @@ namespace storm {
                 // Intentionally left empty.
             }
             
-            template<storm::dd::DdType Type, typename ValueType>
-            void DeterministicModel<Type, ValueType>::reduceToStateBasedRewards() {
-                for (auto& rewardModel : this->getRewardModels()) {
-                    rewardModel.second.reduceToStateBasedRewards(this->getTransitionMatrix(), this->getRowVariables(), this->getColumnVariables(), true);
-                }
-            }
-            
             // Explicitly instantiate the template class.
             template class DeterministicModel<storm::dd::DdType::CUDD>;
             template class DeterministicModel<storm::dd::DdType::Sylvan>;

src/storm/models/symbolic/DeterministicModel.h

@@ -81,7 +81,6 @@ namespace storm {
                                    std::map<std::string, storm::dd::Bdd<Type>> labelToBddMap = std::map<std::string, storm::dd::Bdd<Type>>(),
                                    std::unordered_map<std::string, RewardModelType> const& rewardModels = std::unordered_map<std::string, RewardModelType>());
                 
-                virtual void reduceToStateBasedRewards() override;
             };
             
         } // namespace symbolic

src/storm/models/symbolic/Dtmc.cpp

@@ -43,6 +43,13 @@ namespace storm {
                 // Intentionally left empty.
             }
             
+            template<storm::dd::DdType Type, typename ValueType>
+            void Dtmc<Type, ValueType>::reduceToStateBasedRewards() {
+                for (auto& rewardModel : this->getRewardModels()) {
+                    rewardModel.second.reduceToStateBasedRewards(this->getTransitionMatrix(), this->getRowVariables(), this->getColumnVariables(), true);
+                }
+            }
+            
             template<storm::dd::DdType Type, typename ValueType>
             template<typename NewValueType>
             std::shared_ptr<Dtmc<Type, NewValueType>> Dtmc<Type, ValueType>::toValueType() const {

src/storm/models/symbolic/Dtmc.h

@@ -77,6 +77,8 @@ namespace storm {
                      std::map<std::string, storm::dd::Bdd<Type>> labelToBddMap = std::map<std::string, storm::dd::Bdd<Type>>(),
                      std::unordered_map<std::string, RewardModelType> const& rewardModels = std::unordered_map<std::string, RewardModelType>());
                 
+                virtual void reduceToStateBasedRewards() override;
+
                 template<typename NewValueType>
                 std::shared_ptr<Dtmc<Type, NewValueType>> toValueType() const;
 

src/storm/settings/Argument.cpp

@@ -171,20 +171,28 @@ namespace storm {
             out << std::setw(0) << std::left << "<" << this->getName() << ">";
             if (!this->validators.empty() || this->hasDefaultValue) {
                 out << " (";
+                bool previousEntry = false;
+                if (this->getIsOptional()) {
+                    out << "optional";
+                    previousEntry = true;
+                }
                 if (!this->validators.empty()) {
+                    if (previousEntry) {
+                        out << "; ";
+                    }
                     for (uint64_t i = 0; i < this->validators.size(); ++i) {
                         out << this->validators[i]->toString();
                         if (i + 1 < this->validators.size()) {
                             out << ", ";
                         }
                     }
-                    
-                    if (this->hasDefaultValue) {
-                        out << "; ";
-                    }
+                    previousEntry = true;
                 }
                 
                 if (this->hasDefaultValue) {
+                    if (previousEntry) {
+                        out << "; ";
+                    }
                     out << "default: ";
                     printValue(out, defaultValue);
                 }

src/storm/settings/ArgumentBuilder.h

@@ -89,13 +89,12 @@ namespace storm {
             }
             
             /*!
-             * Sets whether the argument is to be optional.
+             * Make the argument optional.
              *
-             * @param isOptional A flag that indicates whether the argument is to be optional.
              * @return A reference to the argument builder.
              */
-            ArgumentBuilder& setIsOptional(bool isOptional) {
-                this->isOptional = isOptional;
+            ArgumentBuilder& makeOptional() {
+                this->isOptional = true;
                 STORM_LOG_THROW(this->hasDefaultValue, storm::exceptions::IllegalFunctionCallException, "Unable to make argument '" << this->name << "' optional without default value.");
                 return *this;
             }

src/storm/settings/SettingsManager.cpp

@@ -24,6 +24,7 @@
 #include "storm/settings/modules/GmmxxEquationSolverSettings.h"
 #include "storm/settings/modules/NativeEquationSolverSettings.h"
 #include "storm/settings/modules/EliminationSettings.h"
+#include "storm/settings/modules/LongRunAverageSolverSettings.h"
 #include "storm/settings/modules/MinMaxEquationSolverSettings.h"
 #include "storm/settings/modules/GameSolverSettings.h"
 #include "storm/settings/modules/BisimulationSettings.h"
@@ -490,7 +491,7 @@ namespace storm {
             // In case there are optional arguments that were not set, we set them to their default value.
             for (uint_fast64_t i = argumentCache.size(); i < option->getArgumentCount(); ++i) {
                 ArgumentBase& argument = option->getArgument(i);
-                STORM_LOG_THROW(argument.getHasDefaultValue() || argument.getIsOptional(), storm::exceptions::OptionParserException, "Non-optional argument <" << argument.getName() << "> of option:\n" << *option);
+                STORM_LOG_THROW(argument.getIsOptional(), storm::exceptions::OptionParserException, "Non-optional argument <" << argument.getName() << "> of option:\n" << *option);
                 argument.setFromDefaultValue();
             }
             
@@ -654,6 +655,7 @@ namespace storm {
             storm::settings::addModule<storm::settings::modules::EigenEquationSolverSettings>();
             storm::settings::addModule<storm::settings::modules::NativeEquationSolverSettings>();
             storm::settings::addModule<storm::settings::modules::EliminationSettings>();
+            storm::settings::addModule<storm::settings::modules::LongRunAverageSolverSettings>();
             storm::settings::addModule<storm::settings::modules::MinMaxEquationSolverSettings>();
             storm::settings::addModule<storm::settings::modules::GameSolverSettings>();
             storm::settings::addModule<storm::settings::modules::BisimulationSettings>();

src/storm/settings/modules/EliminationSettings.cpp

@@ -28,7 +28,7 @@ namespace storm {
                 
                 this->addOption(storm::settings::OptionBuilder(moduleName, entryStatesLastOptionName, true, "Sets whether the entry states are eliminated last.").setIsAdvanced().build());
                 this->addOption(storm::settings::OptionBuilder(moduleName, maximalSccSizeOptionName, true, "Sets the maximal size of the SCCs for which state elimination is applied.").setIsAdvanced()
-                                .addArgument(storm::settings::ArgumentBuilder::createUnsignedIntegerArgument("maxsize", "The maximal size of an SCC on which state elimination is applied.").setDefaultValueUnsignedInteger(20).setIsOptional(true).build()).build());
+                                .addArgument(storm::settings::ArgumentBuilder::createUnsignedIntegerArgument("maxsize", "The maximal size of an SCC on which state elimination is applied.").setDefaultValueUnsignedInteger(20).build()).build());
                 this->addOption(storm::settings::OptionBuilder(moduleName, useDedicatedModelCheckerOptionName, true, "Sets whether to use the dedicated model elimination checker (only DTMCs).").setIsAdvanced().build());
             }
             

src/storm/settings/modules/GeneralSettings.cpp

@@ -35,7 +35,7 @@ namespace storm {
 
             GeneralSettings::GeneralSettings() : ModuleSettings(moduleName) {
                 this->addOption(storm::settings::OptionBuilder(moduleName, helpOptionName, false, "Shows available options, arguments and descriptions.").setShortName(helpOptionShortName)
-                                .addArgument(storm::settings::ArgumentBuilder::createStringArgument("filter", "'frequent' for frequently used options, 'all' for the complete help, or a regular expression to show help for all matching entities.").setDefaultValueString("frequent").build()).build());
+                                .addArgument(storm::settings::ArgumentBuilder::createStringArgument("filter", "'frequent' for frequently used options, 'all' for the complete help, or a regular expression to show help for all matching entities.").setDefaultValueString("frequent").makeOptional().build()).build());
                 this->addOption(storm::settings::OptionBuilder(moduleName, versionOptionName, false, "Prints the version information.").build());
                 this->addOption(storm::settings::OptionBuilder(moduleName, verboseOptionName, false, "Enables more verbose output.").setShortName(verboseOptionShortName).build());
                 this->addOption(storm::settings::OptionBuilder(moduleName, showProgressOptionName, false, "Sets when additional information (if available) about the progress is printed.").addArgument(storm::settings::ArgumentBuilder::createUnsignedIntegerArgument("delay", "The delay to wait (in seconds) between emitting information (0 means never print progress).").setDefaultValueUnsignedInteger(5).build()).build());

src/storm/settings/modules/IOSettings.cpp

@@ -25,6 +25,7 @@ namespace storm {
             const std::string IOSettings::exportCdfOptionName = "exportcdf";
             const std::string IOSettings::exportCdfOptionShortName = "cdf";
             const std::string IOSettings::exportSchedulerOptionName = "exportscheduler";
+            const std::string IOSettings::exportMonotonicityName = "exportmonotonicity";
             const std::string IOSettings::explicitOptionName = "explicit";
             const std::string IOSettings::explicitOptionShortName = "exp";
             const std::string IOSettings::explicitDrnOptionName = "explicit-drn";
@@ -59,6 +60,7 @@ namespace storm {
                                         .addArgument(storm::settings::ArgumentBuilder::createStringArgument("filename", "The name of the file to which the model is to be written.").build()).build());
                 this->addOption(storm::settings::OptionBuilder(moduleName, exportCdfOptionName, false, "Exports the cumulative density function for reward bounded properties into a .csv file.").setIsAdvanced().setShortName(exportCdfOptionShortName).addArgument(storm::settings::ArgumentBuilder::createStringArgument("directory", "A path to an existing directory where the cdf files will be stored.").build()).build());
                 this->addOption(storm::settings::OptionBuilder(moduleName, exportSchedulerOptionName, false, "Exports the choices of an optimal scheduler to the given file (if supported by engine).").setIsAdvanced().addArgument(storm::settings::ArgumentBuilder::createStringArgument("filename", "The output file.").build()).build());
+                this->addOption(storm::settings::OptionBuilder(moduleName, exportMonotonicityName, false, "Exports the result of monotonicity checking to the given file.").setIsAdvanced().addArgument(storm::settings::ArgumentBuilder::createStringArgument("filename", "The output file.").build()).build());
                 this->addOption(storm::settings::OptionBuilder(moduleName, exportExplicitOptionName, "", "If given, the loaded model will be written to the specified file in the drn format.")
                                 .addArgument(storm::settings::ArgumentBuilder::createStringArgument("filename", "the name of the file to which the model is to be writen.").build()).build());
                 this->addOption(storm::settings::OptionBuilder(moduleName, exportDdOptionName, "", "If given, the loaded model will be written to the specified file in the drdd format.")
@@ -79,7 +81,7 @@ namespace storm {
                 this->addOption(storm::settings::OptionBuilder(moduleName, prismToJaniOptionName, false, "If set, the input PRISM model is transformed to JANI.").setIsAdvanced().build());
                 this->addOption(storm::settings::OptionBuilder(moduleName, propertyOptionName, false, "Specifies the properties to be checked on the model.").setShortName(propertyOptionShortName)
                                         .addArgument(storm::settings::ArgumentBuilder::createStringArgument("property or filename", "The formula or the file containing the formulas.").build())
-                                        .addArgument(storm::settings::ArgumentBuilder::createStringArgument("filter", "The names of the properties to check.").setDefaultValueString("all").build())
+                                        .addArgument(storm::settings::ArgumentBuilder::createStringArgument("filter", "The names of the properties to check.").setDefaultValueString("all").makeOptional().build())
                                         .build());
 
                 this->addOption(storm::settings::OptionBuilder(moduleName, transitionRewardsOptionName, false, "If given, the transition rewards are read from this file and added to the explicit model. Note that this requires the model to be given as an explicit model (i.e., via --" + explicitOptionName + ").")
@@ -91,7 +93,7 @@ namespace storm {
                 this->addOption(storm::settings::OptionBuilder(moduleName, constantsOptionName, false, "Specifies the constant replacements to use in symbolic models. Note that this requires the model to be given as an symbolic model (i.e., via --" + prismInputOptionName + " or --" + janiInputOptionName + ").").setShortName(constantsOptionShortName)
                                 .addArgument(storm::settings::ArgumentBuilder::createStringArgument("values", "A comma separated list of constants and their value, e.g. a=1,b=2,c=3.").setDefaultValueString("").build()).build());
                 this->addOption(storm::settings::OptionBuilder(moduleName, janiPropertyOptionName, false, "Specifies the properties from the jani model (given by --" + janiInputOptionName + ") to be checked.").setShortName(janiPropertyOptionShortName)
-                                .addArgument(storm::settings::ArgumentBuilder::createStringArgument("values", "A comma separated list of properties to be checked").setDefaultValueString("").build()).build());
+                                .addArgument(storm::settings::ArgumentBuilder::createStringArgument("values", "A comma separated list of properties to be checked").setDefaultValueString("").makeOptional().build()).build());
 
                 this->addOption(storm::settings::OptionBuilder(moduleName, qvbsInputOptionName, false, "Selects a model from the Quantitative Verification Benchmark Set.").setShortName(qvbsInputOptionShortName)
                         .addArgument(storm::settings::ArgumentBuilder::createStringArgument("model", "The short model name as in the benchmark set.").build())
@@ -162,6 +164,14 @@ namespace storm {
             std::string IOSettings::getExportSchedulerFilename() const {
                 return this->getOption(exportSchedulerOptionName).getArgumentByName("filename").getValueAsString();
             }
+
+            bool IOSettings::isExportMonotonicitySet() const {
+                return this->getOption(exportMonotonicityName).getHasOptionBeenSet();
+            }
+
+            std::string IOSettings::getExportMonotonicityFilename() const {
+                return this->getOption(exportMonotonicityName).getArgumentByName("filename").getValueAsString();
+            }
             
             bool IOSettings::isExplicitSet() const {
                 return this->getOption(explicitOptionName).getHasOptionBeenSet();

src/storm/settings/modules/IOSettings.h

@@ -105,6 +105,16 @@ namespace storm {
                  * Retrieves a filename to which an optimal scheduler will be exported.
                  */
                  std::string getExportSchedulerFilename() const;
+
+                /*!
+                * Retrieves whether an optimal scheduler is to be exported
+                */
+                bool isExportMonotonicitySet() const;
+
+                /*!
+                 * Retrieves a filename to which an optimal scheduler will be exported.
+                 */
+                std::string getExportMonotonicityFilename() const;
                 
                 /*!
                  * Retrieves whether the explicit option was set.
@@ -339,6 +349,7 @@ namespace storm {
                 static const std::string exportCdfOptionName;
                 static const std::string exportCdfOptionShortName;
                 static const std::string exportSchedulerOptionName;
+                static const std::string exportMonotonicityName;
                 static const std::string explicitOptionName;
                 static const std::string explicitOptionShortName;
                 static const std::string explicitDrnOptionName;

src/storm/settings/modules/LongRunAverageSolverSettings.cpp

@@ -0,0 +1,101 @@
+#include "storm/settings/modules/LongRunAverageSolverSettings.h"
+
+#include "storm/settings/Option.h"
+#include "storm/settings/ArgumentBuilder.h"
+#include "storm/settings/OptionBuilder.h"
+
+#include "storm/utility/macros.h"
+#include "storm/exceptions/IllegalArgumentValueException.h"
+
+namespace storm {
+    namespace settings {
+        namespace modules {
+            
+            const std::string LongRunAverageSolverSettings::moduleName = "lra";
+            const std::string LongRunAverageSolverSettings::detLraMethodOptionName = "detmethod";
+            const std::string LongRunAverageSolverSettings::nondetLraMethodOptionName = "nondetmethod";
+            const std::string LongRunAverageSolverSettings::maximalIterationsOptionName = "maxiter";
+            const std::string LongRunAverageSolverSettings::maximalIterationsOptionShortName = "i";
+            const std::string LongRunAverageSolverSettings::precisionOptionName = "precision";
+            const std::string LongRunAverageSolverSettings::absoluteOptionName = "absolute";
+            const std::string LongRunAverageSolverSettings::aperiodicFactorOptionName = "aperiodicfactor";
+
+            LongRunAverageSolverSettings::LongRunAverageSolverSettings() : ModuleSettings(moduleName) {
+                
+                std::vector<std::string> detLraMethods = {"gb", "gain-bias-equations", "distr", "lra-distribution-equations", "vi", "value-iteration"};
+                this->addOption(storm::settings::OptionBuilder(moduleName, detLraMethodOptionName, true, "Sets which method is preferred for computing long run averages on deterministic models.").setIsAdvanced()
+                                .addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of a long run average computation method.").addValidatorString(ArgumentValidatorFactory::createMultipleChoiceValidator(detLraMethods)).setDefaultValueString("gb").build()).build());
+                
+                std::vector<std::string> nondetLraMethods = {"vi", "value-iteration", "linear-programming", "lp"};
+                this->addOption(storm::settings::OptionBuilder(moduleName, nondetLraMethodOptionName, true, "Sets which method is preferred for computing long run averages on models with nondeterminism.").setIsAdvanced()
+                                .addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of a long run average computation method.").addValidatorString(ArgumentValidatorFactory::createMultipleChoiceValidator(nondetLraMethods)).setDefaultValueString("vi").build()).build());
+                
+                this->addOption(storm::settings::OptionBuilder(moduleName, maximalIterationsOptionName, false, "The maximal number of iterations to perform before iterative solving is aborted.").setShortName(maximalIterationsOptionShortName).setIsAdvanced().addArgument(storm::settings::ArgumentBuilder::createUnsignedIntegerArgument("count", "The maximal iteration count.").build()).build());
+                
+                this->addOption(storm::settings::OptionBuilder(moduleName, precisionOptionName, false, "The precision used for detecting convergence of iterative methods.").setIsAdvanced().addArgument(storm::settings::ArgumentBuilder::createDoubleArgument("value", "The precision to achieve.").setDefaultValueDouble(1e-06).addValidatorDouble(ArgumentValidatorFactory::createDoubleRangeValidatorExcluding(0.0, 1.0)).build()).build());
+
+                this->addOption(storm::settings::OptionBuilder(moduleName, absoluteOptionName, false, "Sets whether the relative or the absolute error is considered for detecting convergence.").setIsAdvanced().build());
+
+                this->addOption(storm::settings::OptionBuilder(moduleName, aperiodicFactorOptionName, true, "If required by the selected method (e.g. vi), this factor controls how the system is made aperiodic").setIsAdvanced().addArgument(storm::settings::ArgumentBuilder::createDoubleArgument("value", "The factor.").setDefaultValueDouble(0.125).addValidatorDouble(ArgumentValidatorFactory::createDoubleRangeValidatorExcluding(0.0, 1.0)).build()).build());
+                
+            }
+            
+            storm::solver::LraMethod LongRunAverageSolverSettings::getDetLraMethod() const {
+                std::string lraMethodString = this->getOption(detLraMethodOptionName).getArgumentByName("name").getValueAsString();
+                if (lraMethodString == "gain-bias-equations" || lraMethodString == "gb") {
+                    return storm::solver::LraMethod::GainBiasEquations;
+                }
+                if (lraMethodString == "lra-distribution-equations" || lraMethodString == "distr") {
+                    return storm::solver::LraMethod::LraDistributionEquations;
+                }
+                if (lraMethodString == "value-iteration" || lraMethodString == "vi") {
+                    return storm::solver::LraMethod::ValueIteration;
+                }
+                STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown lra solving technique for deterministic models:'" << lraMethodString << "'.");
+            }
+            
+            bool LongRunAverageSolverSettings::isDetLraMethodSetFromDefaultValue() const {
+                return !this->getOption(detLraMethodOptionName).getArgumentByName("name").getHasBeenSet() || this->getOption(detLraMethodOptionName).getArgumentByName("name").wasSetFromDefaultValue();
+            }
+            
+            storm::solver::LraMethod LongRunAverageSolverSettings::getNondetLraMethod() const {
+                std::string lraMethodString = this->getOption(nondetLraMethodOptionName).getArgumentByName("name").getValueAsString();
+                if (lraMethodString == "value-iteration" || lraMethodString == "vi") {
+                    return storm::solver::LraMethod::ValueIteration;
+                } else if (lraMethodString == "linear-programming" || lraMethodString == "lp") {
+                    return storm::solver::LraMethod::LinearProgramming;
+                }
+                STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown lra solving technique for nondeterministic models:'" << lraMethodString << "'.");
+            }
+            
+            bool LongRunAverageSolverSettings::isNondetLraMethodSetFromDefaultValue() const {
+                return !this->getOption(nondetLraMethodOptionName).getArgumentByName("name").getHasBeenSet() || this->getOption(nondetLraMethodOptionName).getArgumentByName("name").wasSetFromDefaultValue();
+            }
+            
+            bool LongRunAverageSolverSettings::isMaximalIterationCountSet() const {
+                return this->getOption(maximalIterationsOptionName).getHasOptionBeenSet();
+            }
+            
+            uint_fast64_t LongRunAverageSolverSettings::getMaximalIterationCount() const {
+                return this->getOption(maximalIterationsOptionName).getArgumentByName("count").getValueAsUnsignedInteger();
+            }
+            
+            bool LongRunAverageSolverSettings::isPrecisionSet() const {
+                return this->getOption(precisionOptionName).getHasOptionBeenSet();
+            }
+            
+            double LongRunAverageSolverSettings::getPrecision() const {
+                return this->getOption(precisionOptionName).getArgumentByName("value").getValueAsDouble();
+            }
+            
+            bool LongRunAverageSolverSettings::isRelativePrecision() const {
+                return this->getOption(absoluteOptionName).getHasOptionBeenSet();
+            }
+            
+            double LongRunAverageSolverSettings::getAperiodicFactor() const {
+                return this->getOption(aperiodicFactorOptionName).getArgumentByName("value").getValueAsDouble();
+            }
+            
+        }
+    }
+}

src/storm/settings/modules/LongRunAverageSolverSettings.h

@@ -0,0 +1,92 @@
+#pragma once
+
+#include "storm-config.h"
+#include "storm/settings/modules/ModuleSettings.h"
+
+#include "storm/solver/SolverSelectionOptions.h"
+#include "storm/solver/MultiplicationStyle.h"
+
+namespace storm {
+    namespace settings {
+        namespace modules {
+            
+            /*!
+             * This class represents the LRA solver settings.
+             */
+            class LongRunAverageSolverSettings : public ModuleSettings {
+            public:
+                
+                LongRunAverageSolverSettings();
+              
+                /*!
+                 * Retrieves the selected long run average method for deterministic models.
+                 */
+                storm::solver::LraMethod getDetLraMethod() const;
+                
+                /*!
+                 * Retrieves whether the LraMethod for deterministic models was set from a default value.
+                 */
+                bool isDetLraMethodSetFromDefaultValue() const;
+                
+                /*!
+                 * Retrieves the selected long run average method for nondeterministic models.
+                 */
+                storm::solver::LraMethod getNondetLraMethod() const;
+                
+                /*!
+                 * Retrieves whether the LraMethod for nondeterministic models was set from a default value.
+                 */
+                bool isNondetLraMethodSetFromDefaultValue() const;
+                
+                /*!
+                 * Retrieves whether a maximal iteration count for iterative methods was set.
+                 */
+                bool isMaximalIterationCountSet() const;
+                
+                /*!
+                 * Retrieves the maximal number of iterations to perform until giving up on converging.
+                 *
+                 * @return The maximal iteration count.
+                 */
+                uint_fast64_t getMaximalIterationCount() const;
+                
+                /*!
+                 * Retrieves whether the precision has been set.
+                 *
+                 * @return True iff the precision has been set.
+                 */
+                bool isPrecisionSet() const;
+                
+                /*!
+                 * Retrieves the precision that is used for detecting convergence.
+                 *
+                 * @return The precision to use for detecting convergence.
+                 */
+                double getPrecision() const;
+                
+                /*!
+                 * Retrieves whether the convergence criterion has been set to relative.
+                 */
+                bool isRelativePrecision() const;
+                
+                /*!
+                 * Retrieves a factor that describes how the system is made aperiodic (if necessary by the method)
+                 */
+                double getAperiodicFactor() const;
+                
+                // The name of the module.
+                static const std::string moduleName;
+                
+            private:
+                static const std::string detLraMethodOptionName;
+                static const std::string nondetLraMethodOptionName;
+                static const std::string maximalIterationsOptionName;
+                static const std::string maximalIterationsOptionShortName;
+                static const std::string precisionOptionName;
+                static const std::string absoluteOptionName;
+                static const std::string aperiodicFactorOptionName;
+            };
+            
+        }
+    }
+}

src/storm/settings/modules/MinMaxEquationSolverSettings.cpp

@@ -17,7 +17,6 @@ namespace storm {
             const std::string MinMaxEquationSolverSettings::maximalIterationsOptionShortName = "i";
             const std::string MinMaxEquationSolverSettings::precisionOptionName = "precision";
             const std::string MinMaxEquationSolverSettings::absoluteOptionName = "absolute";
-            const std::string MinMaxEquationSolverSettings::lraMethodOptionName = "lramethod";
             const std::string MinMaxEquationSolverSettings::markovAutomatonBoundedReachabilityMethodOptionName = "mamethod";
             const std::string MinMaxEquationSolverSettings::valueIterationMultiplicationStyleOptionName = "vimult";
             const std::string MinMaxEquationSolverSettings::intervalIterationSymmetricUpdatesOptionName = "symmetricupdates";
@@ -33,10 +32,6 @@ namespace storm {
 
                 this->addOption(storm::settings::OptionBuilder(moduleName, absoluteOptionName, false, "Sets whether the relative or the absolute error is considered for detecting convergence.").setIsAdvanced().build());
 
-                std::vector<std::string> lraMethods = {"vi", "value-iteration", "linear-programming", "lp"};
-                this->addOption(storm::settings::OptionBuilder(moduleName, lraMethodOptionName, false, "Sets which method is preferred for computing long run averages.").setIsAdvanced()
-                                .addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of a long run average computation method.").addValidatorString(ArgumentValidatorFactory::createMultipleChoiceValidator(lraMethods)).setDefaultValueString("vi").build()).build());
-
                 std::vector<std::string> maMethods = {"imca", "unifplus"};
                 this->addOption(storm::settings::OptionBuilder(moduleName, markovAutomatonBoundedReachabilityMethodOptionName, false, "The method to use to solve bounded reachability queries on MAs.").setIsAdvanced().addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of the method to use.").addValidatorString(ArgumentValidatorFactory::createMultipleChoiceValidator(maMethods)).setDefaultValueString("unifplus").build()).build());
                 
@@ -103,20 +98,6 @@ namespace storm {
                 return this->getOption(absoluteOptionName).getHasOptionBeenSet() ? MinMaxEquationSolverSettings::ConvergenceCriterion::Absolute : MinMaxEquationSolverSettings::ConvergenceCriterion::Relative;
             }
             
-            storm::solver::LraMethod MinMaxEquationSolverSettings::getLraMethod() const {
-                std::string lraMethodString = this->getOption(lraMethodOptionName).getArgumentByName("name").getValueAsString();
-                if (lraMethodString == "value-iteration" || lraMethodString == "vi") {
-                    return storm::solver::LraMethod::ValueIteration;
-                } else if (lraMethodString == "linear-programming" || lraMethodString == "lp") {
-                    return storm::solver::LraMethod::LinearProgramming;
-                }
-                STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown lra solving technique '" << lraMethodString << "'.");
-            }
-            
-            bool MinMaxEquationSolverSettings::isLraMethodSetFromDefaultValue() const {
-                return !this->getOption(lraMethodOptionName).getArgumentByName("name").getHasBeenSet() || this->getOption(lraMethodOptionName).getArgumentByName("name").wasSetFromDefaultValue();
-            }
-            
             MinMaxEquationSolverSettings::MarkovAutomatonBoundedReachabilityMethod MinMaxEquationSolverSettings::getMarkovAutomatonBoundedReachabilityMethod() const {
                 std::string techniqueAsString = this->getOption(markovAutomatonBoundedReachabilityMethodOptionName).getArgumentByName("name").getValueAsString();
                 if (techniqueAsString == "imca") {

src/storm/settings/modules/MinMaxEquationSolverSettings.h

@@ -86,18 +86,6 @@ namespace storm {
                  */
                 ConvergenceCriterion getConvergenceCriterion() const;
                 
-                /*!
-                 * Retrieves the selected long run average method.
-                 *
-                 * @return The selected long run average method.
-                 */
-                storm::solver::LraMethod getLraMethod() const;
-                
-                /*!
-                 * Retrieves whether the LraMethod was set from a default value.
-                 */
-                bool isLraMethodSetFromDefaultValue() const;
-                
                 /*!
                  * Retrieves the method to be used for bounded reachability on MAs.
                  *
@@ -126,7 +114,6 @@ namespace storm {
                 static const std::string maximalIterationsOptionShortName;
                 static const std::string precisionOptionName;
                 static const std::string absoluteOptionName;
-                static const std::string lraMethodOptionName;
                 static const std::string markovAutomatonBoundedReachabilityMethodOptionName;
                 static const std::string valueIterationMultiplicationStyleOptionName;
                 static const std::string intervalIterationSymmetricUpdatesOptionName;

src/storm/solver/SolverSelectionOptions.cpp

@@ -49,9 +49,13 @@ namespace storm {
         std::string toString(LraMethod m) {
             switch(m) {
                 case LraMethod::LinearProgramming:
-                    return "linearprogramming";
+                    return "linear-programming";
                 case LraMethod::ValueIteration:
-                    return "valueiteration";
+                    return "value-iteration";
+                case LraMethod::LraDistributionEquations:
+                    return "lra-distribution-equations";
+                case LraMethod::GainBiasEquations:
+                    return "gain-bias-equations";
             }
             return "invalid";
         }

src/storm/solver/SolverSelectionOptions.h

@@ -9,7 +9,7 @@ namespace storm {
         ExtendEnumsWithSelectionField(MinMaxMethod, PolicyIteration, ValueIteration, LinearProgramming, Topological, RationalSearch, IntervalIteration, SoundValueIteration, TopologicalCuda, ViToPi)
         ExtendEnumsWithSelectionField(MultiplierType, Native, Gmmxx)
         ExtendEnumsWithSelectionField(GameMethod, PolicyIteration, ValueIteration)
-        ExtendEnumsWithSelectionField(LraMethod, LinearProgramming, ValueIteration)
+        ExtendEnumsWithSelectionField(LraMethod, LinearProgramming, ValueIteration, GainBiasEquations, LraDistributionEquations)
 
         ExtendEnumsWithSelectionField(LpSolverType, Gurobi, Glpk, Z3)
         ExtendEnumsWithSelectionField(EquationSolverType, Native, Gmmxx, Eigen, Elimination, Topological)

src/storm/storage/SparseMatrix.cpp

@@ -1375,30 +1375,38 @@ namespace storm {
         }
 #endif
         
+        template<typename ValueType>
+        template<typename OtherValueType, typename ResultValueType>
+        ResultValueType SparseMatrix<ValueType>::getPointwiseProductRowSum(storm::storage::SparseMatrix<OtherValueType> const& otherMatrix, index_type const& row) const {
+            typename storm::storage::SparseMatrix<ValueType>::const_iterator it1 = this->begin(row);
+            typename storm::storage::SparseMatrix<ValueType>::const_iterator ite1 = this->end(row);
+            typename storm::storage::SparseMatrix<OtherValueType>::const_iterator it2 = otherMatrix.begin(row);
+            typename storm::storage::SparseMatrix<OtherValueType>::const_iterator ite2 = otherMatrix.end(row);
+            
+            ResultValueType result = storm::utility::zero<ResultValueType>();
+            for (;it1 != ite1 && it2 != ite2; ++it1) {
+                if (it1->getColumn() < it2->getColumn()) {
+                    continue;
+                } else {
+                    // If the precondition of this method (i.e. that the given matrix is a submatrix
+                    // of the current one) was fulfilled, we know now that the two elements are in
+                    // the same column, so we can multiply and add them to the row sum vector.
+                    STORM_LOG_ASSERT(it1->getColumn() == it2->getColumn(), "The given matrix is not a submatrix of this one.");
+                    result += it2->getValue() * OtherValueType(it1->getValue());
+                    ++it2;
+                }
+            }
+            return result;
+        }
+        
         template<typename ValueType>
         template<typename OtherValueType, typename ResultValueType>
         std::vector<ResultValueType> SparseMatrix<ValueType>::getPointwiseProductRowSumVector(storm::storage::SparseMatrix<OtherValueType> const& otherMatrix) const {
-            std::vector<ResultValueType> result(rowCount, storm::utility::zero<ResultValueType>());
-            
-            // Iterate over all elements of the current matrix and either continue with the next element in case the
-            // given matrix does not have a non-zero element at this column position, or multiply the two entries and
-            // add the result to the corresponding position in the vector.
+            std::vector<ResultValueType> result;
+            result.reserve(rowCount);
             for (index_type row = 0; row < rowCount && row < otherMatrix.getRowCount(); ++row) {
-                typename storm::storage::SparseMatrix<OtherValueType>::const_iterator it2 = otherMatrix.begin(row);
-                typename storm::storage::SparseMatrix<OtherValueType>::const_iterator ite2 = otherMatrix.end(row);
-                for (const_iterator it1 = this->begin(row), ite1 = this->end(row); it1 != ite1 && it2 != ite2; ++it1) {
-                    if (it1->getColumn() < it2->getColumn()) {
-                        continue;
-                    } else {
-                        // If the precondition of this method (i.e. that the given matrix is a submatrix
-                        // of the current one) was fulfilled, we know now that the two elements are in
-                        // the same column, so we can multiply and add them to the row sum vector.
-                        result[row] += it2->getValue() * OtherValueType(it1->getValue());
-                        ++it2;
-                    }
-                }
+                result.push_back(getPointwiseProductRowSum<OtherValueType, ResultValueType>(otherMatrix, row));
             }
-            
             return result;
         }
         
@@ -2293,6 +2301,7 @@ namespace storm {
         template class SparseMatrixBuilder<double>;
         template class SparseMatrix<double>;
         template std::ostream& operator<<(std::ostream& out, SparseMatrix<double> const& matrix);
+        template double SparseMatrix<double>::getPointwiseProductRowSum(storm::storage::SparseMatrix<double> const& otherMatrix, typename SparseMatrix<double>::index_type const& row) const;
         template std::vector<double> SparseMatrix<double>::getPointwiseProductRowSumVector(storm::storage::SparseMatrix<double> const& otherMatrix) const;
         template bool SparseMatrix<double>::isSubmatrixOf(SparseMatrix<double> const& matrix) const;
         
@@ -2305,6 +2314,7 @@ namespace storm {
         template class SparseMatrixBuilder<float>;
         template class SparseMatrix<float>;
         template std::ostream& operator<<(std::ostream& out, SparseMatrix<float> const& matrix);
+        template float SparseMatrix<float>::getPointwiseProductRowSum(storm::storage::SparseMatrix<float> const& otherMatrix, typename SparseMatrix<float>::index_type const& row) const;
         template std::vector<float> SparseMatrix<float>::getPointwiseProductRowSumVector(storm::storage::SparseMatrix<float> const& otherMatrix) const;
         template bool SparseMatrix<float>::isSubmatrixOf(SparseMatrix<float> const& matrix) const;
         
@@ -2333,6 +2343,7 @@ namespace storm {
         template class SparseMatrixBuilder<ClnRationalNumber>;
         template class SparseMatrix<ClnRationalNumber>;
         template std::ostream& operator<<(std::ostream& out, SparseMatrix<ClnRationalNumber> const& matrix);
+        template storm::ClnRationalNumber SparseMatrix<storm::ClnRationalNumber>::getPointwiseProductRowSum(storm::storage::SparseMatrix<storm::ClnRationalNumber> const& otherMatrix, typename SparseMatrix<storm::ClnRationalNumber>::index_type const& row) const;
         template std::vector<storm::ClnRationalNumber> SparseMatrix<ClnRationalNumber>::getPointwiseProductRowSumVector(storm::storage::SparseMatrix<storm::ClnRationalNumber> const& otherMatrix) const;
         template bool SparseMatrix<storm::ClnRationalNumber>::isSubmatrixOf(SparseMatrix<storm::ClnRationalNumber> const& matrix) const;
 #endif
@@ -2343,6 +2354,7 @@ namespace storm {
         template class SparseMatrixBuilder<GmpRationalNumber>;
         template class SparseMatrix<GmpRationalNumber>;
         template std::ostream& operator<<(std::ostream& out, SparseMatrix<GmpRationalNumber> const& matrix);
+        template storm::GmpRationalNumber SparseMatrix<storm::GmpRationalNumber>::getPointwiseProductRowSum(storm::storage::SparseMatrix<storm::GmpRationalNumber> const& otherMatrix, typename SparseMatrix<storm::GmpRationalNumber>::index_type const& row) const;
         template std::vector<storm::GmpRationalNumber> SparseMatrix<GmpRationalNumber>::getPointwiseProductRowSumVector(storm::storage::SparseMatrix<storm::GmpRationalNumber> const& otherMatrix) const;
         template bool SparseMatrix<storm::GmpRationalNumber>::isSubmatrixOf(SparseMatrix<storm::GmpRationalNumber> const& matrix) const;
 #endif
@@ -2353,6 +2365,10 @@ namespace storm {
         template class SparseMatrixBuilder<RationalFunction>;
         template class SparseMatrix<RationalFunction>;
         template std::ostream& operator<<(std::ostream& out, SparseMatrix<RationalFunction> const& matrix);
+        template storm::RationalFunction SparseMatrix<storm::RationalFunction>::getPointwiseProductRowSum(storm::storage::SparseMatrix<storm::RationalFunction> const& otherMatrix, typename SparseMatrix<storm::RationalFunction>::index_type const& row) const;
+        template storm::RationalFunction SparseMatrix<double>::getPointwiseProductRowSum(storm::storage::SparseMatrix<storm::RationalFunction> const& otherMatrix, typename SparseMatrix<storm::RationalFunction>::index_type const& row) const;
+        template storm::RationalFunction SparseMatrix<float>::getPointwiseProductRowSum(storm::storage::SparseMatrix<storm::RationalFunction> const& otherMatrix, typename SparseMatrix<storm::RationalFunction>::index_type const& row) const;
+        template storm::RationalFunction SparseMatrix<int>::getPointwiseProductRowSum(storm::storage::SparseMatrix<storm::RationalFunction> const& otherMatrix, typename SparseMatrix<storm::RationalFunction>::index_type const& row) const;
         template std::vector<storm::RationalFunction> SparseMatrix<RationalFunction>::getPointwiseProductRowSumVector(storm::storage::SparseMatrix<storm::RationalFunction> const& otherMatrix) const;
         template std::vector<storm::RationalFunction> SparseMatrix<double>::getPointwiseProductRowSumVector(storm::storage::SparseMatrix<storm::RationalFunction> const& otherMatrix) const;
         template std::vector<storm::RationalFunction> SparseMatrix<float>::getPointwiseProductRowSumVector(storm::storage::SparseMatrix<storm::RationalFunction> const& otherMatrix) const;

src/storm/storage/SparseMatrix.h

@@ -819,6 +819,19 @@ namespace storm {
              */
             std::pair<storm::storage::SparseMatrix<value_type>, std::vector<value_type>> getJacobiDecomposition() const;
             
+            /*!
+             * Performs a pointwise multiplication of the entries in the given row of this matrix and the entries of
+             * the given row of the other matrix and returns the sum.
+             *
+             * @param otherMatrix A reference to the matrix with which to perform the pointwise multiplication. This
+             * matrix must be a submatrix of the current matrix in the sense that it may not have entries at indices
+             * where there is no entry in the current matrix.
+             * @return the sum of the product of the entries in the given row.
+             */
+            template<typename OtherValueType, typename ResultValueType = OtherValueType>
+            ResultValueType getPointwiseProductRowSum(storm::storage::SparseMatrix<OtherValueType> const& otherMatrix, index_type const& row) const;
+            
+            
             /*!
              * Performs a pointwise matrix multiplication of the matrix with the given matrix and returns a vector
              * containing the sum of the entries in each row of the resulting matrix.

src/test/storm-pars/analysis/MonotonicityCheckerTest.cpp

@@ -149,7 +149,7 @@ std::vector<storm::storage::ParameterRegion<storm::RationalFunction>> regions =
     ASSERT_EQ(dtmc->getNumberOfTransitions(), 195ull);
 
     storm::analysis::MonotonicityChecker<storm::RationalFunction> monotonicityChecker = storm::analysis::MonotonicityChecker<storm::RationalFunction>(dtmc, formulas, regions, true);
-    auto result = monotonicityChecker.checkMonotonicity();
+    auto result = monotonicityChecker.checkMonotonicity(std::cout);
     EXPECT_EQ(1ul, result.size());
     EXPECT_EQ(2ul, result.begin()->second.size());
     auto monotone = result.begin()->second.begin();
@@ -197,7 +197,7 @@ TEST(MonotonicityCheckerTest, Brp_with_bisimulation_samples) {
     ASSERT_EQ(dtmc->getNumberOfTransitions(), 195ull);
 
     auto monotonicityChecker = storm::analysis::MonotonicityChecker<storm::RationalFunction>(dtmc, formulas, regions, true, 50);
-    auto result = monotonicityChecker.checkMonotonicity();
+    auto result = monotonicityChecker.checkMonotonicity(std::cout);
     EXPECT_EQ(1ul, result.size());
     EXPECT_EQ(2ul, result.begin()->second.size());
     auto monotone = result.begin()->second.begin();

src/test/storm/modelchecker/csl/CtmcCslModelCheckerTest.cpp

@@ -302,7 +302,6 @@ namespace {
         formulasString += "; P=? [ !\"minimum\" U>=1 \"minimum\"]";
         formulasString += "; P=? [ \"minimum\" U>=1 !\"minimum\"]";
         formulasString += "; R=? [C<=100]";
-        formulasString += "; LRA=? [\"minimum\"]";
         
         auto modelFormulas = this->buildModelFormulas(STORM_TEST_RESOURCES_DIR "/ctmc/cluster2.sm", formulasString);
         auto model = std::move(modelFormulas.first);
@@ -333,9 +332,6 @@ namespace {
         
         result = checker->check(this->env(), tasks[6]);
         EXPECT_NEAR(this->parseNumber("0.8602815057967503"), this->getQuantitativeResultAtInitialState(model, result), this->precision());
-        
-        result = checker->check(this->env(), tasks[7]);
-        EXPECT_NEAR(this->parseNumber("0.99999766034263426"), this->getQuantitativeResultAtInitialState(model, result), this->precision());
     
     }
     
@@ -345,7 +341,6 @@ namespace {
         formulasString += "; P=? [ !\"down\" U<=10000 \"fail_io\"]";
         formulasString += "; P=? [ !\"down\" U<=10000 \"fail_sensors\"]";
         formulasString += "; R{\"up\"}=? [C<=10000]";
-        formulasString += "; LRA=? [\"fail_sensors\"]";
         
         auto modelFormulas = this->buildModelFormulas(STORM_TEST_RESOURCES_DIR "/ctmc/embedded2.sm", formulasString);
         auto model = std::move(modelFormulas.first);
@@ -371,29 +366,6 @@ namespace {
         result = checker->check(this->env(), tasks[4]);
         EXPECT_NEAR(this->parseNumber("2.7745274082080154"), this->getQuantitativeResultAtInitialState(model, result), this->precision());
         
-        result = checker->check(this->env(), tasks[5]);
-        EXPECT_NEAR(this->parseNumber("0.93458866427696596"), this->getQuantitativeResultAtInitialState(model, result), this->precision());
-    }
-    
-    TYPED_TEST(CtmcCslModelCheckerTest, Polling) {
-        std::string formulasString = "P=?[ F<=10 \"target\"]";
-        formulasString += "; LRA=?[\"target\"]";
-        
-        auto modelFormulas = this->buildModelFormulas(STORM_TEST_RESOURCES_DIR "/ctmc/polling2.sm", formulasString);
-        auto model = std::move(modelFormulas.first);
-        auto tasks = this->getTasks(modelFormulas.second);
-        EXPECT_EQ(12ul, model->getNumberOfStates());
-        EXPECT_EQ(22ul, model->getNumberOfTransitions());
-        ASSERT_EQ(model->getType(), storm::models::ModelType::Ctmc);
-        auto checker = this->createModelChecker(model);
-        std::unique_ptr<storm::modelchecker::CheckResult> result;
-        
-        result = checker->check(this->env(), tasks[0]);
-        EXPECT_NEAR(this->parseNumber("1"), this->getQuantitativeResultAtInitialState(model, result), this->precision());
-        
-        result = checker->check(this->env(), tasks[1]);
-        EXPECT_NEAR(this->parseNumber("0.20079750055570736"), this->getQuantitativeResultAtInitialState(model, result), this->precision());
-        
     }
     
     TYPED_TEST(CtmcCslModelCheckerTest, Tandem) {
@@ -403,7 +375,6 @@ namespace {
         formulasString += "; R=? [I=10]";
         formulasString += "; R=? [C<=10]";
         formulasString += "; R=? [F \"first_queue_full\"&\"second_queue_full\"]";
-        formulasString += "; LRA=? [\"first_queue_full\"]";
         
         auto modelFormulas = this->buildModelFormulas(STORM_TEST_RESOURCES_DIR "/ctmc/tandem5.sm", formulasString);
         auto model = std::move(modelFormulas.first);
@@ -431,9 +402,6 @@ namespace {
         
         result = checker->check(this->env(), tasks[5]);
         EXPECT_NEAR(this->parseNumber("262.85103824276308"), this->getQuantitativeResultAtInitialState(model, result), this->precision());
-        
-        result = checker->check(this->env(), tasks[6]);
-        EXPECT_NEAR(this->parseNumber("0.9100373532"), this->getQuantitativeResultAtInitialState(model, result), this->precision());
     }
     
     TYPED_TEST(CtmcCslModelCheckerTest, simple2) {

src/test/storm/modelchecker/csl/LraCtmcCslModelCheckerTest.cpp

@@ -0,0 +1,416 @@
+#include "test/storm_gtest.h"
+#include "test/storm_gtest.h"
+#include "storm-config.h"
+
+#include "storm/api/builder.h"
+#include "storm-parsers/api/model_descriptions.h"
+#include "storm/api/properties.h"
+#include "storm-conv/api/storm-conv.h"
+#include "storm-parsers/api/properties.h"
+#include "storm-parsers/parser/FormulaParser.h"
+#include "storm/logic/Formulas.h"
+#include "storm/models/sparse/StandardRewardModel.h"
+#include "storm/models/symbolic/StandardRewardModel.h"
+#include "storm/models/sparse/Ctmc.h"
+#include "storm/models/symbolic/Ctmc.h"
+#include "storm/modelchecker/csl/SparseCtmcCslModelChecker.h"
+#include "storm/modelchecker/csl/HybridCtmcCslModelChecker.h"
+#include "storm/modelchecker/results/QuantitativeCheckResult.h"
+#include "storm/modelchecker/results/ExplicitQualitativeCheckResult.h"
+#include "storm/modelchecker/results/SymbolicQualitativeCheckResult.h"
+#include "storm/modelchecker/results/QualitativeCheckResult.h"
+#include "storm-parsers/parser/PrismParser.h"
+#include "storm/storage/expressions/ExpressionManager.h"
+#include "storm/settings/modules/CoreSettings.h"
+#include "storm/environment/solver/NativeSolverEnvironment.h"
+#include "storm/environment/solver/GmmxxSolverEnvironment.h"
+#include "storm/environment/solver/EigenSolverEnvironment.h"
+#include "storm/environment/solver/LongRunAverageSolverEnvironment.h"
+
+namespace {
+    
+    enum class CtmcEngine {PrismSparse, JaniSparse, JaniHybrid};
+    
+    class GBSparseGmmxxGmresIluEnvironment {
+    public:
+        static const storm::dd::DdType ddType = storm::dd::DdType::Sylvan; // unused for sparse models
+        static const CtmcEngine engine = CtmcEngine::PrismSparse;
+        static const bool isExact = false;
+        typedef double ValueType;
+        typedef storm::models::sparse::Ctmc<ValueType> ModelType;
+        static storm::Environment createEnvironment() {
+            storm::Environment env;
+            env.solver().setLinearEquationSolverType(storm::solver::EquationSolverType::Gmmxx);
+            env.solver().gmmxx().setMethod(storm::solver::GmmxxLinearEquationSolverMethod::Gmres);
+            env.solver().gmmxx().setPreconditioner(storm::solver::GmmxxLinearEquationSolverPreconditioner::Ilu);
+            env.solver().lra().setDetLraMethod(storm::solver::LraMethod::GainBiasEquations);
+            return env;
+        }
+    };
+    
+    class GBJaniSparseGmmxxGmresIluEnvironment {
+    public:
+        static const storm::dd::DdType ddType = storm::dd::DdType::Sylvan; // unused for sparse models
+        static const CtmcEngine engine = CtmcEngine::JaniSparse;
+        static const bool isExact = false;
+        typedef double ValueType;
+        typedef storm::models::sparse::Ctmc<ValueType> ModelType;
+        static storm::Environment createEnvironment() {
+            storm::Environment env;
+            env.solver().setLinearEquationSolverType(storm::solver::EquationSolverType::Gmmxx);
+            env.solver().gmmxx().setMethod(storm::solver::GmmxxLinearEquationSolverMethod::Gmres);
+            env.solver().gmmxx().setPreconditioner(storm::solver::GmmxxLinearEquationSolverPreconditioner::Ilu);
+            env.solver().lra().setDetLraMethod(storm::solver::LraMethod::GainBiasEquations);
+            return env;
+        }
+    };
+    
+    class GBJaniHybridCuddGmmxxGmresEnvironment {
+    public:
+        static const storm::dd::DdType ddType = storm::dd::DdType::CUDD;
+        static const CtmcEngine engine = CtmcEngine::JaniHybrid;
+        static const bool isExact = false;
+        typedef double ValueType;
+        typedef storm::models::symbolic::Ctmc<ddType, ValueType> ModelType;
+        static storm::Environment createEnvironment() {
+            storm::Environment env;
+            env.solver().setLinearEquationSolverType(storm::solver::EquationSolverType::Gmmxx);
+            env.solver().gmmxx().setMethod(storm::solver::GmmxxLinearEquationSolverMethod::Gmres);
+            env.solver().lra().setDetLraMethod(storm::solver::LraMethod::GainBiasEquations);
+            env.solver().lra().setPrecision(storm::utility::convertNumber<storm::RationalNumber>(1e-8)); // Need to increase precision because eq sys yields incorrect results
+            return env;
+        }
+    };
+    
+    class GBJaniHybridSylvanGmmxxGmresEnvironment {
+    public:
+        static const storm::dd::DdType ddType = storm::dd::DdType::Sylvan;
+        static const CtmcEngine engine = CtmcEngine::JaniHybrid;
+        static const bool isExact = false;
+        typedef double ValueType;
+        typedef storm::models::symbolic::Ctmc<ddType, ValueType> ModelType;
+        static storm::Environment createEnvironment() {
+            storm::Environment env;
+            env.solver().setLinearEquationSolverType(storm::solver::EquationSolverType::Gmmxx);
+            env.solver().gmmxx().setMethod(storm::solver::GmmxxLinearEquationSolverMethod::Gmres);
+            env.solver().lra().setDetLraMethod(storm::solver::LraMethod::GainBiasEquations);
+            env.solver().lra().setPrecision(storm::utility::convertNumber<storm::RationalNumber>(1e-8)); // Need to increase precision because eq sys yields incorrect results
+            return env;
+        }
+    };
+    
+    class GBSparseEigenDGmresEnvironment {
+    public:
+        static const storm::dd::DdType ddType = storm::dd::DdType::Sylvan; // unused for sparse models
+        static const CtmcEngine engine = CtmcEngine::PrismSparse;
+        static const bool isExact = false;
+        typedef double ValueType;
+        typedef storm::models::sparse::Ctmc<ValueType> ModelType;
+        static storm::Environment createEnvironment() {
+            storm::Environment env;
+            env.solver().setLinearEquationSolverType(storm::solver::EquationSolverType::Eigen);
+            env.solver().eigen().setMethod(storm::solver::EigenLinearEquationSolverMethod::DGmres);
+            env.solver().eigen().setPreconditioner(storm::solver::EigenLinearEquationSolverPreconditioner::Ilu);
+            env.solver().lra().setDetLraMethod(storm::solver::LraMethod::GainBiasEquations);
+            return env;
+        }
+    };
+    
+    class GBSparseEigenRationalLUEnvironment {
+    public:
+        static const storm::dd::DdType ddType = storm::dd::DdType::Sylvan; // unused for sparse models
+        static const CtmcEngine engine = CtmcEngine::PrismSparse;
+        static const bool isExact = false;
+        typedef double ValueType;
+        typedef storm::models::sparse::Ctmc<ValueType> ModelType;
+        static storm::Environment createEnvironment() {
+            storm::Environment env;
+            env.solver().setLinearEquationSolverType(storm::solver::EquationSolverType::Eigen);
+            env.solver().eigen().setMethod(storm::solver::EigenLinearEquationSolverMethod::SparseLU);
+            env.solver().lra().setDetLraMethod(storm::solver::LraMethod::GainBiasEquations);
+            return env;
+        }
+    };
+    
+    class GBSparseNativeSorEnvironment {
+    public:
+        static const storm::dd::DdType ddType = storm::dd::DdType::Sylvan; // unused for sparse models
+        static const CtmcEngine engine = CtmcEngine::PrismSparse;
+        static const bool isExact = false;
+        typedef double ValueType;
+        typedef storm::models::sparse::Ctmc<ValueType> ModelType;
+        static storm::Environment createEnvironment() {
+            storm::Environment env;
+            env.solver().setLinearEquationSolverType(storm::solver::EquationSolverType::Native);
+            env.solver().native().setMethod(storm::solver::NativeLinearEquationSolverMethod::SOR);
+            env.solver().native().setSorOmega(storm::utility::convertNumber<storm::RationalNumber>(0.7)); // LRA computation fails for 0.9
+            env.solver().lra().setDetLraMethod(storm::solver::LraMethod::GainBiasEquations);
+            env.solver().lra().setPrecision(storm::utility::convertNumber<storm::RationalNumber>(1e-9));
+            return env;
+        }
+    };
+
+    class DistrSparseGmmxxGmresIluEnvironment {
+    public:
+        static const storm::dd::DdType ddType = storm::dd::DdType::Sylvan; // unused for sparse models
+        static const CtmcEngine engine = CtmcEngine::PrismSparse;
+        static const bool isExact = false;
+        typedef double ValueType;
+        typedef storm::models::sparse::Ctmc<ValueType> ModelType;
+        static storm::Environment createEnvironment() {
+            storm::Environment env;
+            env.solver().setLinearEquationSolverType(storm::solver::EquationSolverType::Gmmxx);
+            env.solver().gmmxx().setMethod(storm::solver::GmmxxLinearEquationSolverMethod::Gmres);
+            env.solver().gmmxx().setPreconditioner(storm::solver::GmmxxLinearEquationSolverPreconditioner::Ilu);
+            env.solver().lra().setDetLraMethod(storm::solver::LraMethod::LraDistributionEquations);
+            env.solver().lra().setPrecision(storm::utility::convertNumber<storm::RationalNumber>(1e-8)); // Need to increase precision because eq sys yields incorrect results
+            return env;
+        }
+    };
+    
+    class DistrSparseEigenDoubleLUEnvironment {
+    public:
+        static const storm::dd::DdType ddType = storm::dd::DdType::Sylvan; // unused for sparse models
+        static const CtmcEngine engine = CtmcEngine::PrismSparse;
+        static const bool isExact = false;
+        typedef double ValueType;
+        typedef storm::models::sparse::Ctmc<ValueType> ModelType;
+        static storm::Environment createEnvironment() {
+            storm::Environment env;
+            env.solver().setLinearEquationSolverType(storm::solver::EquationSolverType::Eigen);
+            env.solver().eigen().setMethod(storm::solver::EigenLinearEquationSolverMethod::DGmres);
+            env.solver().eigen().setPreconditioner(storm::solver::EigenLinearEquationSolverPreconditioner::Ilu);
+            env.solver().eigen().setPrecision(storm::utility::convertNumber<storm::RationalNumber>(1e-8));
+            env.solver().lra().setDetLraMethod(storm::solver::LraMethod::LraDistributionEquations);
+            return env;
+        }
+    };
+    
+    class ValueIterationSparseEnvironment {
+    public:
+        static const storm::dd::DdType ddType = storm::dd::DdType::Sylvan; // unused for sparse models
+        static const CtmcEngine engine = CtmcEngine::PrismSparse;
+        static const bool isExact = false;
+        typedef double ValueType;
+        typedef storm::models::sparse::Ctmc<ValueType> ModelType;
+        static storm::Environment createEnvironment() {
+            storm::Environment env;
+            env.solver().lra().setDetLraMethod(storm::solver::LraMethod::ValueIteration);
+            return env;
+        }
+    };
+    
+    template<typename TestType>
+    class LraCtmcCslModelCheckerTest : public ::testing::Test {
+    public:
+        typedef typename TestType::ValueType ValueType;
+        typedef typename storm::models::sparse::Ctmc<ValueType> SparseModelType;
+        typedef typename storm::models::symbolic::Ctmc<TestType::ddType, ValueType> SymbolicModelType;
+        
+        LraCtmcCslModelCheckerTest() : _environment(TestType::createEnvironment()) {}
+        storm::Environment const& env() const { return _environment; }
+        ValueType parseNumber(std::string const& input) const { return storm::utility::convertNumber<ValueType>(input);}
+        ValueType precision() const { return TestType::isExact ? parseNumber("0") : parseNumber("1e-6");}
+        bool isSparseModel() const { return std::is_same<typename TestType::ModelType, SparseModelType>::value; }
+        bool isSymbolicModel() const { return std::is_same<typename TestType::ModelType, SymbolicModelType>::value; }
+        CtmcEngine getEngine() const { return TestType::engine; }
+        
+        template <typename MT = typename TestType::ModelType>
+        typename std::enable_if<std::is_same<MT, SparseModelType>::value, std::pair<std::shared_ptr<MT>, std::vector<std::shared_ptr<storm::logic::Formula const>>>>::type
+        buildModelFormulas(std::string const& pathToPrismFile, std::string const& formulasAsString, std::string const& constantDefinitionString = "") const {
+            std::pair<std::shared_ptr<MT>, std::vector<std::shared_ptr<storm::logic::Formula const>>> result;
+            storm::prism::Program program = storm::api::parseProgram(pathToPrismFile, true);
+            program = storm::utility::prism::preprocess(program, constantDefinitionString);
+            if (TestType::engine == CtmcEngine::PrismSparse) {
+                result.second = storm::api::extractFormulasFromProperties(storm::api::parsePropertiesForPrismProgram(formulasAsString, program));
+                result.first = storm::api::buildSparseModel<ValueType>(program, result.second)->template as<MT>();
+            } else if (TestType::engine == CtmcEngine::JaniSparse) {
+                auto janiData = storm::api::convertPrismToJani(program, storm::api::parsePropertiesForPrismProgram(formulasAsString, program));
+                janiData.first.substituteFunctions();
+                result.second = storm::api::extractFormulasFromProperties(janiData.second);
+                result.first = storm::api::buildSparseModel<ValueType>(janiData.first, result.second)->template as<MT>();
+            }
+            return result;
+        }
+        
+        template <typename MT = typename TestType::ModelType>
+        typename std::enable_if<std::is_same<MT, SymbolicModelType>::value, std::pair<std::shared_ptr<MT>, std::vector<std::shared_ptr<storm::logic::Formula const>>>>::type
+        buildModelFormulas(std::string const& pathToPrismFile, std::string const& formulasAsString, std::string const& constantDefinitionString = "") const {
+            std::pair<std::shared_ptr<MT>, std::vector<std::shared_ptr<storm::logic::Formula const>>> result;
+            storm::prism::Program program = storm::api::parseProgram(pathToPrismFile, true);
+            program = storm::utility::prism::preprocess(program, constantDefinitionString);
+            if (TestType::engine == CtmcEngine::JaniHybrid) {
+                auto janiData = storm::api::convertPrismToJani(program, storm::api::parsePropertiesForPrismProgram(formulasAsString, program));
+                janiData.first.substituteFunctions();
+                result.second = storm::api::extractFormulasFromProperties(janiData.second);
+                result.first = storm::api::buildSymbolicModel<TestType::ddType, ValueType>(janiData.first, result.second)->template as<MT>();
+            }
+            return result;
+        }
+        
+        std::vector<storm::modelchecker::CheckTask<storm::logic::Formula, ValueType>> getTasks(std::vector<std::shared_ptr<storm::logic::Formula const>> const& formulas) const {
+            std::vector<storm::modelchecker::CheckTask<storm::logic::Formula, ValueType>> result;
+            for (auto const& f : formulas) {
+                result.emplace_back(*f);
+            }
+            return result;
+        }
+        
+        template <typename MT = typename TestType::ModelType>
+        typename std::enable_if<std::is_same<MT, SparseModelType>::value, std::shared_ptr<storm::modelchecker::AbstractModelChecker<MT>>>::type
+        createModelChecker(std::shared_ptr<MT> const& model) const {
+            if (TestType::engine == CtmcEngine::PrismSparse || TestType::engine == CtmcEngine::JaniSparse) {
+                return std::make_shared<storm::modelchecker::SparseCtmcCslModelChecker<SparseModelType>>(*model);
+            }
+            return nullptr;
+        }
+        
+        template <typename MT = typename TestType::ModelType>
+        typename std::enable_if<std::is_same<MT, SymbolicModelType>::value, std::shared_ptr<storm::modelchecker::AbstractModelChecker<MT>>>::type
+        createModelChecker(std::shared_ptr<MT> const& model) const {
+            if (TestType::engine == CtmcEngine::JaniHybrid) {
+                return std::make_shared<storm::modelchecker::HybridCtmcCslModelChecker<SymbolicModelType>>(*model);
+            }
+            return nullptr;
+        }
+        
+        bool getQualitativeResultAtInitialState(std::shared_ptr<storm::models::Model<ValueType>> const& model, std::unique_ptr<storm::modelchecker::CheckResult>& result) {
+            auto filter = getInitialStateFilter(model);
+            result->filter(*filter);
+            return result->asQualitativeCheckResult().forallTrue();
+        }
+        
+        ValueType getQuantitativeResultAtInitialState(std::shared_ptr<storm::models::Model<ValueType>> const& model, std::unique_ptr<storm::modelchecker::CheckResult>& result) {
+            auto filter = getInitialStateFilter(model);
+            result->filter(*filter);
+            return result->asQuantitativeCheckResult<ValueType>().getMin();
+        }
+    
+    private:
+        storm::Environment _environment;
+        
+        std::unique_ptr<storm::modelchecker::QualitativeCheckResult> getInitialStateFilter(std::shared_ptr<storm::models::Model<ValueType>> const& model) const {
+            if (isSparseModel()) {
+                return std::make_unique<storm::modelchecker::ExplicitQualitativeCheckResult>(model->template as<SparseModelType>()->getInitialStates());
+            } else {
+                return std::make_unique<storm::modelchecker::SymbolicQualitativeCheckResult<TestType::ddType>>(model->template as<SymbolicModelType>()->getReachableStates(), model->template as<SymbolicModelType>()->getInitialStates());
+            }
+        }
+    };
+  
+    typedef ::testing::Types<
+            GBSparseGmmxxGmresIluEnvironment,
+            GBJaniSparseGmmxxGmresIluEnvironment,
+            GBJaniHybridCuddGmmxxGmresEnvironment,
+            GBJaniHybridSylvanGmmxxGmresEnvironment,
+            GBSparseEigenDGmresEnvironment,
+            GBSparseEigenRationalLUEnvironment,
+            GBSparseNativeSorEnvironment,
+            DistrSparseGmmxxGmresIluEnvironment,
+            DistrSparseEigenDoubleLUEnvironment,
+            ValueIterationSparseEnvironment
+        > TestingTypes;
+    
+    TYPED_TEST_SUITE(LraCtmcCslModelCheckerTest, TestingTypes,);
+
+    TYPED_TEST(LraCtmcCslModelCheckerTest, Cluster) {
+        std::string formulasString = "LRA=? [\"minimum\"]";
+        
+        auto modelFormulas = this->buildModelFormulas(STORM_TEST_RESOURCES_DIR "/ctmc/cluster2.sm", formulasString);
+        auto model = std::move(modelFormulas.first);
+        auto tasks = this->getTasks(modelFormulas.second);
+        EXPECT_EQ(276ul, model->getNumberOfStates());
+        EXPECT_EQ(1120ul, model->getNumberOfTransitions());
+        ASSERT_EQ(model->getType(), storm::models::ModelType::Ctmc);
+        auto checker = this->createModelChecker(model);
+        std::unique_ptr<storm::modelchecker::CheckResult> result;
+        
+        result = checker->check(this->env(), tasks[0]);
+        EXPECT_NEAR(this->parseNumber("0.99999766034263426"), this->getQuantitativeResultAtInitialState(model, result), this->precision());
+    
+    }
+    
+    TYPED_TEST(LraCtmcCslModelCheckerTest, Embedded) {
+        std::string formulasString = "LRA=? [\"fail_sensors\"]";
+        
+        auto modelFormulas = this->buildModelFormulas(STORM_TEST_RESOURCES_DIR "/ctmc/embedded2.sm", formulasString);
+        auto model = std::move(modelFormulas.first);
+        auto tasks = this->getTasks(modelFormulas.second);
+        EXPECT_EQ(3478ul, model->getNumberOfStates());
+        EXPECT_EQ(14639ul, model->getNumberOfTransitions());
+        ASSERT_EQ(model->getType(), storm::models::ModelType::Ctmc);
+        auto checker = this->createModelChecker(model);
+        std::unique_ptr<storm::modelchecker::CheckResult> result;
+        
+        result = checker->check(this->env(), tasks[0]);
+        EXPECT_NEAR(this->parseNumber("0.93458866427696596"), this->getQuantitativeResultAtInitialState(model, result), this->precision());
+    }
+    
+    TYPED_TEST(LraCtmcCslModelCheckerTest, Polling) {
+        std::string formulasString = "LRA=?[\"target\"]";
+        
+        auto modelFormulas = this->buildModelFormulas(STORM_TEST_RESOURCES_DIR "/ctmc/polling2.sm", formulasString);
+        auto model = std::move(modelFormulas.first);
+        auto tasks = this->getTasks(modelFormulas.second);
+        EXPECT_EQ(12ul, model->getNumberOfStates());
+        EXPECT_EQ(22ul, model->getNumberOfTransitions());
+        ASSERT_EQ(model->getType(), storm::models::ModelType::Ctmc);
+        auto checker = this->createModelChecker(model);
+        std::unique_ptr<storm::modelchecker::CheckResult> result;
+        
+        result = checker->check(this->env(), tasks[0]);
+        EXPECT_NEAR(this->parseNumber("0.20079750055570736"), this->getQuantitativeResultAtInitialState(model, result), this->precision());
+        
+    }
+    
+    TYPED_TEST(LraCtmcCslModelCheckerTest, Tandem) {
+        std::string formulasString = "LRA=? [\"first_queue_full\"]";
+        
+        auto modelFormulas = this->buildModelFormulas(STORM_TEST_RESOURCES_DIR "/ctmc/tandem5.sm", formulasString);
+        auto model = std::move(modelFormulas.first);
+        auto tasks = this->getTasks(modelFormulas.second);
+        EXPECT_EQ(66ul, model->getNumberOfStates());
+        EXPECT_EQ(189ul, model->getNumberOfTransitions());
+        ASSERT_EQ(model->getType(), storm::models::ModelType::Ctmc);
+        auto checker = this->createModelChecker(model);
+        std::unique_ptr<storm::modelchecker::CheckResult> result;
+        
+        result = checker->check(this->env(), tasks[0]);
+        EXPECT_NEAR(this->parseNumber("0.9100373532"), this->getQuantitativeResultAtInitialState(model, result), this->precision());
+    }
+    
+    
+    TYPED_TEST(LraCtmcCslModelCheckerTest, Rewards) {
+        std::string formulasString = "R=? [ LRA ]";
+        
+        auto modelFormulas = this->buildModelFormulas(STORM_TEST_RESOURCES_DIR "/ctmc/lrarewards.sm", formulasString);
+        auto model = std::move(modelFormulas.first);
+        auto tasks = this->getTasks(modelFormulas.second);
+        EXPECT_EQ(4ul, model->getNumberOfStates());
+        EXPECT_EQ(6ul, model->getNumberOfTransitions());
+        ASSERT_EQ(model->getType(), storm::models::ModelType::Ctmc);
+        auto checker = this->createModelChecker(model);
+        std::unique_ptr<storm::modelchecker::CheckResult> result;
+        
+        result = checker->check(this->env(), tasks[0]);
+        EXPECT_NEAR(this->parseNumber("11/15"), this->getQuantitativeResultAtInitialState(model, result), this->precision());
+    }
+    
+    
+    TYPED_TEST(LraCtmcCslModelCheckerTest, kanban) {
+        std::string formulasString = "R{\"throughput\"}=? [ LRA ]";
+        
+        auto modelFormulas = this->buildModelFormulas(STORM_TEST_RESOURCES_DIR "/ctmc/kanban.prism", formulasString);
+        auto model = std::move(modelFormulas.first);
+        auto tasks = this->getTasks(modelFormulas.second);
+        EXPECT_EQ(160ul, model->getNumberOfStates());
+        EXPECT_EQ(616ul, model->getNumberOfTransitions());
+        ASSERT_EQ(model->getType(), storm::models::ModelType::Ctmc);
+        auto checker = this->createModelChecker(model);
+        std::unique_ptr<storm::modelchecker::CheckResult> result;
+        
+        result = checker->check(this->env(), tasks[0]);
+        EXPECT_NEAR(this->parseNumber("113237255213395163953677015242972426399989689654967642609491830216061334090202313396984106738516704120069048184391587092670711590526535239899047608853509681914074220789038015289373871985431257486278/1223067474012215838745994023624181143448435715858923491568712969277184634645504346456117443333632535902774869182827010789201972713368729205674432492059242349591780604188152950845769793378621446766887"), this->getQuantitativeResultAtInitialState(model, result), this->precision());
+    }
+    
+}

src/test/storm/modelchecker/prctl/dtmc/LraDtmcPrctlModelCheckerTest.cpp

@@ -16,80 +16,147 @@
 #include "storm/environment/solver/NativeSolverEnvironment.h"
 #include "storm/environment/solver/EigenSolverEnvironment.h"
 #include "storm/environment/solver/GmmxxSolverEnvironment.h"
+#include "storm/environment/solver/LongRunAverageSolverEnvironment.h"
 
 #include "storm-parsers/parser/AutoParser.h"
 #include "storm/builder/ExplicitModelBuilder.h"
 
 namespace {
     
-    class GmmxxDoubleGmresEnvironment {
+    class GBGmmxxDoubleGmresEnvironment {
     public:
         typedef double ValueType;
         static const bool isExact = false;
         static storm::Environment createEnvironment() {
             storm::Environment env;
+            env.solver().lra().setDetLraMethod(storm::solver::LraMethod::GainBiasEquations);
             env.solver().setLinearEquationSolverType(storm::solver::EquationSolverType::Gmmxx);
             env.solver().gmmxx().setMethod(storm::solver::GmmxxLinearEquationSolverMethod::Gmres);
             env.solver().gmmxx().setPreconditioner(storm::solver::GmmxxLinearEquationSolverPreconditioner::Ilu);
             env.solver().gmmxx().setPrecision(storm::utility::convertNumber<storm::RationalNumber>(1e-8));
+            env.solver().lra().setPrecision(storm::utility::convertNumber<storm::RationalNumber>(1e-8));
             return env;
         }
     };
     
-    class EigenDoubleDGmresEnvironment {
+    class GBEigenDoubleDGmresEnvironment {
     public:
         typedef double ValueType;
         static const bool isExact = false;
         static storm::Environment createEnvironment() {
             storm::Environment env;
+            env.solver().lra().setDetLraMethod(storm::solver::LraMethod::GainBiasEquations);
             env.solver().setLinearEquationSolverType(storm::solver::EquationSolverType::Eigen);
             env.solver().eigen().setMethod(storm::solver::EigenLinearEquationSolverMethod::DGmres);
             env.solver().eigen().setPreconditioner(storm::solver::EigenLinearEquationSolverPreconditioner::Ilu);
+            env.solver().lra().setPrecision(storm::utility::convertNumber<storm::RationalNumber>(1e-8));
             env.solver().eigen().setPrecision(storm::utility::convertNumber<storm::RationalNumber>(1e-8));
             return env;
         }
     };
     
-    class EigenRationalLUEnvironment {
+    class GBEigenRationalLUEnvironment {
     public:
-        typedef double ValueType;
+        typedef storm::RationalNumber ValueType;
         static const bool isExact = true;
         static storm::Environment createEnvironment() {
             storm::Environment env;
+            env.solver().lra().setDetLraMethod(storm::solver::LraMethod::GainBiasEquations);
             env.solver().setLinearEquationSolverType(storm::solver::EquationSolverType::Eigen);
             env.solver().eigen().setMethod(storm::solver::EigenLinearEquationSolverMethod::SparseLU);
             return env;
         }
     };
     
-    class NativeSorEnvironment {
+    class GBNativeSorEnvironment {
     public:
         typedef double ValueType;
         static const bool isExact = false;
         static storm::Environment createEnvironment() {
             storm::Environment env;
+            env.solver().lra().setDetLraMethod(storm::solver::LraMethod::GainBiasEquations);
             env.solver().setLinearEquationSolverType(storm::solver::EquationSolverType::Native);
             env.solver().native().setMethod(storm::solver::NativeLinearEquationSolverMethod::SOR);
-            env.solver().native().setSorOmega(storm::utility::convertNumber<storm::RationalNumber>(0.9));
+            env.solver().native().setSorOmega(storm::utility::convertNumber<storm::RationalNumber>(0.8)); // A test fails if this is set to 0.9...
             env.solver().native().setPrecision(storm::utility::convertNumber<storm::RationalNumber>(1e-8));
+            env.solver().lra().setPrecision(storm::utility::convertNumber<storm::RationalNumber>(1e-8));
+            return env;
+        }
+    };
+    
+    class GBNativeWalkerChaeEnvironment {
+    public:
+        typedef double ValueType;
+        static const bool isExact = false;
+        static storm::Environment createEnvironment() {
+            storm::Environment env;
+            env.solver().lra().setDetLraMethod(storm::solver::LraMethod::GainBiasEquations);
+            env.solver().setLinearEquationSolverType(storm::solver::EquationSolverType::Native);
+            env.solver().native().setMethod(storm::solver::NativeLinearEquationSolverMethod::WalkerChae);
+            env.solver().native().setPrecision(storm::utility::convertNumber<storm::RationalNumber>(1e-8));
+            env.solver().lra().setPrecision(storm::utility::convertNumber<storm::RationalNumber>(1e-8));
+            env.solver().native().setMaximalNumberOfIterations(50000);
+            return env;
+        }
+    };
+    
+    class DistrGmmxxDoubleGmresEnvironment {
+    public:
+        typedef double ValueType;
+        static const bool isExact = false;
+        static storm::Environment createEnvironment() {
+            storm::Environment env;
+            env.solver().lra().setDetLraMethod(storm::solver::LraMethod::LraDistributionEquations);
+            env.solver().setLinearEquationSolverType(storm::solver::EquationSolverType::Gmmxx);
+            env.solver().gmmxx().setMethod(storm::solver::GmmxxLinearEquationSolverMethod::Gmres);
+            env.solver().gmmxx().setPreconditioner(storm::solver::GmmxxLinearEquationSolverPreconditioner::Ilu);
+            env.solver().gmmxx().setPrecision(storm::utility::convertNumber<storm::RationalNumber>(1e-8));
+            env.solver().lra().setPrecision(storm::utility::convertNumber<storm::RationalNumber>(1e-8));
             return env;
         }
     };
     
-    class NativeWalkerChaeEnvironment {
+    class DistrEigenRationalLUEnvironment {
+    public:
+        typedef storm::RationalNumber ValueType;
+        static const bool isExact = true;
+        static storm::Environment createEnvironment() {
+            storm::Environment env;
+            env.solver().lra().setDetLraMethod(storm::solver::LraMethod::LraDistributionEquations);
+            env.solver().setLinearEquationSolverType(storm::solver::EquationSolverType::Eigen);
+            env.solver().eigen().setMethod(storm::solver::EigenLinearEquationSolverMethod::SparseLU);
+            return env;
+        }
+    };
+    
+    class DistrNativeWalkerChaeEnvironment {
     public:
         typedef double ValueType;
         static const bool isExact = false;
         static storm::Environment createEnvironment() {
             storm::Environment env;
+            env.solver().lra().setDetLraMethod(storm::solver::LraMethod::GainBiasEquations);
             env.solver().setLinearEquationSolverType(storm::solver::EquationSolverType::Native);
             env.solver().native().setMethod(storm::solver::NativeLinearEquationSolverMethod::WalkerChae);
+            env.solver().lra().setPrecision(storm::utility::convertNumber<storm::RationalNumber>(1e-8));
             env.solver().native().setPrecision(storm::utility::convertNumber<storm::RationalNumber>(1e-8));
             env.solver().native().setMaximalNumberOfIterations(50000);
             return env;
         }
     };
     
+    class ValueIterationEnvironment {
+    public:
+        typedef double ValueType;
+        static const bool isExact = false;
+        static storm::Environment createEnvironment() {
+            storm::Environment env;
+            env.solver().lra().setDetLraMethod(storm::solver::LraMethod::ValueIteration);
+            env.solver().lra().setPrecision(storm::utility::convertNumber<storm::RationalNumber>(1e-8));
+            return env;
+        }
+    };
+    
     template<typename TestType>
     class LraDtmcPrctlModelCheckerTest : public ::testing::Test {
     public:
@@ -104,11 +171,15 @@ namespace {
     };
   
     typedef ::testing::Types<
-            GmmxxDoubleGmresEnvironment,
-            EigenDoubleDGmresEnvironment,
-            EigenRationalLUEnvironment,
-            NativeSorEnvironment,
-            NativeWalkerChaeEnvironment
+            GBGmmxxDoubleGmresEnvironment,
+            GBEigenDoubleDGmresEnvironment,
+            GBEigenRationalLUEnvironment,
+            GBNativeSorEnvironment,
+            GBNativeWalkerChaeEnvironment,
+            DistrGmmxxDoubleGmresEnvironment,
+            DistrEigenRationalLUEnvironment,
+            DistrNativeWalkerChaeEnvironment,
+            ValueIterationEnvironment
     > TestingTypes;
     
    TYPED_TEST_SUITE(LraDtmcPrctlModelCheckerTest, TestingTypes,);

src/test/storm/modelchecker/prctl/mdp/LraMdpPrctlModelCheckerTest.cpp

@@ -22,7 +22,7 @@
 
 #include "storm/settings/modules/NativeEquationSolverSettings.h"
 #include "storm-parsers/parser/AutoParser.h"
-#include "storm/environment/solver/MinMaxSolverEnvironment.h"
+#include "storm/environment/solver/LongRunAverageSolverEnvironment.h"
 
 namespace {
     
@@ -34,8 +34,8 @@ namespace {
         typedef storm::models::sparse::Mdp<ValueType> ModelType;
         static storm::Environment createEnvironment() {
             storm::Environment env;
-            env.solver().minMax().setLraMethod(storm::solver::LraMethod::ValueIteration);
-            env.solver().minMax().setPrecision(storm::utility::convertNumber<storm::RationalNumber>(1e-10));
+            env.solver().lra().setNondetLraMethod(storm::solver::LraMethod::ValueIteration);
+            env.solver().lra().setPrecision(storm::utility::convertNumber<storm::RationalNumber>(1e-10));
             return env;
         }
     };
@@ -47,8 +47,8 @@ namespace {
         typedef storm::models::sparse::Mdp<ValueType> ModelType;
         static storm::Environment createEnvironment() {
             storm::Environment env;
-            env.solver().minMax().setLraMethod(storm::solver::LraMethod::LinearProgramming);
-            env.solver().minMax().setPrecision(storm::utility::convertNumber<storm::RationalNumber>(1e-10));
+            env.solver().lra().setNondetLraMethod(storm::solver::LraMethod::LinearProgramming);
+            env.solver().lra().setPrecision(storm::utility::convertNumber<storm::RationalNumber>(1e-10));
             return env;
         }
     };
@@ -60,7 +60,7 @@ namespace {
         typedef storm::models::sparse::Mdp<ValueType> ModelType;
         static storm::Environment createEnvironment() {
             storm::Environment env;
-            env.solver().minMax().setLraMethod(storm::solver::LraMethod::LinearProgramming);
+            env.solver().lra().setNondetLraMethod(storm::solver::LraMethod::LinearProgramming);
             return env;
         }
     };

version.cmake

@@ -1,4 +1,4 @@
 set(STORM_VERSION_MAJOR 1)
 set(STORM_VERSION_MINOR 4)
-set(STORM_VERSION_PATCH 0)
+set(STORM_VERSION_PATCH 1)