diff --git a/CMakeLists.txt b/CMakeLists.txt index 82b9289..1ec41fe 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -12,25 +12,20 @@ set(CMAKE_CXX_STANDARD_REQUIRED ON) set(SOURCE_DIR "${PROJECT_SOURCE_DIR}/src") set(INCLUDE_DIR "${PROJECT_SOURCE_DIR}/inc") -file(GLOB SOURCES "${SOURCE_DIR}/*.cpp$") -file(GLOB HEADERS "${SOURCE_INCLUDE}/*.hpp$") +file(GLOB SOURCES "${SOURCE_DIR}/*.cpp") +file(GLOB HEADERS "${INCLUDE_DIR}/*.hpp") #include_directories( # ${GUROBI_INCLUDE_DIRS} #) -message(${SOURCE_DIR}) - aux_source_directory(${SOURCE_DIR} src) aux_source_directory(${INCLUDE_DIR} inc) - -add_executable(Solver - src/main.cpp src/solver.cpp inc/solver.hpp - ) +add_executable(Solver ${SOURCES} ${HEADERS}) #target_link_libraries(Solver -# ${GUROBI_LIBRARIES} +# ${GUROBI_LIBRARIES} #) target_include_directories(Solver PRIVATE ${SOURCE_DIR} ${INCLUDE_DIR}) diff --git a/inc/common.hpp b/inc/common.hpp new file mode 100644 index 0000000..5dc0dbd --- /dev/null +++ b/inc/common.hpp @@ -0,0 +1,163 @@ +#pragma once +#include + +#define MDL_MINIMIZE 1 +#define MDL_MAXIMIZE -1 + +#define LOADED 0 +#define OPTIMAL 1 +#define INFEASIBLE 2 +#define INF_OR_UNBD 3 +#define UNBOUNDED 4 +#define CUTOFF 5 +#define ITERATION_LIMIT 6 +#define NODE_LIMIT 7 +#define TIME_LIMIT 8 +#define SOLUTION_LIMIT 9 +#define INTERRUPTED 10 +#define NUMERIC 11 +#define SUBOPTIMAL 12 +#define INPROGRESS 13 +#define USER_OBJ_LIMIT 14 +#define WORK_LIMIT 15 +#define MEM_LIMIT 16 +#define PIVOT 17 + +namespace sv { + using matrix = std::vector>; + using rtn = int; + class Expr; + class Var; + + enum class IntAttr { + NumConstrs, + NumVars, + NumIntVars, + NumBinVars, + ModelSense, + IsMIP, + IsMultiObj, + Status, + SolCount, + Lazy, + NumObj, + NumCol + }; + + enum class DoubleAttr { + Runtime, + Work, + ObjCon, + LB, + UB, + Obj, + Start, + RHS, + Coeff, + MaxCoeff, + MinCoeff, + MaxBound, + MinBound, + ObjVal, + MIPGap, + IterCount, + NodeCount, + X, + Slack, + }; + + enum class StringAttr { + ModelName, + VarName, + ConstrName, + QCName, + GenConstrName, + ObjNName, + ScenNName, + BatchID, + VTag, + CTag, + QCTag, + BatchErrorMessage + }; + + enum class ConstrOper { + LESS_EQUAL, + GREATER_EQUAL, + EQUAL + }; + + class Var { + public: + Var(double coef = 1) :col(0), val(0), coeffs(coef) {}; + double get(DoubleAttr attr); + int get(IntAttr attr); + friend class Model; + friend class Expr; + private: + double coeffs; + double val; + int col; + }; + + Expr operator+(const Expr& x, const Expr& y); + Expr operator-(const Expr& x, const Expr& y); + Expr operator+(const Expr& x); + Expr operator+(Var x, Var y); + Expr operator+(Var x, double a); + Expr operator+(double a, Var x); + Expr operator-(const Expr& x); + Expr operator-(Var x); + Expr operator-(Var x, Var y); + Expr operator-(Var x, double a); + Expr operator-(double a, Var x); + Expr operator*(double a, Var x); + Expr operator*(Var x, double a); + Expr operator*(const Expr& x, double a); + Expr operator*(double a, const Expr& x); + Expr operator/(Var x, double a); + Expr operator/(const Expr& x, double a); + + + class Expr + { + private: + double constant; + std::vector coeffs; + std::vector vars; + + public: + Expr(const Expr& expr) = default; + Expr(double constant = 0.0); + Expr(Var var, double coeff = 1.0); + + friend class Model; + + friend Expr operator+(const Expr& x, const Expr& y); + friend Expr operator+(const Expr& x); + friend Expr operator+(Var x, Var y); + friend Expr operator+(Var x, double a); + friend Expr operator+(double a, Var x); + friend Expr operator-(const Expr& x, const Expr& y); + friend Expr operator-(const Expr& x); + friend Expr operator-(Var x); + friend Expr operator-(Var x, Var y); + friend Expr operator-(Var x, double a); + friend Expr operator-(double a, Var x); + friend Expr operator*(double a, Var x); + friend Expr operator*(Var x, double a); + friend Expr operator*(const Expr& x, double a); + friend Expr operator*(double a, const Expr& x); + friend Expr operator/(Var x, double a); + friend Expr operator/(const Expr& x, double a); + + + Expr operator=(const Expr& rhs); + void operator+=(const Expr& expr); + void operator-=(const Expr& expr); + void operator*=(double mult); + void operator/=(double a); + Expr operator+(const Expr& rhs); + Expr operator-(const Expr& rhs); + }; +} diff --git a/inc/node.hpp b/inc/node.hpp new file mode 100644 index 0000000..7c06824 --- /dev/null +++ b/inc/node.hpp @@ -0,0 +1,5 @@ +#pragma once + +class Node { + +}; \ No newline at end of file diff --git a/inc/solver.hpp b/inc/solver.hpp index 0fdd6cb..ca974c7 100644 --- a/inc/solver.hpp +++ b/inc/solver.hpp @@ -1,327 +1,34 @@ #pragma once -#include -using std::vector; +#include -class Expr; -class Var; +namespace sv { -#define MDL_MINIMIZE 1 -#define MDL_MAXIMIZE -1 + class Model + { + public: + Model(); + Var* addVars(int col); + void addConstr(const Expr& expr, ConstrOper sense, double rhs); + void setObjective(Expr obje, int sense = MDL_MAXIMIZE); + void print(); + rtn optimize(); + double get(DoubleAttr attr); + int get(IntAttr attr); + private: + double _simplex(); + rtn _pivot(std::pair& p); + void _gaussian(std::pair p); -Expr operator+(const Expr& x, const Expr& y); -Expr operator-(const Expr& x, const Expr& y); -Expr operator+(const Expr& x); -Expr operator+(Var x, Var y); -Expr operator+(Var x, double a); -Expr operator+(double a, Var x); -Expr operator-(const Expr& x); -Expr operator-(Var x); -Expr operator-(Var x, Var y); -Expr operator-(Var x, double a); -Expr operator-(double a, Var x); -Expr operator*(double a, Var x); -Expr operator*(Var x, double a); -Expr operator*(const Expr& x, double a); -Expr operator*(double a, const Expr& x); -Expr operator/(Var x, double a); -Expr operator/(const Expr& x, double a); + Var* vars; -typedef vector> Matrix; - -enum class Rtn { - LOADED, - OPTIMAL, - INFEASIBLE, - INF_OR_UNBD, - UNBOUNDED, - CUTOFF, - ITERATION_LIMIT, - NODE_LIMIT, - TIME_LIMIT, - SOLUTION_LIMIT, - INTERRUPTED, - NUMERIC, - SUBOPTIMAL, - INPROGRESS, - USER_OBJ_LIMIT, - WORK_LIMIT, - MEM_LIMIT, - PIVOT, -}; - -enum class IntAttr { - NumConstrs, - NumVars, - NumSOS, - NumQConstrs, - NumGenConstrs, - NumNZs, - NumQNZs, - NumQCNZs, - NumIntVars, - NumBinVars, - NumPWLObjVars, - ModelSense, - IsMIP, - IsQP, - IsQCP, - IsMultiObj, - Status, - ConcurrentWinMethod, - SolCount, - BarIterCount, - VBasis, - CBasis, - PWLObjCvx, - BranchPriority, - VarPreStat, - BoundVioIndex, - BoundSVioIndex, - ConstrVioIndex, - ConstrSVioIndex, - ConstrResidualIndex, - ConstrSResidualIndex, - DualVioIndex, - DualSVioIndex, - DualResidualIndex, - DualSResidualIndex, - ComplVioIndex, - IntVioIndex, - IISMinimal, - IISLB, - IISUB, - IISConstr, - IISSOS, - IISQConstr, - IISGenConstr, - IISLBForce, - IISUBForce, - IISConstrForce, - IISSOSForce, - IISQConstrForce, - IISGenConstrForce, - TuneResultCount, - Lazy, - VarHintPri, - ObjNPriority, - NumObj, - GenConstrType, - NumStart, - Partition, - LicenseExpiration, - NumScenarios, - FuncPieces, - BatchErrorCode, - BatchStatus, - Fingerprint, - PoolIgnore, - FuncNonlinear, - ModelColumn -}; - -enum class DoubleAttr{ - Runtime, - Work, - ObjCon, - LB, - UB, - Obj, - Start, - PreFixVal, - RHS, - QCRHS, - Coeff, - MaxCoeff, - MinCoeff, - MaxBound, - MinBound, - MaxObjCoeff, - MinObjCoeff, - MaxRHS, - MinRHS, - MaxQCRHS, - MinQCRHS, - MaxQCCoeff, - MinQCCoeff, - MaxQCLCoeff, - MinQCLCoeff, - MaxQObjCoeff, - MinQObjCoeff, - ObjVal, - ObjBound, - ObjBoundC, - MIPGap, - IterCount, - NodeCount, - X, - RC, - Pi, - QCPi, - Slack, - QCSlack, - MaxVio, - BoundVio, - BoundSVio, - BoundVioSum, - BoundSVioSum, - ConstrVio, - ConstrSVio, - ConstrVioSum, - ConstrSVioSum, - ConstrResidual, - ConstrSResidual, - ConstrResidualSum, - ConstrSResidualSum, - DualVio, - DualSVio, - DualVioSum, - DualSVioSum, - DualResidual, - DualSResidual, - DualResidualSum, - DualSResidualSum, - ComplVio, - ComplVioSum, - IntVio, - IntVioSum, - Kappa, - KappaExact, - SAObjLow, - SAObjUp, - SALBLow, - SALBUp, - SARHSLow, - SAUBLow, - SAUBUp, - SARHSUp, - Xn, - FarkasProof, - FarkasDual, - UnbdRay, - PStart, - DStart, - BarX, - VarHintVal, - ObjN, - ObjNCon, - ObjNWeight, - ObjNRelTol, - ObjNAbsTol, - ObjNVal, - OpenNodeCount, - PoolObjBound, - PoolObjVal, - ScenNLB, - ScenNUB, - ScenNObj, - ScenNRHS, - ScenNX, - ScenNObjBound, - ScenNObjVal, - FuncPieceError, - FuncPieceLength, - FuncPieceRatio, - DNumNZs -}; - -enum class StringAttr { - ModelName, - VarName, - ConstrName, - QCName, - GenConstrName, - ObjNName, - ScenNName, - BatchID, - VTag, - CTag, - QCTag, - BatchErrorMessage -}; - -enum class ConstrOper { - LESS_EQUAL, - GREATER_EQUAL, - EQUAL -}; - -class Var { -public: - Var(double coef = 1) :col(0), val(0), coeffs(coef) {}; - double get(DoubleAttr attr); - int get(IntAttr attr); - friend class Model; - friend class Expr; -private: - double coeffs; - double val; - int col; -}; - -class Expr -{ -private: - double constant; - std::vector coeffs; - std::vector vars; - -public: - Expr(const Expr& expr) = default; - Expr(double constant = 0.0); - Expr(Var var, double coeff = 1.0); - - friend class Model; - - friend Expr operator+(const Expr& x, const Expr& y); - friend Expr operator+(const Expr& x); - friend Expr operator+(Var x, Var y); - friend Expr operator+(Var x, double a); - friend Expr operator+(double a, Var x); - friend Expr operator-(const Expr& x, const Expr& y); - friend Expr operator-(const Expr& x); - friend Expr operator-(Var x); - friend Expr operator-(Var x, Var y); - friend Expr operator-(Var x, double a); - friend Expr operator-(double a, Var x); - friend Expr operator*(double a, Var x); - friend Expr operator*(Var x, double a); - friend Expr operator*(const Expr& x, double a); - friend Expr operator*(double a, const Expr& x); - friend Expr operator/(Var x, double a); - friend Expr operator/(const Expr& x, double a); + matrix mt; + size_t cn, bn; + std::vector basic; + rtn rtn_; + double res_; + int sense; + }; +} - Expr operator=(const Expr& rhs); - void operator+=(const Expr& expr); - void operator-=(const Expr& expr); - void operator*=(double mult); - void operator/=(double a); - Expr operator+(const Expr& rhs); - Expr operator-(const Expr& rhs); -}; -class Model -{ -public: - Model(); - Var* addVars(int col); - void addConstr(const Expr& expr, ConstrOper sense, double rhs); - void setObjective(Expr obje, int sense = MDL_MAXIMIZE); - void print(); - Rtn optimize(); - double get(DoubleAttr attr); - int get(IntAttr attr); -private: - double _simplex(); - Rtn _pivot(std::pair& p); - void _gaussian(std::pair p); - - Var* vars; - - Matrix matrix; - size_t cn, bn; - vector basic; - Rtn rtn; - double res; - int sense; -}; \ No newline at end of file diff --git a/src/common.cpp b/src/common.cpp new file mode 100644 index 0000000..2869025 --- /dev/null +++ b/src/common.cpp @@ -0,0 +1,225 @@ +#include "common.hpp" +using namespace sv; + +double Var::get(DoubleAttr attr) +{ + switch (attr) { + case DoubleAttr::Coeff: + return coeffs; + case DoubleAttr::X: + return val; + } + return -1; +} + +int Var::get(IntAttr attr) +{ + return col; +} + +Expr sv::operator+(const Expr& x, const Expr& y) +{ + Expr exp; + exp.coeffs.resize(std::max(x.coeffs.size(), y.coeffs.size()), 0); + for (int c = 0; c < exp.coeffs.size(); c++) { + if (c < x.coeffs.size() && c < y.coeffs.size()) { + exp.coeffs[c] = x.coeffs[c] + y.coeffs[c]; + } + else { + exp.coeffs[c] = c < x.coeffs.size() ? x.coeffs[c] : y.coeffs[c]; + } + } + exp.constant = x.constant + y.constant; + return exp; +} + +Expr sv::operator+(const Expr& x) +{ + return x; +} + +Expr sv::operator+(Var x, Var y) +{ + Expr exp; + exp.coeffs.resize(std::max(x.get(IntAttr::NumCol) + 1, y.get(IntAttr::NumCol) + 1), 0); + exp.coeffs[x.get(IntAttr::NumCol)] = x.get(DoubleAttr::Coeff); + exp.coeffs[y.get(IntAttr::NumCol)] = y.get(DoubleAttr::Coeff); + return exp; +} + +Expr sv::operator+(Var x, double a) +{ + Expr exp; + exp.coeffs.resize(x.get(IntAttr::NumCol) + 1); + exp.constant = a; + return exp; +} + +Expr sv::operator+(double a, Var x) +{ + return x + a; +} + +Expr sv::operator-(const Expr& x, const Expr& y) +{ + Expr exp; + exp.coeffs.resize(std::max(x.coeffs.size(), y.coeffs.size()), 0); + for (int c = 0; c < exp.coeffs.size(); c++) { + if (c < x.coeffs.size() && c < y.coeffs.size()) { + exp.coeffs[c] = x.coeffs[c] - y.coeffs[c]; + } + else { + exp.coeffs[c] = c < x.coeffs.size() ? x.coeffs[c] : y.coeffs[c]; + } + } + exp.constant = x.constant + y.constant; + return exp; +} + +Expr sv::operator-(const Expr& x) +{ + Expr expr(x); + for (int c = 0; c < expr.coeffs.size(); c++) { + expr.coeffs[c] = -expr.coeffs[c]; + } + expr.constant = -expr.constant; + return expr; +} + +Expr sv::operator-(Var x) +{ + return -Expr(x); +} + +Expr sv::operator-(Var x, Var y) +{ + Expr exp; + exp.coeffs.resize(std::max(x.get(IntAttr::NumCol) + 1, y.get(IntAttr::NumCol) + 1), 0); + exp.coeffs[x.get(IntAttr::NumCol)] = x.get(DoubleAttr::Coeff); + exp.coeffs[y.get(IntAttr::NumCol)] = -y.get(DoubleAttr::Coeff); + return exp; +} + +Expr sv::operator-(Var x, double a) +{ + return x - Var(a); +} + +Expr sv::operator-(double a, Var x) +{ + return x - a; +} + +Expr sv::operator*(double a, Var x) +{ + Expr exp; + exp.coeffs.resize(x.get(IntAttr::NumCol) + 1, 0); + exp.coeffs[x.get(IntAttr::NumCol)] = a * x.get(DoubleAttr::Coeff); + return exp; +} + +Expr sv::operator*(Var x, double a) +{ + return a * x; +} + +Expr sv::operator*(const Expr& x, double a) +{ + Expr exp = x; + for (int c = 0; c < exp.coeffs.size(); c++) { + exp.coeffs[c] *= a; + } + exp.constant *= a; + return exp; +} + +Expr sv::operator*(double a, const Expr& x) +{ + return x * a; +} + +Expr sv::operator/(Var x, double a) +{ + return Expr(x) / a; +} + +Expr sv::operator/(const Expr& x, double a) +{ + Expr exp = x; + for (int c = 0; c < exp.coeffs.size(); c++) { + exp.coeffs[c] /= a; + } + exp.constant /= a; + return exp; +} + +Expr::Expr(double constant) + :constant(constant) +{ +} + +Expr::Expr(Var var, double coeff) +{ + this->coeffs.resize(var.col + 1); + this->coeffs[var.col] = coeff; + this->constant = 0; +} + +Expr Expr::operator=(const Expr& rhs) +{ + return *this; +} + +void Expr::operator+=(const Expr& expr) +{ + coeffs.resize(std::max(coeffs.size(), expr.coeffs.size())); + for (int c = 0; c < expr.coeffs.size(); c++) { + coeffs[c] += expr.coeffs[c]; + } + constant += expr.constant; +} + +void Expr::operator-=(const Expr& expr) +{ + coeffs.resize(std::max(coeffs.size(), expr.coeffs.size())); + for (int c = 0; c < expr.coeffs.size(); c++) { + coeffs[c] -= expr.coeffs[c]; + } + constant -= expr.constant; +} + +void Expr::operator*=(double mult) +{ + for (int c = 0; c < coeffs.size(); c++) { + coeffs[c] *= mult; + } + constant *= mult; +} + +void Expr::operator/=(double a) +{ + for (int c = 0; c < coeffs.size(); c++) { + coeffs[c] /= a; + } + constant /= a; +} + +Expr Expr::operator+(const Expr& rhs) +{ + coeffs.resize(std::max(coeffs.size(), rhs.coeffs.size())); + for (int c = 0; c < rhs.coeffs.size(); c++) { + coeffs[c] += rhs.coeffs[c]; + } + constant += rhs.constant; + return *this; +} + +Expr Expr::operator-(const Expr& rhs) +{ + coeffs.resize(std::max(coeffs.size(), rhs.coeffs.size())); + for (int c = 0; c < rhs.coeffs.size(); c++) { + coeffs[c] -= rhs.coeffs[c]; + } + constant -= rhs.constant; + return *this; +} diff --git a/src/main.cpp b/src/main.cpp index 5430942..37e8956 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -4,6 +4,7 @@ #include "solver.hpp" using namespace std; +using namespace sv; int main(int argc, char* argv[]) { @@ -18,16 +19,16 @@ int main(int argc, char* argv[]) mdl.setObjective(vars[0] + 2 * vars[1], MDL_MAXIMIZE); switch(mdl.optimize()) { - case Rtn::OPTIMAL: + case OPTIMAL: cout << "OPTIMAL SOLUTION: " << mdl.get(DoubleAttr::Obj) << endl; for (int i = 0; i < mdl.get(IntAttr::NumVars); i++) { cout << "var[" << i << "] : " << vars[i].get(DoubleAttr::X) << endl; } break; - case Rtn::INFEASIBLE: + case INFEASIBLE: cout << "INFEASIBLE MODEL" << endl; break; - case Rtn::UNBOUNDED: + case UNBOUNDED: cout << "UNBOUNDED SOLUTION" << endl; break; default: diff --git a/src/solver.cpp b/src/solver.cpp index b437bc9..4953777 100644 --- a/src/solver.cpp +++ b/src/solver.cpp @@ -1,148 +1,17 @@ #include "solver.hpp" +#include "common.hpp" + #include #include #include +using namespace sv; + using std::make_pair; using std::pair; using std::cout; using std::endl; - -Expr operator+(const Expr& x, const Expr& y) -{ - Expr exp; - exp.coeffs.resize(std::max(x.coeffs.size(), y.coeffs.size()), 0); - for (int c = 0; c < exp.coeffs.size(); c++) { - if (c < x.coeffs.size() && c < y.coeffs.size()) { - exp.coeffs[c] = x.coeffs[c] + y.coeffs[c]; - } - else { - exp.coeffs[c]= c < x.coeffs.size()? x.coeffs[c] : y.coeffs[c]; - } - } - exp.constant = x.constant + y.constant; - return exp; -} - -Expr operator+(const Expr& x) -{ - return x; -} - -Expr operator+(Var x, Var y) -{ - Expr exp; - exp.coeffs.resize(std::max(x.get(IntAttr::ModelColumn) + 1, y.get(IntAttr::ModelColumn) + 1), 0); - exp.coeffs[x.get(IntAttr::ModelColumn)] = x.get(DoubleAttr::Coeff); - exp.coeffs[y.get(IntAttr::ModelColumn)] = y.get(DoubleAttr::Coeff); - return exp; -} - -Expr operator+(Var x, double a) -{ - Expr exp; - exp.coeffs.resize(x.get(IntAttr::ModelColumn) + 1); - exp.constant = a; - return exp; -} - -Expr operator+(double a, Var x) -{ - return x + a; -} - -Expr operator-(const Expr& x, const Expr& y) -{ - Expr exp; - exp.coeffs.resize(std::max(x.coeffs.size(), y.coeffs.size()), 0); - for (int c = 0; c < exp.coeffs.size(); c++) { - if (c < x.coeffs.size() && c < y.coeffs.size()) { - exp.coeffs[c] = x.coeffs[c] - y.coeffs[c]; - } - else { - exp.coeffs[c] = c < x.coeffs.size() ? x.coeffs[c] : y.coeffs[c]; - } - } - exp.constant = x.constant + y.constant; - return exp; -} - -Expr operator-(const Expr& x) -{ - Expr expr(x); - for (int c = 0; c < expr.coeffs.size(); c++) { - expr.coeffs[c] = -expr.coeffs[c]; - } - expr.constant = -expr.constant; - return expr; -} - -Expr operator-(Var x) -{ - return -Expr(x); -} - -Expr operator-(Var x, Var y) -{ - Expr exp; - exp.coeffs.resize(std::max(x.get(IntAttr::ModelColumn) + 1, y.get(IntAttr::ModelColumn) + 1), 0); - exp.coeffs[x.get(IntAttr::ModelColumn)] = x.get(DoubleAttr::Coeff); - exp.coeffs[y.get(IntAttr::ModelColumn)] = -y.get(DoubleAttr::Coeff); - return exp; -} - -Expr operator-(Var x, double a) -{ - return x - Var(a); -} - -Expr operator-(double a, Var x) -{ - return x - a; -} - -Expr operator*(double a, Var x) -{ - Expr exp; - exp.coeffs.resize(x.get(IntAttr::ModelColumn) + 1, 0); - exp.coeffs[x.get(IntAttr::ModelColumn)] = a * x.get(DoubleAttr::Coeff); - return exp; -} - -Expr operator*(Var x, double a) -{ - return a * x; -} - -Expr operator*(const Expr& x, double a) -{ - Expr exp = x; - for (int c = 0; c < exp.coeffs.size(); c++) { - exp.coeffs[c] *= a; - } - exp.constant *= a; - return exp; -} - -Expr operator*(double a, const Expr& x) -{ - return x * a; -} - -Expr operator/(Var x, double a) -{ - return Expr(x) / a; -} - -Expr operator/(const Expr& x, double a) -{ - Expr exp = x; - for (int c = 0; c < exp.coeffs.size(); c++) { - exp.coeffs[c] /= a; - } - exp.constant /= a; - return exp; -} +using std::vector; Model::Model(): cn(0), @@ -168,15 +37,14 @@ Var* Model::addVars(int num) void Model::addConstr(const Expr& expr, ConstrOper sense, double rhs) { - if (sense == ConstrOper::LESS_EQUAL) { - matrix.push_back(vector(1, rhs - expr.constant)); - matrix.back().insert(matrix.back().end(), expr.coeffs.begin(), expr.coeffs.end()); + mt.push_back(vector(1, rhs - expr.constant)); + mt.back().insert(mt.back().end(), expr.coeffs.begin(), expr.coeffs.end()); } else if (sense == ConstrOper::GREATER_EQUAL) { - matrix.push_back(vector(1, expr.constant - rhs)); + mt.push_back(vector(1, expr.constant - rhs)); for (int coeff : expr.coeffs) { - matrix.back().push_back(-coeff); + mt.back().push_back(-coeff); } } else { @@ -184,8 +52,8 @@ void Model::addConstr(const Expr& expr, ConstrOper sense, double rhs) addConstr(expr, ConstrOper::GREATER_EQUAL, rhs); } - for (int c = matrix.back().size(); c <= cn; c++) { - matrix.back().push_back(0); + for (int c = mt.back().size(); c <= cn; c++) { + mt.back().push_back(0); } } @@ -193,47 +61,48 @@ void Model::setObjective(Expr obje, int _sense) { assert(_sense == 1 || _sense == -1); if (sense == 0) { - matrix.insert(matrix.begin(), obje.coeffs); - matrix.front().insert(matrix.front().begin(), -obje.constant); + mt.insert(mt.begin(), obje.coeffs); + mt.front().insert(mt.front().begin(), -obje.constant); for (int c = obje.coeffs.size() + 1; c <= cn; c++) { - matrix.front().push_back(0); + mt.front().push_back(0); } } else { - matrix.front().front() = -obje.constant; - for (int c = 0; c < obje.coeffs.size(); c++) { - matrix.front()[c + 1] = obje.coeffs[c]; + mt.front().front() = -obje.constant; + for (int c = 0; c < cn; c++) { + if (c < obje.coeffs.size()) { + mt.front()[c + 1] = obje.coeffs[c]; + } + else { + mt.front()[c] = 0; + } } - for (int c = obje.coeffs.size() + 1; c <= cn; c++) { - matrix.front()[c] = 0; - } - } sense = _sense; } void Model::print() { - for (size_t i = 0; i < matrix.size(); i++) { - for (size_t j = 0; j < matrix[0].size(); j++) { - cout << matrix[i][j] << "\t"; + for (size_t i = 0; i < mt.size(); i++) { + for (size_t j = 0; j < mt[0].size(); j++) { + cout << mt[i][j] << "\t"; } cout << endl; } } -Rtn Model::optimize() +rtn Model::optimize() { - bn = matrix.size(); + bn = mt.size(); for (int row = 1; row < bn; row++) { - matrix.front().push_back(0); + mt.front().push_back(0); for (int col = 1; col < bn; col++) { - matrix[row].push_back(col == row ? 1 : 0); + mt[row].push_back(col == row ? 1 : 0); } } - cn = matrix.front().size(); + cn = mt.front().size(); for (int i = 0; i < cn; i++) { - matrix.front()[i] = sense * matrix.front()[i]; + mt.front()[i] = sense * mt.front()[i]; } for (size_t i = 1; i < bn; i++) { @@ -243,7 +112,7 @@ Rtn Model::optimize() // === 判断初始解是否为可行解 === bool initial_feasible = true; for (int row = 1; row < bn; row++) { - if (matrix[row].front() < 0) { + if (mt[row].front() < 0) { initial_feasible = false; break; } @@ -251,29 +120,29 @@ Rtn Model::optimize() // === 构造初始可行解 === if (!initial_feasible) { - vector coeff = matrix.front(); - matrix.front() = vector(cn, .0); - matrix.front().push_back(1); + vector coeff = mt.front(); + mt.front() = vector(cn, .0); + mt.front().push_back(1); pair t = { -1 ,cn }; for (int row = 1; row < bn; row++) { - matrix[row].push_back(-1); - if (t.first == -1 || matrix[row].front() < matrix[t.first].front()) { + mt[row].push_back(-1); + if (t.first == -1 || mt[row].front() < mt[t.first].front()) { t.first = row; } } _gaussian(t); if (fabs(_simplex()) > 1e-10) { - rtn = Rtn::INFEASIBLE; - return rtn; + rtn_ = INFEASIBLE; + return rtn_; } // if the x0 in B, we should pivot it. auto iter = find(basic.begin(), basic.end(), cn); if (iter != basic.end()) { - for (int col = 1; col < matrix.front().size(); col++) { - if (fabs(matrix.front()[col]) > 1e-10) { + for (int col = 1; col < mt.front().size(); col++) { + if (fabs(mt.front()[col]) > 1e-10) { t = make_pair(iter - basic.begin() + 1, col); _gaussian(t); break; @@ -282,33 +151,33 @@ Rtn Model::optimize() } for (int row = 0; row < bn; row++) { - matrix[row].pop_back(); + mt[row].pop_back(); } // recover the coefficient line for (int col = 0; col < cn; col++) { - matrix.front()[col] = coeff[col]; + mt.front()[col] = coeff[col]; } for (int row = 1; row <= basic.size(); row++) { - int norm = matrix.front()[basic[row - 1]]; + int norm = mt.front()[basic[row - 1]]; for (int col = 0; col < cn; col++) { - matrix.front()[col] -= norm * matrix[row][col]; + mt.front()[col] -= norm * mt[row][col]; } } } - res = -sense * _simplex(); + res_ = -sense * _simplex(); for (int row = 1; row < bn; row++) { if (basic[row - 1] - 1 < cn - bn) { - vars[basic[row - 1] - 1].val = matrix[row].front(); + vars[basic[row - 1] - 1].val = mt[row].front(); } } - return rtn; + return rtn_; } double Model::get(DoubleAttr attr) { - return res; + return res_; } int Model::get(IntAttr attr) @@ -320,20 +189,20 @@ double Model::_simplex() { pair t; while (1) { - rtn = _pivot(t); - if (rtn == Rtn::OPTIMAL || rtn == Rtn::UNBOUNDED) { + rtn_ = _pivot(t); + if (rtn_ == OPTIMAL || rtn_ == UNBOUNDED) { break; } _gaussian(t); } - return matrix.front().front(); + return mt.front().front(); } -Rtn Model::_pivot(pair& p) +rtn Model::_pivot(pair& p) { p = make_pair(0, 0); double cmin = INT_MAX; - vector coef = matrix.front(); + vector coef = mt.front(); // === 非主轴元素中找最小值 === for (size_t i = 1; i < coef.size(); i++) { @@ -343,29 +212,29 @@ Rtn Model::_pivot(pair& p) } } if (cmin >= 0) { - return Rtn::OPTIMAL; + return OPTIMAL; } double bmin = INT_MAX; for (size_t row = 1; row < bn; row++) { - double tmp = matrix[row].front() / matrix[row][p.second]; - if (matrix[row][p.second] > 0 && bmin > tmp) { + double tmp = mt[row].front() / mt[row][p.second]; + if (mt[row][p.second] > 0 && bmin > tmp) { bmin = tmp; p.first = row; } } if (abs(bmin - INT_MAX) < 1e-10) { - return Rtn::UNBOUNDED; + return UNBOUNDED; } for (auto iter = basic.begin(); iter != basic.end(); iter++) { - if (matrix[p.first][*iter] != 0) { + if (mt[p.first][*iter] != 0) { *iter = p.second; break; } } assert(basic[p.first - 1] == p.second); - return Rtn::PIVOT; + return PIVOT; } void Model::_gaussian(pair p) @@ -373,9 +242,9 @@ void Model::_gaussian(pair p) size_t x = p.first, y = p.second; // === 主行归一化 === - double norm = matrix[x][y]; - for (size_t col = 0; col < matrix[x].size(); col++) { - matrix[x][col] /= norm; + double norm = mt[x][y]; + for (size_t col = 0; col < mt[x].size(); col++) { + mt[x][col] /= norm; } // === 其余行变换 === @@ -383,100 +252,13 @@ void Model::_gaussian(pair p) if (row == x) { continue; } - if (matrix[row][y] != 0) { - double norm = matrix[row][y]; - for (size_t col = 0; col < matrix[x].size(); col++) { - matrix[row][col] = matrix[row][col] - norm * matrix[x][col]; + if (mt[row][y] != 0) { + double norm = mt[row][y]; + for (size_t col = 0; col < mt[x].size(); col++) { + mt[row][col] = mt[row][col] - norm * mt[x][col]; } } } basic[x - 1] = y; // 换元 } - -Expr::Expr(double constant) - :constant(constant) -{ -} - -Expr::Expr(Var var, double coeff) -{ - this->coeffs.resize(var.col + 1); - this->coeffs[var.col] = coeff; - this->constant = 0; -} - -Expr Expr::operator=(const Expr& rhs) -{ - return *this; -} - -void Expr::operator+=(const Expr& expr) -{ - coeffs.resize(std::max(coeffs.size(), expr.coeffs.size())); - for (int c = 0; c < expr.coeffs.size(); c++) { - coeffs[c] += expr.coeffs[c]; - } - constant += expr.constant; -} - -void Expr::operator-=(const Expr& expr) -{ - coeffs.resize(std::max(coeffs.size(), expr.coeffs.size())); - for (int c = 0; c < expr.coeffs.size(); c++) { - coeffs[c] -= expr.coeffs[c]; - } - constant -= expr.constant; -} - -void Expr::operator*=(double mult) -{ - for (int c = 0; c < coeffs.size(); c++) { - coeffs[c] *= mult; - } - constant *= mult; -} - -void Expr::operator/=(double a) -{ - for (int c = 0; c < coeffs.size(); c++) { - coeffs[c] /= a; - } - constant /= a; -} - -Expr Expr::operator+(const Expr& rhs) -{ - coeffs.resize(std::max(coeffs.size(), rhs.coeffs.size())); - for (int c = 0; c < rhs.coeffs.size(); c++) { - coeffs[c] += rhs.coeffs[c]; - } - constant += rhs.constant; - return *this; -} - -Expr Expr::operator-(const Expr& rhs) -{ - coeffs.resize(std::max(coeffs.size(), rhs.coeffs.size())); - for (int c = 0; c < rhs.coeffs.size(); c++) { - coeffs[c] -= rhs.coeffs[c]; - } - constant -= rhs.constant; - return *this; -} - -double Var::get(DoubleAttr attr) -{ - switch (attr) { - case DoubleAttr::Coeff: - return coeffs; - case DoubleAttr::X: - return val; - } - return val; -} - -int Var::get(IntAttr attr) -{ - return col; -}