From 39abecbad3788d81dfbaa72db32495d94ea29835 Mon Sep 17 00:00:00 2001 From: dehnert Date: Mon, 22 Jun 2015 14:33:29 +0200 Subject: [PATCH] added some tests for LRA in CTMCs Former-commit-id: 3b847d542e52a1ab3c65f551aa83e0d332527049 --- .../csl/SparseCtmcCslModelChecker.cpp | 212 ++++++++++-------- .../GmmxxCtmcCslModelCheckerTest.cpp | 21 ++ 2 files changed, 142 insertions(+), 91 deletions(-) diff --git a/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp b/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp index 62ca9af85..e1dd96467 100644 --- a/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp +++ b/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp @@ -191,7 +191,7 @@ namespace storm { std::vector newSubresult = std::vector(relevantStates.getNumberOfSetBits()); storm::utility::vector::setVectorValues(newSubresult, statesWithProbabilityGreater0NonPsi % relevantStates, subresult); storm::utility::vector::setVectorValues(newSubresult, psiStates % relevantStates, storm::utility::one()); - + // Then compute the transient probabilities of being in such a state after t time units. For this, // we must re-uniformize the CTMC, so we need to compute the second uniformized matrix. uniformizationRate = storm::utility::zero(); @@ -346,7 +346,7 @@ namespace storm { weight = std::get<3>(foxGlynnResult)[index - std::get<0>(foxGlynnResult)]; storm::utility::vector::applyPointwise(result, values, result, addAndScale); } - + return result; } @@ -408,7 +408,7 @@ namespace storm { } uniformizationRate *= 1.02; STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive."); - + storm::storage::SparseMatrix uniformizedMatrix = this->computeUniformizedMatrix(this->getModel().getTransitionMatrix(), storm::storage::BitVector(this->getModel().getNumberOfStates(), true), uniformizationRate, this->getModel().getExitRateVector()); // Compute the total state reward vector. @@ -477,7 +477,7 @@ namespace storm { std::unique_ptr SparseCtmcCslModelChecker::computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative, boost::optional const& optimalityType) { std::unique_ptr subResultPointer = this->check(stateFormula); ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult(); - + storm::storage::SparseMatrix probabilityMatrix = this->computeProbabilityMatrix(this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector()); return std::unique_ptr(new ExplicitQuantitativeCheckResult(computeLongRunAverageHelper(probabilityMatrix, subResult.getTruthValuesVector(), &this->getModel().getExitRateVector(), qualitative, *linearEquationSolverFactory))); } @@ -502,12 +502,13 @@ namespace storm { ValueType one = storm::utility::one(); ValueType zero = storm::utility::zero(); + // Prepare the vector holding the LRA values for each of the BSCCs. + std::vector bsccLra(bsccDecomposition.size(), zero); + // First we check which states are in BSCCs. storm::storage::BitVector statesInBsccs(numOfStates); storm::storage::BitVector firstStatesInBsccs(numOfStates); - std::vector stateToBsccIndexMap(transitionMatrix.getColumnCount()); - for (uint_fast64_t currentBsccIndex = 0; currentBsccIndex < bsccDecomposition.size(); ++currentBsccIndex) { storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[currentBsccIndex]; @@ -518,99 +519,125 @@ namespace storm { if (first) { firstStatesInBsccs.set(state); } - stateToBsccIndexMap[state] = currentBsccIndex; first = false; } } - - firstStatesInBsccs = firstStatesInBsccs % statesInBsccs; storm::storage::BitVector statesNotInBsccs = ~statesInBsccs; - // Then we construct an equation system that yields the steady state probabilities for all states in BSCCs. - storm::storage::SparseMatrix bsccEquationSystem = transitionMatrix.getSubmatrix(false, statesInBsccs, statesInBsccs, true); + // Prepare a vector holding the index within all states that are in BSCCs for every state. + std::vector indexInStatesInBsccs; - // 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); + // Prepare a vector that maps the index within the set of all states in BSCCs to the index of the containing BSCC. + std::vector stateToBsccIndexMap; - // 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; - std::vector indexInStatesInBsccs; - indexInStatesInBsccs.reserve(transitionMatrix.getRowCount()); - for (auto index : statesInBsccs) { - while (lastIndex <= index) { - indexInStatesInBsccs.push_back(currentNumberOfSetBits); - ++lastIndex; + 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 bsccEquationSystem = transitionMatrix.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(transitionMatrix.getRowCount()); + for (auto index : statesInBsccs) { + while (lastIndex <= index) { + indexInStatesInBsccs.push_back(currentNumberOfSetBits); + ++lastIndex; + } + ++currentNumberOfSetBits; } - ++currentNumberOfSetBits; - } - - // Now build the final equation system matrix. - storm::storage::SparseMatrixBuilder builder; - uint_fast64_t currentBsccIndex = 0; - 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. - if (firstStatesInBsccs.get(row)) { + + 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) { - builder.addNextValue(row, indexInStatesInBsccs[state], one); + stateToBsccIndexMap[indexInStatesInBsccs[state]] = currentBsccIndex; } - ++currentBsccIndex; - } else { - // Otherwise, we copy the row, and subtract 1 from the diagonal. - for (auto& entry : bsccEquationSystem.getRow(row)) { - if (entry.getColumn() == row) { - builder.addNextValue(row, entry.getColumn(), entry.getValue() - one); - } else { - builder.addNextValue(row, entry.getColumn(), entry.getValue()); + } + + // Now build the final equation system matrix, the initial guess and the right-hand side in one go. + std::vector bsccEquationSystemRightSide(bsccEquationSystem.getColumnCount(), zero); + storm::storage::SparseMatrixBuilder 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]; + + 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. + for (auto& entry : bsccEquationSystem.getRow(row)) { + if (entry.getColumn() == row) { + builder.addNextValue(row, entry.getColumn(), entry.getValue() - one); + } else { + builder.addNextValue(row, entry.getColumn(), entry.getValue()); + } } } + } - } - - bsccEquationSystem = builder.build(); - - std::vector bsccEquationSystemRightSide(bsccEquationSystem.getColumnCount(), zero); - storm::utility::vector::setVectorValues(bsccEquationSystemRightSide, firstStatesInBsccs, one); - - // As an initial guess, we take the uniform distribution over all states of the containing BSCC. - std::vector bsccEquationSystemSolution(bsccEquationSystem.getColumnCount(), zero); - for (uint_fast64_t currentBsccIndex = 0; currentBsccIndex < bsccDecomposition.size(); ++currentBsccIndex) { - storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[currentBsccIndex]; - for (auto const& state : bscc) { - bsccEquationSystemSolution[indexInStatesInBsccs[state]] = one / bscc.size(); + // Create the initial guess for the LRAs. We take a uniform distribution over all states in a BSCC. + std::vector 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(); + } } - } - - { - std::unique_ptr> solver = linearEquationSolverFactory.create(bsccEquationSystem); - solver->solveEquationSystem(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 bsccTotalValue(bsccDecomposition.size(), zero); - std::size_t i = 0; - for (auto stateIter = statesInBsccs.begin(); stateIter != statesInBsccs.end(); ++i, ++stateIter) { - bsccTotalValue[stateToBsccIndexMap[*stateIter]] += bsccEquationSystemSolution[i] * (one / (*exitRateVector)[*stateIter]); + + bsccEquationSystem = builder.build(); + + { + std::unique_ptr> solver = linearEquationSolverFactory.create(bsccEquationSystem); + solver->solveEquationSystem(bsccEquationSystemSolution, bsccEquationSystemRightSide); } - i = 0; - for (auto stateIter = statesInBsccs.begin(); stateIter != statesInBsccs.end(); ++i, ++stateIter) { - bsccEquationSystemSolution[i] = (bsccEquationSystemSolution[i] * (one / (*exitRateVector)[*stateIter])) / bsccTotalValue[stateToBsccIndexMap[*stateIter]]; + // If exit rates were given, we need to 'fix' the results to also account for the timing behaviour. + if (exitRateVector != nullptr) { + std::vector 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. - std::vector bsccLra(bsccDecomposition.size(), zero); - size_t i = 0; - for (auto stateIter = statesInBsccs.begin(); stateIter != statesInBsccs.end(); ++i, ++stateIter) { - if (psiStates.get(*stateIter)) { - bsccLra[stateToBsccIndexMap[*stateIter]] += bsccEquationSystemSolution[i]; + // 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) { + if (psiStates.get(state)) { + bsccLra[stateToBsccIndexMap[indexInStatesInBsccs[state]]] += bsccEquationSystemSolution[indexInStatesInBsccs[state]]; + } + } + } + } 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). + if (psiStates.get(*bscc.begin())) { + bsccLra[bsccIndex] = 1; + } } } @@ -628,7 +655,7 @@ namespace storm { ValueType reward = zero; for (auto entry : transitionMatrix.getRow(state)) { if (statesInBsccs.get(entry.getColumn())) { - reward += entry.getValue() * bsccLra[stateToBsccIndexMap[entry.getColumn()]]; + reward += entry.getValue() * bsccLra[stateToBsccIndexMap[indexInStatesInBsccs[entry.getColumn()]]]; } } rewardRightSide.push_back(reward); @@ -648,18 +675,21 @@ namespace storm { // Fill the result vector. std::vector result(numOfStates); auto rewardSolutionIter = rewardSolution.begin(); - for (size_t state = 0; state < numOfStates; ++state) { - if (statesInBsccs.get(state)) { - // Assign the value of the bscc the state is in. - result[state] = bsccLra[stateToBsccIndexMap[state]]; - } else { - STORM_LOG_ASSERT(rewardSolutionIter != rewardSolution.end(), "Too few elements in solution."); - // Take the value from the reward computation. Since the n-th state not in any bscc is the n-th - // entry in rewardSolution we can just take the next value from the iterator. - result[state] = *rewardSolutionIter; - ++rewardSolutionIter; + + for (uint_fast64_t bsccIndex = 0; bsccIndex < bsccDecomposition.size(); ++bsccIndex) { + storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[bsccIndex]; + + for (auto const& state : bscc) { + result[state] = bsccLra[bsccIndex]; } } + for (auto state : statesNotInBsccs) { + STORM_LOG_ASSERT(rewardSolutionIter != rewardSolution.end(), "Too few elements in solution."); + // Take the value from the reward computation. Since the n-th state not in any bscc is the n-th + // entry in rewardSolution we can just take the next value from the iterator. + result[state] = *rewardSolutionIter; + ++rewardSolutionIter; + } return result; } diff --git a/test/functional/modelchecker/GmmxxCtmcCslModelCheckerTest.cpp b/test/functional/modelchecker/GmmxxCtmcCslModelCheckerTest.cpp index 4c79aa929..dc3d7a5c7 100644 --- a/test/functional/modelchecker/GmmxxCtmcCslModelCheckerTest.cpp +++ b/test/functional/modelchecker/GmmxxCtmcCslModelCheckerTest.cpp @@ -86,6 +86,13 @@ TEST(GmmxxCtmcCslModelCheckerTest, Cluster) { ASSERT_TRUE(checkResult->isExplicitQuantitativeCheckResult()); storm::modelchecker::ExplicitQuantitativeCheckResult quantitativeCheckResult7 = checkResult->asExplicitQuantitativeCheckResult(); EXPECT_NEAR(0.8602815057967503, quantitativeCheckResult7[initialState], storm::settings::generalSettings().getPrecision()); + + formula = formulaParser.parseFromString("LRA=? [\"minimum\"]"); + checkResult = modelchecker.check(*formula); + + ASSERT_TRUE(checkResult->isExplicitQuantitativeCheckResult()); + storm::modelchecker::ExplicitQuantitativeCheckResult quantitativeCheckResult8 = checkResult->asExplicitQuantitativeCheckResult(); + EXPECT_NEAR(0.99999766034263426, quantitativeCheckResult8[initialState], storm::settings::generalSettings().getPrecision()); } TEST(GmmxxCtmcCslModelCheckerTest, Embedded) { @@ -149,6 +156,13 @@ TEST(GmmxxCtmcCslModelCheckerTest, Embedded) { ASSERT_TRUE(checkResult->isExplicitQuantitativeCheckResult()); storm::modelchecker::ExplicitQuantitativeCheckResult quantitativeCheckResult5 = checkResult->asExplicitQuantitativeCheckResult(); EXPECT_NEAR(2.7745274082080154, quantitativeCheckResult5[initialState], storm::settings::generalSettings().getPrecision()); + + formula = formulaParser.parseFromString("LRA=? [\"fail_sensors\"]"); + checkResult = modelchecker.check(*formula); + + ASSERT_TRUE(checkResult->isExplicitQuantitativeCheckResult()); + storm::modelchecker::ExplicitQuantitativeCheckResult quantitativeCheckResult6 = checkResult->asExplicitQuantitativeCheckResult(); + EXPECT_NEAR(0.93458866427696596, quantitativeCheckResult6[initialState], storm::settings::generalSettings().getPrecision()); } TEST(GmmxxCtmcCslModelCheckerTest, Polling) { @@ -252,4 +266,11 @@ TEST(GmmxxCtmcCslModelCheckerTest, Tandem) { ASSERT_TRUE(checkResult->isExplicitQuantitativeCheckResult()); storm::modelchecker::ExplicitQuantitativeCheckResult quantitativeCheckResult6 = checkResult->asExplicitQuantitativeCheckResult(); EXPECT_NEAR(262.85103661561755, quantitativeCheckResult6[initialState], storm::settings::generalSettings().getPrecision()); + + formula = formulaParser.parseFromString("LRA=? [\"first_queue_full\"]"); + checkResult = modelchecker.check(*formula); + + ASSERT_TRUE(checkResult->isExplicitQuantitativeCheckResult()); + storm::modelchecker::ExplicitQuantitativeCheckResult quantitativeCheckResult7 = checkResult->asExplicitQuantitativeCheckResult(); + EXPECT_NEAR(0.9100373532, quantitativeCheckResult7[initialState], storm::settings::generalSettings().getPrecision()); }