diff options
Diffstat (limited to 'problem.cpp')
| -rw-r--r-- | problem.cpp | 311 |
1 files changed, 311 insertions, 0 deletions
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; +} |