Skip to content

Instantly share code, notes, and snippets.

# avitevet/diet.cpp Created Apr 6, 2017

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 #include #include #include #include using namespace std; typedef double ElementType; typedef vector> 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 &nonbasics, vector &basics, matrix &A, vector &b, vector &c, ElementType &v, int l, int e) { // in CLRS, this is N - {e} vector 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 nonbasics, basics; matrix A; vector b, c; ElementType v; } SlackForm; template std::ostream &operator<<(ostream &out, const vector &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 &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> SimplexIterations(SlackForm &slack) { vector &nonbasics = slack.nonbasics, &basics = slack.basics; matrix &A = slack.A; vector &b = slack.b, &c = slack.c; ElementType &v = slack.v; // this implements lines 2 - 17 of SIMPLEX from CLRS vector delta(basics.size(), std::numeric_limits::infinity()); for (vector::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::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::max(); ElementType min_element = std::numeric_limits::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::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 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 InitializeSimplex(const matrix &A, const vector &b, const vector &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::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 &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 &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> 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[e] != 0. int e = -1; for (int n : aux.nonbasics) { if (aux.A[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 &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> solve_diet_problem( int n, int m, matrix A_initial, vector b_initial, vector c_initial) { // Using simplex algorithm as described in CLRS pair 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(m)); vector b(n); vector c(m); pair> 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; }
to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.