home *** CD-ROM | disk | FTP | other *** search
/ Fish 'n' More 2 / fishmore-publicdomainlibraryvol.ii1991xetec.iso / fish / applications / xlispstat / src / src2.lzh / XLisp-Stat / optimize.c < prev    next >
C/C++ Source or Header  |  1990-10-02  |  8KB  |  314 lines
  1. /* optimimize - basic optimization routines                            */
  2. /* XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney                  */
  3. /* Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz    */
  4. /* You may give out copies of this software; for conditions see the    */
  5. /* file COPYING included with this distribution.                       */
  6.  
  7. /*#ifdef OPTIMIZE  what does this do ?? JKL */
  8.  
  9. #include "xlisp.h"
  10. #include "osdef.h"
  11. #ifdef ANSI
  12. #include "xlproto.h"
  13. #include "xlsproto.h"
  14. #else
  15. #include "xlfun.h"
  16. #include "xlsfun.h"
  17. #endif ANSI
  18. #include "xlvar.h"
  19. #include "xlsvar.h"
  20.  
  21. #ifdef ANSI
  22. double evfun1(LVAL,double);
  23. #else
  24. double evfun1();
  25. #endif ANSI
  26.  
  27. /************************************************************************/
  28. /**                                                                    **/
  29. /**                     One Dimensional Search                         **/
  30. /**                                                                    **/
  31. /************************************************************************/
  32.  
  33. #define GOLD   1.618034
  34. #define GLIMIT 100.0
  35. #ifndef TINY
  36. # define TINY   1.0e-20
  37. #endif
  38. #ifndef HUGE
  39. # define HUGE   1.0e20
  40. #endif
  41. #define SIGN(a,b)     ((b) > 0.0 ? fabs(a) : -fabs(a))
  42. #define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);
  43. #define BRAKMAX 50
  44. #define ZEPS 1.0e-10
  45.  
  46. static double evfun1(f, x)
  47.     LVAL f;
  48.     double x;
  49. {
  50.   LVAL tmp;
  51.   
  52.   xlsave1(tmp);
  53.   tmp = cvflonum((FLOTYPE) x);
  54.   tmp = xsfuncall1(f, tmp);
  55.   xlpop();
  56.   return(makedouble(tmp));
  57. }
  58.  
  59. LVAL xsbracket_search()
  60. {
  61.   LVAL f, result, temp;
  62.   double ax, bx, cx, fa, fb, fc, ulim, u, r, q, fu, dum;
  63.   int i, itmax, verbose;
  64.   
  65.   f = xlgetarg();
  66.   ax = makedouble(xlgetarg());
  67.   bx = makedouble(xlgetarg());
  68.   if (xlgetkeyarg(sk_max_iters, &temp)) itmax = getfixnum(temp);
  69.   else itmax = BRAKMAX;
  70.   if (xlgetkeyarg(k_verbose, &temp)) verbose = (temp != NIL);
  71.   else verbose = FALSE;
  72.   
  73.   fa = evfun1(f, ax);
  74.   fb = evfun1(f, bx);
  75.   if (fb > fa) {
  76.     SHFT(dum, ax, bx, dum)
  77.     SHFT(dum, fb, fa, dum)
  78.   }
  79.   cx = bx + GOLD * (bx - ax);
  80.   fc = evfun1(f, cx);
  81.   for (i = 0; i < itmax && fb > fc; i++) {
  82.     if (verbose) {
  83.       sprintf(buf, "a = %f, b = %f, f(a) = %f, f(b) = %f\n",
  84.                    (ax < cx) ? ax : cx, (ax < cx) ? cx : ax,
  85.                    (ax < cx) ? fa : fc, (ax < cx) ? fc : fa);
  86.       stdputstr(buf);
  87.     }
  88.     r = (bx - ax) * (fb - fc);
  89.     q = (bx - cx) * (fb - fa);
  90.     u = bx - ((bx - cx) * q - (bx - ax) * r)
  91.       / (2.0 * SIGN(Max(fabs(q - r), TINY), q - r));
  92.     ulim = bx + GLIMIT * (cx - bx);
  93.     if ((bx - u) * (u - cx) > 0.0) {
  94.       fu = evfun1(f, u);
  95.       if (fu < fc) {
  96.         ax = bx;
  97.         bx = u;
  98.         fa = fb;
  99.         fb = fu;
  100.         break;
  101.       }
  102.       else if (fu > fb) {
  103.         cx = u;
  104.         fc = fu;
  105.         break;
  106.       }
  107.       u = cx + GOLD * (cx - bx);
  108.       fu = evfun1(f, u);
  109.     }
  110.     else if ((cx - u) * (u - ulim) > 0.0) {
  111.       fu = evfun1(f, u);
  112.       if (fu < fc) {
  113.         SHFT(bx, cx, u, cx + GOLD * (cx - bx))
  114.         SHFT(fb, fc, fu, evfun1(f, u))
  115.       }
  116.     }
  117.     else if ((u - ulim) * (ulim - cx) >= 0.0) {
  118.       u = ulim;
  119.       fu = evfun1(f, u);
  120.     }
  121.     else {
  122.       u = cx + GOLD * (cx - bx);
  123.       fu = evfun1(f, u);
  124.     }
  125.     SHFT(ax, bx, cx, u)
  126.     SHFT(fa, fb, fc, fu)
  127.   }
  128.   
  129.   if (ax > cx) {
  130.       dum = ax; ax = cx; cx = dum;
  131.       dum = fa; fa = fc; fc = dum;
  132.   }
  133.   
  134.   xlsave1(result);
  135.   result = mklist(2, NIL);
  136.   rplaca(result, mklist(3, NIL));
  137.   rplaca(cdr(result), mklist(3, NIL));
  138.   temp = car(result);
  139.   rplaca(temp, cvflonum((FLOTYPE) ax)); temp = cdr(temp);
  140.   rplaca(temp, cvflonum((FLOTYPE) bx)); temp = cdr(temp);
  141.   rplaca(temp, cvflonum((FLOTYPE) cx));
  142.   temp = car(cdr(result));
  143.   rplaca(temp, cvflonum((FLOTYPE) fa)); temp = cdr(temp);
  144.   rplaca(temp, cvflonum((FLOTYPE) fb)); temp = cdr(temp);
  145.   rplaca(temp, cvflonum((FLOTYPE) fc));
  146.   xlpop();
  147.   
  148.   return(result);
  149. }
  150.  
  151. #define R 0.61803399
  152. #define C (1.0 - R)
  153. #define GTOL .000001
  154.  
  155. LVAL xsgolden_search()
  156. {
  157.   double ax, bx, cx, tol, f0, f1, f2, f3, x0, x1, x2, x3, xmin, fmin;
  158.   LVAL f, arg, result;
  159.   int verbose;
  160.   
  161.   f = xlgetarg();
  162.   ax = makedouble(xlgetarg());
  163.   cx = makedouble(xlgetarg());
  164.   
  165.   if (xlgetkeyarg(k_start, &arg)) bx = makedouble(arg);
  166.   else bx = ax + C * (cx - ax);
  167.   
  168.   if (xlgetkeyarg(sk_tolerance, &arg)) tol = makedouble(arg);
  169.   else tol = GTOL;
  170.   
  171.   if (xlgetkeyarg(k_verbose, &arg)) verbose = (arg != NIL);
  172.   else verbose = FALSE;
  173.  
  174.   if (tol < ZEPS) tol = ZEPS;
  175.   
  176.   x0 = ax;
  177.   x3 = cx;
  178.   if (fabs(cx - bx) > fabs(bx - ax)) {
  179.     x1 = bx;
  180.     x2 = bx + C * (cx - bx);
  181.   }
  182.   else {
  183.     x2 = bx;
  184.     x1 = bx - C * (bx - ax);
  185.   }
  186.   f1 = evfun1(f, x1);
  187.   f2 = evfun1(f, x2);
  188.   while (fabs(x3 - x0) > tol * (1.0 + fabs(x1) + fabs(x2))) {
  189.     if (verbose) {
  190.       sprintf(buf, "x = %f, f(x) = %f\n", (f1 < f2) ? x1 : x2, (f1 < f2) ? f1 : f2);
  191.       stdputstr(buf);
  192.     }
  193.     if (f2 < f1) {
  194.       SHFT(x0, x1, x2, R * x1 + C * x3)
  195.       SHFT(f0, f1, f2, evfun1(f, x2)) /* redundant f0 in macro JKL */
  196.     }
  197.     else {
  198.       SHFT(x3, x2, x1, R * x2 + C * x0)
  199.       SHFT(f3, f2, f1, evfun1(f, x1)) /* redundant f3 in macro JKL */
  200.     }
  201.   }
  202.   if (f1 < f2) { xmin = x1; fmin = f1; }
  203.   else { xmin = x2; fmin = f2; }
  204.   
  205.   xlsave1(result);
  206.   xlpop();
  207.   result = mklist(2, NIL);
  208.   rplaca(result, cvflonum((FLOTYPE) xmin));
  209.   rplaca(cdr(result), cvflonum((FLOTYPE) fmin));
  210.   return(result);
  211. }
  212.  
  213. #define PARITMAX 100
  214. #define PARTOL .00001
  215. #define CGOLD 0.3819660
  216.  
  217. LVAL xsparabolic_search()
  218. {
  219.   int iter, itmax, verbose;
  220.   double ax, bx, cx, tol, a, b, d, etemp, fu, fv, fw, fx, p, q, r;
  221.   double tol1, tol2, u, v, w, x, xm;
  222.   double e = 0.0;
  223.   LVAL f, arg, result;
  224.   
  225.   d = 0.0;
  226.   
  227.   f = xlgetarg();
  228.   ax = makedouble(xlgetarg());
  229.   cx = makedouble(xlgetarg());
  230.   
  231.   if (xlgetkeyarg(k_start, &arg)) bx = makedouble(arg);
  232.   else bx = ax + CGOLD * (cx - ax);
  233.   
  234.   if (xlgetkeyarg(sk_tolerance, &arg)) tol = makedouble(arg);
  235.   else tol = PARTOL;
  236.   
  237.   if (xlgetkeyarg(k_verbose, &arg)) verbose = (arg != NIL);
  238.   else verbose = FALSE;
  239.   
  240.   if (xlgetkeyarg(sk_max_iters, &arg)) {
  241.     if (! fixp(arg)) xlerror("not an integer", arg);
  242.     itmax = getfixnum(arg);
  243.   }
  244.   else itmax = PARITMAX;
  245.   
  246.   if (tol < TINY) xlfail("torerance is too small");
  247.   
  248.   a = ((ax < cx) ? ax : cx);
  249.   b = ((ax > cx) ? ax : cx);
  250.   x = w = v = bx;
  251.   fw = fv = fx = evfun1(f, x);
  252.   for (iter = 0; iter < itmax; iter++) {
  253.     xm = 0.5 * (a + b);
  254.     tol2 = 2.0 * (tol1 = tol * (1.0 + fabs(x)) + ZEPS);
  255.     if (fabs(x - xm) <= (tol2 - 0.5 * (b - a))) {
  256.       break;
  257.     }
  258.     if (fabs(e) > tol1) {
  259.       r = (x - w) * (fx - fv);
  260.       q = (x - v) * (fx - fw);
  261.       p = (x - v) * q - (x - w) * r;
  262.       q = 2.0 * (q - r);
  263.       if (q > 0.0) p = -p;
  264.       q = fabs(q);
  265.       etemp = e;
  266.       e = d;
  267.       if (fabs(p) >= fabs(0.5 * q * etemp) || p <= q * (a - x) || p >= q * (b - x))
  268.         d = CGOLD * (e = (x >= xm ? a - x : b - x));
  269.       else {
  270.         d = p / q;
  271.         u = x + d;
  272.         if (u - a < tol2 || b - u < tol2)
  273.           d = SIGN(tol1, xm - x);
  274.       }
  275.     }
  276.     else {
  277.       d = CGOLD * (e = (x >= xm ? a - x : b - x));
  278.     }
  279.     u = (fabs(d) >= tol1 ? x + d : x + SIGN(tol1, d));
  280.     fu = evfun1(f, u);
  281.     if (fu <= fx) {
  282.       if (u >= x) a = x; else b = x;
  283.       SHFT(v, w, x, u)
  284.       SHFT(fv, fw, fx, fu)
  285.     }
  286.     else {
  287.       if (u <  x) a = u; else b = u;
  288.       if (fu <= fw || w == x) {
  289.         v = w;
  290.         w = u;
  291.         fv = fw;
  292.         fw = fu;
  293.       }
  294.       else if (fu <= fv || v == x || v == w) {
  295.         v = u;
  296.         fv = fu;
  297.       }
  298.     }
  299.     if (verbose) {
  300.       sprintf(buf, "x = %f, f(x) = %f, a = %f, b = %f\n", x, fx, a, b);
  301.       stdputstr(buf);
  302.     }
  303.   }
  304.   
  305.   xlsave1(result);
  306.   result = mklist(3, NIL);
  307.   rplaca(result, cvflonum((FLOTYPE) x));
  308.   rplaca(cdr(result), cvflonum((FLOTYPE) fx));
  309.   rplaca(cdr(cdr(result)), cvfixnum((FIXTYPE) iter));
  310.   xlpop();
  311.   
  312.   return(result);
  313. }
  314. /*#endif OPTIMIZE*/