reduce-historical/r37/packages/numeric/steepstd.red

module steepstd;  % Optimization and root finding with method of
                  % steepest descent.

% Author: H. Melenk, ZIB, Berlin.

% Copyright (c): ZIB Berlin 1991, all rights reserved.

fluid '(!*noequiv accuracy!*);
global '(iterations!* !*trnumeric);

symbolic procedure rdmineval u;
  begin scalar e,vars,x,y,z,oldmode,p,!*noequiv;
    oldmode:=steepdecedmode(nil,nil);
    u := for each x in u collect reval x;
    u := accuracycontrol(u,6,40);
    e :=car u; u :=cdr u;
    if eqcar(car u,'list) then
      u := for each x in cdar u collect reval x;
    for each x in u do
     <<if eqcar(x,'equal) then
       <<z:=cadr x; y:=caddr x>> else <<z:=x; y:=random 100>>;
       vars:=z . vars; p := y . p>>;
    vars := reversip vars; p := reversip p;
    x := steepdeceval1(e,vars,p,'num_min);
    steepdecedmode(t,oldmode);
    return list('list, car x, 'list.
       for each p in pair(vars,cdr x) collect
          list('equal,car p,cdr p));
  end;

put('num_min,'psopfn,'rdmineval);

symbolic procedure rdsolvestdeval u;
   % solving a system of equations via minimizing the
   % sum of sqares.
  begin scalar e,vars,x,y,z,oldmode,p,q,!*noequiv;
    oldmode:=steepdecedmode(nil,nil);
    u := for each x in u collect reval x;
    e :=car u; u :=cdr u;
      % construct the function to minimize.
    e:=if eqcar(e,'list) then cdr e else list e;
    q := 'plus . for each f in e collect
       list('expt,if eqexpr f then !*eqn2a f else f,2);
    q := prepsq simp q;
      % starting point & variables.
    if eqcar(car u,'list) then
      u := for each x in cdar u collect reval x;
    for each x in u do
     <<if eqcar(x,'equal) then
       <<z:=cadr x; y:=caddr x>> else <<z:=x; y:=random 100>>;
       vars:=z . vars; p := y . p>>;
    vars := reversip vars; p := reversip p;
    x := steepdeceval1(q,vars,p,'root);
    steepdecedmode(t,oldmode);
    if null x then rederr "no solution found";
    return 'list.
       for each p in pair(vars,cdr x) collect
          list('equal,car p,cdr p);
  end;

symbolic procedure steepdecedmode(m,oldmode);
  if not m then % initial call
   << if !*complex then rederr
      "steepest descent method not applicable under complex";
    if not (dmode!*='!:rd!:)then
       <<oldmode:=t; setdmode('rounded,t)>>;
    {oldmode,precision 0,!*roundbf}>>
  else % reset to old dmode.
    <<if car oldmode then setdmode('rounded,nil);
      precision cadr oldmode;
      !*roundbf:=caddr oldmode>>;

symbolic procedure steepdeceval1(f,vars,x,mode);
  begin scalar e,g,r,acc;
        integer n;
     n := length vars;
    % establish the target function and the gradient.
     e := prepsq simp f;
     if !*trnumeric then lprim "computing symbolic gradient";
     g := for each v in vars collect prepsq simp list('df,f,v);
     !*noequiv:=t;
     if !*trnumeric then
      lprim "starting Fletcher Reeves steepest descent iteration";
     acc := !:!:quotient(1,expt(10,accuracy!*));
     x := for each u in x collect force!-to!-dm numr simp u;
     r:= steepdec2(e,g,vars,acc,x,mode);
     r:= for each x in r collect prepf x;
     return r;
 end;

symbolic procedure steepdec2(f,grad,vars,acc,x,mode);
  % Algorithm for finding the minimum of function f
  % with technique of steepest descent, version Fletcher-Reeves.
  % f:     function to minimize (standard quotient);
  % grad:  symbolic gradient (list of standard quotients);
  % vars:  variables (list of id's);
  % acc:   target precision (e.g. 0.000001)
  % x:     starting point (list of numerical standard forms).
  % mode:  minimum / roots type of operation
 dm!:
  begin scalar e0,e00,e1,e2,a,a1,a1a1,a2,a2a2,x1,x2,g,h,dx,
      delta,limit,gold,gnorm,goldnorm,multi;
    integer count,k,n;
    n:=length grad; multi:=n>1; n:=10*n;
    e00 := e0 := evaluate(f,vars,x);
    a1 := 1;
  init:
    k:=0;
      % evaluate gradient (negative).
    g := for each v in grad collect -evaluate(v,vars,x);
    gnorm := normlist g; h:=g;
  loop:
    count := add1 count; k:=add1 k;
    if count>iterations!* then
    <<lprim "requested accuracy not reached within iteration limit";
      % mode := nil;
      goto ready>>;

      % quadratic interpolation in direction of h in order
      % to find the minimum value of f in that direction.
    a2 := nil;
   l1:
    x1 := list!+list(x, scal!*list(a1,h));
    e1 := evaluate(f,vars,x1);
    if e1>e0 then <<a2:=a1; x2:=x1; e2:=e1;
            a1 := a1/2; goto l1>>;
    if a2 then goto alph;
   l2:
    a2 := a1+a1;
    x2 := list!+list(x, scal!*list(a2,h));
    e2 := evaluate(f,vars,x2);
    if e2<e1 then <<a1:=a2; e1:=e2; goto l2>>;
   alph:   %compute lowest point of parabel
    if abs(e1-e2)<acc then goto ready; % constant?
    a1a1:=a1*a1; a2a2:=a2*a2;
    a:= (a1a1-a2a2)*e0 + a2a2*e1 - a1a1*e2 ;
    a:= a/((a1-a2)*e0 + a2*e1-a1*e2);
    a:= a/2;

     % new point
    dx:=scal!*list(a,h); x:=list!+list(x, dx);
    e0 := evaluate(f,vars,x);
    if e0>e1 then % a1 was better than interpolation.
    <<dx:=scal!*list(a1-a ,h); x:=list!+list(x,dx);
      e0 := e1; dx:=scal!*list(a1,h)>>;
    delta:= normlist dx;
    steepdecprintpoint(count,x,delta,e0,gnorm);
    update!-precision list(delta,e0,gnorm);

      % compute next direction;
    if k>n then goto init; % reinitialize time to time.
    gold := g; goldnorm:= gnorm;
    g := for each v in grad collect -evaluate(v,vars,x);
    gnorm := normlist g;
    if gnorm<limit then goto ready;
    h := if not multi then g else
        list!+list(g,scal!*list(gnorm/goldnorm,h));
    if mode='num_min and gnorm > acc or
       mode='root and e0 > acc then goto loop;

  ready:
    if mode='root and not(abs e0 < acc) then
     <<lprim "probably fallen into local minimum of f^2";
        return nil>>;
    return e0 . x;
  end;

symbolic procedure steepdecprintpoint(count,x,d,e0,ng);
 if !*trnumeric then
   begin integer n;
    writepri(count,'first);
    writepri(". residue=",nil);
    writepri(mkquote prepf e0,nil);
    writepri(", gradient length=",nil);
    writepri(mkquote prepf ng,'last);
    writepri(" at (",nil);
    for each y in x do
       <<if n>0 then writepri(" , ",nil);
         n:=n+1;
         writepri(mkquote prepf y,nil)>>;
    writepri(")",nil);
    writepri(", stepsize=",nil);
    writepri(mkquote prepf d,nil);
    writepri("",'last);
   end;

endmodule;

end;