diff --git a/inc/solver.hpp b/inc/solver.hpp index 5fe27ed..f5f796a 100644 --- a/inc/solver.hpp +++ b/inc/solver.hpp @@ -31,8 +31,8 @@ namespace sv { Var* vars; - matrix mt; - matrix mt_cvt; + matrix table; + matrix ope_table; size_t cn, bn; std::vector basic; rtn rtn_; diff --git a/src/common.cpp b/src/common.cpp index 87bace8..b8b02fe 100644 --- a/src/common.cpp +++ b/src/common.cpp @@ -33,10 +33,10 @@ Expr sv::operator+(const Expr& x, const Expr& y) 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]; + exp.coeffs.at(c) = x.coeffs.at(c) + y.coeffs.at(c); } else { - exp.coeffs[c] = c < x.coeffs.size() ? x.coeffs[c] : y.coeffs[c]; + exp.coeffs.at(c) = c < x.coeffs.size() ? x.coeffs.at(c) : y.coeffs.at(c); } } exp.constant = x.constant + y.constant; @@ -52,8 +52,8 @@ 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); + exp.coeffs.at(x.get(IntAttr::NumCol)) = x.get(DoubleAttr::Coeff); + exp.coeffs.at(y.get(IntAttr::NumCol)) = y.get(DoubleAttr::Coeff); return exp; } @@ -76,10 +76,10 @@ Expr sv::operator-(const Expr& x, const Expr& y) 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]; + exp.coeffs.at(c) = x.coeffs.at(c) - y.coeffs.at(c); } else { - exp.coeffs[c] = c < x.coeffs.size() ? x.coeffs[c] : y.coeffs[c]; + exp.coeffs.at(c) = c < x.coeffs.size() ? x.coeffs.at(c) : y.coeffs.at(c); } } exp.constant = x.constant + y.constant; @@ -90,7 +90,7 @@ 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.coeffs.at(c) = -expr.coeffs.at(c); } expr.constant = -expr.constant; return expr; @@ -105,8 +105,8 @@ 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); + exp.coeffs.at(x.get(IntAttr::NumCol)) = x.get(DoubleAttr::Coeff); + exp.coeffs.at(y.get(IntAttr::NumCol)) = -y.get(DoubleAttr::Coeff); return exp; } @@ -124,7 +124,7 @@ 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); + exp.coeffs.at(x.get(IntAttr::NumCol)) = a * x.get(DoubleAttr::Coeff); return exp; } @@ -137,7 +137,7 @@ 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.coeffs.at(c) *= a; } exp.constant *= a; return exp; @@ -157,7 +157,7 @@ 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.coeffs.at(c) /= a; } exp.constant /= a; return exp; @@ -171,7 +171,7 @@ Expr::Expr(double constant) Expr::Expr(Var var, double coeff) { this->coeffs.resize(var.col + 1); - this->coeffs[var.col] = coeff; + this->coeffs.at(var.col) = coeff; this->constant = 0; } @@ -184,7 +184,7 @@ 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]; + coeffs.at(c) += expr.coeffs.at(c); } constant += expr.constant; } @@ -193,7 +193,7 @@ 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]; + coeffs.at(c) -= expr.coeffs.at(c); } constant -= expr.constant; } @@ -201,7 +201,7 @@ void Expr::operator-=(const Expr& expr) void Expr::operator*=(double mult) { for (int c = 0; c < coeffs.size(); c++) { - coeffs[c] *= mult; + coeffs.at(c) *= mult; } constant *= mult; } @@ -209,7 +209,7 @@ void Expr::operator*=(double mult) void Expr::operator/=(double a) { for (int c = 0; c < coeffs.size(); c++) { - coeffs[c] /= a; + coeffs.at(c) /= a; } constant /= a; } @@ -218,7 +218,7 @@ 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]; + coeffs.at(c) += rhs.coeffs.at(c); } constant += rhs.constant; return *this; @@ -228,7 +228,7 @@ 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]; + coeffs.at(c) -= rhs.coeffs.at(c); } constant -= rhs.constant; return *this; diff --git a/src/main.cpp b/src/main.cpp index 72cb57e..da91f48 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -10,20 +10,14 @@ using namespace sv; int main(int argc, char* argv[]) { Model mdl; - - Var* vars = mdl.addVars(2, VarType::INTEGER); - - - //mdl.addConstr(vars[0] + vars[1], ConstrOper::LESS_EQUAL, 2); - //mdl.addConstr(vars[0] + vars[1], ConstrOper::GREATER_EQUAL, 1); - //mdl.addConstr(5 * vars[0] - 2 * vars[1], ConstrOper::GREATER_EQUAL, -2); - - //mdl.setObjective(vars[0] + 2 * vars[1], MDL_MAXIMIZE); + + Var* vars = mdl.addVars(3, VarType::INTEGER); mdl.addConstr(2 * vars[0] + vars[1], ConstrOper::LESS_EQUAL, 10); mdl.addConstr(3 * vars[0] + 6 * vars[1], ConstrOper::LESS_EQUAL, 40); + mdl.addConstr(3 * vars[0] + 6 * vars[1] + 4 * vars[2], ConstrOper::LESS_EQUAL, 50); + mdl.setObjective(100 * vars[0] + 150 * vars[1] + 120 * vars[2], MDL_MAXIMIZE); - mdl.setObjective(100 * vars[0] + 150 * vars[1], MDL_MAXIMIZE); switch(mdl.optimize()) { case OPTIMAL: cout << "OPTIMAL SOLUTION: " << mdl.get(DoubleAttr::Obj) << endl; diff --git a/src/solver.cpp b/src/solver.cpp index b44ca18..971270a 100644 --- a/src/solver.cpp +++ b/src/solver.cpp @@ -24,11 +24,14 @@ struct Node }; LinSolver::LinSolver() : + obj_(0), + rtn_(LOADED), cn(0), bn(1), sense(0), vars(nullptr) { + } sv::LinSolver::~LinSolver() @@ -40,7 +43,7 @@ sv::LinSolver::LinSolver(const LinSolver& solver) : vars(nullptr), cn(solver.cn), bn(solver.bn), - mt(solver.mt), + table(solver.table), basic(solver.basic), rtn_(solver.rtn_), obj_(solver.obj_), @@ -71,8 +74,7 @@ LinSolver& sv::LinSolver::operator=(const LinSolver& solver) vars[i] = solver.vars[i]; // 深拷贝 } } - - mt = solver.mt; + table = solver.table; cn = solver.cn, bn = solver.bn; basic = solver.basic; rtn_ = solver.rtn_; @@ -91,9 +93,6 @@ Var* LinSolver::addVars(int num, VarType type) vars[c].type = type; } - //for (int i = 0; i < cn; i++) { - // *(vars + i) = *(old + i); - //} memcpy(vars, old, sizeof(Var) * cn); delete[] old; cn += num; @@ -110,14 +109,14 @@ void LinSolver::addConstr(const Expr& expr, ConstrOper sense, double rhs) { if (sense == ConstrOper::LESS_EQUAL) { bn++; - mt.push_back(vector(1, rhs - expr.constant)); - mt.back().insert(mt.back().end(), expr.coeffs.begin(), expr.coeffs.end()); + table.push_back(vector(1, rhs - expr.constant)); + table.back().insert(table.back().end(), expr.coeffs.begin(), expr.coeffs.end()); } else if (sense == ConstrOper::GREATER_EQUAL) { bn++; - mt.push_back(vector(1, expr.constant - rhs)); + table.push_back(vector(1, expr.constant - rhs)); for (int coeff : expr.coeffs) { - mt.back().push_back(-coeff); + table.back().push_back(-coeff); } } else { @@ -125,8 +124,8 @@ void LinSolver::addConstr(const Expr& expr, ConstrOper sense, double rhs) addConstr(expr, ConstrOper::GREATER_EQUAL, rhs); } - for (int c = mt.back().size(); c <= cn; c++) { - mt.back().push_back(0); + for (int c = table.back().size(); c <= cn; c++) { + table.back().push_back(0); } } @@ -134,25 +133,25 @@ void LinSolver::setObjective(Expr obje, int _sense) { assert(_sense == 1 || _sense == -1); if (sense == 0) { - mt.insert(mt.begin(), obje.coeffs); - mt.front().insert(mt.front().begin(), -obje.constant); + table.insert(table.begin(), obje.coeffs); + table.front().insert(table.front().begin(), -obje.constant); for (int c = obje.coeffs.size() + 1; c <= cn; c++) { - mt.front().push_back(0); + table.front().push_back(0); } } else { - mt.front().front() = -obje.constant; - for (int c = 0; c < cn; c++) { - if (c < obje.coeffs.size()) { - mt.front()[c + 1] = obje.coeffs[c]; + table.front().front() = -obje.constant; + for (int col = 0; col < cn; col++) { + if (col < obje.coeffs.size()) { + table.front().at(col + 1) = obje.coeffs.at(col); } else { - mt.front()[c] = 0; + table.front().at(col) = 0; } } } - for (int i = 0; i < mt.front().size(); i++) { - mt.front()[i] = _sense * mt.front()[i]; + for (int row = 0; row < table.front().size(); row++) { + table.front().at(row) = _sense * table.front().at(row); } sense = _sense; } @@ -160,27 +159,29 @@ void LinSolver::setObjective(Expr obje, int _sense) rtn LinSolver::optimize() { assert(sense); - mt_cvt = mt; + ope_table = table; + rtn_ = LOADED; rtn_ = feasible_solution(); if (rtn_ == LOADED) { obj_ = _simplex(); } - cn = mt_cvt.front().size() - bn; - for (int row = 1; row < bn; row++) { - if (basic[row - 1] - 1 <= cn) { - vars[basic[row - 1] - 1].val = mt_cvt[row].front(); + if (rtn_ == OPTIMAL) { + cn = ope_table.front().size() - bn; + for (int row = 1; row < bn; row++) { + if (basic.at(row - 1) - 1 < cn) { + vars[basic.at(row - 1) - 1].val = ope_table.at(row).front(); + } } } - return rtn_; } void LinSolver::print() { - for (size_t i = 0; i < mt_cvt.size(); i++) { - for (size_t j = 0; j < mt_cvt[0].size(); j++) { - cout << mt_cvt[i][j] << "\t"; + for (size_t row = 0; row < ope_table.size(); row++) { + for (size_t col = 0; col < ope_table.front().size(); col++) { + cout << ope_table.at(row).at(col) << "\t"; } cout << endl; } @@ -202,14 +203,11 @@ rtn Model::optimize() std::stack list_; - Node incumbent_node, root_node; - root_node.solver = solver; - root_node.upper_bound = global_upper_bound, root_node.lower_bound = global_lower_bound; + Node root_node(solver, 0, solver.obj_); + Node incumbent_node = root_node; list_.push(root_node); - int cnt = 0; while (list_.size() && global_upper_bound - global_lower_bound > 1e-10) { - cout << ++cnt << endl; Node current_node = list_.top(); list_.pop(); current_node.solver.optimize(); @@ -310,18 +308,18 @@ double LinSolver::_simplex() } _gaussian(t); } - return obj_ = mt_cvt.front().front(); + return obj_ = ope_table.front().front(); } rtn LinSolver::feasible_solution() { for (int row = 1; row < bn; row++) { - mt_cvt.front().push_back(0); + ope_table.front().push_back(0); for (int col = 1; col < bn; col++) { - mt_cvt[row].push_back(col == row ? 1 : 0); + ope_table.at(row).push_back(col == row ? 1 : 0); } } - cn = mt_cvt.front().size(); + cn = ope_table.front().size(); basic.clear(); for (size_t i = 1; i < bn; i++) { basic.push_back(cn - bn + i); @@ -330,7 +328,7 @@ rtn LinSolver::feasible_solution() // === 判断初始解是否为可行解 === bool initial_feasible = true; for (int row = 1; row < bn; row++) { - if (mt_cvt[row].front() < 0) { + if (ope_table.at(row).front() < 0) { initial_feasible = false; break; } @@ -338,65 +336,64 @@ rtn LinSolver::feasible_solution() // === 构造初始可行解 === if (!initial_feasible) { - vector coeff = mt_cvt.front(); - mt_cvt.front() = vector(cn, .0); - mt_cvt.front().push_back(1); + vector coeff = ope_table.front(); + ope_table.front() = vector(cn, .0); + ope_table.front().push_back(1); pair t = { -1 ,cn }; for (int row = 1; row < bn; row++) { - mt_cvt[row].push_back(-1); - if (t.first == -1 || mt_cvt[row].front() < mt_cvt[t.first].front()) { + ope_table.at(row).push_back(-1); + if (t.first == -1 || ope_table.at(row).front() < ope_table.at(t.first).front()) { t.first = row; } } _gaussian(t); if (fabs(_simplex()) > 1e-10) { - rtn_ = INFEASIBLE; - return rtn_; + return rtn_ = INFEASIBLE; } - + rtn_ = LOADED; // 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 < mt_cvt.front().size(); col++) { - if (fabs(mt_cvt.front()[col]) > 1e-10) { + for (int col = 1; col < ope_table.front().size(); col++) { + if (fabs(ope_table.front().at(col)) > 1e-10) { t = make_pair(iter - basic.begin() + 1, col); _gaussian(t); break; } } } - + for (int row = 0; row < bn; row++) { - mt_cvt[row].pop_back(); + ope_table.at(row).pop_back(); } // recover the coefficient line for (int col = 0; col < cn; col++) { - mt_cvt.front()[col] = coeff[col]; + ope_table.front().at(col) = coeff.at(col); } for (int row = 1; row <= basic.size(); row++) { - int norm = mt_cvt.front()[basic[row - 1]]; + int norm = ope_table.front().at(basic.at(row - 1)); for (int col = 0; col < cn; col++) { - mt_cvt.front()[col] -= norm * mt_cvt[row][col]; + ope_table.front().at(col) -= norm * ope_table.at(row).at(col); } } } - return LOADED; + return rtn_; } rtn LinSolver::_pivot(pair& p) { p = make_pair(0, 0); double cmin = INT_MAX; - vector coef = mt_cvt.front(); + vector coef = ope_table.front(); // === 非主轴元素中找最小值 === - for (size_t i = 1; i < coef.size(); i++) { - if (cmin > coef[i] && find(basic.begin(), basic.end(), i) == basic.end()) { - cmin = coef[i]; - p.second = i; + for (size_t col = 1; col < coef.size(); col++) { + if (cmin > coef.at(col) && find(basic.begin(), basic.end(), col) == basic.end()) { + cmin = coef.at(col); + p.second = col; } } if (cmin >= 0) { @@ -404,8 +401,8 @@ rtn LinSolver::_pivot(pair& p) } double bmin = INT_MAX; for (size_t row = 1; row < bn; row++) { - double tmp = mt_cvt[row].front() / mt_cvt[row][p.second]; - if (mt_cvt[row][p.second] > 0 && bmin > tmp) { + double tmp = ope_table.at(row).front() / ope_table.at(row).at(p.second); + if (ope_table.at(row).at(p.second) > 0 && bmin > tmp) { bmin = tmp; p.first = row; } @@ -416,12 +413,12 @@ rtn LinSolver::_pivot(pair& p) } for (auto iter = basic.begin(); iter != basic.end(); iter++) { - if (mt_cvt[p.first][*iter] != 0) { + if (ope_table.at(p.first).at(*iter) != 0) { *iter = p.second; break; } } - assert(basic[p.first - 1] == p.second); + assert(basic.at(p.first - 1) == p.second); return PIVOT; } @@ -430,9 +427,9 @@ void LinSolver::_gaussian(pair p) size_t x = p.first, y = p.second; // === 主行归一化 === - double norm = mt_cvt[x][y]; - for (size_t col = 0; col < mt_cvt[x].size(); col++) { - mt_cvt[x][col] /= norm; + double norm = ope_table.at(x).at(y); + for (size_t col = 0; col < ope_table.at(x).size(); col++) { + ope_table.at(x).at(col) /= norm; } // === 其余行变换 === @@ -440,13 +437,13 @@ void LinSolver::_gaussian(pair p) if (row == x) { continue; } - if (mt_cvt[row][y] != 0) { - double norm = mt_cvt[row][y]; - for (size_t col = 0; col < mt_cvt[x].size(); col++) { - mt_cvt[row][col] = mt_cvt[row][col] - norm * mt_cvt[x][col]; + if (ope_table.at(row).at(y) != 0) { + double norm = ope_table.at(row).at(y); + for (size_t col = 0; col < ope_table.at(x).size(); col++) { + ope_table.at(row).at(col) = ope_table.at(row).at(col) - norm * ope_table.at(x).at(col); } } } - basic[x - 1] = y; // 换元 + basic.at(x - 1) = y; // 换元 }