move stat stuff to statistics, clean up maximize
[CommonLispStat.git] / lib / functions.c
bloba6b42ecbabc16d6166c996357c7f54f12c8bef87
1 /* functions - callbacks for Bayes routines in XLISP-STAT and S */
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. */
7 #include <stdio.h>
8 #include "xmath.h"
10 #define PRINTSTR(s) printf(s)
12 extern void *S_alloc(), *calloc(), *realloc();
13 extern double macheps();
14 char *minresultstring();
16 extern void choldecomp();
18 extern void minsupplyvalues(double , double *, double **,
19 double *, double **);
20 extern void minimize();
21 extern void minresults(double *, double *, double *, double *, double **,
22 double *, double **, int *, int *, double *);
23 extern void minsetup(int n, int k,
24 int *ffun(), int *gfun(),
25 double *x, double typf, double *typx, char *work);
26 /* int (*ffun)(), (*gfun)();*/
27 extern void minsetoptions(double , double, double, int , int , int , int);
29 extern void choldecomp();
33 /* next 2 from cbayes.c */
34 extern void Recover(char *, char *);
35 extern void call_S(char *fun, long narg, char **args, char **mode, long *length,char **names,
36 long nvals, char **values);
38 /* next 2 from derivatives.c */
40 extern void numergrad(int n, double *x, double *grad, double *fsum,
41 int *ffun(), double h, double *typx);
42 extern void numerhess(int n, double *x, double **hess, double f,
43 double *fsum, void ffun(),
44 double h, double *typx);
46 extern void numergrad(int , double *, double *, double *,
47 int ffun(double *, double *, double *, double **),
48 double , double *);
49 extern void numerhess(int , double *, double **, double ,
50 double *,
51 int ffun(double *, double *, double *, double **),
52 double , double *);
54 /* next from minimize */
55 extern int minworkspacesize(int, int);
57 /************************************************************************/
58 /** **/
59 /** Definitions and Globals **/
60 /** **/
61 /************************************************************************/
63 /* #define NULL 0L already defined in stddef */
65 #define nil 0L
66 #define TRUE 1
67 #define FALSE 0
69 #define ROOT2PI 2.50662827463100050241
70 #define PI_INV .31830988618379067153
72 #define GRADTOL_POWER 1.0 / 3.0
73 #define H_POWER 1.0 / 6.0
75 /*typedef double **RMatrix, *RVector;*/
77 typedef struct{
78 char *f, **sf, **g;
79 int n, k;
80 int change_sign, fderivs;
81 int *gderivs;
82 double typf, h, dflt;
83 double *typx, *fsum, *cvals, *ctarget;
84 double **gfsum;
85 } Fundata;
87 static Fundata func, gfuncs, cfuncs;
89 /************************************************************************/
90 /** **/
91 /** Memory Utilities **/
92 /** **/
93 /************************************************************************/
95 /* this function is used to maintain a statically allocated piece of */
96 /* memory of a specified size. If a larger piece is needed the pointer */
97 /* is realloced. This allows functions using memory allocation to be */
98 /* called repeatedly (but not recursively) from within the same call */
99 /* from S. It attempts to avoid the danger of dangling callocs. */
101 void static
102 makespace(pptr, size)
103 char **pptr;
104 int size;
106 if (size <= 0) return;
107 if (*pptr == nil) *pptr = calloc(size, 1);
108 else *pptr = realloc(*pptr, size);
109 if (size > 0 && *pptr == nil) Recover("memory allocation failed", NULL);
112 /************************************************************************/
113 /** **/
114 /** Functions Evaluation Routines **/
115 /** **/
116 /************************************************************************/
119 * All Hessian evaluations by numerical derivatives assume the gradient is
120 * evaluated first at the same location. The results are cached away.
123 /* install log posterior function */
124 static void
125 install_func(char *f, char **sf, int n, int change_sign, /*?? */
126 double typf, double h, double *typx, double dflt)
128 int i;
129 static int inited = FALSE;
131 if (! inited) {
132 func.typx = nil;
133 func.fsum = nil;
134 inited = TRUE;
136 makespace(&func.typx, n * sizeof(double));
137 makespace(&func.fsum, n * sizeof(double));
139 func.f = f;
140 func.sf = sf;
141 func.n = n;
142 func.change_sign = change_sign;
143 func.typf = (typf > 0.0) ? typf : 1.0;
144 func.h = (h > 0.0) ? h : pow(macheps(), H_POWER);
145 for (i = 0; i < n; i++)
146 func.typx[i] = (typx != nil && typx[i] > 0.0) ? typx[i] : 1.0;
147 func.dflt = dflt;
148 func.fderivs = 0;
151 /* install tilt functions */
152 static void
153 install_gfuncs(char **g, int n, int k, int change_sign, double h, double *typx)
155 int i;
156 static int inited = FALSE;
157 static double *gfsumdata = nil;
159 if (! inited) {
160 gfuncs.typx = nil;
161 gfuncs.gfsum = nil;
162 gfuncs.gderivs = nil;
163 inited = TRUE;
165 makespace(&gfuncs.typx, n * sizeof(double));
166 makespace(&gfuncs.gfsum, k * sizeof(double *));
167 makespace(&gfsumdata, k * n * sizeof(double));
168 makespace(&gfuncs.gderivs, k *sizeof(int));
170 gfuncs.g = g;
171 gfuncs.n = n;
172 gfuncs.k = k;
173 gfuncs.change_sign = change_sign;
174 gfuncs.h = (h > 0.0) ? h : pow(macheps(), H_POWER);
175 for (i = 0; i < n; i++)
176 gfuncs.typx[i] = (typx != nil && typx[i] > 0.0) ? typx[i] : 1.0;
177 for (i = 0; i < k; i++) gfuncs.gfsum[i] = gfsumdata + i * n;
178 return;
181 /* install constraint functions */
182 static void
183 install_cfuncs(char **g, int n, int k, double *ctarget, double h, double *typx)
185 int i;
186 static int inited = FALSE;
188 if (! inited) {
189 cfuncs.typx = nil;
190 cfuncs.fsum = nil;
191 cfuncs.gderivs = nil;
192 inited = TRUE;
194 makespace(&cfuncs.typx, n * sizeof(double));
195 makespace(&cfuncs.fsum, n * sizeof(double));
196 makespace(&cfuncs.gderivs, k * sizeof(int));
198 cfuncs.g = g;
199 cfuncs.n = n;
200 cfuncs.k = k;
201 cfuncs.h = (h > 0.0) ? h : pow(macheps(), H_POWER);
202 for (i = 0; i < n; i++)
203 cfuncs.typx[i] = (typx != nil && typx[i] > 0.0) ? typx[i] : 1.0;
204 cfuncs.ctarget = ctarget;
205 return;
208 static int
209 in_support(char **ff, int n, double *x)
211 char *args[1], *values[1];
212 int *result;
213 char *mode[1];
214 long length[1];
216 if (ff == nil || ff[0] == nil) {
217 return(TRUE);
218 } else {
219 mode[0] = "double";
220 length[0] =n;
221 args[0] = (char *) x;
222 call_S(ff[0], 1L, args, mode, length, 0L, 1L, values);
223 result = (int *) values[0];
224 return(result[0]);
228 /* callback for logposterior evaluation */
229 static int
230 evalfunc(double *x, double *pval, double *grad, double **hess)
232 char *args[1], *values[3];
233 double *result, val;
234 char *mode[1];
235 long length[1];
236 int i, j;
238 for (i = 0; i < 3; i++) values[i] = nil;
240 if (in_support(func.sf, func.n, x)) {
241 if (pval != nil || func.fderivs > 0 || hess != nil) {
242 mode[0] = "double";
243 length[0] = func.n;
244 args[0] = (char *) x;
245 call_S(func.f, 1L, args, mode, length, 0L, 3L, values);
246 result = (double *) values[0];
247 val = (! func.change_sign) ? result[0] : -result[0];
248 if (pval != nil) *pval = val;
249 if (values[2] != nil) func.fderivs = 2;
250 else if (values[1] != nil) func.fderivs = 1;
251 else func.fderivs = 0;
253 if (grad != nil) {
254 if (func.fderivs > 0) {
255 result = (double *) values[1];
256 for (i = 0; i < func.n; i++)
257 grad[i] = (! func.change_sign) ? result[i] : -result[i];
259 else {
260 numergrad(func.n, x, grad, func.fsum, evalfunc, func.h, func.typx);
263 if (hess != nil) {
264 if (func.fderivs == 2) {
265 result = (double *) values[2];
266 for (i = 0; i < func.n; i++)
267 for (j = 0; j < func.n; j++)
268 hess[i][j] = (! func.change_sign) ? result[i + j * func.n]
269 : -result[i + j * func.n];
271 else {
272 if (func.fderivs == 1) /* kludge to get fsum for analytic gradients */
273 numergrad(func.n, x, func.fsum, func.fsum,
274 evalfunc, func.h, func.typx);
275 numerhess(func.n, x, hess, val, func.fsum, evalfunc, func.h, func.typx);
278 return(TRUE);
279 } else {
280 if (pval != nil) {
281 *pval = func.dflt;
283 return(FALSE);
285 return(TRUE);
289 /* callback for tilt function evaluation */
290 static int which_gfunc;
292 static int
293 evalgfunc(double *x, double *pval, double *grad, double **hess)
295 char *args[1], *values[3];
296 double *result, val;
297 char *mode[1];
298 long length[1];
299 int i, j;
301 for (i = 0; i < 3; i++) values[i] = nil;
303 if (pval != nil || gfuncs.gderivs[which_gfunc] > 0 || hess != nil) {
304 mode[0] = "double";
305 length[0] = gfuncs.n;
306 args[0] = (char *) x;
307 call_S(gfuncs.g[which_gfunc], 1L, args, mode, length, 0L, 3L, values);
308 result = (double *) values[0];
309 val = result[0];
310 if (pval != nil) *pval = result[0];
311 if (values[2] != nil) gfuncs.gderivs[which_gfunc] = 2;
312 else if (values[1] != nil) gfuncs.gderivs[which_gfunc] = 1;
313 else gfuncs.gderivs[which_gfunc] = 0;
315 if (grad != nil) {
316 if (gfuncs.gderivs[which_gfunc] > 0) {
317 result = (double *) values[1];
318 for (i = 0; i < gfuncs.n; i++) grad[i] = result[i];
320 else {
321 numergrad(gfuncs.n, x, grad, gfuncs.gfsum[which_gfunc], evalgfunc,
322 gfuncs.h, gfuncs.typx);
325 if (hess != nil) {
326 if (gfuncs.gderivs[which_gfunc] == 2) {
327 result = (double *) values[2];
328 for (i = 0; i < gfuncs.n; i++)
329 for (j = 0; j < gfuncs.n; j++)
330 hess[i][j] = result[i + j * gfuncs.n];
332 else {
333 /* kludge to get fsum if analytic gradient used */
334 if (gfuncs.gderivs[which_gfunc] == 1)
335 numergrad(gfuncs.n, x, gfuncs.gfsum[which_gfunc],
336 gfuncs.gfsum[which_gfunc], evalgfunc, gfuncs.h, gfuncs.typx);
337 numerhess(gfuncs.n, x, hess, val, gfuncs.gfsum[which_gfunc], evalgfunc,
338 gfuncs.h, gfuncs.typx);
341 return(TRUE);
344 /* callback for constraint function evaluation */
345 static int which_cfunc;
347 static int
348 evalcfunc(double *x, double *pval, double *grad, double **hess)
350 char *args[1], *values[3];
351 double *result, val;
352 char *mode[1];
353 long length[1];
354 int i, j;
356 if (pval != nil || cfuncs.gderivs[which_cfunc] > 0 || hess != nil) {
357 mode[0] = "double";
358 length[0] = cfuncs.n;
359 args[0] = (char *) x;
360 call_S(cfuncs.g[which_cfunc], 1L, args, mode, length, 0L, 3L, values);
361 result = (double *) values[0];
362 val = result[0];
363 if (pval != nil) {
364 *pval = result[0];
365 if (cfuncs.ctarget != nil) *pval -= cfuncs.ctarget[which_cfunc];
367 if (values[2] != nil) cfuncs.gderivs[which_cfunc] = 2;
368 else if (values[1] != nil) cfuncs.gderivs[which_cfunc] = 1;
369 else cfuncs.gderivs[which_cfunc] = 0;
371 if (grad != nil) {
372 if (cfuncs.gderivs[which_cfunc] > 0) {
373 result = (double *) values[1];
374 for (i = 0; i <cfuncs.n; i++) grad[i] = result[i];
376 else {
377 numergrad(cfuncs.n, x, grad, cfuncs.fsum, evalcfunc,
378 cfuncs.h, cfuncs.typx);
381 if (hess != nil) {
382 if (cfuncs.gderivs[which_cfunc] == 2) {
383 result = (double *) values[2];
384 for (i = 0; i <cfuncs.n; i++)
385 for (j = 0; j <cfuncs.n; j++)
386 hess[i][j] = result[i + j * cfuncs.n];
388 else {
389 /* kludge to get fsum if analytic gradient used */
390 if (cfuncs.gderivs[which_cfunc] == 1)
391 numergrad(cfuncs.n, x, cfuncs.fsum, cfuncs.fsum, evalcfunc,
392 cfuncs.h, cfuncs.typx);
393 numerhess(cfuncs.n, x, hess, val, cfuncs.fsum, evalcfunc,
394 cfuncs.h, cfuncs.typx);
397 return(TRUE);
400 /* S front end for logposterior evaluation */
401 void
402 evalfront(char **ff, int *n, double *x, double *val, double *grad,
403 double *phess, double *h, double *typx)
405 int i;
406 static double **hess = nil;
408 install_func(ff[0], nil, *n, FALSE, 1.0, *h, typx, 0.0);
409 if (phess == nil) hess = nil;
410 else {
411 makespace(&hess, *n * sizeof(double *));
412 for (i = 0; i < *n; i++, phess += *n) hess[i] = phess;
414 evalfunc(x, val, grad, hess);
417 /* S front end for tilt function evaluation */
418 void
419 gevalfront(char **gg, int *n, int *m, double *x, double *h,
420 double *typx, double *val, double *grad)
422 int i;
424 install_gfuncs(gg, *n, *m, FALSE, *h, typx);
425 for (i = 0; i < *m; i++, val++) {
426 which_gfunc = i;
427 evalgfunc(x, val, grad, nil);
428 if (grad != nil) grad += *n;
430 return;
433 /************************************************************************/
434 /** **/
435 /** Derivative Scaling Routines **/
436 /** **/
437 /************************************************************************/
440 void
441 derivscalefront(char **ff, int *n, double *x, double *h, double *typx, double *tol, int *info)
443 int i;
445 if (*tol <= 0.0) *tol = pow(macheps(), GRADTOL_POWER);
447 install_func(ff[0], nil, *n, TRUE, 1.0, *h, typx, 0.0);
448 *info = check_derivs(x, *tol);
450 *h = func.h;
451 for (i = 0; i < *n; i++) typx[i] = func.typx[i];
452 return;
456 static int
457 check_derivs(double *x, double drvtol)
459 static double *grad = nil, work = nil;
460 static double **hess = nil;
461 int i, error;
463 grad = (double *) S_alloc(func.n, sizeof(double));
464 hess = (double **) S_alloc(func.n, sizeof(double *));
465 work = (double *) S_alloc(func.n + func.n * func.n, sizeof(double));
467 for (i = 0; i < func.n; i++) {
468 hess[i] = work;
469 work += func.n;
472 error = derivscale(func.n, x, grad, hess, func.fsum, evalfunc,
473 &func.h, func.typx, drvtol, work);
474 return(error);
478 /************************************************************************/
479 /** **/
480 /** Importance Sampling Routines **/
481 /** **/
482 /************************************************************************/
484 /* joint density of normal-cauchy mixture */
485 static double
486 dncmix(double *x, int n, double p)
488 int i;
489 double dens;
491 for (i = 0, dens = 1.0; i < n; i++) {
492 dens *= p * exp(-0.5 * x[i] * x[i]) / ROOT2PI
493 + (1 - p) * PI_INV / (1.0 + x[i] * x[i]);
495 return(dens);
499 * S front end for computing sample from transformed normal-cauchy
500 * mixture and importance sampling weights
502 void
503 samplefront(char **ff, char **sf, char **rf,
504 double *p, int *n,
505 double *x, double *ch, int *N, double *y, double *w)
507 double val;
508 int i, j, k;
509 char *args[2], *values[1];
510 double *result, mval, c, dens, m;
511 char *mode[2];
512 long length[2];
514 /* get the random variables */
515 mode[0] = "double"; mode[1] = "double";
516 length[0] = 1; length[1] = 1;
517 m = *N * *n; args[0] = (char *) &m; args[1] = (char *) p;
518 call_S(rf[0], 2L, args, mode, length, 0L, 1L, values);
519 result = (double *) values[0];
520 for (i = 0; i < m; i++) y[i] = result[i];
522 /* construct the sample and the weights */
523 install_func(ff[0], sf, *n, FALSE, 1.0, -1.0, nil, 0.0);
524 c = 1.0 / pow(ROOT2PI, (double) *n);
525 evalfunc(x, &mval, nil, nil);
526 for (i = 0; i < *N; i++, y += *n) {
527 dens = dncmix(y, *n, *p);
528 for (j = 0; j < *n; j++) {
529 val = x[j];
530 for (k = j; k < *n; k++) val += y[k] * ch[j + *n * k];
531 y[j] = val;
533 if (evalfunc(y, &val, nil, nil)) w[i] = exp(val - mval) * c / dens;
534 else w[i] = 0.0;
536 return;
540 /************************************************************************/
541 /** **/
542 /** Maximization Routines **/
543 /** **/
544 /************************************************************************/
546 typedef struct {
547 int n, m, k, itnlimit, backtrack, verbose, vals_suppl, exptilt;
548 int count, termcode;
549 } MaxIPars;
551 typedef struct {
552 double typf, h, gradtol, steptol, maxstep, dflt, tilt, newtilt, hessadd;
553 } MaxDPars;
555 struct {
556 double tilt;
557 double *gval;
558 double **ggrad, **ghess;
559 int exptilt;
560 double *tscale;
561 } tiltinfo;
563 static void
564 add_tilt(double *x, double *pval, double *grad, double **hess,
565 double tilt, int exptilt)
567 int i, j, k, n = func.n, m = gfuncs.k;
568 double *gval, *ggrad, **ghess, etilt;
570 if (m == 0) return;
572 if (gfuncs.change_sign) tilt = -tilt;
574 for (k = 0; k < m; k++) {
575 gval = (pval != nil) ? tiltinfo.gval + k : nil;
576 ggrad = (grad != nil) ? tiltinfo.ggrad[k] : nil;
577 ghess = (hess != nil) ? tiltinfo.ghess : nil;
579 which_gfunc = k;
580 evalgfunc(x, gval, ggrad, ghess);
582 if (exptilt) {
583 etilt = (tiltinfo.tscale != nil) ? tilt / tiltinfo.tscale[k] : tilt;
584 if (pval != nil) *pval += etilt * *gval;
585 if (grad != nil)
586 for (i = 0; i < n; i++) grad[i] += etilt * ggrad[i];
587 if (hess != nil)
588 for (i = 0; i < n; i++)
589 for (j = 0; j < n; j++) hess[i][j] += etilt * ghess[i][j];
591 else {
592 gval = tiltinfo.gval;
593 ggrad = tiltinfo.ggrad[k];
594 ghess = tiltinfo.ghess;
595 if (gval[k] <= 0.0) Recover("nonpositive function value", NULL);
596 if (pval != nil) *pval += tilt * log(gval[k]);
597 if (grad != nil)
598 for (i = 0; i < n; i++) grad[i] += tilt * ggrad[i] / gval[k];
599 if (hess != nil)
600 for (i = 0; i < n; i++)
601 for (j = 0; j < n; j++)
602 hess[i][j] +=
603 tilt * (ghess[i][j] / gval[k]
604 - (ggrad[i] / gval[k]) * (ggrad[j] / gval[k]));
609 static void
610 set_tilt_info(int n, int m,
611 double tilt, int exptilt, double *tscale)
613 static double *hessdata = nil, *graddata = nil;
614 int i;
615 static int inited = FALSE;
617 if (! inited) {
618 tiltinfo.gval = nil;
619 tiltinfo.ggrad = nil;
620 tiltinfo.ghess = nil;
621 inited = TRUE;
623 makespace(&tiltinfo.gval, n * sizeof(double));
624 makespace(&tiltinfo.ggrad, m * sizeof(double *));
625 makespace(&tiltinfo.ghess, n * sizeof(double *));
626 makespace(&graddata, n * m * sizeof(double));
627 makespace(&hessdata, n * n * sizeof(double));
629 tiltinfo.tilt = tilt;
630 tiltinfo.exptilt = exptilt;
631 for (i = 0; i < m; i++) tiltinfo.ggrad[i] = graddata + i * n;
632 for (i = 0; i < n; i++) tiltinfo.ghess[i] = hessdata + i * n;
633 tiltinfo.tscale = tscale;
637 static void
638 minfunc(double *x, double *pval, double *grad, double **hess)
640 int k = gfuncs.k;
642 if (evalfunc(x, pval, grad, hess) && (k > 0))
643 add_tilt(x, pval, grad, hess, tiltinfo.tilt, tiltinfo.exptilt);
646 void
647 constfunc(double *x, double *vals, double **jac, double **hess)
649 int i, k = cfuncs.k;
650 double *pvali, *jaci;
652 for (i = 0; i < k; i++) {
653 pvali = (vals != nil) ? vals + i : nil;
654 jaci = (jac != nil) ? jac[i] : nil;
655 which_cfunc = i;
656 evalcfunc(x, pvali, jaci, nil);
660 void
661 maxfront(char **ff, char **gf, char **cf,
662 double *x, double *typx, double *fvals, double *gvals, double *cvals, double *ctarget,
663 MaxIPars *ipars, MaxDPars *dpars,
664 double *tscale, char **msg)
666 static char *work = nil;
667 static double **H = nil, **cJ = nil;
668 double *pf, *grad, *c;
669 int i, n, m, k;
670 void (*cfun)();
672 if (ipars->verbose > 0) PRINTSTR("maximizing...\n");
674 n = ipars->n;
675 m = ipars->m;
676 k = ipars->k;
677 if (k >= n) Recover("too many constraints", NULL);
679 makespace(&H, n * sizeof(double *));
680 makespace(&work, minworkspacesize(n, k));
682 pf = fvals; fvals++;
683 grad = fvals; fvals += n;
684 for (i = 0; i < n; i++, fvals += n) H[i] = fvals;
685 set_tilt_info(n, m, dpars->newtilt, ipars->exptilt, tscale);
687 if (k == 0) {
688 c = nil;
689 cJ = nil;
690 cfun = nil;
691 } else {
692 c = cvals;
693 cvals += k;
694 makespace(&cJ, k * sizeof(double *));
695 for (i = 0; i < k; i++) cJ[i] = cvals + i * n;
696 cfun = &constfunc; /* AJR: pointer to constfunc? */
699 install_func(ff[0], nil, n, TRUE, dpars->typf, dpars->h, typx, dpars->dflt);
700 install_gfuncs(gf, n, m, TRUE, dpars->h, typx);
701 install_cfuncs(cf, n, k, ctarget, dpars->h, typx);
703 minsetup(n, k, &minfunc, &cfun, x, dpars->typf, typx, work); /* AJR: FIXME arg 3,4 */
704 minsetoptions(dpars->gradtol, dpars->steptol, dpars->maxstep,
705 ipars->itnlimit, ipars->verbose, ipars->backtrack, TRUE);
707 if (ipars->vals_suppl) {
708 for (i = 0; i < k; i++) c[i] -= ctarget[i];
709 if (dpars->newtilt != dpars->tilt) {
710 add_tilt(x, pf, grad, H, dpars->newtilt - dpars->tilt, ipars->exptilt);
711 dpars->tilt = dpars->newtilt;
713 minsupplyvalues(*pf, grad, H, c, cJ);
716 minimize();
717 minresults(x, pf, nil, grad, H, c, cJ, &ipars->count, &ipars->termcode,
718 &dpars->hessadd);
719 msg[0] = minresultstring(ipars->termcode);
721 for (i = 0; i < k; i++) c[i] += ctarget[i];
722 ipars->vals_suppl = TRUE;
725 /************************************************************************/
726 /** **/
727 /** Log Laplace Approximation **/
728 /** **/
729 /************************************************************************/
731 void
732 loglapdet(double *fvals, double *cvals, MaxIPars *ipars, MaxDPars *dpars,
733 double *val, int *detonly)
735 int i, j, l, n = ipars->n, k = ipars->k;
736 double f = -fvals[0], *hessdata = fvals + n + 1, *cgraddata = cvals + k;
737 double ldL, ldcv, maxadd;
738 static double **hess = nil, **cgrad = nil;
740 if (k >= n) Recover("too many constraints", NULL);
742 makespace(&hess, n * sizeof(double *));
743 makespace(&cgrad, k * sizeof(double *));
745 for (i = 0; i < n; i++) hess[i] = hessdata + i * n;
746 for (i = 0; i < k; i++) cgrad[i] = cgraddata + i * n;
748 choldecomp(hess, n, 0.0, &maxadd);
749 /**** do something if not pos. definite ****/
751 for (i = 0, ldL = 0.0; i < n; i++) ldL += log(hess[i][i]);
753 if (k > 0) {
754 /* forward solve for (L^-1) cgrad^T */
755 for (l = 0; l < k; l++) {
756 for (i = 0; i < n; i++) {
757 if (hess[i][i] != 0.0) cgrad[l][i] /= hess[i][i];
758 for (j = i + 1; j < n; j++) cgrad[l][j] -= hess[j][i] * cgrad[l][i];
762 /* compute sigma and stdev */
763 for (i = 0; i < k; i++) {
764 for (j = i; j < k; j++) {
765 for (l = 0, hess[i][j] = 0.0; l < n; l++)
766 hess[i][j] += cgrad[i][l] * cgrad[j][l];
767 hess[j][i] = hess[i][j];
771 choldecomp(hess, k, 0.0, &maxadd);
772 /**** do something if not pos. definite ****/
773 for (i = 0, ldcv = 0.0; i < k; i++) ldcv += log(hess[i][i]);
775 else ldcv = 0.0;
777 *val = (n - k) * log(ROOT2PI) - ldL - ldcv;
778 if (! *detonly) *val += f;
781 #ifdef SBAYES
783 loglapfront(fvals, cvals, ipars, dpars, val)
784 double *fvals, *cvals;
785 MaxIPars *ipars;
786 MaxDPars *dpars;
787 double *val;
789 int detonly = FALSE;
791 loglapdet(fvals, cvals, ipars, dpars, val, &detonly);
794 /************************************************************************/
795 /** **/
796 /** First Order Moments **/
797 /** **/
798 /************************************************************************/
800 moms1front(gf, n, m, x, hessdata, mean, stdev, sigmadata, h, typx)
801 char *gf;
802 int *n, *m;
803 double *x, *hessdata, *mean, *stdev, *sigmadata, *h, *typx;
805 int i, j, k;
806 double *hess, *sigma, *delg;
807 double *delgdata, maxadd;
809 hess = (double **) S_alloc(*n, sizeof(double *));
810 sigma = (double **) S_alloc(*m, sizeof(double *));
811 delg = (double **) S_alloc(*m, sizeof(double *));
812 delgdata = (double *) S_alloc(*m * *n, sizeof(double));
814 for (i = 0; i < *n; i++) hess[i] = hessdata + i * *n;
815 for (i = 0; i < *m; i++) sigma[i] = sigmadata + i * *m;
816 for (i = 0; i < *m; i++) delg[i] = delgdata + i * *n;
818 gevalfront(gf, n, m, x, h, typx, mean, delgdata);
820 /* get the cholesky decomposition L of the hessian */
821 choldecomp(hess, *n, 0.0, &maxadd);
823 /* forward solve for (L^-1) delg^T */
824 for (k = 0; k < *m; k++) {
825 for (i = 0; i < *n; i++) {
826 if (hess[i][i] != 0.0) delg[k][i] /= hess[i][i];
827 for (j = i + 1; j < *n; j++) delg[k][j] -= hess[j][i] * delg[k][i];
831 /* compute sigma and stdev */
832 for (i = 0; i < *m; i++) {
833 for (j = i; j < *m; j++) {
834 for (k = 0, sigma[i][j] = 0.0; k < *n; k++)
835 sigma[i][j] += delg[i][k] * delg[j][k];
836 sigma[j][i] = sigma[i][j];
838 stdev[i] = sqrt(sigma[i][i]);
842 /************************************************************************/
843 /** **/
844 /** Second Order Moments **/
845 /** **/
846 /************************************************************************/
848 typedef struct {
849 MaxIPars max;
850 int full, covar;
851 } MomIPars;
853 typedef struct {
854 MaxDPars max;
855 double mgfdel;
856 } MomDPars;
858 moms2front(ff, gf, gnum, x, typx, fvals, gvals, ipars, dpars,
859 mean, stdev, sigmadata)
860 char **ff, **gf;
861 int *gnum;
862 double *x, *typx, *fvals, *gvals, *mean, *stdev, *sigmadata;
863 MomIPars *ipars;
864 MomDPars *dpars;
866 char *msg;
867 double h, loglap0, loglap1, loglap2;
868 double *tilts, *fvals1, *gvals1, *x1;
869 MomDPars dp1, *dpars1 = &dp1;
870 MomIPars ip1, *ipars1 = &ip1;
871 int i, n, m;
873 n = ipars->max.n;
874 m = *gnum;
875 h = dpars->max.h;
877 tilts = (double *) S_alloc(m, sizeof(double));
878 x1 = (double *) S_alloc(n, sizeof(double));
879 fvals1 = (double *) S_alloc(1 + n + n * n, sizeof(double));
880 gvals1 = (double *) S_alloc(m + n * m, sizeof(double));
882 maxfront(ff, nil, nil, x, typx, fvals, gvals, nil, nil,
883 ipars, dpars, nil, &msg);
884 copypars(x, fvals, gvals, ipars, dpars, x1, fvals1, gvals1, ipars1, dpars1);
885 loglapfront(fvals1, nil, ipars1, dpars1, &loglap0);
886 copypars(x, fvals, gvals, ipars, dpars, x1, fvals1, gvals1, ipars1, dpars1);
887 moms1front(gf, &n, &m, x1, fvals1 + n + 1, mean, stdev, sigmadata, &h, typx);
889 if (ipars->full) {
890 for (i = 0; i < m; i++) {
891 copypars(x, fvals, gvals, ipars, dpars,
892 x1, fvals1, gvals1, ipars1, dpars1);
893 ipars1->max.m = 1;
894 ipars1->max.exptilt = FALSE;
895 dpars1->max.newtilt = 1.0;
896 maxfront(ff, gf + i, nil, x1, typx, fvals1, gvals1, nil, nil,
897 ipars1, dpars1, nil, &msg);
898 loglapfront(fvals1, nil, ipars1, dpars1, &loglap1);
899 loglap1 -= loglap0;
901 copypars(x, fvals, gvals, ipars, dpars,
902 x1, fvals1, gvals1, ipars1, dpars1);
903 ipars1->max.m = 1;
904 ipars1->max.exptilt = FALSE;
905 dpars1->max.newtilt = 2.0;
906 maxfront(ff, gf + i, nil, x1, typx, fvals1, gvals1, nil, nil,
907 ipars1, dpars1, nil, &msg);
908 loglapfront(fvals1, nil, ipars1, dpars1, &loglap2);
909 loglap2 -= loglap0;
911 mean[i] = exp(loglap1);
912 stdev[i] = sqrt(exp(loglap2) - exp(2.0 * loglap1));
913 if (ipars->covar) sigmadata[i + i * m] = stdev[i] * stdev[i];
915 if (ipars->covar) {
916 char *cgf[2];
917 int j;
919 for (i = 0; i < m; i++) {
920 for (j = i + 1; j < m; j++) {
921 cgf[0] = gf[i];
922 cgf[1] = gf[j];
923 copypars(x, fvals, gvals, ipars, dpars,
924 x1, fvals1, gvals1, ipars1, dpars1);
925 ipars1->max.m = 2;
926 ipars1->max.exptilt = FALSE;
927 dpars1->max.newtilt = 1.0;
928 maxfront(ff, gf + i, nil, x1, typx, fvals1, gvals1, nil, nil,
929 ipars1, dpars1, nil, &msg);
930 loglapfront(fvals1, nil, ipars1, dpars1, &loglap1);
931 loglap1 -= loglap0;
933 sigmadata[i + j * m] = exp(loglap1) - mean[i] * mean[j];
934 sigmadata[j + i * m] = sigmadata[i + j * m];
939 else {
940 for (i = 0; i < m; i++)
941 tilts[i] = (stdev[i] > 0.0) ? dpars->mgfdel / stdev[i] : dpars->mgfdel;
943 for (i = 0; i < m; i++) {
944 copypars(x, fvals, gvals, ipars, dpars,
945 x1, fvals1, gvals1, ipars1, dpars1);
946 ipars1->max.m = 1;
947 ipars1->max.exptilt = TRUE;
948 dpars1->max.newtilt = tilts[i];
949 maxfront(ff, gf + i, nil, x1, typx, fvals1, gvals1, nil, nil,
950 ipars1, dpars1, nil, &msg);
951 loglapfront(fvals1, nil, ipars1, dpars1, &loglap1);
952 loglap1 -= loglap0;
954 copypars(x, fvals, gvals, ipars, dpars,
955 x1, fvals1, gvals1, ipars1, dpars1);
956 ipars1->max.m = 1;
957 ipars1->max.exptilt = TRUE;
958 dpars1->max.newtilt = -tilts[i];
959 maxfront(ff, gf + i, nil, x1, typx, fvals1, gvals1, nil, nil,
960 ipars1, dpars1, nil, &msg);
961 loglapfront(fvals1, nil, ipars1, dpars1, &loglap2);
962 loglap2 -= loglap0;
964 mean[i] = (loglap1 - loglap2) / (2.0 * tilts[i]);
965 stdev[i] = sqrt((loglap1 + loglap2) / (tilts[i] * tilts[i]));
966 if (ipars->covar) sigmadata[i + i * m] = stdev[i] * stdev[i];
968 if (ipars->covar) {
969 char *cgf[2];
970 double ctilt, tscale[2];
971 int j;
973 ctilt = dpars->mgfdel;
974 for (i = 0; i < m; i++) {
975 for (j = i + 1; j < m; j++) {
976 cgf[0] = gf[i];
977 cgf[1] = gf[j];
978 tscale[0] = stdev[i] > 0 ? stdev[i] : 1.0;
979 tscale[1] = stdev[j] > 0 ? stdev[j] : 1.0;
981 copypars(x, fvals, gvals, ipars, dpars,
982 x1, fvals1, gvals1, ipars1, dpars1);
983 ipars1->max.m = 2;
984 ipars1->max.exptilt = TRUE;
985 dpars1->max.newtilt = ctilt;
986 maxfront(ff, cgf, nil, x1, typx, fvals1, gvals1, nil, nil,
987 ipars1, dpars1, tscale, &msg);
988 loglapfront(fvals1, nil, ipars1, dpars1, &loglap1);
989 loglap1 -= loglap0;
991 copypars(x, fvals, gvals, ipars, dpars,
992 x1, fvals1, gvals1, ipars1, dpars1);
993 ipars1->max.m = 2;
994 ipars1->max.exptilt = TRUE;
995 dpars1->max.newtilt = -ctilt;
996 maxfront(ff, cgf, nil, x1, typx, fvals1, gvals1, nil, nil,
997 ipars1, dpars1, tscale, &msg);
998 loglapfront(fvals1, nil, ipars1, dpars1, &loglap2);
999 loglap2 -= loglap0;
1001 sigmadata[i + j * m] = stdev[i] * stdev[j]
1002 * ((loglap2 + loglap1) /(2 * ctilt * ctilt) - 1.0);
1003 sigmadata[j + i * m] = sigmadata[i + j * m];
1010 static copypars(x, f, g, ip, dp, x1, f1, g1, ip1, dp1)
1011 double *x, *f, *g, *x1, *f1, *g1;
1012 MomIPars *ip, *ip1;
1013 MomDPars *dp, *dp1;
1015 int i, n, m, nf, ng;
1017 n = ip->max.n;
1018 m = ip->max.m;
1019 nf = 1 + n + n * n;
1020 ng = m + n * m;
1022 for (i = 0; i < n; i++) x1[i] = x[i];
1023 for (i = 0; i < nf; i++) f1[i] = f[i];
1024 for (i = 0; i < ng; i++) g1[i] = g[i];
1025 *ip1 = *ip;
1026 *dp1 = *dp;
1029 /************************************************************************/
1030 /** **/
1031 /** Laplace Margins **/
1032 /** **/
1033 /************************************************************************/
1035 lapmar1front(ff, gf, x, typx, fvals, ipars, dpars, xmar, ymar, nmar)
1036 char **ff, **gf;
1037 double *x, *typx, *fvals, *xmar, *ymar;
1038 MaxIPars *ipars;
1039 MaxDPars *dpars;
1040 int *nmar;
1042 char *msg;
1043 int i, n, m, mindex;
1044 double h, loglap0, loglap1, xmode, stdev, sigmadata, ctarget[1];
1045 double *fvals1, *x1, *cvals, *cvals1, *fvals2, *x2, *cvals2;
1046 MaxDPars dp1, dp2, *dpars1 = &dp1, *dpars2 = &dp2;
1047 MaxIPars ip1, ip2, *ipars1 = &ip1, *ipars2 = &ip2;
1049 n = ipars->n;
1050 m = 1;
1051 h = dpars->h;
1053 x1 = (double *) S_alloc(n + 1, sizeof(double));
1054 fvals1 = (double *) S_alloc(1 + n + n * n, sizeof(double));
1055 cvals = (double *) S_alloc(m + n * m, sizeof(double));
1056 cvals1 = (double *) S_alloc(m + n * m, sizeof(double));
1057 x2 = (double *) S_alloc(n + 1, sizeof(double));
1058 fvals2 = (double *) S_alloc(1 + n + n * n, sizeof(double));
1059 cvals2 = (double *) S_alloc(m + n * m, sizeof(double));
1061 maxfront(ff, nil, nil, x, typx, fvals, nil, nil, nil,
1062 ipars, dpars, nil, &msg);
1063 cpmarpars(x, fvals, cvals, ipars, dpars, x1, fvals1, cvals1, ipars1, dpars1);
1064 loglapfront(fvals1, nil, ipars1, dpars1, &loglap0);
1065 cpmarpars(x, fvals, cvals, ipars, dpars, x1, fvals1, cvals1, ipars1, dpars1);
1066 moms1front(gf, &n, &m, x1, fvals1 + n + 1, &xmode, &stdev, &sigmadata,
1067 &h, typx);
1069 ipars->k = 1;
1070 ipars->vals_suppl = FALSE;
1071 ctarget[0] = xmode;
1072 maxfront(ff, nil, gf, x, typx, fvals, nil, cvals, ctarget,
1073 ipars, dpars, nil, &msg);
1075 for (mindex = 0; mindex < *nmar && xmar[mindex] <= xmode; mindex++);
1077 cpmarpars(x, fvals, cvals, ipars, dpars, x1, fvals1, cvals1, ipars1, dpars1);
1078 for (i = mindex; i >= 0; i--) {
1079 ctarget[0] = xmar[i];
1080 maxfront(ff, nil, gf, x1, typx, fvals1, nil, cvals1, ctarget,
1081 ipars1, dpars1, nil, &msg);
1082 cpmarpars(x1, fvals1, cvals1, ipars1, dpars1, x2,
1083 fvals2, cvals2, ipars2, dpars2);
1084 loglapfront(fvals2, cvals2, ipars2, dpars2, &loglap1);
1085 ymar[i] = exp(loglap1 - loglap0);
1087 cpmarpars(x, fvals, cvals, ipars, dpars, x1, fvals1, cvals1, ipars1, dpars1);
1088 for (i = mindex + 1; i < *nmar; i++) {
1089 ctarget[0] = xmar[i];
1090 maxfront(ff, nil, gf, x1, typx, fvals1, nil, cvals1, ctarget,
1091 ipars1, dpars1, nil, &msg);
1092 cpmarpars(x1, fvals1, cvals1, ipars1, dpars1, x2,
1093 fvals2, cvals2, ipars2, dpars2);
1094 loglapfront(fvals2, cvals2, ipars2, dpars2, &loglap1);
1095 ymar[i] = exp(loglap1 - loglap0);
1099 static cpmarpars(x, f, g, ip, dp, x1, f1, g1, ip1, dp1)
1100 double *x, *f, *g, *x1, *f1, *g1;
1101 MaxIPars *ip, *ip1;
1102 MaxDPars *dp, *dp1;
1104 int i, n, k, nf, ng;
1106 n = ip->n;
1107 k = ip->k;
1108 nf = 1 + n + n * n;
1109 ng = k + n * k;
1111 for (i = 0; i < n; i++) x1[i] = x[i];
1112 for (i = 0; i < nf; i++) f1[i] = f[i];
1113 for (i = 0; i < ng; i++) g1[i] = g[i];
1114 *ip1 = *ip;
1115 *dp1 = *dp;
1117 #endif /* SBAYES */
1120 TODO
1122 get hessian from gradiant for analytical gradiants
1124 avoid repeated derivative calls in mimimize.
1126 2d margins
1128 use pos. definiteness info in margins