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 #ifdef SBAYES
   8 # include <math.h>
   9 #define PRINTSTR(s) printf(s)
  10 #else
  11 # include "xmath.h"
  12 #define PRINTSTR(s) stdputstr(s)
  13 #endif SBAYES
  14
  15 extern char *S_alloc(), *calloc(), *realloc();
  16 extern double macheps();
  17 char *minresultstring();
  18
  19 /************************************************************************/
  20 /**                                                                    **/
  21 /**                      Definitions and Globals                       **/
  22 /**                                                                    **/
  23 /************************************************************************/
  24
  25 #define nil 0L
  26 #define NULL 0L
  27 #define TRUE 1
  28 #define FALSE 0
  29
  30 #define ROOT2PI 2.50662827463100050241
  31 #define PI_INV .31830988618379067153
  32
  33 #define GRADTOL_POWER 1.0 / 3.0
  34 #define H_POWER 1.0 / 6.0
  35
  36 typedef double **RMatrix, *RVector;
  37
  38 typedef struct{
  39   char *f, **sf, **g;
  40   int n, k;
  41   int change_sign, fderivs;
  42   int *gderivs;
  43   double typf, h, dflt;
  44   RVector typx, fsum, cvals, ctarget;
  45   RMatrix gfsum;
  46 } Fundata;
  47
  48 static Fundata func, gfuncs, cfuncs;
  49
  50 /************************************************************************/
  51 /**                                                                    **/
  52 /**                         Memory Utilities                           **/
  53 /**                                                                    **/
  54 /************************************************************************/
  55
  56 /* this function is used to maintain a statically allocated piece of    */
  57 /* memory of a specified size. If a larger piece is needed the pointer  */
  58 /* is realloced. This allows functions using memory allocation to be    */
  59 /* called repeatedly (but not recursively) from within the same call    */
  60 /* from S. It attempts to avoid the danger of dangling callocs.         */
  61
  62 static makespace(pptr, size)
  63      char **pptr;
  64      int size;
  65 {
  66   if (size <= 0) return;
  67   if (*pptr == nil) *pptr = calloc(size, 1);
  68   else *pptr = realloc(*pptr, size);
  69   if (size > 0 && *pptr == nil) Recover("memory allocation failed", NULL);
  70 }
  71
  72 /************************************************************************/
  73 /**                                                                    **/
  74 /**                    Functions Evaluation Routines                   **/
  75 /**                                                                    **/
  76 /************************************************************************/
  77
  78 /*
  79  * All Hessianevaluations by numerical derivatives assume the gradient is
  80  * evaluated first at the same location. The results are cached away.
  81  */
  82
  83 /* install log posterior function */
  84 static install_func(f, sf, n, change_sign, typf, h, typx, dflt)
  85      char *f, **sf;
  86      int n;
  87      double typf, h, dflt;
  88      RVector typx;
  89 {
  90   int i;
  91   static int inited = FALSE;
  92
  93   if (! inited) {
  94     func.typx = nil;
  95     func.fsum = nil;
  96     inited = TRUE;
  97   }
  98   makespace(&func.typx, n * sizeof(double));
  99   makespace(&func.fsum, n * sizeof(double));
 100
 101   func.f = f;
 102   func.sf = sf;
 103   func.n = n;
 104   func.change_sign = change_sign;
 105   func.typf = (typf > 0.0) ? typf : 1.0;
 106   func.h = (h > 0.0) ? h : pow(macheps(), H_POWER);
 107   for (i = 0; i < n; i++)
 108     func.typx[i] = (typx != nil && typx[i] > 0.0) ? typx[i] : 1.0;
 109   func.dflt = dflt;
 110   func.fderivs = 0;
 111 }
 112
 113 /* install tilt functions */
 114 static install_gfuncs(g, n, k, change_sign, h, typx)
 115      char **g;
 116      int n, k, change_sign;
 117      double h;
 118      RVector typx;
 119 {
 120   int i;
 121   static int inited = FALSE;
 122   static double *gfsumdata = nil;
 123
 124   if (! inited) {
 125     gfuncs.typx = nil;
 126     gfuncs.gfsum = nil;
 127     gfuncs.gderivs = nil;
 128     inited = TRUE;
 129   }
 130   makespace(&gfuncs.typx, n * sizeof(double));
 131   makespace(&gfuncs.gfsum, k * sizeof(double *));
 132   makespace(&gfsumdata, k * n * sizeof(double));
 133   makespace(&gfuncs.gderivs, k *sizeof(int));
 134
 135   gfuncs.g = g;
 136   gfuncs.n = n;
 137   gfuncs.k = k;
 138   gfuncs.change_sign = change_sign;
 139   gfuncs.h = (h > 0.0) ? h : pow(macheps(), H_POWER);
 140   for (i = 0; i < n; i++)
 141     gfuncs.typx[i] = (typx != nil && typx[i] > 0.0) ? typx[i] : 1.0;
 142   for (i = 0; i < k; i++) gfuncs.gfsum[i] = gfsumdata + i * n;
 143 }
 144
 145 /* install constraint functions */
 146 static install_cfuncs(g, n, k, ctarget, h, typx)
 147      char **g;
 148      int n, k;
 149      double h;
 150      RVector typx, ctarget;
 151 {
 152   int i;
 153   static int inited = FALSE;
 154
 155   if (! inited) {
 156     cfuncs.typx = nil;
 157     cfuncs.fsum = nil;
 158     cfuncs.gderivs = nil;
 159     inited = TRUE;
 160   }
 161   makespace(&cfuncs.typx, n * sizeof(double));
 162   makespace(&cfuncs.fsum, n * sizeof(double));
 163   makespace(&cfuncs.gderivs, k * sizeof(int));
 164
 165   cfuncs.g = g;
 166   cfuncs.n = n;
 167   cfuncs.k = k;
 168   cfuncs.h = (h > 0.0) ? h : pow(macheps(), H_POWER);
 169   for (i = 0; i < n; i++)
 170     cfuncs.typx[i] = (typx != nil && typx[i] > 0.0) ? typx[i] : 1.0;
 171   cfuncs.ctarget = ctarget;
 172 }
 173
 174 /* callback to test if x is in the support of the posterior */
 175 static in_support(ff, n, x)
 176      char **ff;
 177      int n;
 178      double *x;
 179 {
 180   char *args[1], *values[1];
 181   int *result;
 182   char *mode[1];
 183   long length[1];
 184
 185   if (ff == nil || ff[0] == nil) return(TRUE);
 186   else {
 187     mode[0] = "double";
 188     length[0] =n;
 189     args[0] = (char *) x;
 190     call_S(ff[0], 1L, args, mode, length, 0L, 1L, values);
 191     result = (int *) values[0];
 192     return(result[0]);
 193   }
 194 }
 195
 196 /* callback for logposterior evaluation */
 197 static evalfunc(x, pval, grad, hess)
 198      RVector x, grad;
 199      double *pval;
 200      RMatrix hess;
 201 {
 202   char *args[1], *values[3];
 203   double *result, val;
 204   char *mode[1];
 205   long length[1];
 206   int i, j;
 207
 208   for (i = 0; i < 3; i++) values[i] = nil;
 209
 210   if (in_support(func.sf, func.n, x)) {
 211     if (pval != nil || func.fderivs > 0 || hess != nil) {
 212       mode[0] = "double";
 213       length[0] = func.n;
 214       args[0] = (char *) x;
 215       call_S(func.f, 1L, args, mode, length, 0L, 3L, values);
 216       result = (double *) values[0];
 217           val = (! func.change_sign) ? result[0] : -result[0];
 218       if (pval != nil) *pval = val;
 219       if (values[2] != nil) func.fderivs = 2;
 220       else if (values[1] != nil) func.fderivs = 1;
 221       else func.fderivs = 0;
 222     }
 223     if (grad != nil) {
 224       if (func.fderivs > 0) {
 225         result = (double *) values[1];
 226         for (i = 0; i < func.n; i++)
 227           grad[i] = (! func.change_sign) ? result[i] : -result[i];
 228       }
 229       else {
 230         numergrad(func.n, x, grad, func.fsum, evalfunc, func.h, func.typx);
 231       }
 232     }
 233     if (hess != nil) {
 234       if (func.fderivs == 2) {
 235         result = (double *) values[2];
 236         for (i = 0; i < func.n; i++)
 237           for (j = 0; j < func.n; j++)
 238             hess[i][j] = (! func.change_sign) ? result[i + j * func.n]
 239                                               : -result[i + j * func.n];
 240       }
 241       else {
 242         if (func.fderivs == 1) /* kludge to get fsum for analytic gradients */
 243           numergrad(func.n, x, func.fsum, func.fsum,
 244                     evalfunc, func.h, func.typx);
 245         numerhess(func.n, x, hess, val, func.fsum, evalfunc, func.h, func.typx);
 246       }
 247     }
 248     return(TRUE);
 249   }
 250   else {
 251     if (pval != nil) *pval = func.dflt;
 252     return(FALSE);
 253   }
 254 }
 255
 256
 257 /* callback for tilt function evaluation */
 258 static int which_gfunc;
 259
 260 static evalgfunc(x, pval, grad, hess)
 261      RVector x, grad;
 262      double *pval;
 263      RMatrix hess;
 264 {
 265   char *args[1], *values[3];
 266   double *result, val;
 267   char *mode[1];
 268   long length[1];
 269   int i, j;
 270
 271   for (i = 0; i < 3; i++) values[i] = nil;
 272
 273   if (pval != nil || gfuncs.gderivs[which_gfunc] > 0 || hess != nil) {
 274     mode[0] = "double";
 275     length[0] = gfuncs.n;
 276     args[0] = (char *) x;
 277     call_S(gfuncs.g[which_gfunc], 1L, args, mode, length, 0L, 3L, values);
 278     result = (double *) values[0];
 279         val = result[0];
 280     if (pval != nil) *pval = result[0];
 281     if (values[2] != nil) gfuncs.gderivs[which_gfunc] = 2;
 282     else if (values[1] != nil) gfuncs.gderivs[which_gfunc] = 1;
 283     else gfuncs.gderivs[which_gfunc] = 0;
 284   }
 285   if (grad != nil) {
 286     if (gfuncs.gderivs[which_gfunc] > 0) {
 287       result = (double *) values[1];
 288       for (i = 0; i < gfuncs.n; i++) grad[i] = result[i];
 289     }
 290     else {
 291       numergrad(gfuncs.n, x, grad, gfuncs.gfsum[which_gfunc], evalgfunc,
 292                 gfuncs.h, gfuncs.typx);
 293     }
 294   }
 295   if (hess != nil) {
 296     if (gfuncs.gderivs[which_gfunc] == 2) {
 297       result = (double *) values[2];
 298       for (i = 0; i < gfuncs.n; i++)
 299         for (j = 0; j < gfuncs.n; j++)
 300           hess[i][j] = result[i + j * gfuncs.n];
 301     }
 302     else {
 303       /* kludge to get fsum if analytic gradient used */
 304       if (gfuncs.gderivs[which_gfunc] == 1)
 305         numergrad(gfuncs.n, x, gfuncs.gfsum[which_gfunc],
 306                   gfuncs.gfsum[which_gfunc], evalgfunc, gfuncs.h, gfuncs.typx);
 307       numerhess(gfuncs.n, x, hess, val, gfuncs.gfsum[which_gfunc], evalgfunc,
 308                 gfuncs.h, gfuncs.typx);
 309     }
 310   }
 311 }
 312
 313 /* callback for constraint function evaluation */
 314 static int which_cfunc;
 315
 316 static evalcfunc(x, pval, grad, hess)
 317      RVector x, grad;
 318      double *pval;
 319      RMatrix hess;
 320 {
 321   char *args[1], *values[3];
 322   double *result, val;
 323   char *mode[1];
 324   long length[1];
 325   int i, j;
 326
 327   if (pval != nil || cfuncs.gderivs[which_cfunc] > 0 || hess != nil) {
 328     mode[0] = "double";
 329     length[0] = cfuncs.n;
 330     args[0] = (char *) x;
 331     call_S(cfuncs.g[which_cfunc], 1L, args, mode, length, 0L, 3L, values);
 332     result = (double *) values[0];
 333         val = result[0];
 334     if (pval != nil) {
 335       *pval = result[0];
 336       if (cfuncs.ctarget != nil) *pval -= cfuncs.ctarget[which_cfunc];
 337     }
 338     if (values[2] != nil) cfuncs.gderivs[which_cfunc] = 2;
 339     else if (values[1] != nil) cfuncs.gderivs[which_cfunc] = 1;
 340     else cfuncs.gderivs[which_cfunc] = 0;
 341   }
 342   if (grad != nil) {
 343     if (cfuncs.gderivs[which_cfunc] > 0) {
 344       result = (double *) values[1];
 345       for (i = 0; i <cfuncs.n; i++) grad[i] = result[i];
 346     }
 347     else {
 348       numergrad(cfuncs.n, x, grad, cfuncs.fsum, evalcfunc,
 349                 cfuncs.h, cfuncs.typx);
 350     }
 351   }
 352   if (hess != nil) {
 353     if (cfuncs.gderivs[which_cfunc] == 2) {
 354       result = (double *) values[2];
 355       for (i = 0; i <cfuncs.n; i++)
 356         for (j = 0; j <cfuncs.n; j++)
 357           hess[i][j] = result[i + j * cfuncs.n];
 358     }
 359     else {
 360       /* kludge to get fsum if analytic gradient used */
 361       if (cfuncs.gderivs[which_cfunc] == 1)
 362         numergrad(cfuncs.n, x, cfuncs.fsum, cfuncs.fsum, evalcfunc,
 363                   cfuncs.h, cfuncs.typx);
 364       numerhess(cfuncs.n, x, hess, val, cfuncs.fsum, evalcfunc,
 365                 cfuncs.h, cfuncs.typx);
 366     }
 367   }
 368 }
 369
 370 /* S front end for logposterior evaluation */
 371 evalfront(ff, n, x, val, grad, phess, h, typx)
 372      char **ff;
 373      int *n;
 374      double *x, *val, *grad, *phess, *typx, *h;
 375 {
 376   int i;
 377   static RMatrix hess = nil;
 378
 379   install_func(ff[0], nil, *n, FALSE, 1.0, *h, typx, 0.0);
 380   if (phess == nil) hess = nil;
 381   else {
 382     makespace(&hess, *n * sizeof(double *));
 383     for (i = 0; i < *n; i++, phess += *n) hess[i] = phess;
 384   }
 385   evalfunc(x, val, grad, hess);
 386 }
 387
 388 #ifdef SBAYES
 389 /* S front end for tilt function evaluation */
 390 gevalfront(gg, n, m, x, h, typx, val, grad)
 391      char **gg;
 392      int *n, *m;
 393      double *x, *h, *typx, *val, *grad;
 394 {
 395   int i;
 396
 397   install_gfuncs(gg, *n, *m, FALSE, *h, typx);
 398   for (i = 0; i < *m; i++, val++) {
 399     which_gfunc = i;
 400     evalgfunc(x, val, grad, nil);
 401     if (grad != nil) grad += *n;
 402   }
 403 }
 404
 405 /************************************************************************/
 406 /**                                                                    **/
 407 /**                     Derivative Scaling Routines                    **/
 408 /**                                                                    **/
 409 /************************************************************************/
 410
 411 static check_derivs(x, drvtol)
 412      RVector x;
 413      double drvtol;
 414 {
 415   static RVector grad = nil, work = nil;
 416   static RMatrix hess = nil;
 417   int i, error;
 418
 419   grad = (RVector) S_alloc(func.n, sizeof(double));
 420   hess = (RMatrix) S_alloc(func.n, sizeof(double *));
 421   work = (RVector) S_alloc(func.n + func.n * func.n, sizeof(double));
 422
 423   for (i = 0; i < func.n; i++) {
 424     hess[i] = work;
 425     work += func.n;
 426   }
 427
 428   error = derivscale(func.n, x, grad, hess, func.fsum, evalfunc,
 429                      &func.h, func.typx, drvtol, work);
 430   return(error);
 431 }
 432
 433 derivscalefront(ff, n, x, h, typx, tol, info)
 434      char **ff;
 435      int *n, *info;
 436      double *x, *h, *typx, *tol;
 437 {
 438   int i;
 439
 440   if (*tol <= 0.0) *tol = pow(macheps(), GRADTOL_POWER);
 441
 442   install_func(ff[0], nil, *n, TRUE, 1.0, *h, typx, 0.0);
 443   *info = check_derivs(x, *tol);
 444
 445   *h = func.h;
 446   for (i = 0; i < *n; i++) typx[i] = func.typx[i];
 447 }
 448
 449 /************************************************************************/
 450 /**                                                                    **/
 451 /**                    Importance Sampling Routines                    **/
 452 /**                                                                    **/
 453 /************************************************************************/
 454
 455 /* joint density of normal-cauchy mixture */
 456 static double dncmix(x, n, p)
 457      double *x, p;
 458      int n;
 459 {
 460   int i;
 461   double dens;
 462
 463   for (i = 0, dens = 1.0; i < n; i++) {
 464     dens *= p * exp(-0.5 * x[i] * x[i]) / ROOT2PI
 465       + (1 - p) * PI_INV / (1.0 + x[i] * x[i]);
 466   }
 467   return(dens);
 468 }
 469
 470 /*
 471  * S front end for computing sample from transformed normal-cauchy
 472  * mixture and importance sampling weights
 473  */
 474 samplefront(ff, sf, rf, p, n, x, ch, N, y, w)
 475      char **ff, **sf, **rf;
 476      int *n, *N;
 477      double *p,  *x, *ch, *y, *w;
 478 {
 479   double val;
 480   int i, j, k;
 481   char *args[2], *values[1];
 482   double *result, mval, c, dens, m;
 483   char *mode[2];
 484   long length[2];
 485
 486   /* get the random variables */
 487   mode[0] = "double"; mode[1] = "double";
 488   length[0] = 1; length[1] = 1;
 489   m = *N * *n; args[0] = (char *) &m; args[1] = (char *) p;
 490   call_S(rf[0], 2L, args, mode, length, 0L, 1L, values);
 491   result = (double *) values[0];
 492   for (i = 0; i < m; i++) y[i] = result[i];
 493
 494   /* construct the sample and the weights */
 495   install_func(ff[0], sf, *n, FALSE, 1.0, -1.0, nil, 0.0);
 496   c = 1.0 / pow(ROOT2PI, (double) *n);
 497   evalfunc(x, &mval, nil, nil);
 498   for (i = 0; i < *N; i++, y += *n) {
 499     dens = dncmix(y, *n, *p);
 500     for (j = 0; j < *n; j++) {
 501       val = x[j];
 502       for (k = j; k < *n; k++) val += y[k] * ch[j + *n * k];
 503       y[j] = val;
 504     }
 505     if (evalfunc(y, &val, nil, nil)) w[i] = exp(val - mval) * c / dens;
 506     else w[i] = 0.0;
 507   }
 508 }
 509 #endif SBAYES
 510 /************************************************************************/
 511 /**                                                                    **/
 512 /**                       Maximization Routines                        **/
 513 /**                                                                    **/
 514 /************************************************************************/
 515
 516 typedef struct {
 517   int n, m, k, itnlimit, backtrack, verbose, vals_suppl, exptilt;
 518   int count, termcode;
 519 } MaxIPars;
 520
 521 typedef struct {
 522   double typf, h, gradtol, steptol, maxstep, dflt, tilt, newtilt, hessadd;
 523 } MaxDPars;
 524
 525 struct {
 526   double tilt;
 527   RVector gval;
 528   RMatrix  ggrad, ghess;
 529   int exptilt;
 530   RVector tscale;
 531 } tiltinfo;
 532
 533 static set_tilt_info(n, m, tilt, exptilt, tscale)
 534      int n, m;
 535      double tilt, *tscale;
 536      int exptilt;
 537 {
 538   static double *hessdata = nil, *graddata = nil;
 539   int i;
 540   static int inited = FALSE;
 541
 542   if (! inited) {
 543     tiltinfo.gval = nil;
 544     tiltinfo.ggrad = nil;
 545     tiltinfo.ghess = nil;
 546     inited = TRUE;
 547   }
 548   makespace(&tiltinfo.gval, n * sizeof(double));
 549   makespace(&tiltinfo.ggrad, m * sizeof(double *));
 550   makespace(&tiltinfo.ghess, n * sizeof(double *));
 551   makespace(&graddata, n * m * sizeof(double));
 552   makespace(&hessdata, n * n * sizeof(double));
 553
 554   tiltinfo.tilt = tilt;
 555   tiltinfo.exptilt = exptilt;
 556   for (i = 0; i < m; i++) tiltinfo.ggrad[i] = graddata + i * n;
 557   for (i = 0; i < n; i++) tiltinfo.ghess[i] = hessdata + i * n;
 558   tiltinfo.tscale = tscale;
 559 }
 560
 561 static minfunc(x, pval, grad, hess)
 562      RVector x, grad;
 563      double *pval;
 564      RMatrix hess;
 565 {
 566   int k = gfuncs.k;
 567
 568   if (evalfunc(x, pval, grad, hess) && (k > 0))
 569     add_tilt(x, pval, grad, hess, tiltinfo.tilt, tiltinfo.exptilt);
 570 }
 571
 572 constfunc(x, vals, jac, hess)
 573      RVector x, vals;
 574      RMatrix jac, hess;
 575 {
 576   int i, k = cfuncs.k;
 577   double *pvali, *jaci;
 578
 579   for (i = 0; i < k; i++) {
 580     pvali = (vals != nil) ? vals + i : nil;
 581     jaci = (jac != nil) ? jac[i] : nil;
 582     which_cfunc = i;
 583     evalcfunc(x, pvali, jaci, nil);
 584   }
 585 }
 586
 587 static add_tilt(x, pval, grad, hess, tilt, exptilt)
 588      RVector x, grad;
 589      double *pval, tilt;
 590      RMatrix hess;
 591      int exptilt;
 592 {
 593   int i, j, k, n = func.n, m = gfuncs.k;
 594   double *gval, *ggrad, **ghess, etilt;
 595
 596   if (m == 0) return;
 597
 598   if (gfuncs.change_sign) tilt = -tilt;
 599
 600   for (k = 0; k < m; k++) {
 601     gval = (pval != nil) ? tiltinfo.gval + k : nil;
 602     ggrad = (grad != nil) ? tiltinfo.ggrad[k] : nil;
 603     ghess = (hess != nil) ? tiltinfo.ghess : nil;
 604
 605     which_gfunc = k;
 606     evalgfunc(x, gval, ggrad, ghess);
 607
 608     if (exptilt) {
 609       etilt = (tiltinfo.tscale != nil) ? tilt / tiltinfo.tscale[k] : tilt;
 610       if (pval != nil) *pval += etilt * *gval;
 611       if (grad != nil)
 612         for (i = 0; i < n; i++) grad[i] += etilt * ggrad[i];
 613       if (hess != nil)
 614         for (i = 0; i < n; i++)
 615           for (j = 0; j < n; j++) hess[i][j] += etilt * ghess[i][j];
 616     }
 617     else {
 618       gval = tiltinfo.gval;
 619       ggrad = tiltinfo.ggrad[k];
 620       ghess = tiltinfo.ghess;
 621       if (gval[k] <= 0.0) Recover("nonpositive function value", NULL);
 622       if (pval != nil) *pval += tilt * log(gval[k]);
 623       if (grad != nil)
 624         for (i = 0; i < n; i++) grad[i] += tilt * ggrad[i] / gval[k];
 625       if (hess != nil)
 626         for (i = 0; i < n; i++)
 627           for (j = 0; j < n; j++)
 628             hess[i][j] +=
 629               tilt * (ghess[i][j] / gval[k]
 630                       - (ggrad[i] / gval[k]) * (ggrad[j] / gval[k]));
 631     }
 632   }
 633 }
 634
 635 maxfront(ff, gf, cf, x, typx, fvals, gvals, cvals, ctarget, ipars, dpars,
 636          tscale, msg)
 637      char **ff, **gf, **cf;
 638      double *x, *typx, *fvals, *gvals, *cvals, *ctarget, *tscale;
 639      MaxIPars *ipars;
 640      MaxDPars *dpars;
 641      char **msg;
 642 {
 643   static char *work = nil;
 644   static RMatrix H = nil, cJ = nil;
 645   double *pf, *grad, *c;
 646   int i, n, m, k;
 647   int (*cfun)();
 648
 649   if (ipars->verbose > 0) PRINTSTR("maximizing...\n");
 650
 651   n = ipars->n;
 652   m = ipars->m;
 653   k = ipars->k;
 654   if (k >= n) Recover("too many constraints", NULL);
 655
 656   makespace(&H, n * sizeof(double *));
 657   makespace(&work, minworkspacesize(n, k));
 658
 659   pf = fvals; fvals++;
 660   grad = fvals; fvals += n;
 661   for (i = 0; i < n; i++, fvals += n) H[i] = fvals;
 662   set_tilt_info(n, m, dpars->newtilt, ipars->exptilt, tscale);
 663
 664   if (k == 0) {
 665     c = nil;
 666     cJ = nil;
 667     cfun = nil;
 668   }
 669   else {
 670     c = cvals;
 671     cvals += k;
 672     makespace(&cJ, k * sizeof(double *));
 673     for (i = 0; i < k; i++) cJ[i] = cvals + i * n;
 674     cfun = constfunc;
 675   }
 676
 677   install_func(ff[0], nil, n, TRUE, dpars->typf, dpars->h, typx, dpars->dflt);
 678   install_gfuncs(gf, n, m, TRUE, dpars->h, typx);
 679   install_cfuncs(cf, n, k, ctarget, dpars->h, typx);
 680
 681   minsetup(n, k, minfunc, cfun, x, dpars->typf, typx, work);
 682   minsetoptions(dpars->gradtol, dpars->steptol, dpars->maxstep,
 683                 ipars->itnlimit, ipars->verbose, ipars->backtrack, TRUE);
 684
 685   if (ipars->vals_suppl) {
 686     for (i = 0; i < k; i++) c[i] -= ctarget[i];
 687     if (dpars->newtilt != dpars->tilt) {
 688       add_tilt(x, pf, grad, H, dpars->newtilt - dpars->tilt, ipars->exptilt);
 689       dpars->tilt = dpars->newtilt;
 690     }
 691     minsupplyvalues(*pf, grad, H, c, cJ);
 692   }
 693
 694   minimize();
 695   minresults(x, pf, nil, grad, H, c, cJ, &ipars->count, &ipars->termcode,
 696              &dpars->hessadd);
 697   msg[0] = minresultstring(ipars->termcode);
 698
 699   for (i = 0; i < k; i++) c[i] += ctarget[i];
 700   ipars->vals_suppl = TRUE;
 701 }
 702
 703 /************************************************************************/
 704 /**                                                                    **/
 705 /**                     Log Laplace Approximation                      **/
 706 /**                                                                    **/
 707 /************************************************************************/
 708
 709 loglapdet(fvals, cvals, ipars, dpars, val, detonly)
 710      double *fvals, *cvals;
 711      MaxIPars *ipars;
 712      MaxDPars *dpars;
 713      double *val;
 714      int *detonly;
 715 {
 716   int i, j, l, n = ipars->n, k = ipars->k;
 717   double f = -fvals[0], *hessdata = fvals + n + 1, *cgraddata = cvals + k;
 718   double ldL, ldcv, maxadd;
 719   static RMatrix hess = nil, cgrad = nil;
 720
 721   if (k >= n) Recover("too many constraints", NULL);
 722
 723   makespace(&hess, n * sizeof(double *));
 724   makespace(&cgrad, k * sizeof(double *));
 725
 726   for (i = 0; i < n; i++) hess[i] = hessdata + i * n;
 727   for (i = 0; i < k; i++) cgrad[i] = cgraddata + i * n;
 728
 729   choldecomp(hess, n, 0.0, &maxadd);
 730   /**** do something if not pos. definite ****/
 731
 732   for (i = 0, ldL = 0.0; i < n; i++) ldL += log(hess[i][i]);
 733
 734   if (k > 0) {
 735     /* forward solve for (L^-1) cgrad^T */
 736     for (l = 0; l < k; l++) {
 737       for (i = 0; i < n; i++) {
 738         if (hess[i][i] != 0.0) cgrad[l][i] /= hess[i][i];
 739         for (j = i + 1; j < n; j++) cgrad[l][j] -= hess[j][i] * cgrad[l][i];
 740       }
 741     }
 742
 743     /* compute sigma and stdev */
 744     for (i = 0; i < k; i++) {
 745       for (j = i; j < k; j++) {
 746         for (l = 0, hess[i][j] = 0.0; l < n; l++)
 747           hess[i][j] += cgrad[i][l] * cgrad[j][l];
 748         hess[j][i] = hess[i][j];
 749       }
 750     }
 751
 752     choldecomp(hess, k, 0.0, &maxadd);
 753     /**** do something if not pos. definite ****/
 754     for (i = 0, ldcv = 0.0; i < k; i++) ldcv += log(hess[i][i]);
 755   }
 756   else ldcv = 0.0;
 757
 758   *val = (n - k) * log(ROOT2PI) - ldL - ldcv;
 759   if (! *detonly) *val += f;
 760 }
 761
 762 #ifdef SBAYES
 763
 764 loglapfront(fvals, cvals, ipars, dpars, val)
 765      double *fvals, *cvals;
 766      MaxIPars *ipars;
 767      MaxDPars *dpars;
 768      double *val;
 769 {
 770   int detonly = FALSE;
 771
 772   loglapdet(fvals, cvals, ipars, dpars, val, &detonly);
 773 }
 774
 775 /************************************************************************/
 776 /**                                                                    **/
 777 /**                        First Order Moments                         **/
 778 /**                                                                    **/
 779 /************************************************************************/
 780
 781 moms1front(gf, n, m, x, hessdata, mean, stdev, sigmadata, h, typx)
 782      char *gf;
 783      int *n, *m;
 784      double *x, *hessdata, *mean, *stdev, *sigmadata, *h, *typx;
 785 {
 786   int i, j, k;
 787   RMatrix hess, sigma, delg;
 788   double *delgdata, maxadd;
 789
 790   hess = (RMatrix) S_alloc(*n, sizeof(double *));
 791   sigma = (RMatrix) S_alloc(*m, sizeof(double *));
 792   delg = (RMatrix) S_alloc(*m, sizeof(double *));
 793   delgdata = (double *) S_alloc(*m * *n, sizeof(double));
 794
 795   for (i = 0; i < *n; i++) hess[i] = hessdata + i * *n;
 796   for (i = 0; i < *m; i++) sigma[i] = sigmadata + i * *m;
 797   for (i = 0; i < *m; i++) delg[i] = delgdata + i * *n;
 798
 799   gevalfront(gf, n, m, x, h, typx, mean, delgdata);
 800
 801   /* get the cholesky decomposition L of the hessian */
 802   choldecomp(hess, *n, 0.0, &maxadd);
 803
 804   /* forward solve for (L^-1) delg^T */
 805   for (k = 0; k < *m; k++) {
 806     for (i = 0; i < *n; i++) {
 807       if (hess[i][i] != 0.0) delg[k][i] /= hess[i][i];
 808       for (j = i + 1; j < *n; j++) delg[k][j] -= hess[j][i] * delg[k][i];
 809     }
 810   }
 811
 812   /* compute sigma and stdev */
 813   for (i = 0; i < *m; i++) {
 814     for (j = i; j < *m; j++) {
 815       for (k = 0, sigma[i][j] = 0.0; k < *n; k++)
 816         sigma[i][j] += delg[i][k] * delg[j][k];
 817       sigma[j][i] = sigma[i][j];
 818     }
 819     stdev[i] = sqrt(sigma[i][i]);
 820   }
 821 }
 822
 823 /************************************************************************/
 824 /**                                                                    **/
 825 /**                        Second Order Moments                        **/
 826 /**                                                                    **/
 827 /************************************************************************/
 828
 829 typedef struct {
 830   MaxIPars max;
 831   int full, covar;
 832 } MomIPars;
 833
 834 typedef struct {
 835   MaxDPars max;
 836   double mgfdel;
 837 } MomDPars;
 838
 839 moms2front(ff, gf, gnum, x, typx, fvals, gvals, ipars, dpars,
 840            mean, stdev, sigmadata)
 841      char **ff, **gf;
 842      int *gnum;
 843      double *x, *typx, *fvals, *gvals, *mean, *stdev, *sigmadata;
 844      MomIPars *ipars;
 845      MomDPars *dpars;
 846 {
 847   char *msg;
 848   double h, loglap0, loglap1, loglap2;
 849   double *tilts, *fvals1, *gvals1, *x1;
 850   MomDPars dp1, *dpars1 = &dp1;
 851   MomIPars ip1, *ipars1 = &ip1;
 852   int i, n, m;
 853
 854   n = ipars->max.n;
 855   m = *gnum;
 856   h = dpars->max.h;
 857
 858   tilts = (double *) S_alloc(m, sizeof(double));
 859   x1 = (double *) S_alloc(n, sizeof(double));
 860   fvals1 = (double *) S_alloc(1 + n + n * n, sizeof(double));
 861   gvals1 = (double *) S_alloc(m + n * m, sizeof(double));
 862
 863   maxfront(ff, nil, nil, x, typx, fvals, gvals, nil, nil,
 864            ipars, dpars, nil, &msg);
 865   copypars(x, fvals, gvals, ipars, dpars, x1, fvals1, gvals1, ipars1, dpars1);
 866   loglapfront(fvals1, nil, ipars1, dpars1, &loglap0);
 867   copypars(x, fvals, gvals, ipars, dpars, x1, fvals1, gvals1, ipars1, dpars1);
 868   moms1front(gf, &n, &m, x1, fvals1 + n + 1, mean, stdev, sigmadata, &h, typx);
 869
 870   if (ipars->full) {
 871     for (i = 0; i < m; i++) {
 872       copypars(x, fvals, gvals, ipars, dpars,
 873                x1, fvals1, gvals1, ipars1, dpars1);
 874       ipars1->max.m = 1;
 875       ipars1->max.exptilt = FALSE;
 876       dpars1->max.newtilt = 1.0;
 877       maxfront(ff, gf + i, nil, x1, typx, fvals1, gvals1, nil, nil,
 878                ipars1, dpars1, nil, &msg);
 879       loglapfront(fvals1, nil, ipars1, dpars1, &loglap1);
 880       loglap1 -= loglap0;
 881
 882       copypars(x, fvals, gvals, ipars, dpars,
 883                x1, fvals1, gvals1, ipars1, dpars1);
 884       ipars1->max.m = 1;
 885       ipars1->max.exptilt = FALSE;
 886       dpars1->max.newtilt = 2.0;
 887       maxfront(ff, gf + i, nil, x1, typx, fvals1, gvals1, nil, nil,
 888                ipars1, dpars1, nil, &msg);
 889       loglapfront(fvals1, nil, ipars1, dpars1, &loglap2);
 890       loglap2 -= loglap0;
 891
 892       mean[i] = exp(loglap1);
 893       stdev[i] = sqrt(exp(loglap2) - exp(2.0 * loglap1));
 894       if (ipars->covar) sigmadata[i + i * m] = stdev[i] * stdev[i];
 895     }
 896     if (ipars->covar) {
 897       char *cgf[2];
 898       int j;
 899
 900       for (i = 0; i < m; i++) {
 901         for (j = i + 1; j < m; j++) {
 902           cgf[0] = gf[i];
 903           cgf[1] = gf[j];
 904           copypars(x, fvals, gvals, ipars, dpars,
 905                    x1, fvals1, gvals1, ipars1, dpars1);
 906           ipars1->max.m = 2;
 907           ipars1->max.exptilt = FALSE;
 908           dpars1->max.newtilt = 1.0;
 909           maxfront(ff, gf + i, nil, x1, typx, fvals1, gvals1, nil, nil,
 910                    ipars1, dpars1, nil, &msg);
 911           loglapfront(fvals1, nil, ipars1, dpars1, &loglap1);
 912           loglap1 -= loglap0;
 913
 914           sigmadata[i + j * m] = exp(loglap1) - mean[i] * mean[j];
 915           sigmadata[j + i * m] = sigmadata[i + j * m];
 916         }
 917       }
 918     }
 919   }
 920   else {
 921     for (i = 0; i < m; i++)
 922       tilts[i] = (stdev[i] > 0.0) ? dpars->mgfdel / stdev[i] : dpars->mgfdel;
 923
 924     for (i = 0; i < m; i++) {
 925       copypars(x, fvals, gvals, ipars, dpars,
 926                x1, fvals1, gvals1, ipars1, dpars1);
 927       ipars1->max.m = 1;
 928       ipars1->max.exptilt = TRUE;
 929       dpars1->max.newtilt = tilts[i];
 930       maxfront(ff, gf + i, nil, x1, typx, fvals1, gvals1, nil, nil,
 931                ipars1, dpars1, nil, &msg);
 932       loglapfront(fvals1, nil, ipars1, dpars1, &loglap1);
 933       loglap1 -= loglap0;
 934
 935       copypars(x, fvals, gvals, ipars, dpars,
 936                x1, fvals1, gvals1, ipars1, dpars1);
 937       ipars1->max.m = 1;
 938       ipars1->max.exptilt = TRUE;
 939       dpars1->max.newtilt = -tilts[i];
 940       maxfront(ff, gf + i, nil, x1, typx, fvals1, gvals1, nil, nil,
 941                ipars1, dpars1, nil, &msg);
 942       loglapfront(fvals1, nil, ipars1, dpars1, &loglap2);
 943       loglap2 -= loglap0;
 944
 945       mean[i] = (loglap1 - loglap2) / (2.0 * tilts[i]);
 946       stdev[i] = sqrt((loglap1 + loglap2) / (tilts[i] * tilts[i]));
 947       if (ipars->covar) sigmadata[i + i * m] = stdev[i] * stdev[i];
 948     }
 949     if (ipars->covar) {
 950       char *cgf[2];
 951       double ctilt, tscale[2];
 952       int j;
 953
 954       ctilt = dpars->mgfdel;
 955       for (i = 0; i < m; i++) {
 956         for (j = i + 1; j < m; j++) {
 957           cgf[0] = gf[i];
 958           cgf[1] = gf[j];
 959           tscale[0] = stdev[i] > 0 ? stdev[i] : 1.0;
 960           tscale[1] = stdev[j] > 0 ? stdev[j] : 1.0;
 961
 962           copypars(x, fvals, gvals, ipars, dpars,
 963                    x1, fvals1, gvals1, ipars1, dpars1);
 964           ipars1->max.m = 2;
 965           ipars1->max.exptilt = TRUE;
 966           dpars1->max.newtilt = ctilt;
 967           maxfront(ff, cgf, nil, x1, typx, fvals1, gvals1, nil, nil,
 968                    ipars1, dpars1, tscale, &msg);
 969           loglapfront(fvals1, nil, ipars1, dpars1, &loglap1);
 970           loglap1 -= loglap0;
 971
 972           copypars(x, fvals, gvals, ipars, dpars,
 973                    x1, fvals1, gvals1, ipars1, dpars1);
 974           ipars1->max.m = 2;
 975           ipars1->max.exptilt = TRUE;
 976           dpars1->max.newtilt = -ctilt;
 977           maxfront(ff, cgf, nil, x1, typx, fvals1, gvals1, nil, nil,
 978                    ipars1, dpars1, tscale, &msg);
 979           loglapfront(fvals1, nil, ipars1, dpars1, &loglap2);
 980           loglap2 -= loglap0;
 981
 982           sigmadata[i + j * m] = stdev[i] * stdev[j]
 983             * ((loglap2 + loglap1) /(2 * ctilt * ctilt) - 1.0);
 984           sigmadata[j + i * m] = sigmadata[i + j * m];
 985         }
 986       }
 987     }
 988   }
 989 }
 990
 991 static copypars(x, f, g, ip, dp, x1, f1, g1, ip1, dp1)
 992      double *x, *f, *g, *x1, *f1, *g1;
 993      MomIPars *ip, *ip1;
 994      MomDPars *dp, *dp1;
 995 {
 996   int i, n, m, nf, ng;
 997
 998   n = ip->max.n;
 999   m = ip->max.m;
1000   nf = 1 + n + n * n;
1001   ng = m + n * m;
1002
1003   for (i = 0; i < n; i++) x1[i] = x[i];
1004   for (i = 0; i < nf; i++) f1[i] = f[i];
1005   for (i = 0; i < ng; i++) g1[i] = g[i];
1006   *ip1 = *ip;
1007   *dp1 = *dp;
1008 }
1009
1010 /************************************************************************/
1011 /**                                                                    **/
1012 /**                          Laplace Margins                           **/
1013 /**                                                                    **/
1014 /************************************************************************/
1015
1016 lapmar1front(ff, gf, x, typx, fvals, ipars, dpars, xmar, ymar, nmar)
1017      char **ff, **gf;
1018      double *x, *typx, *fvals, *xmar, *ymar;
1019      MaxIPars *ipars;
1020      MaxDPars *dpars;
1021      int *nmar;
1022 {
1023   char *msg;
1024   int i, n, m, mindex;
1025   double h, loglap0, loglap1, xmode, stdev, sigmadata, ctarget[1];
1026   double *fvals1, *x1, *cvals, *cvals1, *fvals2, *x2, *cvals2;
1027   MaxDPars dp1, dp2, *dpars1 = &dp1, *dpars2 = &dp2;
1028   MaxIPars ip1, ip2, *ipars1 = &ip1, *ipars2 = &ip2;
1029
1030   n = ipars->n;
1031   m = 1;
1032   h = dpars->h;
1033
1034   x1 = (double *) S_alloc(n + 1, sizeof(double));
1035   fvals1 = (double *) S_alloc(1 + n + n * n, sizeof(double));
1036   cvals = (double *) S_alloc(m + n * m, sizeof(double));
1037   cvals1 = (double *) S_alloc(m + n * m, sizeof(double));
1038   x2 = (double *) S_alloc(n + 1, sizeof(double));
1039   fvals2 = (double *) S_alloc(1 + n + n * n, sizeof(double));
1040   cvals2 = (double *) S_alloc(m + n * m, sizeof(double));
1041
1042   maxfront(ff, nil, nil, x, typx, fvals, nil, nil, nil,
1043            ipars, dpars, nil, &msg);
1044   cpmarpars(x, fvals, cvals, ipars, dpars, x1, fvals1, cvals1, ipars1, dpars1);
1045   loglapfront(fvals1, nil, ipars1, dpars1, &loglap0);
1046   cpmarpars(x, fvals, cvals, ipars, dpars, x1, fvals1, cvals1, ipars1, dpars1);
1047   moms1front(gf, &n, &m, x1, fvals1 + n + 1, &xmode, &stdev, &sigmadata,
1048              &h, typx);
1049
1050   ipars->k = 1;
1051   ipars->vals_suppl = FALSE;
1052   ctarget[0] = xmode;
1053   maxfront(ff, nil, gf, x, typx, fvals, nil, cvals, ctarget,
1054            ipars, dpars, nil, &msg);
1055
1056   for (mindex = 0; mindex < *nmar && xmar[mindex] <= xmode; mindex++);
1057
1058   cpmarpars(x, fvals, cvals, ipars, dpars, x1, fvals1, cvals1, ipars1, dpars1);
1059   for (i = mindex; i >= 0; i--) {
1060     ctarget[0] = xmar[i];
1061     maxfront(ff, nil, gf, x1, typx, fvals1, nil, cvals1, ctarget,
1062              ipars1, dpars1, nil, &msg);
1063     cpmarpars(x1, fvals1, cvals1, ipars1, dpars1, x2,
1064               fvals2, cvals2, ipars2, dpars2);
1065     loglapfront(fvals2, cvals2, ipars2, dpars2, &loglap1);
1066     ymar[i] = exp(loglap1 - loglap0);
1067   }
1068   cpmarpars(x, fvals, cvals, ipars, dpars, x1, fvals1, cvals1, ipars1, dpars1);
1069   for (i = mindex + 1; i < *nmar; i++) {
1070     ctarget[0] = xmar[i];
1071     maxfront(ff, nil, gf, x1, typx, fvals1, nil, cvals1, ctarget,
1072              ipars1, dpars1, nil, &msg);
1073     cpmarpars(x1, fvals1, cvals1, ipars1, dpars1, x2,
1074               fvals2, cvals2, ipars2, dpars2);
1075     loglapfront(fvals2, cvals2, ipars2, dpars2, &loglap1);
1076     ymar[i] = exp(loglap1 - loglap0);
1077   }
1078 }
1079
1080 static cpmarpars(x, f, g, ip, dp, x1, f1, g1, ip1, dp1)
1081      double *x, *f, *g, *x1, *f1, *g1;
1082      MaxIPars *ip, *ip1;
1083      MaxDPars *dp, *dp1;
1084 {
1085   int i, n, k, nf, ng;
1086
1087   n = ip->n;
1088   k = ip->k;
1089   nf = 1 + n + n * n;
1090   ng = k + n * k;
1091
1092   for (i = 0; i < n; i++) x1[i] = x[i];
1093   for (i = 0; i < nf; i++) f1[i] = f[i];
1094   for (i = 0; i < ng; i++) g1[i] = g[i];
1095   *ip1 = *ip;
1096   *dp1 = *dp;
1097 }
1098 #endif /* SBAYES */
1099
1100 #ifdef TODO
1101 get hessian from gradiant for analytical gradiants
1102 avoid repeated derivative calls in mimimize.
1103 2d margins
1104 use pos. definiteness info in margins
1105 #endif TODO