lib/cbayes.c

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