lib/cbayes.c

   1 /* cbayes - Lisp interface to laplace approximation stuff              */
   2 /* Copyright (c) 1990, by Luke Tierney                                 */
   3
   4 #include "linalg.h"
   5
   6 #ifdef INTPTR
   7 typedef int PTR;
   8 #else
   9 typedef char *PTR;
  10 #endif
  11
  12 extern char *calloc(), *realloc();
  13
  14 /************************************************************************/
  15 /**                                                                    **/
  16 /**                      Definitions and Globals                       **/
  17 /**                                                                    **/
  18 /************************************************************************/
  19
  20 #define MAXALLOC 20
  21
  22 static char  *mem[MAXALLOC], memcount;
  23
  24 typedef struct {
  25   int n, m, k, itnlimit, backtrack, verbose, vals_suppl, exptilt;
  26   int count, termcode;
  27 } MaxIPars;
  28
  29 typedef struct {
  30   double typf, h, gradtol, steptol, maxstep, dflt, tilt, newtilt, hessadd;
  31 } MaxDPars;
  32
  33 typedef struct {
  34   MaxIPars max;
  35   int full, covar;
  36 } MomIPars;
  37
  38 typedef struct {
  39   MaxDPars max;
  40   double mgfdel;
  41 } MomDPars;
  42
  43 /************************************************************************/
  44 /**                                                                    **/
  45 /**                 Fake Replacements for S Interface                  **/
  46 /**                                                                    **/
  47 /************************************************************************/
  48
  49 static meminit()
  50 {
  51   static inited = FALSE;
  52   int i;
  53
  54   if (! inited) {
  55     for (i = 0; i < MAXALLOC; i++) mem[i] = nil;
  56     inited = TRUE;
  57   }
  58
  59   memcount = 0;
  60 }
  61
  62 static makespace(pptr, size)
  63      char **pptr;
  64      int size;
  65 {
  66   if (size <= 0) return;
  67   if (*pptr == nil) *pptr = calloc(size, 1);
  68   else *pptr = realloc(*pptr, size);
  69   if (size > 0 && *pptr == nil) xlfail("memory allocation failed");
  70 }
  71
  72 call_S(fun, narg, args, mode, length, names, nvals, values)
  73      char *fun, **args, **mode, **names, **values;
  74      long narg, nvals, *length;
  75 {
  76   int n = length[0], derivs;
  77   static double *fval = nil, *grad = nil, *hess = nil;
  78
  79   makespace(&fval, 1 * sizeof(double));
  80   makespace(&grad, n * sizeof(double));
  81   makespace(&hess, n * n * sizeof(double));
  82
  83   callLminfun(n, args[0], fval, grad, hess, &derivs);
  84
  85   values[0] = (char *) fval;
  86   values[1] = (derivs > 0) ? (char *) grad : nil;
  87   values[2] = (derivs > 1) ? (char *) hess : nil;
  88 }
  89
  90 Recover(s, w)
  91      char *s, *w;
  92 {
  93   xlfail(s);
  94 }
  95
  96 /************************************************************************/
  97 /**                                                                    **/
  98 /**                         Callback Function                          **/
  99 /**                                                                    **/
 100 /************************************************************************/
 101
 102 static callLminfun(n, x, fval, grad, hess, derivs)
 103      int n, *derivs;
 104      RVector x, grad, hess;
 105      double *fval;
 106 {
 107   maximize_callback(n, (PTR) x,
 108                     (PTR) fval, (PTR) grad, (PTR) hess, (PTR) derivs);
 109 }
 110
 111 /************************************************************************/
 112 /**                                                                    **/
 113 /**                      Numerical Derivatives                         **/
 114 /**                                                                    **/
 115 /************************************************************************/
 116
 117 numgrad_front(n, px, pgrad, h, pscale)
 118      int n;
 119      PTR px, pgrad, pscale;
 120      double h;
 121 {
 122   LVAL f = nil;
 123   double fval;
 124
 125   evalfront(&f, &n, (double *) px,
 126             &fval, (double *) pgrad, nil, &h, (double *) pscale);
 127 }
 128
 129 numhess_front(n, px, pf, pgrad, phess, h, pscale)
 130      int n;
 131      PTR px, pf, pgrad, phess, pscale;
 132      double h;
 133 {
 134   LVAL f = nil;
 135
 136   evalfront(&f, &n, (double *) px,
 137             (double *) pf, (double *) pgrad, (double *) phess,
 138             &h, (double *) pscale);
 139 }
 140
 141 /************************************************************************/
 142 /**                                                                    **/
 143 /**                      Maximization Interface                        **/
 144 /**                                                                    **/
 145 /************************************************************************/
 146
 147 /* internals array information */
 148 #define INTLEN 12
 149 #define F_POS 0
 150 #define G_POS 1
 151 #define C_POS 2
 152 #define X_POS 3
 153 #define SCALE_POS 4
 154 #define FVALS_POS 5
 155 #define CVALS_POS 6
 156 #define CTARG_POS 7
 157 #define IPARS_POS 8
 158 #define DPARS_POS 9
 159 #define TSCAL_POS 10
 160 #define MULT_POS 11
 161
 162 static MaxIPars getMaxIPars(ipars)
 163      int *ipars;
 164 {
 165   MaxIPars ip;
 166
 167   ip.n = ipars[0];
 168   ip.m = ipars[1];
 169   ip.k = ipars[2];
 170   ip.itnlimit = ipars[3];
 171   ip.backtrack = ipars[4];
 172   ip.verbose = ipars[5];
 173   ip.vals_suppl = ipars[6];
 174   ip.exptilt = ipars[7];
 175   ip.count = ipars[8];
 176   ip.termcode = ipars[9];
 177
 178   return(ip);
 179 }
 180
 181 static MaxDPars getMaxDPars(dpars)
 182      double *dpars;
 183 {
 184   MaxDPars dp;
 185
 186   dp.typf = dpars[0];
 187   dp.h = dpars[1];
 188   dp.gradtol = dpars[2];
 189   dp.steptol = dpars[3];
 190   dp.maxstep = dpars[4];
 191   dp.dflt = dpars[5];
 192   dp.tilt = dpars[6];
 193   dp.newtilt = dpars[7];
 194   dp.hessadd = dpars[8];
 195
 196   return(dp);
 197 }
 198
 199 minfo_maximize(px, pfvals, pscale, pip, pdp, verbose)
 200      PTR px, pfvals, pscale, pip, pdp;
 201      int verbose;
 202 {
 203   LVAL f = nil;
 204   MaxIPars ip;
 205   MaxDPars dp;
 206   int n, m, k;
 207   static double *dx, *typx, *fvals;
 208   char *msg;
 209
 210   dx = (double *) px;
 211   typx = (double *) pscale;
 212   fvals = (double *) pfvals;
 213
 214   ip = getMaxIPars((int *) pip);
 215   dp = getMaxDPars((double *) pdp);
 216
 217   m = 0;
 218   k = 0;
 219   n = ip.n;
 220   if (verbose >= 0) ip.verbose = verbose;
 221
 222   meminit();
 223   maxfront(&f, nil, nil, dx, typx, fvals, nil, nil, nil,
 224            &ip, &dp, nil, &msg);
 225
 226   bufputstr(msg);
 227 }