Using linear programming/simplex method to solve the diet problem: given a set of foods with nutritional data, and a set of nutritional constraints, find the cheapest combination of foods that meets your nutritional needs.

diet.cpp

#include <algorithm>
#include <iostream>
#include <vector>
#include <cstdio>
#include <cassert>
using namespace std;

typedef double ElementType;
typedef vector<vector<ElementType>> matrix;

const int NO_SOLUTION = -1;
const int BOUNDED_SOLUTION = 0;
const int INFINITY_SOLUTION = 1;

const ElementType C_THRESHOLD = 0.0000001;
const ElementType DELTA_RATIO_THRESHOLD = 0; //  0.0001;
const ElementType DELTA_COMPARE_EPSILON = 0.0000001;

// implements PIVOT from CLRS
void pivot(
	vector<int> &nonbasics,
	vector<int> &basics,
	matrix &A,
	vector<ElementType> &b,
	vector<ElementType> &c,
	ElementType &v,
	int l,
	int e) {

	// in CLRS, this is N - {e}
	vector<int> otherNonbasics;
	for (int n : nonbasics) {
		if (n != e) {
			otherNonbasics.push_back(n);
		}
	}

	// the e & l variables provide the indices of the entering and leaving variables, but
	// they are not the same as the row in A that will be processed.  row_l provides that mapping
	// (aka, the row in A that currently represents the basic variable x[l])
	int row_l = -1;
	for (size_t i = 0; i < basics.size(); ++i) {
		if (basics[i] == l) {
			row_l = i;
			break;
		}
	}

	// compute the coefficients of the equation for new basic variable x[e].
	// in other words, we are solving for x[e] using the constraint indexed by l.
	// to do so, we scale the constraint by dividing by A[l][e], which sets the coefficient
	// in A in constraint l for x[e]  to 1.0 and effectively solves for it
	b[row_l] = b[row_l] / A[row_l][e];
	for (int j : otherNonbasics) {
		A[row_l][j] = A[row_l][j] / A[row_l][e];
	}
	A[row_l][l] = 1.0 / A[row_l][e];
	A[row_l][e] = 0.0;

	// compute the coefficients of the remaining constraints.
	// Effectively, this updates the constraints not indexed by l
	// by substituting the RHS of equation for x[e] into each constraint
	// for each occurrence of x[e]
	for (size_t i = 0; i < basics.size(); ++i) {
		if (i == row_l) {
			continue;
		}

		b[i] -= A[i][e] * b[row_l];
		for (int j : otherNonbasics) {
			A[i][j] -= A[i][e] * A[row_l][j];
		}

		A[i][l] = -A[i][e] * A[row_l][l];
		A[i][e] = 0.0;
	}

	// compute the objective function
	// by substituting in constraint l (solved for x[e])
	v += c[e] * b[row_l];
	for (int j : otherNonbasics) {
		c[j] -= c[e] * A[row_l][j];
	}
	c[l] = -c[e] * A[row_l][l];
	c[e] = 0.0;

	// compute new sets of basic & nonbasic variables by swapping e & l
	for (size_t n = 0; n < nonbasics.size(); ++n) {
		if (nonbasics[n] == e) {
			nonbasics[n] = l;
			break;
		}
	}
	for (size_t n = 0; n < basics.size(); ++n) {
		if (basics[n] == l) {
			basics[n] = e;
			break;
		}
	}

	return;
}

typedef struct SlackForm {
	vector<int> nonbasics, basics;
	matrix A;
	vector<ElementType> b, c;
	ElementType v;
} SlackForm;

template<typename T>
std::ostream &operator<<(ostream &out, const vector<T> &v) {
	out << "{ ";
	for (size_t i = 0; i < v.size(); ++i) {
		out << v[i];
		if (i != v.size() - 1) {
			out << ", ";
		}
	}
	out << " }";

	return out;
}

std::ostream &operator<<(ostream &out, SlackForm &s) {
	out << "N = " << s.nonbasics << endl;
	out << "B = " << s.basics << endl;
	out << "A = { ";
	for (vector<ElementType> &row : s.A) {
		out << row << endl;
	}
	out << " }" << endl;
	out << "b = " << s.b << endl;
	out << "c = " << s.c << endl;
	out << "v = " << s.v << endl;
	return out;
}

pair<int, vector<ElementType>> SimplexIterations(SlackForm &slack) {
	vector<int>
		&nonbasics = slack.nonbasics,
		&basics = slack.basics;
	matrix &A = slack.A;
	vector<ElementType>
		&b = slack.b,
		&c = slack.c;
	ElementType &v = slack.v;

	// this implements lines 2 - 17 of SIMPLEX from CLRS

	vector<ElementType> delta(basics.size(), std::numeric_limits<ElementType>::infinity());

	for (vector<ElementType>::iterator j = std::max_element(c.begin(), c.end());
		*j > C_THRESHOLD;
		j = std::max_element(c.begin(), c.end())) {
		int e = std::distance(c.begin(), j);

		// choose l, the index of the variables that will be the leaving variable.
		// do this by seeing which of the constraints allows the largest value of
		// x[e] while not violating the non-negativity constraints on x
		for (size_t i = 0; i < basics.size(); ++i) {
			delta[i] = (A[i][e] > DELTA_RATIO_THRESHOLD) ? b[i] / A[i][e] : std::numeric_limits<ElementType>::infinity();
		}

		// now find "the smallest" delta, but there is a fudge factor for "ties".  If delta[i] =~ delta[j] then choose l = min(basics[i], basics[j])
		int l = std::numeric_limits<int>::max();
		ElementType min_element = std::numeric_limits<ElementType>::infinity();
		for (size_t i = 0; i < delta.size(); ++i) {
			bool lt = delta[i] < min_element - DELTA_COMPARE_EPSILON;
			bool appxEq = (min_element - DELTA_COMPARE_EPSILON <= delta[i]) && (delta[i] < min_element + DELTA_COMPARE_EPSILON);
			if (lt) {
				min_element = delta[i];
				l = basics[i];
			}
			else if (appxEq && (basics[i] < l)) {
				min_element = delta[i];
				l = basics[i];
			}
		}

		if (min_element == std::numeric_limits<ElementType>::infinity()) {
			return{ INFINITY_SOLUTION, {} };
		}

		pivot(nonbasics, basics, A, b, c, v, l, e);
		
		//cout << "POST-pivot on (l, e) = (" << l << ", " << e << "): " << endl << slack;
	}

	// form p, the optimal vertex.  The nonbasics are stored in 0..(c.size()), and
	// the basics are stored in c.size..c.size+b.size()-1.
	int n = c.size();
	vector<ElementType> p(n, 0);

	for (int i = 0; i < n; ++i) {
		// if i is in basics:
		auto it = std::find(basics.begin(), basics.end(), i);
		if (it != basics.end()) {
			int row = std::distance(basics.begin(), it);
			p[i] = b[row];
		}
		else {
			p[i] = 0;
		}
	}

	return{ BOUNDED_SOLUTION, p };
}

// initialize-simplex from CLRS
pair<int, SlackForm> InitializeSimplex(const matrix &A, const vector<ElementType> &b, const vector<ElementType> &c) {
	// So in this function, we are given a program in standard form: 
	//   maximize     cx
	//   subject to   Ax <= b; x >= 0
	// and we must convert it to slack form:
	//   maximize     cx
	//   subject to   b - Ax = s; x,s >= 0
	// 
	// A is given as an m x n matrix, where m = number of x variables; n = number of linear inequalities.
	//   |c| = |x| = m, and |b| = n
	//
	// In slack form,
	//   |s| = |b| = n
	//
	// To convert to slack form, A must be converted to a (m + s) x n matrix.  Also c must be changed so
	//   |c| = |x| + |s|


	// if all the elements of b are non-negative, then just return 
	// the inputs, expanded so they include the slack variables
	vector<ElementType>::const_iterator min_b = std::min_element(b.begin(), b.end());
	if (*min_b >= 0) {
		SlackForm ret;

		ret.A = A;
		ret.b = b;
		ret.c = c;
		ret.v = 0;
		
		// c, and each row of A, need to be expanded to include the slack variables.  
		// The slack variable values will be initialized to 0
		ret.c.resize(ret.c.size() + ret.b.size(), 0.0);
		for (vector<ElementType> &row : ret.A) {
			row.resize(row.size() + ret.b.size(), 0.0);
		}

		// nonbasics = {0, 1, .., n - 1}
		for (size_t i = 0; i < c.size(); ++i) {
			ret.nonbasics.push_back(i);
		}
		// basics = {n, n + 1, .., n + m -1}
		for (size_t i = 0; i < b.size(); ++i) {
			ret.basics.push_back(c.size() + i);
		}

		return{ BOUNDED_SOLUTION, ret };
	}

	// if there is a negative value in b, then we need to correct that by introducing yet another variable.
	// In CLRS, the variable is nonbasic x0, but
	// rather than using x0 as in CLRS, I'll use xN, because I'm using 0 based arrays everywhere else already.
	// Form Laux by adding -xN to the lhs of each constraint and setting objective function to -1.0 * xN.
	// Also expand c & each row of A to include the slack variables, which are initialized to 0.0
	SlackForm aux;
	aux.A = A;
	for (vector<ElementType> &row : aux.A) {
		row.push_back(-1.0);
		row.resize(row.size() + b.size(), 0.0);
	}
	aux.b = b;
	// c[0..n-1] = 0.0, c[n] = -1, c[n+1..n+m-1] = 0.0
	aux.c.resize(c.size() + b.size() + 1, 0.0);
	aux.c[c.size()] = -1.0;

	// aux.nonbasics = {0, 1, .., n} - the original n nonbasic variables + xN
	for (size_t i = 0; i < c.size() + 1; ++i) {
		aux.nonbasics.push_back(i);
	}

	// aux.basics = {n + 1, n + 2, .., n + m} - only the original basic variables
	for (size_t i = 0; i < b.size(); ++i) {
		aux.basics.push_back(c.size() + 1 + i);
	}
	aux.v = 0.0;

	int l = c.size() + 1 + std::distance(b.begin(), min_b);

	pivot(aux.nonbasics, aux.basics, aux.A, aux.b, aux.c, aux.v, l, c.size());
	// the basic solution is now feasible for Laux

	pair<int, vector<ElementType>> solution = SimplexIterations(aux);
	if (solution.first != BOUNDED_SOLUTION) {
		return{ solution.first, SlackForm() };
	}

	// if optimal solution to Laux does not set xN == 0, this is infeasible
	if (solution.second[c.size()] == 0) {
		// if xN is basic, ie if N is a member of aux.basics
		if (std::find(aux.basics.begin(), aux.basics.end(), c.size()) != aux.basics.end()) {
			// perform a degenerate pivot to make it nonbasic.  The leaving variable l
			// is N from xN (== c.size()).  The entering variable e can be any of the nonbasics
			// such that A[0][e] != 0.
			int e = -1;
			for (int n : aux.nonbasics) {
				if (aux.A[0][n] != 0) {
					e = n;
					break;
				}
			}
			pivot(aux.nonbasics, aux.basics, aux.A, aux.b, aux.c, aux.v, c.size(), e);
		}

		// since xN is 0, it can be removed
		
		// remove xN from the constraints
		for (vector<ElementType> &row : aux.A) {
			row.erase(row.begin() + c.size());
		}
		// remove xN from nonbasics
		aux.nonbasics.erase(std::remove(aux.nonbasics.begin(), aux.nonbasics.end(), c.size()), aux.nonbasics.end());
		// adjust the indices of the basic & nonbasic vars
		for (int &b : aux.basics) {
			if (b > c.size()) {
				b--;
			}
		}
		for (int &nb : aux.nonbasics) {
			if (nb > c.size()) {
				nb--;
			}
		}

		// remove xN from c - this will be done in the next step

		// restore original objective function of L, and increase its size 
		// to accommodate the |b| slack variables
		aux.c = c;
		aux.c.resize(c.size() + b.size(), 0.0);


		// Restore the original objective function L, but
		// replace each basic variable in this objective function by the RHS of its associated constraint
		for (size_t i = 0; i < aux.c.size(); ++i) {
			// if the coefficient of the variable in the objective function is not zero,
			// and it's a basic variable, then it needs to be replaced with the RHS of its associated constraint
			ElementType objVarCoefficient = aux.c[i];
			auto ib = std::find(aux.basics.begin(), aux.basics.end(), i);
			bool isBasic = ib != aux.basics.end();
			if ((objVarCoefficient != 0.0) && isBasic) {
				size_t basicRow = std::distance(aux.basics.begin(), ib);
				for (size_t j = 0; j < aux.c.size(); ++j) {
					// use -= here because when we replace the basic var in the objective function with the RHS of its
					// associated constraint, we must remember that the form is b - Ax = s
					aux.c[j] -= objVarCoefficient * aux.A[basicRow][j];
				}
				aux.c[i] = 0.0;
				aux.v += objVarCoefficient * aux.b[basicRow];
			}
		}

		return{ BOUNDED_SOLUTION, aux };

	}

	return{ NO_SOLUTION, SlackForm() };
}

pair<int, vector<ElementType>> solve_diet_problem(
    int n, 
    int m, 
    matrix A_initial, 
    vector<ElementType> b_initial, 
    vector<ElementType> c_initial) {

	// Using simplex algorithm as described in CLRS
	pair<int, SlackForm> simplexValues = InitializeSimplex(A_initial, b_initial, c_initial);

	if (simplexValues.first != BOUNDED_SOLUTION) {
		return { simplexValues.first, {} };
	}

	return SimplexIterations(simplexValues.second);
}

// written for https://www.coursera.org/learn/advanced-algorithms-and-complexity, week 2, based on the simplex method as described in 
// Introduction to Algorithms, CLRS.  Data input code has been removed as minimal protection against use by unscrupulous students.
int main(){
  const int n = 10, m = 20;
  matrix A(n, vector<ElementType>(m));
  vector<ElementType> b(n);
  vector<ElementType> c(m);

  pair<int, vector<ElementType>> ans = solve_diet_problem(n, m, A, b, c);

  switch (ans.first) {
    case NO_SOLUTION: 
      printf("No solution\n");
      break;
    case BOUNDED_SOLUTION: 
      printf("Bounded solution\n");
      for (int i = 0; i < m; i++) {
        printf("%.18f%c", ans.second[i], " \n"[i + 1 == m]);
      }
      break;
    case INFINITY_SOLUTION:
      printf("Infinity\n");
      break;      
  }
  return 0;
}