diff options
Diffstat (limited to '05-advanced_algorithms_and_complexity/02-linear_programming/02-diet')
| -rw-r--r-- | 05-advanced_algorithms_and_complexity/02-linear_programming/02-diet/diet.cpp | 205 |
1 files changed, 196 insertions, 9 deletions
diff --git a/05-advanced_algorithms_and_complexity/02-linear_programming/02-diet/diet.cpp b/05-advanced_algorithms_and_complexity/02-linear_programming/02-diet/diet.cpp index efe4e89..7e1056a 100644 --- a/05-advanced_algorithms_and_complexity/02-linear_programming/02-diet/diet.cpp +++ b/05-advanced_algorithms_and_complexity/02-linear_programming/02-diet/diet.cpp @@ -1,21 +1,208 @@ #include <algorithm> +#include <limits> #include <iostream> #include <vector> #include <cstdio> using namespace std; +// #define DEBUG + typedef vector<vector<double>> matrix; -pair<int, vector<double>> solve_diet_problem( - int n, - int m, - matrix A, - vector<double> b, - vector<double> c) { +void print(matrix A, vector<double> b) +{ +#ifdef DEBUG + for (int i = 0; i < A.size(); i++) { + for (int j = 0; j < A[0].size(); j++) { + cerr << A[i][j] << " "; + } + cerr << ": " << b[i] << endl; + } + cerr << endl; +#endif +} + +#define DELTA 0.000001 +static inline bool same(double a, double b) +{ + double d = (a - b); + return (d < 0 ? -d : d) < DELTA; +} + +void gaussian_elimination(matrix& m) +{ + int rows = m.size(); + int cols = m[0].size(); + int np_row = 0; + + // each column + for (int col = 0; col < (cols - 1); col++) { + // each non-pivot row + for (int row = np_row; row < rows; row++) { + // leftmost non-zero + if (m[row][col] != 0) { + // swap row with np_row + if (row != np_row) + std::swap(m[np_row], m[row]); + + vector<double>& pivot_row = m[np_row]; + // scale pivot row + if (pivot_row[col] != 1.0) { + double f = 1.0 / pivot_row[col]; + for (int i = 0; i < cols; i++) + pivot_row[i] *= f; + } + + // scale then add to non-zero other rows + for (int i = 0; i < rows; i++) { + if (i != np_row && m[i][col] != 0) { + vector<double>& r = m[i]; + double f = r[col] / pivot_row[col]; + for (int j = 0; j < cols; j++) + r[j] -= pivot_row[j] * f; + } + } + + np_row++; + break; + } + } + } +} + +#define BIG_NUM 1000000000 +#define BEST -std::numeric_limits<double>::max() + +pair<int, vector<double>> solve( + matrix A, + vector<double> b, + vector<double> c) +{ + int cstrs = A.size(); // n + int vars = A[0].size(); // m + + matrix M(vars, vector<double>(vars + 1)); + double best = BEST; + vector<double> sol(vars); + + // build permutations of A + int nbits = 0; + int v = 0; + int max = (1 << cstrs) - 1; + vector<bool> bits(cstrs, false); + + bool augmented = false; // will be set to true when the last constraint is used + + // for each permutation of size vars + for (int v = 0; v < max; v++) { + // build next value + for (int i = 0; i < bits.size(); i++) { + if (!bits[i]) { + bits[i] = true; + nbits++; + break; + } + bits[i] = !bits[i]; + nbits--; + } - // Write your code here +#ifdef DEBUG + cerr << "####" << v + 1 << " : "; + for (int i = 0; i < bits.size(); i++) + cerr << bits[i] << " "; + cerr << endl; +#endif + + if (nbits != vars) + continue; + + if (!augmented && bits[bits.size() - 1]) { + // no solution found so far + if (same(best, BEST)) + return {-1, sol}; // no solution + augmented = true; + } + + // build matrix to solve + int x = 0; + for (int i = 0; i < bits.size(); i++) { + if (bits[i]) { + for (int j = 0; j < vars; j++) + M[x][j] = A[i][j]; + M[x][vars] = b[i]; + x++; + } + } + + print(M, b); + gaussian_elimination(M); + print(M, b); + + bool ok = true; + + // check that solution matrix pivots are 1's + for (int i = 0; i < vars; i++) { + if (!same(M[i][i], 1.0)) { + ok = false; + break; + } + } + if (!ok) + continue; + + // check that other inequalities holds + for (int i = 0; i < bits.size(); i++) { + if (!bits[i]) { + double r = 0; + for (int j = 0; j < vars; j++) + r += (M[j][vars] * A[i][j]); + if (r > (b[i] + DELTA)) { + ok = false; + break; + } + } + } + if (!ok) + continue; + + // compute fitness + double r = 0; + for (int j = 0; j < vars; j++) + r += (M[j][vars] * c[j]); + if (r > best) { + if (augmented) + return { 1, sol}; // infinity, see last paragraph in What To Do + best = r; + for (int j = 0; j < vars; j++) + sol[j] = M[j][vars]; + } + } + + if (same(best, BEST)) + return {-1, sol}; // no solution + + return {0, sol}; +} + +pair<int, vector<double>> solve_diet_problem( + int n, + int m, + matrix A, + vector<double> b, + vector<double> c) +{ + // add m constraints : -var <= 0 + for (int i = 0; i < m; i++) { + vector<double> v(m, 0); + v[i] = -1; + A.push_back(v); + b.push_back(0); + } + // add augmented constraint : sum vars < BIG_NUM + A.push_back(vector<double>(m, 1)); + b.push_back(BIG_NUM); - return {0, vector<double>(m, 0)}; + return solve(A, b, c); } int main(){ @@ -45,7 +232,7 @@ int main(){ case 0: printf("Bounded solution\n"); for (int i = 0; i < m; i++) { - printf("%.18f%c", ans.second[i], " \n"[i + 1 == m]); + printf("%.14f%c", ans.second[i], " \n"[i + 1 == m]); } break; case 1: |