CodeArchive Simplex

Simplex

With some feasibility code

#include <iostream>
#include <vector>
#include <string>
#include <cmath>
using namespace std;

typedef long double ld;
typedef vector < ld > vd;
typedef vector < vd > vvd;
const ld EPS = 1e-9;

void pivot(vvd& A, int r, int c) { // A is A|b
   int m = A.size(), n = A[0].size();
   ld tmp = A[r][c];
   for (int i = 0; i < n; ++i) A[r][i]/=tmp;
   for (int i = 0; i < m; ++i)
      if (i != r) {
         ld k = A[i][c];
         for(int j = 0; j < n; ++j)
            A[i][j] -= A[r][j]*k;
      }
}

void gaussian(vvd& A) {
}

vvd create_system(const vvd& coef, const vd &cons, const vd& obj) {
   vvd A(coef.size() + 1, vd(obj.size() + cons.size() + 1, 0));
   for(int i = 0; i < (int)coef.size(); i++)
      for(int j = 0; j < (int)coef[i].size(); j++)
         A[i][j] = coef[i][j];
   for(int i = 0; i < (int)obj.size(); i++) A.back()[i] = -obj[i];
   for(int i = 0; i < (int)cons.size(); i++) A[i][obj.size() + i] = 1;
   for(int i = 0; i < (int)cons.size(); i++) A[i].back() = cons[i];
   return A;
}

// Maximizes objective (SUM(var*obj))
// st. SUM(var*coef) <= constrain, var >= 0
// coef.size() = cons.size() = number of inequalities
// coef[i].size() =  obj.size() = number of variables
bool maximize(vvd& A, int skip = 0){
   int m = A.size(), n = A[0].size();

   while(true){
      int c = -1; ld mn = 0;
      for(int i = 0; i + 1 < n; i++)
         if(A.back()[i] < mn) mn = A.back()[c=i];

      if(mn > -EPS) return true;;

      int r = -1; mn = 1e99;
      for(int i = 0; i + 1 + skip < m; i++)
         if(A[i][c] > EPS && A[i].back()/A[i][c] < mn)
            mn = A[i].back()/A[r=i][c];
      if(mn > 1e98) return false; // Unbounded
      pivot(A, r, c);
   }
   return false;
}

bool feasible(vvd & A) {
   int m = A.size(), n = A[0].size();
   ld tmp;
   for (int i = 0; i < m; ++i) {
      tmp = A[i].back();
      A[i].back() = -1;
      A[i].push_back(tmp);
   }
   A.push_back(vd(n+1, 0));
   A.back()[n-1] = 1; // Objective
   int r = -1; ld mn = 0;
   for (int i = 0; i < m; ++i)
      if (A[i].back() < mn) mn = A[r=i].back();
   if (mn < -EPS) {
    pivot(A, r, n-1);
      maximize(A, 1);
      if (A.back().back() < -EPS) return false;
   }
   A.pop_back();
   for (int i = 0; i < m; ++i) {
      tmp = A[i].back();
      A[i].pop_back();
      A[i].back() = tmp;
   }
   return true;
}

ld x[128], y[128];
ld dist[128][128];

int main(){
  cout.setf(ios::fixed, ios::floatfield);
  cout.precision(9);
   int n;
   cin >> n;
   for (int i = 0; i < n; ++i) cin >> x[i] >> y[i];
   int num_cons = n * (n-1)/2;

   for (int i = 0; i < n; ++i)
      for (int j = 0; j < n; ++j)
         dist[i][j] = sqrt( (x[i] - x[j]) * (x[i] - x[j]) + (y[i] - y[j]) * (y[i] - y[j]) );

   vvd coef;
   vd cons, obj;
   //for (int i = 0; i < n; ++i) obj.push_back(-1);
   for (int i = 0; i < n; ++i) obj.push_back(-1);

   for (int i = 0; i < n; ++i)
      for (int j = i+1; j < n; ++j) {
         // Since i and j are not adjacent
         if ( ((j - i + n)%n) == 1 ) continue;
         // Add constraints
         coef.push_back(vd(n, 0));
         cons.push_back(dist[i][j] + 1e-5);
         coef.back()[i] = 1;
         coef.back()[j] = 1;
      }

   for (int i = 0; i < n; ++i) {
      int j = (i+1)%n;
      coef.push_back(vd(n,0));
      coef.back()[i] = -1;
      coef.back()[j] = -1;
    cons.push_back(-dist[i][j] + 1e-5);
  }

  coef = create_system(coef, cons, obj);

  /*
  cerr << "Before" << endl;
  for (int i = 0; i < coef.size(); ++i) {
    for (int j = 0; j < coef[i].size(); ++j)
      cerr << ' ' << coef[i][j];
    cerr << endl;
  }
  */

  //cerr << feasible(coef) << endl;
  feasible(coef);

  /*
  cerr << "After" << endl;
  for (int i = 0; i < coef.size(); ++i) {
    for (int j = 0; j < coef[i].size(); ++j)
      cerr << ' ' << coef[i][j];
    cerr << endl;
  }
  */

  for (int i = 0; i < n; ++i) {
    bool found = false;
    ld x = 0;
    int cnt = 0;
    for (int j = 0; j < coef.size(); ++j)
      if ( fabs(coef[j][i]) > EPS ) ++cnt;
    if (cnt != 1) { cout << "0.00" << endl; continue; }

    for (int j = 0; !found && j < coef.size(); ++j) {
      if (fabs(coef[j][i] - 1) < EPS) {
        cout << coef[j].back() << endl;
        break;
      }
    }
  }
  return 0;
}

Simplified version, without feasibility testing

Give praise to Yury.

#include <iostream>
#include <vector>
#include <string>
#include <cmath>
using namespace std;

typedef vector < double > vd;
typedef vector < vd > vvd;
const double EPS = 1e-9;

// Maximizes objective (SUM(var*obj))
// st. SUM(var*coef) <= constrain, var >= 0
// coef.size() = cons.size() = number of inequalities
// coef[i].size() =  obj.size() = number of variables
double simplex(const vvd &coef, const vd &cons, const vd & obj){

 vvd table(coef.size() + 1, vd(obj.size() + cons.size() + 1, 0));
 for(int i = 0; i < (int)coef.size(); i++)
   for(int j = 0; j < (int)coef[i].size(); j++)
     table[i][j] = coef[i][j];
 for(int i = 0; i < (int)obj.size(); i++) table.back()[i] = -obj[i];
 for(int i = 0; i < (int)cons.size(); i++) table[i][obj.size() + i] = 1;
 for(int i = 0; i < (int)cons.size(); i++) table[i].back() = cons[i];

 while(true){
   int c = -1; double mn = 0;
   for(int i = 0; i + 1 < (int)table.back().size(); i++)
     if(table.back()[i] < mn) mn = table.back()[c=i];

   if(mn > -EPS) break;

   int r = -1; mn = 1e99;
   for(int i = 0; i + 1 < (int)table.size(); i++)
     if(table[i][c] > EPS && table[i].back()/table[i][c] < mn)
       mn = table[i].back()/table[r=i][c];

   if(mn > 1e98) break;

   double temp = table[r][c];
   for(int i = 0; i < (int)table[r].size(); i++)
     table[r][i] /= temp;

   for(int i = 0; i < (int)table.size(); i++)
     if(i != r){
       double k = table[i][c]/table[r][c];
       for(int j = 0; j < (int)table[i].size(); j++)
         table[i][j] -= table[r][j]*k;
     }
 }
 return table.back().back();
}

int main(){


 return 0;
}

-- Main.MatthewChan - 08 Dec 2005