Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
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.
#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;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment