home *** CD-ROM | disk | FTP | other *** search
/ Collection of Education / collectionofeducationcarat1997.iso / SCIENCE / PRAXIS.ZIP / PRAXIS.INC < prev    next >
Text File  |  1987-07-15  |  24KB  |  565 lines
  1. (****************************************************************************)
  2. (*              P R O C E D U R E      P R A X I S                          *)
  3. (*                                                                          *)
  4. (*   P R A X I S  is for the minimization of a function in several          *)
  5. (*   variables. The algorithm used is a modification of a conjugate         *)
  6. (*   gradient method developed by Powell. Changes are due to Brent,         *)
  7. (*   who gives an ALGOL W program.                                          *)
  8. (*   Users who are interested in more of the details should read            *)
  9. (*          - Powell, M. J. D., 1962. An efficient method for finding       *)
  10. (*            the minimum of a function of several variables without        *)
  11. (*            calculating derivatives, Computer Journal 7, 155-162.         *)
  12. (*          - Brent, R. P., 1973. Algorithms for minimization without       *)
  13. (*            derivatives. Prentice Hall, Englewood Cliffs.                 *)
  14. (*   If you have any further comments, questions, etc. please feel free     *)
  15. (*   and contact me.                                                        *)
  16. (*                                                                          *)
  17. (*                           Karl Gegenfurtner     02/17/86                 *)
  18. (*                                                 TURBO - P Version 1.0    *)
  19. (****************************************************************************)
  20. (*  The use of PRAXIS is fairly simple. There are only two parameters:      *)
  21. (*      - X is of type vector and contains on input an initial guess        *)
  22. (*        of the solution. on output X holds the final solution of the      *)
  23. (*        system of equations.                                              *)
  24. (*      - N is of type integer and gives the number of unknown parameters   *)
  25. (*  The result of PRAXIS is the least calculated value of the function F    *)
  26. (*  F is one of the global parameters used by PRAXIS:                       *)
  27. (*    - F(X,N) is declared forward and is the function to be minimized      *)
  28. (*  All other globals are optional, i.e. you can change them or else        *)
  29. (*  PRAXIS will use some default values, which are adequate for most        *)
  30. (*  problems under consideration.                                           *)
  31. (*    - Prin controls the printout from PRAXIS                              *)
  32. (*           0:  no printout at all                                         *)
  33. (*           1:  only initial and final values                              *)
  34. (*           2:  detailed map of the minimization process                   *)
  35. (*           3:  also prints eigenvalues and eigenvectors of the            *)
  36. (*               direction matices in use (for insiders only).              *)
  37. (*    - T is the tolerance for the precision of the solution                *)
  38. (*      PRAXIS returns if the criterion                                     *)
  39. (*             2 * ||x(k)-x(k-1)|| <= Sqrt(MachEps) * ||x(k)|| + T          *)
  40. (*             is fulfilled more than KTM times, where                      *)
  41. (*    - KTM is another global parameter. It's default value is 1 and        *)
  42. (*      a value of 4 leads to a very(!) cautious stopping criterion         *)
  43. (*    - MachEps is the relative machine precision and is                    *)
  44. (*      1.0 E-15 using an 8087 and TURBO-87 and                             *)
  45. (*      1.0 E-06 without 8087.                                              *)
  46. (*    - H is a steplength parameter and shoul be set to the expected        *)
  47. (*      distance to the solution. An exceptional large or small value       *)
  48. (*      of H leads to slower convergence on the first few iterations.       *)
  49. (*    - Scbd is a scaling parameter and should be set to about 10.          *)
  50. (*      The default is 1 and with that value no scaling is done at all      *)
  51. (*      It's only necessary when the parameters are scaled very different   *)
  52. (*    - IllC is a Boolean variable and should be set to TRUE if the         *)
  53. (*      problem is known to be illconditioned.                              *)
  54. (****************************************************************************)
  55.  
  56.  
  57. CONST MaxPar = 20;
  58.  
  59. TYPE  Vector = ARRAY[1..MaxPar] OF REAL;
  60.       Matrix = ARRAY[1..MaxPar] OF Vector;
  61.       PrString = STRING[80];
  62.  
  63. CONST Macheps: REAL = 1.0E-08;  (* set to 1.0E-13 for TURBO-87 *)
  64.       H:       REAL = 1.0E+00;
  65.       T:       REAL = 1.0E-05;
  66.       Prin:    INTEGER = 2;
  67.       IllC:    BOOLEAN = FALSE;
  68.       Scbd:    REAL = 1;
  69.       Ktm:     INTEGER = 2;
  70.  
  71.  
  72. FUNCTION F(VAR x:Vector;n:INTEGER):REAL; FORWARD;
  73.  
  74. FUNCTION Power(a,b: REAL):REAL;
  75. BEGIN    Power := Exp(b*Ln(a));
  76. END;     (* Power *)
  77.  
  78.  
  79. PROCEDURE MinFit(N:INTEGER;Eps,Tol:REAL;VAR ab:Matrix;VAR q:Vector);
  80. LABEL     TestFsplitting,
  81.           TestFconvergence,
  82.           Convergence,
  83.           Cancellation;
  84. VAR       l, kt,
  85.           l2,
  86.           i, j, k: INTEGER;
  87.           c, f, g,
  88.           h, s, x,
  89.           y, z:    REAL;
  90.           e:       Vector;
  91. BEGIN     (* Householders reduction to bidiagonal form *)
  92.           x:= 0; g:= 0;
  93.           FOR i:= 1 TO N DO
  94.           BEGIN e[i]:= g; s:= 0; l:= i+1;
  95.                 FOR j:= i TO N DO
  96.                     s:= s + Sqr(ab[j,i]);
  97.                 IF s<Tol THEN g:= 0 ELSE
  98.                 BEGIN f:= ab[i,i];
  99.                       IF f<0
  100.                          THEN g:= Sqrt(s)
  101.                          ELSE g:= - Sqrt(s);
  102.                       h:= f*g-s;  ab[i,i]:= f - g;
  103.                       FOR j:= l TO N DO
  104.                       BEGIN f:= 0;
  105.                             FOR k:= i TO N DO
  106.                                 f:= f + ab[k,i]*ab[k,j];
  107.                             f:= f/h;
  108.                             FOR k:= i TO N DO
  109.                                 ab[k,j]:= ab[k,j] + f*ab[k,i];
  110.                       END; (* j *)
  111.                 END; (* IF *)
  112.                 q[i]:= g; s:= 0;
  113.                 IF i<=N
  114.                    THEN FOR j:= l TO N DO
  115.                             s:= s + Sqr(ab[i,j]);
  116.                 IF s<Tol THEN g:= 0 ELSE
  117.                 BEGIN f:= ab[i,i+1];
  118.                       IF f<0
  119.                          THEN g:= Sqrt(s)
  120.                          ELSE g:= - Sqrt(s);
  121.                       h:= f*g-s;  ab[i,i+1]:= f-g;
  122.                       FOR j:= l TO N DO e[j]:= ab[i,j]/h;
  123.                       FOR j:= l TO N DO
  124.                       BEGIN s:= 0;
  125.                             FOR k:= l TO N DO s:= s + ab[j,k]*ab[i,k];
  126.                             FOR k:= l TO N DO ab[j,k]:= ab[j,k] + s*e[k];
  127.                       END; (* J *)
  128.                 END; (* IF *)
  129.                 y:= Abs(q[i])+Abs(e[i]);
  130.                 IF y > x THEN x:= y;
  131.           END; (* I *)
  132.           (* accumulation of right hand transformations *)
  133.           FOR i:= N DOWNTO 1 DO
  134.           BEGIN IF g<>0.0 THEN
  135.                 BEGIN h:= ab[i,i+1]*g;
  136.                       FOR j:= l TO N DO ab[j,i]:= ab[i,j]/h;
  137.                       FOR j:= l TO N DO
  138.                       BEGIN s:= 0;
  139.                             FOR k:= l TO N DO s:= s + ab[i,k]*ab[k,j];
  140.                             FOR k:= l TO N DO ab[k,j]:= ab[k,j] + s*ab[k,i];
  141.                       END; (* J *)
  142.                 END; (* IF *)
  143.                 FOR j:= l TO N DO
  144.                 BEGIN ab[j,i]:= 0;
  145.                       ab[i,j]:= 0;
  146.                 END;
  147.                 ab[i,i]:= 1; g:= e[i]; l:= i;
  148.           END; (* I *)
  149.           (* diagonalization to bidiagonal form *)
  150.           eps:= eps*x;
  151.           FOR k:= N DOWNTO 1 DO
  152.           BEGIN kt:= 0;
  153. TestFsplitting: kt:= kt + 1;
  154.                 IF kt > 30 THEN
  155.                 BEGIN e[k]:= 0;
  156.                       WriteLn('+++ QR failed');
  157.                 END;
  158.                 FOR l2:= k DOWNTO 1 DO
  159.                 BEGIN l:= l2;
  160.                       IF Abs(e[l])<=Eps THEN Goto TestFconvergence;
  161.                       IF Abs(q[l-1])<=Eps THEN Goto Cancellation;
  162.                 END; (* L2 *)
  163. C