diff options
| -rw-r--r-- | .gitignore | 4 | ||||
| -rw-r--r-- | Makefile | 24 | ||||
| -rw-r--r-- | lp.cpp | 19 | ||||
| -rw-r--r-- | problem.cpp | 311 | ||||
| -rw-r--r-- | problem.h | 126 |
5 files changed, 484 insertions, 0 deletions
diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..70014ba --- /dev/null +++ b/.gitignore @@ -0,0 +1,4 @@ +*.o +lp +*.dSYM +.DS_Store diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..a8bf246 --- /dev/null +++ b/Makefile @@ -0,0 +1,24 @@ +CXX = g++ +CXXFLAGS = -Wall -Wextra -std=c++11 -fwrapv +ifneq ($(DEBUG),) + CXXFLAGS += -g +else + CXXFLAGS += -O2 +endif +BIN = lp + +.PHONY: all clean remake + +all: $(BIN) + +clean: + rm -rf $(BIN) *.o *.dSYM + +remake: clean all + + +$(BIN): $(patsubst %.cpp,%.o,$(wildcard *.cpp)) + $(CXX) -o $@ $^ + +%.o: %.cpp $(wildcard *.h) + $(CXX) $(CXXFLAGS) -c -o $@ $< @@ -0,0 +1,19 @@ +#include <iostream> +#include "problem.h" + +using namespace std; + +int main(){ + Problem problem(cin); + cout<<problem<<endl<<endl; + + problem.solve(); + cout<<problem<<endl; + + cout<<"solution value: "<<problem.solutionValue()<<endl; + + vector<double> vec=problem.solutionVars(); + cout<<"variable values:"; + for(double v : vec)cout<<' '<<v; + cout<<endl; +} diff --git a/problem.cpp b/problem.cpp new file mode 100644 index 0000000..3eb7ecc --- /dev/null +++ b/problem.cpp @@ -0,0 +1,311 @@ +#include <iomanip> +#include <stdexcept> +#include <cassert> +#include "problem.h" + +#define DEBUG + +using namespace std; + +Problem::Problem(istream &in){ + in>>*this; +} + +template <typename T,typename U> +void rowopadd(int n,T &dst,const U &src,double times){ + for(int i=0;i<n;i++)dst[i]+=src[i]*times; +} + +template <typename T> +void rowopmult(int n,T &dst,double times){ + for(int i=0;i<n;i++)dst[i]*=times; +} + +template <typename T> +int min_index(const vector<T> &v){ + if(v.size()==0)throw logic_error("No minimum element in empty vector"); + int mi=0; + int size=v.size(); + for(int i=1;i<size;i++){ + if(v[i]<v[mi]){ + mi=i; + } + } + return mi; +} + +void Problem::removefrombasis(int varidx){ + if(zfunc[varidx]==0){ +#ifdef DEBUG + cerr<<"Not needing to remove variable "<<varidx<<" from basis"<<endl; +#endif + return; + } + + const int nrestr=restr.height(),nvars=restr.width(); + + int restridx=-1; + for(int i=0;i<nrestr;i++){ + if(basis[i]==varidx){ + restridx=i; + break; + } + } + throw logic_error("Variable to remove from basis is not in basis"); + + if(restr[restridx][varidx]==0)throw logic_error("Basis variable column not in identity form"); + + zvalue-=zfunc[varidx]*bvec[restridx]; + rowopadd(nvars,zfunc,restr[restridx],-zfunc[varidx]); + +#ifdef DEBUG + cerr<<"Removed variable "<<varidx<<" from basis:"<<endl; + cerr<<*this<<endl<<endl; +#endif +} + +void Problem::solve(){ + const int nvars=restr.width(); + + bool haveart=false; + for(VarType type : vartype){ + if(type==VT_ART){ + haveart=true; + break; + } + } + if(haveart){ +#ifdef DEBUG + cerr<<" === SOLVING ARTIFICIAL PROBLEM ==="<<endl; +#endif + vector<double> origzfunc=zfunc; + for(int j=0;j<nvars;j++){ + zfunc[j]=vartype[j]==VT_ART; + } +#ifdef DEBUG + cerr<<*this<<endl<<endl; +#endif + solve_noart(); +#ifdef DEBUG + cerr<<*this<<endl<<endl; +#endif + if(zvalue!=0)throw invalid_argument("No feasible solution (no zero value of artificial variables)"); + zfunc=origzfunc; + + int nart=0; + for(int j=nvars-1;j>=0&&vartype[j]==VT_ART;j--){ + removefrombasis(j); + nart++; + } + vartype.resize(nvars-nart); + zfunc.resize(nvars-nart); + restr.resize(nvars-nart,restr.height()); + +#ifdef DEBUG + cerr<<" === SOLVING ORIGINAL PROBLEM ==="<<endl; +#endif + solve_noart(); + } +} + +void Problem::solve_noart(){ + const int nrestr=restr.height(),nvars=restr.width(); + + basis.clear(); + basis.resize(nrestr,-1); + for(int j=0;j<nvars;j++){ //find basis variables + int i; + for(i=0;i<nrestr;i++)if(restr[i][j]!=0)break; + if(i==nrestr)throw invalid_argument("Variable with zero matrix column"); + if(restr[i][j]<=0)continue; + int onerestr=i; + if(basis[onerestr]!=-1)continue; + i++; + for(;i<nrestr;i++)if(restr[i][j]!=0)break; + if(i!=nrestr)continue; + basis[onerestr]=j; + double value=restr[onerestr][j]; + assert(value!=0); + for(int j2=0;j2<nvars;j2++)restr[onerestr][j2]/=value; + bvec[onerestr]/=value; + } + + for(int j : basis){ + if(j==-1)throw invalid_argument("No feasible solution (overspecified system)"); + } + +#ifdef DEBUG + cerr<<"basis ="; + for(int i=0;i<nrestr;i++)cerr<<' '<<basis[i]; + cerr<<endl; + cerr<<*this<<endl<<endl; +#endif + + for(int i=0;i<nrestr;i++){ //express z in terms of non-basis variables + if(zfunc[basis[i]]!=0){ + zvalue-=zfunc[basis[i]]*bvec[i]; + rowopadd(nvars,zfunc,restr[i],-zfunc[basis[i]]); + } + } + +#ifdef DEBUG + cerr<<"Expressed z in terms of non-basis variables:"<<endl; + cerr<<*this<<endl<<endl; +#endif + + while(true){ + int pivotvar=min_index(zfunc); //find pivot column + if(zfunc[pivotvar]>=0)break; //optimal solution found + int pivotrestr=-1; + for(int i=0;i<nrestr;i++){ //find pivot row + if(restr[i][pivotvar]<=0)continue; + if(pivotrestr==-1||bvec[i]/restr[i][pivotvar]<bvec[pivotrestr]/restr[pivotrestr][pivotvar]){ + pivotrestr=i; + } + } + if(pivotrestr==-1){ + throw invalid_argument("Unbounded problem (no positive entry in pivot column)"); + } + if(restr[pivotrestr][pivotvar]!=1){ //normalise row + bvec[pivotrestr]/=restr[pivotrestr][pivotvar]; + auto row=restr[pivotrestr]; + rowopmult(nvars,row,1/restr[pivotrestr][pivotvar]); + } + for(int i=0;i<nrestr;i++){ //zero this column in other rows + if(i==pivotrestr)continue; + if(restr[i][pivotvar]==0)continue; + bvec[i]-=restr[i][pivotvar]*bvec[pivotrestr]; + auto row=restr[i]; + rowopadd(nvars,row,restr[pivotrestr],-restr[i][pivotvar]); + } + if(zfunc[pivotvar]!=0){ //and the z row + zvalue-=zfunc[pivotvar]*bvec[pivotrestr]; + rowopadd(nvars,zfunc,restr[pivotrestr],-zfunc[pivotvar]); + } + basis[pivotrestr]=pivotvar; //replace variable in basis + +#ifdef DEBUG + cerr<<"basis ="; + for(int i=0;i<nrestr;i++)cerr<<' '<<basis[i]; + cerr<<endl; +#endif + } +} + +double Problem::solutionValue() const { + return zvalue==0?0:-zvalue; //prevent -0 +} + +vector<double> Problem::solutionVars() const { + vector<double> normal; + int nrestr=restr.height(),nvars=restr.width(); + int nnormal; + for(nnormal=0;nnormal<nvars;nnormal++)if(vartype[nnormal]!=VT_NORMAL)break; + normal.resize(nnormal); + for(int i=0;i<nrestr;i++){ + if(vartype[basis[i]]==VT_NORMAL){ + normal[basis[i]]=bvec[i]; + if(normal[basis[i]]==0)normal[basis[i]]=0; //prevent -0 + } + } + return normal; +} + +template <typename T> +void readthrow(istream &in,T &dst){ + in>>dst; + if(in.fail())throw invalid_argument("Failure reading input"); +} + +istream& operator>>(istream &in,Problem &prob){ + string type; + in>>type; + bool negate=false; + if(type=="max")negate=true; + else if(type!="min")throw invalid_argument("Invalid LP problem type"); + int nvars,nrestr; + in>>nvars>>nrestr; + if(nvars<=0)throw invalid_argument("Invalid number of variables"); + if(nrestr<=0)throw invalid_argument("Invalid number of restrictions"); + prob.vartype.clear(); prob.vartype.resize(nvars); + prob.zfunc.clear(); prob.zfunc.resize(nvars); + prob.restr.clear(); prob.restr.resize(nvars,nrestr); + prob.bvec.clear(); prob.bvec.resize(nrestr); + prob.basis.clear(); prob.basis.resize(nrestr,-1); + for(int i=0;i<nvars;i++){ + readthrow(in,prob.zfunc[i]); + if(negate)prob.zfunc[i]=-prob.zfunc[i]; + } + vector<double> addslack(nrestr,0),addart(nrestr,0); + for(int i=0;i<nrestr;i++){ + for(int j=0;j<nvars;j++){ + readthrow(in,prob.restr[i][j]); + } + readthrow(in,type); + readthrow(in,prob.bvec[i]); + if(prob.bvec[i]<0){ + prob.bvec[i]=-prob.bvec[i]; + for(int j=0;j<nvars;j++)prob.restr[i][j]=-prob.restr[i][j]; + if(type=="<=")type=">="; + else if(type==">=")type="<="; + } + if(type=="<="){ + addslack[i]=1; + } else if(type==">="){ + addslack[i]=-1; + addart[i]=1; + } else if(type=="="){ + addart[i]=1; + } else { + throw invalid_argument("Invalid restriction type"); + } + } + for(int i=0;i<nrestr;i++){ + if(addslack[i]==0)continue; + prob.vartype.push_back(Problem::VT_SLACK); + prob.zfunc.emplace_back(); + prob.restr.addcolumn(); + prob.restr[i].back()=addslack[i]; + } + for(int i=0;i<nrestr;i++){ + if(addart[i]==0)continue; + prob.vartype.push_back(Problem::VT_ART); + prob.zfunc.emplace_back(); + prob.restr.addcolumn(); + prob.restr[i].back()=addart[i]; + } + return in; +} + +ostream& operator<<(ostream &os,const Problem &prob){ + const int cellwid=4; + os<<"min."; + for(Problem::VarType type : prob.vartype){ + switch(type){ + case Problem::VT_NORMAL: os<<' '<<setw(cellwid)<<'N'; break; + case Problem::VT_SLACK: os<<' '<<setw(cellwid)<<'S'; break; + case Problem::VT_ART: os<<' '<<setw(cellwid)<<'A'; break; + default: assert(false); + } + } + os<<endl; + + os<<"z = "; + for(double c : prob.zfunc){ + os<<' '<<setw(cellwid)<<c; + } + os<<" = "<<setw(cellwid)<<prob.zvalue<<endl; + + for(int i=0;i<prob.restr.height();i++){ + if(i==0)os<<"s.t."; + else os<<" "; + for(int j=0;j<prob.restr.width();j++){ + os<<' '<<setw(cellwid)<<prob.restr[i][j]; + } + os<<" = "<<setw(cellwid)<<prob.bvec[i]; + os<<" (var "<<prob.basis[i]<<')'; + if(i!=prob.restr.height()-1)os<<endl; + } + + return os; +} diff --git a/problem.h b/problem.h new file mode 100644 index 0000000..678ab59 --- /dev/null +++ b/problem.h @@ -0,0 +1,126 @@ +#pragma once + +#include <iostream> +#include <vector> +#include <utility> + +using namespace std; + +template <typename T> +class Matrix{ + vector<vector<T>> mat; //list of rows + int W=0,H=0; + +public: + template <typename U> + class RowAccess{ + vector<T> &row; + public: + RowAccess(vector<T> &row):row(row){} + T& operator[](int x){return row[x];} + const T& operator[](int x) const {return row[x];} + T& front(){return row.front();} + const T& front() const {return row.front();} + T& back(){return row.back();} + const T& back() const {return row.back();} + }; + + template <typename U> + class ConstRowAccess{ + const vector<T> &row; + public: + ConstRowAccess(const vector<T> &row):row(row){} + const T& operator[](int x) const {return row[x];} + const T& front() const {return row.front();} + const T& back() const {return row.back();} + }; + + Matrix() = default; + + Matrix(int W,int H) + :W(W),H(H){ + mat.resize(H); + for(auto &row : mat){ + row.resize(W); + } + } + + void addrow(){ + mat.emplace_back(W,T()); + H++; + } + + void addcolumn(){ + for(auto &row : mat){ + row.emplace_back(); + } + W++; + } + + void resize(int w,int h){ + if(h<H){ + H=h; + mat.resize(H); + } + if(w!=W){ + W=w; + for(auto &row : mat){ + row.resize(w); + } + } + if(h>H){ + mat.resize(h); + for(int i=H;i<h;i++){ + mat[i].resize(W); + } + H=h; + } + } + + void clear(){ + mat.clear(); + } + + RowAccess<T> operator[](int y){ + return RowAccess<T>(mat[y]); + } + + ConstRowAccess<T> operator[](int y) const { + return ConstRowAccess<T>(mat[y]); + } + + int width() const {return W;} + + int height() const {return H;} +}; + +class Problem{ + enum VarType{ + VT_NORMAL=0, + VT_SLACK, + VT_ART, + }; + vector<VarType> vartype; + vector<double> zfunc; + double zvalue=0; + Matrix<double> restr; + vector<double> bvec; + vector<int> basis; + + void removefrombasis(int varidx); + void solve_noart(); + +public: + Problem(const Problem&) = default; + Problem(istream &in); + + void solve(); + double solutionValue() const; + vector<double> solutionVars() const; + + friend istream& operator>>(istream &in,Problem &prob); + friend ostream& operator<<(ostream &os,const Problem &prob); +}; + +istream& operator>>(istream &in,Problem &prob); +ostream& operator<<(ostream &os,const Problem &prob); |