From 3b36e04bc4e0ef7ac3c968eb854f266946516b52 Mon Sep 17 00:00:00 2001
From: tomsmeding <tom.smeding@gmail.com>
Date: Sat, 5 Nov 2016 18:28:07 +0100
Subject: Initial

---
 .gitignore  |   4 +
 Makefile    |  24 +++++
 lp.cpp      |  19 ++++
 problem.cpp | 311 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 problem.h   | 126 ++++++++++++++++++++++++
 5 files changed, 484 insertions(+)
 create mode 100644 .gitignore
 create mode 100644 Makefile
 create mode 100644 lp.cpp
 create mode 100644 problem.cpp
 create mode 100644 problem.h

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);
-- 
cgit v1.2.3-70-g09d2