path: root/problem.cpp

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