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