Algorithms : complete 05-advanced_algorithms_and_complexity 02-linear_programming

8 files changed, 634 insertions, 74 deletions

diff --git a/05-advanced_algorithms_and_complexity/02-linear_programming/01-energy_values/energy_values.cpp b/05-advanced_algorithms_and_complexity/02-linear_programming/01-energy_values/energy_values.cpp
index b7956a7..253c961 100644
--- a/05-advanced_algorithms_and_complexity/02-linear_programming/01-energy_values/energy_values.cpp
+++ b/05-advanced_algorithms_and_complexity/02-linear_programming/01-energy_values/energy_values.cpp
@@ -3,7 +3,7 @@
 #include <vector>
 
 const double EPS = 1e-6;
-const int PRECISION = 20;
+const int PRECISION = 5;
 
 typedef std::vector<double> Column;
 typedef std::vector<double> Row;
@@ -19,16 +19,6 @@ struct Equation {
     Column b;
 };
 
-struct Position {
-    Position(int column, int row):
-        column(column),
-        row(row)
-    {}
-
-    int column;
-    int row;
-};
-
 Equation ReadEquation() {
     int size;
     std::cin >> size;
@@ -42,48 +32,50 @@ Equation ReadEquation() {
     return Equation(a, b);
 }
 
-Position SelectPivotElement(
-  const Matrix &a,
-  std::vector <bool> &used_rows,
-  std::vector <bool> &used_columns) {
-    // This algorithm selects the first free element.
-    // You'll need to improve it to pass the problem.
-    Position pivot_element(0, 0);
-    while (used_rows[pivot_element.row])
-        ++pivot_element.row;
-    while (used_columns[pivot_element.column])
-        ++pivot_element.column;
-    return pivot_element;
-}
-
-void SwapLines(Matrix &a, Column &b, std::vector <bool> &used_rows, Position &pivot_element) {
-    std::swap(a[pivot_element.column], a[pivot_element.row]);
-    std::swap(b[pivot_element.column], b[pivot_element.row]);
-    std::swap(used_rows[pivot_element.column], used_rows[pivot_element.row]);
-    pivot_element.row = pivot_element.column;
-}
-
-void ProcessPivotElement(Matrix &a, Column &b, const Position &pivot_element) {
-    // Write your code here
-}
-
-void MarkPivotElementUsed(const Position &pivot_element, std::vector <bool> &used_rows, std::vector <bool> &used_columns) {
-    used_rows[pivot_element.row] = true;
-    used_columns[pivot_element.column] = true;
-}
-
 Column SolveEquation(Equation equation) {
     Matrix &a = equation.a;
     Column &b = equation.b;
     int size = a.size();
-
-    std::vector <bool> used_columns(size, false);
-    std::vector <bool> used_rows(size, false);
-    for (int step = 0; step < size; ++step) {
-        Position pivot_element = SelectPivotElement(a, used_rows, used_columns);
-        SwapLines(a, b, used_rows, pivot_element);
-        ProcessPivotElement(a, b, pivot_element);
-        MarkPivotElementUsed(pivot_element, used_rows, used_columns);
+    int np_row = 0;
+
+    if (size == 0)
+        return b;
+
+    // Gaussian Elimination
+    // each column
+    for (int col = 0; col < a[0].size(); col++) {
+        // each non-pivot row
+        for (int row = np_row; row < size; row++) {
+            // leftmost non-zero
+            if (a[row][col] != 0) {
+                // swap row with np_row
+                if (row != np_row) {
+                    std::swap(a[np_row], a[row]);
+                    std::swap(b[np_row], b[row]);
+                }
+                Row& pivot_row = a[np_row];
+                // scale pivot row
+                if (pivot_row[col] != 1) {
+                    double f = 1 / pivot_row[col];
+                    for (int i = 0; i < a[0].size(); i++)
+                        pivot_row[i] *= f;
+                    b[np_row] *= f;
+                }
+                // scale then add to non-zero other rows
+                for (int i = 0; i < size; i++) {
+                    if (i != np_row && a[i][col] != 0) {
+                        Row& r = a[i];
+                        double f = r[col] / pivot_row[col];
+                        for (int j = 0; j < a[0].size(); j++)
+                            r[j] -= pivot_row[j] * f;
+                        b[i] -= b[np_row] * f;
+                    }
+                }
+
+                np_row++;
+                break;
+            }
+        }
     }
 
     return b;
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:
diff --git a/05-advanced_algorithms_and_complexity/02-linear_programming/03-SimplexMethod.pdf b/05-advanced_algorithms_and_complexity/02-linear_programming/03-SimplexMethod.pdf
new file mode 100644
index 0000000..1fe3c42
--- /dev/null
+++ b/05-advanced_algorithms_and_complexity/02-linear_programming/03-SimplexMethod.pdf
diff --git a/05-advanced_algorithms_and_complexity/02-linear_programming/03-SolvingLinearPrograms.pdf b/05-advanced_algorithms_and_complexity/02-linear_programming/03-SolvingLinearPrograms.pdf
new file mode 100644
index 0000000..f227818
--- /dev/null
+++ b/05-advanced_algorithms_and_complexity/02-linear_programming/03-SolvingLinearPrograms.pdf
diff --git a/05-advanced_algorithms_and_complexity/02-linear_programming/03-ad_allocation/ad_allocation.cpp b/05-advanced_algorithms_and_complexity/02-linear_programming/03-ad_allocation/ad_allocation.cpp
index 88bd8e6..b9566a0 100644
--- a/05-advanced_algorithms_and_complexity/02-linear_programming/03-ad_allocation/ad_allocation.cpp
+++ b/05-advanced_algorithms_and_complexity/02-linear_programming/03-ad_allocation/ad_allocation.cpp
@@ -1,42 +1,333 @@
 #include <algorithm>
+#include <cmath>
+#include <limits>
 #include <iostream>
+#include <iomanip>
 #include <vector>
 #include <cstdio>
 using namespace std;
 
-typedef vector<vector<double>> matrix;
+// #define DEBUG 1
 
-pair<int, vector<double>> allocate_ads(
-    int n,
-    int m,
-    matrix A,
-    vector<double> b,
-    vector<double> c) {
+#ifdef DEBUG
+# define say_dbg(more) cerr << "INF " << more << std::endl
+#else
+# define say_dbg(more) {}
+#endif
 
-  // Write your code here
+using my_type = long double;
+using matrix = vector<vector<my_type>>;
+using row = vector<my_type>;
 
-  return {0, vector<double>(m, 0)};
-}
+const my_type EPS = numeric_limits<my_type>::epsilon();
+
+struct position {
+    int row;
+    int column;
+    bool is_optimal() { return row == -1 || column == -1; }
+};
+
+enum class solution_state {
+  optimal,
+  infeasible,
+  unbounded
+};
+
+enum class method_phase {
+  one,
+  two
+};
+
+struct simplex_method {
+  int eqs, vars;
+  matrix A;
+  vector<my_type> rhs, obj, obj_phase1;
+  solution_state state;
+  method_phase phase;
+
+  static inline bool equals(my_type a, my_type b, my_type epsilon = 0.00001) {
+    return abs(a - b) < epsilon;
+  }
+
+  static inline bool equals_zero(my_type a, my_type epsilon = 0.001) {
+    return abs(a) < epsilon;
+  }
+
+  static inline my_type fmt(my_type v) { return (abs(v) < 0.0000005f) ? 0 : v; }
+
+  void print() const {
+#ifdef DEBUG
+    #define PR 2
+    cerr << endl;
+    for (size_t i = 0; i < A.size(); i++) {
+      for (size_t j = 0; j < A[i].size(); j++) {
+        cerr << fixed << setw(7) << setprecision(PR) << fmt(A[i][j]) << ' ';
+      }
+      cerr << fixed << " = " << setw(5) << setprecision(PR) << fmt(rhs[i]) << endl;
+    }
+    for (const auto &v : obj) {
+      cerr << setw(7) << fmt(v) << ' ';
+    }
+    cerr << fixed << " = " << setw(5) << setprecision(PR) << fmt(rhs[A.size()]) << endl;
+    if (phase == method_phase::one) {
+      for (const auto &v : obj_phase1) {
+        cerr << setw(7) << fmt(v) << ' ';
+      }
+      cerr << fixed << " = " << setw(5) << setprecision(PR) << fmt(rhs[A.size() + 1]) << endl;
+    }
+    cerr << endl;
+#endif
+  }
+
+  bool prepare() {
+    int negative_rhs = count_if(rhs.cbegin(), rhs.cend(), [](my_type v){return v < 0.0;});
+    // FIXME remove from negs the rows that have at least 1 basic variables (they won't need an artificial variable)
+    int delta = eqs - negative_rhs;
+    if (delta > 0) {
+      for (auto &v : A) {
+        v.resize(v.size() - delta);
+      }
+    }
+    for(size_t i = 0; i < negative_rhs; i++) {
+      obj.push_back(0.0);
+    }
+    // add slacks | surplus (will be inverted)
+    for (size_t i = 0; i < A.size(); i++) {
+      A[i][i + vars] = 1;
+    }
+    if (negative_rhs > 0) {
+      obj_phase1 = vector<my_type>(vars + eqs + negative_rhs, 0);
+      for (size_t i = 0; i < rhs.size(); i++) {
+        if (rhs[i] < 0.0) {
+          // invert row and rhs
+          rhs[i] = -rhs[i];
+          transform(A[i].begin(), A[i].end(), A[i].begin(), negate<my_type>());
+          // sum row to obj_phase1 and last rhs
+          for (size_t j = 0; j < obj_phase1.size(); j++) {
+            obj_phase1[j] += A[i][j];
+          }
+          rhs[rhs.size() -1] += rhs[i];
+          // add artificial variable
+          A[i][A[i].size() - negative_rhs--] = 1;
+        }
+      }
+      return true;
+    }
+    return false;
+  }
+
+  bool is_column_basic(int col) { // FIXME not used
+    bool basic = false;
+    for (size_t row = 0; row < A.size(); row++) {
+      double v = A[row][col];
+      if (!equals_zero(v, EPS)) {
+        if (!basic && equals(v, 1.0f)) {
+          basic = true;
+        } else {
+          basic = false;
+          break;
+        }
+      }
+    }
+    return basic;
+  }
+
+  position find_pivot(row &cw) {
+    // Optimality test : search largest  positive coefficient in objective function
+    // Bland's rule : Choose the leftmost nonbasic column with a negative (reduced) cost  TODO
+    int j = distance(cw.begin(), max_element(cw.begin(), cw.end()));
+    if (cw[j] <= 0.000000000001L) {                                             // FIXME
+      say_dbg("pivot col not found" << endl);
+      return {-1, -1};  // is_optimal
+    }
+
+    // Unboundedness test : search a positive coefficient
+    // Bland's rule : Chose the lowest ratio (RHS vs coefficient) where the coefficient is greater than zero.
+    // Bland's rule : If the minimum ratio is shared by several rows, choose the row with the lowest-numbered column (variable) basic in it.
+    int i = -1;
+    vector<int> pivot_rows;
+    my_type ratio = numeric_limits<my_type>::max() - 1;
+    for (size_t k = 0; k < A.size(); k++) {
+      my_type v = A[k][j];
+      if (v > 0.0) {
+        my_type r = (rhs[k] / v);
+        bool tiny =  (fabs(r - ratio) < 0.000000001L);                          // FIXME could divert with time ;) || (tiny && r > 0) ????
+        if (r < ratio || tiny) {
+          if (!tiny) {
+            ratio = r;
+            pivot_rows.clear();
+          }
+          pivot_rows.push_back(k);
+        }
+      }
+    }
+    int s = pivot_rows.size();
+    if (s == 0) {
+      state = solution_state::unbounded;
+    } else if (s > 1) {
+      say_dbg("chose 1 pivot amongs " << s);
+      size_t first = 0;
+      for (auto it = pivot_rows.begin(); it != pivot_rows.end(); it++) {
+        int r = *it;
+        for (size_t k = 0; k < A[0].size(); k++) {
+          if (equals(A[r][k], 1.0f)) {
+            if (is_column_basic(k)) {
+              if (k >= first) {  // choose the right most !!!
+                first = k;
+                i = r;
+              }
+              break;
+            }
+          }
+        }
+      }
+    } else {
+      i = pivot_rows[0];
+    }
+    say_dbg("pivot : [ " << i << " ; " << j << " ] " << setprecision(15) <<  A[i][j]<< " " << ratio);
+    return {i, j};
+  }
+
+  void scale_pivot(position p) {
+    auto d = A[p.row][p.column];                                                // optimization with 1/factor fails !!!
+    rhs[p.row] /= d;
+    for (auto &n : A[p.row]) {
+      n /= d;
+    }
+  }
+
+  void process_pivot(position p, row &w) {
+    for (int i = 0, s = A.size(); i < s; i++) {
+      if (p.row != i && !equals_zero(A[i][p.column], EPS)) {
+        my_type f = A[i][p.column];
+        for (size_t j = 0; j < A[i].size(); j++) {
+          A[i][j] -= A[p.row][j] * f;
+        }
+        rhs[i] -= rhs[p.row] * f;
+      }
+    }
+
+    if (phase == method_phase::one) {
+      auto o = obj[p.column];
+      auto o1 = obj_phase1[p.column];
+      rhs[rhs.size() - 2] -= rhs[p.row] * o;
+      rhs[rhs.size() - 1] -= rhs[p.row] * o1;
+      for (size_t i = 0; i < w.size(); i++) {
+        obj[i] -= A[p.row][i] * o;
+        obj_phase1[i] -= A[p.row][i] * o1;
+      }
+    } else {
+      auto o = obj[p.column];
+      rhs[rhs.size() - 1] -= rhs[p.row] * o;
+      for (size_t i = 0; i < w.size(); i++) {
+        obj[i] -= A[p.row][i] * o;
+      }
+    }
+  }
+
+  void Gauss_Jordam_eliminations(row &obj) {
+    for(;;) {
+      position p = find_pivot(obj);
+      if (p.is_optimal() || state == solution_state::unbounded) {
+        break;
+      }
+      scale_pivot(p);
+      process_pivot(p, obj);
+      print();
+    }
+  }
+
+  void trim_table_from_avars() {
+    obj.resize(vars + eqs);
+    rhs.pop_back();
+    for (auto &r : A) {
+      r.resize(vars + eqs);
+    }
+  }
+
+  void phase_I() {
+    say_dbg("*** Phase I" << endl);
+    phase = method_phase::one;
+    Gauss_Jordam_eliminations(obj_phase1);
+    state = equals_zero(rhs.back()) ? solution_state::optimal : solution_state::infeasible;
+  }
 
-int main(){
+  void phase_II() {
+    say_dbg("*** Phase II");
+    phase = method_phase::two;
+    trim_table_from_avars();
+    print();
+    Gauss_Jordam_eliminations(obj);
+  }
+
+  void round_values(vector<my_type> &result) {
+    for (auto it = result.begin(); it != result.end(); it++) {
+      if (equals(*it, ((int)*it)+1)) {        // slightly under int value
+        *it = (my_type)((int)*it+1);
+      } else if (equals(*it, ((int)*it))) {   // slightly over int value
+        *it = (my_type)((int)*it);
+      }
+    }
+  }
+
+  pair<int, vector<my_type> > solve() {
+    state = solution_state::optimal;
+    bool mixed_constraints = prepare();
+    print();
+    if (mixed_constraints) {
+      phase_I();
+      if (state == solution_state::infeasible) {
+        return {-1, {}};
+      }
+    }
+    phase_II();
+    if (state == solution_state::unbounded) {
+      return {1, {}};
+    }
+
+    vector<my_type> result(vars, 0);
+    for (int i = 0; i < vars; i++) {
+      int k = -1;
+      my_type sum = 0.0;
+      for (int j = 0; j < eqs; j++) {
+        sum += fabs(A[j][i]);
+        if (equals(A[j][i], 1.0)) {
+          k = j;
+        }
+      }
+      if (k != -1) {
+        result[i] = (sum - 1.0 > EPS) ? 0.0 : rhs[k];  // use RHS unless there more then 1 cumulated factor
+      }
+    }
+
+    round_values(result);
+
+    return {0, move(result)};
+  }
+
+};
+
+int main() {
   int n, m;
   cin >> n >> m;
-  matrix A(n, vector<double>(m));
+  matrix A(n, vector<my_type>(m + n + n, 0.0));  // variables[m], slacks|surplus[n], artificials[0;n]
   for (int i = 0; i < n; i++) {
     for (int j = 0; j < m; j++) {
       cin >> A[i][j];
     }
   }
-  vector<double> b(n);
+  vector<my_type> b(n + 2, 0);  // rhs + 2 objective
   for (int i = 0; i < n; i++) {
     cin >> b[i];
   }
-  vector<double> c(m);
+  vector<my_type> c(m + n);  // variables + slacks
   for (int i = 0; i < m; i++) {
     cin >> c[i];
   }
 
-  pair<int, vector<double>> ans = allocate_ads(n, m, A, b, c);
+  simplex_method sm{n, m, move(A), move(b), move(c)};
+
+  pair<int, vector<my_type>> ans = sm.solve();
 
   switch (ans.first) {
     case -1:
@@ -45,7 +336,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("%.18Lf%c", ans.second[i], " \n"[i + 1 == m]); // format must match my_type !!!
       }
       break;
     case 1:
diff --git a/05-advanced_algorithms_and_complexity/02-linear_programming/03-ad_allocation/check b/05-advanced_algorithms_and_complexity/02-linear_programming/03-ad_allocation/check
new file mode 100755
index 0000000..1f7160a
--- /dev/null
+++ b/05-advanced_algorithms_and_complexity/02-linear_programming/03-ad_allocation/check
@@ -0,0 +1,40 @@
+#! /bin/bash
+
+RESET="\033[0m"
+RED="\033[0;31m"
+GREEN="\033[0;32m"
+BROWN="\033[0;33m"
+
+GPP_OPTS="-std=c++11 -O0 -ggdb"
+#GPP_OPTS="-std=c++11 -O2"
+
+binaries=()
+
+for src in $(find . -name \*.cpp | sort)
+do
+  bin=/tmp/${src##*/}.bin
+  sed -i 's/^\/\/ #define DEBUG 1/#define DEBUG 1/' $src
+  echo -e "${RED}compile $GREEN${src}$RESET" && g++ $GPP_OPTS "$src" -o "$bin" || exit 1
+  binaries+=("$bin")
+done
+
+for t in "$@"
+do
+  tf="./tests/$t"
+  echo -en "$GREEN*****$RESET try $RED$t$RESET ...\n$BROWN>$RESET\n"
+  # cat "$tf"
+  cat "$tf.a"
+  cat "$tf.a" | sed 's/\([0-9]\+\.[0-9]\{4\}\)[0-9]\+/\1/g' > "/tmp/$t.a"
+  for bin in ${binaries[@]}
+  do
+    echo -en "$BROWN> $bin$RESET\n"
+    cat "$tf" | "$bin" 2>"${bin}-${t}-steps" | sed 's/\([0-9]\+\.[0-9]\{4\}\)[0-9]\+/\1/g' > "${bin}-${t}-out"
+    cmp "/tmp/$t.a" "${bin}-${t}-out" >/dev/null
+    if [ $? -ne 0 ]; then
+      echo -e "     $BROWN$t$RESET is ${RED}KO$RESET"
+      # cat "${bin}-${t}-out"
+    else
+      echo -e "     $BROWN$t$RESET is ${GREEN}ok$RESET"
+    fi
+  done
+done
diff --git a/05-advanced_algorithms_and_complexity/02-linear_programming/URLs b/05-advanced_algorithms_and_complexity/02-linear_programming/URLs
new file mode 100644
index 0000000..6944ff1
--- /dev/null
+++ b/05-advanced_algorithms_and_complexity/02-linear_programming/URLs
@@ -0,0 +1,3 @@
+http://reshmat.ru/simplex_method_lpp.html
+http://www.phpsimplex.com/simplex/simplex.htm?l=en
+http://cbom.atozmath.com/CBOM/Simplex.aspx
diff --git a/05-advanced_algorithms_and_complexity/02-linear_programming/check b/05-advanced_algorithms_and_complexity/02-linear_programming/check
new file mode 100755
index 0000000..b449883
--- /dev/null
+++ b/05-advanced_algorithms_and_complexity/02-linear_programming/check
@@ -0,0 +1,47 @@
+#! /bin/bash
+
+RESET="\033[0m"
+RED="\033[0;31m"
+GREEN="\033[0;32m"
+BROWN="\033[0;33m"
+
+BIN=/tmp/bin
+OUTA=/tmp/_outa
+OUTB=/tmp/_outb
+# GPP_OPTS="-std=c++11 -O0 -ggdb"
+GPP_OPTS="-std=c++11 -O2"
+
+for path in $(find . -name \*.cpp | sort); do
+    src=${path##*/}
+    dir=${path%/*}
+    echo -e "${RED}validate $BROWN$dir$RESET/$GREEN$src$RESET"
+    pushd $dir >/dev/null || exit 1
+    sed -i 's/^#define DEBUG 1/\/\/ #define DEBUG 1/' $src
+    echo -e "   ${RED}compile $GREEN$src$RESET" && g++ $GPP_OPTS $src -o $BIN || exit 1
+    if [ -d tests ]; then
+        start=`date +%s.%N`
+        echo -e "   ${RED}check $GREEN$src$RESET"
+        # if [ "$src" == "good.cpp" -o "$src" == "ad_allocation-new.cpp" ]; then
+        for t in $(find ./tests -name "*[^a~]"|sort); do
+            if [ -f $t -a -f "$t.a" ]; then
+                cat $t | $BIN 2>/dev/null | sed 's/\([0-9]\+\.[0-9]\{4\}\)[0-9]\+/\1/g' > $OUTA || echo "segfault"
+                cat $t.a  | sed 's/\([0-9]\+\.[0-9]\{4\}\)[0-9]\+/\1/g' > $OUTB
+                cmp $OUTA $OUTB >/dev/null
+                if [ $? -ne 0 ]; then
+                    echo -e "     $BROWN$t$RESET is ${RED}KO$RESET"
+                    # cat $OUTA
+                    # cat $OUTB
+                # else
+                #     echo -e "     $BROWN$t$RESET is ${GREEN}ok$RESET"
+                fi
+            fi
+        done
+        end=`date +%s.%N`
+        echo "runtime : $(echo "$end - $start" | bc -l)"
+        # fi
+    else
+        echo -e "   ${RED}no tests$RESET"
+    fi
+    popd > /dev/null
+done
+