summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--.gitignore4
-rw-r--r--Makefile24
-rw-r--r--lp.cpp19
-rw-r--r--problem.cpp311
-rw-r--r--problem.h126
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 $@ $<
diff --git a/lp.cpp b/lp.cpp
new file mode 100644
index 0000000..2a4a283
--- /dev/null
+++ b/lp.cpp
@@ -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);