diff --git a/src/storm/modelchecker/prctl/helper/HybridMdpPrctlHelper.cpp b/src/storm/modelchecker/prctl/helper/HybridMdpPrctlHelper.cpp index 171a7bf70..4d6f3426b 100644 --- a/src/storm/modelchecker/prctl/helper/HybridMdpPrctlHelper.cpp +++ b/src/storm/modelchecker/prctl/helper/HybridMdpPrctlHelper.cpp @@ -65,6 +65,49 @@ namespace storm { return result; } + template + void eliminateExtendedStatesFromExplicitRepresentation(std::pair, std::vector>& explicitRepresentation, std::vector& scheduler, storm::storage::BitVector const& properMaybeStates) { + // Eliminate superfluous entries from the scheduler. + uint64_t position = 0; + for (auto state : properMaybeStates) { + scheduler[position] = scheduler[state]; + position++; + } + scheduler.resize(properMaybeStates.getNumberOfSetBits()); + scheduler.shrink_to_fit(); + + // Treat the matrix. + explicitRepresentation.first = explicitRepresentation.first.getSubmatrix(true, properMaybeStates, properMaybeStates); + } + + template + void eliminateEndComponentsAndExtendedStatesUntilProbabilities(std::pair, std::vector>& explicitRepresentation, SolverRequirementsData& solverRequirementsData, storm::storage::BitVector const& targetStates) { + + // Get easier handles to the data. + auto& transitionMatrix = explicitRepresentation.first; + auto& oneStepProbabilities = explicitRepresentation.second; + + bool doDecomposition = !solverRequirementsData.properMaybeStates.empty(); + + storm::storage::MaximalEndComponentDecomposition endComponentDecomposition; + if (doDecomposition) { + // Then compute the states that are in MECs with zero reward. + endComponentDecomposition = storm::storage::MaximalEndComponentDecomposition(transitionMatrix, transitionMatrix.transpose(), solverRequirementsData.properMaybeStates); + } + + // Only do more work if there are actually end-components. + if (doDecomposition && !endComponentDecomposition.empty()) { + STORM_LOG_DEBUG("Eliminating " << endComponentDecomposition.size() << " EC(s)."); + std::vector subvector; + SparseMdpEndComponentInformation::eliminateEndComponents(endComponentDecomposition, transitionMatrix, solverRequirementsData.properMaybeStates, &targetStates, nullptr, nullptr, explicitRepresentation.first, &subvector, nullptr); + oneStepProbabilities = std::move(subvector); + } else { + STORM_LOG_DEBUG("Not eliminating ECs as there are none."); + eliminateExtendedStatesFromExplicitRepresentation(explicitRepresentation, solverRequirementsData.initialScheduler.get(), solverRequirementsData.properMaybeStates); + oneStepProbabilities = explicitRepresentation.first.getConstrainedRowGroupSumVector(solverRequirementsData.properMaybeStates, targetStates); + } + } + template std::unique_ptr HybridMdpPrctlHelper::computeUntilProbabilities(OptimizationDirection dir, storm::models::symbolic::NondeterministicModel const& model, storm::dd::Add const& transitionMatrix, storm::dd::Bdd const& phiStates, storm::dd::Bdd const& psiStates, bool qualitative, storm::solver::MinMaxLinearEquationSolverFactory const& linearEquationSolverFactory) { // We need to identify the states which have to be taken out of the matrix, i.e. all states that have @@ -91,7 +134,8 @@ namespace storm { storm::solver::MinMaxLinearEquationSolverRequirements requirements = linearEquationSolverFactory.getRequirements(storm::solver::EquationSystemType::UntilProbabilities, dir); storm::solver::MinMaxLinearEquationSolverRequirements clearedRequirements = requirements; SolverRequirementsData solverRequirementsData; - bool extendMaybeStates = true; + bool extendMaybeStates = false; + if (!clearedRequirements.empty()) { if (clearedRequirements.requiresNoEndComponents()) { STORM_LOG_DEBUG("Scheduling EC elimination, because the solver requires it."); @@ -116,31 +160,50 @@ namespace storm { // Create the ODD for the translation between symbolic and explicit storage. storm::dd::Odd odd = extendedMaybeStates.createOdd(); - // Create the matrix and the vector for the equation system. - storm::dd::Add extendedMaybeStatesAdd = extendedMaybeStates.template toAdd(); + // Convert the maybe states BDD to an ADD. + storm::dd::Add maybeStatesAdd = maybeStates.template toAdd(); // Start by cutting away all rows that do not belong to maybe states. Note that this leaves columns targeting // non-maybe states in the matrix. - storm::dd::Add submatrix = transitionMatrix * extendedMaybeStatesAdd; + storm::dd::Add submatrix = transitionMatrix * maybeStatesAdd; - // Then compute the vector that contains the one-step probabilities to a state with probability 1 for all - // maybe states. - storm::dd::Add prob1StatesAsColumn = statesWithProbability01.second.template toAdd(); - prob1StatesAsColumn = prob1StatesAsColumn.swapVariables(model.getRowColumnMetaVariablePairs()); - storm::dd::Add subvector = submatrix * prob1StatesAsColumn; - subvector = subvector.sumAbstract(model.getColumnVariables()); - - // Before cutting the non-maybe columns, we need to compute the sizes of the row groups. - std::vector rowGroupSizes = submatrix.notZero().existsAbstract(model.getColumnVariables()).template toAdd().sumAbstract(model.getNondeterminismVariables()).toVector(odd); - - // Finally cut away all columns targeting non-maybe states. - submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs()); - - // Translate the symbolic matrix/vector to their explicit representations and solve the equation system. - std::pair, std::vector> explicitRepresentation = submatrix.toMatrixVector(subvector, std::move(rowGroupSizes), model.getNondeterminismVariables(), odd, odd); - - if (requirements.requiresValidInitialScheduler()) { - solverRequirementsData.initialScheduler = computeValidInitialSchedulerForUntilProbabilities(explicitRepresentation.first, explicitRepresentation.second); + // If the maybe states were extended, we generate the explicit representation slightly differently. + std::pair, std::vector> explicitRepresentation; + if (extendMaybeStates) { + // Eliminate all transitions to non-extended-maybe states. + submatrix *= extendedMaybeStates.template toAdd().swapVariables(model.getRowColumnMetaVariablePairs()); + + // Only translate the matrix for now. + explicitRepresentation.first = submatrix.toMatrix(model.getNondeterminismVariables(), odd, odd); + + // Get all original maybe states in the extended matrix. + solverRequirementsData.properMaybeStates = maybeStates.toVector(odd); + + // Compute the target states within the set of extended maybe states. + storm::storage::BitVector targetStates = (extendedMaybeStates && statesWithProbability01.second).toVector(odd); + + // Eliminate the end components and remove the states that are not interesting (target or non-filter). + eliminateEndComponentsAndExtendedStatesUntilProbabilities(explicitRepresentation, solverRequirementsData, targetStates); + } else { + // Then compute the vector that contains the one-step probabilities to a state with probability 1 for all + // maybe states. + storm::dd::Add prob1StatesAsColumn = statesWithProbability01.second.template toAdd(); + prob1StatesAsColumn = prob1StatesAsColumn.swapVariables(model.getRowColumnMetaVariablePairs()); + storm::dd::Add subvector = submatrix * prob1StatesAsColumn; + subvector = subvector.sumAbstract(model.getColumnVariables()); + + // Before cutting the non-maybe columns, we need to compute the sizes of the row groups. + std::vector rowGroupSizes = submatrix.notZero().existsAbstract(model.getColumnVariables()).template toAdd().sumAbstract(model.getNondeterminismVariables()).toVector(odd); + + // Finally cut away all columns targeting non-maybe states. + submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs()); + + // Translate the symbolic matrix/vector to their explicit representations and solve the equation system. + explicitRepresentation = submatrix.toMatrixVector(subvector, std::move(rowGroupSizes), model.getNondeterminismVariables(), odd, odd); + + if (requirements.requiresValidInitialScheduler()) { + solverRequirementsData.initialScheduler = computeValidInitialSchedulerForUntilProbabilities(explicitRepresentation.first, explicitRepresentation.second); + } } // Create the solution vector. @@ -154,6 +217,18 @@ namespace storm { solver->setRequirementsChecked(); solver->solveEquations(dir, x, explicitRepresentation.second); + // If we included some target and non-filter states in the ODD, we need to expand the result from the solver. + if (requirements.requiresNoEndComponents() && solverRequirementsData.ecInformation) { + std::vector extendedVector(solverRequirementsData.properMaybeStates.getNumberOfSetBits()); + solverRequirementsData.ecInformation.get().setValues(extendedVector, solverRequirementsData.properMaybeStates, x); + x = std::move(extendedVector); + } + + // If we extended the maybe states, we create a new ODD containing only the propery maybe states. + if (extendMaybeStates) { + odd = maybeStates.createOdd(); + } + // Return a hybrid check result that stores the numerical values explicitly. return std::unique_ptr(new storm::modelchecker::HybridQuantitativeCheckResult(model.getReachableStates(), model.getReachableStates() && !maybeStates, statesWithProbability01.second.template toAdd(), maybeStates, odd, x)); } else { @@ -310,34 +385,6 @@ namespace storm { return result; } - template - void eliminateTargetStatesFromExplicitRepresentation(std::pair, std::vector>& explicitRepresentation, std::vector& scheduler, storm::storage::BitVector const& properMaybeStates) { - // Eliminate superfluous entries from the scheduler. - uint64_t position = 0; - for (auto state : properMaybeStates) { - scheduler[position] = scheduler[state]; - position++; - } - scheduler.resize(properMaybeStates.getNumberOfSetBits()); - scheduler.shrink_to_fit(); - - // Treat the matrix. - explicitRepresentation.first = explicitRepresentation.first.getSubmatrix(true, properMaybeStates, properMaybeStates); - } - - template - std::vector insertTargetStateValuesIntoExplicitRepresentation(std::vector const& x, storm::storage::BitVector const& properMaybeStates) { - std::vector expandedResult(properMaybeStates.size(), storm::utility::zero()); - - uint64_t position = 0; - for (auto state : properMaybeStates) { - expandedResult[state] = x[position]; - position++; - } - - return expandedResult; - } - template void eliminateEndComponentsAndTargetStatesReachabilityRewards(std::pair, std::vector>& explicitRepresentation, SolverRequirementsData& solverRequirementsData) { @@ -389,6 +436,7 @@ namespace storm { rewardVector = std::move(subvector); } else { STORM_LOG_DEBUG("Not eliminating ECs as there are none."); + eliminateExtendedStatesFromExplicitRepresentation(explicitRepresentation, solverRequirementsData.initialScheduler.get(), solverRequirementsData.properMaybeStates); } } @@ -487,7 +535,7 @@ namespace storm { // Since we needed the transitions to target states to be translated as well for the computation // of the scheduler, we have to get rid of them now. - eliminateTargetStatesFromExplicitRepresentation(explicitRepresentation, solverRequirementsData.initialScheduler.get(), solverRequirementsData.properMaybeStates); + eliminateExtendedStatesFromExplicitRepresentation(explicitRepresentation, solverRequirementsData.initialScheduler.get(), solverRequirementsData.properMaybeStates); } } @@ -505,20 +553,19 @@ namespace storm { solver->setLowerBound(storm::utility::zero()); solver->setRequirementsChecked(); solver->solveEquations(dir, x, explicitRepresentation.second); - - // If we included the target states in the ODD, we need to expand the result from the solver. + + // If we eliminated end components, we need to extend the solution vector. + if (requirements.requiresNoEndComponents() && solverRequirementsData.ecInformation) { + std::vector extendedVector(solverRequirementsData.properMaybeStates.getNumberOfSetBits()); + solverRequirementsData.ecInformation.get().setValues(extendedVector, solverRequirementsData.properMaybeStates, x); + x = std::move(extendedVector); + } + + // If we extended the maybe states, we create a new ODD that only contains proper maybe states. if (extendMaybeStates) { - if (requirements.requiresNoEndComponents()) { - if (solverRequirementsData.ecInformation) { - std::vector extendedVector(solverRequirementsData.properMaybeStates.size()); - solverRequirementsData.ecInformation.get().setValues(extendedVector, solverRequirementsData.properMaybeStates, x); - x = std::move(extendedVector); - } - } else { - x = insertTargetStateValuesIntoExplicitRepresentation(x, solverRequirementsData.properMaybeStates); - } + odd = maybeStates.createOdd(); } - + // Return a hybrid check result that stores the numerical values explicitly. return std::unique_ptr(new storm::modelchecker::HybridQuantitativeCheckResult(model.getReachableStates(), model.getReachableStates() && !maybeStates, infinityStates.ite(model.getManager().getConstant(storm::utility::infinity()), model.getManager().template getAddZero()), maybeStates, odd, x)); } else { diff --git a/src/storm/modelchecker/prctl/helper/SparseMdpEndComponentInformation.cpp b/src/storm/modelchecker/prctl/helper/SparseMdpEndComponentInformation.cpp index fb00c8243..f9d46ae94 100644 --- a/src/storm/modelchecker/prctl/helper/SparseMdpEndComponentInformation.cpp +++ b/src/storm/modelchecker/prctl/helper/SparseMdpEndComponentInformation.cpp @@ -161,6 +161,7 @@ namespace storm { if (stateActions.second.find(row) != stateActions.second.end()) { continue; } + // If the choices is not in the selected ones, drop it. if (selectedChoices && !selectedChoices->get(row)) { continue; @@ -198,7 +199,7 @@ namespace storm { } } - submatrix = builder.build(); + submatrix = builder.build(currentRow); return result; } diff --git a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp index 8f08a586c..e849f7b9e 100644 --- a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp +++ b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp @@ -489,6 +489,7 @@ namespace storm { this->createLowerBoundsVector(*lowerX); this->createUpperBoundsVector(this->auxiliaryRowGroupVector, this->A->getRowGroupCount()); std::vector* upperX = this->auxiliaryRowGroupVector.get(); + std::vector* tmp = nullptr; if (!useGaussSeidelMultiplication) { auxiliaryRowGroupVector2 = std::make_unique>(lowerX->size()); @@ -618,7 +619,8 @@ namespace storm { reportStatus(status, iterations); // We take the means of the lower and upper bound so we guarantee the desired precision. - storm::utility::vector::applyPointwise(*lowerX, *upperX, *lowerX, [] (ValueType const& a, ValueType const& b) { return (a + b) / storm::utility::convertNumber(2.0); }); + ValueType two = storm::utility::convertNumber(2.0); + storm::utility::vector::applyPointwise(*lowerX, *upperX, *lowerX, [&two] (ValueType const& a, ValueType const& b) -> ValueType { return (a + b) / two; }); // Since we shuffled the pointer around, we need to write the actual results to the input/output vector x. if (&x == tmp) { diff --git a/src/storm/storage/SparseMatrix.cpp b/src/storm/storage/SparseMatrix.cpp index b37472baf..928e190a4 100644 --- a/src/storm/storage/SparseMatrix.cpp +++ b/src/storm/storage/SparseMatrix.cpp @@ -231,16 +231,13 @@ namespace storm { rowIndications.push_back(currentEntryCount); } - // If there are no rows, we need to erase the start index of the current (non-existing) row. - if (rowCount == 0) { - rowIndications.pop_back(); - } - // We put a sentinel element at the last position of the row indices array. This eases iteration work, // as now the indices of row i are always between rowIndications[i] and rowIndications[i + 1], also for // the first and last row. - rowIndications.push_back(currentEntryCount); - STORM_LOG_ASSERT(rowCount == rowIndications.size() - 1, "Wrong sizes of vectors."); + if (rowCount > 0) { + rowIndications.push_back(currentEntryCount); + } + STORM_LOG_ASSERT(rowCount == rowIndications.size() - 1, "Wrong sizes of vectors: " << rowCount << " != " << (rowIndications.size() - 1) << "."); uint_fast64_t columnCount = hasEntries ? highestColumn + 1 : 0; if (initialColumnCountSet && forceInitialDimensions) { STORM_LOG_THROW(columnCount <= initialColumnCount, storm::exceptions::InvalidStateException, "Expected not more than " << initialColumnCount << " columns, but got " << columnCount << "."); @@ -1478,7 +1475,8 @@ namespace storm { template ValueType SparseMatrix::multiplyRowWithVector(index_type row, std::vector const& vector) const { ValueType result = storm::utility::zero(); - for(auto const& entry : this->getRow(row)){ + + for (auto const& entry : this->getRow(row)){ result += entry.getValue() * vector[entry.getColumn()]; } return result; diff --git a/src/storm/storage/dd/Add.cpp b/src/storm/storage/dd/Add.cpp index 0a1d25553..f01680610 100644 --- a/src/storm/storage/dd/Add.cpp +++ b/src/storm/storage/dd/Add.cpp @@ -826,7 +826,7 @@ namespace storm { rowIndications[i] = rowIndications[i - 1]; } rowIndications[0] = 0; - + return std::make_pair(storm::storage::SparseMatrix(columnOdd.getTotalOffset(), std::move(rowIndications), std::move(columnsAndValues), std::move(rowGroupIndices)), std::move(explicitVector)); }