lib/derivatives.c

   1 /* derivatives - for Bayes code 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 /*
   8 #ifdef __APPLE__
   9 #include <stdlib.h>
  10 #endif
  11 */
  12
  13 #include <stdlib.h>
  14
  15 # include "xmath.h"
  16 extern double macheps();
  17
  18 # define nil 0L
  19 # define FALSE 0
  20 # define TRUE 1
  21
  22 void
  23 numergrad(size_t n, double *x, double *grad, double *fsum,
  24           int ffun(double *, double *, double *, double **),
  25           double h, double *typx)
  26      /*     int (*ffun)();*/
  27 {
  28   size_t i;
  29   double old_xi, f1, f2, hi;
  30
  31   for (i = 0; i < n; i++) {
  32     old_xi = x[i];
  33     hi = (typx != nil) ? typx[i] * h : h;
  34     x[i] = old_xi + hi;
  35     (*ffun)(x, &f1, nil, nil);
  36     x[i] = old_xi - hi;
  37     (*ffun)(x, &f2, nil, nil);
  38     x[i] = old_xi;
  39     grad[i] = (f1 - f2) / (2.0 * hi);
  40     fsum[i] = f1 + f2;
  41   }
  42 }
  43
  44 void
  45 numerhess(size_t n, double *x, double **hess, double f, double *fsum,
  46           int *ffun(double *, double *, double *, double *),
  47           double h, double *typx)
  48 {
  49   size_t i, j;
  50   double old_xi, old_xj, f1, f2, hi, hj;
  51
  52   for (i = 0; i < n; i++) {
  53     hi = (typx != nil) ? typx[i] * h : h;
  54     hess[i][i] = (fsum[i] - 2 * f) / (hi * hi);
  55     for (j = i + 1; j < n; j++) {
  56       hj = (typx != nil) ? typx[j] * h : h;
  57       old_xi = x[i];
  58       old_xj = x[j];
  59       x[i] = old_xi + hi;
  60       x[j] = old_xj + hj;
  61       (*ffun)(x, &f1, nil, nil);
  62       x[i] = old_xi - hi;
  63       x[j] = old_xj - hj;
  64       (*ffun)(x, &f2, nil, nil);
  65       x[i] = old_xi;
  66       x[j] = old_xj;
  67       hess[i][j] = (2 * f + f1 + f2 - fsum[i] - fsum[j]) / (2.0 * hi * hj);
  68       hess[j][i] = hess[i][j];
  69     }
  70   }
  71 }