lib/functions.c

   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.                       */
   6
   7 #include <stdio.h>
   8 #include "xmath.h"
   9
  10 #define PRINTSTR(s) printf(s)
  11
  12 extern void *S_alloc(), *calloc(), *realloc();
  13 extern double macheps();
  14 char *minresultstring();
  15
  16 extern void choldecomp();
  17
  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);
  28
  29 extern void choldecomp();
  30
  31
  32
  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);
  37
  38 /* next 2 from derivatives.c */
  39 /*
  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);
  45 */
  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 *);
  53
  54 /* next from minimize */
  55 extern int minworkspacesize(int, int);
  56
  57 /************************************************************************/
  58 /**                                                                    **/
  59 /**                      Definitions and Globals                       **/
  60 /**                                                                    **/
  61 /************************************************************************/
  62
  63 /* #define NULL 0L  already defined in stddef */
  64
  65 #define nil 0L
  66 #define TRUE 1
  67 #define FALSE 0
  68
  69 #define ROOT2PI 2.50662827463100050241
  70 #define PI_INV .31830988618379067153
  71
  72 #define GRADTOL_POWER 1.0 / 3.0
  73 #define H_POWER 1.0 / 6.0
  74
  75 /*typedef double **RMatrix, *RVector;*/
  76
  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;
  86
  87 static Fundata func, gfuncs, cfuncs;
  88
  89 /************************************************************************/
  90 /**                                                                    **/
  91 /**                         Memory Utilities                           **/
  92 /**                                                                    **/
  93 /************************************************************************/
  94
  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.         */
 100
 101 void static
 102 makespace(pptr, size)
 103      char **pptr;
 104      int size;
 105 {
 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);
 110 }
 111
 112 /************************************************************************/
 113 /**                                                                    **/
 114 /**                    Functions Evaluation Routines                   **/
 115 /**                                                                    **/
 116 /************************************************************************/
 117
 118 /*
 119  * All Hessian evaluations by numerical derivatives assume the gradient is
 120  * evaluated first at the same location. The results are cached away.
 121  */
 122
 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)
 127 {
 128   int i;
 129   static int inited = FALSE;
 130
 131   if (! inited) {
 132     func.typx = nil;
 133     func.fsum = nil;
 134     inited = TRUE;
 135   }
 136   makespace(&func.typx, n * sizeof(double));
 137   makespace(&func.fsum, n * sizeof(double));
 138
 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;
 149 }
 150
 151 /* install tilt functions */
 152 static void
 153 install_gfuncs(char **g, int n, int k, int change_sign, double h, double *typx)
 154 {
 155   int i;
 156   static int inited = FALSE;
 157   static double *gfsumdata = nil;
 158
 159   if (! inited) {
 160     gfuncs.typx = nil;
 161     gfuncs.gfsum = nil;
 162     gfuncs.gderivs = nil;
 163     inited = TRUE;
 164   }
 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));
 169
 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;
 179 }
 180
 181 /* install constraint functions */
 182 static void
 183 install_cfuncs(char **g, int n, int k, double *ctarget, double h, double *typx)
 184 {
 185   int i;
 186   static int inited = FALSE;
 187
 188   if (! inited) {
 189     cfuncs.typx = nil;
 190     cfuncs.fsum = nil;
 191     cfuncs.gderivs = nil;
 192     inited = TRUE;
 193   }
 194   makespace(&cfuncs.typx, n * sizeof(double));
 195   makespace(&cfuncs.fsum, n * sizeof(double));
 196   makespace(&cfuncs.gderivs, k * sizeof(int));
 197
 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;
 206 }
 207
 208 static int
 209 in_support(char **ff, int n, double *x)
 210 {
 211   char *args[1], *values[1];
 212   int *result;
 213   char *mode[1];
 214   long length[1];
 215
 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]);
 225   }
 226 }
 227
 228 /* callback for logposterior evaluation */
 229 static int
 230 evalfunc(double *x, double *pval, double *grad, double **hess)
 231 {
 232   char *args[1], *values[3];
 233   double *result, val;
 234   char *mode[1];
 235   long length[1];
 236   int i, j;
 237
 238   for (i = 0; i < 3; i++) values[i] = nil;
 239
 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;
 252     }
 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];
 258       }
 259       else {
 260         numergrad(func.n, x, grad, func.fsum, evalfunc, func.h, func.typx);
 261       }
 262     }
 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];
 270       }
 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);
 276       }
 277     }
 278     return(TRUE);
 279   } else {
 280     if (pval != nil) {
 281       *pval = func.dflt;
 282     }
 283     return(FALSE);
 284   }
 285   return(TRUE);
 286 }
 287
 288
 289 /* callback for tilt function evaluation */
 290 static int which_gfunc;
 291
 292 static int
 293 evalgfunc(double *x, double *pval, double *grad, double **hess)
 294 {
 295   char *args[1], *values[3];
 296   double *result, val;
 297   char *mode[1];
 298   long length[1];
 299   int i, j;
 300
 301   for (i = 0; i < 3; i++) values[i] = nil;
 302
 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;
 314   }
 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];
 319     }
 320     else {
 321       numergrad(gfuncs.n, x, grad, gfuncs.gfsum[which_gfunc], evalgfunc,
 322                 gfuncs.h, gfuncs.typx);
 323     }
 324   }
 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];
 331     }
 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);
 339     }
 340   }
 341   return(TRUE);
 342 }
 343
 344 /* callback for constraint function evaluation */
 345 static int which_cfunc;
 346
 347 static int
 348 evalcfunc(double *x, double *pval, double *grad, double **hess)
 349 {
 350   char *args[1], *values[3];
 351   double *result, val;
 352   char *mode[1];
 353   long length[1];
 354   int i, j;
 355
 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];
 366     }
 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;
 370   }
 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];
 375     }
 376     else {
 377       numergrad(cfuncs.n, x, grad, cfuncs.fsum, evalcfunc,
 378                 cfuncs.h, cfuncs.typx);
 379     }
 380   }
 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];
 387     }
 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);
 395     }
 396   }
 397   return(TRUE);
 398 }
 399
 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)
 404 {
 405   int i;
 406   static double **hess = nil;
 407
 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;
 413   }
 414   evalfunc(x, val, grad, hess);
 415 }
 416
 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)
 421 {
 422   int i;
 423
 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;
 429   }
 430   return;
 431 }
 432
 433 /************************************************************************/
 434 /**                                                                    **/
 435 /**                     Derivative Scaling Routines                    **/
 436 /**                                                                    **/
 437 /************************************************************************/
 438
 439 /*
 440 void
 441 derivscalefront(char **ff, int *n, double *x, double *h, double *typx, double *tol, int *info)
 442 {
 443   int i;
 444
 445   if (*tol <= 0.0) *tol = pow(macheps(), GRADTOL_POWER);
 446
 447   install_func(ff[0], nil, *n, TRUE, 1.0, *h, typx, 0.0);
 448   *info = check_derivs(x, *tol);
 449
 450   *h = func.h;
 451   for (i = 0; i < *n; i++) typx[i] = func.typx[i];
 452   return;
 453 }
 454
 455
 456 static int
 457 check_derivs(double *x, double drvtol)
 458 {
 459   static double *grad = nil, work = nil;
 460   static double **hess = nil;
 461   int i, error;
 462
 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));
 466
 467   for (i = 0; i < func.n; i++) {
 468     hess[i] = work;
 469     work += func.n;
 470   }
 471
 472   error = derivscale(func.n, x, grad, hess, func.fsum, evalfunc,
 473                      &func.h, func.typx, drvtol, work);
 474   return(error);
 475 }
 476 */
 477
 478 /************************************************************************/
 479 /**                                                                    **/
 480 /**                    Importance Sampling Routines                    **/
 481 /**                                                                    **/
 482 /************************************************************************/
 483
 484 /* joint density of normal-cauchy mixture */
 485 static double
 486 dncmix(double *x, int n, double p)
 487 {
 488   int i;
 489   double dens;
 490
 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]);
 494   }
 495   return(dens);
 496 }
 497
 498 /*
 499  * S front end for computing sample from transformed normal-cauchy
 500  * mixture and importance sampling weights
 501  */
 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)
 506 {
 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];
 513
 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];
 521
 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;
 532     }
 533     if (evalfunc(y, &val, nil, nil)) w[i] = exp(val - mval) * c / dens;
 534     else w[i] = 0.0;
 535   }
 536   return;
 537 }
 538
 539
 540 /************************************************************************/
 541 /**                                                                    **/
 542 /**                       Maximization Routines                        **/
 543 /**                                                                    **/
 544 /************************************************************************/
 545
 546 typedef struct {
 547   int n, m, k, itnlimit, backtrack, verbose, vals_suppl, exptilt;
 548   int count, termcode;
 549 } MaxIPars;
 550
 551 typedef struct {
 552   double typf, h, gradtol, steptol, maxstep, dflt, tilt, newtilt, hessadd;
 553 } MaxDPars;
 554
 555 struct {
 556   double tilt;
 557   double *gval;
 558   double  **ggrad, **ghess;
 559   int exptilt;
 560   double *tscale;
 561 } tiltinfo;
 562
 563 static void
 564 add_tilt(double *x, double *pval, double *grad, double **hess,
 565          double tilt, int exptilt)
 566 {
 567   int i, j, k, n = func.n, m = gfuncs.k;
 568   double *gval, *ggrad, **ghess, etilt;
 569
 570   if (m == 0) return;
 571
 572   if (gfuncs.change_sign) tilt = -tilt;
 573
 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;
 578
 579     which_gfunc = k;
 580     evalgfunc(x, gval, ggrad, ghess);
 581
 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];
 590     }
 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]));
 605     }
 606   }
 607 }
 608
 609 static void
 610 set_tilt_info(int n, int m,
 611               double tilt, int exptilt, double *tscale)
 612 {
 613   static double *hessdata = nil, *graddata = nil;
 614   int i;
 615   static int inited = FALSE;
 616
 617   if (! inited) {
 618     tiltinfo.gval = nil;
 619     tiltinfo.ggrad = nil;
 620     tiltinfo.ghess = nil;
 621     inited = TRUE;
 622   }
 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));
 628
 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;
 634 }
 635
 636
 637 static void
 638 minfunc(double *x, double *pval, double *grad, double **hess)
 639 {
 640   int k = gfuncs.k;
 641
 642   if (evalfunc(x, pval, grad, hess) && (k > 0))
 643     add_tilt(x, pval, grad, hess, tiltinfo.tilt, tiltinfo.exptilt);
 644 }
 645
 646 void
 647 constfunc(double *x, double *vals, double **jac, double **hess)
 648 {
 649   int i, k = cfuncs.k;
 650   double *pvali, *jaci;
 651
 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);
 657   }
 658 }
 659
 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)
 665 {
 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)();
 671
 672   if (ipars->verbose > 0) PRINTSTR("maximizing...\n");
 673
 674   n = ipars->n;
 675   m = ipars->m;
 676   k = ipars->k;
 677   if (k >= n) Recover("too many constraints", NULL);
 678
 679   makespace(&H, n * sizeof(double *));
 680   makespace(&work, minworkspacesize(n, k));
 681
 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);
 686
 687   if (k == 0) {
 688     c = nil;
 689     cJ = nil;
 690     cfun = nil;
 691   }
 692   else {
 693     c = cvals;
 694     cvals += k;
 695     makespace(&cJ, k * sizeof(double *));
 696     for (i = 0; i < k; i++) cJ[i] = cvals + i * n;
 697     cfun = &constfunc; /* AJR: pointer to constfunc? */
 698   }
 699
 700   install_func(ff[0], nil, n, TRUE, dpars->typf, dpars->h, typx, dpars->dflt);
 701   install_gfuncs(gf, n, m, TRUE, dpars->h, typx);
 702   install_cfuncs(cf, n, k, ctarget, dpars->h, typx);
 703
 704   minsetup(n, k, &minfunc, &cfun, x, dpars->typf, typx, work); /* AJR: FIXME arg 3,4 */
 705   minsetoptions(dpars->gradtol, dpars->steptol, dpars->maxstep,
 706                 ipars->itnlimit, ipars->verbose, ipars->backtrack, TRUE);
 707
 708   if (ipars->vals_suppl) {
 709     for (i = 0; i < k; i++) c[i] -= ctarget[i];
 710     if (dpars->newtilt != dpars->tilt) {
 711       add_tilt(x, pf, grad, H, dpars->newtilt - dpars->tilt, ipars->exptilt);
 712       dpars->tilt = dpars->newtilt;
 713     }
 714     minsupplyvalues(*pf, grad, H, c, cJ);
 715   }
 716
 717   minimize();
 718   minresults(x, pf, nil, grad, H, c, cJ, &ipars->count, &ipars->termcode,
 719              &dpars->hessadd);
 720   msg[0] = minresultstring(ipars->termcode);
 721
 722   for (i = 0; i < k; i++) c[i] += ctarget[i];
 723   ipars->vals_suppl = TRUE;
 724 }
 725
 726 /************************************************************************/
 727 /**                                                                    **/
 728 /**                     Log Laplace Approximation                      **/
 729 /**                                                                    **/
 730 /************************************************************************/
 731
 732 void
 733 loglapdet(double *fvals, double *cvals, MaxIPars *ipars, MaxDPars *dpars,
 734           double *val, int *detonly)
 735 {
 736   int i, j, l, n = ipars->n, k = ipars->k;
 737   double f = -fvals[0], *hessdata = fvals + n + 1, *cgraddata = cvals + k;
 738   double ldL, ldcv, maxadd;
 739   static double **hess = nil, **cgrad = nil;
 740
 741   if (k >= n) Recover("too many constraints", NULL);
 742
 743   makespace(&hess, n * sizeof(double *));
 744   makespace(&cgrad, k * sizeof(double *));
 745
 746   for (i = 0; i < n; i++) hess[i] = hessdata + i * n;
 747   for (i = 0; i < k; i++) cgrad[i] = cgraddata + i * n;
 748
 749   choldecomp(hess, n, 0.0, &maxadd);
 750   /**** do something if not pos. definite ****/
 751
 752   for (i = 0, ldL = 0.0; i < n; i++) ldL += log(hess[i][i]);
 753
 754   if (k > 0) {
 755     /* forward solve for (L^-1) cgrad^T */
 756     for (l = 0; l < k; l++) {
 757       for (i = 0; i < n; i++) {
 758         if (hess[i][i] != 0.0) cgrad[l][i] /= hess[i][i];
 759         for (j = i + 1; j < n; j++) cgrad[l][j] -= hess[j][i] * cgrad[l][i];
 760       }
 761     }
 762
 763     /* compute sigma and stdev */
 764     for (i = 0; i < k; i++) {
 765       for (j = i; j < k; j++) {
 766         for (l = 0, hess[i][j] = 0.0; l < n; l++)
 767           hess[i][j] += cgrad[i][l] * cgrad[j][l];
 768         hess[j][i] = hess[i][j];
 769       }
 770     }
 771
 772     choldecomp(hess, k, 0.0, &maxadd);
 773     /**** do something if not pos. definite ****/
 774     for (i = 0, ldcv = 0.0; i < k; i++) ldcv += log(hess[i][i]);
 775   }
 776   else ldcv = 0.0;
 777
 778   *val = (n - k) * log(ROOT2PI) - ldL - ldcv;
 779   if (! *detonly) *val += f;
 780 }
 781
 782 #ifdef SBAYES
 783
 784 loglapfront(fvals, cvals, ipars, dpars, val)
 785      double *fvals, *cvals;
 786      MaxIPars *ipars;
 787      MaxDPars *dpars;
 788      double *val;
 789 {
 790   int detonly = FALSE;
 791
 792   loglapdet(fvals, cvals, ipars, dpars, val, &detonly);
 793 }
 794
 795 /************************************************************************/
 796 /**                                                                    **/
 797 /**                        First Order Moments                         **/
 798 /**                                                                    **/
 799 /************************************************************************/
 800
 801 moms1front(gf, n, m, x, hessdata, mean, stdev, sigmadata, h, typx)
 802      char *gf;
 803      int *n, *m;
 804      double *x, *hessdata, *mean, *stdev, *sigmadata, *h, *typx;
 805 {
 806   int i, j, k;
 807   double *hess, *sigma, *delg;
 808   double *delgdata, maxadd;
 809
 810   hess = (double **) S_alloc(*n, sizeof(double *));
 811   sigma = (double **) S_alloc(*m, sizeof(double *));
 812   delg = (double **) S_alloc(*m, sizeof(double *));
 813   delgdata = (double *) S_alloc(*m * *n, sizeof(double));
 814
 815   for (i = 0; i < *n; i++) hess[i] = hessdata + i * *n;
 816   for (i = 0; i < *m; i++) sigma[i] = sigmadata + i * *m;
 817   for (i = 0; i < *m; i++) delg[i] = delgdata + i * *n;
 818
 819   gevalfront(gf, n, m, x, h, typx, mean, delgdata);
 820
 821   /* get the cholesky decomposition L of the hessian */
 822   choldecomp(hess, *n, 0.0, &maxadd);
 823
 824   /* forward solve for (L^-1) delg^T */
 825   for (k = 0; k < *m; k++) {
 826     for (i = 0; i < *n; i++) {
 827       if (hess[i][i] != 0.0) delg[k][i] /= hess[i][i];
 828       for (j = i + 1; j < *n; j++) delg[k][j] -= hess[j][i] * delg[k][i];
 829     }
 830   }
 831
 832   /* compute sigma and stdev */
 833   for (i = 0; i < *m; i++) {
 834     for (j = i; j < *m; j++) {
 835       for (k = 0, sigma[i][j] = 0.0; k < *n; k++)
 836         sigma[i][j] += delg[i][k] * delg[j][k];
 837       sigma[j][i] = sigma[i][j];
 838     }
 839     stdev[i] = sqrt(sigma[i][i]);
 840   }
 841 }
 842
 843 /************************************************************************/
 844 /**                                                                    **/
 845 /**                        Second Order Moments                        **/
 846 /**                                                                    **/
 847 /************************************************************************/
 848
 849 typedef struct {
 850   MaxIPars max;
 851   int full, covar;
 852 } MomIPars;
 853
 854 typedef struct {
 855   MaxDPars max;
 856   double mgfdel;
 857 } MomDPars;
 858
 859 moms2front(ff, gf, gnum, x, typx, fvals, gvals, ipars, dpars,
 860            mean, stdev, sigmadata)
 861      char **ff, **gf;
 862      int *gnum;
 863      double *x, *typx, *fvals, *gvals, *mean, *stdev, *sigmadata;
 864      MomIPars *ipars;
 865      MomDPars *dpars;
 866 {
 867   char *msg;
 868   double h, loglap0, loglap1, loglap2;
 869   double *tilts, *fvals1, *gvals1, *x1;
 870   MomDPars dp1, *dpars1 = &dp1;
 871   MomIPars ip1, *ipars1 = &ip1;
 872   int i, n, m;
 873
 874   n = ipars->max.n;
 875   m = *gnum;
 876   h = dpars->max.h;
 877
 878   tilts = (double *) S_alloc(m, sizeof(double));
 879   x1 = (double *) S_alloc(n, sizeof(double));
 880   fvals1 = (double *) S_alloc(1 + n + n * n, sizeof(double));
 881   gvals1 = (double *) S_alloc(m + n * m, sizeof(double));
 882
 883   maxfront(ff, nil, nil, x, typx, fvals, gvals, nil, nil,
 884            ipars, dpars, nil, &msg);
 885   copypars(x, fvals, gvals, ipars, dpars, x1, fvals1, gvals1, ipars1, dpars1);
 886   loglapfront(fvals1, nil, ipars1, dpars1, &loglap0);
 887   copypars(x, fvals, gvals, ipars, dpars, x1, fvals1, gvals1, ipars1, dpars1);
 888   moms1front(gf, &n, &m, x1, fvals1 + n + 1, mean, stdev, sigmadata, &h, typx);
 889
 890   if (ipars->full) {
 891     for (i = 0; i < m; i++) {
 892       copypars(x, fvals, gvals, ipars, dpars,
 893                x1, fvals1, gvals1, ipars1, dpars1);
 894       ipars1->max.m = 1;
 895       ipars1->max.exptilt = FALSE;
 896       dpars1->max.newtilt = 1.0;
 897       maxfront(ff, gf + i, nil, x1, typx, fvals1, gvals1, nil, nil,
 898                ipars1, dpars1, nil, &msg);
 899       loglapfront(fvals1, nil, ipars1, dpars1, &loglap1);
 900       loglap1 -= loglap0;
 901
 902       copypars(x, fvals, gvals, ipars, dpars,
 903                x1, fvals1, gvals1, ipars1, dpars1);
 904       ipars1->max.m = 1;
 905       ipars1->max.exptilt = FALSE;
 906       dpars1->max.newtilt = 2.0;
 907       maxfront(ff, gf + i, nil, x1, typx, fvals1, gvals1, nil, nil,
 908                ipars1, dpars1, nil, &msg);
 909       loglapfront(fvals1, nil, ipars1, dpars1, &loglap2);
 910       loglap2 -= loglap0;
 911
 912       mean[i] = exp(loglap1);
 913       stdev[i] = sqrt(exp(loglap2) - exp(2.0 * loglap1));
 914       if (ipars->covar) sigmadata[i + i * m] = stdev[i] * stdev[i];
 915     }
 916     if (ipars->covar) {
 917       char *cgf[2];
 918       int j;
 919
 920       for (i = 0; i < m; i++) {
 921         for (j = i + 1; j < m; j++) {
 922           cgf[0] = gf[i];
 923           cgf[1] = gf[j];
 924           copypars(x, fvals, gvals, ipars, dpars,
 925                    x1, fvals1, gvals1, ipars1, dpars1);
 926           ipars1->max.m = 2;
 927           ipars1->max.exptilt = FALSE;
 928           dpars1->max.newtilt = 1.0;
 929           maxfront(ff, gf + i, nil, x1, typx, fvals1, gvals1, nil, nil,
 930                    ipars1, dpars1, nil, &msg);
 931           loglapfront(fvals1, nil, ipars1, dpars1, &loglap1);
 932           loglap1 -= loglap0;
 933
 934           sigmadata[i + j * m] = exp(loglap1) - mean[i] * mean[j];
 935           sigmadata[j + i * m] = sigmadata[i + j * m];
 936         }
 937       }
 938     }
 939   }
 940   else {
 941     for (i = 0; i < m; i++)
 942       tilts[i] = (stdev[i] > 0.0) ? dpars->mgfdel / stdev[i] : dpars->mgfdel;
 943
 944     for (i = 0; i < m; i++) {
 945       copypars(x, fvals, gvals, ipars, dpars,
 946                x1, fvals1, gvals1, ipars1, dpars1);
 947       ipars1->max.m = 1;
 948       ipars1->max.exptilt = TRUE;
 949       dpars1->max.newtilt = tilts[i];
 950       maxfront(ff, gf + i, nil, x1, typx, fvals1, gvals1, nil, nil,
 951                ipars1, dpars1, nil, &msg);
 952       loglapfront(fvals1, nil, ipars1, dpars1, &loglap1);
 953       loglap1 -= loglap0;
 954
 955       copypars(x, fvals, gvals, ipars, dpars,
 956                x1, fvals1, gvals1, ipars1, dpars1);
 957       ipars1->max.m = 1;
 958       ipars1->max.exptilt = TRUE;
 959       dpars1->max.newtilt = -tilts[i];
 960       maxfront(ff, gf + i, nil, x1, typx, fvals1, gvals1, nil, nil,
 961                ipars1, dpars1, nil, &msg);
 962       loglapfront(fvals1, nil, ipars1, dpars1, &loglap2);
 963       loglap2 -= loglap0;
 964
 965       mean[i] = (loglap1 - loglap2) / (2.0 * tilts[i]);
 966       stdev[i] = sqrt((loglap1 + loglap2) / (tilts[i] * tilts[i]));
 967       if (ipars->covar) sigmadata[i + i * m] = stdev[i] * stdev[i];
 968     }
 969     if (ipars->covar) {
 970       char *cgf[2];
 971       double ctilt, tscale[2];
 972       int j;
 973
 974       ctilt = dpars->mgfdel;
 975       for (i = 0; i < m; i++) {
 976         for (j = i + 1; j < m; j++) {
 977           cgf[0] = gf[i];
 978           cgf[1] = gf[j];
 979           tscale[0] = stdev[i] > 0 ? stdev[i] : 1.0;
 980           tscale[1] = stdev[j] > 0 ? stdev[j] : 1.0;
 981
 982           copypars(x, fvals, gvals, ipars, dpars,
 983                    x1, fvals1, gvals1, ipars1, dpars1);
 984           ipars1->max.m = 2;
 985           ipars1->max.exptilt = TRUE;
 986           dpars1->max.newtilt = ctilt;
 987           maxfront(ff, cgf, nil, x1, typx, fvals1, gvals1, nil, nil,
 988                    ipars1, dpars1, tscale, &msg);
 989           loglapfront(fvals1, nil, ipars1, dpars1, &loglap1);
 990           loglap1 -= loglap0;
 991
 992           copypars(x, fvals, gvals, ipars, dpars,
 993                    x1, fvals1, gvals1, ipars1, dpars1);
 994           ipars1->max.m = 2;
 995           ipars1->max.exptilt = TRUE;
 996           dpars1->max.newtilt = -ctilt;
 997           maxfront(ff, cgf, nil, x1, typx, fvals1, gvals1, nil, nil,
 998                    ipars1, dpars1, tscale, &msg);
 999           loglapfront(fvals1, nil, ipars1, dpars1, &loglap2);
1000           loglap2 -= loglap0;
1001
1002           sigmadata[i + j * m] = stdev[i] * stdev[j]
1003             * ((loglap2 + loglap1) /(2 * ctilt * ctilt) - 1.0);
1004           sigmadata[j + i * m] = sigmadata[i + j * m];
1005         }
1006       }
1007     }
1008   }
1009 }
1010
1011 static copypars(x, f, g, ip, dp, x1, f1, g1, ip1, dp1)
1012      double *x, *f, *g, *x1, *f1, *g1;
1013      MomIPars *ip, *ip1;
1014      MomDPars *dp, *dp1;
1015 {
1016   int i, n, m, nf, ng;
1017
1018   n = ip->max.n;
1019   m = ip->max.m;
1020   nf = 1 + n + n * n;
1021   ng = m + n * m;
1022
1023   for (i = 0; i < n; i++) x1[i] = x[i];
1024   for (i = 0; i < nf; i++) f1[i] = f[i];
1025   for (i = 0; i < ng; i++) g1[i] = g[i];
1026   *ip1 = *ip;
1027   *dp1 = *dp;
1028 }
1029
1030 /************************************************************************/
1031 /**                                                                    **/
1032 /**                          Laplace Margins                           **/
1033 /**                                                                    **/
1034 /************************************************************************/
1035
1036 lapmar1front(ff, gf, x, typx, fvals, ipars, dpars, xmar, ymar, nmar)
1037      char **ff, **gf;
1038      double *x, *typx, *fvals, *xmar, *ymar;
1039      MaxIPars *ipars;
1040      MaxDPars *dpars;
1041      int *nmar;
1042 {
1043   char *msg;
1044   int i, n, m, mindex;
1045   double h, loglap0, loglap1, xmode, stdev, sigmadata, ctarget[1];
1046   double *fvals1, *x1, *cvals, *cvals1, *fvals2, *x2, *cvals2;
1047   MaxDPars dp1, dp2, *dpars1 = &dp1, *dpars2 = &dp2;
1048   MaxIPars ip1, ip2, *ipars1 = &ip1, *ipars2 = &ip2;
1049
1050   n = ipars->n;
1051   m = 1;
1052   h = dpars->h;
1053
1054   x1 = (double *) S_alloc(n + 1, sizeof(double));
1055   fvals1 = (double *) S_alloc(1 + n + n * n, sizeof(double));
1056   cvals = (double *) S_alloc(m + n * m, sizeof(double));
1057   cvals1 = (double *) S_alloc(m + n * m, sizeof(double));
1058   x2 = (double *) S_alloc(n + 1, sizeof(double));
1059   fvals2 = (double *) S_alloc(1 + n + n * n, sizeof(double));
1060   cvals2 = (double *) S_alloc(m + n * m, sizeof(double));
1061
1062   maxfront(ff, nil, nil, x, typx, fvals, nil, nil, nil,
1063            ipars, dpars, nil, &msg);
1064   cpmarpars(x, fvals, cvals, ipars, dpars, x1, fvals1, cvals1, ipars1, dpars1);
1065   loglapfront(fvals1, nil, ipars1, dpars1, &loglap0);
1066   cpmarpars(x, fvals, cvals, ipars, dpars, x1, fvals1, cvals1, ipars1, dpars1);
1067   moms1front(gf, &n, &m, x1, fvals1 + n + 1, &xmode, &stdev, &sigmadata,
1068              &h, typx);
1069
1070   ipars->k = 1;
1071   ipars->vals_suppl = FALSE;
1072   ctarget[0] = xmode;
1073   maxfront(ff, nil, gf, x, typx, fvals, nil, cvals, ctarget,
1074            ipars, dpars, nil, &msg);
1075
1076   for (mindex = 0; mindex < *nmar && xmar[mindex] <= xmode; mindex++);
1077
1078   cpmarpars(x, fvals, cvals, ipars, dpars, x1, fvals1, cvals1, ipars1, dpars1);
1079   for (i = mindex; i >= 0; i--) {
1080     ctarget[0] = xmar[i];
1081     maxfront(ff, nil, gf, x1, typx, fvals1, nil, cvals1, ctarget,
1082              ipars1, dpars1, nil, &msg);
1083     cpmarpars(x1, fvals1, cvals1, ipars1, dpars1, x2,
1084               fvals2, cvals2, ipars2, dpars2);
1085     loglapfront(fvals2, cvals2, ipars2, dpars2, &loglap1);
1086     ymar[i] = exp(loglap1 - loglap0);
1087   }
1088   cpmarpars(x, fvals, cvals, ipars, dpars, x1, fvals1, cvals1, ipars1, dpars1);
1089   for (i = mindex + 1; i < *nmar; i++) {
1090     ctarget[0] = xmar[i];
1091     maxfront(ff, nil, gf, x1, typx, fvals1, nil, cvals1, ctarget,
1092              ipars1, dpars1, nil, &msg);
1093     cpmarpars(x1, fvals1, cvals1, ipars1, dpars1, x2,
1094               fvals2, cvals2, ipars2, dpars2);
1095     loglapfront(fvals2, cvals2, ipars2, dpars2, &loglap1);
1096     ymar[i] = exp(loglap1 - loglap0);
1097   }
1098 }
1099
1100 static cpmarpars(x, f, g, ip, dp, x1, f1, g1, ip1, dp1)
1101      double *x, *f, *g, *x1, *f1, *g1;
1102      MaxIPars *ip, *ip1;
1103      MaxDPars *dp, *dp1;
1104 {
1105   int i, n, k, nf, ng;
1106
1107   n = ip->n;
1108   k = ip->k;
1109   nf = 1 + n + n * n;
1110   ng = k + n * k;
1111
1112   for (i = 0; i < n; i++) x1[i] = x[i];
1113   for (i = 0; i < nf; i++) f1[i] = f[i];
1114   for (i = 0; i < ng; i++) g1[i] = g[i];
1115   *ip1 = *ip;
1116   *dp1 = *dp;
1117 }
1118 #endif /* SBAYES */
1119
1120 /*
1121   TODO
1122
1123   get hessian from gradiant for analytical gradiants
1124
1125   avoid repeated derivative calls in mimimize.
1126
1127   2d margins
1128
1129   use pos. definiteness info in margins
1130
1131 */