initial work towards getting xarray support into CLS.
[CommonLispStat.git] / lib / cbayes.c
blob8a842a4e0d59ae4df8cb8742cf72b56f12567edb
1 /* cbayes - Lisp interface to laplace approximation stuff */
2 /* Copyright (c) 1990, by Luke Tierney */
4 #include <stdlib.h> /* for calloc/realloc */
5 #include "linalg.h"
7 #ifdef INTPTR
8 typedef int PTR;
9 #else
10 typedef char *PTR;
11 #endif
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 *);
17 extern void maxfront();
18 extern void bufputstr(char *);
20 /************************************************************************/
21 /** **/
22 /** Definitions and Globals **/
23 /** **/
24 /************************************************************************/
26 #define MAXALLOC 20
28 static char *mem[MAXALLOC], memcount;
30 typedef struct {
31 int n, m, k, itnlimit, backtrack, verbose, vals_suppl, exptilt;
32 int count, termcode;
33 } MaxIPars;
35 typedef struct {
36 double typf, h, gradtol, steptol, maxstep, dflt, tilt, newtilt, hessadd;
37 } MaxDPars;
39 typedef struct {
40 MaxIPars max;
41 int full, covar;
42 } MomIPars;
44 typedef struct {
45 MaxDPars max;
46 double mgfdel;
47 } MomDPars;
49 /************************************************************************/
50 /** **/
51 /** Fake Replacements for S Interface **/
52 /** **/
53 /************************************************************************/
55 static void
56 meminit()
58 static int inited = FALSE;
59 int i;
61 if (! inited) {
62 for (i = 0; i < MAXALLOC; i++) mem[i] = nil;
63 inited = TRUE;
66 memcount = 0;
69 static int
70 makespace(void **pptr, size_t size) /* why are we using **char? */
72 if (size <= 0) {
73 return(1); /* we've done, by default, what we asked for. */
75 if (*pptr == nil) {
76 *pptr = calloc(size, 1);
77 } else {
78 *pptr = realloc(*pptr, size);
80 if (size > 0 && *pptr == nil) {
81 return(0); /* xlfail("memory allocation failed"); FIXME:AJR xlfail redef. */
83 return(1);
87 /************************************************************************/
88 /** **/
89 /** Callback Function **/
90 /** **/
91 /************************************************************************/
93 static void
94 callLminfun(size_t n,
95 double *x, double *fval, double *grad, double *hess,
96 long *derivs)
98 maximize_callback(n, (PTR) x,
99 (PTR) fval, (PTR) grad, (PTR) hess, (PTR) derivs);
102 void
103 call_S(char *fun, long narg, char **args, char **mode, long *length,char **names,
104 long nvals, char **values)
106 long n = length[0], derivs;
107 static double *fval = nil, *grad = nil, *hess = nil;
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));
116 callLminfun(n,(double *)args[0], fval, grad, hess, &derivs);
118 values[0] = (char *) fval;
119 values[1] = (derivs > 0) ? (char *) grad : nil;
120 values[2] = (derivs > 1) ? (char *) hess : nil;
123 void
124 Recover(s, w)
125 char *s, *w;
127 /* FIXME:AJR: xlfail(s); */
128 return;
131 /************************************************************************/
132 /** **/
133 /** Numerical Derivatives **/
134 /** **/
135 /************************************************************************/
137 void
138 numgrad_front(size_t n, PTR px, PTR pgrad, double h, PTR pscale)
140 LVAL f = nil;
141 double fval;
143 evalfront((char **)&f, &n, (double *) px,
144 &fval, (double *) pgrad, nil, &h, (double *) pscale);
147 void
148 numhess_front(size_t n, PTR px, PTR pf, PTR pgrad, PTR phess, double h, PTR pscale)
150 LVAL f = nil;
152 evalfront((char **)&f, &n, (double *) px,
153 (double *) pf, (double *) pgrad, (double *) phess,
154 &h, (double *) pscale);
157 /************************************************************************/
158 /** **/
159 /** Maximization Interface **/
160 /** **/
161 /************************************************************************/
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
178 static MaxIPars getMaxIPars(int *ipars)
180 MaxIPars ip;
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];
193 return(ip);
196 static MaxDPars getMaxDPars(dpars)
197 double *dpars;
199 MaxDPars dp;
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];
211 return(dp);
214 void
215 minfo_maximize(PTR px, PTR pfvals, PTR pscale, PTR pip, PTR pdp, int verbose)
217 LVAL f = nil;
218 MaxIPars ip;
219 MaxDPars dp;
220 int n, m, k;
221 static double *dx, *typx, *fvals;
222 char *msg;
224 dx = (double *) px;
225 typx = (double *) pscale;
226 fvals = (double *) pfvals;
228 ip = getMaxIPars((int *) pip);
229 dp = getMaxDPars((double *) pdp);
231 m = 0;
232 k = 0;
233 n = ip.n;
234 if (verbose >= 0) ip.verbose = verbose;
236 meminit();
237 maxfront(&f, nil, nil, dx, typx, fvals, nil, nil, nil,
238 &ip, &dp, nil, &msg);
240 bufputstr(msg);