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(size_t, PTR, PTR, PTR, PTR, PTR);
  14 extern void evalfront(char **, size_t *, 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, size_t 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(size_t n,
  95             double *x, double *fval, double *grad, double *hess,
  96             long *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   long 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(size_t n, PTR px, PTR pgrad, double h, PTR pscale)
 139 {
 140   LVAL f = nil;
 141   double fval;
 142
 143   evalfront((char **)&f, &n, (double *) px,
 144             &fval, (double *) pgrad, nil, &h, (double *) pscale);
 145 }
 146
 147 void
 148 numhess_front(size_t n, PTR px, PTR pf, PTR pgrad, PTR phess, double h, PTR pscale)
 149 {
 150   LVAL f = nil;
 151
 152   evalfront((char **)&f, &n, (double *) px,
 153             (double *) pf, (double *) pgrad, (double *) phess,
 154             &h, (double *) pscale);
 155 }
 156
 157 /************************************************************************/
 158 /**                                                                    **/
 159 /**                      Maximization Interface                        **/
 160 /**                                                                    **/
 161 /************************************************************************/
 162
 163 /* internals array information */
 164 #define INTLEN 12
 165 #define F_POS 0
 166 #define G_POS 1
 167 #define C_POS 2
 168 #define X_POS 3
 169 #define SCALE_POS 4
 170 #define FVALS_POS 5
 171 #define CVALS_POS 6
 172 #define CTARG_POS 7
 173 #define IPARS_POS 8
 174 #define DPARS_POS 9
 175 #define TSCAL_POS 10
 176 #define MULT_POS 11
 177
 178 static MaxIPars getMaxIPars(int *ipars)
 179 {
 180   MaxIPars ip;
 181
 182   ip.n = ipars[0];
 183   ip.m = ipars[1];
 184   ip.k = ipars[2];
 185   ip.itnlimit = ipars[3];
 186   ip.backtrack = ipars[4];
 187   ip.verbose = ipars[5];
 188   ip.vals_suppl = ipars[6];
 189   ip.exptilt = ipars[7];
 190   ip.count = ipars[8];
 191   ip.termcode = ipars[9];
 192
 193   return(ip);
 194 }
 195
 196 static MaxDPars getMaxDPars(dpars)
 197      double *dpars;
 198 {
 199   MaxDPars dp;
 200
 201   dp.typf = dpars[0];
 202   dp.h = dpars[1];
 203   dp.gradtol = dpars[2];
 204   dp.steptol = dpars[3];
 205   dp.maxstep = dpars[4];
 206   dp.dflt = dpars[5];
 207   dp.tilt = dpars[6];
 208   dp.newtilt = dpars[7];
 209   dp.hessadd = dpars[8];
 210
 211   return(dp);
 212 }
 213
 214 void
 215 minfo_maximize(PTR px, PTR pfvals, PTR pscale, PTR pip, PTR pdp, int verbose)
 216 {
 217   LVAL f = nil;
 218   MaxIPars ip;
 219   MaxDPars dp;
 220   int n, m, k;
 221   static double *dx, *typx, *fvals;
 222   char *msg;
 223
 224   dx = (double *) px;
 225   typx = (double *) pscale;
 226   fvals = (double *) pfvals;
 227
 228   ip = getMaxIPars((int *) pip);
 229   dp = getMaxDPars((double *) pdp);
 230
 231   m = 0;
 232   k = 0;
 233   n = ip.n;
 234   if (verbose >= 0) ip.verbose = verbose;
 235
 236   meminit();
 237   maxfront(&f, nil, nil, dx, typx, fvals, nil, nil, nil,
 238            &ip, &dp, nil, &msg);
 239
 240   bufputstr(msg);
 241 }