diff --git a/src/storm/solver/LpMinMaxLinearEquationSolver.cpp b/src/storm/solver/LpMinMaxLinearEquationSolver.cpp new file mode 100644 index 000000000..2787a1b6f --- /dev/null +++ b/src/storm/solver/LpMinMaxLinearEquationSolver.cpp @@ -0,0 +1,141 @@ +#include "storm/solver/LpMinMaxLinearEquationSolver.h" + +#include "storm/utility/vector.h" +#include "storm/utility/macros.h" +#include "storm/exceptions/InvalidSettingsException.h" +#include "storm/exceptions/InvalidOperationException.h" +#include "storm/exceptions/UnexpectedException.h" + +namespace storm { + namespace solver { + + template<typename ValueType> + LpMinMaxLinearEquationSolver<ValueType>::LpMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, std::unique_ptr<LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory, std::unique_ptr<storm::utility::solver::LpSolverFactory>&& lpSolverFactory) : StandardMinMaxLinearEquationSolver<ValueType>(A, std::move(linearEquationSolverFactory)), lpSolverFactory(std::move(lpSolverFactory)) { + // Intentionally left empty. + } + + template<typename ValueType> + LpMinMaxLinearEquationSolver<ValueType>::LpMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType>&& A, std::unique_ptr<LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory, std::unique_ptr<storm::utility::solver::LpSolverFactory>&& lpSolverFactory) : StandardMinMaxLinearEquationSolver<ValueType>(std::move(A), std::move(linearEquationSolverFactory)), lpSolverFactory(std::move(lpSolverFactory)) { + // Intentionally left empty. + } + + template<typename ValueType> + bool LpMinMaxLinearEquationSolver<ValueType>::solveEquations(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const { + + // Set up the LP solver + std::unique_ptr<storm::solver::LpSolver> solver = lpSolverFactory->create("MinMaxLinearEquationSolver"); + solver->setOptimizationDirection(invert(dir)); + + // Create a variable for each row group + std::vector<storm::expressions::Variable> variables; + variables.reserve(this->A.getRowGroupCount()); + for (uint64_t rowGroup = 0; rowGroup < this->A.getRowGroupCount(); ++rowGroup) { + if (this->lowerBound) { + if (this->upperBound) { + variables.push_back(solver->addBoundedContinuousVariable("x" + std::to_string(rowGroup), storm::utility::convertNumber<double>(this->lowerBound.get()), storm::utility::convertNumber<double>(this->upperBound.get()), 1.0)); + } else { + variables.push_back(solver->addLowerBoundedContinuousVariable("x" + std::to_string(rowGroup), storm::utility::convertNumber<double>(this->lowerBound.get()), 1.0)); + } + } else { + if (this->upperBound) { + variables.push_back(solver->addUpperBoundedContinuousVariable("x" + std::to_string(rowGroup), storm::utility::convertNumber<double>(this->upperBound.get()), 1.0)); + } else { + variables.push_back(solver->addUnboundedContinuousVariable("x" + std::to_string(rowGroup), 1.0)); + } + } + } + solver->update(); + + // Add a constraint for each row + for (uint64_t rowGroup = 0; rowGroup < this->A.getRowGroupCount(); ++rowGroup) { + for (uint64_t row = this->A.getRowGroupIndices()[rowGroup]; row < this->A.getRowGroupIndices()[rowGroup + 1]; ++row) { + storm::expressions::Expression rowConstraint = solver->getConstant(storm::utility::convertNumber<double>(b[row])); + for (auto const& entry : this->A.getRow(row)) { + rowConstraint = rowConstraint + (solver->getConstant(storm::utility::convertNumber<double>(entry.getValue())) * variables[entry.getColumn()].getExpression()); + } + if (minimize(dir)) { + rowConstraint = variables[rowGroup].getExpression() <= rowConstraint; + } else { + rowConstraint = variables[rowGroup].getExpression() >= rowConstraint; + } + solver->addConstraint("row" + std::to_string(row), rowConstraint); + } + } + + // Invoke optimization + solver->optimize(); + STORM_LOG_THROW(!solver->isInfeasible(), storm::exceptions::UnexpectedException, "The MinMax equation system is infeasible."); + STORM_LOG_THROW(!solver->isUnbounded(), storm::exceptions::UnexpectedException, "The MinMax equation system is unbounded."); + STORM_LOG_THROW(solver->isOptimal(), storm::exceptions::UnexpectedException, "Unable to find optimal solution for MinMax equation system."); + + // write the solution into the solution vector + STORM_LOG_ASSERT(x.size() == variables.size(), "Dimension of x-vector does not match number of varibales."); + auto xIt = x.begin(); + auto vIt = variables.begin(); + for (; xIt != x.end(); ++xIt, ++vIt) { + *xIt = storm::utility::convertNumber<ValueType>(solver->getContinuousValue(*vIt)); + } + return true; + } + + template<typename ValueType> + void LpMinMaxLinearEquationSolver<ValueType>::clearCache() const { + StandardMinMaxLinearEquationSolver<ValueType>::clearCache(); + } + + template<typename ValueType> + LpMinMaxLinearEquationSolverFactory<ValueType>::LpMinMaxLinearEquationSolverFactory(bool trackScheduler) : StandardMinMaxLinearEquationSolverFactory<ValueType>(MinMaxMethodSelection::LinearProgramming, trackScheduler), lpSolverFactory(std::make_unique<storm::utility::solver::LpSolverFactory>()) { + // Intentionally left empty + } + + template<typename ValueType> + LpMinMaxLinearEquationSolverFactory<ValueType>::LpMinMaxLinearEquationSolverFactory(std::unique_ptr<storm::utility::solver::LpSolverFactory>&& lpSolverFactory, bool trackScheduler) : StandardMinMaxLinearEquationSolverFactory<ValueType>(MinMaxMethodSelection::LinearProgramming, trackScheduler), lpSolverFactory(std::move(lpSolverFactory)) { + // Intentionally left empty + } + + template<typename ValueType> + LpMinMaxLinearEquationSolverFactory<ValueType>::LpMinMaxLinearEquationSolverFactory(std::unique_ptr<LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory, std::unique_ptr<storm::utility::solver::LpSolverFactory>&& lpSolverFactory, bool trackScheduler) : StandardMinMaxLinearEquationSolverFactory<ValueType>(std::move(linearEquationSolverFactory), MinMaxMethodSelection::LinearProgramming, trackScheduler), lpSolverFactory(std::move(lpSolverFactory)) { + // Intentionally left empty + } + + template<typename ValueType> + std::unique_ptr<MinMaxLinearEquationSolver<ValueType>> LpMinMaxLinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType> const& matrix) const { + STORM_LOG_ASSERT(this->linearEquationSolverFactory, "Linear equation solver factory not initialized."); + STORM_LOG_ASSERT(this->lpSolverFactory, "Lp solver factory not initialized."); + + std::unique_ptr<MinMaxLinearEquationSolver<ValueType>> result = std::make_unique<LpMinMaxLinearEquationSolver<ValueType>>(std::move(matrix), this->linearEquationSolverFactory->clone(), this->lpSolverFactory->clone()); + result->setTrackScheduler(this->isTrackSchedulerSet()); + return result; + } + + template<typename ValueType> + std::unique_ptr<MinMaxLinearEquationSolver<ValueType>> LpMinMaxLinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType>&& matrix) const { + STORM_LOG_ASSERT(this->linearEquationSolverFactory, "Linear equation solver factory not initialized."); + STORM_LOG_ASSERT(this->lpSolverFactory, "Lp solver factory not initialized."); + + std::unique_ptr<MinMaxLinearEquationSolver<ValueType>> result = std::make_unique<LpMinMaxLinearEquationSolver<ValueType>>(std::move(matrix), this->linearEquationSolverFactory->clone(), this->lpSolverFactory->clone()); + result->setTrackScheduler(this->isTrackSchedulerSet()); + return result; + } + + template<typename ValueType> + void LpMinMaxLinearEquationSolverFactory<ValueType>::setMinMaxMethod(MinMaxMethodSelection const& newMethod) { + STORM_LOG_THROW(newMethod == MinMaxMethodSelection::LinearProgramming, storm::exceptions::InvalidOperationException, "The factory can only create linear programming based MinMax solvers."); + MinMaxLinearEquationSolverFactory<ValueType>::setMinMaxMethod(newMethod); + } + + template<typename ValueType> + void LpMinMaxLinearEquationSolverFactory<ValueType>::setMinMaxMethod(MinMaxMethod const& newMethod) { + STORM_LOG_THROW(newMethod == MinMaxMethod::LinearProgramming, storm::exceptions::InvalidOperationException, "The factory can only create linear programming based MinMax solvers."); + MinMaxLinearEquationSolverFactory<ValueType>::setMinMaxMethod(newMethod); + } + + template class LpMinMaxLinearEquationSolver<double>; + template class LpMinMaxLinearEquationSolverFactory<double>; + +#ifdef STORM_HAVE_CARL + template class LpMinMaxLinearEquationSolver<storm::RationalNumber>; + template class LpMinMaxLinearEquationSolverFactory<storm::RationalNumber>; +#endif + } +} diff --git a/src/storm/solver/LpMinMaxLinearEquationSolver.h b/src/storm/solver/LpMinMaxLinearEquationSolver.h new file mode 100644 index 000000000..bd15326cb --- /dev/null +++ b/src/storm/solver/LpMinMaxLinearEquationSolver.h @@ -0,0 +1,41 @@ +#pragma once + +#include "storm/solver/LpSolver.h" +#include "storm/solver/StandardMinMaxLinearEquationSolver.h" +#include "storm/utility/solver.h" + +namespace storm { + namespace solver { + + template<typename ValueType> + class LpMinMaxLinearEquationSolver : public StandardMinMaxLinearEquationSolver<ValueType> { + public: + LpMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, std::unique_ptr<LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory, std::unique_ptr<storm::utility::solver::LpSolverFactory>&& lpSolverFactory); + LpMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType>&& A, std::unique_ptr<LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory, std::unique_ptr<storm::utility::solver::LpSolverFactory>&& lpSolverFactory); + + virtual bool solveEquations(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const override; + + virtual void clearCache() const override; + + private: + std::unique_ptr<storm::utility::solver::LpSolverFactory> lpSolverFactory; + }; + + template<typename ValueType> + class LpMinMaxLinearEquationSolverFactory : public StandardMinMaxLinearEquationSolverFactory<ValueType> { + public: + LpMinMaxLinearEquationSolverFactory(bool trackScheduler = false); + LpMinMaxLinearEquationSolverFactory(std::unique_ptr<storm::utility::solver::LpSolverFactory>&& lpSolverFactory, bool trackScheduler = false); + LpMinMaxLinearEquationSolverFactory(std::unique_ptr<LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory, std::unique_ptr<storm::utility::solver::LpSolverFactory>&& lpSolverFactory, bool trackScheduler = false); + + virtual std::unique_ptr<MinMaxLinearEquationSolver<ValueType>> create(storm::storage::SparseMatrix<ValueType> const& matrix) const override; + virtual std::unique_ptr<MinMaxLinearEquationSolver<ValueType>> create(storm::storage::SparseMatrix<ValueType>&& matrix) const override; + + virtual void setMinMaxMethod(MinMaxMethodSelection const& newMethod) override; + virtual void setMinMaxMethod(MinMaxMethod const& newMethod) override; + + private: + std::unique_ptr<storm::utility::solver::LpSolverFactory> lpSolverFactory; + }; + } +} diff --git a/src/storm/utility/solver.cpp b/src/storm/utility/solver.cpp index fcee52548..bc632ce9d 100644 --- a/src/storm/utility/solver.cpp +++ b/src/storm/utility/solver.cpp @@ -46,18 +46,34 @@ namespace storm { return LpSolverFactory::create(name, storm::solver::LpSolverTypeSelection::FROMSETTINGS); } + std::unique_ptr<LpSolverFactory> LpSolverFactory::clone() const { + return std::make_unique<LpSolverFactory>(*this); + } + std::unique_ptr<storm::solver::LpSolver> GlpkLpSolverFactory::create(std::string const& name) const { return LpSolverFactory::create(name, storm::solver::LpSolverTypeSelection::Glpk); } + std::unique_ptr<LpSolverFactory> GlpkLpSolverFactory::clone() const { + return std::make_unique<GlpkLpSolverFactory>(*this); + } + std::unique_ptr<storm::solver::LpSolver> GurobiLpSolverFactory::create(std::string const& name) const { return LpSolverFactory::create(name, storm::solver::LpSolverTypeSelection::Gurobi); } - + + std::unique_ptr<LpSolverFactory> GurobiLpSolverFactory::clone() const { + return std::make_unique<GurobiLpSolverFactory>(*this); + } + std::unique_ptr<storm::solver::LpSolver> Z3LpSolverFactory::create(std::string const& name) const { return LpSolverFactory::create(name, storm::solver::LpSolverTypeSelection::Z3); } + std::unique_ptr<LpSolverFactory> Z3LpSolverFactory::clone() const { + return std::make_unique<Z3LpSolverFactory>(*this); + } + std::unique_ptr<storm::solver::LpSolver> getLpSolver(std::string const& name, storm::solver::LpSolverTypeSelection solvType) { std::unique_ptr<storm::utility::solver::LpSolverFactory> factory(new LpSolverFactory()); return factory->create(name, solvType); diff --git a/src/storm/utility/solver.h b/src/storm/utility/solver.h index 4fbc0745e..3f2bd6f28 100644 --- a/src/storm/utility/solver.h +++ b/src/storm/utility/solver.h @@ -70,21 +70,27 @@ namespace storm { */ virtual std::unique_ptr<storm::solver::LpSolver> create(std::string const& name) const; virtual std::unique_ptr<storm::solver::LpSolver> create(std::string const& name, storm::solver::LpSolverTypeSelection solvType) const; + + virtual std::unique_ptr<LpSolverFactory> clone() const; + }; class GlpkLpSolverFactory : public LpSolverFactory { public: virtual std::unique_ptr<storm::solver::LpSolver> create(std::string const& name) const override; + virtual std::unique_ptr<LpSolverFactory> clone() const override; }; class GurobiLpSolverFactory : public LpSolverFactory { public: virtual std::unique_ptr<storm::solver::LpSolver> create(std::string const& name) const override; + virtual std::unique_ptr<LpSolverFactory> clone() const override; }; class Z3LpSolverFactory : public LpSolverFactory { public: virtual std::unique_ptr<storm::solver::LpSolver> create(std::string const& name) const override; + virtual std::unique_ptr<LpSolverFactory> clone() const override; }; std::unique_ptr<storm::solver::LpSolver> getLpSolver(std::string const& name, storm::solver::LpSolverTypeSelection solvType = storm::solver::LpSolverTypeSelection::FROMSETTINGS) ;