CodeArchive Simplex
From Ubcacm
Revision as of 22:06, 29 January 2007 by 142.103.22.72 (Talk)
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