More tests, bugfixes: All tests pass.

Former-commit-id: f37c02a9d7

src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp

@@ -341,18 +341,19 @@ namespace storm {
 
 			// Get some data members for convenience.
 			typename storm::storage::SparseMatrix<ValueType> const& transitionMatrix = this->getModel().getTransitionMatrix();
+			std::vector<uint_fast64_t> const& nondeterministicChoiceIndices = this->getModel().getNondeterministicChoiceIndices();
 			ValueType one = storm::utility::one<ValueType>();
 			ValueType zero = storm::utility::zero<ValueType>();
 
 			//first calculate LRA for the Maximal End Components.
 			storm::storage::BitVector statesInMecs(numOfStates);
 			std::vector<uint_fast64_t> stateToMecIndexMap(transitionMatrix.getColumnCount());
-			std::vector<ValueType> mecLra(mecDecomposition.size(), zero);
+			std::vector<ValueType> lraValuesForEndComponents(mecDecomposition.size(), zero);
 
 			for (uint_fast64_t currentMecIndex = 0; currentMecIndex < mecDecomposition.size(); ++currentMecIndex) {
 				storm::storage::MaximalEndComponent const& mec = mecDecomposition[currentMecIndex];
 
-				mecLra[currentMecIndex] = computeLraForMaximalEndComponent(minimize, transitionMatrix, psiStates, mec);
+				lraValuesForEndComponents[currentMecIndex] = computeLraForMaximalEndComponent(minimize, transitionMatrix, psiStates, mec);
 
 				// Gather information for later use.
 				for (auto const& stateChoicesPair : mec) {
@@ -361,66 +362,122 @@ namespace storm {
 				}
 			}
 
-			//calculate LRA for states not in MECs as expected reachability rewards
-			//we add an auxiliary target state, every state in any MEC has a choice to move to that state with prob 1.
-			//This transitions have the LRA of the MEC as reward.
-			//The expected reward corresponds to sum of LRAs in MEC weighted by the reachability probability of the MEC
-
-			//we now build the submatrix of the transition matrix of the system with the auxiliary state, that only contains the states from
-			//the original state, i.e. all "maybe-states"
-			storm::storage::SparseMatrixBuilder<ValueType> rewardEquationSystemMatrixBuilder(transitionMatrix.getRowCount() + statesInMecs.getNumberOfSetBits(),
-				transitionMatrix.getColumnCount(),
-				transitionMatrix.getEntryCount(),
-				false,
-				true,
-				transitionMatrix.getRowGroupCount());
-
-			std::vector<ValueType> rewardRightSide(transitionMatrix.getRowCount() + statesInMecs.getNumberOfSetBits(), zero);
-
-			uint_fast64_t rowIndex = 0;
-			uint_fast64_t oldRowIndex = 0;
-			for (uint_fast64_t rowGroupIndex = 0; rowGroupIndex < transitionMatrix.getRowGroupCount(); ++rowGroupIndex) {
-				rewardEquationSystemMatrixBuilder.newRowGroup(rowIndex);
-				for (uint_fast64_t i = 0; i < transitionMatrix.getRowGroupSize(rowGroupIndex); ++i) {
-					//we have to make sure that an entry exists for all diagonal elements, even if it is zero. Other wise the call to convertToEquationSystem will produce wrong results or fail.
-					bool foundDiagonal = false;
-					for (auto entry : transitionMatrix.getRow(oldRowIndex)) {
-						if (!foundDiagonal) {
-							if (entry.getColumn() > rowGroupIndex) {
-								foundDiagonal = true;
-								rewardEquationSystemMatrixBuilder.addNextValue(rowIndex, rowGroupIndex, zero);
-							} else if (entry.getColumn() == rowGroupIndex) {
-								foundDiagonal = true;
-							}
+			// For fast transition rewriting, we build some auxiliary data structures.
+			storm::storage::BitVector statesNotContainedInAnyMec = ~statesInMecs;
+			uint_fast64_t firstAuxiliaryStateIndex = statesNotContainedInAnyMec.getNumberOfSetBits();
+			uint_fast64_t lastStateNotInMecs = 0;
+			uint_fast64_t numberOfStatesNotInMecs = 0;
+			std::vector<uint_fast64_t> statesNotInMecsBeforeIndex;
+			statesNotInMecsBeforeIndex.reserve(this->getModel().getNumberOfStates());
+			for (auto state : statesNotContainedInAnyMec) {
+				while (lastStateNotInMecs <= state) {
+					statesNotInMecsBeforeIndex.push_back(numberOfStatesNotInMecs);
+					++lastStateNotInMecs;
+				}
+				++numberOfStatesNotInMecs;
+			}
+
+			// Finally, we are ready to create the SSP matrix and right-hand side of the SSP.
+			std::vector<ValueType> b;
+			typename storm::storage::SparseMatrixBuilder<ValueType> sspMatrixBuilder(0, 0, 0, false, true, numberOfStatesNotInMecs + mecDecomposition.size());
+
+			// If the source state is not contained in any MEC, we copy its choices (and perform the necessary modifications).
+			uint_fast64_t currentChoice = 0;
+			for (auto state : statesNotContainedInAnyMec) {
+				sspMatrixBuilder.newRowGroup(currentChoice);
+
+				for (uint_fast64_t choice = nondeterministicChoiceIndices[state]; choice < nondeterministicChoiceIndices[state + 1]; ++choice, ++currentChoice) {
+					std::vector<ValueType> auxiliaryStateToProbabilityMap(mecDecomposition.size());
+					b.push_back(storm::utility::zero<ValueType>());
+
+					for (auto element : transitionMatrix.getRow(choice)) {
+						if (statesNotContainedInAnyMec.get(element.getColumn())) {
+							// If the target state is not contained in an MEC, we can copy over the entry.
+							sspMatrixBuilder.addNextValue(currentChoice, statesNotInMecsBeforeIndex[element.getColumn()], element.getValue());
+						} else {
+							// If the target state is contained in MEC i, we need to add the probability to the corresponding field in the vector
+							// so that we are able to write the cumulative probability to the MEC into the matrix.
+							auxiliaryStateToProbabilityMap[stateToMecIndexMap[element.getColumn()]] += element.getValue();
 						}
-						//copy over values from transition matrix of the actual system
-						rewardEquationSystemMatrixBuilder.addNextValue(rowIndex, entry.getColumn(), entry.getValue());
 					}
-					if (!foundDiagonal) {
-						rewardEquationSystemMatrixBuilder.addNextValue(rowIndex, rowGroupIndex, zero);
+
+					// Now insert all (cumulative) probability values that target an MEC.
+					for (uint_fast64_t mecIndex = 0; mecIndex < auxiliaryStateToProbabilityMap.size(); ++mecIndex) {
+						if (auxiliaryStateToProbabilityMap[mecIndex] != 0) {
+							sspMatrixBuilder.addNextValue(currentChoice, firstAuxiliaryStateIndex + mecIndex, auxiliaryStateToProbabilityMap[mecIndex]);
+						}
 					}
-					++oldRowIndex;
-					++rowIndex;
-				}
-				if (statesInMecs.get(rowGroupIndex)) {
-					//put the transition-reward on the right side
-					rewardRightSide[rowIndex] = mecLra[stateToMecIndexMap[rowGroupIndex]];
-					//add the choice where we go to the auxiliary state, which is a row with all zeros in the submatrix we build
-					rewardEquationSystemMatrixBuilder.addNextValue(rowIndex, rowGroupIndex, zero);
-					++rowIndex;
 				}
 			}
 
-			storm::storage::SparseMatrix<ValueType> rewardEquationSystemMatrix = rewardEquationSystemMatrixBuilder.build();
-			rewardEquationSystemMatrix.convertToEquationSystem();
+			// Now we are ready to construct the choices for the auxiliary states.
+			for (uint_fast64_t mecIndex = 0; mecIndex < mecDecomposition.size(); ++mecIndex) {
+				storm::storage::MaximalEndComponent const& mec = mecDecomposition[mecIndex];
+				sspMatrixBuilder.newRowGroup(currentChoice);
+
+				for (auto const& stateChoicesPair : mec) {
+					uint_fast64_t state = stateChoicesPair.first;
+					boost::container::flat_set<uint_fast64_t> const& choicesInMec = stateChoicesPair.second;
+
+					for (uint_fast64_t choice = nondeterministicChoiceIndices[state]; choice < nondeterministicChoiceIndices[state + 1]; ++choice) {
+						std::vector<ValueType> auxiliaryStateToProbabilityMap(mecDecomposition.size());
+
+						// If the choice is not contained in the MEC itself, we have to add a similar distribution to the auxiliary state.
+						if (choicesInMec.find(choice) == choicesInMec.end()) {
+							b.push_back(storm::utility::zero<ValueType>());
+
+							for (auto element : transitionMatrix.getRow(choice)) {
+								if (statesNotContainedInAnyMec.get(element.getColumn())) {
+									// If the target state is not contained in an MEC, we can copy over the entry.
+									sspMatrixBuilder.addNextValue(currentChoice, statesNotInMecsBeforeIndex[element.getColumn()], element.getValue());
+								} else {
+									// If the target state is contained in MEC i, we need to add the probability to the corresponding field in the vector
+									// so that we are able to write the cumulative probability to the MEC into the matrix.
+									auxiliaryStateToProbabilityMap[stateToMecIndexMap[element.getColumn()]] += element.getValue();
+								}
+							}
+
+							// Now insert all (cumulative) probability values that target an MEC.
+							for (uint_fast64_t targetMecIndex = 0; targetMecIndex < auxiliaryStateToProbabilityMap.size(); ++targetMecIndex) {
+								if (auxiliaryStateToProbabilityMap[targetMecIndex] != 0) {
+									// If the target MEC is the same as the current one, instead of adding a transition, we need to add the weighted reward
+									// to the right-hand side vector of the SSP.
+									if (mecIndex == targetMecIndex) {
+										b.back() += auxiliaryStateToProbabilityMap[mecIndex] * lraValuesForEndComponents[mecIndex];
+									} else {
+										// Otherwise, we add a transition to the auxiliary state that is associated with the target MEC.
+										sspMatrixBuilder.addNextValue(currentChoice, firstAuxiliaryStateIndex + targetMecIndex, auxiliaryStateToProbabilityMap[targetMecIndex]);
+									}
+								}
+							}
 
-			std::vector<ValueType> result(rewardEquationSystemMatrix.getColumnCount(), one);
+							++currentChoice;
+						}
+					}
+				}
 
-			{
-				auto solver = this->MinMaxLinearEquationSolverFactory->create(rewardEquationSystemMatrix);
-				solver->solveEquationSystem(minimize, result, rewardRightSide);
+				// For each auxiliary state, there is the option to achieve the reward value of the LRA associated with the MEC.
+				++currentChoice;
+				b.push_back(lraValuesForEndComponents[mecIndex]);
 			}
 
+			// Finalize the matrix and solve the corresponding system of equations.
+			storm::storage::SparseMatrix<ValueType> sspMatrix = sspMatrixBuilder.build(currentChoice);
+
+			std::vector<ValueType> sspResult(numberOfStatesNotInMecs + mecDecomposition.size());
+			std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = MinMaxLinearEquationSolverFactory->create(sspMatrix);
+			solver->solveEquationSystem(minimize, sspResult, b);
+
+			// Prepare result vector.
+			std::vector<ValueType> result(this->getModel().getNumberOfStates());
+
+			// Set the values for states not contained in MECs.
+			storm::utility::vector::setVectorValues(result, statesNotContainedInAnyMec, sspResult);
+
+			// Set the values for all states in MECs.
+			for (auto state : statesInMecs) {
+				result[state] = sspResult[firstAuxiliaryStateIndex + stateToMecIndexMap[state]];
+			}
 
 			return result;
 		}
@@ -455,16 +512,16 @@ namespace storm {
 
 				// Now, based on the type of the state, create a suitable constraint.
 				for (auto choice : stateChoicesPair.second) {
-					storm::expressions::Expression constraint = solver->getConstant(1);
-					ValueType w = 0;
+					storm::expressions::Expression constraint = -lambda;
+					ValueType r = 0;
 
 					for (auto element : transitionMatrix.getRow(choice)) {
 						constraint = constraint + stateToVariableMap.at(element.getColumn()) * solver->getConstant(element.getValue());
 						if (psiStates.get(element.getColumn())) {
-							w += element.getValue();
+							r += element.getValue();
 						}
 					}
-					constraint = constraint - solver->getConstant(w) * lambda;
+					constraint = solver->getConstant(r) + constraint;
 
 					if (minimize) {
 						constraint = stateToVariableMap.at(state) <= constraint;

test/functional/modelchecker/SparseMdpPrctlModelCheckerTest.cpp

@@ -195,7 +195,7 @@ TEST(SparseMdpPrctlModelCheckerTest, AsynchronousLeader) {
 	EXPECT_NEAR(4.285689611, quantitativeResult6[0], storm::settings::nativeEquationSolverSettings().getPrecision());
 }
 
-TEST(SparseMdpPrctlModelCheckerTest, LRA) {
+TEST(SparseMdpPrctlModelCheckerTest, LRA_SingleMec) {
 	storm::storage::SparseMatrixBuilder<double> matrixBuilder;
 	std::shared_ptr<storm::models::sparse::Mdp<double>> mdp;
 
@@ -216,10 +216,314 @@ TEST(SparseMdpPrctlModelCheckerTest, LRA) {
 		auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("a");
 		auto lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(storm::logic::OptimalityType::Maximize, labelFormula);
 
-		std::unique_ptr<storm::modelchecker::CheckResult> result = checker.check(*lraFormula);
-		storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult = result->asExplicitQuantitativeCheckResult<double>();
+		std::unique_ptr<storm::modelchecker::CheckResult> result = std::move(checker.check(*lraFormula));
+		storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult1 = result->asExplicitQuantitativeCheckResult<double>();
+
+		EXPECT_NEAR(.5, quantitativeResult1[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(.5, quantitativeResult1[1], storm::settings::nativeEquationSolverSettings().getPrecision());
+
+		lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(storm::logic::OptimalityType::Minimize, labelFormula);
+
+		result = std::move(checker.check(*lraFormula));
+		storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult2 = result->asExplicitQuantitativeCheckResult<double>();
+
+		EXPECT_NEAR(.5, quantitativeResult2[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(.5, quantitativeResult2[1], storm::settings::nativeEquationSolverSettings().getPrecision());
+	}
+	{
+		matrixBuilder = storm::storage::SparseMatrixBuilder<double>(2, 2, 4);
+		matrixBuilder.addNextValue(0, 0, .5);
+		matrixBuilder.addNextValue(0, 1, .5);
+		matrixBuilder.addNextValue(1, 0, .5);
+		matrixBuilder.addNextValue(1, 1, .5);
+		storm::storage::SparseMatrix<double> transitionMatrix = matrixBuilder.build();
+
+		storm::models::sparse::StateLabeling ap(2);
+		ap.addLabel("a");
+		ap.addLabelToState("a", 1);
+
+		mdp.reset(new storm::models::sparse::Mdp<double>(transitionMatrix, ap, boost::none, boost::none, boost::none));
+
+		storm::modelchecker::SparseMdpPrctlModelChecker<double> checker(*mdp, std::unique_ptr<storm::utility::solver::MinMaxLinearEquationSolverFactory<double>>(new storm::utility::solver::NativeMinMaxLinearEquationSolverFactory<double>()));
+
+		auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("a");
+		auto lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(storm::logic::OptimalityType::Maximize, labelFormula);
+
+		std::unique_ptr<storm::modelchecker::CheckResult> result = std::move(checker.check(*lraFormula));
+		storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult1 = result->asExplicitQuantitativeCheckResult<double>();
+
+		EXPECT_NEAR(.5, quantitativeResult1[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(.5, quantitativeResult1[1], storm::settings::nativeEquationSolverSettings().getPrecision());
+
+		lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(storm::logic::OptimalityType::Minimize, labelFormula);
+
+		result = std::move(checker.check(*lraFormula));
+		storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult2 = result->asExplicitQuantitativeCheckResult<double>();
+
+		EXPECT_NEAR(.5, quantitativeResult2[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(.5, quantitativeResult2[1], storm::settings::nativeEquationSolverSettings().getPrecision());
+	}
+
+	{
+		matrixBuilder = storm::storage::SparseMatrixBuilder<double>(4, 3, 4, true, true, 3);
+		matrixBuilder.newRowGroup(0);
+		matrixBuilder.addNextValue(0, 1, 1);
+		matrixBuilder.newRowGroup(1);
+		matrixBuilder.addNextValue(1, 0, 1);
+		matrixBuilder.addNextValue(2, 2, 1);
+		matrixBuilder.newRowGroup(3);
+		matrixBuilder.addNextValue(3, 0, 1);
+		storm::storage::SparseMatrix<double> transitionMatrix = matrixBuilder.build();
+
+		storm::models::sparse::StateLabeling ap(3);
+		ap.addLabel("a");
+		ap.addLabelToState("a", 2);
+		ap.addLabel("b");
+		ap.addLabelToState("b", 0);
+		ap.addLabel("c");
+		ap.addLabelToState("c", 0);
+		ap.addLabelToState("c", 2);
+
+		mdp.reset(new storm::models::sparse::Mdp<double>(transitionMatrix, ap, boost::none, boost::none, boost::none));
+
+		storm::modelchecker::SparseMdpPrctlModelChecker<double> checker(*mdp, std::unique_ptr<storm::utility::solver::MinMaxLinearEquationSolverFactory<double>>(new storm::utility::solver::NativeMinMaxLinearEquationSolverFactory<double>()));
+
+		auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("a");
+		auto lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(storm::logic::OptimalityType::Maximize, labelFormula);
+
+		std::unique_ptr<storm::modelchecker::CheckResult> result = std::move(checker.check(*lraFormula));
+		storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult1 = result->asExplicitQuantitativeCheckResult<double>();
+
+		EXPECT_NEAR(1. / 3., quantitativeResult1[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(1. / 3., quantitativeResult1[1], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(1. / 3., quantitativeResult1[2], storm::settings::nativeEquationSolverSettings().getPrecision());
+
+		lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(storm::logic::OptimalityType::Minimize, labelFormula);
+
+		result = std::move(checker.check(*lraFormula));
+		storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult2 = result->asExplicitQuantitativeCheckResult<double>();
+
+		EXPECT_NEAR(0.0, quantitativeResult2[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(0.0, quantitativeResult2[1], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(0.0, quantitativeResult2[2], storm::settings::nativeEquationSolverSettings().getPrecision());
+
+		labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("b");
+		lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(storm::logic::OptimalityType::Maximize, labelFormula);
+
+		result = std::move(checker.check(*lraFormula));
+		storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult3 = result->asExplicitQuantitativeCheckResult<double>();
+
+		EXPECT_NEAR(0.5, quantitativeResult3[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(0.5, quantitativeResult3[1], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(0.5, quantitativeResult3[2], storm::settings::nativeEquationSolverSettings().getPrecision());
+
+		lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(storm::logic::OptimalityType::Minimize, labelFormula);
+
+		result = std::move(checker.check(*lraFormula));
+		storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult4 = result->asExplicitQuantitativeCheckResult<double>();
+
+		EXPECT_NEAR(1. / 3., quantitativeResult4[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(1. / 3., quantitativeResult4[1], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(1. / 3., quantitativeResult4[2], storm::settings::nativeEquationSolverSettings().getPrecision());
+
+		labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("c");
+		lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(storm::logic::OptimalityType::Maximize, labelFormula);
+
+		result = std::move(checker.check(*lraFormula));
+		storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult5 = result->asExplicitQuantitativeCheckResult<double>();
+
+		EXPECT_NEAR(2. / 3., quantitativeResult5[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(2. / 3., quantitativeResult5[1], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(2. / 3., quantitativeResult5[2], storm::settings::nativeEquationSolverSettings().getPrecision());
+
+		lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(storm::logic::OptimalityType::Minimize, labelFormula);
+
+		result = std::move(checker.check(*lraFormula));
+		storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult6 = result->asExplicitQuantitativeCheckResult<double>();
+
+		EXPECT_NEAR(0.5, quantitativeResult6[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(0.5, quantitativeResult6[1], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(0.5, quantitativeResult6[2], storm::settings::nativeEquationSolverSettings().getPrecision());
+	}
+}
+
+TEST(SparseMdpPrctlModelCheckerTest, LRA) {
+	storm::storage::SparseMatrixBuilder<double> matrixBuilder;
+	std::shared_ptr<storm::models::sparse::Mdp<double>> mdp;
+
+	{
+		matrixBuilder = storm::storage::SparseMatrixBuilder<double>(4, 3, 4, true, true, 3);
+		matrixBuilder.newRowGroup(0);
+		matrixBuilder.addNextValue(0, 1, 1);
+		matrixBuilder.newRowGroup(1);
+		matrixBuilder.addNextValue(1, 1, 1);
+		matrixBuilder.addNextValue(2, 2, 1);
+		matrixBuilder.newRowGroup(3);
+		matrixBuilder.addNextValue(3, 2, 1);
+		storm::storage::SparseMatrix<double> transitionMatrix = matrixBuilder.build();
+
+		storm::models::sparse::StateLabeling ap(3);
+		ap.addLabel("a");
+		ap.addLabelToState("a", 0);
+		ap.addLabel("b");
+		ap.addLabelToState("b", 1);
+		ap.addLabel("c");
+		ap.addLabelToState("c", 2);
+
+		mdp.reset(new storm::models::sparse::Mdp<double>(transitionMatrix, ap, boost::none, boost::none, boost::none));
+
+		storm::modelchecker::SparseMdpPrctlModelChecker<double> checker(*mdp, std::unique_ptr<storm::utility::solver::MinMaxLinearEquationSolverFactory<double>>(new storm::utility::solver::NativeMinMaxLinearEquationSolverFactory<double>()));
+
+		auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("a");
+		auto lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(storm::logic::OptimalityType::Maximize, labelFormula);
+
+		std::unique_ptr<storm::modelchecker::CheckResult> result = std::move(checker.check(*lraFormula));
+		storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult1 = result->asExplicitQuantitativeCheckResult<double>();
+
+		EXPECT_NEAR(0.0, quantitativeResult1[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(0.0, quantitativeResult1[1], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(0.0, quantitativeResult1[2], storm::settings::nativeEquationSolverSettings().getPrecision());
+
+		lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(storm::logic::OptimalityType::Minimize, labelFormula);
+
+		result = std::move(checker.check(*lraFormula));
+		storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult2 = result->asExplicitQuantitativeCheckResult<double>();
+
+		EXPECT_NEAR(0.0, quantitativeResult2[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(0.0, quantitativeResult2[1], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(0.0, quantitativeResult2[2], storm::settings::nativeEquationSolverSettings().getPrecision());
+
+		labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("b");
+		lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(storm::logic::OptimalityType::Maximize, labelFormula);
+
+		result = std::move(checker.check(*lraFormula));
+		storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult3 = result->asExplicitQuantitativeCheckResult<double>();
+
+		EXPECT_NEAR(1.0, quantitativeResult3[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(1.0, quantitativeResult3[1], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(0.0, quantitativeResult3[2], storm::settings::nativeEquationSolverSettings().getPrecision());
+
+		lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(storm::logic::OptimalityType::Minimize, labelFormula);
+
+		result = std::move(checker.check(*lraFormula));
+		storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult4 = result->asExplicitQuantitativeCheckResult<double>();
+
+		EXPECT_NEAR(0.0, quantitativeResult4[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(0.0, quantitativeResult4[1], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(0.0, quantitativeResult4[2], storm::settings::nativeEquationSolverSettings().getPrecision());
+
+		labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("c");
+		lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(storm::logic::OptimalityType::Maximize, labelFormula);
+
+		result = std::move(checker.check(*lraFormula));
+		storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult5 = result->asExplicitQuantitativeCheckResult<double>();
+
+		EXPECT_NEAR(1.0, quantitativeResult5[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(1.0, quantitativeResult5[1], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(1.0, quantitativeResult5[2], storm::settings::nativeEquationSolverSettings().getPrecision());
+
+		lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(storm::logic::OptimalityType::Minimize, labelFormula);
+
+		result = std::move(checker.check(*lraFormula));
+		storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult6 = result->asExplicitQuantitativeCheckResult<double>();
+
+		EXPECT_NEAR(0.0, quantitativeResult6[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(0.0, quantitativeResult6[1], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(1.0, quantitativeResult6[2], storm::settings::nativeEquationSolverSettings().getPrecision());
+	}
+	{
+		matrixBuilder = storm::storage::SparseMatrixBuilder<double>(22, 15, 28, true, true, 15);
+		matrixBuilder.newRowGroup(0);
+		matrixBuilder.addNextValue(0, 1, 1);
+		matrixBuilder.newRowGroup(1);
+		matrixBuilder.addNextValue(1, 0, 1);
+		matrixBuilder.addNextValue(2, 2, 1);
+		matrixBuilder.addNextValue(3, 4, 0.7);
+		matrixBuilder.addNextValue(3, 6, 0.3);
+		matrixBuilder.newRowGroup(4);
+		matrixBuilder.addNextValue(4, 0, 1);
+
+		matrixBuilder.newRowGroup(5);
+		matrixBuilder.addNextValue(5, 4, 1);
+		matrixBuilder.addNextValue(6, 5, 0.8);
+		matrixBuilder.addNextValue(6, 9, 0.2);
+		matrixBuilder.newRowGroup(7);
+		matrixBuilder.addNextValue(7, 3, 1);
+		matrixBuilder.addNextValue(8, 5, 1);
+		matrixBuilder.newRowGroup(9);
+		matrixBuilder.addNextValue(9, 3, 1);
+
+		matrixBuilder.newRowGroup(10);
+		matrixBuilder.addNextValue(10, 7, 1);
+		matrixBuilder.newRowGroup(11);
+		matrixBuilder.addNextValue(11, 6, 1);
+		matrixBuilder.addNextValue(12, 8, 1);
+		matrixBuilder.newRowGroup(13);
+		matrixBuilder.addNextValue(13, 6, 1);
+
+		matrixBuilder.newRowGroup(14);
+		matrixBuilder.addNextValue(14, 10, 1);
+		matrixBuilder.newRowGroup(15);
+		matrixBuilder.addNextValue(15, 9, 1);
+		matrixBuilder.addNextValue(16, 11, 1);
+		matrixBuilder.newRowGroup(17);
+		matrixBuilder.addNextValue(17, 9, 1);
+
+		matrixBuilder.newRowGroup(18);
+		matrixBuilder.addNextValue(18, 5, 0.4);
+		matrixBuilder.addNextValue(18, 8, 0.3);
+		matrixBuilder.addNextValue(18, 11, 0.3);
+
+		matrixBuilder.newRowGroup(19);
+		matrixBuilder.addNextValue(19, 7, 0.7);
+		matrixBuilder.addNextValue(19, 12, 0.3);
+
+		matrixBuilder.newRowGroup(20);
+		matrixBuilder.addNextValue(20, 12, 0.1);
+		matrixBuilder.addNextValue(20, 13, 0.9);
+		matrixBuilder.addNextValue(21, 12, 1);
+
+		storm::storage::SparseMatrix<double> transitionMatrix = matrixBuilder.build();
+
+		storm::models::sparse::StateLabeling ap(15);
+		ap.addLabel("a");
+		ap.addLabelToState("a", 1);
+		ap.addLabelToState("a", 4);
+		ap.addLabelToState("a", 5);
+		ap.addLabelToState("a", 7);
+		ap.addLabelToState("a", 11);
+		ap.addLabelToState("a", 13);
+		ap.addLabelToState("a", 14);
+
+		mdp.reset(new storm::models::sparse::Mdp<double>(transitionMatrix, ap, boost::none, boost::none, boost::none));
+
+		storm::modelchecker::SparseMdpPrctlModelChecker<double> checker(*mdp, std::unique_ptr<storm::utility::solver::MinMaxLinearEquationSolverFactory<double>>(new storm::utility::solver::NativeMinMaxLinearEquationSolverFactory<double>()));
+
+		auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("a");
+		auto lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(storm::logic::OptimalityType::Maximize, labelFormula);
 
-		EXPECT_NEAR(.5, quantitativeResult[0], storm::settings::nativeEquationSolverSettings().getPrecision());
-		EXPECT_NEAR(.5, quantitativeResult[1], storm::settings::nativeEquationSolverSettings().getPrecision());
+		std::unique_ptr<storm::modelchecker::CheckResult> result = std::move(checker.check(*lraFormula));
+		storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult1 = result->asExplicitQuantitativeCheckResult<double>();
+
+		EXPECT_NEAR(37./60., quantitativeResult1[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(2./3., quantitativeResult1[3], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(0.5, quantitativeResult1[6], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(1./3., quantitativeResult1[9], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(31./60., quantitativeResult1[12], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(101./200., quantitativeResult1[13], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(31./60., quantitativeResult1[14], storm::settings::nativeEquationSolverSettings().getPrecision());
+
+		lraFormula = std::make_shared<storm::logic::LongRunAverageOperatorFormula>(storm::logic::OptimalityType::Minimize, labelFormula);
+
+		result = std::move(checker.check(*lraFormula));
+		storm::modelchecker::ExplicitQuantitativeCheckResult<double>& quantitativeResult2 = result->asExplicitQuantitativeCheckResult<double>();
+
+		EXPECT_NEAR(1./3., quantitativeResult2[0], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(0.4, quantitativeResult2[3], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(1./3., quantitativeResult2[6], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(0.0, quantitativeResult2[9], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(.26, quantitativeResult2[12], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(467./1500., quantitativeResult2[13], storm::settings::nativeEquationSolverSettings().getPrecision());
+		EXPECT_NEAR(.26, quantitativeResult2[14], storm::settings::nativeEquationSolverSettings().getPrecision());
 	}
 }