调整表达式运算错误，重新排版代码，优化单纯形搜索过程

CMakeLists.txt

@@ -1,6 +1,6 @@
 cmake_minimum_required(VERSION 3.16.0)
 
-project(Solver)
+project(solver)
 
 set(CMAKE_CXX_STANDARD 20)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
@@ -22,10 +22,10 @@ file(GLOB HEADERS "${INCLUDE_DIR}/*.hpp")
 aux_source_directory(${SOURCE_DIR} src)
 aux_source_directory(${INCLUDE_DIR} inc)
 
-add_executable(Solver ${SOURCES} ${HEADERS})
+add_executable(solver ${SOURCES} ${HEADERS})
 
-#target_link_libraries(Solver 
+#target_link_libraries(solver 
 #  ${GUROBI_LIBRARIES}
 #)
 
-target_include_directories(Solver PRIVATE ${SOURCE_DIR} ${INCLUDE_DIR})
+target_include_directories(solver PRIVATE ${SOURCE_DIR} ${INCLUDE_DIR})

TodoList.md

@@ -0,0 +1,19 @@
+## P0 基础修复
+
+1. 补齐变量约束能力：完善变量 LB/UB 上下界定义与校验；对 BINARY 类型严格强制 0-1 取值
+2. 增加测试与回归用例：搭建基础测试基线，覆盖整数、二进制、简单 MIP 算例，建立回归机制
+
+## P1 核心能力
+
+3. 完善终止条件与工程可控性：实现 TimeLimit 时间限制、MIPGap 收敛间隙终止条件；完善运行日志输出
+4. 迭代 MIP 分支与节点选择策略：按顺序落地：best-bound 选节点 → pseudocost 分支策略
+
+## P2 进阶算法
+
+5. 增加启发式、割平面与预处理能力：实现一种rounding舍入启发式、root层Gomory割平面、简单预处理逻辑
+
+## P3 底层求解架构
+
+6. 学习稀疏单纯形，对接HiGHS等成熟求解器架构：调研学习稀疏单纯形、稀疏基分解、对偶单纯形热启动机制
+
+## P4 工程收口与标准化（待定）

inc/common.hpp

@@ -24,12 +24,12 @@
 #define PIVOT 17
 
 namespace sv {
-	using matrix = std::vector<std::vector<double>>;
-	using rtn = int;
-	class Expr;
-	class Var;
+using matrix = std::vector<std::vector<double>>;
+using rtn = int;
+class Expr;
+class Var;
 
-	enum class IntAttr {
+enum class IntAttr {
 	NumConstrs,
 	NumVars,
 	NumIntVars,
@@ -42,9 +42,9 @@ namespace sv {
 	Lazy,
 	NumObj,
 	NumCol
-	};
+};
 
-	enum class DoubleAttr {
+enum class DoubleAttr {
 	Runtime,
 	Work,
 	ObjCon,
@@ -64,9 +64,9 @@ namespace sv {
 	NodeCount,
 	X,
 	Slack,
-	};
+};
 
-	enum class StringAttr {
+enum class StringAttr {
 	ModelName,
 	VarName,
 	ConstrName,
@@ -79,23 +79,18 @@ namespace sv {
 	CTag,
 	QCTag,
 	BatchErrorMessage
-	};
+};
 
-	enum class ConstrOper {
-		LESS_EQUAL,
-		GREATER_EQUAL,
-		EQUAL
-	};
+enum class ConstrOper { LESS_EQUAL, GREATER_EQUAL, EQUAL };
 
-
-	enum class VarType {
+enum class VarType {
 	CONTINUOUS,
 	BINARY,
 	INTEGER,
-	};
+};
 
-	class Var {
-	public:
+class Var {
+public:
 	friend class Model;
 	friend class LinSolver;
 	friend class Expr;
@@ -104,71 +99,68 @@ namespace sv {
 	double get(DoubleAttr attr);
 	int get(IntAttr attr);
 
-	private:
+private:
 	double coeffs;
 	double val;
 	int col;
 	VarType type;
-	};
+};
 
-	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);
+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:
+class Expr {
+private:
 	double constant;
 	std::vector<double> coeffs;
 	std::vector<Var> vars;
 
-	public:
-		Expr(const Expr& expr) = default;
+public:
+	Expr(const Expr &expr) = default;
 	Expr(double constant = 0.0);
 	Expr(Var var, double coeff = 1.0);
 
 	friend class LinSolver;
 
-		friend Expr operator+(const Expr& x, const Expr& y);
-		friend Expr operator+(const Expr& x);
+	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-(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*(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);
+	friend Expr operator/(const Expr &x, double a);
 
-
-		Expr operator=(const Expr& rhs);
-		void operator+=(const Expr& expr);
-		void operator-=(const Expr& expr);
+	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);
-	};
-}
+	Expr operator+(const Expr &rhs);
+	Expr operator-(const Expr &rhs);
+};
+} // namespace sv

inc/solver.hpp

@@ -3,18 +3,18 @@
 
 namespace sv {
 
-	class LinSolver {
-	public:
+class LinSolver {
+public:
 	friend class Model;
 	LinSolver();
 	~LinSolver();
 
-		LinSolver(const LinSolver& solver);
-		LinSolver& operator=(const LinSolver& solver);
+	LinSolver(const LinSolver &solver);
+	LinSolver &operator=(const LinSolver &solver);
 
-		Var* addVars(int col, VarType type);
-		Var& getVar(int idx);
-		void addConstr(const Expr& expr, ConstrOper sense, double rhs);
+	Var *addVars(int col, VarType type);
+	Var &getVar(int idx);
+	void addConstr(const Expr &expr, ConstrOper sense, double rhs);
 	void setObjective(Expr obje, int sense = MDL_MAXIMIZE);
 	void print();
 
@@ -22,14 +22,14 @@ namespace sv {
 	int get(IntAttr attr);
 
 	rtn optimize();
-	protected:
 
+protected:
 	double _simplex();
-		rtn _pivot(std::pair<size_t, size_t>& p);
+	rtn _pivot(std::pair<size_t, size_t> &p);
 	rtn feasible_solution();
 	void _gaussian(std::pair<size_t, size_t> p);
 
-		std::vector<std::pair<int, Var*>> vars;
+	std::vector<std::pair<int, Var *>> vars;
 
 	matrix table;
 	matrix ope_table;
@@ -38,25 +38,21 @@ namespace sv {
 	rtn rtn_;
 	double obj_;
 	int sense;
-	};
+};
 
-	class Model
-	{
-	public:
+class Model {
+public:
 	Model();
 	rtn optimize();
 
-		Var* addVars(int col, VarType type);
-		void addConstr(const Expr& expr, ConstrOper sense, double rhs);
+	Var *addVars(int col, VarType type);
+	void addConstr(const Expr &expr, ConstrOper sense, double rhs);
 	void setObjective(Expr obje, int sense = MDL_MAXIMIZE);
 
 	double get(DoubleAttr attr);
 	int get(IntAttr attr);
 
-	private:
+private:
 	LinSolver solver;
-	};
-}
-
-
-
+};
+} // namespace sv

src/common.cpp

@@ -2,17 +2,9 @@
 #include <iostream>
 using namespace sv;
 
+Var::Var(double coef, VarType type_) : col(0), val(0), coeffs(coef), type(type_){};
 
-Var::Var(double coef, VarType type_) :
-	col(0), 
-	val(0), 
-	coeffs(coef), 
-	type(type_) 
-{
-};
-
-double Var::get(DoubleAttr attr)
-{
+double Var::get(DoubleAttr attr) {
 	switch (attr) {
 	case DoubleAttr::Coeff:
 		return coeffs;
@@ -22,13 +14,11 @@ double Var::get(DoubleAttr attr)
 	return -1;
 }
 
-int Var::get(IntAttr attr)
-{
+int Var::get(IntAttr attr) {
 	return col;
 }
 
-Expr sv::operator+(const Expr& x, const Expr& y)
-{
+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++) {
@@ -43,13 +33,11 @@ Expr sv::operator+(const Expr& x, const Expr& y)
 	return exp;
 }
 
-Expr sv::operator+(const Expr& x)
-{
+Expr sv::operator+(const Expr &x) {
 	return x;
 }
 
-Expr sv::operator+(Var x, Var y)
-{
+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.at(x.get(IntAttr::NumCol)) = x.get(DoubleAttr::Coeff);
@@ -57,21 +45,19 @@ Expr sv::operator+(Var x, Var y)
 	return exp;
 }
 
-Expr sv::operator+(Var x, double a)
-{
+Expr sv::operator+(Var x, double a) {
 	Expr exp;
-	exp.coeffs.resize(x.get(IntAttr::NumCol) + 1);
+	exp.coeffs.resize(x.get(IntAttr::NumCol) + 1, 0);
+	exp.coeffs.at(x.get(IntAttr::NumCol)) = x.get(DoubleAttr::Coeff);
 	exp.constant = a;
 	return exp;
 }
 
-Expr sv::operator+(double a, Var x)
-{
+Expr sv::operator+(double a, Var x) {
 	return x + a;
 }
 
-Expr sv::operator-(const Expr& x, const Expr& y)
-{
+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++) {
@@ -79,15 +65,14 @@ Expr sv::operator-(const Expr& x, const Expr& y)
 			exp.coeffs.at(c) = x.coeffs.at(c) - y.coeffs.at(c);
 		}
 		else {
-			exp.coeffs.at(c) = c < x.coeffs.size() ? x.coeffs.at(c) : y.coeffs.at(c);
+			exp.coeffs.at(c) = c < x.coeffs.size() ? x.coeffs.at(c) : -y.coeffs.at(c);
 		}
 	}
-	exp.constant = x.constant + y.constant;
+	exp.constant = x.constant - y.constant;
 	return exp;
 }
 
-Expr sv::operator-(const Expr& x)
-{
+Expr sv::operator-(const Expr &x) {
 	Expr expr(x);
 	for (int c = 0; c < expr.coeffs.size(); c++) {
 		expr.coeffs.at(c) = -expr.coeffs.at(c);
@@ -96,13 +81,11 @@ Expr sv::operator-(const Expr& x)
 	return expr;
 }
 
-Expr sv::operator-(Var x)
-{
+Expr sv::operator-(Var x) {
 	return -Expr(x);
 }
 
-Expr sv::operator-(Var x, Var y)
-{
+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.at(x.get(IntAttr::NumCol)) = x.get(DoubleAttr::Coeff);
@@ -110,31 +93,30 @@ Expr sv::operator-(Var x, Var y)
 	return exp;
 }
 
-Expr sv::operator-(Var x, double a)
-{
-	return x - Var(a);
+Expr sv::operator-(Var x, double a) {
+	Expr exp;
+	exp.coeffs.resize(x.get(IntAttr::NumCol) + 1, 0);
+	exp.coeffs.at(x.get(IntAttr::NumCol)) = x.get(DoubleAttr::Coeff);
+	exp.constant = -a;
+	return exp;
 }
 
-Expr sv::operator-(double a, Var x)
-{
-	return x - a;
+Expr sv::operator-(double a, Var x) {
+	return a + (-x);
 }
 
-Expr sv::operator*(double a, Var x)
-{
+Expr sv::operator*(double a, Var x) {
 	Expr exp;
 	exp.coeffs.resize(x.get(IntAttr::NumCol) + 1, 0);
 	exp.coeffs.at(x.get(IntAttr::NumCol)) = a * x.get(DoubleAttr::Coeff);
 	return exp;
 }
 
-Expr sv::operator*(Var x, double a)
-{
+Expr sv::operator*(Var x, double a) {
 	return a * x;
 }
 
-Expr sv::operator*(const Expr& x, double a)
-{
+Expr sv::operator*(const Expr &x, double a) {
 	Expr exp = x;
 	for (int c = 0; c < exp.coeffs.size(); c++) {
 		exp.coeffs.at(c) *= a;
@@ -143,18 +125,15 @@ Expr sv::operator*(const Expr& x, double a)
 	return exp;
 }
 
-Expr sv::operator*(double a, const Expr& x)
-{
+Expr sv::operator*(double a, const Expr &x) {
 	return x * a;
 }
 
-Expr sv::operator/(Var x, double a)
-{
+Expr sv::operator/(Var x, double a) {
 	return Expr(x) / a;
 }
 
-Expr sv::operator/(const Expr& x, double a)
-{
+Expr sv::operator/(const Expr &x, double a) {
 	Expr exp = x;
 	for (int c = 0; c < exp.coeffs.size(); c++) {
 		exp.coeffs.at(c) /= a;
@@ -163,25 +142,24 @@ Expr sv::operator/(const Expr& x, double a)
 	return exp;
 }
 
-Expr::Expr(double constant)
-	:constant(constant)
-{
-}
+Expr::Expr(double constant) : constant(constant) {}
 
-Expr::Expr(Var var, double coeff)
-{
+Expr::Expr(Var var, double coeff) {
 	this->coeffs.resize(var.col + 1);
 	this->coeffs.at(var.col) = coeff;
 	this->constant = 0;
 }
 
-Expr Expr::operator=(const Expr& rhs)
-{
+Expr &Expr::operator=(const Expr &rhs) {
+	if (this != &rhs) {
+		constant = rhs.constant;
+		coeffs = rhs.coeffs;
+		vars = rhs.vars;
+	}
 	return *this;
 }
 
-void Expr::operator+=(const Expr& expr)
-{
+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.at(c) += expr.coeffs.at(c);
@@ -189,8 +167,7 @@ void Expr::operator+=(const Expr& expr)
 	constant += expr.constant;
 }
 
-void Expr::operator-=(const Expr& expr)
-{
+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.at(c) -= expr.coeffs.at(c);
@@ -198,24 +175,21 @@ void Expr::operator-=(const Expr& expr)
 	constant -= expr.constant;
 }
 
-void Expr::operator*=(double mult)
-{
+void Expr::operator*=(double mult) {
 	for (int c = 0; c < coeffs.size(); c++) {
 		coeffs.at(c) *= mult;
 	}
 	constant *= mult;
 }
 
-void Expr::operator/=(double a)
-{
+void Expr::operator/=(double a) {
 	for (int c = 0; c < coeffs.size(); c++) {
 		coeffs.at(c) /= a;
 	}
 	constant /= a;
 }
 
-Expr Expr::operator+(const Expr& rhs)
-{
+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.at(c) += rhs.coeffs.at(c);
@@ -224,8 +198,7 @@ Expr Expr::operator+(const Expr& rhs)
 	return *this;
 }
 
-Expr Expr::operator-(const Expr& rhs)
-{
+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.at(c) -= rhs.coeffs.at(c);

src/main.cpp

@@ -6,29 +6,28 @@
 using namespace std;
 using namespace sv;
 
-
-int main(int argc, char* argv[])
-{
+int main(int argc, char *argv[]) {
 	Model mdl;
 
-	Var* int_var = mdl.addVars(3, VarType::INTEGER);
-	Var* con_var = mdl.addVars(2, VarType::CONTINUOUS);
+	Var *x = mdl.addVars(3, VarType::INTEGER);
+	Var *y = mdl.addVars(2, VarType::CONTINUOUS);
 
-	mdl.addConstr(2 * int_var[0] + int_var[1], ConstrOper::LESS_EQUAL, 10);
-	mdl.addConstr(3 * int_var[0] + 6 * int_var[1], ConstrOper::LESS_EQUAL, 40);
-	mdl.addConstr(3 * int_var[0] + 6 * int_var[1] + 4 * int_var[2], ConstrOper::LESS_EQUAL, 50);
-	mdl.addConstr(2.3 * con_var[0] + 2.6 * con_var[1], ConstrOper::LESS_EQUAL, 80);
-	mdl.addConstr(con_var[0] + 2 * con_var[1], ConstrOper::LESS_EQUAL, 70);
-	mdl.setObjective(100 * int_var[0] + 150 * int_var[1] + 120 * int_var[2] + 82.6 * con_var[0] + 90.4 * con_var[1], MDL_MAXIMIZE);
+	mdl.addConstr(3 * x[0] + 6 * x[1], ConstrOper::LESS_EQUAL, 28);
+	mdl.addConstr(3 * x[0] + 6 * x[1] + 4 * x[2], ConstrOper::LESS_EQUAL, 30);
+	mdl.addConstr(3 * y[0] + 2 * y[1], ConstrOper::LESS_EQUAL, 80);
+	mdl.addConstr(y[0] + 2 * y[1], ConstrOper::LESS_EQUAL, 70);
+	mdl.addConstr(2 * x[1] + x[2], ConstrOper::EQUAL, 7);
+	mdl.setObjective(100 * x[0] + 150 * x[1] + 120 * x[2] + 80 * y[0] + 90 * y[1] + 15,
+	                 MDL_MAXIMIZE);
 
-	switch(mdl.optimize()) {
+	switch (mdl.optimize()) {
 	case OPTIMAL:
 		cout << "OPTIMAL SOLUTION: " << mdl.get(DoubleAttr::Obj) << endl;
 		for (int i = 0; i < 3; i++) {
-			cout << "integer var [" << i << "] : " << int_var[i].get(DoubleAttr::X) << endl;
+			cout << "integer var [" << i << "] : " << x[i].get(DoubleAttr::X) << endl;
 		}
 		for (int i = 0; i < 2; i++) {
-			cout << "continuous var [" << i << "] : " << con_var[i].get(DoubleAttr::X) << endl;
+			cout << "continuous var [" << i << "] : " << y[i].get(DoubleAttr::X) << endl;
 		}
 		break;
 	case INFEASIBLE:

src/solver.cpp

@@ -5,68 +5,53 @@
 #include <cassert>
 #include <algorithm>
 #include <stack>
-#include <cassert>
-#include <cstring>
+#include <cmath>
 #include <cfloat>
-#include <math.h>
+#include <limits>
+
 using namespace sv;
 
-using std::make_pair;
-using std::pair;
 using std::cout;
 using std::endl;
+using std::make_pair;
+using std::pair;
 using std::vector;
 
-
-struct Node
-{
+struct Node {
 	LinSolver solver;
-	double lower_bound;
-	double upper_bound;
-
+	double bound; // LP bound in true-objective space
 };
 
-LinSolver::LinSolver() :
-	obj_(0),
-	rtn_(LOADED),
-	cn(0),
-	bn(1),
-	sense(0)
-{
+LinSolver::LinSolver() : obj_(0), rtn_(LOADED), cn(0), bn(1), sense(0) {}
 
-}
-
-sv::LinSolver::~LinSolver()
-{
-	for (auto& var : vars) {
+sv::LinSolver::~LinSolver() {
+	for (auto &var : vars) {
 		delete[] var.second;
 	}
 }
 
-sv::LinSolver::LinSolver(const LinSolver& solver)
-{
+sv::LinSolver::LinSolver(const LinSolver &solver) {
 	*this = solver;
 }
 
-LinSolver& sv::LinSolver::operator=(const LinSolver& solver)
-{
+LinSolver &sv::LinSolver::operator=(const LinSolver &solver) {
 	if (this == &solver) {
 		return *this;
 	}
-	for (auto& var : vars) {
+	for (auto &var : vars) {
 		delete[] var.second;
 	}
 	vars.clear();
 	vars.reserve(solver.vars.size());
-	for (auto& var : solver.vars) {
+	for (auto &var : solver.vars) {
 		vars.push_back(std::make_pair(var.first, new Var[var.first]));
 		for (int i = 0; i < var.first; i++) {
 			vars.back().second[i] = var.second[i];
 		}
 	}
 	cn = solver.cn;
+	bn = solver.bn;
 	table = solver.table;
-	cn = solver.cn, bn = solver.bn;
 	basic = solver.basic;
 	rtn_ = solver.rtn_;
 	obj_ = solver.obj_;
@@ -75,85 +60,73 @@ LinSolver& sv::LinSolver::operator=(const LinSolver& solver)
 	return *this;
 }
 
-Var* LinSolver::addVars(int num, VarType type)
-{
-	Var* var = new Var[num];
+Var *LinSolver::addVars(int num, VarType type) {
+	Var *var = new Var[num];
 	for (int c = 0; c < num; c++) {
-		var[c].col = c + cn, var[c].type = type;
+		var[c].col = c + cn;
+		var[c].type = type;
 	}
 	vars.push_back(std::make_pair(num, var));
 	cn += num;
 	return vars.back().second;
 }
 
-Var& sv::LinSolver::getVar(int idx)
-{
-	assert(idx >= 0 && idx < cn);
+Var &sv::LinSolver::getVar(int idx) {
+	assert(idx >= 0 && idx < static_cast<int>(cn));
 	static Var err_var;
-	for (auto& var : vars) {
-		if (var.first <= idx) {
-			idx -= var.first;
-		}
-		else {
-			return var.second[idx];
+	int offset = 0;
+	for (auto &var : vars) {
+		if (idx < offset + var.first) {
+			return var.second[idx - offset];
 		}
+		offset += var.first;
 	}
 	return err_var;
 }
 
-void LinSolver::addConstr(const Expr& expr, ConstrOper sense, double rhs)
-{
-	if (sense == ConstrOper::LESS_EQUAL) {
+void LinSolver::addConstr(const Expr &expr, ConstrOper sense, double rhs) {
+	if (sense == ConstrOper::EQUAL) {
+		addConstr(expr, ConstrOper::LESS_EQUAL, rhs);
+		addConstr(expr, ConstrOper::GREATER_EQUAL, rhs);
+		return;
+	}
+
 	bn++;
+	if (sense == ConstrOper::LESS_EQUAL) {
 		table.push_back(vector<double>(1, rhs - expr.constant));
 		table.back().insert(table.back().end(), expr.coeffs.begin(), expr.coeffs.end());
 	}
-	else if (sense == ConstrOper::GREATER_EQUAL) {
-		bn++;
+	else {
 		table.push_back(vector<double>(1, expr.constant - rhs));
-		for (int coeff : expr.coeffs) {
+		for (double coeff : expr.coeffs) {
 			table.back().push_back(-coeff);
 		}
 	}
-	else {
-		addConstr(expr, ConstrOper::LESS_EQUAL, rhs);
-		addConstr(expr, ConstrOper::GREATER_EQUAL, rhs);
-	}
 
-	for (int c = table.back().size(); c <= cn; c++) {
+	for (size_t c = table.back().size(); c <= cn; c++) {
 		table.back().push_back(0);
 	}
 }
 
-void LinSolver::setObjective(Expr obje, int _sense)
-{
+void LinSolver::setObjective(Expr obje, int _sense) {
 	assert(_sense == 1 || _sense == -1);
 	if (sense == 0) {
 		table.insert(table.begin(), obje.coeffs);
+	}
+	else {
+		table.front() = obje.coeffs;
+	}
 	table.front().insert(table.front().begin(), -obje.constant);
-		for (int c = obje.coeffs.size() + 1; c <= cn; c++) {
+	for (size_t c = obje.coeffs.size() + 1; c <= cn; c++) {
 		table.front().push_back(0);
 	}
-	}
-	else {
-		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 {
-				table.front().at(col) = 0;
-			}
-		}
-	}
-	for (int row = 0; row < table.front().size(); row++) {
-		table.front().at(row) = _sense * table.front().at(row);
+	for (size_t col = 0; col < table.front().size(); col++) {
+		table.front().at(col) = _sense * table.front().at(col);
 	}
 	sense = _sense;
 }
 
-rtn LinSolver::optimize()
-{
+rtn LinSolver::optimize() {
 	assert(sense);
 	ope_table = table;
 	rtn_ = LOADED;
@@ -163,18 +136,22 @@ rtn LinSolver::optimize()
 	}
 
 	if (rtn_ == OPTIMAL) {
-		cn = ope_table.front().size() - bn;
-		for (int row = 1; row < bn; row++) {
-			if (basic.at(row - 1) - 1 < cn) {
-				getVar(basic.at(row - 1) - 1).val = ope_table.at(row).front();
+		const size_t num_vars = ope_table.front().size() - bn;
+		cn = num_vars;
+		for (size_t i = 0; i < num_vars; i++) {
+			getVar(static_cast<int>(i)).val = 0;
+		}
+		for (size_t row = 1; row < bn; row++) {
+			int var_idx = basic.at(row - 1) - 1;
+			if (var_idx >= 0 && static_cast<size_t>(var_idx) < num_vars) {
+				getVar(var_idx).val = ope_table.at(row).front();
 			}
 		}
 	}
 	return rtn_;
 }
 
-void LinSolver::print()
-{
+void LinSolver::print() {
 	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";
@@ -183,121 +160,163 @@ void LinSolver::print()
 	}
 }
 
-Model::Model()
-{
+Model::Model() {}
 
+namespace {
+
+bool is_integer_type(VarType type) {
+	return type == VarType::INTEGER || type == VarType::BINARY;
 }
 
-rtn Model::optimize()
-{
+bool is_fractional(double val, double eps = 1e-6) {
+	return fabs(val - std::round(val)) > eps;
+}
+
+bool is_better(double candidate, double incumbent, int sense, double eps = 1e-10) {
+	if (sense == MDL_MAXIMIZE) {
+		return candidate > incumbent + eps;
+	}
+	return candidate < incumbent - eps;
+}
+
+} // namespace
+
+rtn Model::optimize() {
 	solver.optimize();
 	if (solver.rtn_ != OPTIMAL) {
 		return solver.rtn_;
 	}
 
-	double global_upper_bound = solver.obj_, global_lower_bound = 0;
+	const int sense = solver.sense;
+	const int num_vars = solver.get(IntAttr::NumVars);
 
-	std::stack<Node> list_;
-	
-	Node root_node = { solver, 0, solver.obj_ };
-	Node incumbent_node = root_node;
-
-	list_.push(root_node);
-	while (list_.size() && global_upper_bound - global_lower_bound > 1e-10) {
-		Node current_node = list_.top();
-		list_.pop();
-		current_node.solver.optimize();
-		
-		if (current_node.solver.get(IntAttr::Status) == OPTIMAL) {
-			int branch_var_index = -1;
-
-			for (int i = 0; i < current_node.solver.get(IntAttr::NumVars); i++) {
-				if (current_node.solver.getVar(i).type == VarType::INTEGER) {
-					if (fabs(int(current_node.solver.getVar(i).val) - current_node.solver.getVar(i).val) > 1e-10) {
-						branch_var_index = i;
+	bool has_integer = false;
+	for (int i = 0; i < num_vars; i++) {
+		if (is_integer_type(solver.getVar(i).type)) {
+			has_integer = true;
 			break;
 		}
 	}
+	if (!has_integer) {
+		return solver.rtn_;
+	}
+
+	double best_obj = (sense == MDL_MAXIMIZE) ? -std::numeric_limits<double>::infinity()
+	                                          : std::numeric_limits<double>::infinity();
+	bool found_integer = false;
+	Node incumbent_node{solver, solver.get(DoubleAttr::Obj)};
+
+	std::stack<Node> open_nodes;
+	open_nodes.push(Node{solver, solver.get(DoubleAttr::Obj)});
+	while (!open_nodes.empty()) {
+		Node current_node = std::move(open_nodes.top());
+		open_nodes.pop();
+
+		if (found_integer && !is_better(current_node.bound, best_obj, sense)) {
+			continue;
+		}
+
+		current_node.solver.optimize();
+		if (current_node.solver.get(IntAttr::Status) != OPTIMAL) {
+			continue;
+		}
+
+		const double lp_obj = current_node.solver.get(DoubleAttr::Obj);
+		current_node.bound = lp_obj;
+
+		if (found_integer && !is_better(lp_obj, best_obj, sense)) {
+			continue;
+		}
+
+		int branch_var_index = -1;
+		for (int i = 0; i < current_node.solver.get(IntAttr::NumVars); i++) {
+			Var &var = current_node.solver.getVar(i);
+			if (!is_integer_type(var.type)) {
+				continue;
+			}
+			if (is_fractional(var.val)) {
+				if (branch_var_index < 0) {
+					branch_var_index = i;
+				}
+			}
+			else {
+				// Snap near-integer values so later checks / output stay clean
+				var.val = std::round(var.val);
+			}
 		}
 
 		if (branch_var_index == -1) {
-				current_node.lower_bound = current_node.solver.obj_;
-				current_node.upper_bound = current_node.solver.obj_;
-				if (current_node.lower_bound > global_lower_bound) {
-					global_lower_bound = current_node.lower_bound;
+			if (!found_integer || is_better(lp_obj, best_obj, sense)) {
+				best_obj = lp_obj;
+				found_integer = true;
 				incumbent_node = current_node;
 			}
+			continue;
 		}
-			else {
-				if (current_node.upper_bound >= global_lower_bound) {
-					const Var& branch_var = current_node.solver.getVar(branch_var_index);
-					int left_var_bound = branch_var.val;
-					int right_var_bound = branch_var.val + 1;
+
+		const Var &branch_var = current_node.solver.getVar(branch_var_index);
+		const int left_bound = static_cast<int>(std::floor(branch_var.val));
+		const int right_bound = left_bound + 1;
 
 		Node left_node = current_node;
-					left_node.solver.addConstr(branch_var, ConstrOper::LESS_EQUAL, left_var_bound);
-					list_.push(left_node);
+		left_node.solver.addConstr(branch_var, ConstrOper::LESS_EQUAL, left_bound);
+		open_nodes.push(std::move(left_node));
 
 		Node right_node = current_node;
-					right_node.solver.addConstr(branch_var, ConstrOper::GREATER_EQUAL, right_var_bound);
-					list_.push(right_node);
-				}
-			}
+		right_node.solver.addConstr(branch_var, ConstrOper::GREATER_EQUAL, right_bound);
+		open_nodes.push(std::move(right_node));
 	}
+
+	if (!found_integer) {
+		return solver.rtn_ = INFEASIBLE;
 	}
+
 	solver.rtn_ = incumbent_node.solver.rtn_;
 	solver.obj_ = incumbent_node.solver.obj_;
-	for (int i = 0; i < solver.cn; i++) {
+	for (int i = 0; i < solver.get(IntAttr::NumVars); i++) {
 		solver.getVar(i).val = incumbent_node.solver.getVar(i).val;
 	}
 	return solver.rtn_;
 }
 
-Var* sv::Model::addVars(int col, VarType type)
-{
+Var *sv::Model::addVars(int col, VarType type) {
 	return solver.addVars(col, type);
 }
 
-void sv::Model::addConstr(const Expr& expr, ConstrOper sense, double rhs)
-{
-	return solver.addConstr(expr, sense, rhs);
+void sv::Model::addConstr(const Expr &expr, ConstrOper sense, double rhs) {
+	solver.addConstr(expr, sense, rhs);
 }
 
-void sv::Model::setObjective(Expr obje, int sense)
-{
-	return solver.setObjective(obje, sense);
+void sv::Model::setObjective(Expr obje, int sense) {
+	solver.setObjective(obje, sense);
 }
 
-double sv::Model::get(DoubleAttr attr)
-{
+double sv::Model::get(DoubleAttr attr) {
 	return solver.get(attr);
 }
 
-int sv::Model::get(IntAttr attr)
-{
+int sv::Model::get(IntAttr attr) {
 	return solver.get(attr);
 }
 
-double LinSolver::get(DoubleAttr attr)
-{
+double LinSolver::get(DoubleAttr attr) {
 	return -sense * obj_;
 }
 
-int LinSolver::get(IntAttr attr)
-{
+int LinSolver::get(IntAttr attr) {
 	switch (attr) {
 	case IntAttr::NumVars:
-		return cn;
+		return static_cast<int>(cn);
 	case IntAttr::Status:
 		return rtn_;
-	}
+	default:
 		return -1;
+	}
 }
 
-double LinSolver::_simplex()
-{
+double LinSolver::_simplex() {
 	pair<size_t, size_t> t;
-	while (1) {
+	while (true) {
 		rtn_ = _pivot(t);
 		if (rtn_ == OPTIMAL || rtn_ == UNBOUNDED) {
 			break;
@@ -307,39 +326,40 @@ double LinSolver::_simplex()
 	return obj_ = ope_table.front().front();
 }
 
-rtn LinSolver::feasible_solution()
-{
-	for (int row = 1; row < bn; row++) {
+rtn LinSolver::feasible_solution() {
+	for (size_t row = 1; row < bn; row++) {
 		ope_table.front().push_back(0);
-		for (int col = 1; col < bn; col++) {
-			ope_table.at(row).push_back(col == row ? 1 : 0);
+		for (size_t col = 1; col < bn; col++) {
+			ope_table.at(row).push_back(col == row ? 1.0 : 0.0);
 		}
 	}
 	cn = ope_table.front().size();
 	basic.clear();
+	basic.reserve(bn - 1);
 	for (size_t i = 1; i < bn; i++) {
-		basic.push_back(cn - bn + i);
+		basic.push_back(static_cast<int>(cn - bn + i));
 	}
 
-	// === 判断初始解是否为可行解 ===
+	// Check whether the initial basic solution is feasible
 	bool initial_feasible = true;
-	for (int row = 1; row < bn; row++) {
+	for (size_t row = 1; row < bn; row++) {
 		if (ope_table.at(row).front() < 0) {
 			initial_feasible = false;
 			break;
 		}
 	}
 
-	// === 构造初始可行解 === 
+	// Two-phase method when the initial basis is infeasible
 	if (!initial_feasible) {
 		vector<double> coeff = ope_table.front();
-		ope_table.front() = vector<double>(cn, .0);
+		ope_table.front() = vector<double>(cn, 0.0);
 		ope_table.front().push_back(1);
-		pair<size_t, size_t> t = { -1 ,cn };
+		pair<size_t, size_t> t = {static_cast<size_t>(-1), cn};
 
-		for (int row = 1; row < bn; row++) {
+		for (size_t row = 1; row < bn; row++) {
 			ope_table.at(row).push_back(-1);
-			if (t.first == -1 || ope_table.at(row).front() < ope_table.at(t.first).front()) {
+			if (t.first == static_cast<size_t>(-1) ||
+			    ope_table.at(row).front() < ope_table.at(t.first).front()) {
 				t.first = row;
 			}
 		}
@@ -349,29 +369,30 @@ rtn LinSolver::feasible_solution()
 			return rtn_ = INFEASIBLE;
 		}
 		rtn_ = LOADED;
-		// if the x0 in B, we should pivot it.
-		auto iter = find(basic.begin(), basic.end(), cn);
+
+		// If artificial variable remains basic, pivot it out
+		auto iter = find(basic.begin(), basic.end(), static_cast<int>(cn));
 		if (iter != basic.end()) {
-			for (int col = 1; col < ope_table.front().size(); col++) {
+			for (size_t 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);
+					t = make_pair(static_cast<size_t>(iter - basic.begin() + 1), col);
 					_gaussian(t);
 					break;
 				}
 			}
 		}
 
-		for (int row = 0; row < bn; row++) {
+		for (size_t row = 0; row < bn; row++) {
 			ope_table.at(row).pop_back();
 		}
 
-		// recover the coefficient line
-		for (int col = 0; col < cn; col++) {
+		// Restore original objective row and re-express in current basis
+		for (size_t col = 0; col < cn; col++) {
 			ope_table.front().at(col) = coeff.at(col);
 		}
-		for (int row = 1; row <= basic.size(); row++) {
-			int norm = ope_table.front().at(basic.at(row - 1));
-			for (int col = 0; col < cn; col++) {
+		for (size_t row = 1; row <= basic.size(); row++) {
+			double norm = ope_table.front().at(basic.at(row - 1));
+			for (size_t col = 0; col < cn; col++) {
 				ope_table.front().at(col) -= norm * ope_table.at(row).at(col);
 			}
 		}
@@ -379,15 +400,15 @@ rtn LinSolver::feasible_solution()
 	return rtn_;
 }
 
-rtn LinSolver::_pivot(pair<size_t, size_t>& p)
-{
+rtn LinSolver::_pivot(pair<size_t, size_t> &p) {
 	p = make_pair(0, 0);
 	double cmin = DBL_MAX;
-	vector<double> coef = ope_table.front();
+	const vector<double> &coef = ope_table.front();
 
-	// === 非主轴元素中找最小值 ===
+	// Entering variable: most negative reduced cost
 	for (size_t col = 1; col < coef.size(); col++) {
-		if (cmin > coef.at(col) && find(basic.begin(), basic.end(), col) == basic.end()) {
+		if (cmin > coef.at(col) &&
+		    find(basic.begin(), basic.end(), static_cast<int>(col)) == basic.end()) {
 			cmin = coef.at(col);
 			p.second = col;
 		}
@@ -395,51 +416,47 @@ rtn LinSolver::_pivot(pair<size_t, size_t>& p)
 	if (cmin >= 0) {
 		return OPTIMAL;
 	}
+
 	double bmin = DBL_MAX;
 	for (size_t row = 1; row < bn; row++) {
-		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) {
+		const double pivot_col = ope_table.at(row).at(p.second);
+		if (pivot_col > 0) {
+			const double tmp = ope_table.at(row).front() / pivot_col;
+			if (bmin > tmp) {
 				bmin = tmp;
 				p.first = row;
 			}
 		}
+	}
 
-	if (abs(bmin - DBL_MAX) < 1e-10) {
+	if (bmin >= DBL_MAX / 2) {
 		return UNBOUNDED;
 	}
 
-	for (auto iter = basic.begin(); iter != basic.end(); iter++) {
-		if (ope_table.at(p.first).at(*iter) != 0) {
-			*iter = p.second;
-			break;
-		}
-	}
-	assert(basic.at(p.first - 1) == p.second);
+	basic.at(p.first - 1) = static_cast<int>(p.second);
 	return PIVOT;
 }
 
-void LinSolver::_gaussian(pair<size_t, size_t> p)
-{
+void LinSolver::_gaussian(pair<size_t, size_t> p) {
 	size_t x = p.first, y = p.second;
 
-	// === 主行归一化 ===
+	// Normalize pivot row
 	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;
 	}
 
-	// === 其余行变换 ===
+	// Eliminate pivot column in other rows
 	for (size_t row = 0; row < bn; row++) {
 		if (row == x) {
 			continue;
 		}
 		if (ope_table.at(row).at(y) != 0) {
-			double norm = ope_table.at(row).at(y);
+			double row_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);
+				ope_table.at(row).at(col) -= row_norm * ope_table.at(x).at(col);
 			}
 		}
 	}
-	basic.at(x - 1) = y;	// 换元
+	basic.at(x - 1) = static_cast<int>(y);
 }
-