sp
/
tempest
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

We now avoid the renaming of 'Eigen' to 'StormEigen' as this is (hopefully) not needed anymore.

tempestpy_adaptions
Tim Quatmann 5 years ago
parent
739ef0b8b4
commit
9f1487392e
15 changed files with 106 additions and 104 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Unified View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 30
      src/storm/adapters/EigenAdapter.cpp
  2. 12
      src/storm/adapters/EigenAdapter.h
  3. 62
      src/storm/solver/EigenLinearEquationSolver.cpp
  4. 2
      src/storm/solver/EigenLinearEquationSolver.h
  5. 42
      src/storm/storage/geometry/NativePolytope.cpp
  6. 4
      src/storm/storage/geometry/NativePolytope.h
  7. 2
      src/storm/storage/geometry/nativepolytopeconversion/HyperplaneCollector.cpp
  8. 4
      src/storm/storage/geometry/nativepolytopeconversion/HyperplaneCollector.h
  9. 14
      src/storm/storage/geometry/nativepolytopeconversion/HyperplaneEnumeration.cpp
  10. 4
      src/storm/storage/geometry/nativepolytopeconversion/HyperplaneEnumeration.h
  11. 10
      src/storm/storage/geometry/nativepolytopeconversion/QuickHull.cpp
  12. 4
      src/storm/storage/geometry/nativepolytopeconversion/QuickHull.h
  13. 8
      src/storm/storage/geometry/nativepolytopeconversion/SubsetEnumerator.cpp
  14. 8
      src/storm/utility/eigen.h
  15. 4
      src/test/storm/modelchecker/csl/LraCtmcCslModelCheckerTest.cpp

30
src/storm/adapters/EigenAdapter.cpp
View File

@ -4,9 +4,9 @@ namespace storm {
namespace adapters { namespace adapters {
template<typename ValueType> template<typename ValueType>
std::unique_ptr<StormEigen::SparseMatrix<ValueType>> EigenAdapter::toEigenSparseMatrix(storm::storage::SparseMatrix<ValueType> const& matrix) {
std::unique_ptr<Eigen::SparseMatrix<ValueType>> EigenAdapter::toEigenSparseMatrix(storm::storage::SparseMatrix<ValueType> const& matrix) {
// Build a list of triplets and let Eigen care about the insertion. // Build a list of triplets and let Eigen care about the insertion.
std::vector<StormEigen::Triplet<ValueType>> triplets;
std::vector<Eigen::Triplet<ValueType>> triplets;
triplets.reserve(matrix.getNonzeroEntryCount()); triplets.reserve(matrix.getNonzeroEntryCount());
for (uint64_t row = 0; row < matrix.getRowCount(); ++row) { for (uint64_t row = 0; row < matrix.getRowCount(); ++row) {
@ -15,33 +15,33 @@ namespace storm {
} }
} }
std::unique_ptr<StormEigen::SparseMatrix<ValueType>> result = std::make_unique<StormEigen::SparseMatrix<ValueType>>(matrix.getRowCount(), matrix.getColumnCount());
std::unique_ptr<Eigen::SparseMatrix<ValueType>> result = std::make_unique<Eigen::SparseMatrix<ValueType>>(matrix.getRowCount(), matrix.getColumnCount());
result->setFromTriplets(triplets.begin(), triplets.end()); result->setFromTriplets(triplets.begin(), triplets.end());
return result; return result;
} }
template <typename ValueType> template <typename ValueType>
std::vector<ValueType> EigenAdapter::toStdVector(StormEigen::Matrix<ValueType, StormEigen::Dynamic, 1> const& v){
std::vector<ValueType> EigenAdapter::toStdVector(Eigen::Matrix<ValueType, Eigen::Dynamic, 1> const& v){
return std::vector<ValueType>(v.data(), v.data() + v.rows()); return std::vector<ValueType>(v.data(), v.data() + v.rows());
} }
template <typename ValueType> template <typename ValueType>
StormEigen::Matrix<ValueType, StormEigen::Dynamic, 1> EigenAdapter::toEigenVector(std::vector<ValueType> const& v){
return StormEigen::Matrix<ValueType, StormEigen::Dynamic, 1>::Map(v.data(), v.size());
Eigen::Matrix<ValueType, Eigen::Dynamic, 1> EigenAdapter::toEigenVector(std::vector<ValueType> const& v){
return Eigen::Matrix<ValueType, Eigen::Dynamic, 1>::Map(v.data(), v.size());
} }
template std::unique_ptr<StormEigen::SparseMatrix<double>> EigenAdapter::toEigenSparseMatrix(storm::storage::SparseMatrix<double> const& matrix);
template std::vector<double> EigenAdapter::toStdVector(StormEigen::Matrix<double, StormEigen::Dynamic, 1> const& v);
template StormEigen::Matrix<double, StormEigen::Dynamic, 1> EigenAdapter::toEigenVector(std::vector<double> const& v);
template std::unique_ptr<Eigen::SparseMatrix<double>> EigenAdapter::toEigenSparseMatrix(storm::storage::SparseMatrix<double> const& matrix);
template std::vector<double> EigenAdapter::toStdVector(Eigen::Matrix<double, Eigen::Dynamic, 1> const& v);
template Eigen::Matrix<double, Eigen::Dynamic, 1> EigenAdapter::toEigenVector(std::vector<double> const& v);
#ifdef STORM_HAVE_CARL #ifdef STORM_HAVE_CARL
template std::unique_ptr<StormEigen::SparseMatrix<storm::RationalNumber>> EigenAdapter::toEigenSparseMatrix(storm::storage::SparseMatrix<storm::RationalNumber> const& matrix);
template std::vector<storm::RationalNumber> EigenAdapter::toStdVector(StormEigen::Matrix<storm::RationalNumber, StormEigen::Dynamic, 1> const& v);
template StormEigen::Matrix<storm::RationalNumber, StormEigen::Dynamic, 1> EigenAdapter::toEigenVector(std::vector<storm::RationalNumber> const& v);
template std::unique_ptr<Eigen::SparseMatrix<storm::RationalNumber>> EigenAdapter::toEigenSparseMatrix(storm::storage::SparseMatrix<storm::RationalNumber> const& matrix);
template std::vector<storm::RationalNumber> EigenAdapter::toStdVector(Eigen::Matrix<storm::RationalNumber, Eigen::Dynamic, 1> const& v);
template Eigen::Matrix<storm::RationalNumber, Eigen::Dynamic, 1> EigenAdapter::toEigenVector(std::vector<storm::RationalNumber> const& v);
template std::unique_ptr<StormEigen::SparseMatrix<storm::RationalFunction>> EigenAdapter::toEigenSparseMatrix(storm::storage::SparseMatrix<storm::RationalFunction> const& matrix);
template std::vector<storm::RationalFunction> EigenAdapter::toStdVector(StormEigen::Matrix<storm::RationalFunction, StormEigen::Dynamic, 1> const& v);
template StormEigen::Matrix<storm::RationalFunction, StormEigen::Dynamic, 1> EigenAdapter::toEigenVector(std::vector<storm::RationalFunction> const& v);
template std::unique_ptr<Eigen::SparseMatrix<storm::RationalFunction>> EigenAdapter::toEigenSparseMatrix(storm::storage::SparseMatrix<storm::RationalFunction> const& matrix);
template std::vector<storm::RationalFunction> EigenAdapter::toStdVector(Eigen::Matrix<storm::RationalFunction, Eigen::Dynamic, 1> const& v);
template Eigen::Matrix<storm::RationalFunction, Eigen::Dynamic, 1> EigenAdapter::toEigenVector(std::vector<storm::RationalFunction> const& v);
#endif #endif
} }
} }

12
src/storm/adapters/EigenAdapter.h
View File

@ -18,13 +18,13 @@ namespace storm {
* @return A pointer to a row-major sparse matrix in gmm++ format. * @return A pointer to a row-major sparse matrix in gmm++ format.
*/ */
template<class ValueType> template<class ValueType>
static std::unique_ptr<StormEigen::SparseMatrix<ValueType>> toEigenSparseMatrix(storm::storage::SparseMatrix<ValueType> const& matrix);
static std::unique_ptr<Eigen::SparseMatrix<ValueType>> toEigenSparseMatrix(storm::storage::SparseMatrix<ValueType> const& matrix);
template <typename ValueType> template <typename ValueType>
static std::vector<ValueType> toStdVector(StormEigen::Matrix<ValueType, StormEigen::Dynamic, 1> const& v);
static std::vector<ValueType> toStdVector(Eigen::Matrix<ValueType, Eigen::Dynamic, 1> const& v);
template <typename ValueType> template <typename ValueType>
static StormEigen::Matrix<ValueType, StormEigen::Dynamic, 1> toEigenVector(std::vector<ValueType> const& v);
static Eigen::Matrix<ValueType, Eigen::Dynamic, 1> toEigenVector(std::vector<ValueType> const& v);
}; };
} }
@ -33,8 +33,8 @@ namespace storm {
namespace std { namespace std {
template<class ValueType> template<class ValueType>
struct hash<StormEigen::Matrix<ValueType, StormEigen::Dynamic, 1>> {
std::size_t operator()(StormEigen::Matrix<ValueType, StormEigen::Dynamic, 1> const &vector) const {
struct hash<Eigen::Matrix<ValueType, Eigen::Dynamic, 1>> {
std::size_t operator()(Eigen::Matrix<ValueType, Eigen::Dynamic, 1> const &vector) const {
size_t seed = 0; size_t seed = 0;
for (uint_fast64_t i = 0; i < static_cast<uint_fast64_t>(vector.rows()); ++i) { for (uint_fast64_t i = 0; i < static_cast<uint_fast64_t>(vector.rows()); ++i) {
carl::hash_add(seed, std::hash<ValueType>()(vector(i))); carl::hash_add(seed, std::hash<ValueType>()(vector(i)));
@ -44,7 +44,7 @@ namespace std {
}; };
} }
namespace StormEigen {
namespace Eigen {
template<> struct NumTraits<storm::RationalNumber> : GenericNumTraits<storm::RationalNumber> template<> struct NumTraits<storm::RationalNumber> : GenericNumTraits<storm::RationalNumber>
{ {
typedef storm::RationalNumber Real; typedef storm::RationalNumber Real;

62
src/storm/solver/EigenLinearEquationSolver.cpp
View File

@ -72,14 +72,14 @@ namespace storm {
STORM_LOG_INFO("Solving linear equation system (" << x.size() << " rows) with with rational numbers using LU factorization (Eigen library)."); STORM_LOG_INFO("Solving linear equation system (" << x.size() << " rows) with with rational numbers using LU factorization (Eigen library).");
// Map the input vectors to Eigen's format. // Map the input vectors to Eigen's format.
auto eigenX = StormEigen::Matrix<storm::RationalNumber, StormEigen::Dynamic, 1>::Map(x.data(), x.size());
auto eigenB = StormEigen::Matrix<storm::RationalNumber, StormEigen::Dynamic, 1>::Map(b.data(), b.size());
auto eigenX = Eigen::Matrix<storm::RationalNumber, Eigen::Dynamic, 1>::Map(x.data(), x.size());
auto eigenB = Eigen::Matrix<storm::RationalNumber, Eigen::Dynamic, 1>::Map(b.data(), b.size());
StormEigen::SparseLU<StormEigen::SparseMatrix<storm::RationalNumber>, StormEigen::COLAMDOrdering<int>> solver;
Eigen::SparseLU<Eigen::SparseMatrix<storm::RationalNumber>, Eigen::COLAMDOrdering<int>> solver;
solver.compute(*eigenA); solver.compute(*eigenA);
solver._solve_impl(eigenB, eigenX); solver._solve_impl(eigenB, eigenX);
return solver.info() == StormEigen::ComputationInfo::Success;
return solver.info() == Eigen::ComputationInfo::Success;
} }
// Specialization for storm::RationalFunction // Specialization for storm::RationalFunction
@ -90,13 +90,13 @@ namespace storm {
STORM_LOG_INFO("Solving linear equation system (" << x.size() << " rows) with rational functions using LU factorization (Eigen library)."); STORM_LOG_INFO("Solving linear equation system (" << x.size() << " rows) with rational functions using LU factorization (Eigen library).");
// Map the input vectors to Eigen's format. // Map the input vectors to Eigen's format.
auto eigenX = StormEigen::Matrix<storm::RationalFunction, StormEigen::Dynamic, 1>::Map(x.data(), x.size());
auto eigenB = StormEigen::Matrix<storm::RationalFunction, StormEigen::Dynamic, 1>::Map(b.data(), b.size());
auto eigenX = Eigen::Matrix<storm::RationalFunction, Eigen::Dynamic, 1>::Map(x.data(), x.size());
auto eigenB = Eigen::Matrix<storm::RationalFunction, Eigen::Dynamic, 1>::Map(b.data(), b.size());
StormEigen::SparseLU<StormEigen::SparseMatrix<storm::RationalFunction>, StormEigen::COLAMDOrdering<int>> solver;
Eigen::SparseLU<Eigen::SparseMatrix<storm::RationalFunction>, Eigen::COLAMDOrdering<int>> solver;
solver.compute(*eigenA); solver.compute(*eigenA);
solver._solve_impl(eigenB, eigenX); solver._solve_impl(eigenB, eigenX);
return solver.info() == StormEigen::ComputationInfo::Success;
return solver.info() == Eigen::ComputationInfo::Success;
} }
#endif #endif
@ -104,19 +104,19 @@ namespace storm {
template<typename ValueType> template<typename ValueType>
bool EigenLinearEquationSolver<ValueType>::internalSolveEquations(Environment const& env, std::vector<ValueType>& x, std::vector<ValueType> const& b) const { bool EigenLinearEquationSolver<ValueType>::internalSolveEquations(Environment const& env, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
// Map the input vectors to Eigen's format. // Map the input vectors to Eigen's format.
auto eigenX = StormEigen::Matrix<ValueType, StormEigen::Dynamic, 1>::Map(x.data(), x.size());
auto eigenB = StormEigen::Matrix<ValueType, StormEigen::Dynamic, 1>::Map(b.data(), b.size());
auto eigenX = Eigen::Matrix<ValueType, Eigen::Dynamic, 1>::Map(x.data(), x.size());
auto eigenB = Eigen::Matrix<ValueType, Eigen::Dynamic, 1>::Map(b.data(), b.size());
auto solutionMethod = getMethod(env, env.solver().isForceExact()); auto solutionMethod = getMethod(env, env.solver().isForceExact());
if (solutionMethod == EigenLinearEquationSolverMethod::SparseLU) { if (solutionMethod == EigenLinearEquationSolverMethod::SparseLU) {
STORM_LOG_INFO("Solving linear equation system (" << x.size() << " rows) with sparse LU factorization (Eigen library)."); STORM_LOG_INFO("Solving linear equation system (" << x.size() << " rows) with sparse LU factorization (Eigen library).");
StormEigen::SparseLU<StormEigen::SparseMatrix<ValueType>, StormEigen::COLAMDOrdering<int>> solver;
Eigen::SparseLU<Eigen::SparseMatrix<ValueType>, Eigen::COLAMDOrdering<int>> solver;
solver.compute(*this->eigenA); solver.compute(*this->eigenA);
solver._solve_impl(eigenB, eigenX); solver._solve_impl(eigenB, eigenX);
} else { } else {
bool converged = false; bool converged = false;
uint64_t numberOfIterations = 0; uint64_t numberOfIterations = 0;
StormEigen::Index maxIter = std::numeric_limits<StormEigen::Index>::max();
Eigen::Index maxIter = std::numeric_limits<Eigen::Index>::max();
if (env.solver().eigen().getMaximalNumberOfIterations() < static_cast<uint64_t>(maxIter)) { if (env.solver().eigen().getMaximalNumberOfIterations() < static_cast<uint64_t>(maxIter)) {
maxIter = env.solver().eigen().getMaximalNumberOfIterations(); maxIter = env.solver().eigen().getMaximalNumberOfIterations();
} }
@ -127,102 +127,102 @@ namespace storm {
if (preconditioner == EigenLinearEquationSolverPreconditioner::Ilu) { if (preconditioner == EigenLinearEquationSolverPreconditioner::Ilu) {
STORM_LOG_INFO("Solving linear equation system (" << x.size() << " rows) with BiCGSTAB with Ilu preconditioner (Eigen library)."); STORM_LOG_INFO("Solving linear equation system (" << x.size() << " rows) with BiCGSTAB with Ilu preconditioner (Eigen library).");
StormEigen::BiCGSTAB<StormEigen::SparseMatrix<ValueType>, StormEigen::IncompleteLUT<ValueType>> solver;
Eigen::BiCGSTAB<Eigen::SparseMatrix<ValueType>, Eigen::IncompleteLUT<ValueType>> solver;
solver.compute(*this->eigenA); solver.compute(*this->eigenA);
solver.setTolerance(precision); solver.setTolerance(precision);
solver.setMaxIterations(maxIter); solver.setMaxIterations(maxIter);
eigenX = solver.solveWithGuess(eigenB, eigenX); eigenX = solver.solveWithGuess(eigenB, eigenX);
converged = solver.info() == StormEigen::ComputationInfo::Success;
converged = solver.info() == Eigen::ComputationInfo::Success;
numberOfIterations = solver.iterations(); numberOfIterations = solver.iterations();
} else if (preconditioner == EigenLinearEquationSolverPreconditioner::Diagonal) { } else if (preconditioner == EigenLinearEquationSolverPreconditioner::Diagonal) {
STORM_LOG_INFO("Solving linear equation system (" << x.size() << " rows) with BiCGSTAB with Diagonal preconditioner (Eigen library)."); STORM_LOG_INFO("Solving linear equation system (" << x.size() << " rows) with BiCGSTAB with Diagonal preconditioner (Eigen library).");
StormEigen::BiCGSTAB<StormEigen::SparseMatrix<ValueType>, StormEigen::DiagonalPreconditioner<ValueType>> solver;
Eigen::BiCGSTAB<Eigen::SparseMatrix<ValueType>, Eigen::DiagonalPreconditioner<ValueType>> solver;
solver.setTolerance(precision); solver.setTolerance(precision);
solver.setMaxIterations(maxIter); solver.setMaxIterations(maxIter);
solver.compute(*this->eigenA); solver.compute(*this->eigenA);
eigenX = solver.solveWithGuess(eigenB, eigenX); eigenX = solver.solveWithGuess(eigenB, eigenX);
converged = solver.info() == StormEigen::ComputationInfo::Success;
converged = solver.info() == Eigen::ComputationInfo::Success;
numberOfIterations = solver.iterations(); numberOfIterations = solver.iterations();
} else { } else {
STORM_LOG_INFO("Solving linear equation system (" << x.size() << " rows) with BiCGSTAB with identity preconditioner (Eigen library)."); STORM_LOG_INFO("Solving linear equation system (" << x.size() << " rows) with BiCGSTAB with identity preconditioner (Eigen library).");
StormEigen::BiCGSTAB<StormEigen::SparseMatrix<ValueType>, StormEigen::IdentityPreconditioner> solver;
Eigen::BiCGSTAB<Eigen::SparseMatrix<ValueType>, Eigen::IdentityPreconditioner> solver;
solver.setTolerance(precision); solver.setTolerance(precision);
solver.setMaxIterations(maxIter); solver.setMaxIterations(maxIter);
solver.compute(*this->eigenA); solver.compute(*this->eigenA);
eigenX = solver.solveWithGuess(eigenB, eigenX); eigenX = solver.solveWithGuess(eigenB, eigenX);
numberOfIterations = solver.iterations(); numberOfIterations = solver.iterations();
converged = solver.info() == StormEigen::ComputationInfo::Success;
converged = solver.info() == Eigen::ComputationInfo::Success;
} }
} else if (solutionMethod == EigenLinearEquationSolverMethod::DGmres) { } else if (solutionMethod == EigenLinearEquationSolverMethod::DGmres) {
if (preconditioner == EigenLinearEquationSolverPreconditioner::Ilu) { if (preconditioner == EigenLinearEquationSolverPreconditioner::Ilu) {
STORM_LOG_INFO("Solving linear equation system (" << x.size() << " rows) with DGMRES with Ilu preconditioner (Eigen library)."); STORM_LOG_INFO("Solving linear equation system (" << x.size() << " rows) with DGMRES with Ilu preconditioner (Eigen library).");
StormEigen::DGMRES<StormEigen::SparseMatrix<ValueType>, StormEigen::IncompleteLUT<ValueType>> solver;
Eigen::DGMRES<Eigen::SparseMatrix<ValueType>, Eigen::IncompleteLUT<ValueType>> solver;
solver.setTolerance(precision); solver.setTolerance(precision);
solver.setMaxIterations(maxIter); solver.setMaxIterations(maxIter);
solver.set_restart(restartThreshold); solver.set_restart(restartThreshold);
solver.compute(*this->eigenA); solver.compute(*this->eigenA);
eigenX = solver.solveWithGuess(eigenB, eigenX); eigenX = solver.solveWithGuess(eigenB, eigenX);
converged = solver.info() == StormEigen::ComputationInfo::Success;
converged = solver.info() == Eigen::ComputationInfo::Success;
numberOfIterations = solver.iterations(); numberOfIterations = solver.iterations();
} else if (preconditioner == EigenLinearEquationSolverPreconditioner::Diagonal) { } else if (preconditioner == EigenLinearEquationSolverPreconditioner::Diagonal) {
STORM_LOG_INFO("Solving linear equation system (" << x.size() << " rows) with DGMRES with Diagonal preconditioner (Eigen library)."); STORM_LOG_INFO("Solving linear equation system (" << x.size() << " rows) with DGMRES with Diagonal preconditioner (Eigen library).");
StormEigen::DGMRES<StormEigen::SparseMatrix<ValueType>, StormEigen::DiagonalPreconditioner<ValueType>> solver;
Eigen::DGMRES<Eigen::SparseMatrix<ValueType>, Eigen::DiagonalPreconditioner<ValueType>> solver;
solver.setTolerance(precision); solver.setTolerance(precision);
solver.setMaxIterations(maxIter); solver.setMaxIterations(maxIter);
solver.set_restart(restartThreshold); solver.set_restart(restartThreshold);
solver.compute(*this->eigenA); solver.compute(*this->eigenA);
eigenX = solver.solveWithGuess(eigenB, eigenX); eigenX = solver.solveWithGuess(eigenB, eigenX);
converged = solver.info() == StormEigen::ComputationInfo::Success;
converged = solver.info() == Eigen::ComputationInfo::Success;
numberOfIterations = solver.iterations(); numberOfIterations = solver.iterations();
} else { } else {
STORM_LOG_INFO("Solving linear equation system (" << x.size() << " rows) with DGMRES with identity preconditioner (Eigen library)."); STORM_LOG_INFO("Solving linear equation system (" << x.size() << " rows) with DGMRES with identity preconditioner (Eigen library).");
StormEigen::DGMRES<StormEigen::SparseMatrix<ValueType>, StormEigen::IdentityPreconditioner> solver;
Eigen::DGMRES<Eigen::SparseMatrix<ValueType>, Eigen::IdentityPreconditioner> solver;
solver.setTolerance(precision); solver.setTolerance(precision);
solver.setMaxIterations(maxIter); solver.setMaxIterations(maxIter);
solver.set_restart(restartThreshold); solver.set_restart(restartThreshold);
solver.compute(*this->eigenA); solver.compute(*this->eigenA);
eigenX = solver.solveWithGuess(eigenB, eigenX); eigenX = solver.solveWithGuess(eigenB, eigenX);
converged = solver.info() == StormEigen::ComputationInfo::Success;
converged = solver.info() == Eigen::ComputationInfo::Success;
numberOfIterations = solver.iterations(); numberOfIterations = solver.iterations();
} }
} else if (solutionMethod == EigenLinearEquationSolverMethod::Gmres) { } else if (solutionMethod == EigenLinearEquationSolverMethod::Gmres) {
if (preconditioner == EigenLinearEquationSolverPreconditioner::Ilu) { if (preconditioner == EigenLinearEquationSolverPreconditioner::Ilu) {
STORM_LOG_INFO("Solving linear equation system (" << x.size() << " rows) with GMRES with Ilu preconditioner (Eigen library)."); STORM_LOG_INFO("Solving linear equation system (" << x.size() << " rows) with GMRES with Ilu preconditioner (Eigen library).");
StormEigen::GMRES<StormEigen::SparseMatrix<ValueType>, StormEigen::IncompleteLUT<ValueType>> solver;
Eigen::GMRES<Eigen::SparseMatrix<ValueType>, Eigen::IncompleteLUT<ValueType>> solver;
solver.setTolerance(precision); solver.setTolerance(precision);
solver.setMaxIterations(maxIter); solver.setMaxIterations(maxIter);
solver.set_restart(restartThreshold); solver.set_restart(restartThreshold);
solver.compute(*this->eigenA); solver.compute(*this->eigenA);
eigenX = solver.solveWithGuess(eigenB, eigenX); eigenX = solver.solveWithGuess(eigenB, eigenX);
converged = solver.info() == StormEigen::ComputationInfo::Success;
converged = solver.info() == Eigen::ComputationInfo::Success;
numberOfIterations = solver.iterations(); numberOfIterations = solver.iterations();
} else if (preconditioner == EigenLinearEquationSolverPreconditioner::Diagonal) { } else if (preconditioner == EigenLinearEquationSolverPreconditioner::Diagonal) {
STORM_LOG_INFO("Solving linear equation system (" << x.size() << " rows) with GMRES with Diagonal preconditioner (Eigen library)."); STORM_LOG_INFO("Solving linear equation system (" << x.size() << " rows) with GMRES with Diagonal preconditioner (Eigen library).");
StormEigen::GMRES<StormEigen::SparseMatrix<ValueType>, StormEigen::DiagonalPreconditioner<ValueType>> solver;
Eigen::GMRES<Eigen::SparseMatrix<ValueType>, Eigen::DiagonalPreconditioner<ValueType>> solver;
solver.setTolerance(precision); solver.setTolerance(precision);
solver.setMaxIterations(maxIter); solver.setMaxIterations(maxIter);
solver.set_restart(restartThreshold); solver.set_restart(restartThreshold);
solver.compute(*this->eigenA); solver.compute(*this->eigenA);
eigenX = solver.solveWithGuess(eigenB, eigenX); eigenX = solver.solveWithGuess(eigenB, eigenX);
converged = solver.info() == StormEigen::ComputationInfo::Success;
converged = solver.info() == Eigen::ComputationInfo::Success;
numberOfIterations = solver.iterations(); numberOfIterations = solver.iterations();
} else { } else {
STORM_LOG_INFO("Solving linear equation system (" << x.size() << " rows) with GMRES with identity preconditioner (Eigen library)."); STORM_LOG_INFO("Solving linear equation system (" << x.size() << " rows) with GMRES with identity preconditioner (Eigen library).");
StormEigen::GMRES<StormEigen::SparseMatrix<ValueType>, StormEigen::IdentityPreconditioner> solver;
Eigen::GMRES<Eigen::SparseMatrix<ValueType>, Eigen::IdentityPreconditioner> solver;
solver.setTolerance(precision); solver.setTolerance(precision);
solver.setMaxIterations(maxIter); solver.setMaxIterations(maxIter);
solver.set_restart(restartThreshold); solver.set_restart(restartThreshold);
solver.compute(*this->eigenA); solver.compute(*this->eigenA);
eigenX = solver.solveWithGuess(eigenB, eigenX); eigenX = solver.solveWithGuess(eigenB, eigenX);
converged = solver.info() == StormEigen::ComputationInfo::Success;
converged = solver.info() == Eigen::ComputationInfo::Success;
numberOfIterations = solver.iterations(); numberOfIterations = solver.iterations();
} }
} }

2
src/storm/solver/EigenLinearEquationSolver.h
View File

@ -32,7 +32,7 @@ namespace storm {
virtual uint64_t getMatrixColumnCount() const override; virtual uint64_t getMatrixColumnCount() const override;
// The (eigen) matrix associated with this equation solver. // The (eigen) matrix associated with this equation solver.
std::unique_ptr<StormEigen::SparseMatrix<ValueType>> eigenA;
std::unique_ptr<Eigen::SparseMatrix<ValueType>> eigenA;
}; };

42
src/storm/storage/geometry/NativePolytope.cpp
View File

@ -22,12 +22,12 @@ namespace storm {
// The polytope is universal // The polytope is universal
emptyStatus = EmptyStatus::Nonempty; emptyStatus = EmptyStatus::Nonempty;
} else { } else {
StormEigen::Index maxCol = halfspaces.front().normalVector().size();
StormEigen::Index maxRow = halfspaces.size();
Eigen::Index maxCol = halfspaces.front().normalVector().size();
Eigen::Index maxRow = halfspaces.size();
A = EigenMatrix(maxRow, maxCol); A = EigenMatrix(maxRow, maxCol);
b = EigenVector(maxRow); b = EigenVector(maxRow);
for (StormEigen::Index row = 0; row < A.rows(); ++row ){
assert((StormEigen::Index) halfspaces[row].normalVector().size() == maxCol);
for (Eigen::Index row = 0; row < A.rows(); ++row ){
assert((Eigen::Index) halfspaces[row].normalVector().size() == maxCol);
b(row) = halfspaces[row].offset(); b(row) = halfspaces[row].offset();
A.row(row) = storm::adapters::EigenAdapter::toEigenVector(halfspaces[row].normalVector()); A.row(row) = storm::adapters::EigenAdapter::toEigenVector(halfspaces[row].normalVector());
} }
@ -107,7 +107,7 @@ namespace storm {
std::vector<Halfspace<ValueType>> result; std::vector<Halfspace<ValueType>> result;
result.reserve(A.rows()); result.reserve(A.rows());
for (StormEigen::Index row=0; row < A.rows(); ++row){
for (Eigen::Index row=0; row < A.rows(); ++row){
result.emplace_back(storm::adapters::EigenAdapter::toStdVector(EigenVector(A.row(row))), b(row)); result.emplace_back(storm::adapters::EigenAdapter::toStdVector(EigenVector(A.row(row))), b(row));
} }
return result; return result;
@ -145,7 +145,7 @@ namespace storm {
template <typename ValueType> template <typename ValueType>
bool NativePolytope<ValueType>::contains(Point const& point) const{ bool NativePolytope<ValueType>::contains(Point const& point) const{
EigenVector x = storm::adapters::EigenAdapter::toEigenVector(point); EigenVector x = storm::adapters::EigenAdapter::toEigenVector(point);
for (StormEigen::Index row=0; row < A.rows(); ++row){
for (Eigen::Index row=0; row < A.rows(); ++row){
if ((A.row(row) * x)(0) > b(row)){ if ((A.row(row) * x)(0) > b(row)){
return false; return false;
} }
@ -252,7 +252,7 @@ namespace storm {
resultConstraints.reserve(A.rows() + nativeRhs.A.rows()); resultConstraints.reserve(A.rows() + nativeRhs.A.rows());
// evaluation of rhs in directions of lhs // evaluation of rhs in directions of lhs
for (StormEigen::Index i = 0; i < A.rows(); ++i) {
for (Eigen::Index i = 0; i < A.rows(); ++i) {
auto optimizationRes = nativeRhs.optimize(A.row(i)); auto optimizationRes = nativeRhs.optimize(A.row(i));
if ( optimizationRes.second ) { if ( optimizationRes.second ) {
resultConstraints.emplace_back(A.row(i), b(i) + (A.row(i) * optimizationRes.first)(0)); resultConstraints.emplace_back(A.row(i), b(i) + (A.row(i) * optimizationRes.first)(0));
@ -261,7 +261,7 @@ namespace storm {
} }
// evaluation of lhs in directions of rhs // evaluation of lhs in directions of rhs
for (StormEigen::Index i = 0; i < nativeRhs.A.rows(); ++i) {
for (Eigen::Index i = 0; i < nativeRhs.A.rows(); ++i) {
auto optimizationRes = optimize(nativeRhs.A.row(i)); auto optimizationRes = optimize(nativeRhs.A.row(i));
if ( optimizationRes.second ) { if ( optimizationRes.second ) {
resultConstraints.emplace_back(nativeRhs.A.row(i), nativeRhs.b(i) + (nativeRhs.A.row(i) * optimizationRes.first)(0)); resultConstraints.emplace_back(nativeRhs.A.row(i), nativeRhs.b(i) + (nativeRhs.A.row(i) * optimizationRes.first)(0));
@ -274,7 +274,7 @@ namespace storm {
} else { } else {
EigenMatrix newA(resultConstraints.size(), resultConstraints.front().first.rows()); EigenMatrix newA(resultConstraints.size(), resultConstraints.front().first.rows());
EigenVector newb(resultConstraints.size()); EigenVector newb(resultConstraints.size());
for (StormEigen::Index i = 0; i < newA.rows(); ++i) {
for (Eigen::Index i = 0; i < newA.rows(); ++i) {
newA.row(i) = resultConstraints[i].first; newA.row(i) = resultConstraints[i].first;
newb(i) = resultConstraints[i].second; newb(i) = resultConstraints[i].second;
} }
@ -285,15 +285,15 @@ namespace storm {
template <typename ValueType> template <typename ValueType>
std::shared_ptr<Polytope<ValueType>> NativePolytope<ValueType>::affineTransformation(std::vector<Point> const& matrix, Point const& vector) const { std::shared_ptr<Polytope<ValueType>> NativePolytope<ValueType>::affineTransformation(std::vector<Point> const& matrix, Point const& vector) const {
STORM_LOG_THROW(!matrix.empty(), storm::exceptions::InvalidArgumentException, "Invoked affine transformation with a matrix without rows."); STORM_LOG_THROW(!matrix.empty(), storm::exceptions::InvalidArgumentException, "Invoked affine transformation with a matrix without rows.");
StormEigen::Index rows = matrix.size();
StormEigen::Index columns = matrix.front().size();
Eigen::Index rows = matrix.size();
Eigen::Index columns = matrix.front().size();
EigenMatrix eigenMatrix(rows, columns); EigenMatrix eigenMatrix(rows, columns);
for (StormEigen::Index row = 0; row < rows; ++row) {
for (Eigen::Index row = 0; row < rows; ++row) {
eigenMatrix.row(row) = storm::adapters::EigenAdapter::toEigenVector(matrix[row]); eigenMatrix.row(row) = storm::adapters::EigenAdapter::toEigenVector(matrix[row]);
} }
EigenVector eigenVector = storm::adapters::EigenAdapter::toEigenVector(vector); EigenVector eigenVector = storm::adapters::EigenAdapter::toEigenVector(vector);
StormEigen::FullPivLU<EigenMatrix> luMatrix( eigenMatrix );
Eigen::FullPivLU<EigenMatrix> luMatrix( eigenMatrix );
STORM_LOG_THROW(luMatrix.isInvertible(), storm::exceptions::NotImplementedException, "Affine Transformation of native polytope only implemented if the transformation matrix is invertable"); STORM_LOG_THROW(luMatrix.isInvertible(), storm::exceptions::NotImplementedException, "Affine Transformation of native polytope only implemented if the transformation matrix is invertable");
if (isUniversal()) { if (isUniversal()) {
return std::make_shared<NativePolytope<ValueType>>(std::vector<Halfspace<ValueType>>()); return std::make_shared<NativePolytope<ValueType>>(std::vector<Halfspace<ValueType>>());
@ -312,7 +312,7 @@ namespace storm {
storm::solver::Z3LpSolver<ValueType> solver(storm::solver::OptimizationDirection::Maximize); storm::solver::Z3LpSolver<ValueType> solver(storm::solver::OptimizationDirection::Maximize);
std::vector<storm::expressions::Variable> variables; std::vector<storm::expressions::Variable> variables;
variables.reserve(A.cols()); variables.reserve(A.cols());
for (StormEigen::Index i = 0; i < A.cols(); ++i) {
for (Eigen::Index i = 0; i < A.cols(); ++i) {
variables.push_back(solver.addUnboundedContinuousVariable("x" + std::to_string(i), direction[i])); variables.push_back(solver.addUnboundedContinuousVariable("x" + std::to_string(i), direction[i]));
} }
std::vector<storm::expressions::Expression> constraints = getConstraints(solver.getManager(), variables); std::vector<storm::expressions::Expression> constraints = getConstraints(solver.getManager(), variables);
@ -343,7 +343,7 @@ namespace storm {
storm::solver::Z3LpSolver<ValueType> solver(storm::solver::OptimizationDirection::Maximize); storm::solver::Z3LpSolver<ValueType> solver(storm::solver::OptimizationDirection::Maximize);
std::vector<storm::expressions::Variable> variables; std::vector<storm::expressions::Variable> variables;
variables.reserve(A.cols()); variables.reserve(A.cols());
for (StormEigen::Index i = 0; i < A.cols(); ++i) {
for (Eigen::Index i = 0; i < A.cols(); ++i) {
variables.push_back(solver.addUnboundedContinuousVariable("x" + std::to_string(i), static_cast<ValueType>(direction(i)))); variables.push_back(solver.addUnboundedContinuousVariable("x" + std::to_string(i), static_cast<ValueType>(direction(i))));
} }
std::vector<storm::expressions::Expression> constraints = getConstraints(solver.getManager(), variables); std::vector<storm::expressions::Expression> constraints = getConstraints(solver.getManager(), variables);
@ -354,7 +354,7 @@ namespace storm {
solver.optimize(); solver.optimize();
if (solver.isOptimal()) { if (solver.isOptimal()) {
auto result = std::make_pair(EigenVector(A.cols()), true); auto result = std::make_pair(EigenVector(A.cols()), true);
for (StormEigen::Index i = 0; i < A.cols(); ++i) {
for (Eigen::Index i = 0; i < A.cols(); ++i) {
result.first(i) = solver.getContinuousValue(variables[i]); result.first(i) = solver.getContinuousValue(variables[i]);
} }
return result; return result;
@ -379,7 +379,7 @@ namespace storm {
std::vector<storm::expressions::Variable> NativePolytope<ValueType>::declareVariables(storm::expressions::ExpressionManager& manager, std::string const& namePrefix) const { std::vector<storm::expressions::Variable> NativePolytope<ValueType>::declareVariables(storm::expressions::ExpressionManager& manager, std::string const& namePrefix) const {
std::vector<storm::expressions::Variable> result; std::vector<storm::expressions::Variable> result;
result.reserve(A.cols()); result.reserve(A.cols());
for (StormEigen::Index col=0; col < A.cols(); ++col){
for (Eigen::Index col=0; col < A.cols(); ++col){
result.push_back(manager.declareVariable(namePrefix + std::to_string(col), manager.getRationalType())); result.push_back(manager.declareVariable(namePrefix + std::to_string(col), manager.getRationalType()));
} }
return result; return result;
@ -388,9 +388,9 @@ namespace storm {
template <typename ValueType> template <typename ValueType>
std::vector<storm::expressions::Expression> NativePolytope<ValueType>::getConstraints(storm::expressions::ExpressionManager const& manager, std::vector<storm::expressions::Variable> const& variables) const { std::vector<storm::expressions::Expression> NativePolytope<ValueType>::getConstraints(storm::expressions::ExpressionManager const& manager, std::vector<storm::expressions::Variable> const& variables) const {
std::vector<storm::expressions::Expression> result; std::vector<storm::expressions::Expression> result;
for (StormEigen::Index row = 0; row < A.rows(); ++row) {
for (Eigen::Index row = 0; row < A.rows(); ++row) {
storm::expressions::Expression lhs = manager.rational(A(row,0)) * variables[0].getExpression(); storm::expressions::Expression lhs = manager.rational(A(row,0)) * variables[0].getExpression();
for (StormEigen::Index col=1; col < A.cols(); ++col) {
for (Eigen::Index col=1; col < A.cols(); ++col) {
lhs = lhs + manager.rational(A(row,col)) * variables[col].getExpression(); lhs = lhs + manager.rational(A(row,col)) * variables[col].getExpression();
} }
result.push_back(lhs <= manager.rational(b(row))); result.push_back(lhs <= manager.rational(b(row)));
@ -411,9 +411,9 @@ namespace storm {
solver->add(constraint); solver->add(constraint);
} }
storm::storage::BitVector keptConstraints(A.rows(), false); storm::storage::BitVector keptConstraints(A.rows(), false);
for (StormEigen::Index row = 0; row < A.rows(); ++row) {
for (Eigen::Index row = 0; row < A.rows(); ++row) {
storm::expressions::Expression lhs = manager->rational(A(row,0)) * variables[0].getExpression(); storm::expressions::Expression lhs = manager->rational(A(row,0)) * variables[0].getExpression();
for (StormEigen::Index col=1; col < A.cols(); ++col) {
for (Eigen::Index col=1; col < A.cols(); ++col) {
lhs = lhs + manager->rational(A(row,col)) * variables[col].getExpression(); lhs = lhs + manager->rational(A(row,col)) * variables[col].getExpression();
} }
solver->push(); solver->push();

4
src/storm/storage/geometry/NativePolytope.h
View File

@ -14,8 +14,8 @@ namespace storm {
class NativePolytope : public Polytope<ValueType> { class NativePolytope : public Polytope<ValueType> {
public: public:
typedef StormEigen::Matrix<ValueType, StormEigen::Dynamic, StormEigen::Dynamic> EigenMatrix;
typedef StormEigen::Matrix<ValueType, StormEigen::Dynamic, 1> EigenVector;
typedef Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic> EigenMatrix;
typedef Eigen::Matrix<ValueType, Eigen::Dynamic, 1> EigenVector;
enum class EmptyStatus{ enum class EmptyStatus{
Unknown, //It is unknown whether the polytope is empty or not Unknown, //It is unknown whether the polytope is empty or not

2
src/storm/storage/geometry/nativepolytopeconversion/HyperplaneCollector.cpp
View File

@ -15,7 +15,7 @@ namespace storm {
template< typename ValueType> template< typename ValueType>
bool HyperplaneCollector<ValueType>::insert(EigenVector && normal, ValueType && offset, std::vector<uint_fast64_t>const* indexList) { bool HyperplaneCollector<ValueType>::insert(EigenVector && normal, ValueType && offset, std::vector<uint_fast64_t>const* indexList) {
//Normalize //Normalize
ValueType infinityNorm = normal.template lpNorm<StormEigen::Infinity>();
ValueType infinityNorm = normal.template lpNorm<Eigen::Infinity>();
if(infinityNorm != (ValueType)0 ){ if(infinityNorm != (ValueType)0 ){
normal /= infinityNorm; normal /= infinityNorm;
offset /= infinityNorm; offset /= infinityNorm;

4
src/storm/storage/geometry/nativepolytopeconversion/HyperplaneCollector.h
View File

@ -17,8 +17,8 @@ namespace storm {
class HyperplaneCollector { class HyperplaneCollector {
public: public:
typedef StormEigen::Matrix<ValueType, StormEigen::Dynamic, StormEigen::Dynamic> EigenMatrix;
typedef StormEigen::Matrix<ValueType, StormEigen::Dynamic, 1> EigenVector;
typedef Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic> EigenMatrix;
typedef Eigen::Matrix<ValueType, Eigen::Dynamic, 1> EigenVector;
HyperplaneCollector() = default; HyperplaneCollector() = default;
HyperplaneCollector(const HyperplaneCollector& orig) = default; HyperplaneCollector(const HyperplaneCollector& orig) = default;

14
src/storm/storage/geometry/nativepolytopeconversion/HyperplaneEnumeration.cpp
View File

@ -11,7 +11,7 @@ namespace storm {
template<typename ValueType> template<typename ValueType>
void HyperplaneEnumeration<ValueType>::generateVerticesFromConstraints(EigenMatrix const& constraintMatrix, EigenVector const& constraintVector, bool generateRelevantHyperplanesAndVertexSets) { void HyperplaneEnumeration<ValueType>::generateVerticesFromConstraints(EigenMatrix const& constraintMatrix, EigenVector const& constraintVector, bool generateRelevantHyperplanesAndVertexSets) {
STORM_LOG_DEBUG("Invoked Hyperplane enumeration with " << constraintMatrix.rows() << " constraints."); STORM_LOG_DEBUG("Invoked Hyperplane enumeration with " << constraintMatrix.rows() << " constraints.");
StormEigen::Index dimension = constraintMatrix.cols();
Eigen::Index dimension = constraintMatrix.cols();
if (dimension == 0) { if (dimension == 0) {
// No halfspaces means no vertices // No halfspaces means no vertices
resultVertices.clear(); resultVertices.clear();
@ -28,14 +28,14 @@ namespace storm {
EigenMatrix subMatrix(dimension, dimension); EigenMatrix subMatrix(dimension, dimension);
EigenVector subVector(dimension); EigenVector subVector(dimension);
for (StormEigen::Index i = 0; i < dimension; ++i){
for (Eigen::Index i = 0; i < dimension; ++i){
subMatrix.row(i) = constraintMatrix.row(subset[i]); subMatrix.row(i) = constraintMatrix.row(subset[i]);
subVector(i) = constraintVector(subset[i]); subVector(i) = constraintVector(subset[i]);
} }
EigenVector point = subMatrix.fullPivLu().solve(subVector); EigenVector point = subMatrix.fullPivLu().solve(subVector);
bool pointContained = true; bool pointContained = true;
for (StormEigen::Index row=0; row < constraintMatrix.rows(); ++row){
for (Eigen::Index row=0; row < constraintMatrix.rows(); ++row){
if ((constraintMatrix.row(row) * point)(0) > constraintVector(row)){ if ((constraintMatrix.row(row) * point)(0) > constraintVector(row)){
pointContained = false; pointContained = false;
break; break;
@ -53,7 +53,7 @@ namespace storm {
if (generateRelevantHyperplanesAndVertexSets){ if (generateRelevantHyperplanesAndVertexSets){
//For each hyperplane, get the number of (unique) vertices that lie on it. //For each hyperplane, get the number of (unique) vertices that lie on it.
std::vector<StormEigen::Index> verticesOnHyperplaneCounter(constraintMatrix.rows(), 0);
std::vector<Eigen::Index> verticesOnHyperplaneCounter(constraintMatrix.rows(), 0);
for (auto const& mapEntry : vertexCollector){ for (auto const& mapEntry : vertexCollector){
for (auto const& hyperplaneIndex : mapEntry.second){ for (auto const& hyperplaneIndex : mapEntry.second){
++verticesOnHyperplaneCounter[hyperplaneIndex]; ++verticesOnHyperplaneCounter[hyperplaneIndex];
@ -91,7 +91,7 @@ namespace storm {
for (auto const& mapEntry : vertexCollector){ for (auto const& mapEntry : vertexCollector){
for (auto const& oldHyperplaneIndex : mapEntry.second){ for (auto const& oldHyperplaneIndex : mapEntry.second){
//ignore the hyperplanes which are redundant, i.e. for which there is no new index //ignore the hyperplanes which are redundant, i.e. for which there is no new index
if ((StormEigen::Index) oldToNewIndexMapping[oldHyperplaneIndex] < relevantVector.rows()){
if ((Eigen::Index) oldToNewIndexMapping[oldHyperplaneIndex] < relevantVector.rows()){
vertexSets[oldToNewIndexMapping[oldHyperplaneIndex]].insert(resultVertices.size()); vertexSets[oldToNewIndexMapping[oldHyperplaneIndex]].insert(resultVertices.size());
} }
} }
@ -120,10 +120,10 @@ namespace storm {
bool HyperplaneEnumeration<ValueType>::linearDependenciesFilter(std::vector<uint_fast64_t> const& subset, uint_fast64_t const& item, EigenMatrix const& A) { bool HyperplaneEnumeration<ValueType>::linearDependenciesFilter(std::vector<uint_fast64_t> const& subset, uint_fast64_t const& item, EigenMatrix const& A) {
EigenMatrix subMatrix(subset.size() + 1, A.cols()); EigenMatrix subMatrix(subset.size() + 1, A.cols());
for (uint_fast64_t i = 0; i < subset.size(); ++i){ for (uint_fast64_t i = 0; i < subset.size(); ++i){
subMatrix.row((StormEigen::Index) i) = A.row(StormEigen::Index (subset[i]));
subMatrix.row((Eigen::Index) i) = A.row(Eigen::Index (subset[i]));
} }
subMatrix.row(subset.size()) = A.row(item); subMatrix.row(subset.size()) = A.row(item);
StormEigen::FullPivLU<EigenMatrix> lUMatrix( subMatrix );
Eigen::FullPivLU<EigenMatrix> lUMatrix( subMatrix );
if (lUMatrix.rank() < subMatrix.rows()){ if (lUMatrix.rank() < subMatrix.rows()){
//Linear dependent! //Linear dependent!
return false; return false;

4
src/storm/storage/geometry/nativepolytopeconversion/HyperplaneEnumeration.h
View File

@ -12,8 +12,8 @@ namespace storm{
class HyperplaneEnumeration { class HyperplaneEnumeration {
public: public:
typedef StormEigen::Matrix<ValueType, StormEigen::Dynamic, StormEigen::Dynamic> EigenMatrix;
typedef StormEigen::Matrix<ValueType, StormEigen::Dynamic, 1> EigenVector;
typedef Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic> EigenMatrix;
typedef Eigen::Matrix<ValueType, Eigen::Dynamic, 1> EigenVector;
HyperplaneEnumeration() = default; HyperplaneEnumeration() = default;
virtual ~HyperplaneEnumeration() = default; virtual ~HyperplaneEnumeration() = default;

10
src/storm/storage/geometry/nativepolytopeconversion/QuickHull.cpp
View File

@ -117,7 +117,7 @@ namespace storm {
vectorMatrix.col(i) << points[subset[i]], storm::utility::one<ValueType>(); vectorMatrix.col(i) << points[subset[i]], storm::utility::one<ValueType>();
} }
vectorMatrix.col(subset.size()) << points[item], storm::utility::one<ValueType>(); vectorMatrix.col(subset.size()) << points[item], storm::utility::one<ValueType>();
return (vectorMatrix.fullPivLu().rank() > (StormEigen::Index) subset.size());
return (vectorMatrix.fullPivLu().rank() > (Eigen::Index) subset.size());
} }
template<typename ValueType> template<typename ValueType>
@ -172,7 +172,7 @@ namespace storm {
// clear the additionally added points. Note that the order of the points might have changed // clear the additionally added points. Note that the order of the points might have changed
storm::storage::BitVector keptPoints(points.size(), true); storm::storage::BitVector keptPoints(points.size(), true);
for (uint_fast64_t pointIndex = 0; pointIndex < points.size(); ++pointIndex) { for (uint_fast64_t pointIndex = 0; pointIndex < points.size(); ++pointIndex) {
for (StormEigen::Index row = 0; row < resultMatrix.rows(); ++row) {
for (Eigen::Index row = 0; row < resultMatrix.rows(); ++row) {
if ((resultMatrix.row(row) * points[pointIndex])(0) > resultVector(row)) { if ((resultMatrix.row(row) * points[pointIndex])(0) > resultVector(row)) {
keptPoints.set(pointIndex, false); keptPoints.set(pointIndex, false);
break; break;
@ -184,7 +184,7 @@ namespace storm {
if (generateRelevantVerticesAndVertexSets) { if (generateRelevantVerticesAndVertexSets) {
storm::storage::BitVector keptVertices(relevantVertices.size(), true); storm::storage::BitVector keptVertices(relevantVertices.size(), true);
for (uint_fast64_t vertexIndex = 0; vertexIndex < relevantVertices.size(); ++vertexIndex) { for (uint_fast64_t vertexIndex = 0; vertexIndex < relevantVertices.size(); ++vertexIndex) {
for (StormEigen::Index row = 0; row < resultMatrix.rows(); ++row) {
for (Eigen::Index row = 0; row < resultMatrix.rows(); ++row) {
if ((resultMatrix.row(row) * relevantVertices[vertexIndex])(0) > resultVector(row)) { if ((resultMatrix.row(row) * relevantVertices[vertexIndex])(0) > resultVector(row)) {
keptVertices.set(vertexIndex, false); keptVertices.set(vertexIndex, false);
break; break;
@ -389,7 +389,7 @@ namespace storm {
std::unordered_map<EigenVector, std::vector<uint_fast64_t>> relevantVerticesMap; std::unordered_map<EigenVector, std::vector<uint_fast64_t>> relevantVerticesMap;
relevantVerticesMap.reserve(points.size()); relevantVerticesMap.reserve(points.size());
for (uint_fast64_t vertexIndex = 0; vertexIndex < hyperplanesOnVertexCounter.size(); ++vertexIndex){ for (uint_fast64_t vertexIndex = 0; vertexIndex < hyperplanesOnVertexCounter.size(); ++vertexIndex){
if ((StormEigen::Index) hyperplanesOnVertexCounter[vertexIndex] >= points.front().rows()){
if ((Eigen::Index) hyperplanesOnVertexCounter[vertexIndex] >= points.front().rows()){
auto mapEntry = relevantVerticesMap.insert(typename std::unordered_map<EigenVector, std::vector<uint_fast64_t>>::value_type(points[vertexIndex], std::vector<uint_fast64_t>())).first; auto mapEntry = relevantVerticesMap.insert(typename std::unordered_map<EigenVector, std::vector<uint_fast64_t>>::value_type(points[vertexIndex], std::vector<uint_fast64_t>())).first;
mapEntry->second.push_back(vertexIndex); mapEntry->second.push_back(vertexIndex);
} }
@ -408,7 +408,7 @@ namespace storm {
for (auto& vertexVector : vertexSets){ for (auto& vertexVector : vertexSets){
std::set<uint_fast64_t> vertexSet; std::set<uint_fast64_t> vertexSet;
for (auto const& oldIndex : vertexVector){ for (auto const& oldIndex : vertexVector){
if ((StormEigen::Index) hyperplanesOnVertexCounter[oldIndex] >= points.front().rows()){
if ((Eigen::Index) hyperplanesOnVertexCounter[oldIndex] >= points.front().rows()){
vertexSet.insert(oldToNewIndexMapping[oldIndex]); vertexSet.insert(oldToNewIndexMapping[oldIndex]);
} }
} }

4
src/storm/storage/geometry/nativepolytopeconversion/QuickHull.h
View File

@ -14,8 +14,8 @@ namespace storm {
class QuickHull { class QuickHull {
public: public:
typedef StormEigen::Matrix<ValueType, StormEigen::Dynamic, StormEigen::Dynamic> EigenMatrix;
typedef StormEigen::Matrix<ValueType, StormEigen::Dynamic, 1> EigenVector;
typedef Eigen::Matrix<ValueType, Eigen::Dynamic, Eigen::Dynamic> EigenMatrix;
typedef Eigen::Matrix<ValueType, Eigen::Dynamic, 1> EigenVector;
QuickHull() = default; QuickHull() = default;
~QuickHull() = default; ~QuickHull() = default;

8
src/storm/storage/geometry/nativepolytopeconversion/SubsetEnumerator.cpp
View File

@ -89,10 +89,10 @@ namespace storm {
} }
template class SubsetEnumerator<std::nullptr_t>; template class SubsetEnumerator<std::nullptr_t>;
template class SubsetEnumerator<std::vector<StormEigen::Matrix<double, StormEigen::Dynamic, 1>>>;
template class SubsetEnumerator<std::vector<StormEigen::Matrix<storm::RationalNumber, StormEigen::Dynamic, 1>>>;
template class SubsetEnumerator<StormEigen::Matrix<double, StormEigen::Dynamic, StormEigen::Dynamic>>;
template class SubsetEnumerator<StormEigen::Matrix<storm::RationalNumber, StormEigen::Dynamic, StormEigen::Dynamic>>;
template class SubsetEnumerator<std::vector<Eigen::Matrix<double, Eigen::Dynamic, 1>>>;
template class SubsetEnumerator<std::vector<Eigen::Matrix<storm::RationalNumber, Eigen::Dynamic, 1>>>;
template class SubsetEnumerator<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>>;
template class SubsetEnumerator<Eigen::Matrix<storm::RationalNumber, Eigen::Dynamic, Eigen::Dynamic>>;
} }
} }
} }

8
src/storm/utility/eigen.h
View File

@ -1,6 +1,6 @@
#pragma once #pragma once
// Include this utility header so we can access utility function from Eigen.
// Include these utility headers so we can access utility function from Eigen.
#include "storm/utility/constants.h" #include "storm/utility/constants.h"
#if defined(__clang__) #if defined(__clang__)
@ -14,9 +14,9 @@
#endif #endif
// Finally include the parts of Eigen we need. // Finally include the parts of Eigen we need.
#include <StormEigen/Dense>
#include <StormEigen/Sparse>
#include <unsupported/StormEigen/IterativeSolvers>
#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <unsupported/Eigen/IterativeSolvers>
#if defined(__clang__) #if defined(__clang__)
#pragma clang diagnostic pop #pragma clang diagnostic pop

4
src/test/storm/modelchecker/csl/LraCtmcCslModelCheckerTest.cpp
View File

@ -44,6 +44,7 @@ namespace {
env.solver().gmmxx().setMethod(storm::solver::GmmxxLinearEquationSolverMethod::Gmres); env.solver().gmmxx().setMethod(storm::solver::GmmxxLinearEquationSolverMethod::Gmres);
env.solver().gmmxx().setPreconditioner(storm::solver::GmmxxLinearEquationSolverPreconditioner::Ilu); env.solver().gmmxx().setPreconditioner(storm::solver::GmmxxLinearEquationSolverPreconditioner::Ilu);
env.solver().lra().setDetLraMethod(storm::solver::LraMethod::GainBiasEquations); 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; return env;
} }
}; };
@ -61,6 +62,7 @@ namespace {
env.solver().gmmxx().setMethod(storm::solver::GmmxxLinearEquationSolverMethod::Gmres); env.solver().gmmxx().setMethod(storm::solver::GmmxxLinearEquationSolverMethod::Gmres);
env.solver().gmmxx().setPreconditioner(storm::solver::GmmxxLinearEquationSolverPreconditioner::Ilu); env.solver().gmmxx().setPreconditioner(storm::solver::GmmxxLinearEquationSolverPreconditioner::Ilu);
env.solver().lra().setDetLraMethod(storm::solver::LraMethod::GainBiasEquations); 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; return env;
} }
}; };
@ -343,7 +345,7 @@ namespace {
std::unique_ptr<storm::modelchecker::CheckResult> result; std::unique_ptr<storm::modelchecker::CheckResult> result;
result = checker->check(this->env(), tasks[0]); result = checker->check(this->env(), tasks[0]);
EXPECT_NEAR(this->parseNumber("0.93458866427696596"), this->getQuantitativeResultAtInitialState(model, result), this->precision());
EXPECT_NEAR(this->parseNumber("6201111489217/6635130141055"), this->getQuantitativeResultAtInitialState(model, result), this->precision());
} }
TYPED_TEST(LraCtmcCslModelCheckerTest, Polling) { TYPED_TEST(LraCtmcCslModelCheckerTest, Polling) {

Write Preview
Loading…
Cancel
Save