diff options
Diffstat (limited to 'Master/texmf/asymptote/simplex.asy')
-rw-r--r-- | Master/texmf/asymptote/simplex.asy | 272 |
1 files changed, 272 insertions, 0 deletions
diff --git a/Master/texmf/asymptote/simplex.asy b/Master/texmf/asymptote/simplex.asy new file mode 100644 index 00000000000..bd9f69c9324 --- /dev/null +++ b/Master/texmf/asymptote/simplex.asy @@ -0,0 +1,272 @@ +/***** + * simplex.asy + * Andy Hammerlindl 2004/07/27 + * + * Solves the two-variable linear programming problem using the simplex method. + * This problem is specialized in that the second variable, "b", does not have + * a non-negativity condition, and the first variable, "a", is the quantity + * being maximized. + * Correct execution of the algorithm also assumes that the coefficient of "b" + * will be +1 or -1 in every added restriction, and that the problem can be + * initialized to a valid state by pivoting b with one of the slack + * variables. This assumption may in fact be incorrect. + *****/ + +private real infinity=sqrt(0.25*realMax); + +struct problem { + typedef int var; + static var VAR_A = 0; + static var VAR_B = 1; + + static int OPTIMAL = -1; + static var UNBOUNDED = -2; + static int INVALID = -3; + + struct row { + real c, t[]; + } + + // The variables of the rows. + // Initialized for the two variable problem. + var[] v = {VAR_A, VAR_B}; + + // The rows of equalities. + row rowA() { + row r = new row; + r.c = 0; + r.t = new real[] {1, 0}; + return r; + } + row rowB() { + row r = new row; + r.c = 0; + r.t = new real[] {0, 1}; + return r; + } + row[] rows = {rowA(), rowB()}; + + // The number of original variables. + int n = rows.length; + + // Pivot the variable v[col] with vp. + void pivot(int col, var vp) + { + var vc=v[col]; + + // Recalculate rows v[col] and vp for the pivot-swap. + row rvc = rows[vc], rvp = rows[vp]; + real factor=1/rvp.t[col]; // NOTE: Handle rvp.t[col] == 0 case. + rvc.c=-rvp.c*factor; + rvp.c=0; + rvc.t=-rvp.t*factor; + rvp.t *= 0; + rvc.t[col]=factor; + rvp.t[col]=1; + + var a=min(vc,vp); + var b=max(vc,vp); + + // Recalculate the rows other than the two used for the above pivot. + for (var i = 0; i < a; ++i) { + row r=rows[i]; + real m = r.t[col]; + r.c += m*rvc.c; + r.t += m*rvc.t; + r.t[col]=m*factor; + } + for (var i = a+1; i < b; ++i) { + row r=rows[i]; + real m = r.t[col]; + r.c += m*rvc.c; + r.t += m*rvc.t; + r.t[col]=m*factor; + } + for (var i = b+1; i < rows.length; ++i) { + row r=rows[i]; + real m = r.t[col]; + r.c += m*rvc.c; + r.t += m*rvc.t; + r.t[col]=m*factor; + } + + // Relabel the vars. + v[col] = vp; + } + + // As b does not have a non-negativity condition, it must initially be + // pivoted out for a variable that does. This selects the initial + // variable to pivot with b. It also assumes that there is a valid + // solution with a == 0 to the linear programming problem, and if so, it + // picks a pivot to get to that state. In our case, a == 0 corresponds to + // a picture with the user coordinates shrunk down to zero, and if that + // doesn't fit, nothing will. + var initVar() + { + real min=infinity, max=-infinity; + var argmin=0, argmax=0; + + for (var i = 2; i < rows.length; ++i) { + row r=rows[i]; + if (r.t[VAR_B] > 0) { + real val=r.c/r.t[VAR_B]; + if (val < min) { + min=val; + argmin=i; + } + } else if (r.t[VAR_B] < 0) { + real val=r.c/r.t[VAR_B]; + if (val > max) { + max=val; + argmax=i; + } + } + } + + // If b has a minimal value, choose a pivot that will give b its minimal + // value. Otherwise, if b has maximal value, choose a pivot to give b its + // maximal value. + return argmin != 0 ? argmin : + argmax != 0 ? argmax : + UNBOUNDED; + } + + // Initialize the linear program problem by moving into an acceptable state + // this assumes that b is unrestrained and is the second variable. + // NOTE: Works in limited cases, may be bug-ridden. + void init() + { + // Find the lowest constant term in the equations. + var lowest = 0; + for (var i = 2; i < rows.length; ++i) { + if (rows[i].c < rows[lowest].c) + lowest = i; + } + + // Pivot if necessary. + if (lowest != 0) + pivot(VAR_B, lowest); + } + + // Selects a column to pivot on. Returns OPTIMAL if the current state is + // optimal. Assumes we are optimizing the first row. + int selectColumn() + { + int i=find(rows[0].t > 0,1); + return (i >= 0) ? i : OPTIMAL; + } + + // Select the new variable associated with a pivot on the column given. + // Returns UNBOUNDED if the space is unbounded. + var selectVar(int col) + { + // We assume that the first two vars (a and b) once swapped out, won't be + // swapped back in. This finds the variable which gives the tightest + // non-negativity condition restricting our optimization. This turns + // out to be the max of c/t[col]. Note that as c is positive, and + // t[col] is negative, all c/t[col] will be negative, so we are finding + // the smallest in magnitude. + var vp=UNBOUNDED; + real max=-infinity; + for (int i = 2; i < rows.length; ++i) { + row r=rows[i]; + if(r.c < max*r.t[col]) { + max=r.c/r.t[col]; vp=i; + } + } + + return vp; + } + + // Checks that the rows are in a valid state. + bool valid() + { + // Checks that constants are valid. + bool validConstants() { + for (int i = 0; i < rows.length; ++i) + if (rows[i].c < 0) + return false; + return true; + } + + // Check a variable to see if its row is simple. + // NOTE: Simple rows could be optimized out, since they are not really + // used. + bool validVar(int col) { + + var vc = v[col]; + row rvc = rows[vc]; + + if (rvc.c != 0) + return false; + for (int i = 0; i < n; ++i) + if (rvc.t[i] != (i == col ? 1 : 0)) + return false; + + return true; + } + + if (!validConstants()) { + return false; + } + for (int i = 0; i < n; ++i) + if (!validVar(i)) { + return false; + } + + return true; + } + + + // Perform the algorithm to find the optimal solution. Returns OPTIMAL, + // UNBOUNDED, or INVALID (if no solution is possible). + int optimize() + { + // Put into a valid state to begin and pivot b out. + var iv=initVar(); + if (iv == UNBOUNDED) + return iv; + pivot(VAR_B, iv); + + if (!valid()) + return INVALID; + + while(true) { + int col = selectColumn(); + + if (col == OPTIMAL) + return col; + var vp = selectVar(col); + + if (vp == UNBOUNDED) + return vp; + + pivot(col, vp); + } + + // Shouldn't reach here. + return INVALID; + } + + // Add a restriction to the problem: + // t1*a + t2*b + c >= 0 + void addRestriction(real t1, real t2, real c) + { + row r = new row; + r.c = c; + r.t = new real[] {t1, t2}; + rows.push(r); + } + + // Return the value of a computed. + real a() + { + return rows[VAR_A].c; + } + + // Return the value of b computed. + real b() + { + return rows[VAR_B].c; + } +} |