home *** CD-ROM | disk | FTP | other *** search
/ Fish 'n' More 2 / fishmore-publicdomainlibraryvol.ii1991xetec.iso / fish / applications / xlispstat / src / src2.lzh / XLisp-Stat / xsbayes.c < prev    next >
C/C++ Source or Header  |  1990-10-04  |  23KB  |  798 lines
  1. /* xsbayes - Lisp interface to laplace approximation stuff             */
  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. #include <stdlib.h>
  8. #include "xlisp.h"
  9. #include "osdef.h"
  10. #ifdef ANSI
  11. #include "xlproto.h"
  12. #include "xlsproto.h"
  13. #else
  14. #include "xlfun.h"
  15. #include "xlsfun.h"
  16. #endif ANSI
  17. #include "xlsvar.h"
  18.  
  19. #ifdef ANSI
  20. void callLminfun(LVAL,int,RVector,double *,RVector,RVector,int *),meminit(void),
  21.      double2list(int,double *,LVAL),seq2double(int,LVAL,RVector),
  22.      makespace(char **,int),seq2pchar(int,LVAL,LVAL *),setinternals(LVAL,LVAL),
  23.      call_scaled_Lfun(RVector,double *,RVector,RVector);
  24. LVAL newipars(MaxIPars),newdpars(MaxDPars),makecallarg(int),numderiv(int),
  25.      getinternals(LVAL),newinternals(LVAL,LVAL,LVAL,double),scaled_eval(int);
  26. MaxIPars getMaxIPars(LVAL);
  27. MaxDPars getMaxDPars(LVAL);
  28. #else
  29. void callLminfun(),meminit(),double2list(),seq2double(),
  30.      makespace(),seq2pchar(),setinternals(),call_scaled_Lfun();
  31. LVAL newipars(),newdpars(),makecallarg(),numderiv(),
  32.      getinternals(),newinternals(),scaled_eval();
  33. MaxIPars getMaxIPars();
  34. MaxDPars getMaxDPars();
  35. #endif ANSI
  36.  
  37. /************************************************************************/
  38. /**                                                                    **/
  39. /**                      Definitions and Globals                       **/
  40. /**                                                                    **/
  41. /************************************************************************/
  42.  
  43. #define MAXALLOC 20
  44.  
  45. static char  *mem[MAXALLOC], memcount;
  46. static LVAL arg;
  47. /* in xlsdef.h JKL
  48. typedef struct {
  49.   int n, m, k, itnlimit, backtrack, verbose, vals_suppl, exptilt;
  50.   int count, termcode;
  51. } MaxIPars;
  52.  
  53. typedef struct {
  54.   double typf, h, gradtol, steptol, maxstep, dflt, tilt, newtilt, hessadd;
  55. } MaxDPars;
  56.  
  57. typedef struct {
  58.   MaxIPars max;
  59.   int full, covar;
  60. } MomIPars;
  61.  
  62. typedef struct {
  63.   MaxDPars max;
  64.   double mgfdel;
  65. } MomDPars;
  66. */
  67. /************************************************************************/
  68. /**                                                                    **/
  69. /**                 Fake Replacements for S Interface                  **/
  70. /**                                                                    **/
  71. /************************************************************************/
  72.  
  73. static void meminit()
  74. {
  75.   static inited = FALSE;
  76.   int i;
  77.  
  78.   if (! inited) {
  79.     for (i = 0; i < MAXALLOC; i++) mem[i] = nil;
  80.     inited = TRUE;
  81.   }
  82.  
  83.   memcount = 0;
  84. }
  85.  
  86. static void makespace(pptr, size)
  87.      char **pptr;
  88.      int size;
  89. {
  90.   if (size <= 0) return;
  91.   if (*pptr == nil) *pptr = calloc(size, 1);
  92.   else *pptr = realloc(*pptr, size);
  93.   if (size > 0 && *pptr == nil) xlfail("memory allocation failed");
  94. }
  95.  
  96. #ifdef SBAYES
  97. char *S_alloc(n, m)
  98.      int n, m;
  99. {
  100.   char *result;
  101.   if (memcount >= MAXALLOC) result = nil;
  102.   else {
  103.     makespace(mem + memcount, n * m);
  104.     result = mem[memcount];
  105.     memcount++;
  106.   }
  107.   return(result);
  108. }
  109. #endif SBAYES
  110.  
  111. void call_S(fun, narg, args, mode, length, names, nvals, values)
  112.      LVAL fun; /* type changed  JKL */
  113.      char /* *fun,*/ **args, **mode, **names, **values;
  114.      long narg, nvals, *length;
  115. {
  116.   int n = length[0], derivs;
  117.   static double *fval = nil, *grad = nil, *hess = nil;
  118.  
  119.   makespace((char **)&fval, 1 * sizeof(double));  /* casts added  JKL */
  120.   makespace((char **)&grad, n * sizeof(double));
  121.   makespace((char **)&hess, n * n * sizeof(double));
  122.  
  123.   callLminfun(fun, n,(RVector)args[0], fval, grad, hess, &derivs);
  124.   
  125.   values[0] = (char *) fval;
  126.   values[1] = (derivs > 0) ? (char *) grad : nil;
  127.   values[2] = (derivs > 1) ? (char *) hess : nil;
  128. }
  129.  
  130. void Recover(s, w)
  131.      char *s, *w;
  132. {
  133.   xlfail(s);
  134. }
  135.  
  136. /************************************************************************/
  137. /**                                                                    **/
  138. /**                   Callback and Copying Function                    **/
  139. /**                                                                    **/
  140. /************************************************************************/
  141.  
  142. static void callLminfun(fun, n, x, fval, grad, hess, derivs)
  143.      LVAL fun;
  144.      int n, *derivs;
  145.      RVector x, grad, hess;
  146.      double *fval;
  147. {
  148.   LVAL result, seq;
  149.   
  150.   /* arg does not realy have to containt the function; this is a carryover */
  151.   double2list(n, x, car(cdr(arg)));
  152.   rplaca(arg, fun);
  153.   result = xlapply(pushargs(car(arg), cdr(arg)));
  154.  
  155.   if (realp(result)) {
  156.     *fval = makedouble(result);
  157.     *derivs = 0;
  158.   }
  159.   else if (consp(result)) {
  160.     *fval = makedouble(car(result));
  161.     *derivs = 0;
  162.     if (consp(cdr(result))) {
  163.       seq = compounddataseq(car(cdr(result)));
  164.       if (seqlen(seq) != n) xlerror("bad gradient", car(cdr(result)));
  165.       seq2double(n, seq, grad);
  166.       *derivs = 1;
  167.       if (consp(cdr(cdr(result)))) {
  168.     seq = compounddataseq(car(cdr(cdr(result))));
  169.     if (seqlen(seq) != n * n)
  170.       xlerror("bad hessian", car(cdr(cdr(result))));
  171.     seq2double(n * n, seq, hess);
  172.     *derivs = 2;
  173.       }
  174.     }
  175.   }
  176.   else xlerror("bad function result", result);
  177. }       
  178.   
  179. static LVAL makecallarg(n)
  180.      int n;
  181. {
  182.   LVAL carg;
  183.  
  184.   xlsave1(carg);
  185.   carg = mklist(n, NIL);
  186.   carg = consa(carg);
  187.   carg = cons(NIL, carg);
  188.   xlpop();
  189.   return(carg);
  190. }
  191.  
  192. static void seq2pchar(n, x, p)
  193.      int n;
  194.      LVAL x;
  195.      /*char **/LVAL *p;/* changed JKL */
  196. {
  197.   int i;
  198.  
  199.   for (i = 0; i < n; i++) 
  200.     p[i] = /*(char *)*/ getnextelement(&x, i);/* changed JKL */
  201. }
  202.  
  203. static void seq2double(n, x, dx)
  204.      int n;
  205.      LVAL x;
  206.      /*double **/RVector dx; /* changed JKL */
  207. {
  208.   int i;
  209.  
  210.   for (i = 0; i < n; i++) 
  211.     dx[i] = makedouble(getnextelement(&x, i));
  212. }
  213.  
  214. static void double2list(n, dx, x)
  215.      int n;
  216.      double *dx;
  217.      LVAL x;
  218. {
  219.   int i;
  220.  
  221.   for (i = 0; i < n && consp(x) ; i++, x = cdr(x))
  222.     rplaca(x, cvflonum((FLOTYPE) dx[i]));
  223. }
  224.  
  225. /************************************************************************/
  226. /**                                                                    **/
  227. /**                      Numerical Derivatives                         **/
  228. /**                                                                    **/
  229. /************************************************************************/
  230.  
  231. static LVAL numderiv(order)
  232.      int order;
  233. {
  234.   LVAL f, x, scale, result, data;
  235.   double h, fval;
  236.   static double *dx = nil, *grad = nil, *hess = nil, *typx = nil;
  237.   int n, i, all;
  238.  
  239.   f = xlgetarg();
  240.   x = xsgetsequence();
  241.   scale = (moreargs()) ? xsgetsequence() : NIL;
  242.   h = (moreargs()) ? makedouble(xlgetarg()) : -1.0;
  243.   all = (moreargs() && xlgetarg() != NIL) ? TRUE : FALSE;
  244.  
  245.   n = seqlen(x);
  246.   makespace((char **)&dx, n * sizeof(double));  /* casts added JKL */
  247.   makespace((char **)&grad, n * sizeof(double));
  248.   if (order > 1) makespace((char **)&hess, n * n * sizeof(double));
  249.   makespace((char **)&typx, n * sizeof(double));
  250.  
  251.   seq2double(n, x, dx);
  252.   if (scale != NIL) seq2double(n, scale, typx);
  253.   else for (i = 0; i < n; i++) typx[i] = 1.0;
  254.  
  255.   xlstkcheck(2);
  256.   xlsave(arg);
  257.   xlsave(result);
  258.   arg = makecallarg(n);
  259.  
  260.   evalfront(&f, &n, dx, &fval, grad, (order > 1) ? hess : nil, &h, typx);
  261.   
  262.   if (order == 1) {
  263.     result = mklist(n, NIL);
  264.     double2list(n, grad, result);
  265.   }
  266.   else {
  267.     result = integer_list_2(n, n);
  268.     result = newarray(result, NIL, NIL);
  269.     data = arraydata(result);
  270.     for (i = 0; i < n * n; i++)
  271.       setelement(data, i, cvflonum((FLOTYPE) hess[i]));
  272.     if (all) {
  273.       result = consa(result);
  274.       result = cons(mklist(n, NIL), result);
  275.       double2list(n, grad, car(result));
  276.       result = cons(cvflonum((FLOTYPE) fval), result);
  277.     }
  278.   }
  279.   xlpopn(2);
  280.  
  281.   return(result);
  282. }
  283.  
  284. LVAL xsnumgrad() { return(numderiv(1)); }
  285. LVAL xsnumhess() { return(numderiv(2)); }
  286.  
  287. /************************************************************************/
  288. /**                                                                    **/
  289. /**                      Maximization Interface                        **/
  290. /**                                                                    **/
  291. /************************************************************************/
  292.  
  293. #define INTLEN 12
  294. #define F_POS 0
  295. #define G_POS 1
  296. #define C_POS 2
  297. #define X_POS 3
  298. #define SCALE_POS 4
  299. #define FVALS_POS 5
  300. #define CVALS_POS 6
  301. #define CTARG_POS 7
  302. #define IPARS_POS 8
  303. #define DPARS_POS 9
  304. #define TSCAL_POS 10
  305. #define MULT_POS 11
  306.  
  307. #define getffun(i) getelement(i, F_POS)
  308. #define getgfuns(i) getelement(i, G_POS)
  309. #define getcfuns(i) getelement(i, C_POS)
  310. #define getx(i) getelement(i, X_POS)
  311. #define getscale(i) getelement(i, SCALE_POS)
  312. #define getfvals(i) getelement(i, FVALS_POS)
  313. #define getcvals(i) getelement(i, CVALS_POS)
  314. #define getctarget(i) getelement(i, CTARG_POS)
  315. #define getipars(i) getelement(i, IPARS_POS)
  316. #define getdpars(i) getelement(i, DPARS_POS)
  317. #define gettscale(i) getelement(i, TSCAL_POS)
  318. #define getmults(i) getelement(i, MULT_POS)
  319. #define getderivstep(i) getflonum(getelement(getdpars(i), 1))
  320.  
  321. #define setx(i, v) setelement(i, X_POS, v)
  322. #define setscale(i, v) setelement(i, SCALE_POS, v)
  323. #define setmults(i, v) setelement(i, MULT_POS, v)
  324. #define setfvals(i, v) setelement(i, FVALS_POS, v)
  325. #define setcvals(i, v) setelement(i, CVALS_POS, v)
  326. #define setipars(i, v) setelement(i, IPARS_POS, v)
  327. #define setdpars(i, v) setelement(i, DPARS_POS, v)
  328. #define setderivstep(i, v) setelement(getdpars(i), 1, cvflonum((FLOTYPE) (v)))
  329.  
  330. static LVAL getinternals(minfo)
  331.      LVAL minfo;
  332. {
  333.   return(slot_value(minfo, xlenter("INTERNALS")));
  334. }
  335.  
  336. static void setinternals(minfo, internals)
  337.      LVAL minfo, internals;
  338. {
  339.   set_slot_value(minfo, xlenter("INTERNALS"), internals);
  340. }
  341.  
  342. static LVAL newipars(ip)
  343.      MaxIPars ip;
  344. {
  345.   LVAL result;
  346.  
  347.   xlsave1(result);
  348.   result = newvector(10);
  349.   setelement(result, 0, cvfixnum((FIXTYPE) ip.n));
  350.   setelement(result, 1, cvfixnum((FIXTYPE) ip.m));
  351.   setelement(result, 2, cvfixnum((FIXTYPE) ip.k));
  352.   setelement(result, 3, cvfixnum((FIXTYPE) ip.itnlimit));
  353.   setelement(result, 4, cvfixnum((FIXTYPE) ip.backtrack));
  354.   setelement(result, 5, cvfixnum((FIXTYPE) ip.verbose));
  355.   setelement(result, 6, cvfixnum((FIXTYPE) ip.vals_suppl));
  356.   setelement(result, 7, cvfixnum((FIXTYPE) ip.exptilt));
  357.   setelement(result, 8, cvfixnum((FIXTYPE) ip.count));
  358.   setelement(result, 9, cvfixnum((FIXTYPE) ip.termcode));
  359.   xlpop();
  360.   return(result);
  361. }
  362.  
  363. static LVAL newdpars(dp)
  364.      MaxDPars dp;
  365. {
  366.   LVAL result;
  367.  
  368.   xlsave1(result);
  369.   result = newvector(9);
  370.   setelement(result, 0, cvflonum((FLOTYPE) dp.typf));
  371.   setelement(result, 1, cvflonum((FLOTYPE) dp.h));
  372.   setelement(result, 2, cvflonum((FLOTYPE) dp.gradtol));
  373.   setelement(result, 3, cvflonum((FLOTYPE) dp.steptol));
  374.   setelement(result, 4, cvflonum((FLOTYPE) dp.maxstep));
  375.   setelement(result, 5, cvflonum((FLOTYPE) dp.dflt));
  376.   setelement(result, 6, cvflonum((FLOTYPE) dp.tilt));
  377.   setelement(result, 7, cvflonum((FLOTYPE) dp.newtilt));
  378.   setelement(result, 8, cvflonum((FLOTYPE) dp.hessadd));
  379.   xlpop();
  380.   return(result);
  381. }
  382.  
  383. static MaxIPars getMaxIPars(internals)
  384.      LVAL internals;
  385. {
  386.   LVAL ipars = getipars(internals);
  387.   MaxIPars ip;
  388.  
  389.   ip.n = getfixnum(getelement(ipars, 0));
  390.   ip.m = getfixnum(getelement(ipars, 1));
  391.   ip.k = getfixnum(getelement(ipars, 2));
  392.   ip.itnlimit = getfixnum(getelement(ipars, 3));
  393.   ip.backtrack = getfixnum(getelement(ipars, 4));
  394.   ip.verbose = getfixnum(getelement(ipars, 5));
  395.   ip.vals_suppl = getfixnum(getelement(ipars, 6));
  396.   ip.exptilt = getfixnum(getelement(ipars, 7));
  397.   ip.count = getfixnum(getelement(ipars, 8));
  398.   ip.termcode = getfixnum(getelement(ipars, 9));
  399.  
  400.   return(ip);
  401. }
  402.  
  403. static MaxDPars getMaxDPars(internals)
  404.      LVAL internals;
  405. {
  406.   LVAL dpars = getdpars(internals);
  407.   MaxDPars dp;
  408.  
  409.   dp.typf = getflonum(getelement(dpars, 0));
  410.   dp.h = getflonum(getelement(dpars, 1));
  411.   dp.gradtol = getflonum(getelement(dpars, 2));
  412.   dp.steptol = getflonum(getelement(dpars, 3));
  413.   dp.maxstep = getflonum(getelement(dpars, 4));
  414.   dp.dflt = getflonum(getelement(dpars, 5));
  415.   dp.tilt = getflonum(getelement(dpars, 6));
  416.   dp.newtilt = getflonum(getelement(dpars, 7));
  417.   dp.hessadd = getflonum(getelement(dpars, 8));
  418.  
  419.   return(dp);
  420. }
  421.  
  422. static LVAL newinternals(f, x, scale, h)
  423.      LVAL f, x, scale;
  424.      double h;
  425. {
  426.   LVAL result;
  427.   int n, i;
  428.   MaxIPars ip;
  429.   MaxDPars dp;
  430.  
  431.   n = seqlen(x);
  432.   ip.n = n; ip.k = 0; ip.m = 0; ip.itnlimit = -1; ip.backtrack = TRUE;
  433.   ip.verbose = 0; ip.vals_suppl = FALSE; ip.exptilt = TRUE;
  434.   ip.count = 0; ip.termcode = 0;
  435.  
  436.   dp.typf = 1.0; dp.h = h; dp.gradtol = -1.0; dp.steptol = -1.0;
  437.   dp.maxstep = -1.0; dp.dflt = 0.0; dp.tilt = 0.0; dp.newtilt = 0.0;
  438.   dp.hessadd = 0.0;
  439.  
  440.   xlsave1(result);
  441.   result = newvector(INTLEN);
  442.   setelement(result, F_POS, f);
  443.   setelement(result, X_POS, x);
  444.   if (scale != NIL) setelement(result, SCALE_POS, scale);
  445.   else {
  446.     scale = newvector(n);
  447.     setelement(result, SCALE_POS, scale);
  448.     for (i = 0; i < n; i++) setelement(scale, i, cvflonum((FLOTYPE) 1.0));
  449.   }
  450.   setelement(result, IPARS_POS, newipars(ip));
  451.   setelement(result, DPARS_POS, newdpars(dp));
  452.   xlpop();
  453.   
  454.   return(result);
  455. }
  456.  
  457. LVAL xsminfo_isnew()
  458. {
  459.   LVAL myinfo, f, x, scale, step;
  460.   double h;
  461.  
  462.   myinfo = xlgaobject();
  463.   f = xlgetarg();
  464.   x = xsgetsequence();
  465.   if (! xlgetkeyarg(xlenter(":SCALE"), &scale)) scale = NIL;
  466.   if (xlgetkeyarg(xlenter(":DERIVSTEP"), &step)) h = makedouble(step);
  467.   else h = -1.0;
  468.  
  469.   setinternals(myinfo, newinternals(f, x, scale, h));
  470.   return(NIL);
  471. }
  472.  
  473. LVAL xsminfo_maximize()
  474. {
  475.   LVAL minfo, internals, f, g, c, x, scale, verbose, tiltscale, result, mults; 
  476.   MaxIPars ip;
  477.   MaxDPars dp;
  478.   int n, m, k;
  479.   RVector dx = nil, typx = nil; /* changed JKL */
  480.   static double /* *dx = nil, *typx = nil, */ *fvals = nil, *tscale = nil;
  481.   static double *cvals = nil, *ctarget = nil;
  482. /*  static char *gg = nil, *cc = nil;*/
  483.   static LVAL gg = nil, cc = nil;/* changed JKL */
  484.   char *msg;
  485.  
  486.   minfo = xlgaobject();
  487.   verbose = (moreargs()) ? xlgafixnum() : NIL;
  488.  
  489.   internals = getinternals(minfo);
  490.   f = getffun(internals);
  491.   g = getgfuns(internals);
  492.   c = getcfuns(internals);
  493.   x = getx(internals);
  494.   scale = getscale(internals);
  495.   ip = getMaxIPars(internals);
  496.   dp = getMaxDPars(internals);
  497.   tiltscale = gettscale(internals);
  498.   mults = getmults(internals);
  499.  
  500.   m = ip.m;
  501.   k = ip.k;
  502.   if (verbose != NIL) ip.verbose = getfixnum(verbose);
  503.   n = seqlen(x);
  504.   ip.n = n;
  505.   makespace((char **)&gg, m * sizeof(char *));  /* casts added  JKL */
  506.   makespace((char **)&cc, k * sizeof(char *));
  507.   makespace((char **)&dx, (n + k) * sizeof(double));
  508.   makespace((char **)&typx, n * sizeof(double));
  509.   makespace((char **)&fvals, (1 + n + n * n) * sizeof(double));
  510.   makespace((char **)&tscale, m * sizeof(double));
  511.   makespace((char **)&cvals, (k + k * n) * sizeof(double));
  512.   makespace((char **)&ctarget, k * sizeof(double));
  513.  
  514.   seq2pchar(m, g, &gg);  /* pointers added JKL */
  515.   seq2pchar(k, c, &cc);
  516.   seq2double(n, x, dx);
  517.   seq2double(n, scale, typx);
  518.   if (ip.vals_suppl) {
  519.     seq2double(1 + n + n * n, getfvals(internals), fvals);
  520.     seq2double(k + k * n, getcvals(internals), cvals);
  521.   }
  522.   seq2double(m, tiltscale, tscale);
  523.   seq2double(k, getctarget(internals), ctarget);
  524.   seq2double(k, mults, dx + n);
  525.  
  526.   xlsave1(arg);
  527.   arg = makecallarg(n);
  528.  
  529.   meminit();  /* pointers added JKL */
  530.   maxfront(&f, &gg, &cc, dx, typx, fvals, nil, cvals, ctarget, &ip, &dp, tscale, &msg);
  531.  
  532.   result = mklist(n, NIL);
  533.   setx(internals, result);
  534.   double2list(n, dx, result);
  535.   result = mklist(1 + n + n * n, NIL);
  536.   setfvals(internals, result);
  537.   double2list(1 + n + n * n, fvals, result);
  538.   if (k > 0) {
  539.     result = mklist(k + k * n, NIL);
  540.     setcvals(internals, result);
  541.     double2list(k + k * n, cvals, result);
  542.     result = mklist(k, NIL);
  543.     setmults(internals, result);
  544.     double2list(k, dx + n, result);
  545.   }
  546.   setipars(internals, newipars(ip));
  547.   setdpars(internals, newdpars(dp));
  548.  
  549.   xlpop();
  550.   
  551.   return(make_string(msg));
  552. }
  553.  
  554. /************************************************************************/
  555. /**                                                                    **/
  556. /**                        Laplace Interface                           **/
  557. /**                                                                    **/
  558. /************************************************************************/
  559.  
  560. LVAL xsminfo_loglap()
  561. {
  562.   LVAL minfo, internals;
  563.   MaxIPars ip;
  564.   MaxDPars dp;
  565.   int n, k, detonly;
  566.   double val;
  567.   static double *fvals = nil, *cvals = nil;
  568.   
  569.   minfo = xlgaobject();
  570.   detonly = (moreargs()) ? (xlgetarg() != NIL) : FALSE;
  571.  
  572.   internals = getinternals(minfo);
  573.   ip = getMaxIPars(internals);
  574.   dp = getMaxDPars(internals);
  575.  
  576.   n = ip.n;
  577.   k = ip.k;
  578.             /* casts added JKL */
  579.   makespace((char **)&fvals, (1 + n + n * n) * sizeof(double));
  580.   makespace((char **)&cvals, (k + k * n) * sizeof(double));
  581.   seq2double(1 + n + n * n, getfvals(internals), fvals);
  582.   seq2double(k + k * n, getcvals(internals), cvals);
  583.  
  584.   loglapdet(fvals, cvals, &ip, &dp, &val, &detonly);
  585.  
  586.   return(cvflonum((FLOTYPE) val));
  587. }
  588.  
  589. /************************************************************************/
  590. /**                                                                    **/
  591. /**            Scaled Evaluation and Numerical Derivatives             **/
  592. /**                                                                    **/
  593. /************************************************************************/
  594.  
  595. #ifdef DODO
  596. /**** function objects; exact derivatives; NIL values ****/
  597.  
  598. static struct scaled_Lfun_info {
  599.   LVAL arglist, value;
  600.   RVector scaling, z;
  601.   int dim, in_range, num_derivs, cached;
  602.   double cached_value;
  603. } Lfun_info;
  604.  
  605. static LVAL scaled_eval(order)
  606.     int order;
  607. {
  608.   struct scaled_Lfun_info old_info;
  609.   LVAL Lfun, Lx, Lscaling, Lvalue, space, hspace, Ldata;
  610.   RVector x, fsum, grad;
  611.   RMatrix hess;
  612.   double value, h;
  613.   int k, i, j, all;
  614.   
  615.   old_info = Lfun_info;
  616.   
  617.   Lfun = xlgetarg();
  618.   Lx = xsgetsequence();
  619.   Lscaling = xsgetsequence();
  620.   if (order > 0) h = makedouble(xlgetarg());
  621.   if (order == 2) all = (moreargs() && xlgetarg() != NIL) ? TRUE : FALSE;
  622.   xllastarg();
  623.   
  624.   xlstkcheck(5);
  625.   xlsave(space);
  626.   xlsave(hspace);
  627.   xlsave(Lfun_info.arglist);
  628.   xlsave(Lfun_info.value);
  629.   xlsave(Lvalue);
  630.   
  631.   k = Lfun_info.dim = seqlen(Lx);
  632.   space = newadata(sizeof(double), 2 + k * (2 * k + 5), FALSE);
  633.   Lfun_info.scaling = (RVector) getadaddr(space);
  634.   Lfun_info.z = Lfun_info.scaling + 2 + k * (k + 1);
  635.   x = Lfun_info.scaling + 2 + k * (k + 2);
  636.   fsum = Lfun_info.scaling + 2 + k * (k + 3);
  637.   grad = Lfun_info.scaling + 2 + k * (k + 4);
  638.   Lfun_info.arglist = mklist(2, NIL);
  639.   rplaca(Lfun_info.arglist, Lfun);
  640.   rplaca(cdr(Lfun_info.arglist), mklist(k, NIL));
  641.   Lfun_info.in_range = TRUE;
  642.   if (order == 2) {
  643.     hspace = newadata(sizeof(double *), k, FALSE);
  644.     hess = (RMatrix) getadaddr(hspace);
  645.     for (i = 0; i < k; i++)
  646.       hess[i] = Lfun_info.scaling + 2 + k * (k + 5 + i);
  647.   }
  648.   else hess = nil;
  649.   Lfun_info.cached = FALSE;
  650.   Lfun_info.num_derivs = 0;
  651.   
  652.   seq2double(2 + k * (k + 1), Lscaling, Lfun_info.scaling);
  653.   seq2double(k, Lx, x);
  654.  
  655.   call_scaled_Lfun(x, &value, grad, hess);
  656.   switch (order) {
  657.   case 0:
  658.     if (Lfun_info.in_range) Lvalue = cvflonum((FLOTYPE) value);
  659.     else Lvalue = Lfun_info.value;
  660.     break;
  661.   case 1:
  662.     if (Lfun_info.num_derivs < 1)
  663.       numergrad(k, x, grad, fsum, call_scaled_Lfun, h, nil);
  664.     if (Lfun_info.in_range) {
  665.       Lvalue = mklist(k, NIL);
  666.       double2list(k, grad, Lvalue);
  667.     }
  668.     else Lvalue = NIL;
  669.     break;
  670.   case 2:
  671.     /* currently uses second differences even for analytic gradient */
  672.     if (Lfun_info.num_derivs < 2) {
  673.       numergrad(k, x, grad, fsum, call_scaled_Lfun, h, nil);
  674.       /* Lfun_info.cached_value = value;
  675.       Lfun_info.cached = TRUE;*/
  676.       numerhess(k, x, hess, value, fsum, call_scaled_Lfun, h, nil);
  677.     }
  678.     if (Lfun_info.in_range) {
  679.       Lvalue = integer_list_2(k, k);
  680.       Lvalue = newarray(Lvalue, NIL, NIL);
  681.       Ldata = arraydata(Lvalue);
  682.       for (i = 0; i < k; i++)
  683.         for (j = 0; j < k; j++)
  684.           setelement(Ldata, i * k + j, cvflonum((FLOTYPE) hess[i][j]));
  685.       if (all) {
  686.         Lvalue = consa(Lvalue);
  687.         Lvalue = cons(mklist(k, NIL), Lvalue);
  688.         double2list(k, grad, car(Lvalue));
  689.         Lvalue = cons(cvflonum((FLOTYPE) value), Lvalue);
  690.       }
  691.     }
  692.     else Lvalue = NIL;
  693.     break;
  694.   }
  695.   
  696.   xlpopn(5);
  697.   Lfun_info = old_info;
  698.   return(Lvalue);
  699. }
  700.  
  701. static void call_scaled_Lfun(x, value, grad, hess)
  702.     RVector x, grad, hess;
  703.     double *value;
  704. {
  705.   LVAL temp, Lvalue, Lgrad, Lhess;
  706.   double center, scale;
  707.   RVector mean, sigma, z;
  708.   int i, j, k;
  709.   
  710.   /* cheat to avoid recalculating function value in Hessian evaluation */
  711.   if (Lfun_info.cached && Lfun_info.num_derivs == 0) {
  712.     Lfun_info.cached = FALSE;
  713.     *value = Lfun_info.cached_value;
  714.     return;
  715.   }
  716.   
  717.   k = Lfun_info.dim;
  718.   center = Lfun_info.scaling[0];
  719.   scale = Lfun_info.scaling[1];
  720.   if (scale == 0.0) xlfail("division by zero");
  721.   mean = Lfun_info.scaling + 2;
  722.   sigma = Lfun_info.scaling + 2 + k;
  723.   z = Lfun_info.z;
  724.   
  725.   for (i = 0; i < k; i++) {
  726.     z[i] = mean[i];
  727.     for (j = 0; j <= i; j++)
  728.       z[i] += sigma[i * k + j] * x[j];
  729.   }
  730.   
  731.   double2list(k, z, car(cdr(Lfun_info.arglist)));
  732.  
  733.   temp = xlapply(pushargs(car(Lfun_info.arglist), cdr(Lfun_info.arglist)));
  734.   
  735.   Lfun_info.num_derivs = (consp(temp)) ? seqlen(temp) - 1 : 0;
  736.   
  737.   Lvalue = (consp(temp)) ? car(Lvalue) : temp;
  738.   Lgrad = (consp(temp) && consp(cdr(temp))) ? car(cdr(temp)) : NIL;
  739.   Lhess = (consp(temp) && consp(cdr(temp)) && consp(cdr(cdr(temp)))) 
  740.         ? car(cdr(cdr(temp))) : NIL;
  741.   
  742.   if (fixp(Lvalue) || floatp(Lvalue)) {
  743.       *value = (makedouble(Lvalue) - center) / scale;
  744.     if (grad != nil && Lgrad != NIL) {
  745.       seq2double(k, Lgrad, grad);
  746.       for (i = 0; i < k; i++) {
  747.         z[i] = 0.0;
  748.         for (j = i; j < k; j++) z[i] += sigma[j * k + i] * grad[j];
  749.         grad[i] = z[i] / scale;
  750.       }
  751.     }
  752.     if (hess != nil && Lhess != NIL) {
  753.       if (! matrixp(Lhess)) xlerror("not an a matrix", Lhess);
  754.       seq2double(k * k, arraydata(Lhess), hess);
  755.       /* rescale */
  756.     }
  757.   }
  758.   else {
  759.     *value = 0.0;
  760.     Lfun_info.in_range = FALSE;
  761.     Lfun_info.value = Lvalue;
  762.   }
  763. }
  764.  
  765. LVAL xsscaled_c2_eval() { return(scaled_eval(0)); }
  766. LVAL xsscaled_c2_grad() { return(scaled_eval(1)); }
  767. LVAL xsscaled_c2_hess() { return(scaled_eval(2)); }
  768. #endif DODO
  769.  
  770. LVAL xsaxpy()
  771. {
  772.   LVAL result, next, tx, a, x, y;
  773.   int i, j, n, start, end, lower;
  774.   double val;
  775.   
  776.   a = arraydata(xsgetmatrix());
  777.   x = xsgetsequence();
  778.   y = xsgetsequence();
  779.   lower = (moreargs() && xlgetarg() != NIL) ? TRUE : FALSE;
  780.   
  781.   n = seqlen(x);
  782.   if (n * n != seqlen(a) || n != seqlen(y)) xlfail("dimensions do not match");
  783.   
  784.   xlsave(result);
  785.   result = mklist(n, NIL);
  786.   for (i = 0, start = 0, next = result; i < n; i++, start += n, next = cdr(next)) {
  787.       val = makedouble(getnextelement(&y, i));
  788.     end = (lower) ? i : n - 1;
  789.     for (j = 0, tx = x; j <= end; j++) {
  790.       val += makedouble(getnextelement(&tx, j)) 
  791.            * makedouble(getelement(a, start + j));
  792.     }
  793.     rplaca(next, cvflonum((FLOTYPE) val));
  794.   }
  795.   xlpop();
  796.   return(result);
  797. }
  798.