reduce-historical/r35/lib/linineq.red

module linineq; % Linear inequalities and linear optimization.

% Author:       Herbert Melenk <melenk@sc.zib-berlin.dbp.de>

% Version 1     January 1990
% Version 1.1   February 1990
%               added parameter "record=t"
% Version 2     May 1991
%               added Branch-and-Bound for Integer Prgramming

%
% Solution of linear inequalities & equations with numerical
% coefficients.
%
%   Fourier(1826) /Motzkin(1936): George. B. Dantzig,
%                  Linear Programming and Extensions.
%   Branch and Bound: Beightler, Phillips, Wilde,
%                  Foundations of Optimization
%                  Prentice Hall, 1979
%

put('linineq,'psopfn,'linineqeval);

global '(!*trlinineq !*trlinineqint !*prlinineq);
switch trlinineq,prlinineq,trlinineqint;
fluid '(linineqinterval!* linineqrecord!*);

algebraic operator <=,>=;

symbolic procedure linineqeval(u);
  % Interface for algebraic mode.
   begin scalar prob,equa,requa,vars,oldorder,res,u1,x,y,p,e,msg;
         scalar direction,rec,linineqrecord!*,r,intvars;
      msg := !*prlinineq or !*trlinineq;
      u1 := getrlist reval car u; u := cdr u;
      if u then <<x:= reval car u; u := cdr u>>;
      if eqcar(x,'LIST) then vars := getrlist x else
                       u := x . u;
      while u do <<x := reval car u; u := cdr u;
                  if eqcar(x,'EQUAL) and 
                     ((cadr x='RECORD and (rec:=t)) or
                      (cadr x='INT and (intvars:=getrlist caddr x)))
                  then t else
                  <<print x; rederr "illegal parameter">>;
                >>;
      x := vars := for each u in vars collect
       <<u := reval u;
         if eqcar(u,'equal) then
           if  member(caddr u,'(min max)) then
            <<direction:=(cadr u. caddr u) . direction;
              u := cadr u>> else
                 rederr "illegal form in 2nd parameter";
         u>>;
      for each u in u1 collect
      <<u := reval u;
        if eqcar(u,'geq)then
           prob := (simp cadr u . simp caddr u) . prob else
        if eqcar(u,'leq)then
           prob := (simp caddr u . simp cadr u) . prob else
        if eqcar(u,'equal) then
           equa := (simp caddr u . simp cadr u) . equa else
           rederr "relational operator '>=','<='or'='missing" >>;
        % control the linearity
      for each p in append(equa,prob) do
      <<if not domainp denr car p or not domainp denr cdr p
          then rederr "unable to process nonlinear system";
        vars := linineqevaltest(numr car p,
                        linineqevaltest(numr cdr p,vars));
      >>;
      if msg then <<prin2 "variables:"; prin2t vars;>>;
      oldorder := setkorder vars;
      prob := for each p in prob collect
       (reorder numr car p./denr car p).
             (reorder numr cdr p./denr cdr p);
      equa:= for each p in equa collect
       (reorder numr car p./denr car p).
             (reorder numr cdr p./denr cdr p);
       % eliminate variables from equations
      while equa do
      <<e := car equa; equa := cdr equa;
        e := addsq(car e,negsq cdr e);
        if domainp numr e then
        <<if numr e then  % nonzero constant equated to 0
          rederr "equation part inconsistent">>
         else
        <<u := list((x := mvar numr e).
                  prepsq(y:=multsq(negf red numr e ./ 1,
                                   invsq(lc numr e ./ 1))));
                 if member(x,intvars) then 
                          % Dont eliminate integer variables;
                          % represent equation by double inequality instead.
                 <<x:=simp x; prob:=append({x.y,y.x},prob)>> 
                 else
                 <<
          prob := for each p in prob collect
                        subsq(car p,u).subsq(cdr p,u);
          equa := for each p in equa collect
                        subsq(car p,u).subsq(cdr p,u);
          requa := append(u,requa);
          if msg then
            <<prin2 "         ";prin2 x;
              prin2 " eliminated by equation";
              terpri()>>;
          vars := delete(x,vars);
                 >>
      >> >>;
      res := if intvars then linineqint(prob,vars,msg,direction,rec,intvars)
             else linineq1(prob,vars,msg,direction,rec);
          % backsubstitution in equations;
      for each e in requa do
      <<x := prepsq subsq(y:=simp cdr e,res);
        res := (car e . x) . res;
        if rec then
        <<x := prepsq y;
          linineqrecord!* := list(x,x) . linineqrecord!*>>;
      >>;
      setkorder oldorder;
      r := if rec then for each p in pair(res,linineqrecord!*) collect
               list('LIST,list('EQUAL,caar p,cdar p),cadr p,caddr p)
            else
           for each p in res collect list('EQUAL,car p,cdr p);
      return 'LIST . r;
   end;

symbolic procedure linineqevaltest(f,v);
   % Collect the variables in standard form f and control linearity.
     if domainp f then v else
     if not(ldeg f=1) then
           rederr "unable to process nonlinear system" else
     if member(mvar f,v) then linineqevaltest(red f,v) else
         linineqevaltest(red f,mvar f.v);

symbolic procedure linineq0(prob,vars,dir,rec);
  % Interface for symbolic mode.
  % Prob is a list (e1,e2,..) of algebraic expressions without
  % relational operators, which are interpreted as
  % set of inequalities ei >= 0. They are linear in the
  % variables vars.
  % Silent operation: result=nil if the system is inconsistent.
   begin scalar oldorder,res;
      linineqrecord!* := nil;
      oldorder := setkorder vars;
      prob := for each u in prob collect simp u . (nil ./ 1);
      res := linineq1(prob,vars,nil,dir,rec);
      setkorder oldorder;
      return res;
   end;

symbolic procedure linineqint(prob,vars,msg,dir,rec,intvars);
  begin scalar x,x0,y,y0,y1,z,w,problems,best,z,z0,zbest,zf,bestr;
       % test integer variables and adjust order;
    for each x in vars do
      if member(x,intvars) then<<w:=x.w;intvars:=delete(x,intvars)>>;
    if intvars then typerr('list.intvars,"int variables");
    intvars:=reversip w;
       % select primary optimization principle.
        if dir then<<z:=caar dir;zf:=if cdar dir='MAX then 1 else -1>>;
    problems:=list (nil.prob);
       % macro loop.
    while problems do
    <<z0:=caar problems; prob:=cdar problems; problems:=cdr problems;
      if msg or !*trlinineqint 
        then linineqprint2("=== next integer subproblem",prob);
      w:=if best and not evalgreaterp({'times,zf,z0},{'times,zf,zbest})
          then nil  % skip problem with suboptimal bound.
         else linineq1(prob,vars,msg,dir,rec);
      if !*trlinineqint then linineqprint3("=== subresult",w);
      if w and dir then
      <<% is better than best so far?
        z0:=cdr assoc(z,w);
        if best and evalgreaterp({'times,zf,zbest},{'times,zf,z0})
           then w:=nil;
      >>;
      if w then
      <<% test feasability;
        y:=list prob;
        for each x in intvars do
        <<x0:=cdr assoc(x,w);
          if not fixp x0 then  % branch and bound
          <<x:=simp x; y0:=simp{'ceiling,x0}; y1:=simp {'floor,x0};
            y:= for each q in y join {(x.y0).q, (y1.x).q};
            if msg or !*trlinineqint then
            <<writepri("branch and bound with",'first);
              writepri(mkquote{'list,{'geq,x:=prepsq x,prepsq y0},
                                     {'leq,x,prepsq y1}},'last);
           >>;
          >>;
        >>;
        if cdr y then
         problems:=append(problems,for each q in y collect z0.q)
        else
         <<zbest:=z0; best:=w; bestr:=linineqrecord!*;
           if !*trlinineqint then prin2t "===>  is feasable";>>;
       >>;  % if w
               % without target dont need additional result.
           if best and null dir then problems:=nil; 
     >>;  % while problems
   linineqrecord!*:=bestr;
   return best;
 end;

symbolic procedure linineq1(prob,vars,msg,dir,rec);
  % Algebraic evaluation of a set of inequalities:
  % prob is a list of pairs of standard quotients,
  % (( p1 . q1)(p2 . q2) .. (pn . qn))
  % which are interpreted as inequalities:
  %     pi >= qi ;
  % vars is the list of (linear) variables.
  % dir  the direction of final optimization
  % rec  switch; if t, the record of inequatlities is produced
  % Result is NIL if the system has no solution; otherwise
  % the solution has the form of an association list
  %  ((v1 . val1)(v2 . val2) ... (vn . valn)),
  % where vi are the variables and vali are values in algebraic
  % form. NIL if the system has no solution.
  %
   begin scalar v,vq,lh,rh,x,y,z,prob1,prob2,prob3,prob4,nprob,sw,sol;
      v := car vars; vars := cdr vars;
      vq := mksq(v,1);
      if !*trlinineq then
       linineqprint2(list("next variable:",v,"; initial system:"),prob);
      prob := linineqnormalize prob;
      for each p in prob do
       <<lh := car p; rh := cdr p;
          % if v appears on the lhs, isolate it
         if not domainp numr lh and mvar numr lh = v then
         <<x :=  invsq(lc numr lh ./ 1);
           sw := (numr x < 0);
           lh := multsq(lh,x); rh := multsq(rh,x);
           rh := addsq(rh,negf red numr lh ./ denr lh);
           if not sw then prob1 := (vq . rh) . prob1 else
                          prob2 := (rh . vq) . prob2;
         >>else if domainp numr rh and domainp numr lh then
                prob4 := (lh . rh) . prob4 else
                prob3 := (lh . rh) . prob3;
       >>;
      if null prob1 and null prob2 and vars then
      << sol := linineq1(prob,vars,msg,dir,rec);
         if rec then linineqrecord!* :=
             append(linineqrecord!*,'(((MINUS INF),'INF)));
         return if sol then (v . 0) . sol else nil>>;
      if !*trlinineq then
      <<linineqprint2("class 1:",prob1);
        linineqprint2("class 2:",prob2);
        linineqprint2("class 3:",prob3);
        linineqprint2("class 4:",prob4);
      >>;
      if rec then
      << x := for each u in prob1 collect prepsq cdr u;
         y := for each u in prob2 collect prepsq car u;
         x := if null x then '(MINUS INF) else
              if null cdr x then car x else 'MAX . x;
         y := if null y then ' INF else
              if null cdr y then car y else 'MIN . y;
         linineqrecord!* := append(linineqrecord!*, list list(x,y))
      >>;
      if not linineq2(prob4,msg) then return nil;
      nprob := append(prob3,
         for each x in prob1 join
           for each y in prob2 collect
             car y . cdr x);
      if vars then
       << if null (sol := linineq1(nprob,vars,msg,dir,rec)) then
                         return nil>>
        else if not linineq2(nprob,msg) then return nil;
         % lower bound:
      x := if null prob1 then nil else
         linineqevalmax for each p in prob1 collect
                subsq(cdr p,sol);
         % upper bound:
      y := if null prob2 then nil else
         linineqevalmin for each p in prob2 collect
                subsq(car p,sol);
      if (z:=assoc(v,dir)) then z:= cdr z;
      if msg then
      <<writepri("         ",'first);
        writepri(mkquote if x then prepsq x else '(minus inf),nil);
        writepri(" <= ",nil);
        writepri(mkquote v,nil);
        writepri(" <= ",nil);
        writepri(mkquote if y then prepsq y else 'inf,nil);
        writepri(";   ",nil)>>;
     linineqinterval!* := x . y;
     if z='min and null x or z='max and null y then
      <<if msg then writepri( " max/min cannot be resolved",'last);
        return nil>>;
      if not x=y then
        if z='min then y:=nil else if z='max then x:=nil;
      if msg then
      << writepri(
        if null x and null y then " completely free: " else
        if null y then " minimum: " else
        if null x then " maximum: " else
        if x=y then " zero length interval: " else
                " middle: ",nil);
       >>;
      if null x and null y then x := 0 else % completely free
      if null x then x := prepsq y else
      if null y then x := prepsq x else
      if sqlessp(y,x) then
        <<prin2 "system inconsistent:";
          prin2 prepsq x; prin2 " not <= "; prin2t prepsq y;
          return nil>> else
        x := list('quotient,list('plus,prepsq x,prepsq y),2);
      x := aeval x;
      if msg then
        writepri(mkquote list('EQUAL,v,x),'last);
      return (v . x) . sol;
   end;

symbolic procedure linineq2(prob,msg);
   % All variables are elimitated. Control, if the
   % remaining numerical inequalities are consistent.
   begin scalar rh,lh;
 loop: if null prob then return t;
      lh := caar prob; rh := cdar prob;
      if not domainp numr rh or not domainp numr lh then
        <<prin2t " non numeric:"; print rh; print lh;
          rederr "linineq";>>;
      if sqlessp(lh,rh) then
        <<if msg then <<writepri("system inconsistent: ",'first);
                        writepri(mkquote prepsq lh,nil);
                        writepri(" not >= ",nil);
                        writepri(mkquote prepsq rh,'last);>>;
          return nil>>;
      prob := cdr prob;
      goto loop;
    end;

symbolic procedure linineqnormalize prob;
    % Normalize system: reform all inequalities such that they have
    % the canonical form %     polynomial >= constant
    %      (canonical restriction: absolute term of lhs=0,
    %                              denominator of lhs = 1).
    % and remove those, which have same lhs, but smaller rhs
    % (the latter are superfluous).
   begin scalar r,lh,rh,d,ab,x;
     for each p in prob do
     <<lh := car p; rh := cdr p;
         % arithmetic normalizaion
       lh := addsq(lh,negsq rh);
       d := denr lh;
       lh := numr lh;
       ab := lh; x := if domainp lh then 1 else lc ab;
       while not domainp ab do <<x := gcdf(x,lc ab);ab := red ab>>;
       ab := negf ab;
       lh := multsq(addf(lh,ab)./1,1 ./ x);
       rh := multsq(ab ./ 1, 1 ./ x);
         % removal of redundant elements
       x := assoc(lh,r);
       if null x then r:=(lh.rh) . r else
         if sqlessp(cdr x,rh) then rplacd(x,rh);
     >>;
     if !*trlinineq then
         linineqprint2("normalized and reduced:",r);
     return r;
   end;

symbolic procedure linineqevalmin u;
   % Compute the minimum among the list u with sq's.
     linineqevalmin1(car u,cdr u);

symbolic procedure linineqevalmin1(q,u);
   if null u then q else
  (linineqevalmin1( if x and !:minusp x then q else car u, cdr u)
       ) where x=numr addsq(q,negsq car u);

symbolic procedure linineqevalmax u;
   % compute the maximum among the list u with sq's
     linineqevalmax1(car u,cdr u);

symbolic procedure linineqevalmax1(q,u);
   if null u then q else
  (linineqevalmax1(
     if x and !:minusp x then car u else q, cdr u)
        ) where x=numr addsq(q,negsq car u);

symbolic procedure sqlessp(q1,q2);
   (x and !:minusp x) where x=numr addsq(q1,negsq q2);


symbolic procedure linineqprint1(text,lh,rh);
         <<writepri(text,'first);
           writepri(mkquote prepsq lh,nil);
           writepri(" >= ",nil);
           writepri(mkquote prepsq rh,'last)>>;


symbolic procedure linineqprint2(text,prob);
    <<prin2t "--------------------------------";
          if atom text then text:={text};
      for each u in text do prin2 u; terpri();
          writepri(mkquote('list .
         for each p in prob collect
                    {'geq,prepsq car p,prepsq cdr p}),'last)
    >>;

symbolic procedure linineqprint3(text,res);
    <<writepri(text,'first);
          writepri(mkquote('list . for each p in res collect
                     {'equal,car p,cdr p}), 'last);
        >>;

endmodule;

end;