Artifact fb49c6d742de6e37ae32eb41e0835657d631b535f131c9f348231bf34a8743f4:
- Executable file
r38/packages/linalg/simplex.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 26241) [annotate] [blame] [check-ins using] [more...]
module simplex; % added a patch by Thomas Sturm and Andreas Dolzmann for Simplex1 % WN 21-Oct-1999 %**********************************************************************% % % % Computation of the optimal value of an objective function given a % % number of linear inequalities using the SIMPLEX algorithm. % % % % Author: Matt Rebbeck, June 1994. % % % % Many of the ideas were taken from "Linear Programming" by % % M.J.Best & K. Ritter % % % % Minor changes: Herbert Melenk, Jan 1995 % % % % replacing first, second etc. by car, cadr % % converted big smacros to ordinary procedures % % % %**********************************************************************% if not get('leq,'simpfn) then << algebraic operator <=; algebraic operator >=; >>; symbolic smacro procedure smplx_prepsq u; % % If u in (!*sq) standard quotient form then get !:rd!: part. % if atom u then u else if car u = 'minus and atom cadr u then u else if car u = 'minus and caadr u = '!*sq then {'minus,car cadadr u} else if car u = 'minus and caadr u = '!:rd!: then u else if car u = '!:rd!: then u else if car u = '!*sq then prepsq cadr u; symbolic smacro procedure fast_row_dim(in_mat); % % Finds row dimension of a matrix with no error checking. % length in_mat #- 1; symbolic smacro procedure fast_column_dim(in_mat); % % Finds column dimension of a matrix with no error checking. % length cadr in_mat; symbolic smacro procedure fast_stack_rows(in_mat,row_list); % % row_list is always an integer in simplex. % 'mat.{nth(cdr in_mat,row_list)}; symbolic smacro procedure fast_getmat(matri,i,j); % % Get matrix element (i,j). % fast_unchecked_getmatelem list(matri,i,j); symbolic smacro procedure fast_my_letmtr(u,v,y); rplaca(pnth(nth(cdr y,car my_revlis cdr u),cadr my_revlis cdr u),v); put('simplex,'psopfn,'simplex1); symbolic procedure simplex1(input); % % The simplex problem is: % % min {c'x | Ax = b, x>=0}, % % where Ax = b describes the linear equations and c is the function % that is to be optimised. % % This code implements the basic phaseI-phaseII revised simplex % algorithm. (phaseI checks for feasibility and phaseII finds the % optimal solution). % % A general idea of tha algorithm is as follows: % % 1: create phase 1 data. % % Add slack and artificial variables to equations to create matrix % A1. The initial basis (ib) consists of the numbers of the columns % relating to the artificial variables. The basic feasible solution % (xb) (if one exists) is B^(-1)*b where b is the r.h.s. of the % equations. Throughout, cb denotes the columns of the objective % matrix corresponding to ib. % This data goes to the revised simplex algorithm(2). % % 2: revised simplex: % % step 1: Computation of search direction sb. % % Compute u = -(B^(-1))'*cb, the smallest index k, and vk s.t. % % vk = min{c(i) + A(i)'u | i not elt of ib}. % % If vk>=0, stop with optimal solution xb = B^(-1)*b. % If vk<0, set sb = B^(-1)*A(k) and go to step 2. % % step 2: Computation of maximum feasible step size Ob. % % If sb<=0 then rederr "Problem unbounded from below". % If (sb)(i) >0 for at least one i, compute Ob and the smallest % index l s.t. % % (xb)(l) { (xb)(i) | } % Ob = ------- = min { ------ | all i with (sb)(i)>0 }, % (sb)(l) { (sb)(i) | } % % and go to step 3. % % step3: Update. % % Replace B^(-1) with [phiprm((B^(-1)',A(k),l)]', xb with B^(-1)*b % and the l'th elt in ib with k. Compute cb'xb and go to step 1. % % 3: If we get this far (ie: we are feasible) then apply revised % simplex to A (equations with slacks added), and the new xb, % Binv, and ib. % % % Further details and far more advanced algorithms can be found % in almost any linear programming book. The above was adapted from % "Linear Programming" M.J.Best and K. Ritter. To date, the code % contains no anti_cycling algorithm. % begin scalar max_or_min,objective,equation_list,tmp,A,b,obj_mat,X,A1, phase1_obj,ib,xb,Binv,simp_calc,phase1_obj_value,big,sum, stop,work,I_turned_rounded_on,ans_list,optimal_value; integer m,n,k,i,ell,no_coeffs,no_variables,equation_variables; max_or_min := reval car input; objective := reval cadr input; equation_list := normalize!-equationl reval caddr input; % <--- patch if not !*rounded then << I_turned_rounded_on := t; on rounded; >>; if (max_or_min neq 'max) and (max_or_min neq 'min) then rederr "Error in simplex(first argument): must be either max or min."; if pairp equation_list and car equation_list = 'list then equation_list := cdr equation_list else rederr "Error in simplex(third argument}: must be a list."; % get rid of any two (or more) equal equations! Probably makes % things faster in general. tmp := unique_equation_list(equation_list); equation_list := car tmp; equation_variables := cadr tmp; % If r.h.s. and l.h.s. of inequalities are both <0 then multiply % by -1. equation_list := make_equations_positive(equation_list); % If there are variables in the objective function that are not in % the equation_list then add these variables to the equation list. % (They are added as variable>=0). equation_list := add_not_defined_variables(objective,equation_list, equation_variables); tmp := initialise(max_or_min,objective,equation_list); A := car tmp; b := cadr tmp; obj_mat := caddr tmp; X := cadddr tmp; % no_variables is the number of variables in the input equation % list. this is used in make_answer_list. no_variables := car cddddr tmp; % r.h.s. (i.e. b matrix) must be positive. tmp := check_minus_b(A,b); A := car tmp; b := cadr tmp; m := fast_row_dim(A); n := no_coeffs := fast_column_dim(A); tmp := create_phase1_A1_and_obj_and_ib(A); A1 := car tmp; phase1_obj := cadr tmp; ib := caddr tmp; xb := copy_mat(b); Binv := fast_make_identity(fast_row_dim(A)); % Phase 1 data is now ready to go... simp_calc := simplex_calculation(phase1_obj,A1,b,ib,Binv,xb); phase1_obj_value := car simp_calc; xb := cadr simp_calc; Binv := cadddr(simp_calc); if get_num_part(phase1_obj_value) neq 0 then rederr "Error in simplex: Problem has no feasible solution."; % Are any artificials still basic? for ell:=1:m do if nth(ib,ell) <= n then <<>> else << % so here, an artificial is basic in row ell. big := -1; k := 0; stop := nil; i := 1; while i<=n and not stop do << sum := get_num_part(smplx_prepsq(fast_getmat(reval {'times,fast_stack_rows(Binv,ell), fast_augment_columns(A,i)},1,1))); if abs(sum) leq big then stop := t else << big := abs(sum); k := i; >>; i := i+1; >>; if big geq 0 then <<>> else rederr {"Error in simplex: constraint",k," is redundant."}; work := fast_augment_columns(A,k); Binv := phiprm(Binv,work,ell); nth(ib,ell) := k; >>; % Into Phase 2. simp_calc := simplex_calculation(obj_mat,A,b,ib,Binv,xb); optimal_value := car simp_calc; xb := cadr simp_calc; ib := caddr simp_calc; ans_list := make_answer_list(xb,ib,no_coeffs,X,no_variables); if I_turned_rounded_on then off rounded; if max_or_min = 'max then optimal_value := my_reval{'minus,optimal_value}; return {'list,optimal_value,'list.ans_list}; end; flag('(simplex1),'opfn); procedure normalize!-equationl(eql); if eqcar(eql , 'list) then 'list . for each equ in cdr eql collect normalize!-equation equ else eql; procedure normalize!-equation(equ); begin scalar lhs,b; lhs := numr subtrsq(simp cadr equ,simp caddr equ); b := negf abssummand lhs; return {car equ,prepf addf(lhs,b),prepf b} end; procedure abssummand(f); if domainp f then f else if red f then abssummand red f; symbolic procedure unique_equation_list(equation_list); % % Removes repititions in input. Also returns coeffecients in equation % list. % begin scalar unique_equation_list,coeff_list; for each equation in equation_list do << if not intersection({equation},unique_equation_list) then << unique_equation_list := append(unique_equation_list,{equation}); coeff_list := union(coeff_list,get_coeffs(cadr equation)); >>; >>; return {unique_equation_list,coeff_list}; end; symbolic procedure make_equations_positive(equation_list); % % If r.h.s. and l.h.s. of inequality are <0 then mult. both sides by % -1. % for each equation in equation_list collect if pairp cadr equation and caadr equation = 'minus and pairp caddr equation and caaddr equation = 'minus then {car equation,my_minus(cadr equation),my_minus(caddr equation)} else equation; symbolic procedure add_not_defined_variables (objective,equation_list,equation_variables); % % If variables in the objective have not been defined in the % inequalities(equation_list) then add them. They are added as % variable >= 0. % begin scalar obj_variables; obj_variables := get_coeffs(objective); if length obj_variables = length equation_variables then return equation_list; for each variable in obj_variables do << if not intersection({variable},equation_variables) then << prin2 "*** Warning: variable "; prin2 variable; prin2t " not defined in input. Has been defined as >=0."; equation_list := append(equation_list,{{'geq,variable,0}}); >>; >>; return equation_list; end; symbolic procedure initialise(max_or_min,objective,equation_list); % % Creates A (with slack variables included), b (r.h.s. of equations), % the objective matrix (obj_mat) and X s.t. AX=b and % obj_mat * X = objective function. % Also returns the number of equations in the equation_list so we know % where to stop making answers in make_answer_list. % begin scalar more_init,A,b,obj_mat,X; integer no_variables; if max_or_min = 'max then objective := reval{'times,objective,-1}; more_init := more_initialise(objective,equation_list); A := car more_init; b := cadr more_init; obj_mat := caddr more_init; X := cadddr more_init; no_variables := car cddddr more_init; return {A,b,obj_mat,X,no_variables}; end; symbolic procedure more_initialise(objective,equation_list); begin scalar objective,equation_list,non_slack_variable_list,obj_mat, no_of_non_slacks,tmp,variable_list,slack_equations,A,b,X; non_slack_variable_list := get_preliminary_variable_list(equation_list); no_of_non_slacks := length non_slack_variable_list; tmp := add_slacks_to_equations(equation_list); slack_equations := car tmp; b := cadr tmp; variable_list := union(non_slack_variable_list,caddr tmp); tmp := get_X_and_obj_mat(objective,variable_list); X := car tmp; obj_mat := cadr tmp; A := simp_get_A(slack_equations,variable_list); return {A,b,obj_mat,X,no_of_non_slacks}; end; symbolic procedure check_minus_b(A,b); % % The algorithm requires the r.h.s. (i.e. the b matrix) to contain % only positive entries. % begin for i:=1:row_dim(b) do << if get_num_part( reval getmat(b,i,1) ) < 0 then << b := mult_rows(b,i,-1); A := mult_rows(A,i,-1); >>; >>; return {A,b}; end; symbolic procedure create_phase1_A1_and_obj_and_ib(A); begin scalar phase1_obj,A1,ib; integer column_dim_A1,column_dim_A,row_dim_A1,i; column_dim_A := fast_column_dim(A); % Add artificials to A. A1 := fast_matrix_augment({A,fast_make_identity(fast_row_dim(A))}); column_dim_A1 := fast_column_dim(A1); row_dim_A1 := fast_row_dim(A1); phase1_obj := mkmatrix(1,fast_column_dim(A1)); for i:=column_dim_A+1:fast_column_dim(A1) do fast_setmat(phase1_obj,1,i,1); ib := for i:=column_dim_A+1:fast_column_dim(A1) collect i; return {A1,phase1_obj,ib}; end; symbolic procedure simplex_calculation(obj_mat,A,b,ib,Binv,xb); % % Applies the revised simplex algorithm. See above for details. % begin scalar rs1,sb,rs2,rs3,u,continue,obj_value; integer k,iter,ell; obj_value := compute_objective(get_cb(obj_mat,ib),xb); while continue neq 'optimal do << rs1 := rstep1(A,obj_mat,Binv,ib); sb := car rs1; k := cadr rs1; u := caddr rs1; continue := cadddr rs1; if continue neq 'optimal then << rs2 := rstep2(xb,sb); ell := cadr rs2; rs3 := rstep3(xb,obj_mat,b,Binv,A,ib,k,ell); iter := iter + 1; Binv := car rs3; obj_value := cadr rs3; xb := caddr rs3; >>; >>; return {obj_value,xb,ib,Binv}; end; symbolic procedure get_preliminary_variable_list(equation_list); % % Gets all variables before slack variables are added. % begin scalar variable_list; for each equation in equation_list do variable_list := union(variable_list,get_coeffs(cadr equation)); return variable_list; end; symbolic procedure add_slacks_to_equations(equation_list); % % Takes list of equations (=, <=, >=) and adds required slack % variables. Also returns all the rhs integers in a column matrix, % and a list of the added slack variables. % begin scalar slack_list,rhs_mat,slack_variable,slack_variable_list; integer i,row; rhs_mat := mkmatrix(length equation_list,1); row := 1; for each equation in equation_list do << if not numberp reval caddr equation then << prin2 "***** Error in simplex(third argument): "; rederr "right hand side of each inequality must be a number"; >> else fast_setmat(rhs_mat,row,1,caddr equation); row := row+1; % % Put in slack/surplus variables where required. % if car equation = 'geq then << i := i+1; slack_variable := mkid('sl_var,i); equation := {'plus,{'minus,mkid('sl_var,i)},cadr equation}; slack_variable_list := append(slack_variable_list, {slack_variable}); >> else if car equation = 'leq then << i := i+1; slack_variable := mkid('sl_var,i); equation := {'plus,mkid('sl_var,i),cadr equation}; slack_variable_list := append(slack_variable_list, {slack_variable}); >> else if car equation = 'equal then equation := cadr equation else << prin2 "***** Error in simplex(third argument):"; rederr "inequalities must contain either >=, <=, or =."; >>; slack_list := append(slack_list,{equation}); >>; return {slack_list,rhs_mat,slack_variable_list}; end; flag('(add_slacks_to_list),'opfn); symbolic procedure simp_get_A(slack_equations,variable_list); % % Extracts the A matrix from the slack equations. % begin scalar A,slack_elt,var_elt; integer row,col,length_slack_equations,length_variable_list; length_slack_equations := length slack_equations; length_variable_list := length variable_list; A := mkmatrix(length slack_equations,length variable_list); for row := 1:length_slack_equations do << for col := 1:length_variable_list do << slack_elt := nth(slack_equations,row); var_elt := nth(variable_list,col); fast_setmat(A,row,col,smplx_prepsq( algebraic coeffn(slack_elt,var_elt,1))); >>; >>; return A; end; symbolic procedure get_X_and_obj_mat(objective,variable_list); % % Converts the variable list into a matrix and creates the objective % matrix. This is the matrix s.t. obj_mat * X = objective function. % begin scalar X,obj_mat; integer i,length_variable_list,tmp; length_variable_list := length variable_list; X := mkmatrix(length_variable_list,1); obj_mat := mkmatrix(1,length_variable_list); for i := 1:length variable_list do << fast_setmat(X,i,1,nth(variable_list,i)); tmp := nth(variable_list,i); fast_setmat(obj_mat,1,i,algebraic coeffn(objective,tmp,1)); >>; return {X,obj_mat}; end; symbolic procedure get_cb(obj_mat,ib); % % Gets hold of the columns of the objective matrix that are pointed % at in ib. % fast_augment_columns(obj_mat,ib); symbolic procedure compute_objective(cb,xb); % % Objective value := cb * xb (both are matrices) % fast_getmat(reval {'times,cb,xb},1,1); symbolic procedure rstep1(A,obj_mat,Binv,ib); % % Implements step 1 of the revised simplex algorithm. % ie: Computation of search direction sb. % % See above for details. (comments in simplex). % begin scalar u,sb,sum,i_in_ib; integer i,j,m,n,k,vkmin; m := fast_row_dim(A); n := fast_column_dim(A); u := mkmatrix(m,1); sb := mkmatrix(m,1); % Compute u. u := reval {'times,{'minus,algebraic (tp(Binv))}, algebraic tp(symbolic get_cb(obj_mat,ib))}; k := 0; vkmin := 10^10; i := 1; for i:=1:n do << i_in_ib := nil; % Check if i is in ib. for j:=1:m do << if i = nth(ib,j) then i_in_ib := t; >>; if not i_in_ib then << sum := specrd!:plus(smplx_prepsq(fast_getmat(obj_mat,1,i)), two_column_scalar_product(fast_augment_columns(A,i),u)); % if i is not in ib. %sum := fast_getmat(obj_mat,1,i); %for p:=1:m do %<< %sum := reval % {'plus,sum,{'times,fast_getmat(A,p,i),fast_getmat(u,p,1)}}; %>>; if get_num_part(sum) geq get_num_part(vkmin) then <<>> else << vkmin := sum; k := i; >>; >>; >>; % Do we need a tolerance here? if get_num_part(vkmin) < 0 then << % Form sb. for i:=1:m do << sum := 0; for j:=1:m do sum := specrd!:plus(sum, specrd!:times(fast_getmat(Binv,i,j),fast_getmat(A,j,k))); fast_setmat(sb,i,1,sum); >>; return {sb,k,u,nil}; >> else return {sb,k,u,'optimal}; end; symbolic procedure rstep2(xb,sb); % % step 2: Computation of maximum feasible step size Ob. % % see above for details. (comments in simplex). % begin scalar ratio; integer ell,sigb; sigb := 1*10^30; for i:=1:fast_row_dim(sb) do << if get_num_part(my_reval fast_getmat(sb,i,1)) leq 0 then <<>> else << ratio := specrd!:quotient(smplx_prepsq(fast_getmat(xb,i,1)), smplx_prepsq(fast_getmat(sb,i,1))); if get_num_part(ratio) geq get_num_part(sigb) then <<>> else << sigb := ratio; ell := i; >>; >>; >>; if ell= 0 then rederr "Error in simplex: The problem is unbounded."; return {sigb,ell}; end; symbolic procedure rstep3(xb,obj_mat,b,Binv,A,ib,k,ell); % % step3: Update. % % see above for details. (comments in simplex). % begin scalar work,Binv; work := fast_augment_columns(A,k); Binv := phiprm(Binv,work,ell); xb := reval{'times,Binv,b}; nth(ib,ell) := k; obj_mat := compute_objective(get_cb(obj_mat,ib),xb); return {Binv,obj_mat,xb}; end; symbolic procedure phiprm(Binv,D,ell); % % Replaces B^(-1) with [phi((B^(-1)',A(k),l)]'. % begin scalar sum,temp; integer m,j; m := fast_column_dim(Binv); sum := scalar_product(fast_stack_rows(Binv,ell),D); % if get_num_part(sum) = 0 then % rederr %{"Error in simplex: new matrix would be singular. Inner product = 0."}; if not zerop get_num_part(sum) then sum := specrd!:quotient(1,sum); Binv := fast_mult_rows(Binv,ell,sum); for j:=1:m do << if j = ell then <<>> else << temp := fast_getmat(reval{'times,fast_stack_rows(Binv,j),D}, 1,1); Binv := fast_add_rows(Binv,ell,j,{'minus,temp}); >>; >>; return Binv; end; symbolic procedure make_answer_list(xb,ib,no_coeffs,X,no_variables); % % Creates a list of the values of the variables at the optimal % solution. % begin scalar x_mat,ans_list; integer i; x_mat := mkmatrix(1,no_coeffs); i := 1; for each elt in ib do << if fast_getmat(xb,i,1) neq 0 then fast_setmat(x_mat,1,elt,fast_getmat(xb,i,1)); i := i+1; >>; ans_list := for i:=1:no_variables collect {'equal,my_reval fast_getmat(X,i,1), get_num_part(my_reval fast_getmat(x_mat,1,i))}; return ans_list; end; % Speed functions symbolic procedure fast_add_rows(in_mat,r1,r2,mult1); % % Replaces row2 (r2) by mult1*r1 + r2 without messing around. % begin scalar new_mat,fast_getmatel; integer i,coldim; coldim := fast_column_dim(in_mat); new_mat := copy_mat(in_mat); if (my_reval mult1) = 0 then return new_mat; for i:=1:coldim do << if not((fast_getmatel :=my_reval fast_getmat(new_mat,r1,i)) = 0) then fast_setmat(new_mat,r2,i,specrd!:plus(specrd!:times( smplx_prepsq(mult1),smplx_prepsq(fast_getmatel)),smplx_prepsq( fast_getmat(in_mat,r2,i)))); >>; return new_mat; end; symbolic procedure fast_augment_columns(in_mat,col_list); % % Quickly augments columns of in_mat specified in col_list. % if atom col_list then 'mat.for i:=1:fast_row_dim(in_mat) collect {fast_getmat(in_mat,i,col_list)} else 'mat.for each row in cdr in_mat collect for each elt in col_list collect nth(row,elt); symbolic procedure fast_matrix_augment(mat_list); % % As in linear_algebra package but doesn't produce !*sq output. % begin scalar ll,new_list; if length mat_list = 1 then return mat_list else << new_list := {}; for i:=1:fast_row_dim(car mat_list) do << ll := {}; for each mat1 in mat_list do ll := append(ll,nth(cdr mat1,i)); new_list := append(new_list,{ll}); >>; return 'mat.new_list; >>; end; symbolic procedure fast_setmat(matri,i,j,val); % % Set matrix element (i,j) to val. % fast_my_letmtr(list(matri,i,j),val,matri); symbolic procedure fast_unchecked_getmatelem u; nth(nth(cdr car u,cadr u),caddr u); symbolic procedure fast_mult_rows(in_mat,row_list,mult1); % % In simplex row_list is always an integer. % begin scalar new_list,new_row; integer row_no; row_no := 1; for each row in cdr in_mat do << if row_no neq row_list then new_list := append(new_list,{row}) else << new_row := for each elt in row collect my_reval{'times,mult1,elt}; new_list := append(new_list,{new_row}); >>; row_no := row_no+1; >>; return 'mat.new_list; end; symbolic procedure fast_make_identity(sq_size); % % Creates identity matrix. % 'mat. (for i:=1:sq_size collect for j:=1:sq_size collect if i=j then 1 else 0); symbolic procedure two_column_scalar_product(col1,col2); % % Calculates tp(col1)*col2. % % Uses sparsity efficiently. % begin scalar sum; sum := 0; for i:=1:length cdr col1 do << if car nth(cdr col1,i)=0 or car nth(cdr col2,i)=0 then <<>> else sum := specrd!:plus(sum,specrd!:times(smplx_prepsq( car nth(cdr col1,i)),smplx_prepsq( car nth(cdr col2,i)))); >>; return sum; end; symbolic procedure scalar_product(row,col); % % Calculates row*col. % % Uses sparsity efficiently. % begin scalar sum; sum := 0; for i:=1:length cadr row do << if nth(cadr row,i)=0 or car nth(cdr col,i)=0 then <<>> else sum := specrd!:plus(sum, specrd!:times(smplx_prepsq(nth(cadr row,i)), smplx_prepsq(car nth(cdr col,i)))); >>; return sum; end; endmodule; % simplex. end;