LU/general-solver
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

调整目录结构

Browse Source
This commit is contained in:
hit_lu
2025-05-19 14:02:28 +08:00
parent 0ae1785e6c
commit 26ac3ddebd
7 changed files with 494 additions and 616 deletions
Download Patch File Download Diff File Expand all files Collapse all files

11
CMakeLists.txt
View File

@ -12,22 +12,17 @@ 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}

163
inc/common.hpp Normal file
View File

@ -0,0 +1,163 @@
#pragma once
#include <vector>
#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<std::vector<double>>;
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<double> coeffs;
std::vector<Var> 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);
};
}

5
inc/node.hpp Normal file
View File

@ -0,0 +1,5 @@
#pragma once
class Node {
};

317
inc/solver.hpp
View File

@ -1,304 +1,7 @@
#pragma once
#include <vector>
using std::vector;
#include <common.hpp>
class Expr;
class Var;
#define MDL_MINIMIZE 1
#define MDL_MAXIMIZE -1
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);
typedef vector<vector<double>> 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<double> coeffs;
std::vector<Var> 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);
};
namespace sv {
class Model
{
@ -308,20 +11,24 @@ public:
void addConstr(const Expr& expr, ConstrOper sense, double rhs);
void setObjective(Expr obje, int sense = MDL_MAXIMIZE);
void print();
Rtn optimize();
rtn optimize();
double get(DoubleAttr attr);
int get(IntAttr attr);
private:
double _simplex();
Rtn _pivot(std::pair<size_t, size_t>& p);
rtn _pivot(std::pair<size_t, size_t>& p);
void _gaussian(std::pair<size_t, size_t> p);
Var* vars;
Matrix matrix;
matrix mt;
size_t cn, bn;
vector<int> basic;
Rtn rtn;
double res;
std::vector<int> basic;
rtn rtn_;
double res_;
int sense;
};
}

225
src/common.cpp Normal file
View File

@ -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;
}

7
src/main.cpp
View File

@ -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:

350
src/solver.cpp
View File

@ -1,148 +1,17 @@
#include "solver.hpp"
#include "common.hpp"
#include <iostream>
#include <cassert>
#include <algorithm>
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<double>(1, rhs - expr.constant));
matrix.back().insert(matrix.back().end(), expr.coeffs.begin(), expr.coeffs.end());
mt.push_back(vector<double>(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<double>(1, expr.constant - rhs));
mt.push_back(vector<double>(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()
// === <20>жϳ<D0B6>ʼ<EFBFBD><CABC><EFBFBD>Ƿ<EFBFBD>Ϊ<EFBFBD><CEAA><EFBFBD>н<EFBFBD> ===
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()
// === <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʼ<EFBFBD><CABC><EFBFBD>н<EFBFBD> ===
if (!initial_feasible) {
vector<double> coeff = matrix.front();
matrix.front() = vector<double>(cn, .0);
matrix.front().push_back(1);
vector<double> coeff = mt.front();
mt.front() = vector<double>(cn, .0);
mt.front().push_back(1);
pair<size_t, size_t> 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<size_t, size_t> 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<size_t, size_t>& p)
rtn Model::_pivot(pair<size_t, size_t>& p)
{
p = make_pair(0, 0);
double cmin = INT_MAX;
vector<double> coef = matrix.front();
vector<double> coef = mt.front();
// === <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ԫ<EFBFBD><D4AA><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Сֵ ===
for (size_t i = 1; i < coef.size(); i++) {
@ -343,29 +212,29 @@ Rtn Model::_pivot(pair<size_t, size_t>& 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<size_t, size_t> p)
@ -373,9 +242,9 @@ void Model::_gaussian(pair<size_t, size_t> p)
size_t x = p.first, y = p.second;
// === <20><><EFBFBD>й<EFBFBD>һ<EFBFBD><D2BB> ===
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;
}
// === <20><><EFBFBD><EFBFBD><EFBFBD>б任 ===
@ -383,100 +252,13 @@ void Model::_gaussian(pair<size_t, size_t> 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; // <20><>Ԫ
}
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;
}