Algorithms : complete 05-advanced_algorithms_and_complexity 02-linear_programming

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: