Linear Programming

linearprogram.hpp

#include <vector>
#include <iostream>
#include <algorithm>
#include <limits>
#include <utility>

class StandardForm
{
	friend std::ostream &operator<<(std::ostream &out, const StandardForm &standrdform);
	friend class SlackForm;
	template <class T> using vector = std::vector<T>;
	vector<vector<double>> a;
	vector<double> b, c;
public:
	StandardForm(const vector<vector<double>> a, 
		const vector<double> &b,
		const vector<double> &c): a(a), b(b), c(c) {}
};

std::ostream &operator<<(std::ostream &out, const StandardForm &standrdform)
{
	using std::endl;
	out << "object:" << endl;
	for (int i = 0; i < standrdform.c.size(); i++) {
		if (i) out << ' ';
		out << standrdform.c[i] << "*x[" << i << ']';
	}
	out << std::endl << "constraints:" << endl;
	for (int i = 0; i < standrdform.a.size(); i++) {
		for (int j = 0; j < standrdform.a[i].size(); j++) {
			if (j) out << ' ';
			out << standrdform.a[i][j] << "*x[" << i << ']';
		}
		out << "<=" << standrdform.b[i] << endl;
	}
}

class SlackForm
{
	friend std::ostream &operator<<(std::ostream &out, const SlackForm &slackform);

	template <class T> using vector = std::vector<T>;
	template <class T1, class T2> using pair = std::pair<T1, T2>;

	vector<int> nonbasics, basics;
	vector<vector<double>> a;
	vector<double> b, c;
	double v;

	void pivot(int leaving, int entering)
	{
		// compute the coefficients of the equation for new basic variable x[entreing]
		b[entering] = b[leaving] / a[leaving][entering];
		for (int j : nonbasics)
			if (j != entering)
				a[entering][j] = a[leaving][j] / a[leaving][entering];
		a[entering][leaving] = 1.0 / a[leaving][entering];
		// compute the cofficients of the remaining constraints
		for (int i : basics)
			if (i != leaving) {
				b[i] = b[i] - a[i][entering] * b[entering];
				for (int j : nonbasics)
					if (j != entering)
						a[i][j] = a[i][j] - a[i][entering] * a[entering][j];
				a[i][leaving] = -a[i][entering] * a[entering][leaving];
			}
		// compute the object function
		v += c[entering] * b[entering];
		for (int i : nonbasics)
			if (i != entering)
				c[i] -= c[entering]*a[entering][i];
		c[leaving] = -c[entering]*a[entering][leaving];
		// compute new set of basic and nonbasic variables
		for (int& i : nonbasics)
			if (i == entering) {
				i = leaving;
				break;
			}
		for (int& i : basics)
			if (i == leaving) {
				i = entering;
				break;
			}
		std::sort(nonbasics.begin(), nonbasics.end());
		std::sort(basics.begin(), basics.end());
	}

	ReturnCode initializeSimplex()
	{
		// find 'leaving' that minimize b[leaving]
		int leaving = basics[0];
		for (int i = 1; i < basics.size(); i++)
			if (b[basics[i]] < b[leaving])
				leaving = basics[i];
		// is initialize feasible ?
		if (b[leaving] >= 0) return OK;
		// construct a helper linear programming 
		const int NB = nonbasics.size() + basics.size();
		nonbasics.push_back(NB);
		vector<double> tempC(NB+1);
		tempC[NB] = -1;
		std::swap(c, tempC);
		b.push_back(0);
		for (vector<double> &rows : a)
			rows.push_back(-1);
		a.push_back(vector<double>(NB+1));
		pivot(leaving, NB);
		// solve helper linear programming
		vector<double> deltas(basics. size());
		while (true) {
			// find 'entering' that c[entering] > 0
			int entering = -1;
			for (int i : nonbasics)
				if (c[i] > 0) {
					entering = i;
					break;
				}
			if (entering == -1) break;
			// find 'leaving' that minimize b[leaving]/a[leaving][entering]
			for (int i = 0; i < basics.size(); i++)
				if (a[basics[i]][entering] > 0)
					deltas[i] = b[basics[i]]/a[basics[i]][entering];
				else
					deltas[i] = std::numeric_limits<double>::infinity();
			int leavingIndex = 0;
			for (int i = 1; i < deltas.size(); i++)
				if (deltas[i] < deltas[leavingIndex])
					leavingIndex = i;
			if (deltas[leavingIndex] == std::numeric_limits<double>::infinity())
				return UNBOUNDED;
			else
				pivot(basics[leavingIndex], entering);
		}
		double NBValue;
		bool isBasic = false;
		for (int i : nonbasics)
			if (i == NB) {
				NBValue = 0;
				break;
			}
		for (int i : basics)
			if (i == NB) {
				NBValue = b[i];
				isBasic = true;
				break;
			}
		if (NBValue == 0) {
			if (isBasic)
				pivot(NB, nonbasics[0]);
			// remove x[NB]
			nonbasics.pop_back();
			std::swap(tempC, c);
			v = 0;
			for (int i : basics) {
				for (int j : nonbasics)
					c[j] -= a[i][j] * c[i];
				v += b[i] * c[i];
			}
			b.pop_back();
			a.pop_back();
			for (vector<double> &row : a)
				row.pop_back();
			return OK;
		} else return INFEASIBLE;
	}

public:

	enum ReturnCode { OK, UNBOUNDED, INFEASIBLE };

	SlackForm(const StandardForm &standrdform) {
		// convert standard form into slack form
		const int N = standrdform.c.size();
		const int B = standrdform.b.size();
		nonbasics = vector<int>(N);
		for (int i = 0; i < N; i++)
			nonbasics[i] = i;
		basics = vector<int>(B);
		for (int i = 0; i < B; i++)
			basics[i] = N+i;
		a = vector<vector<double>>(N+B, vector<double>(N+B));
		for (int i : basics)
			for (int j : nonbasics)
				a[i][j] = standrdform.a[i-N][j];
		b = vector<double>(N+B);
		for (int i : basics)
			b[i] = standrdform.b[i-N];
		c = vector<double>(N+B);
		for (int i : nonbasics)
			c[i] = standrdform.c[i];
		v = 0;
	}

	pair<ReturnCode, vector<double>> simplex()
	{
		ReturnCode rc = initializeSimplex();
		if (rc != OK) 
			return make_pair(rc, vector<double>());
		vector<double> deltas(basics. size());
		while (true) {
			// find 'entering' that c[entering] > 0
			int entering = -1;
			for (int i : nonbasics)
				if (c[i] > 0) {
					entering = i;
					break;
				}
			if (entering == -1) break;
			// find 'leaving' that minimize b[leaving]/a[leaving][entering]
			for (int i = 0; i < basics.size(); i++)
				if (a[basics[i]][entering] > 0)
					deltas[i] = b[basics[i]]/a[basics[i]][entering];
				else
					deltas[i] = std::numeric_limits<double>::infinity();
			int leavingIndex = 0;
			for (int i = 1; i < deltas.size(); i++)
				if (deltas[i] < deltas[leavingIndex])
					leavingIndex = i;
			if (deltas[leavingIndex] == std::numeric_limits<double>::infinity())
				return make_pair(UNBOUNDED, vector<double>());
			else
				pivot(basics[leavingIndex], entering);
		}
		// construct answer
		vector<double> ans(nonbasics.size());
		for (int i : nonbasics)
			if (i < nonbasics.size())
				ans[i] = 0;
		for (int i : basics)
			if (i < basics.size())
				ans[i] = b[i];
		return make_pair(OK, ans);
	}
};

std::ostream &operator<<(std::ostream &out, const SlackForm &slackform) {
	using std::endl;
	out << "object:" << endl;
	out << slackform.v;
	for (int i : slackform.nonbasics)
		out << ' '<< slackform.c[i] << "*x[" << i << ']';
	out << endl << "constraints:" << endl;
	for (int i : slackform.basics) {
		out << "x[" << i << "]=" << slackform.b[i];
		for (int j : slackform.nonbasics)
			out << ' ' << -slackform.a[i][j] << "*x[" << j << ']';
		out << endl;
	}
	return out;
}

main.cpp

#define CATCH_CONFIG_MAIN

#include "catch.hpp"
#include "linearprogram.hpp"

using namespace std;

TEST_CASE("simplex algorithm", "[simplex]") {

	SECTION("linear programming 1") {
		SlackForm slack(StandardForm({{1,1,3},{2,2,5},{4,1,2}}, {30,24,36}, {3,1,2}));
		pair<SlackForm::ReturnCode, vector<double>> ret = slack.simplex();
		REQUIRE(ret.first == SlackForm::OK);
		REQUIRE(ret.second == (vector<double>{8,4,0}));
	}

	SECTION("linear programming 2") {
		SlackForm slack(StandardForm({{1,1},{1,0},{0,1}}, {20,12,16}, {18,12.5}));
		pair<SlackForm::ReturnCode, vector<double>> ret = slack.simplex();
		REQUIRE(ret.first == SlackForm::OK);
		REQUIRE(ret.second == (vector<double>{12,8}));
	}

	SECTION("linear programming 3") {
		SlackForm slack(StandardForm({{1,-1},{2,1}}, {1,2}, {5,-3}));
		pair<SlackForm::ReturnCode, vector<double>> ret = slack.simplex();
		REQUIRE(ret.first == SlackForm::OK);
		REQUIRE(ret.second == (vector<double>{1,0}));
	}

	SECTION("linear programming 4") {
		SlackForm slack(StandardForm({{1,-1},{-1,-1},{-1,4}}, {8,-3,2}, {1,3}));
		pair<SlackForm::ReturnCode, vector<double>> ret = slack.simplex();
		REQUIRE(ret.first == SlackForm::OK);
	}

	SECTION("infeasible linear programming 1") {
		SlackForm slack(StandardForm({{1,1},{-2,-2}}, {2,-10}, {3,-2}));
		pair<SlackForm::ReturnCode, vector<double>> ret = slack.simplex();
		REQUIRE(ret.first == SlackForm::INFEASIBLE);
	}

	SECTION("infeasible linear programming 2") {
		SlackForm slack(StandardForm({{1,2},{-2,-6},{0,1}}, {4,-12,1}, {1,-2}));
		pair<SlackForm::ReturnCode, vector<double>> ret = slack.simplex();
		REQUIRE(ret.first == SlackForm::INFEASIBLE);
	}

	SECTION("unbounded linear programming 1") {
		SlackForm slack(StandardForm({{-2,1},{-1,-2}}, {-1,-2}, {1,-1}));
		pair<SlackForm::ReturnCode, vector<double>> ret = slack.simplex();
		REQUIRE(ret.first == SlackForm::UNBOUNDED);
	}

	SECTION("unbounded linear programming 2") {
		SlackForm slack(StandardForm({{-1,1},{-1,-1},{-1,4}}, {-1,-3,2}, {1,3}));
		pair<SlackForm::ReturnCode, vector<double>> ret = slack.simplex();
		REQUIRE(ret.first == SlackForm::UNBOUNDED);
	}
}