diff --git a/CMakeLists.txt b/CMakeLists.txt index f338ee334..fd7d430ca 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -57,6 +57,7 @@ set(TBB_INSTALL_DIR "${PROJECT_SOURCE_DIR}/resources/3rdparty/tbb42_20140122_mer find_package(Boost REQUIRED) find_package(Doxygen REQUIRED) +find_package(Gurobi) find_package(TBB) find_package(Threads REQUIRED) find_package(GMP) @@ -69,7 +70,7 @@ else() endif() message(STATUS "StoRM - Building ${CMAKE_BUILD_TYPE} version.") -if ("${GUROBI_ROOT}" STREQUAL "") +if ("${GUROBI_ROOT}" STREQUAL "" AND NOT GUROBI_FOUND) set(ENABLE_GUROBI OFF) else() set(ENABLE_GUROBI ON) @@ -338,7 +339,10 @@ endif(ADDITIONAL_LINK_DIRS) ## ############################################################# if (ENABLE_GUROBI) - link_directories("${GUROBI_ROOT}/lib") + if (NOT GUROBI_FOUND) + message(FATAL_ERROR "Gurobi was requested, but not found!") + endif() + #link_directories("${GUROBI_ROOT}/lib") endif() if (ENABLE_Z3) link_directories("${Z3_ROOT}/bin") @@ -425,11 +429,11 @@ target_link_libraries(storm ltl2dstar) ## ############################################################# if (ENABLE_GUROBI) - message (STATUS "StoRM - Linking with Gurobi") - include_directories("${GUROBI_ROOT}/include") - target_link_libraries(storm "gurobi56") - target_link_libraries(storm-functional-tests "gurobi56") - target_link_libraries(storm-performance-tests "gurobi56") + message (STATUS "StoRM - Linking with Gurobi (include: ${GUROBI_INCLUDE_DIRS})") + include_directories(${GUROBI_INCLUDE_DIRS}) + target_link_libraries(storm ${GUROBI_LIBRARY}) + target_link_libraries(storm-functional-tests ${GUROBI_LIBRARY}) + target_link_libraries(storm-performance-tests ${GUROBI_LIBRARY}) endif(ENABLE_GUROBI) ############################################################# diff --git a/resources/cmake/FindGurobi.cmake b/resources/cmake/FindGurobi.cmake new file mode 100644 index 000000000..93de74af7 --- /dev/null +++ b/resources/cmake/FindGurobi.cmake @@ -0,0 +1,66 @@ +#### Taken from http://www.openflipper.org/svnrepo/CoMISo/trunk/CoMISo/cmake/FindGUROBI.cmake + + +# - Try to find GUROBI +# Once done this will define +# GUROBI_FOUND - System has Gurobi +# GUROBI_INCLUDE_DIRS - The Gurobi include directories +# GUROBI_LIBRARIES - The libraries needed to use Gurobi + +if (GUROBI_INCLUDE_DIR) + # in cache already + set(GUROBI_FOUND TRUE) + set(GUROBI_INCLUDE_DIRS "${GUROBI_INCLUDE_DIR}" ) + set(GUROBI_LIBRARIES "${GUROBI_LIBRARY};${GUROBI_CXX_LIBRARY}" ) +else (GUROBI_INCLUDE_DIR) + +find_path(GUROBI_INCLUDE_DIR + NAMES gurobi_c++.h + PATHS "$ENV{GUROBI_HOME}/include" + "/Library/gurobi502/mac64/include" + "C:\\libs\\gurobi502\\include" + "C:\\gurobi600\\win64\\include" + "${GUROBI_ROOT}/include" + ) + +find_library( GUROBI_LIBRARY + NAMES gurobi + gurobi45 + gurobi46 + gurobi50 + gurobi51 + gurobi52 + gurobi55 + gurobi56 + gurobi60 + PATHS "$ENV{GUROBI_HOME}/lib" + "/Library/gurobi502/mac64/lib" + "C:\\libs\\gurobi502\\lib" + "C:\\gurobi600\\win64\\lib" + "${GUROBI_ROOT}/lib" + ) + +find_library( GUROBI_CXX_LIBRARY + NAMES gurobi_c++ + PATHS "$ENV{GUROBI_HOME}/lib" + "/Library/gurobi502/mac64/lib" + "C:\\libs\\gurobi502\\lib" + "C:\\gurobi600\\win64\\lib" + "${GUROBI_ROOT}/lib" + ) + +set(GUROBI_INCLUDE_DIRS "${GUROBI_INCLUDE_DIR}" ) +set(GUROBI_LIBRARIES "${GUROBI_LIBRARY};${GUROBI_CXX_LIBRARY}" ) + +# use c++ headers as default +# set(GUROBI_COMPILER_FLAGS "-DIL_STD" CACHE STRING "Gurobi Compiler Flags") + +include(FindPackageHandleStandardArgs) +# handle the QUIETLY and REQUIRED arguments and set LIBCPLEX_FOUND to TRUE +# if all listed variables are TRUE +find_package_handle_standard_args(GUROBI DEFAULT_MSG + GUROBI_LIBRARY GUROBI_CXX_LIBRARY GUROBI_INCLUDE_DIR) + +mark_as_advanced(GUROBI_INCLUDE_DIR GUROBI_LIBRARY GUROBI_CXX_LIBRARY) + +endif(GUROBI_INCLUDE_DIR) \ No newline at end of file diff --git a/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp b/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp index e7ac8c425..a2af58046 100644 --- a/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp +++ b/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp @@ -59,8 +59,6 @@ namespace storm { } } else if (formula.isPropositionalFormula()) { return true; - } else { - std::cout << formula << " and type " << typeid(formula).name() << std::endl; } return false; } @@ -203,7 +201,7 @@ namespace storm { storm::storage::BitVector trueStates(model.getNumberOfStates(), true); // Do some sanity checks to establish some required properties. - STORM_LOG_THROW(storm::settings::sparseDtmcEliminationModelCheckerSettings().getEliminationMethod() == storm::settings::modules::SparseDtmcEliminationModelCheckerSettings::EliminationMethod::State, storm::exceptions::InvalidArgumentException, "Unsupported elimination method for conditional probabilities."); + STORM_LOG_WARN_COND(storm::settings::sparseDtmcEliminationModelCheckerSettings().getEliminationMethod() == storm::settings::modules::SparseDtmcEliminationModelCheckerSettings::EliminationMethod::State, "The chosen elimination method is not available for computing conditional probabilities. Falling back to regular state elimination."); STORM_LOG_THROW(model.getInitialStates().getNumberOfSetBits() == 1, storm::exceptions::IllegalArgumentException, "Input model is required to have exactly one initial state."); storm::storage::sparse::state_type initialState = *model.getInitialStates().begin(); @@ -244,7 +242,8 @@ namespace storm { // Determine the set of initial states of the sub-DTMC. storm::storage::BitVector newInitialStates = model.getInitialStates() % maybeStates; - + STORM_LOG_DEBUG("Found new initial states: " << newInitialStates << " (old: " << model.getInitialStates() << ")"); + // Create a dummy vector for the one-step probabilities. std::vector oneStepProbabilities(maybeStates.getNumberOfSetBits(), storm::utility::zero()); @@ -292,14 +291,20 @@ namespace storm { } STORM_LOG_INFO("Eliminated " << states.size() << " states." << std::endl); - // Eliminate the transitions going into the initial state. - eliminateState(flexibleMatrix, oneStepProbabilities, *newInitialStates.begin(), flexibleBackwardTransitions, missingStateRewards, false); + // Eliminate the transitions going into the initial state (if there are any). + if (!flexibleBackwardTransitions.getRow(*newInitialStates.begin()).empty()) { + eliminateState(flexibleMatrix, oneStepProbabilities, *newInitialStates.begin(), flexibleBackwardTransitions, missingStateRewards, false); + } // Now we need to basically eliminate all chains of not-psi states after phi states and chains of not-phi // states after psi states. for (auto const& trans1 : flexibleMatrix.getRow(*newInitialStates.begin())) { auto initialStateSuccessor = trans1.getColumn(); + + STORM_LOG_DEBUG("Exploring successor " << initialStateSuccessor << " of the initial state."); + if (phiStates.get(initialStateSuccessor)) { + STORM_LOG_DEBUG("Is a phi state."); // If the state is both a phi and a psi state, we do not need to eliminate chains. if (psiStates.get(initialStateSuccessor)) { @@ -319,8 +324,13 @@ namespace storm { for (auto const& element : currentRow) { // If any of the successors is a phi state, we eliminate it (wrt. all its phi predecessors). if (!psiStates.get(element.getColumn())) { - eliminateState(flexibleMatrix, oneStepProbabilities, element.getColumn(), flexibleBackwardTransitions, missingStateRewards, false, true, phiStates); - hasNonPsiSuccessor = true; + typename FlexibleSparseMatrix::row_type const& successorRow = flexibleMatrix.getRow(element.getColumn()); + // Eliminate the successor only if there possibly is a psi state reachable through it. + if (successorRow.size() > 1 || (!successorRow.empty() && successorRow.front().getColumn() != element.getColumn())) { + STORM_LOG_DEBUG("Found non-psi successor " << element.getColumn() << " that needs to be eliminated."); + eliminateState(flexibleMatrix, oneStepProbabilities, element.getColumn(), flexibleBackwardTransitions, missingStateRewards, false, true, phiStates); + hasNonPsiSuccessor = true; + } } } STORM_LOG_ASSERT(!flexibleMatrix.getRow(initialStateSuccessor).empty(), "(1) New transitions expected to be non-empty."); @@ -328,7 +338,8 @@ namespace storm { } } else { STORM_LOG_ASSERT(psiStates.get(initialStateSuccessor), "Expected psi state."); - + STORM_LOG_DEBUG("Is a psi state."); + // At this point, we know that the state satisfies psi and not phi. // This means, we must compute the probability to reach phi states, which in turn means that we need // to eliminate all chains of non-phi states between the current state and phi states. @@ -343,8 +354,12 @@ namespace storm { for (auto const& element : currentRow) { // If any of the successors is a psi state, we eliminate it (wrt. all its psi predecessors). if (!phiStates.get(element.getColumn())) { - eliminateState(flexibleMatrix, oneStepProbabilities, element.getColumn(), flexibleBackwardTransitions, missingStateRewards, false, true, psiStates); - hasNonPhiSuccessor = true; + typename FlexibleSparseMatrix::row_type const& successorRow = flexibleMatrix.getRow(element.getColumn()); + if (successorRow.size() > 1 || (!successorRow.empty() && successorRow.front().getColumn() != element.getColumn())) { + STORM_LOG_DEBUG("Found non-phi successor " << element.getColumn() << " that needs to be eliminated."); + eliminateState(flexibleMatrix, oneStepProbabilities, element.getColumn(), flexibleBackwardTransitions, missingStateRewards, false, true, psiStates); + hasNonPhiSuccessor = true; + } } } } @@ -364,8 +379,9 @@ namespace storm { } else { ValueType additiveTerm = storm::utility::zero(); for (auto const& trans2 : flexibleMatrix.getRow(initialStateSuccessor)) { - STORM_LOG_ASSERT(psiStates.get(trans2.getColumn()), "Expected " << trans2.getColumn() << " to be a psi state."); - additiveTerm += trans2.getValue(); + if (psiStates.get(trans2.getColumn())) { + additiveTerm += trans2.getValue(); + } } additiveTerm *= trans1.getValue(); numerator += additiveTerm; @@ -376,8 +392,9 @@ namespace storm { denominator += trans1.getValue(); ValueType additiveTerm = storm::utility::zero(); for (auto const& trans2 : flexibleMatrix.getRow(initialStateSuccessor)) { - STORM_LOG_ASSERT(phiStates.get(trans2.getColumn()), "Expected " << trans2.getColumn() << " to be a phi state."); - additiveTerm += trans2.getValue(); + if (phiStates.get(trans2.getColumn())) { + additiveTerm += trans2.getValue(); + } } numerator += trans1.getValue() * additiveTerm; } @@ -734,6 +751,11 @@ namespace storm { typename FlexibleSparseMatrix::row_type& currentStatePredecessors = backwardTransitions.getRow(state); std::size_t numberOfPredecessors = currentStatePredecessors.size(); std::size_t predecessorForwardTransitionCount = 0; + + // In case we have a constrained elimination, we need to keep track of the new predecessors. + typename FlexibleSparseMatrix::row_type newCurrentStatePredecessors; + + // Now go through the predecessors and eliminate the ones (satisfying the constraint if given). for (auto const& predecessorEntry : currentStatePredecessors) { uint_fast64_t predecessor = predecessorEntry.getColumn(); @@ -745,8 +767,11 @@ namespace storm { // Skip the state if the elimination is constrained, but the predecessor is not in the constraint. if (constrained && !predecessorConstraint.get(predecessor)) { + newCurrentStatePredecessors.emplace_back(predecessor, storm::utility::one()); + STORM_LOG_DEBUG("Not eliminating predecessor " << predecessor << ", because it does not fit the filter."); continue; } + STORM_LOG_DEBUG("Eliminating predecessor " << predecessor << "."); // First, find the probability with which the predecessor can move to the current state, because // the other probabilities need to be scaled with this factor. @@ -826,9 +851,11 @@ namespace storm { for (auto const& successorEntry : currentStateSuccessors) { typename FlexibleSparseMatrix::row_type& successorBackwardTransitions = backwardTransitions.getRow(successorEntry.getColumn()); - // Delete the current state as a predecessor of the successor state. - typename FlexibleSparseMatrix::row_type::iterator elimIt = std::find_if(successorBackwardTransitions.begin(), successorBackwardTransitions.end(), [&](storm::storage::MatrixEntry const& a) { return a.getColumn() == state; }); - if (elimIt != successorBackwardTransitions.end()) { + // Delete the current state as a predecessor of the successor state only if we are going to remove the + // current state's forward transitions. + if (removeForwardTransitions) { + typename FlexibleSparseMatrix::row_type::iterator elimIt = std::find_if(successorBackwardTransitions.begin(), successorBackwardTransitions.end(), [&](storm::storage::MatrixEntry const& a) { return a.getColumn() == state; }); + STORM_LOG_ASSERT(elimIt != successorBackwardTransitions.end(), "Expected a proper backward transition, but found none."); successorBackwardTransitions.erase(elimIt); } @@ -841,28 +868,49 @@ namespace storm { newPredecessors.reserve((last1 - first1) + (last2 - first2)); std::insert_iterator result(newPredecessors, newPredecessors.end()); - - for (; first1 != last1; ++result) { - if (first2 == last2) { - std::copy_if(first1, last1, result, [&] (storm::storage::MatrixEntry const& a) { return a.getColumn() != state; }); - break; - } - if (first2->getColumn() < first1->getColumn()) { - if (first2->getColumn() != state) { - *result = *first2; + if (!constrained) { + for (; first1 != last1; ++result) { + if (first2 == last2) { + std::copy(first1, last1, result); + break; } - ++first2; - } else { - if (first1->getColumn() != state) { + if (first2->getColumn() < first1->getColumn()) { + if (first2->getColumn() != state) { + *result = *first2; + } + ++first2; + } else { *result = *first1; + if (first1->getColumn() == first2->getColumn()) { + ++first2; + } + ++first1; } - if (first1->getColumn() == first2->getColumn()) { + } + std::copy_if(first2, last2, result, [&] (storm::storage::MatrixEntry const& a) { return a.getColumn() != state; }); + } else { + // If the elimination is constrained, we need to be more selective when we set the new predecessors + // of the successor state. + for (; first1 != last1; ++result) { + if (first2 == last2) { + std::copy(first1, last1, result); + break; + } + if (first2->getColumn() < first1->getColumn()) { + if (first2->getColumn() != state) { + *result = *first2; + } ++first2; + } else { + *result = *first1; + if (first1->getColumn() == first2->getColumn()) { + ++first2; + } + ++first1; } - ++first1; } + std::copy_if(first2, last2, result, [&] (storm::storage::MatrixEntry const& a) { return a.getColumn() != state && (!constrained || predecessorConstraint.get(a.getColumn())); }); } - std::copy_if(first2, last2, result, [&] (storm::storage::MatrixEntry const& a) { return a.getColumn() != state; }); // Now move the new predecessors in place. successorBackwardTransitions = std::move(newPredecessors); @@ -875,9 +923,10 @@ namespace storm { currentStateSuccessors.shrink_to_fit(); } if (!constrained) { - // FIXME: is this safe? If the elimination is constrained, we might have to repair the predecessor relation. currentStatePredecessors.clear(); currentStatePredecessors.shrink_to_fit(); + } else { + currentStatePredecessors = std::move(newCurrentStatePredecessors); } auto eliminationEnd = std::chrono::high_resolution_clock::now(); diff --git a/src/utility/cli.h b/src/utility/cli.h index ce4597096..c4ae034b1 100644 --- a/src/utility/cli.h +++ b/src/utility/cli.h @@ -387,8 +387,8 @@ namespace storm { std::unique_ptr result; if (model->getType() == storm::models::DTMC) { std::shared_ptr> dtmc = model->as>(); - storm::modelchecker::SparseDtmcPrctlModelChecker modelchecker(*dtmc); -// storm::modelchecker::SparseDtmcEliminationModelChecker modelchecker(*dtmc); +// storm::modelchecker::SparseDtmcPrctlModelChecker modelchecker(*dtmc); + storm::modelchecker::SparseDtmcEliminationModelChecker modelchecker(*dtmc); result = modelchecker.check(*formula.get()); } else if (model->getType() == storm::models::MDP) { std::shared_ptr> mdp = model->as>();