From 0cb8750e8e4f20da09c390ee34022a8dc9dd51e7 Mon Sep 17 00:00:00 2001 From: AJ Rossini Date: Sun, 18 Mar 2007 16:33:48 +0100 Subject: [PATCH] ansi-fication of C glue; .gitignore rejects so'd objects; more doc cleanup --- .gitignore | 1 + README | 6 +- lib/Makefile | 19 +- lib/betabase.c | 9 +- lib/bivnor.c | 17 +- lib/cbayes.c | 95 ++++++---- lib/cdists.c | 365 ++++++++++++++++++++------------------ lib/{mclglue.c => cffi-glue.c} | 342 ++++++++++++++++++------------------ lib/cfft.c | 307 ++++++++++++++++---------------- lib/cholesky.c | 10 +- lib/clinalg.c | 109 ++++++------ lib/complex.c | 102 ++++++----- lib/complex.h | 14 +- lib/derivatives.c | 22 +-- lib/exclglue.c | 385 ---------------------------------------- lib/functions.c | 390 ++++++++++++++++++++++------------------- lib/gammabase.c | 11 +- lib/kernel.c | 16 +- lib/linalg.h | 3 +- lib/linalgdata.c | 49 +++--- lib/lowess.c | 154 ++++++++-------- lib/ludecomp.c | 14 +- lib/makerotation.c | 7 +- lib/minimize.c | 284 ++++++++++++++++-------------- lib/nor.c | 66 ++++--- lib/qrdecomp.c | 46 ++--- lib/rcondest.c | 64 +++---- lib/splines.c | 18 +- lib/studentbase.c | 27 +-- lib/svdecomp.c | 17 +- 30 files changed, 1366 insertions(+), 1603 deletions(-) rename lib/{mclglue.c => cffi-glue.c} (53%) delete mode 100644 lib/exclglue.c diff --git a/.gitignore b/.gitignore index dace02a..30c2d42 100644 --- a/.gitignore +++ b/.gitignore @@ -2,3 +2,4 @@ *.fasl semantic.cache fasl +*.so diff --git a/README b/README index 8e320b4..7f6051a 100644 --- a/README +++ b/README @@ -8,7 +8,11 @@ proceed. First, we update CLS. Then, continue with TSL. Separation is needed, so that we can make everything separate. We are working with SBCL and CLISP implementations to keep a wide -range of target platforms. +range of target platforms, as well as to challenge the system by the +two more common CL implementations (compiled, and byte-compiled). + +This is done to provide a point of comparison as well as a means for +building upon. diff --git a/lib/Makefile b/lib/Makefile index a0b3b8e..c0062bf 100644 --- a/lib/Makefile +++ b/lib/Makefile @@ -1,3 +1,5 @@ +LISPSTAT-SO = liblispstat.so + CDISTOBJS = betabase.o bivnor.o cdists.o gamln.o gammabase.o nor.o ppnd.o \ studentbase.o @@ -7,18 +9,27 @@ CLINALGOBJS = cholesky.o linalgdata.o ludecomp.o complex.o svdecomp.o \ CBAYESOBJS = cbayes.o functions.o derivatives.o minimize.o -CFLAGS = -O -G 0 -DINTPTR +CGLUE = cffi-glue.o + +CFLAGS = -g -DINTPTR -fPIC -Wall + +cffi : $(LISPSTAT-SO) + +sbcl : clib.a kcl: clib.a excl: clib.a exclglue.o -clib.a: ${CDISTOBJS} ${CLINALGOBJS} ${CBAYESOBJS} - ar rcv clib.a ${CDISTOBJS} ${CLINALGOBJS} ${CBAYESOBJS} +clib.a: ${CDISTOBJS} ${CLINALGOBJS} ${CBAYESOBJS} $(CGLUE) + ar rcv clib.a ${CDISTOBJS} ${CLINALGOBJS} ${CBAYESOBJS} $(CGLUE) ranlib clib.a +$(LISPSTAT-SO) :${CDISTOBJS} ${CLINALGOBJS} ${CBAYESOBJS} $(CGLUE) + gcc -shared -o $(LISPSTAT-SO) ${CDISTOBJS} ${CLINALGOBJS} ${CBAYESOBJS} $(CGLUE) + clean: - rm -f ${CDISTOBJS} ${CLINALGOBJS} ${CBAYESOBJS} *~ + rm -f ${CDISTOBJS} ${CLINALGOBJS} ${CBAYESOBJS} $(CGLUE) *~ cleanall: make clean diff --git a/lib/betabase.c b/lib/betabase.c index 27193cd..d7b29e8 100644 --- a/lib/betabase.c +++ b/lib/betabase.c @@ -9,16 +9,13 @@ extern double macheps(), gamma(); /* forward declarations */ static double logbeta(), betai(), xinbta(); -betabase(x, a, b, gia, gib, cdf) - int *gia, *gib; - double *x, *a, *b, *cdf; +void +betabase(double *x, double *a, double *b, int *gia, int *gib, double *cdf) { *cdf = betai(*x, *a, *b); } -double ppbeta(p, a, b, ifault) - double p, a, b; - int *ifault; +double ppbeta(double p, double a, double b, int *ifault) { double lbeta; diff --git a/lib/bivnor.c b/lib/bivnor.c index c0c3ed8..d10f573 100644 --- a/lib/bivnor.c +++ b/lib/bivnor.c @@ -3,8 +3,9 @@ #define twopi 6.283195307179587 #define con (twopi / 2.0) * 10.0e-10 -double bivnor(ah, ak, r) - double ah, ak, r; +extern void normbase(double *, double *); + +double bivnor(double ah, double ak, double r) { /* based on alg 4628 comm. acm oct 73 @@ -64,26 +65,28 @@ double bivnor(ah, ak, r) sgn = -1; t = 0; if (wk!=0) { - if (fabs(wk)>=1) + if (fabs(wk)>=1) { if (fabs(wk)==1) { t = wk*gw*(1-gw)/2; goto label40; - } - else { + } else { sgn = -sgn; wh = wh*wk; normbase(&wh, &g2); wk = 1/wk; - if (wk<0) + if (wk<0) { b = b+.5; + } b = b-(gw+g2)/2+gw*g2; } + } h2 = wh*wh; a2 = wk*wk; h4 = h2*.5; ex = 0; - if (h4<150.0) + if (h4<150.0) { ex = exp(-h4); + } w2 = h4*ex; ap = 1; s2 = ap-ex; diff --git a/lib/cbayes.c b/lib/cbayes.c index 40d1d6d..f6402ff 100644 --- a/lib/cbayes.c +++ b/lib/cbayes.c @@ -1,6 +1,7 @@ /* cbayes - Lisp interface to laplace approximation stuff */ /* Copyright (c) 1990, by Luke Tierney */ - + +#include /* for calloc/realloc */ #include "linalg.h" #ifdef INTPTR @@ -9,7 +10,12 @@ typedef int PTR; typedef char *PTR; #endif -extern char *calloc(), *realloc(); +extern void maximize_callback(int, PTR, PTR, PTR, PTR, PTR); +extern void evalfront(char **, int *, double *, double *, double *, + double *, double *, double *); + +extern void maxfront(); +extern void bufputstr(char *); /************************************************************************/ /** **/ @@ -46,9 +52,10 @@ typedef struct { /** **/ /************************************************************************/ -static meminit() +static void +meminit() { - static inited = FALSE; + static int inited = FALSE; int i; if (! inited) { @@ -59,54 +66,67 @@ static meminit() memcount = 0; } -static makespace(pptr, size) - char **pptr; - int size; +static int +makespace(void **pptr, int size) /* why are we using **char? */ +{ + if (size <= 0) { + return(1); /* we've done, by default, what we asked for. */ + } + if (*pptr == nil) { + *pptr = calloc(size, 1); + } else { + *pptr = realloc(*pptr, size); + } + if (size > 0 && *pptr == nil) { + return(0); /* xlfail("memory allocation failed"); FIXME:AJR xlfail redef. */ + } + return(1); +} + + +/************************************************************************/ +/** **/ +/** Callback Function **/ +/** **/ +/************************************************************************/ + +static void +callLminfun(int n, + double *x, double *fval, double *grad, double *hess, + int *derivs) { - if (size <= 0) return; - if (*pptr == nil) *pptr = calloc(size, 1); - else *pptr = realloc(*pptr, size); - if (size > 0 && *pptr == nil) xlfail("memory allocation failed"); + maximize_callback(n, (PTR) x, + (PTR) fval, (PTR) grad, (PTR) hess, (PTR) derivs); } -call_S(fun, narg, args, mode, length, names, nvals, values) - char *fun, **args, **mode, **names, **values; - long narg, nvals, *length; +void +call_S(char *fun, long narg, char **args, char **mode, long *length,char **names, + long nvals, char **values) { int n = length[0], derivs; static double *fval = nil, *grad = nil, *hess = nil; - makespace(&fval, 1 * sizeof(double)); - makespace(&grad, n * sizeof(double)); - makespace(&hess, n * n * sizeof(double)); + makespace((void **)&fval, 1 * sizeof(double)); /* probably should test the + result of this and the next + 2 to make sure that they are + appropriate */ + makespace((void **)&grad, n * sizeof(double)); + makespace((void **)&hess, n * n * sizeof(double)); - callLminfun(n, args[0], fval, grad, hess, &derivs); + callLminfun(n,(double *)args[0], fval, grad, hess, &derivs); values[0] = (char *) fval; values[1] = (derivs > 0) ? (char *) grad : nil; values[2] = (derivs > 1) ? (char *) hess : nil; } +void Recover(s, w) char *s, *w; { - xlfail(s); + /* FIXME:AJR: xlfail(s); */ + return; } - -/************************************************************************/ -/** **/ -/** Callback Function **/ -/** **/ -/************************************************************************/ - -static callLminfun(n, x, fval, grad, hess, derivs) - int n, *derivs; - RVector x, grad, hess; - double *fval; -{ - maximize_callback(n, (PTR) x, - (PTR) fval, (PTR) grad, (PTR) hess, (PTR) derivs); -} /************************************************************************/ /** **/ @@ -114,6 +134,7 @@ static callLminfun(n, x, fval, grad, hess, derivs) /** **/ /************************************************************************/ +void numgrad_front(n, px, pgrad, h, pscale) int n; PTR px, pgrad, pscale; @@ -122,10 +143,11 @@ numgrad_front(n, px, pgrad, h, pscale) LVAL f = nil; double fval; - evalfront(&f, &n, (double *) px, + evalfront((char **)&f, &n, (double *) px, &fval, (double *) pgrad, nil, &h, (double *) pscale); } +void numhess_front(n, px, pf, pgrad, phess, h, pscale) int n; PTR px, pf, pgrad, phess, pscale; @@ -133,7 +155,7 @@ numhess_front(n, px, pf, pgrad, phess, h, pscale) { LVAL f = nil; - evalfront(&f, &n, (double *) px, + evalfront((char **)&f, &n, (double *) px, (double *) pf, (double *) pgrad, (double *) phess, &h, (double *) pscale); } @@ -196,6 +218,7 @@ static MaxDPars getMaxDPars(dpars) return(dp); } +void minfo_maximize(px, pfvals, pscale, pip, pdp, verbose) PTR px, pfvals, pscale, pip, pdp; int verbose; diff --git a/lib/cdists.c b/lib/cdists.c index a79f11a..582b950 100644 --- a/lib/cdists.c +++ b/lib/cdists.c @@ -1,5 +1,14 @@ #include "xmath.h" +extern void normbase(double *,double *); +extern void gammabase(double *,double *, double *); +extern void studentbase(double *,double *, double *); +extern void betabase(double *, double *,double *,int *, int *,double *); +extern double max(double, double); +extern double min(double, double); + +extern void xlfail(char *); + extern double ppnd(), gamma(), bivnor(), uni(), ppgamma(), ppbeta(), ppstudent(); @@ -8,7 +17,7 @@ extern double tdens(); #ifndef PI #define PI 3.14159265358979323846 -#endif PI +#endif /* PI*/ #define TRUE 1 #define FALSE 0 @@ -36,57 +45,61 @@ double floor(x) #endif #endif -static checkflag(flag) - int flag; +static void +checkflag(int flag) { /* do nothing for now */ } -static checkexp(a) - double a; +static void +checkexp(double a) { - if (a <= 0.0) xlfail("non-positive gamma or beta exponent"); + if (a <= 0.0) { + xlfail("non-positive gamma or beta exponent"); + } } -static checkdf(df) - double df; +static void +checkdf(double df) { - if (df <= 0.0) xlfail("non-positive degrees of freedom"); + if (df <= 0.0) { + xlfail("non-positive degrees of freedom"); + } } -static checkprob(p, zerostrict, onestrict) - double p; - int zerostrict, onestrict; +static void +checkprob(double p, int zerostrict, int onestrict) { if (zerostrict) { if (p <= 0.0) xlfail("non-positive probability argument"); - } - else { + } else { if (p < 0.0) xlfail("negative probability argument"); } if (onestrict) { if (p >= 1.0) xlfail("probability argument not less than one"); - } - else { + } else { if (p > 1.0) xlfail("probability argument greater than one"); } } -static checkrho(r) - double r; +static void +checkrho(double r) { - if (r < -1 || r > 1) xlfail("correlation out of range"); + if (r < -1 || r > 1) { + xlfail("correlation out of range"); + } } -static checkpoisson(L) - double L; +static void +checkpoisson(double L) { - if (L < 0.0) xlfail("negative Poisson mean"); + if (L < 0.0) { + xlfail("negative Poisson mean"); + } } -static checkbinomial(n, p) - int n; - double p; +static void +checkbinomial(int n, double p) { if (p < 0.0 || p > 1.0) xlfail("binomial p out of range"); if (n < 1) xlfail("non-positive binomial n"); @@ -182,30 +195,28 @@ double bnormcdf(x, y, r) /*****************************************************************************/ /*****************************************************************************/ -/* Cauchy cdf */ -double cauchycdf(dx) - double dx; +double +cauchycdf(double dx) { return((atan(dx) + PI / 2) / PI); } -/* Cauchy quantile function */ -double cauchyquant(p) - double p; +double +cauchyquant(double p) { checkprob(p, TRUE, TRUE); return(tan(PI * (p - 0.5))); } -/* cauchy density */ -double cauchydens(dx) - double dx; +double +cauchydens(double dx) { return(tdens(dx, 1.0)); } /* cauchy generator */ -double cauchyrand() +double +cauchyrand() { double u1, u2, v1, v2; @@ -220,16 +231,13 @@ double cauchyrand() } /*****************************************************************************/ -/*****************************************************************************/ /** **/ /** Gamma Distribution **/ /** **/ /*****************************************************************************/ -/*****************************************************************************/ -/* gamma cdf */ -double gammacdf(x, a) - double x, a; +double +gammacdf(double x, double a) { double p; @@ -239,8 +247,8 @@ double gammacdf(x, a) return(p); } -double gammaquant(p, a) - double p, a; +double +gammaquant(double p, double a) { int flag; double x; @@ -253,20 +261,22 @@ double gammaquant(p, a) } /* gamma density */ -double gammadens(x, a) - double x, a; +double gammadens(double x, double a) { double dens; checkexp(a); - if (x <= 0.0) dens = 0.0; - else dens = exp(log(x) * (a - 1) - x - gamma(a)); + if (x <= 0.0) { + dens = 0.0; + } else { + dens = exp(log(x) * (a - 1) - x - gamma(a)); + } return(dens); } /* gamma generator */ -double gammarand(a) - double a; +double +gammarand(double a) { double x, u0, u1, u2, v, w, c, c1, c2, c3, c4, c5; static double e = -1.0; @@ -315,27 +325,28 @@ double gammarand(a) } /*****************************************************************************/ -/*****************************************************************************/ /** **/ /** Chi-Square Distribution **/ /** **/ /*****************************************************************************/ -/*****************************************************************************/ -double chisqcdf(x, df) - double x, df; +double chisqcdf(double x, double df) { double p, a; checkdf(df); - a = 0.5 * df; x = 0.5 * x; - if (x <= 0.0) p = 0.0; - else gammabase(&x, &a, &p); + a = 0.5 * df; + x = 0.5 * x; + if (x <= 0.0) { + p = 0.0; + } else { + gammabase(&x, &a, &p); + } return(p); } -double chisqquant(p, df) - double p, df; +double +chisqquant(double p, double df) { double x, a; int flag; @@ -348,9 +359,8 @@ double chisqquant(p, df) return(x); } -/* chi-square density */ -double chisqdens(dx, da) - double dx, da; +double +chisqdens(double dx, double da) { checkdf(da); da = 0.5 * da; @@ -358,51 +368,55 @@ double chisqdens(dx, da) return(0.5 * gammadens(dx, da)); } -/* chi-square generator */ -double chisqrand(df) - double df; +double +chisqrand(double df) { checkdf(df); return(2.0 * gammarand(df / 2.0)); } /*****************************************************************************/ -/*****************************************************************************/ /** **/ /** Beta Distribution **/ /** **/ /*****************************************************************************/ -/*****************************************************************************/ -double betacdf(x, a, b) - double x, a, b; +double +betacdf(double x, double a, double b) { double p; int ia, ib; checkexp(a); checkexp(b); ia = a; ib = b; - if (x <= 0.0) p = 0.0; - else if (x >= 1.0) p = 1.0; - else betabase(&x, &a, &b, &ia, &ib, &p); + if (x <= 0.0) { + p = 0.0; + } else { + if (x >= 1.0) { + p = 1.0; + } else { + betabase(&x, &a, &b, &ia, &ib, &p); + } + } return(p); } -double betaquant(p, a, b) - double p, a, b; +double +betaquant(double p, double a, double b) { double x; int flag; - checkexp(a); checkexp(b); + checkexp(a); + checkexp(b); checkprob(p, FALSE, FALSE); x = ppbeta(p, a, b, &flag); checkflag(flag); return(x); } -static double logbeta(a, b) - double a, b; +static double +logbeta(double a, double b) { static double da = 0.0, db = 0.0, lbeta = 0.0; @@ -414,21 +428,24 @@ static double logbeta(a, b) } /* beta density */ -double betadens(x, a, b) - double x, a, b; +double +betadens(double x, double a, double b) { double dens; checkexp(a); checkexp(b); - if (x <= 0.0 || x >= 1.0) dens = 0.0; - else dens = exp(log(x) * (a - 1) + log(1 - x) * (b - 1) - logbeta(a, b)); + if (x <= 0.0 || x >= 1.0) { + dens = 0.0; + } else { + dens = exp(log(x) * (a - 1) + log(1 - x) * (b - 1) - logbeta(a, b)); + } return(dens); } /* beta generator */ -double betarand(a, b) - double a, b; +double +betarand(double a, double b) { double x, y; @@ -440,16 +457,14 @@ double betarand(a, b) } /*****************************************************************************/ -/*****************************************************************************/ /** **/ /** t Distribution **/ /** **/ /*****************************************************************************/ -/*****************************************************************************/ /* t cdf */ -double tcdf(x, df) - double x, df; +double +tcdf(double x, double df) { double p; @@ -458,9 +473,8 @@ double tcdf(x, df) return(p); } -/* t quantile function */ -double tquant(p, df) - double p, df; +double +tquant(double p, double df) { double x; int flag; @@ -472,9 +486,8 @@ double tquant(p, df) return(x); } -/* t density */ -double tdens(x, a) - double x, a; +double +tdens(double x, double a) { double dens; @@ -485,42 +498,35 @@ double tdens(x, a) return(dens); } -/* t generator */ -double trand(df) - double df; +double trand(double df) { checkdf(df); return(normalrand() / sqrt(chisqrand(df) / df)); } /*****************************************************************************/ -/*****************************************************************************/ /** **/ /** F Distribution **/ /** **/ /*****************************************************************************/ -/*****************************************************************************/ -/* f cdf */ -double fcdf(x, ndf, ddf) - double x, ndf, ddf; +double fcdf(double x, double ndf, double ddf) { double p, a, b; checkdf(ndf); checkdf(ddf); a = 0.5 * ddf; b = 0.5 * ndf; - if (x <= 0.0) p = 0.0; - else { + if (x <= 0.0) { + p = 0.0; + } else { x = a / (a + b * x); p = 1.0 - betacdf(x, a, b); } return(p); } -/* f quantile function */ -double fquant(p, ndf, ddf) - double p, ndf, ddf; +double fquant(double p, double ndf, double ddf) { double x, a, b; int flag; @@ -529,8 +535,9 @@ double fquant(p, ndf, ddf) checkprob(p, FALSE, TRUE); a = 0.5 * ddf; b = 0.5 * ndf; - if (p == 0.0) x = 0.0; - else { + if (p == 0.0) { + x = 0.0; + } else { p = 1.0 - p; x = ppbeta(p, a, b, &flag); checkflag(flag); @@ -539,25 +546,26 @@ double fquant(p, ndf, ddf) return(x); } -/* f density */ -double fdens(dx, da, db) - double dx, da, db; +double +fdens(double dx, double da, double db) { double dens; checkdf(da); checkdf(db); - if (dx <= 0.0) dens = 0.0; - else dens = exp(0.5 * da * log(da) + 0.5 * db *log(db) - + (0.5 * da - 1.0) * log(dx) - - logbeta(0.5 * da, 0.5 * db) - - 0.5 * (da + db) * log(db + da * dx)); + if (dx <= 0.0) { + dens = 0.0; + } else { + dens = exp(0.5 * da * log(da) + 0.5 * db *log(db) + + (0.5 * da - 1.0) * log(dx) + - logbeta(0.5 * da, 0.5 * db) + - 0.5 * (da + db) * log(db + da * dx)); + } return(dens); } /* f generator */ -double frand(ndf, ddf) - double ndf, ddf; +double frand(double ndf, double ddf) { checkdf(ndf); checkdf(ddf); @@ -565,39 +573,39 @@ double frand(ndf, ddf) } /*****************************************************************************/ -/*****************************************************************************/ /** **/ /** Poisson Distribution **/ /** **/ /*****************************************************************************/ -/*****************************************************************************/ -static double poisson_cdf(k, L) - int k; - double L; +static double +poisson_cdf(int k, double L) { double dp, dx; - if (k < 0) dp = 0.0; - else if (L == 0.0) dp = (k < 0) ? 0.0 : 1.0; - else { - dx = k + 1.0; - gammabase(&L, &dx, &dp); - dp = 1.0 - dp; + if (k < 0) { + dp = 0.0; + } else { + if (L == 0.0) { + dp = (k < 0) ? 0.0 : 1.0; + } else { + dx = k + 1.0; + gammabase(&L, &dx, &dp); + dp = 1.0 - dp; + } } return(dp); } -double poissoncdf(k, L) - double k; - double L; +double +poissoncdf(double k, double L) { checkpoisson(L); return(poisson_cdf((int) floor(k), L)); } -int poissonquant(x, L) - double x, L; +int +poissonquant(double x, double L) { int k, k1, k2, del, ia; double m, s, p1, p2, pk; @@ -608,15 +616,22 @@ int poissonquant(x, L) s = sqrt(L); del = max(1, (int) (0.2 * s)); - if (x == 0.0) k = 0.0; - else k = m + s * ppnd(x, &ia); - k1 = k; k2 = k; + if (x == 0.0) { + k = 0.0; + } else { + k = m + s * ppnd(x, &ia); + } + k1 = k; + k2 = k; do { k1 = k1 - del; k1 = max(0, k1); p1 = poisson_cdf(k1, L); } while (k1 > 0 && p1 > x); - if (k1 == 0 && p1 >= x) return(k1); + + if (k1 == 0 && p1 >= x) { + return(k1); + } do { k2 = k2 + del; @@ -626,29 +641,36 @@ int poissonquant(x, L) while (k2 - k1 > 1) { k = (k1 + k2) / 2; pk = poisson_cdf(k, L); - if (pk < x) { k1 = k; p1 = pk; } - else { k2 = k; p2 = pk; } + if (pk < x) { + k1 = k; p1 = pk; + } else { + k2 = k; p2 = pk; + } } return(k2); } -double poissonpmf(k, L) - int k; - double L; +double +poissonpmf(int k, double L) { double dx, dp; checkpoisson(L); dx = k; - if (L == 0.0) dp = (k == 0) ? 1.0 : 0.0; - else if (dx < 0.0) dp = 0.0; - else dp = exp(dx * log(L) - L - gamma(dx + 1.0)); + if (L == 0.0) { + dp = (k == 0) ? 1.0 : 0.0; + } else { + if (dx < 0.0) { + dp = 0.0; + } else { + dp = exp(dx * log(L) - L - gamma(dx + 1.0)); + } + } return(dp); } /* poisson random generator from Numerical Recipes */ -int poissonrand(xm) - double xm; +int poissonrand(double xm) { static double sqrt2xm, logxm, expxm, g, oldxm = -1.0; double t, y; @@ -663,8 +685,7 @@ int poissonrand(xm) k++; t *= uni(); } while (t > expxm); - } - else { + } else { if (xm != oldxm) { oldxm = xm; sqrt2xm = sqrt(2.0 * xm); @@ -683,46 +704,52 @@ int poissonrand(xm) } /*****************************************************************************/ -/*****************************************************************************/ /** **/ /** Binomial Distribution **/ /** **/ /*****************************************************************************/ -/*****************************************************************************/ -static double binomial_cdf(k, n, p) - int k, n; - double p; +static double +binomial_cdf(int k, int n, double p) { double da, db, dp; int ia, ib; - if (k < 0) dp = 0.0; - else if (k >= n) dp = 1.0; - else if (p == 0.0) dp = (k < 0) ? 0.0 : 1.0; - else if (p == 1.0) dp = (k < n) ? 0.0 : 1.0; - else { - da = k + 1.0; - db = n - k; - ia = floor(da); ib = floor(db); - betabase(&p, &da, &db, &ia, &ib, &dp); - dp = 1.0 - dp; + if (k < 0) { + dp = 0.0; + } else { + if (k >= n) { + dp = 1.0; + } else { + if (p == 0.0) { + dp = (k < 0) ? 0.0 : 1.0; + } else { + if (p == 1.0) { + dp = (k < n) ? 0.0 : 1.0; + } else { + da = k + 1.0; + db = n - k; + ia = floor(da); + ib = floor(db); + betabase(&p, &da, &db, &ia, &ib, &dp); + dp = 1.0 - dp; + } + } + } } return(dp); } -double binomialcdf(k, n, p) - double k, p; - int n; +double +binomialcdf(double k, int n, double p) { checkbinomial(n, p); return(binomial_cdf((int) floor(k), n, p)); } -int binomialquant(x, n, p) - double x, p; - int n; +int +binomialquant(double x, int n, double p) { int k, k1, k2, del, ia; double m, s, p1, p2, pk; diff --git a/lib/mclglue.c b/lib/cffi-glue.c similarity index 53% rename from lib/mclglue.c rename to lib/cffi-glue.c index e3c9735..19cff4c 100644 --- a/lib/mclglue.c +++ b/lib/cffi-glue.c @@ -1,3 +1,6 @@ +#include +#include + #include "linalg.h" typedef char *PTR; @@ -16,31 +19,44 @@ extern int poissonquant(), poissonrand(); extern double binomialcdf(), binomialpmf(); extern int binomialquant(), binomialrand(); +void xlfail(char *); + +extern int chol_decomp_front(); +extern int lu_decomp_front(); +extern int lu_solve_front(); +extern int lu_inverse_front(); +extern int sv_decomp_front(); +extern int qr_decomp_front(); +extern int rcondest_front(); +extern int make_rotation_front(); +extern int eigen_front(); +extern int range_to_rseq(); +extern int spline_front(); +extern int kernel_dens_front(); +extern int kernel_smooth_front(); +extern int base_lowess_front(); +extern int fft_front(); +extern int numgrad_front(); +extern int numhess_front(); +extern int minfo_maximize(); + /***********************************************************************/ -/***********************************************************************/ -/**** ****/ /**** Basic Utilities ****/ -/**** ****/ -/***********************************************************************/ /***********************************************************************/ /***********************************************************************/ -/** **/ /** Callback Support Functions **/ -/** **/ /***********************************************************************/ static int ccl_integer_value; static double ccl_double_value; static PTR ccl_ptr_value; -ccl_store_integer(x) int x; { ccl_integer_value = x; } -ccl_store_double(x) double x; { ccl_double_value = x; } -ccl_store_ptr(x) PTR x; { ccl_ptr_value = x; } +void ccl_store_integer(int x) { ccl_integer_value = x; } +void ccl_store_double(double x) { ccl_double_value = x; } +void ccl_store_ptr(PTR x) { ccl_ptr_value = x; } /***************************************************************************/ -/** **/ /** Lisp-Managed Calloc/Free **/ -/** **/ /***************************************************************************/ #ifdef DODO @@ -70,7 +86,7 @@ void free(p) { (*free_ptr)(p); } -#endif DODO +#endif /* DODO*/ /***************************************************************************/ /** **/ @@ -78,23 +94,23 @@ void free(p) /** **/ /***************************************************************************/ -PTR la_base_allocate(n, m) - unsigned n, m; +PTR +la_base_allocate(unsigned n, unsigned m) { char *p = calloc(n, m); if (p == nil) xlfail("allocation failed"); return((PTR) p); } -int la_base_free_alloc(p) - PTR p; +int +la_base_free_alloc(PTR p) { if (p) free((char *) p); return(0); } -int la_mode_size(mode) - int mode; +int +la_mode_size(int mode) { switch (mode) { case IN: return(sizeof(int)); @@ -112,18 +128,18 @@ int la_mode_size(mode) int (*ccl_la_allocate)(), (*ccl_la_free_alloc)(); -register_la_allocate(f) int (*f)(); { ccl_la_allocate = f; } -register_la_free_alloc(f) int (*f)(); { ccl_la_free_alloc = f; } +void register_la_allocate(f) int (*f)(); { ccl_la_allocate = f; } +void register_la_free_alloc(f) int (*f)(); { ccl_la_free_alloc = f; } -PTR la_allocate(n, m) - int n, m; +PTR +la_allocate(int n, int m) { (*ccl_la_allocate)(n, m); return(ccl_ptr_value); } -la_free_alloc(p) - PTR p; +void +la_free_alloc(PTR p) { (*ccl_la_free_alloc)(p); } @@ -134,39 +150,34 @@ la_free_alloc(p) /** **/ /***************************************************************************/ -int la_get_integer(p, i) - PTR p; - int i; +int +la_get_integer(PTR p, int i) { return(*(((int *) p) + i)); } -double la_get_double(p, i) - PTR p; - int i; +double +la_get_double(PTR p, int i) { return(*(((double *) p) + i)); } -double la_get_complex_real(p, i) - PTR p; - int i; +double +la_get_complex_real(PTR p, int i) { Complex *c = ((Complex *) p) + i; return(c->real); } -double la_get_complex_imag(p, i) - PTR p; - int i; +double +la_get_complex_imag(PTR p, int i) { Complex *c = ((Complex *) p) + i; return(c->imag); } -PTR la_get_pointer(p, i) - PTR p; - int i; +PTR +la_get_pointer(PTR p, int i) { return(*(((PTR *) p) + i)); } @@ -177,27 +188,21 @@ PTR la_get_pointer(p, i) /** **/ /***************************************************************************/ -int la_put_integer(p, i, x) - PTR p; - int i, x; +int +la_put_integer(PTR p, int i, int x) { *(((int *) p) + i) = x; return(0); } -int la_put_double(p, i, x) - PTR p; - int i; - double x; +int la_put_double(PTR p, int i, double x) { *(((double *) p) + i) = x; return(0); } -int la_put_complex(p, i, x, y) - PTR p; - int i; - double x, y; +int +la_put_complex(PTR p, int i, double x, double y) { Complex *c = ((Complex *) p) + i; c->real = x; @@ -205,9 +210,8 @@ int la_put_complex(p, i, x, y) return(0); } -int la_put_pointer(p, i, x) - PTR p, x; - int i; +int +la_put_pointer(PTR p, int i, PTR x) { *(((PTR *) p) + i) = x; return(0); @@ -222,19 +226,23 @@ int la_put_pointer(p, i, x) char buf[1000]; static int (*ccl_set_buf_char_fptr)(); -register_set_buf_char(f) int (*f)(); { ccl_set_buf_char_fptr = f; } -set_buf_char(n, c) int n, c; { (*ccl_set_buf_char_fptr)(n, c); } +void register_set_buf_char(f) int (*f)(); { ccl_set_buf_char_fptr = f; } +void set_buf_char(int n, int c) { (*ccl_set_buf_char_fptr)(n, c); } static int (*ccl_print_buffer)(); -register_print_buffer(f) int (*f)(); { ccl_print_buffer = f; } -print_buffer(n, m) int n, m; { (*ccl_print_buffer)(n, m); } +void register_print_buffer(f) int (*f)(); { ccl_print_buffer = f; } +void print_buffer(int n, int m) { (*ccl_print_buffer)(n, m); } static int bufpos = 0; -static resetbuf() { bufpos = 0; } +static void +resetbuf() +{ + bufpos = 0; +} -static prbuf(s) - char *s; +static void +prbuf(char *s) { int i, n; @@ -243,174 +251,156 @@ static prbuf(s) set_buf_char(bufpos, 0); } -xlfail(s) - char *s; +void +xlfail(char *s) { resetbuf(); prbuf(s); print_buffer(bufpos, 1); } -stdputstr(s) - char *s; +void +stdputstr(char *s) { resetbuf(); prbuf(s); print_buffer(bufpos, 0); } -bufputstr(s) - char *s; +void +bufputstr(char *s) { resetbuf(); prbuf(s); } /***************************************************************************/ -/***************************************************************************/ /**** ****/ /***** Lisp Interfaces to Linear Algebra Routines ****/ /**** ****/ /***************************************************************************/ -/***************************************************************************/ -ccl_chol_decomp_front(mat, n, dpars) - PTR mat, dpars; - int n; +int +ccl_chol_decomp_front(PTR mat, int n, PTR dpars) { return(chol_decomp_front(mat, n, dpars)); } -ccl_lu_decomp_front(mat, n, iv, mode, dp) - PTR mat, iv, dp; - int n, mode; +int +ccl_lu_decomp_front(PTR mat, int n, PTR iv, int mode, PTR dp) { return(lu_decomp_front(mat, n, iv, mode, dp)); } -ccl_lu_solve_front(a, n, indx, b, mode) - PTR a, indx, b; - int n, mode; +int +ccl_lu_solve_front(PTR a, int n, PTR indx, PTR b, int mode) { return(lu_solve_front(a, n, indx, b, mode)); } -ccl_lu_inverse_front(pmat, n, piv, pv, mode, pinv) - PTR pmat, piv, pv, pinv; - int n, mode; +int +ccl_lu_inverse_front(PTR pmat, int n, PTR piv, PTR pv, int mode, PTR pinv) { return(lu_inverse_front(pmat, n, piv, pv, mode, pinv)); } -ccl_sv_decomp_front(mat, m, n, w, v) - PTR mat, w, v; - int m, n; +int +ccl_sv_decomp_front(PTR mat, int m, int n, PTR w, PTR v) { return(sv_decomp_front(mat, m, n, w, v)); } -ccl_qr_decomp_front(mat, m, n, v, jpvt, pivot) - PTR mat, v, jpvt; - int m, n, pivot; + +int +ccl_qr_decomp_front(PTR mat, int m, int n, PTR v, PTR jpvt, int pivot) { return(qr_decomp_front(mat, m, n, v, jpvt, pivot)); } -double ccl_rcondest_front(mat, n) - PTR mat; - int n; +double +ccl_rcondest_front(PTR mat, int n) { return(rcondest_front(mat, n)); } -ccl_make_rotation_front(n, rot, x, y, use_alpha, alpha) - int n, use_alpha; - PTR rot, x, y; - double alpha; +int +ccl_make_rotation_front(int n, PTR rot, PTR x, PTR y, int use_alpha, double alpha) { return(make_rotation_front(n, rot, x, y, use_alpha, alpha)); } -ccl_eigen_front(a, n, w, z, fv1) - PTR a, w, z, fv1; - int n; +int +ccl_eigen_front(PTR a, int n, PTR w, PTR z, PTR fv1) { return(eigen_front(a, n, w, z, fv1)); } -ccl_range_to_rseq(n, px, ns, pxs) - int n, ns; - PTR px, pxs; +int +ccl_range_to_rseq(int n, PTR px, int ns, PTR pxs) { return(range_to_rseq(n, px, ns, pxs)); } -ccl_spline_front(n, x, y, ns, xs, ys, work) - PTR x, y, xs, ys, work; - int n, ns; +int +ccl_spline_front(int n, PTR x, PTR y, int ns, PTR xs, PTR ys, PTR work) { return(spline_front(n, x, y, ns, xs, ys, work)); } -ccl_kernel_dens_front(x, n, width, xs, ys, ns, code) - PTR x, xs, ys; - int n, ns, code; - double width; +int +ccl_kernel_dens_front(PTR x, int n, double width, PTR xs, PTR ys, int ns, int code) { return(kernel_dens_front(x, n, width, xs, ys, ns, code)); } -ccl_kernel_smooth_front(x, y, n, width, xs, ys, ns, code) - PTR x, y, xs, ys; - int n, ns, code; - double width; +int +ccl_kernel_smooth_front(PTR x, PTR y, int n, double width, PTR xs, PTR ys, int ns, int code) { return(kernel_smooth_front(x, y, n, width, xs, ys, ns, code)); } -ccl_base_lowess_front(x, y, n, f, nsteps, delta, ys, rw, res) - PTR x, y, ys, rw, res; - int n, nsteps; - double f, delta; +int +ccl_base_lowess_front(PTR x, PTR y, int n, double f, + int nsteps, double delta, PTR ys, PTR rw, PTR res) { return(base_lowess_front(x, y, n, f, nsteps, delta, ys, rw, res)); } -ccl_fft_front(n, x, work, isign) - int n, isign; - PTR x, work; +int +ccl_fft_front(int n, PTR x, PTR work, int isign) { return(fft_front(n, x, work, isign)); } static int (*ccl_maximize_callback)(); -register_maximize_callback(f) int (*f)(); { ccl_maximize_callback = f; } -maximize_callback(n, px, pfval, pgrad, phess, pderivs) - int n; - PTR px, pfval, pgrad, phess, pderivs; + +void +register_maximize_callback(f) + int (*f)(); +{ + ccl_maximize_callback = f; +} + +void +maximize_callback(int n, PTR px, PTR pfval, PTR pgrad, PTR phess, PTR pderivs) { (*ccl_maximize_callback)(n, px, pfval, pgrad, phess, pderivs); } -ccl_numgrad_front(n, px, pgrad, h, pscale) - int n; - PTR px, pgrad, pscale; - double h; +int +ccl_numgrad_front(int n, PTR px, PTR pgrad, double h, PTR pscale) { return(numgrad_front(n, px, pgrad, h, pscale)); } -ccl_numhess_front(n, px, pf, pgrad, phess, h, pscale) - int n; - PTR px, pf, pgrad, phess, pscale; - double h; +int +ccl_numhess_front(int n, PTR px, PTR pf, PTR pgrad, PTR phess, double h, PTR pscale) { return(numhess_front(n, px, pf, pgrad, phess, h, pscale)); } -ccl_minfo_maximize(px, pfvals, pscale, pip, pdp, verbose) - PTR px, pfvals, pscale, pip, pdp; - int verbose; +int +ccl_minfo_maximize(PTR px, PTR pfvals, PTR pscale, PTR pip, PTR pdp, int verbose) { return(minfo_maximize(px, pfvals, pscale, pip, pdp, verbose)); } @@ -424,53 +414,55 @@ ccl_minfo_maximize(px, pfvals, pscale, pip, pdp, verbose) /***********************************************************************/ static int (*ccl_uni_fptr)(); -register_uni(f) int (*f)(); { ccl_uni_fptr = f; } +void register_uni(f) int (*f)(); { ccl_uni_fptr = f; } double uni() { (*ccl_uni_fptr)(); return(ccl_double_value); } -double ccl_gamma(x) double x; { return(gamma(x)); } +double ccl_gamma(double x) { return(gamma(x)); } -double ccl_normalcdf(x) double x; { return(normalcdf(x)); } -double ccl_normalquant(x) double x; { return(normalquant(x)); } -double ccl_normaldens(x) double x; { return(normaldens(x)); } +double ccl_normalcdf(double x) { return(normalcdf(x)); } +double ccl_normalquant(double x) { return(normalquant(x)); } +double ccl_normaldens(double x) { return(normaldens(x)); } double ccl_normalrand() { return(normalrand()); } -double ccl_bnormcdf(x, y, z) double x, y, z; { return(bnormcdf(x,y,z)); } +double ccl_bnormcdf(double x, double y, double z) { + return(bnormcdf(x,y,z)); +} -double ccl_cauchycdf(x) double x; { return(cauchycdf(x)); } -double ccl_cauchyquant(x) double x; { return(cauchyquant(x)); } -double ccl_cauchydens(x) double x; { return(cauchydens(x)); } +double ccl_cauchycdf(double x) { return(cauchycdf(x)); } +double ccl_cauchyquant(double x) { return(cauchyquant(x)); } +double ccl_cauchydens(double x) { return(cauchydens(x)); } double ccl_cauchyrand() { return(cauchyrand()); } -double ccl_gammacdf(x, y) double x, y; { return(gammacdf(x, y)); } -double ccl_gammaquant(x, y) double x, y; { return(gammaquant(x, y)); } -double ccl_gammadens(x, y) double x, y; { return(gammadens(x, y)); } -double ccl_gammarand(x) double x; { return(gammarand(x)); } - -double ccl_chisqcdf(x, y) double x, y; { return(chisqcdf(x, y)); } -double ccl_chisqquant(x, y) double x, y; { return(chisqquant(x, y)); } -double ccl_chisqdens(x, y) double x, y; { return(chisqdens(x, y)); } -double ccl_chisqrand(x) double x; { return(chisqrand(x)); } - -double ccl_betacdf(x, y, z) double x, y, z; { return(betacdf(x, y, z)); } -double ccl_betaquant(x, y, z) double x, y, z; { return(betaquant(x, y, z)); } -double ccl_betadens(x, y, z) double x, y, z; { return(betadens(x, y, z)); } -double ccl_betarand(x, y) double x, y; { return(betarand(x, y)); } - -double ccl_tcdf(x, y) double x, y; { return(tcdf(x, y)); } -double ccl_tquant(x, y) double x, y; { return(tquant(x, y)); } -double ccl_tdens(x, y) double x, y; { return(tdens(x, y)); } -double ccl_trand(x) double x; { return(trand(x)); } - -double ccl_fcdf(x, y, z) double x, y, z; { return(fcdf(x, y, z)); } -double ccl_fquant(x, y, z) double x, y, z; { return(fquant(x, y, z)); } -double ccl_fdens(x, y, z) double x, y, z; { return(fdens(x, y, z)); } -double ccl_frand(x, y) double x, y; { return(frand(x, y)); } - -double ccl_poissoncdf(x, y) double x, y; { return(poissoncdf(x, y)); } -int ccl_poissonquant(x, y) double x, y; { return(poissonquant(x, y)); } -double ccl_poissonpmf(x, y) int x; double y; { return(poissonpmf(x, y)); } -int ccl_poissonrand(x) double x; { return(poissonrand(x)); } - -double ccl_binomialcdf(x, y, z) double x, z; int y; { return(binomialcdf(x, y, z)); } -int ccl_binomialquant(x, y, z) double x, z; int y; { return(binomialquant(x, y, z)); } -double ccl_binomialpmf(x, y, z) int x, y; double z; { return(binomialpmf(x, y, z)); } -int ccl_binomialrand(x, y) int x; double y; { return(binomialrand(x, y)); } +double ccl_gammacdf(double x, double y) { return(gammacdf(x, y)); } +double ccl_gammaquant(double x, double y) { return(gammaquant(x, y)); } +double ccl_gammadens(double x, double y) { return(gammadens(x, y)); } +double ccl_gammarand(double x) { return(gammarand(x)); } + +double ccl_chisqcdf(double x, double y) { return(chisqcdf(x, y)); } +double ccl_chisqquant(double x, double y) { return(chisqquant(x, y)); } +double ccl_chisqdens(double x, double y) { return(chisqdens(x, y)); } +double ccl_chisqrand(double x) { return(chisqrand(x)); } + +double ccl_betacdf(double x, double y, double z) { return(betacdf(x, y, z)); } +double ccl_betaquant(double x, double y, double z) { return(betaquant(x, y, z)); } +double ccl_betadens(double x, double y, double z) { return(betadens(x, y, z)); } +double ccl_betarand(double x, double y) { return(betarand(x, y)); } + +double ccl_tcdf(double x, double y) { return(tcdf(x, y)); } +double ccl_tquant(double x, double y) { return(tquant(x, y)); } +double ccl_tdens(double x, double y) { return(tdens(x, y)); } +double ccl_trand(double x) { return(trand(x)); } + +double ccl_fcdf(double x, double y, double z) { return(fcdf(x, y, z)); } +double ccl_fquant(double x, double y, double z) { return(fquant(x, y, z)); } +double ccl_fdens(double x, double y, double z) { return(fdens(x, y, z)); } +double ccl_frand(double x, double y) { return(frand(x, y)); } + +double ccl_poissoncdf(double x, double y) { return(poissoncdf(x, y)); } +int ccl_poissonquant(double x, double y) { return(poissonquant(x, y)); } +double ccl_poissonpmf(int x, double y) { return(poissonpmf(x, y)); } +int ccl_poissonrand(double x) { return(poissonrand(x)); } + +double ccl_binomialcdf(double x, int y, double z) { return(binomialcdf(x, y, z)); } +int ccl_binomialquant(double x, int y, double z) { return(binomialquant(x, y, z)); } +double ccl_binomialpmf(int x, int y, double z) { return(binomialpmf(x, y, z)); } +int ccl_binomialrand(int x, double y) { return(binomialrand(x, y)); } diff --git a/lib/cfft.c b/lib/cfft.c index bbe818e..80d97a0 100644 --- a/lib/cfft.c +++ b/lib/cfft.c @@ -3,162 +3,9 @@ #include "xmath.h" /* - Public Routine -*/ - -/* - * cfft computes the forward or backward complex discrete fourier transform. - * - * Input Parameters: - * - * n The length of the complex sequence c. The method is - * more efficient when n is the product of small primes. - * - * c A complex array of length n which contains the sequence - * - * wsave a real work array which must be dimensioned at least 4n+15 - * in the program that calls cfft. - * isign 1 for transform, -1 for inverse transform. - * A call of cfft with isign = 1 followed by a call of cfft with - * isign = -1 will multiply the sequence by n. - * - * Output Parameters: - * - * c For j=1,...,n - * - * c(j)=the sum from k=1,...,n of - * - * c(k)*exp(-isign*i*(j-1)*(k-1)*2*pi/n) - * - * where i=sqrt(-1) - */ - -int cfft(n, c, wsave, isign) - int n; - double *c, *wsave; - int isign; -{ - int iw1, iw2; - - /* Parameter adjustments */ - --c; - --wsave; - - /* Function Body */ - if (n != 1) { - iw1 = n + n + 1; - iw2 = iw1 + n + n; - cffti1_(&n, &wsave[iw1], &wsave[iw2]); - cfft1_(&n, &c[1], &wsave[1], &wsave[iw1], &wsave[iw2], &isign); - } - return 0; -} - -/* Internal Routines */ -static int cfft1_(n, c, ch, wa, ifac, isign) - int *n; - double *c, *ch, *wa; - int *ifac; - int *isign; -{ - /* System generated locals */ - int i_1; - - /* Local variables */ - int idot; - int i, k1, l1, l2, n2, na, nf, ip, iw, ix2, ix3, ix4, nac, ido, idl1; - - /* Parameter adjustments */ - --c; - --ch; - --wa; - --ifac; - - /* Function Body */ - nf = ifac[2]; - na = 0; - l1 = 1; - iw = 1; - i_1 = nf; - for (k1 = 1; k1 <= i_1; ++k1) { - ip = ifac[k1 + 2]; - l2 = ip * l1; - ido = *n / l2; - idot = ido + ido; - idl1 = idot * l1; - if (ip == 4) { - ix2 = iw + idot; - ix3 = ix2 + idot; - if (na == 0) { - pass4_(&idot, &l1, &c[1], &ch[1], &wa[iw], &wa[ix2], &wa[ix3], isign); - } - else { - pass4_(&idot, &l1, &ch[1], &c[1], &wa[iw], &wa[ix2], &wa[ix3], isign); - } - na = 1 - na; - } - else if (ip == 2) { - if (na == 0) { - pass2_(&idot, &l1, &c[1], &ch[1], &wa[iw], isign); - } - else { - pass2_(&idot, &l1, &ch[1], &c[1], &wa[iw], isign); - } - na = 1 - na; - } - else if (ip == 3) { - ix2 = iw + idot; - if (na == 0) { - pass3_(&idot, &l1, &c[1], &ch[1], &wa[iw], &wa[ix2], isign); - } - else { - pass3_(&idot, &l1, &ch[1], &c[1], &wa[iw], &wa[ix2], isign); - } - na = 1 - na; - } - else if (ip == 5) { - ix2 = iw + idot; - ix3 = ix2 + idot; - ix4 = ix3 + idot; - if (na == 0) { - pass5_(&idot, &l1, &c[1], &ch[1], &wa[iw], &wa[ix2], &wa[ix3], - &wa[ix4], isign); - } - else { - pass5_(&idot, &l1, &ch[1], &c[1], &wa[iw], &wa[ix2], &wa[ix3], - &wa[ix4], isign); - } - na = 1 - na; - } - else { - if (na == 0) { - pass_(&nac, &idot, &ip, &l1, &idl1, &c[1], &c[1], &c[1], &ch[1], - &ch[1], &wa[iw], isign); - } - else { - pass_(&nac, &idot, &ip, &l1, &idl1, &ch[1], &ch[1], &ch[1], &c[1], - &c[1], &wa[iw], isign); - } - if (nac != 0) { - na = 1 - na; - } - } - l1 = l2; - iw += (ip - 1) * idot; - } - if (na != 0) { - n2 = *n + *n; - i_1 = n2; - for (i = 1; i <= i_1; ++i) { - c[i] = ch[i]; - } - } - return 0; -} - static int cffti1_(n, wa, ifac) int *n; double *wa; @@ -825,3 +672,157 @@ static int pass5_(ido, l1, cc, ch, wa1, wa2, wa3, wa4, isign) } return 0; } + + +/*static */ +int +cfft1_(int *n, double *c, double *ch, double *wa, + double /* int??*/ *ifac, int *isign) +{ + /* System generated locals */ + int i_1; + + /* Local variables */ + int idot; + int i, k1, l1, l2, n2, na, nf, ip, iw, ix2, ix3, ix4, nac, ido, idl1; + + /* Parameter adjustments */ + --c; + --ch; + --wa; + --ifac; + + /* Function Body */ + nf = ifac[2]; + na = 0; + l1 = 1; + iw = 1; + i_1 = nf; + for (k1 = 1; k1 <= i_1; ++k1) { + ip = ifac[k1 + 2]; + l2 = ip * l1; + ido = *n / l2; + idot = ido + ido; + idl1 = idot * l1; + if (ip == 4) { + ix2 = iw + idot; + ix3 = ix2 + idot; + if (na == 0) { + pass4_(&idot, &l1, &c[1], &ch[1], &wa[iw], &wa[ix2], &wa[ix3], isign); + } + else { + pass4_(&idot, &l1, &ch[1], &c[1], &wa[iw], &wa[ix2], &wa[ix3], isign); + } + na = 1 - na; + } + else if (ip == 2) { + if (na == 0) { + pass2_(&idot, &l1, &c[1], &ch[1], &wa[iw], isign); + } + else { + pass2_(&idot, &l1, &ch[1], &c[1], &wa[iw], isign); + } + na = 1 - na; + } + else if (ip == 3) { + ix2 = iw + idot; + if (na == 0) { + pass3_(&idot, &l1, &c[1], &ch[1], &wa[iw], &wa[ix2], isign); + } + else { + pass3_(&idot, &l1, &ch[1], &c[1], &wa[iw], &wa[ix2], isign); + } + na = 1 - na; + } + else if (ip == 5) { + ix2 = iw + idot; + ix3 = ix2 + idot; + ix4 = ix3 + idot; + if (na == 0) { + pass5_(&idot, &l1, &c[1], &ch[1], &wa[iw], &wa[ix2], &wa[ix3], + &wa[ix4], isign); + } + else { + pass5_(&idot, &l1, &ch[1], &c[1], &wa[iw], &wa[ix2], &wa[ix3], + &wa[ix4], isign); + } + na = 1 - na; + } + else { + if (na == 0) { + pass_(&nac, &idot, &ip, &l1, &idl1, &c[1], &c[1], &c[1], &ch[1], + &ch[1], &wa[iw], isign); + } + else { + pass_(&nac, &idot, &ip, &l1, &idl1, &ch[1], &ch[1], &ch[1], &c[1], + &c[1], &wa[iw], isign); + } + if (nac != 0) { + na = 1 - na; + } + } + l1 = l2; + iw += (ip - 1) * idot; + } + if (na != 0) { + n2 = *n + *n; + i_1 = n2; + for (i = 1; i <= i_1; ++i) { + c[i] = ch[i]; + } + } + return 0; +} + + +/* + Public Routine +*/ + +/* + * cfft computes the forward or backward complex discrete fourier transform. + * + * Input Parameters: + * + * n The length of the complex sequence c. The method is + * more efficient when n is the product of small primes. + * + * c A complex array of length n which contains the sequence + * + * wsave a real work array which must be dimensioned at least 4n+15 + * in the program that calls cfft. + * isign 1 for transform, -1 for inverse transform. + * A call of cfft with isign = 1 followed by a call of cfft with + * isign = -1 will multiply the sequence by n. + * + * Output Parameters: + * + * c For j=1,...,n + * + * c(j)=the sum from k=1,...,n of + * + * c(k)*exp(-isign*i*(j-1)*(k-1)*2*pi/n) + * + * where i=sqrt(-1) + */ + +int cfft(n, c, wsave, isign) + int n; + double *c, *wsave; + int isign; +{ + int iw1, iw2; + + /* Parameter adjustments */ + --c; + --wsave; + + /* Function Body */ + if (n != 1) { + iw1 = n + n + 1; + iw2 = iw1 + n + n; + cffti1_(&n, &wsave[iw1], &wsave[iw2]); + cfft1_(&n, &c[1], &wsave[1], &wsave[iw1], &wsave[iw2], &isign); + } + return 0; +} diff --git a/lib/cholesky.c b/lib/cholesky.c index 21d1e18..8d649fc 100644 --- a/lib/cholesky.c +++ b/lib/cholesky.c @@ -28,16 +28,14 @@ choldecomp(a, n) } */ -static double Max(a, b) - double a, b; +static double +Max(double a, double b) { return(a > b ? a : b); } -choldecomp(a, n, maxoffl, maxadd) - RMatrix a; - int n; - double maxoffl, *maxadd; +void +choldecomp(double **a, int n, double maxoffl, double *maxadd) { double minl, minljj, minl2; int i, j, k; diff --git a/lib/clinalg.c b/lib/clinalg.c index 77bfe70..cafe65f 100644 --- a/lib/clinalg.c +++ b/lib/clinalg.c @@ -11,8 +11,31 @@ typedef char *PTR; extern double rcondest(); -int min (x, y) int x, y; { return((x < y) ? x : y); } -int max (x, y) int x, y; { return((x > y) ? x : y); } +extern void choldecomp(); +extern int crludcmp(); +extern int crlubksb(); +extern int svdcmp(); +extern void qrdecomp(); +extern void make_rotation(); +extern void eigen(); +extern void cfft(); +extern int lowess(); +extern int fit_spline(); +extern int kernel_smooth(); + + + +int +min (int x, int y) +{ + return((x < y) ? x : y); +} + +int +max (int x, int y) +{ + return((x > y) ? x : y); +} /************************************************************************/ /** **/ @@ -37,31 +60,27 @@ double macheps() /** **/ /************************************************************************/ -chol_decomp_front(mat, n, dpars) - PTR mat, dpars; - int n; +void +chol_decomp_front(PTR mat, int n, PTR dpars) { double *dp = (double *) dpars; choldecomp((double **) mat, n, *dp, dp + 1); } -int lu_decomp_front(mat, n, iv, mode, dp) - PTR mat, iv, dp; - int n, mode; +int +lu_decomp_front(PTR mat, int n, PTR iv, int mode, PTR dp) { return(crludcmp((char **) mat, n, (int *) iv, mode, (double *) dp)); } -int lu_solve_front(a, n, indx, b, mode) - PTR a, indx, b; - int n, mode; +int +lu_solve_front(PTR a, int n, PTR indx, PTR b, int mode) { return(crlubksb((char **) a, n, (int *) indx, (char *) b, mode)); } -int lu_inverse_front(pmat, n, piv, pv, mode, pinv) - PTR pmat, piv, pv, pinv; - int n, mode; +int +lu_inverse_front(PTR pmat, int n, PTR piv, PTR pv, int mode, PTR pinv) { Matrix mat = (Matrix) pmat, inv = (Matrix) pinv; IVector iv = (IVector) piv; @@ -98,38 +117,32 @@ int lu_inverse_front(pmat, n, piv, pv, mode, pinv) return(singular); } -sv_decomp_front(mat, m, n, w, v) - PTR mat, w, v; - int m, n; +int +sv_decomp_front(PTR mat, int m, int n, PTR w, PTR v) { return(svdcmp((char **) mat, m, n, (char *) w, (char **) v)); } -qr_decomp_front(mat, m, n, v, jpvt, pivot) - PTR mat, v, jpvt; - int m, n, pivot; +void +qr_decomp_front(PTR mat, int m, int n, PTR v, PTR jpvt, int pivot) { qrdecomp((char **) mat, m, n, (char **) v, (char *) jpvt, pivot); } -double rcondest_front(mat, n) - PTR mat; - int n; +double +rcondest_front(PTR mat, int n) { return(rcondest((char **) mat, n)); } -make_rotation_front(n, rot, x, y, use_alpha, alpha) - int n, use_alpha; - PTR rot, x, y; - double alpha; +void +make_rotation_front(int n, PTR rot, PTR x, PTR y, int use_alpha, double alpha) { make_rotation(n, (char **) rot, (char *) x, (char *) y, use_alpha, alpha); } -int eigen_front(a, n, w, z, fv1) - PTR a, w, z, fv1; - int n; +int +eigen_front(PTR a, int n, PTR w, PTR z, PTR fv1) { int ierr; @@ -137,25 +150,22 @@ int eigen_front(a, n, w, z, fv1) return(ierr); } -fft_front(n, x, work, isign) - int n, isign; - PTR x, work; +void +fft_front(int n, PTR x, PTR work, int isign) { cfft(n, (char *) x, (char *) work, isign); } -int base_lowess_front(x, y, n, f, nsteps, delta, ys, rw, res) - PTR x, y, ys, rw, res; - int n, nsteps; - double f, delta; +int +base_lowess_front(PTR x, PTR y, int n, double f, + int nsteps, double delta, PTR ys, PTR rw, PTR res) { return(lowess((char *) x, (char *) y, n, f, nsteps, delta, (char *) ys, (char *) rw, (char *) res)); } -range_to_rseq(n, px, ns, pxs) - int n, ns; - PTR px, pxs; +void +range_to_rseq(int n, PTR px, int ns, PTR pxs) { int i; double xmin, xmax, *x, *xs; @@ -170,18 +180,16 @@ range_to_rseq(n, px, ns, pxs) xs[i] = xmin + (xmax - xmin) * ((double) i) / ((double) (ns - 1)); } -int spline_front(n, x, y, ns, xs, ys, work) - PTR x, y, xs, ys, work; - int n, ns; +int +spline_front(int n, PTR x, PTR y, int ns, + PTR xs, PTR ys, PTR work) { return(fit_spline(n, (char *) x, (char *) y, ns, (char *) xs, (char *) ys, (char *) work)); } -kernel_dens_front(x, n, width, xs, ys, ns, code) - PTR x, xs, ys; - int n, ns, code; - double width; +int +kernel_dens_front(PTR x, int n, double width, PTR xs, PTR ys, int ns, int code) { int ktype; @@ -194,10 +202,9 @@ kernel_dens_front(x, n, width, xs, ys, ns, code) (char *) xs, (char *) ys, ns, ktype)); } -int kernel_smooth_front(x, y, n, width, xs, ys, ns, code) - PTR x, y, xs, ys; - int n, ns, code; - double width; +int +kernel_smooth_front(PTR x, PTR y, int n, double width, + PTR xs, PTR ys, int ns, int code) { int ktype; diff --git a/lib/complex.c b/lib/complex.c index 10711b6..4056b2e 100644 --- a/lib/complex.c +++ b/lib/complex.c @@ -1,11 +1,16 @@ /* complex - Complex number functions */ /* Copyright (c) 1990, by Luke Tierney */ +/* patched up and semi-ansified, (c) 2006, AJ Rossini, blindglobe@gmail.com */ + #include "xmath.h" #include "complex.h" -static double phase(c) - Complex c; +extern void xlfail(char *); + +/*static */ +double +phase(Complex c) { double phi; @@ -13,26 +18,31 @@ static double phase(c) if (c.imag > 0.0) phi = PI / 2; else if (c.imag == 0.0) phi = 0.0; else phi = -PI / 2; - } - else { + } else { phi = atan(c.imag / c.real); if (c.real < 0.0) { - if (c.imag > 0.0) phi += PI; - else if (c.imag < 0.0) phi -= PI; - else phi = PI; + if (c.imag > 0.0) { + phi += PI; + } else { + if (c.imag < 0.0) { + phi -= PI; + } else { + phi = PI; + } + } } } return(phi); } -double modulus(c) - Complex c; +double +modulus(Complex c) { return(sqrt(c.real * c.real + c.imag * c.imag)); } -Complex cart2complex(real, imag) - double real, imag; +Complex +cart2complex(double real, double imag) { Complex val; val.real = real; @@ -40,8 +50,8 @@ Complex cart2complex(real, imag) return(val); } -static Complex polar2complex(mod, phi) - double mod, phi; +Complex +polar2complex(double mod, double phi) { Complex val; double cs, sn; @@ -49,42 +59,46 @@ static Complex polar2complex(mod, phi) if (phi == 0) { cs = 1.0; sn = 0.0; - } - else if (phi == PI / 2) { - cs = 0.0; - sn = 1.0; - } - else if (phi == PI) { - cs = -1.0; - sn = 0.0; - } - else if (phi == -PI / 2) { - cs = 0.0; - sn = -1.0; - } - else { - cs = cos(phi); - sn = sin(phi); + } else { + if (phi == PI / 2) { + cs = 0.0; + sn = 1.0; + } else { + if (phi == PI) { + cs = -1.0; + sn = 0.0; + } else { + if (phi == -PI / 2) { + cs = 0.0; + sn = -1.0; + } else { + cs = cos(phi); + sn = sin(phi); + } + } + } } val.real = mod * cs; val.imag = mod * sn; return(val); } -Complex cadd(c1, c2) - Complex c1, c2; +Complex +cadd(Complex c1, + Complex c2) { return(cart2complex(c1.real + c2.real, c1.imag + c2.imag)); } -Complex csub(c1, c2) - Complex c1, c2; +Complex +csub(Complex c1, + Complex c2) { return(cart2complex(c1.real - c2.real, c1.imag - c2.imag)); } -Complex cmul(c1, c2) - Complex c1, c2; +Complex +cmul(Complex c1, Complex c2) { double m1, m2, p1, p2; @@ -95,8 +109,16 @@ Complex cmul(c1, c2) return(polar2complex(m1 * m2, p1 + p2)); } -Complex cdiv(c1, c2) - Complex c1, c2; +static void +checkfzero(double x) +{ + if (x == 0.0) { + xlfail("division by zero"); + } +} + +Complex +cdiv(Complex c1, Complex c2) { double m1, m2, p1, p2; @@ -107,9 +129,3 @@ Complex cdiv(c1, c2) checkfzero(m2); return(polar2complex(m1 / m2, p1 - p2)); } - -static checkfzero(x) - double x; -{ - if (x == 0.0) xlfail("division by zero"); -} diff --git a/lib/complex.h b/lib/complex.h index 54ade1e..3658a8f 100644 --- a/lib/complex.h +++ b/lib/complex.h @@ -1,14 +1,18 @@ +/* I hate C programming... */ + typedef struct { double real, imag; } Complex; #ifndef PI #define PI 3.141592653589793 -#endif PI +#endif /* PI */ extern Complex makecomplex(); -extern Complex cart2complex(), polar2complex(); -extern Complex csqrt(), cexp(), clog(), cexpt(), csin(), ccos(), ctan(); -extern Complex casin(), cacos(), catan(); +extern Complex cart2complex(); +extern Complex polar2complex(double, double); +/* extern Complex csqrt(), cexp(), clog(), cexpt(), csin(), ccos(), ctan(); + extern Complex casin(), cacos(), catan(); */ extern Complex cadd(), csub(), cmul(), cdiv(); -extern double phase(), modulus(); +extern /* static */ double phase(Complex); +extern double modulus(Complex); diff --git a/lib/derivatives.c b/lib/derivatives.c index c4852c5..6a70383 100644 --- a/lib/derivatives.c +++ b/lib/derivatives.c @@ -7,17 +7,15 @@ # include "xmath.h" extern double macheps(); -typedef double **RMatrix, *RVector; - # define nil 0L # define FALSE 0 # define TRUE 1 -numergrad(n, x, grad, fsum, ffun, h, typx) - int n; - RVector x, grad, fsum, typx; - int (*ffun)(); - double h; +void +numergrad(int n, double *x, double *grad, double *fsum, + int ffun(double *, double *, double *, double **), + double h, double *typx) + /* int (*ffun)();*/ { int i; double old_xi, f1, f2, hi; @@ -35,12 +33,10 @@ numergrad(n, x, grad, fsum, ffun, h, typx) } } -numerhess(n, x, hess, f, fsum, ffun, h, typx) - int n; - RVector x, fsum, typx; - RMatrix hess; - int (*ffun)(); - double h, f; +void +numerhess(int n, double *x, double **hess, double f, double *fsum, + int *ffun(double *, double *, double *, double *), + double h, double *typx) { int i, j; double old_xi, old_xj, f1, f2, hi, hj; diff --git a/lib/exclglue.c b/lib/exclglue.c deleted file mode 100644 index ab795a9..0000000 --- a/lib/exclglue.c +++ /dev/null @@ -1,385 +0,0 @@ -#include "linalg.h" - -extern double rcondest_front(); - -extern double unirand(), gamma(); -extern double normalcdf(), normalquant(), normaldens(), normalrand(), - bnormcdf(); -extern double cauchycdf(), cauchyquant(), cauchydens(), cauchyrand(); -extern double gammacdf(), gammaquant(), gammadens(), gammarand(); -extern double chisqcdf(), chisqquant(), chisqdens(), chisqrand(); -extern double betacdf(), betaquant(), betadens(), betarand(); -extern double tcdf(), tquant(), tdens(), trand(); -extern double fcdf(), fquant(), fdens(), frand(); -extern double poissoncdf(), poissonpmf(); -extern int poissonquant(), poissonrand(); -extern double binomialcdf(), binomialpmf(); -extern int binomialquant(), binomialrand(); - -/***************************************************************************/ -/***************************************************************************/ -/**** ****/ -/**** Basic Utilities ****/ -/**** ****/ -/***************************************************************************/ -/***************************************************************************/ - -/***************************************************************************/ -/** **/ -/** Callback Value Storage **/ -/** **/ -/***************************************************************************/ - -static int excl_integer_value; -static double excl_double_value; - -excl_set_integer_value(x) - int x; -{ - excl_integer_value = x; -} - -excl_set_double_value(x) - double x; -{ - excl_double_value = x; -} - -/***************************************************************************/ -/** **/ -/** Storage Allocation Functions **/ -/** **/ -/***************************************************************************/ - -int la_base_allocate(n, m) - unsigned n, m; -{ - char *p = calloc(n, m); - if (p == nil) xlfail("allocation failed"); - return((int) p); -} - -int la_base_free_alloc(p) - int p; -{ - if (p) free((char *) p); - return(0); -} - -int la_mode_size(mode) - int mode; -{ - switch (mode) { - case IN: return(sizeof(int)); - case RE: return(sizeof(double)); - case CX: return(sizeof(Complex)); - } - return(0); -} - -/***************************************************************************/ -/** **/ -/** Callbacks for Internal Storage **/ -/** **/ -/***************************************************************************/ - -int la_allocate_index, la_free_alloc_index; - -excl_register_la_allocate(f) int f; { la_allocate_index = f; } -excl_register_la_free_alloc(f) int f; { la_free_alloc_index = f; } - -int la_allocate(n, m) - int n, m; -{ - lisp_call(la_allocate_index, n, m); - return(excl_integer_value); -} - -la_free_alloc(p) - int p; -{ - lisp_call(la_free_alloc_index, p); -} - -/***************************************************************************/ -/** **/ -/** Storage Access Functions **/ -/** **/ -/***************************************************************************/ - -int la_get_integer(p, i) - int p, i; -{ - return(*(((int *) p) + i)); -} - -double la_get_double(p, i) - int p, i; -{ - return(*(((double *) p) + i)); -} - -double la_get_complex_real(p, i) - int p, i; -{ - Complex *c = ((Complex *) p) + i; - return(c->real); -} - -double la_get_complex_imag(p, i) - int p, i; -{ - Complex *c = ((Complex *) p) + i; - return(c->imag); -} - -/***************************************************************************/ -/** **/ -/** Storage Mutation Functions **/ -/** **/ -/***************************************************************************/ - -int la_put_integer(p, i, x) - int p, i, x; -{ - *(((int *) p) + i) = x; - return(0); -} - -int la_put_double(p, i, x) - int p, i; - double x; -{ - *(((double *) p) + i) = x; - return(0); -} - -int la_put_complex(p, i, x, y) - int p, i; - double x, y; -{ - Complex *c = ((Complex *) p) + i; - c->real = x; - c->imag = y; - return(0); -} - -/***************************************************************************/ -/** **/ -/** XLISP internal error message emulation **/ -/** **/ -/***************************************************************************/ - -char buf[1000]; - -static int excl_set_buf_char_index; -excl_register_set_buf_char(f) int f; { excl_set_buf_char_index = f; } -set_buf_char(n, c) int n, c; { lisp_call(excl_set_buf_char_index, n, c); } - -static int excl_print_buffer_index; -excl_register_print_buffer(f) int f; { excl_print_buffer_index = f; } -print_buffer(n, m) int n, m; { lisp_call(excl_print_buffer_index, n, m); } - -static int bufpos = 0; - -static resetbuf() { bufpos = 0; } - -static prbuf(s) - char *s; -{ - int i, n; - - n = strlen(s); - for (i = 0; i +#include +#include "xmath.h" + #define PRINTSTR(s) printf(s) -#else -# include "xmath.h" -#define PRINTSTR(s) stdputstr(s) -#endif SBAYES -extern char *S_alloc(), *calloc(), *realloc(); +extern void *S_alloc(), *calloc(), *realloc(); extern double macheps(); char *minresultstring(); +extern void choldecomp(); + +extern void minsupplyvalues(double , double *, double **, + double *, double **); +extern void minimize(); +extern void minresults(double *, double *, double *, double *, double **, + double *, double **, int *, int *, double *); +extern void minsetup(int n, int k, + int *ffun(), int *gfun(), + double *x, double typf, double *typx, char *work); + /* int (*ffun)(), (*gfun)();*/ +extern void minsetoptions(double , double, double, int , int , int , int); + +extern void choldecomp(); + + + +/* next 2 from cbayes.c */ +extern void Recover(char *, char *); +extern void call_S(char *fun, long narg, char **args, char **mode, long *length,char **names, + long nvals, char **values); + +/* next 2 from derivatives.c */ +/* +extern void numergrad(int n, double *x, double *grad, double *fsum, + int *ffun(), double h, double *typx); +extern void numerhess(int n, double *x, double **hess, double f, + double *fsum, void ffun(), + double h, double *typx); +*/ +extern void numergrad(int , double *, double *, double *, + int ffun(double *, double *, double *, double **), + double , double *); +extern void numerhess(int , double *, double **, double , + double *, + int ffun(double *, double *, double *, double **), + double , double *); + +/* next from minimize */ +extern int minworkspacesize(int, int); + /************************************************************************/ /** **/ /** Definitions and Globals **/ /** **/ /************************************************************************/ +/* #define NULL 0L already defined in stddef */ + #define nil 0L -#define NULL 0L #define TRUE 1 #define FALSE 0 @@ -33,7 +72,7 @@ char *minresultstring(); #define GRADTOL_POWER 1.0 / 3.0 #define H_POWER 1.0 / 6.0 -typedef double **RMatrix, *RVector; +/*typedef double **RMatrix, *RVector;*/ typedef struct{ char *f, **sf, **g; @@ -41,8 +80,8 @@ typedef struct{ int change_sign, fderivs; int *gderivs; double typf, h, dflt; - RVector typx, fsum, cvals, ctarget; - RMatrix gfsum; + double *typx, *fsum, *cvals, *ctarget; + double **gfsum; } Fundata; static Fundata func, gfuncs, cfuncs; @@ -59,7 +98,8 @@ static Fundata func, gfuncs, cfuncs; /* called repeatedly (but not recursively) from within the same call */ /* from S. It attempts to avoid the danger of dangling callocs. */ -static makespace(pptr, size) +void static +makespace(pptr, size) char **pptr; int size; { @@ -76,16 +116,14 @@ static makespace(pptr, size) /************************************************************************/ /* - * All Hessianevaluations by numerical derivatives assume the gradient is + * All Hessian evaluations by numerical derivatives assume the gradient is * evaluated first at the same location. The results are cached away. */ /* install log posterior function */ -static install_func(f, sf, n, change_sign, typf, h, typx, dflt) - char *f, **sf; - int n; - double typf, h, dflt; - RVector typx; +static void +install_func(char *f, char **sf, int n, int change_sign, /*?? */ + double typf, double h, double *typx, double dflt) { int i; static int inited = FALSE; @@ -111,11 +149,8 @@ static install_func(f, sf, n, change_sign, typf, h, typx, dflt) } /* install tilt functions */ -static install_gfuncs(g, n, k, change_sign, h, typx) - char **g; - int n, k, change_sign; - double h; - RVector typx; +static void +install_gfuncs(char **g, int n, int k, int change_sign, double h, double *typx) { int i; static int inited = FALSE; @@ -140,14 +175,12 @@ static install_gfuncs(g, n, k, change_sign, h, typx) for (i = 0; i < n; i++) gfuncs.typx[i] = (typx != nil && typx[i] > 0.0) ? typx[i] : 1.0; for (i = 0; i < k; i++) gfuncs.gfsum[i] = gfsumdata + i * n; + return; } /* install constraint functions */ -static install_cfuncs(g, n, k, ctarget, h, typx) - char **g; - int n, k; - double h; - RVector typx, ctarget; +static void +install_cfuncs(char **g, int n, int k, double *ctarget, double h, double *typx) { int i; static int inited = FALSE; @@ -169,21 +202,20 @@ static install_cfuncs(g, n, k, ctarget, h, typx) for (i = 0; i < n; i++) cfuncs.typx[i] = (typx != nil && typx[i] > 0.0) ? typx[i] : 1.0; cfuncs.ctarget = ctarget; + return; } -/* callback to test if x is in the support of the posterior */ -static in_support(ff, n, x) - char **ff; - int n; - double *x; +static int +in_support(char **ff, int n, double *x) { char *args[1], *values[1]; int *result; char *mode[1]; long length[1]; - if (ff == nil || ff[0] == nil) return(TRUE); - else { + if (ff == nil || ff[0] == nil) { + return(TRUE); + } else { mode[0] = "double"; length[0] =n; args[0] = (char *) x; @@ -194,10 +226,8 @@ static in_support(ff, n, x) } /* callback for logposterior evaluation */ -static evalfunc(x, pval, grad, hess) - RVector x, grad; - double *pval; - RMatrix hess; +static int +evalfunc(double *x, double *pval, double *grad, double **hess) { char *args[1], *values[3]; double *result, val; @@ -246,21 +276,21 @@ static evalfunc(x, pval, grad, hess) } } return(TRUE); - } - else { - if (pval != nil) *pval = func.dflt; + } else { + if (pval != nil) { + *pval = func.dflt; + } return(FALSE); } + return(TRUE); } /* callback for tilt function evaluation */ static int which_gfunc; -static evalgfunc(x, pval, grad, hess) - RVector x, grad; - double *pval; - RMatrix hess; +static int +evalgfunc(double *x, double *pval, double *grad, double **hess) { char *args[1], *values[3]; double *result, val; @@ -308,15 +338,14 @@ static evalgfunc(x, pval, grad, hess) gfuncs.h, gfuncs.typx); } } + return(TRUE); } /* callback for constraint function evaluation */ static int which_cfunc; -static evalcfunc(x, pval, grad, hess) - RVector x, grad; - double *pval; - RMatrix hess; +static int +evalcfunc(double *x, double *pval, double *grad, double **hess) { char *args[1], *values[3]; double *result, val; @@ -365,16 +394,16 @@ static evalcfunc(x, pval, grad, hess) cfuncs.h, cfuncs.typx); } } + return(TRUE); } /* S front end for logposterior evaluation */ -evalfront(ff, n, x, val, grad, phess, h, typx) - char **ff; - int *n; - double *x, *val, *grad, *phess, *typx, *h; +void +evalfront(char **ff, int *n, double *x, double *val, double *grad, + double *phess, double *h, double *typx) { int i; - static RMatrix hess = nil; + static double **hess = nil; install_func(ff[0], nil, *n, FALSE, 1.0, *h, typx, 0.0); if (phess == nil) hess = nil; @@ -385,12 +414,10 @@ evalfront(ff, n, x, val, grad, phess, h, typx) evalfunc(x, val, grad, hess); } -#ifdef SBAYES /* S front end for tilt function evaluation */ -gevalfront(gg, n, m, x, h, typx, val, grad) - char **gg; - int *n, *m; - double *x, *h, *typx, *val, *grad; +void +gevalfront(char **gg, int *n, int *m, double *x, double *h, + double *typx, double *val, double *grad) { int i; @@ -400,6 +427,7 @@ gevalfront(gg, n, m, x, h, typx, val, grad) evalgfunc(x, val, grad, nil); if (grad != nil) grad += *n; } + return; } /************************************************************************/ @@ -408,17 +436,33 @@ gevalfront(gg, n, m, x, h, typx, val, grad) /** **/ /************************************************************************/ -static check_derivs(x, drvtol) - RVector x; - double drvtol; +/* +void +derivscalefront(char **ff, int *n, double *x, double *h, double *typx, double *tol, int *info) +{ + int i; + + if (*tol <= 0.0) *tol = pow(macheps(), GRADTOL_POWER); + + install_func(ff[0], nil, *n, TRUE, 1.0, *h, typx, 0.0); + *info = check_derivs(x, *tol); + + *h = func.h; + for (i = 0; i < *n; i++) typx[i] = func.typx[i]; + return; +} + + +static int +check_derivs(double *x, double drvtol) { - static RVector grad = nil, work = nil; - static RMatrix hess = nil; + static double *grad = nil, work = nil; + static double **hess = nil; int i, error; - grad = (RVector) S_alloc(func.n, sizeof(double)); - hess = (RMatrix) S_alloc(func.n, sizeof(double *)); - work = (RVector) S_alloc(func.n + func.n * func.n, sizeof(double)); + grad = (double *) S_alloc(func.n, sizeof(double)); + hess = (double **) S_alloc(func.n, sizeof(double *)); + work = (double *) S_alloc(func.n + func.n * func.n, sizeof(double)); for (i = 0; i < func.n; i++) { hess[i] = work; @@ -429,22 +473,7 @@ static check_derivs(x, drvtol) &func.h, func.typx, drvtol, work); return(error); } - -derivscalefront(ff, n, x, h, typx, tol, info) - char **ff; - int *n, *info; - double *x, *h, *typx, *tol; -{ - int i; - - if (*tol <= 0.0) *tol = pow(macheps(), GRADTOL_POWER); - - install_func(ff[0], nil, *n, TRUE, 1.0, *h, typx, 0.0); - *info = check_derivs(x, *tol); - - *h = func.h; - for (i = 0; i < *n; i++) typx[i] = func.typx[i]; -} +*/ /************************************************************************/ /** **/ @@ -453,9 +482,8 @@ derivscalefront(ff, n, x, h, typx, tol, info) /************************************************************************/ /* joint density of normal-cauchy mixture */ -static double dncmix(x, n, p) - double *x, p; - int n; +static double +dncmix(double *x, int n, double p) { int i; double dens; @@ -471,10 +499,10 @@ static double dncmix(x, n, p) * S front end for computing sample from transformed normal-cauchy * mixture and importance sampling weights */ -samplefront(ff, sf, rf, p, n, x, ch, N, y, w) - char **ff, **sf, **rf; - int *n, *N; - double *p, *x, *ch, *y, *w; +void +samplefront(char **ff, char **sf, char **rf, + double *p, int *n, + double *x, double *ch, int *N, double *y, double *w) { double val; int i, j, k; @@ -505,8 +533,10 @@ samplefront(ff, sf, rf, p, n, x, ch, N, y, w) if (evalfunc(y, &val, nil, nil)) w[i] = exp(val - mval) * c / dens; else w[i] = 0.0; } + return; } -#endif SBAYES + + /************************************************************************/ /** **/ /** Maximization Routines **/ @@ -524,71 +554,15 @@ typedef struct { struct { double tilt; - RVector gval; - RMatrix ggrad, ghess; + double *gval; + double **ggrad, **ghess; int exptilt; - RVector tscale; + double *tscale; } tiltinfo; -static set_tilt_info(n, m, tilt, exptilt, tscale) - int n, m; - double tilt, *tscale; - int exptilt; -{ - static double *hessdata = nil, *graddata = nil; - int i; - static int inited = FALSE; - - if (! inited) { - tiltinfo.gval = nil; - tiltinfo.ggrad = nil; - tiltinfo.ghess = nil; - inited = TRUE; - } - makespace(&tiltinfo.gval, n * sizeof(double)); - makespace(&tiltinfo.ggrad, m * sizeof(double *)); - makespace(&tiltinfo.ghess, n * sizeof(double *)); - makespace(&graddata, n * m * sizeof(double)); - makespace(&hessdata, n * n * sizeof(double)); - - tiltinfo.tilt = tilt; - tiltinfo.exptilt = exptilt; - for (i = 0; i < m; i++) tiltinfo.ggrad[i] = graddata + i * n; - for (i = 0; i < n; i++) tiltinfo.ghess[i] = hessdata + i * n; - tiltinfo.tscale = tscale; -} - -static minfunc(x, pval, grad, hess) - RVector x, grad; - double *pval; - RMatrix hess; -{ - int k = gfuncs.k; - - if (evalfunc(x, pval, grad, hess) && (k > 0)) - add_tilt(x, pval, grad, hess, tiltinfo.tilt, tiltinfo.exptilt); -} - -constfunc(x, vals, jac, hess) - RVector x, vals; - RMatrix jac, hess; -{ - int i, k = cfuncs.k; - double *pvali, *jaci; - - for (i = 0; i < k; i++) { - pvali = (vals != nil) ? vals + i : nil; - jaci = (jac != nil) ? jac[i] : nil; - which_cfunc = i; - evalcfunc(x, pvali, jaci, nil); - } -} - -static add_tilt(x, pval, grad, hess, tilt, exptilt) - RVector x, grad; - double *pval, tilt; - RMatrix hess; - int exptilt; +static void +add_tilt(double *x, double *pval, double *grad, double **hess, + double tilt, int exptilt) { int i, j, k, n = func.n, m = gfuncs.k; double *gval, *ggrad, **ghess, etilt; @@ -632,19 +606,68 @@ static add_tilt(x, pval, grad, hess, tilt, exptilt) } } -maxfront(ff, gf, cf, x, typx, fvals, gvals, cvals, ctarget, ipars, dpars, - tscale, msg) - char **ff, **gf, **cf; - double *x, *typx, *fvals, *gvals, *cvals, *ctarget, *tscale; - MaxIPars *ipars; - MaxDPars *dpars; - char **msg; +static void +set_tilt_info(int n, int m, + double tilt, int exptilt, double *tscale) +{ + static double *hessdata = nil, *graddata = nil; + int i; + static int inited = FALSE; + + if (! inited) { + tiltinfo.gval = nil; + tiltinfo.ggrad = nil; + tiltinfo.ghess = nil; + inited = TRUE; + } + makespace(&tiltinfo.gval, n * sizeof(double)); + makespace(&tiltinfo.ggrad, m * sizeof(double *)); + makespace(&tiltinfo.ghess, n * sizeof(double *)); + makespace(&graddata, n * m * sizeof(double)); + makespace(&hessdata, n * n * sizeof(double)); + + tiltinfo.tilt = tilt; + tiltinfo.exptilt = exptilt; + for (i = 0; i < m; i++) tiltinfo.ggrad[i] = graddata + i * n; + for (i = 0; i < n; i++) tiltinfo.ghess[i] = hessdata + i * n; + tiltinfo.tscale = tscale; +} + + +static void +minfunc(double *x, double *pval, double *grad, double **hess) +{ + int k = gfuncs.k; + + if (evalfunc(x, pval, grad, hess) && (k > 0)) + add_tilt(x, pval, grad, hess, tiltinfo.tilt, tiltinfo.exptilt); +} + +void +constfunc(double *x, double *vals, double **jac, double **hess) +{ + int i, k = cfuncs.k; + double *pvali, *jaci; + + for (i = 0; i < k; i++) { + pvali = (vals != nil) ? vals + i : nil; + jaci = (jac != nil) ? jac[i] : nil; + which_cfunc = i; + evalcfunc(x, pvali, jaci, nil); + } +} + +void +maxfront(char **ff, char **gf, char **cf, + double *x, double *typx, double *fvals, double *gvals, double *cvals, double *ctarget, + MaxIPars *ipars, MaxDPars *dpars, + double *tscale, char **msg) { static char *work = nil; - static RMatrix H = nil, cJ = nil; + static double **H = nil, **cJ = nil; double *pf, *grad, *c; int i, n, m, k; - int (*cfun)(); + void (*cfun)(); if (ipars->verbose > 0) PRINTSTR("maximizing...\n"); @@ -671,14 +694,14 @@ maxfront(ff, gf, cf, x, typx, fvals, gvals, cvals, ctarget, ipars, dpars, cvals += k; makespace(&cJ, k * sizeof(double *)); for (i = 0; i < k; i++) cJ[i] = cvals + i * n; - cfun = constfunc; + cfun = &constfunc; /* AJR: pointer to constfunc? */ } install_func(ff[0], nil, n, TRUE, dpars->typf, dpars->h, typx, dpars->dflt); install_gfuncs(gf, n, m, TRUE, dpars->h, typx); install_cfuncs(cf, n, k, ctarget, dpars->h, typx); - minsetup(n, k, minfunc, cfun, x, dpars->typf, typx, work); + minsetup(n, k, &minfunc, &cfun, x, dpars->typf, typx, work); /* AJR: FIXME arg 3,4 */ minsetoptions(dpars->gradtol, dpars->steptol, dpars->maxstep, ipars->itnlimit, ipars->verbose, ipars->backtrack, TRUE); @@ -706,17 +729,14 @@ maxfront(ff, gf, cf, x, typx, fvals, gvals, cvals, ctarget, ipars, dpars, /** **/ /************************************************************************/ -loglapdet(fvals, cvals, ipars, dpars, val, detonly) - double *fvals, *cvals; - MaxIPars *ipars; - MaxDPars *dpars; - double *val; - int *detonly; +void +loglapdet(double *fvals, double *cvals, MaxIPars *ipars, MaxDPars *dpars, + double *val, int *detonly) { int i, j, l, n = ipars->n, k = ipars->k; double f = -fvals[0], *hessdata = fvals + n + 1, *cgraddata = cvals + k; double ldL, ldcv, maxadd; - static RMatrix hess = nil, cgrad = nil; + static double **hess = nil, **cgrad = nil; if (k >= n) Recover("too many constraints", NULL); @@ -784,12 +804,12 @@ moms1front(gf, n, m, x, hessdata, mean, stdev, sigmadata, h, typx) double *x, *hessdata, *mean, *stdev, *sigmadata, *h, *typx; { int i, j, k; - RMatrix hess, sigma, delg; + double *hess, *sigma, *delg; double *delgdata, maxadd; - hess = (RMatrix) S_alloc(*n, sizeof(double *)); - sigma = (RMatrix) S_alloc(*m, sizeof(double *)); - delg = (RMatrix) S_alloc(*m, sizeof(double *)); + hess = (double **) S_alloc(*n, sizeof(double *)); + sigma = (double **) S_alloc(*m, sizeof(double *)); + delg = (double **) S_alloc(*m, sizeof(double *)); delgdata = (double *) S_alloc(*m * *n, sizeof(double)); for (i = 0; i < *n; i++) hess[i] = hessdata + i * *n; @@ -1097,9 +1117,15 @@ static cpmarpars(x, f, g, ip, dp, x1, f1, g1, ip1, dp1) } #endif /* SBAYES */ -#ifdef TODO -get hessian from gradiant for analytical gradiants -avoid repeated derivative calls in mimimize. -2d margins -use pos. definiteness info in margins -#endif TODO +/* + TODO + + get hessian from gradiant for analytical gradiants + + avoid repeated derivative calls in mimimize. + + 2d margins + + use pos. definiteness info in margins + +*/ diff --git a/lib/gammabase.c b/lib/gammabase.c index 5faf343..66f833e 100644 --- a/lib/gammabase.c +++ b/lib/gammabase.c @@ -3,18 +3,19 @@ #define TRUE 1 #define FALSE 0 +extern void normbase(double *,double *); extern double gamma(), ppnd(); static double gammp(), gser(), gcf(), gnorm(), ppchi2(); -gammabase(x, a, p) - double *x, *a, *p; + +void +gammabase(double *x, double *a, double *p) { *p = gammp(*a, *x); } -double ppgamma(p, a, ifault) - double p, a; - int *ifault; +double +ppgamma(double p, double a, int *ifault) { double x, v, g; diff --git a/lib/kernel.c b/lib/kernel.c index cdee735..e00550f 100644 --- a/lib/kernel.c +++ b/lib/kernel.c @@ -2,15 +2,14 @@ #ifndef ROOT2PI #define ROOT2PI 2.50662827463100050241 -#endif ROOT2PI +#endif /* ROOT2PI*/ #ifndef nil #define nil 0L -#endif nil +#endif /*nil */ -static double kernel(x, y, w, type) - double x, y, w; - int type; +static double +kernel(double x, double y, double w, int type) { double z, k; @@ -50,9 +49,10 @@ static double kernel(x, y, w, type) return(k); } -kernel_smooth(x, y, n, width, wts, wds, xs, ys, ns, ktype) - double *x, *y, width, *wts, *wds, *xs, *ys; - int n, ns, ktype; +int +kernel_smooth(double *x, double *y, int n, + double width, double *wts, double *wds, + double *xs, double *ys, int ns, int ktype) { int i, j; double wsum, ysum, lwidth, lwt, xmin, xmax; diff --git a/lib/linalg.h b/lib/linalg.h index 6cf505c..f8e65c9 100644 --- a/lib/linalg.h +++ b/lib/linalg.h @@ -2,7 +2,7 @@ # include "xlisp.h" # include "complex.h" -extern char *calloc(); +extern void *calloc(); /*conflict? was (char *) */ extern double macheps(); #define nil 0L @@ -24,3 +24,4 @@ extern IMatrix imatrix(); extern RMatrix rmatrix(); extern CMatrix cmatrix(); +extern void free_vector(double *); diff --git a/lib/linalgdata.c b/lib/linalgdata.c index f81f25b..85a9e68 100644 --- a/lib/linalgdata.c +++ b/lib/linalgdata.c @@ -9,7 +9,10 @@ typedef int PTR; typedef char *PTR; #endif -extern PTR la_allocate(); +extern PTR la_allocate(int, int); +extern void la_free_alloc(PTR); +extern void xlfail(char *); + /************************************************************************/ /** **/ @@ -17,30 +20,30 @@ extern PTR la_allocate(); /** **/ /************************************************************************/ -static char *allocate(n, m) - unsigned n, m; +static char +*allocate(unsigned n, unsigned m) { char *p = (char *) la_allocate(n, m); if (p == nil) xlfail("allocation failed"); return(p); } -static free_alloc(p) - char *p; +static void +free_alloc(char *p) { if (p != nil) la_free_alloc((PTR) p); } -IVector ivector(n) - unsigned n; +IVector +ivector(unsigned n) { return((IVector) allocate(n, sizeof(int))); } -RVector rvector(n) - unsigned n; +double +*rvector(unsigned n) { - return((RVector) allocate(n, sizeof(double))); + return((double *) allocate(n, sizeof(double))); } CVector cvector(n) @@ -49,10 +52,13 @@ CVector cvector(n) return((CVector) allocate(n, sizeof(Complex))); } -free_vector(v) Vector v; { free_alloc(v); } +void +free_vector(double *v) +{ + free_alloc((char *)v); +} -IMatrix imatrix(n, m) - unsigned n, m; +IMatrix imatrix(unsigned n, unsigned m) { int i; IMatrix mat = (IMatrix) allocate(n, sizeof(IVector)); @@ -60,8 +66,8 @@ IMatrix imatrix(n, m) return(mat); } -RMatrix rmatrix(n, m) - unsigned n, m; +RMatrix +rmatrix(unsigned n, unsigned m) { int i; RMatrix mat = (RMatrix) allocate(n, sizeof(RVector)); @@ -69,8 +75,8 @@ RMatrix rmatrix(n, m) return(mat); } -CMatrix cmatrix(n, m) - unsigned n, m; +CMatrix +cmatrix(unsigned n, unsigned m) { int i; CMatrix mat = (CMatrix) allocate(n, sizeof(CVector)); @@ -78,12 +84,11 @@ CMatrix cmatrix(n, m) return(mat); } -free_matrix(mat, n) - Matrix mat; - int n; +void +free_matrix(double **mat, int n) /* Matrix?? Not RMatrix?? */ { int i; - if (mat != nil) for (i = 0; i < n; i++) free_alloc(mat[i]); - free_alloc(mat); + if (mat != nil) for (i = 0; i < n; i++) free_alloc((char *)mat[i]); + free_alloc((char *)mat); } diff --git a/lib/lowess.c b/lib/lowess.c index fb3875f..9d229fe 100644 --- a/lib/lowess.c +++ b/lib/lowess.c @@ -7,13 +7,88 @@ #define FALSE 0 #define TRUE 1 +extern int max(int, int); +extern int min(int, int); + static double pow2(x) double x; { return(x * x); } static double pow3(x) double x; { return(x * x * x); } static double fmax(x,y) double x, y; { return (x > y ? x : y); } -lowess(x, y, n, f, nsteps, delta, ys, rw, res) - double *x, *y, f, delta, *ys, *rw, *res; - int n, nsteps; +int +static compar(double *a, double *b) +{ + if (*a < *b) return(-1); + else if (*a > *b) return(1); + else return(0); +} + +static void +lowest(double *x, double *y, int n, double xs, double *ys, int nleft, int nright, + double *w, int userw, double *rw, int *ok) +{ + double range, h, h1, h9, a, b, c, r; + int j, nrt; + + range = x[n - 1] - x[0]; + h = fmax(xs - x[nleft], x[nright] - xs); + h9 = .999 * h; + h1 = .001 * h; + + /* compute weights (pick up all ties on right) */ + a = 0.0; /* sum of weights */ + for(j = nleft; j < n; j++) { + w[j]=0.0; + r = fabs(x[j] - xs); + if (r <= h9) { /* small enough for non-zero weight */ + if (r > h1) w[j] = pow3(1.0-pow3(r/h)); + else w[j] = 1.0; + if (userw) w[j] = rw[j] * w[j]; + a += w[j]; + } + else if (x[j] > xs) break; /* get out at first zero wt on right */ + } + nrt = j - 1; /* rightmost pt (may be greater than nright because of ties) */ + if (a <= 0.0) *ok = FALSE; + else { /* weighted least squares */ + *ok = TRUE; + + /* make sum of w[j] == 1 */ + for (j = nleft; j <= nrt; j++) w[j] = w[j] / a; + + if (h > 0.0) { /* use linear fit */ + + /* find weighted center of x values */ + for (j = nleft, a = 0.0; j <= nrt; j++) a += w[j] * x[j]; + + b = xs - a; + for (j = nleft, c = 0.0; j <= nrt; j++) + c += w[j] * (x[j] - a) * (x[j] - a); + + if(sqrt(c) > .001 * range) { + /* points are spread out enough to compute slope */ + b = b/c; + for (j = nleft; j <= nrt; j++) + w[j] = w[j] * (1.0 + b*(x[j] - a)); + } + } + for (j = nleft, *ys = 0.0; j <= nrt; j++) *ys += w[j] * y[j]; + } +} + +static void +sort(double *x, int n) +{ + extern void qsort(); + + qsort(x, n, sizeof(double), compar); +} + + + +int +lowess(double *x, double *y, int n, + double f, int nsteps, + double delta, double *ys, double *rw, double *res) { int iter, ns, ok, nleft, nright, i, j, last, m1, m2; double d1, d2, denom, alpha, cut, cmad, c9, c1, r; @@ -35,7 +110,11 @@ lowess(x, y, n, f, nsteps, delta, ys, rw, res) nleft++; nright++; } - lowest(x, y, n, x[i], &ys[i], nleft, nright, res, (iter > 1), rw, &ok); + lowest(x, y, + n, x[i], + &ys[i], + nleft, nright, + res, (iter > 1), rw, &ok); /* fitted value at x[i] */ if (! ok) ys[i] = y[i]; /* all weights zero - copy over value (all rw==0) */ @@ -78,71 +157,4 @@ lowess(x, y, n, f, nsteps, delta, ys, rw, res) } -static lowest(x, y, n, xs, ys, nleft, nright, w, userw, rw, ok) - double *x, *y, *w, *rw, xs, *ys; - int n, nleft, nright, userw, *ok; -{ - double range, h, h1, h9, a, b, c, r; - int j, nrt; - - range = x[n - 1] - x[0]; - h = fmax(xs - x[nleft], x[nright] - xs); - h9 = .999 * h; - h1 = .001 * h; - - /* compute weights (pick up all ties on right) */ - a = 0.0; /* sum of weights */ - for(j = nleft; j < n; j++) { - w[j]=0.0; - r = fabs(x[j] - xs); - if (r <= h9) { /* small enough for non-zero weight */ - if (r > h1) w[j] = pow3(1.0-pow3(r/h)); - else w[j] = 1.0; - if (userw) w[j] = rw[j] * w[j]; - a += w[j]; - } - else if (x[j] > xs) break; /* get out at first zero wt on right */ - } - nrt = j - 1; /* rightmost pt (may be greater than nright because of ties) */ - if (a <= 0.0) *ok = FALSE; - else { /* weighted least squares */ - *ok = TRUE; - - /* make sum of w[j] == 1 */ - for (j = nleft; j <= nrt; j++) w[j] = w[j] / a; - - if (h > 0.0) { /* use linear fit */ - - /* find weighted center of x values */ - for (j = nleft, a = 0.0; j <= nrt; j++) a += w[j] * x[j]; - - b = xs - a; - for (j = nleft, c = 0.0; j <= nrt; j++) - c += w[j] * (x[j] - a) * (x[j] - a); - - if(sqrt(c) > .001 * range) { - /* points are spread out enough to compute slope */ - b = b/c; - for (j = nleft; j <= nrt; j++) - w[j] = w[j] * (1.0 + b*(x[j] - a)); - } - } - for (j = nleft, *ys = 0.0; j <= nrt; j++) *ys += w[j] * y[j]; - } -} -static compar(a, b) - double *a, *b; -{ - if (*a < *b) return(-1); - else if (*a > *b) return(1); - else return(0); -} - -static sort(x, n) - double *x; - int n; -{ - qsort(x, n, sizeof(double), compar); -} - diff --git a/lib/ludecomp.c b/lib/ludecomp.c index 02ad677..19fa1d0 100644 --- a/lib/ludecomp.c +++ b/lib/ludecomp.c @@ -7,11 +7,8 @@ #include "linalg.h" -crludcmp(mat, n, indx, mode, d) - Matrix mat; - IVector indx; - int n, mode; - double *d; +int +crludcmp(Matrix mat, int n, IVector indx, int mode, double *d) { int i, imax, j, k, singular = FALSE; double big, temp; @@ -107,11 +104,8 @@ crludcmp(mat, n, indx, mode, d) return(singular); } -crlubksb(a, n, indx, b, mode) - Matrix a; - IVector indx; - Vector b; - int n, mode; +int +crlubksb(double **a, int n, int *indx, double *b, int mode) { int i, ii, ip, j, singular = FALSE; CMatrix ca = (CMatrix) a; diff --git a/lib/makerotation.c b/lib/makerotation.c index 1e665be..006d53b 100644 --- a/lib/makerotation.c +++ b/lib/makerotation.c @@ -18,11 +18,8 @@ static double inner_product(n, x, y) #define NORM(n, x) (sqrt(inner_product(n, x, x))) -make_rotation(n, rot, x, y, use_alpha, alpha) - int n, use_alpha; - RMatrix rot; - RVector x, y; - double alpha; +void +make_rotation(int n, double **rot, double *x, double *y, int use_alpha, double alpha) { double nx, ny, xy, c, s; int i, j; diff --git a/lib/minimize.c b/lib/minimize.c index 454aeed..e6440c8 100644 --- a/lib/minimize.c +++ b/lib/minimize.c @@ -10,18 +10,53 @@ * Equations." */ -#ifdef SBAYES -# include -static char buf[200]; +#include +#include "xmath.h" +extern char buf[200]; #define PRINTSTR(s) printf(s) -#else -# include "xmath.h" -extern char buf[]; -#define PRINTSTR(s) stdputstr(s) -#endif SBAYES + +typedef struct { + int n, k; + void (*ffun)(), (*gfun)(); + double f, typf, new_f; + double crit, new_crit; + double *x, *new_x, *sx, *delf, *new_delf, *qnstep, *F; + double **hessf, **H, **L, **new_delg; + double gradtol, steptol, maxstep; + int itncount, itnlimit, maxtaken, consecmax, retcode, termcode; + int use_line_search, verbose, values_supplied, change_sign; + double diagadd; +} Iteration; + extern double macheps(); + +extern double choldecomp(); +static void eval_gradient(Iteration *); +static void eval_next_gradient(Iteration *); + + +/* +This is a little more advanced section on functions, but is very useful. Take this for example: + int applyeqn(int F(int), int max, int min) { + int itmp; + + itmp = F(int) + min; + itmp = itmp - max; + return itmp; + } + + What does this function do if we call it with applyeqn(square(x), + y, z);? What happens is that the int F(int) is a reference to the + function that is passed in as a parameter. Thus inside applyeqn where + there is a call to F, it actually is a call to square! This is very + useful if we have one set function, but wish to vary the input + according to a particular function. So if we had a different function + called cube we could change how we call applyeqn by calling the + function by applyeqn(cube(x), y, z);. */ + + /************************************************************************/ /** **/ /** Definitions and Globals **/ @@ -44,18 +79,6 @@ extern double macheps(); typedef double **RMatrix, *RVector; -typedef struct { - int n, k; - int (*ffun)(), (*gfun)(); - double f, typf, new_f; - double crit, new_crit; - RVector x, new_x, sx, delf, new_delf, qnstep, F; - RMatrix hessf, H, L, new_delg; - double gradtol, steptol, maxstep; - int itncount, itnlimit, maxtaken, consecmax, retcode, termcode; - int use_line_search, verbose, values_supplied, change_sign; - double diagadd; -} Iteration; static char *termcodes[] = {"not yet terminated", "gradient size is less than gradient tolerance", @@ -88,28 +111,10 @@ static double Min(a, b) /** **/ /************************************************************************/ -/* solve (L L^T) s = -g for s */ -static cholsolve(n, g, L, s) - int n; - RVector g, s; - RMatrix L; -{ - int i; - - /* solve Ly = g */ - lsolve(n, g, L, s); - - /* solve L^Ts = y */ - ltsolve(n, s, L, s); - - for (i = 0; i < n; i++) s[i] = -s[i]; -} /* solve Ly = b for y */ -static lsolve(n, b, L, y) - int n; - RVector b, y; - RMatrix L; +static void +lsolve(int n, double *b, double **L, double *y) { int i, j; @@ -121,10 +126,8 @@ static lsolve(n, b, L, y) } /* solve (L^T)x = y for x */ -static ltsolve(n, y, L, x) - int n; - RVector y, x; - RMatrix L; +static void +ltsolve(int n, double *y, double **L, double *x) { int i, j; @@ -135,11 +138,25 @@ static ltsolve(n, y, L, x) } } -static modelhess(n, sx, H, L, diagadd) - int n; - RVector sx; - RMatrix H, L; - double *diagadd; +/* solve (L L^T) s = -g for s */ +static void +cholsolve(int n, + double *g, double **L, + double *s) +{ + int i; + + /* solve Ly = g */ + lsolve(n, g, L, s); + + /* solve L^Ts = y */ + ltsolve(n, s, L, s); + + for (i = 0; i < n; i++) s[i] = -s[i]; +} + +static void +modelhess(int n, double *sx, double **H, double **L, double *diagadd) { int i, j; double sqrteps, maxdiag, mindiag, maxoff, maxoffl, maxposdiag, mu, @@ -232,13 +249,12 @@ static modelhess(n, sx, H, L, diagadd) /** **/ /************************************************************************/ -static double gradsize(iter, new) - Iteration *iter; - int new; +static double +gradsize(Iteration *iter, int new) { int n, i; double size, term, crit, typf; - RVector x, delf, sx; + double *x, *delf, *sx; n = iter->n + iter->k; crit = iter->crit; @@ -254,12 +270,12 @@ static double gradsize(iter, new) return(size); } -static double incrsize(iter) - Iteration *iter; +static double +incrsize(Iteration *iter) { int n, i; double size, term; - RVector x, new_x, sx; + double *x, *new_x, *sx; new_x = iter->new_x; n = iter->n + iter->k; @@ -273,8 +289,8 @@ static double incrsize(iter) return(size); } -static stoptest0(iter) - Iteration *iter; +static int /* AJR:TESTME guess! */ +stoptest0(Iteration *iter) { iter->consecmax = 0; @@ -285,8 +301,8 @@ static stoptest0(iter) return(iter->termcode); } -static stoptest(iter) - Iteration *iter; +static int +stoptest(Iteration *iter) { int termcode, retcode, itncount, itnlimit, maxtaken, consecmax; double gradtol, steptol; @@ -322,8 +338,8 @@ static stoptest(iter) /** **/ /************************************************************************/ -static eval_funval(iter) - Iteration *iter; +static void +eval_funval(Iteration *iter) { int i; @@ -336,8 +352,8 @@ static eval_funval(iter) } } -static eval_next_funval(iter) - Iteration *iter; +static void +eval_next_funval(Iteration *iter) { int i; @@ -350,8 +366,8 @@ static eval_next_funval(iter) } } -static eval_gradient(iter) - Iteration *iter; +static void +eval_gradient(Iteration *iter) { int i, j, n, k; @@ -373,8 +389,8 @@ static eval_gradient(iter) } } -static eval_next_gradient(iter) - Iteration *iter; +static void +eval_next_gradient(Iteration *iter) { int i, j, n, k; @@ -391,8 +407,7 @@ static eval_next_gradient(iter) } } -static eval_hessian(iter) - Iteration *iter; +static void eval_hessian(Iteration *iter) { int i, j, n, k; @@ -409,14 +424,14 @@ static eval_hessian(iter) /** **/ /************************************************************************/ -static linesearch(iter) - Iteration *iter; +static void +linesearch(Iteration *iter) { int i, n; double newtlen, maxstep, initslope, rellength, lambda, minlambda, lambdatemp, lambdaprev, a, b, disc, critprev, f1, f2, a11, a12, a21, a22, del; - RVector qnstep, delf, x, new_x, sx; + double *qnstep, *delf, *x, *new_x, *sx; n = iter->n + iter->k; if (! iter->use_line_search) { @@ -510,8 +525,8 @@ static linesearch(iter) /** **/ /************************************************************************/ -static print_header(iter) - Iteration *iter; +static void +print_header(Iteration *iter) { if (iter->verbose > 0) { sprintf(buf, "Iteration %d.\n", iter->itncount); @@ -519,8 +534,8 @@ static print_header(iter) } } -static print_status(iter) - Iteration *iter; +static void +print_status(Iteration *iter) { int i, j; @@ -566,8 +581,8 @@ static print_status(iter) /** **/ /************************************************************************/ -static findqnstep(iter) - Iteration *iter; +static void +findqnstep(Iteration *iter) { int i, j, N, l; @@ -589,8 +604,7 @@ static findqnstep(iter) } } -static iterupdate(iter) - Iteration *iter; +static void iterupdate(Iteration *iter) { int i, j, n, k; @@ -610,8 +624,8 @@ static iterupdate(iter) } } -static mindriver(iter) - Iteration *iter; +static void +mindriver(Iteration *iter) { iter->consecmax = 0; iter->itncount = 0; @@ -646,8 +660,8 @@ static mindriver(iter) static Iteration myiter; -minworkspacesize(n, k) - int n, k; +int +minworkspacesize(int n, int k) { int size; @@ -664,19 +678,19 @@ minworkspacesize(n, k) return(size); } -char *minresultstring(code) - int code; +char * +minresultstring(int code) { if (code <= 0) return("bad input data"); else if (code <= 5) return(termcodes[code]); else return("unknown return code"); } -minsetup(n, k, ffun, gfun, x, typf, typx, work) - int n, k, (*ffun)(), (*gfun)(); - RVector x, typx; - double typf; - char *work; +void +minsetup(int n, int k, + void (*ffun)(), void (*gfun)(), + double *x, double typf, double *typx, char *work) + /* int (*ffun)(), (*gfun)();*/ { Iteration *iter = &myiter; int i, j; @@ -687,16 +701,16 @@ minsetup(n, k, ffun, gfun, x, typf, typx, work) iter->n = n; iter->k = k; - iter->ffun = ffun; + iter->ffun = ffun; /* FIXME:AJR */ iter->gfun = gfun; - iter->x = (RVector) work; work += sizeof(double) * (n + k); - iter->new_x = (RVector) work; work += sizeof(double) * (n + k); - iter->sx = (RVector) work; work += sizeof(double) * (n + k); - iter->delf = (RVector) work; work += sizeof(double) * (n + k); - iter->new_delf = (RVector) work; work += sizeof(double) * (n + k); - iter->qnstep = (RVector) work; work += sizeof(double) * (n + k); - iter->F = (RVector) work; work += sizeof(double) * (n + k); + iter->x = (double *) work; work += sizeof(double) * (n + k); + iter->new_x = (double *) work; work += sizeof(double) * (n + k); + iter->sx = (double *) work; work += sizeof(double) * (n + k); + iter->delf = (double *) work; work += sizeof(double) * (n + k); + iter->new_delf = (double *) work; work += sizeof(double) * (n + k); + iter->qnstep = (double *) work; work += sizeof(double) * (n + k); + iter->F = (double *) work; work += sizeof(double) * (n + k); for (i = 0; i < n; i++) { iter->x[i] = x[i]; iter->sx[i] = (typx != nil && typx[i] > 0.0) ? 1.0 / typx[i] : 1.0; @@ -706,28 +720,28 @@ minsetup(n, k, ffun, gfun, x, typf, typx, work) iter->sx[n + i] = 1.0; } - iter->hessf = (RMatrix) work; work += sizeof(double *) * (n + k); + iter->hessf = (double **) work; work += sizeof(double *) * (n + k); for (i = 0; i < n + k; i++) { - iter->hessf[i] = (RVector) work; + iter->hessf[i] = (double *) work; work += sizeof(double) * (n + k); } for (i = 0; i < n + k; i++) for (j = 0; j < n + k; j++) iter->hessf[i][j] = 0.0; - iter->L = (RMatrix) work; work += sizeof(double *) * (n + k); + iter->L = (double **) work; work += sizeof(double *) * (n + k); for (i = 0; i < n + k; i++) { - iter->L[i] = (RVector) work; + iter->L[i] = (double *) work; work += sizeof(double) * (n + k); } - iter->H = (RMatrix) work; work += sizeof(double *) * (n + k); + iter->H = (double **) work; work += sizeof(double *) * (n + k); for (i = 0; i < n + k; i++) { - iter->H[i] = (RVector) work; + iter->H[i] = (double *) work; work += sizeof(double) * (n + k); } - iter->new_delg = (k > 0) ? (RMatrix) work : nil; + iter->new_delg = (k > 0) ? (double **) work : nil; work += sizeof(double *) * k; for (i = 0; i < k; i++) { - iter->new_delg[i] = (RVector) work; + iter->new_delg[i] = (double *) work; work += sizeof(double) * n; } @@ -747,9 +761,9 @@ minsetup(n, k, ffun, gfun, x, typf, typx, work) iter->values_supplied = FALSE; } -minsetoptions(gradtol, steptol, maxstep, itnlimit, verbose, use_search, change_sign) - double gradtol, steptol, maxstep; - int itnlimit, verbose, use_search, change_sign; +void +minsetoptions(double gradtol, double steptol, double maxstep, + int itnlimit, int verbose, int use_search, int change_sign) { Iteration *iter = &myiter; @@ -762,10 +776,9 @@ minsetoptions(gradtol, steptol, maxstep, itnlimit, verbose, use_search, change_s iter->change_sign = change_sign; } -minsupplyvalues(f, delf, hessf, g, delg) - double f; - RVector delf, g; - RMatrix hessf, delg; +void +minsupplyvalues(double f, double *delf, double **hessf, + double *g, double **delg) { Iteration *iter = &myiter; int i, j, n, k; @@ -795,13 +808,15 @@ minsupplyvalues(f, delf, hessf, g, delg) iter->values_supplied = TRUE; } -minimize() { mindriver(&myiter); } +void +minimize() +{ + mindriver(&myiter); +} -minresults(x, pf, pcrit, delf, hessf, g, delg, pcount, ptermcode, diagadd) - RVector x, delf, g; - double *pf, *pcrit, *diagadd; - RMatrix hessf, delg; - int *pcount, *ptermcode; +void +minresults(double *x, double *pf, double *pcrit, double *delf, double **hessf, + double *g, double **delg, int *pcount, int *ptermcode, double *diagadd) { Iteration *iter = &myiter; int i, j, n, k; @@ -827,10 +842,8 @@ minresults(x, pf, pcrit, delf, hessf, g, delg, pcount, ptermcode, diagadd) if (diagadd != nil) *diagadd = iter->diagadd; } -#ifdef SBAYES -double pdlogdet(n, H) - int n; - RMatrix H; + +double pdlogdet(int n, double **H) { int i; double logdet, maxadd; @@ -839,10 +852,17 @@ double pdlogdet(n, H) for (i = 0, logdet = 0.0; i < n; i++) logdet += 2.0 * log(H[i][i]); return(logdet); } -#endif /* SBAYES */ -#ifdef TODO -return amount added to make pos definite -scaling for constraints -alternate global strategies -callback function for verbose mode -#endif TODO + + +/* + TODO + + return amount added to make pos definite + + scaling for constraints + + alternate global strategies + + callback function for verbose mode +*/ + diff --git a/lib/nor.c b/lib/nor.c index 67fca0c..3e789cb 100644 --- a/lib/nor.c +++ b/lib/nor.c @@ -41,8 +41,8 @@ #define SQRT2 1.414213562373095049 #define SQRTPI 1.772453850905516027 -normbase(x, phi) -double *x, *phi; +void +normbase(double *x, double *phi) { int sn; double R1, R2, y, y2, y3, y4, y5, y6, y7, erf, erfc, z, z2, z3, z4; @@ -51,8 +51,9 @@ double *x, *phi; if (y < 0) { y = -y; sn = -1; + } else { + sn = 1; } - else sn = 1; y2 = y * y; y4 = y2 * y2; y6 = y4 * y2; @@ -61,30 +62,39 @@ double *x, *phi; R1 = P10 + P11 * y2 + P12 * y4 + P13 * y6; R2 = Q10 + Q11 * y2 + Q12 * y4 + Q13 * y6; erf = y * R1 / R2; - if(sn == 1) *phi = 0.5 + 0.5*erf; - else *phi = 0.5 - 0.5*erf; - } - else if(y < 4.0) { - y3 = y2 * y; - y5 = y4 * y; - y7 = y6 * y; - R1 = P20 + P21 * y + P22 * y2 + P23 * y3 + P24 * y4 + P25 * y5 - + P26 * y6 + P27 * y7; - R2 = Q20 + Q21 * y + Q22 * y2 + Q23 * y3 + Q24 * y4 + Q25 * y5 - + Q26 * y6 + Q27 * y7; - erfc = exp(-y2) * R1 / R2; - if(sn == 1) *phi = 1.0 - 0.5*erfc; - else *phi = 0.5*erfc; - } - else { - z = y4; - z2 = z * z; - z3 = z2 * z; - z4 = z2 * z2; - R1 = P30 + P31 * z + P32 * z2 + P33 * z3 + P34 * z4; - R2 = Q30 + Q31 * z + Q32 * z2 + Q33 * z3 + Q34 * z4; - erfc = (exp(-y2)/y) * (1.0 / SQRTPI + R1 / (R2 * y2)); - if(sn == 1) *phi = 1.0 - 0.5*erfc; - else *phi = 0.5*erfc; + if(sn == 1) { + *phi = 0.5 + 0.5*erf; + } else { + *phi = 0.5 - 0.5*erf; + } + } else { + if(y < 4.0) { + y3 = y2 * y; + y5 = y4 * y; + y7 = y6 * y; + R1 = P20 + P21 * y + P22 * y2 + P23 * y3 + P24 * y4 + P25 * y5 + + P26 * y6 + P27 * y7; + R2 = Q20 + Q21 * y + Q22 * y2 + Q23 * y3 + Q24 * y4 + Q25 * y5 + + Q26 * y6 + Q27 * y7; + erfc = exp(-y2) * R1 / R2; + if (sn == 1) { + *phi = 1.0 - (0.5 * erfc); + } else { + *phi = 0.5 * erfc; + } + } else { + z = y4; + z2 = z * z; + z3 = z2 * z; + z4 = z2 * z2; + R1 = P30 + P31 * z + P32 * z2 + P33 * z3 + P34 * z4; + R2 = Q30 + Q31 * z + Q32 * z2 + Q33 * z3 + Q34 * z4; + erfc = (exp(-y2)/y) * (1.0 / SQRTPI + R1 / (R2 * y2)); + if(sn == 1) { + *phi = 1.0 - ( 0.5 * erfc); + } else { + *phi = (0.5 * erfc); + } + } } } diff --git a/lib/qrdecomp.c b/lib/qrdecomp.c index b7c9708..8e04db8 100644 --- a/lib/qrdecomp.c +++ b/lib/qrdecomp.c @@ -4,9 +4,8 @@ #define SIGN(a, b) ((b) >= 0.0 ? fabs(a) : -fabs(a)) -static double NORM2(i, j, n, x) - int i, j, n; - double **x; +static double +NORM2(int i, int j, int n, double **x) { int k; double maxx, sum, temp; @@ -25,9 +24,8 @@ static double NORM2(i, j, n, x) } } -static double DOT(i, j, k, n, x) - int i, j, k; - double **x; +static double +DOT(int i, int j, int k, int n, double **x) { int l; double sum; @@ -36,25 +34,27 @@ static double DOT(i, j, k, n, x) return(sum); } -static AXPY(i, j, k, n, a, x) - int i, j, k, n; - double a, **x; +static void +AXPY(int i, int j, int k, int n, + double a, double **x) { int l; - for (l = i; l < n; l++) x[l][k] = a * x[l][j] + x[l][k]; + for (l = i; l < n; l++) { + x[l][k] = a * x[l][j] + x[l][k]; + } } -static SCALE(i, j, n, a, x) - int i, j, n; - double a, **x; +static void +SCALE(int i, int j, int n, double a, double **x) { int k; - for (k = i; k < n; k++) x[k][j] *= a; + for (k = i; k < n; k++) { + x[k][j] *= a; + } } -static SWAP(i, j, n, a) - int i, j, n; - double **a; +static void +SWAP(int i, int j, int n, double **a) { int k; double temp; @@ -65,16 +65,17 @@ static SWAP(i, j, n, a) } } -qrdecomp(x,n,p,v,jpvt,pivot) - int n, p, pivot; - int *jpvt; - double **x, **v; +void +qrdecomp(double **x, int n, int p, + double **v, int *jpvt, int pivot) { int i,j,k,jp,l,lp1,lup,maxj; double maxnrm,tt,*qraux,*work; double nrmxl,t; - if (n < 0) return; + if (n < 0) { + return; + } work = v[0]; qraux = rvector(p); @@ -181,6 +182,5 @@ qrdecomp(x,n,p,v,jpvt,pivot) } } } - free_vector(qraux); } diff --git a/lib/rcondest.c b/lib/rcondest.c index 8361e3b..ccb0e31 100644 --- a/lib/rcondest.c +++ b/lib/rcondest.c @@ -6,15 +6,42 @@ #include "linalg.h" -static double Max(a, b) - double a, b; +static +double Max(double a, double b) { return(a > b ? a : b); } -double rcondest(a, n) - RMatrix a; - int n; + +void +backsolve(double **a, double *b, int n, int mode) +{ + int i, j; + RMatrix ra = (RMatrix) a; + RVector rb = (RVector) b; + CMatrix ca = (CMatrix) a; + CVector cb = (CVector) b; + + switch (mode) { + case RE: + for (i = n - 1; i >= 0; i--) { + if (ra[i][i] != 0.0) rb[i] = rb[i] / ra[i][i]; + for (j = i + 1; j < n; j++) rb[i] -= ra[i][j] * rb[j]; + } + break; + case CX: + for (i = n - 1; i >= 0; i--) { + if (modulus(ca[i][i]) != 0.0) cb[i] = cdiv(cb[i], ca[i][i]); + for (j = i + 1; j < n; j++) + cb[i] = csub(cb[i], cmul(ca[i][j], cb[j])); + } + break; + } +} + +/** Exported function **/ + +double rcondest(double **a, int n) { RVector p, pm, x; double est, xp, xm, temp, tempm, xnorm; @@ -71,30 +98,3 @@ double rcondest(a, n) return(est); } -backsolve(a, b, n, mode) - Matrix a; - Vector b; - int n; -{ - int i, j; - RMatrix ra = (RMatrix) a; - RVector rb = (RVector) b; - CMatrix ca = (CMatrix) a; - CVector cb = (CVector) b; - - switch (mode) { - case RE: - for (i = n - 1; i >= 0; i--) { - if (ra[i][i] != 0.0) rb[i] = rb[i] / ra[i][i]; - for (j = i + 1; j < n; j++) rb[i] -= ra[i][j] * rb[j]; - } - break; - case CX: - for (i = n - 1; i >= 0; i--) { - if (modulus(ca[i][i]) != 0.0) cb[i] = cdiv(cb[i], ca[i][i]); - for (j = i + 1; j < n; j++) - cb[i] = csub(cb[i], cmul(ca[i][j], cb[j])); - } - break; - } -} diff --git a/lib/splines.c b/lib/splines.c index 2e63933..b9c95b4 100644 --- a/lib/splines.c +++ b/lib/splines.c @@ -3,9 +3,9 @@ /* natural cubic spline interpolation based on Numerical Recipes in C */ /* calculate second derivatives; assumes strictly increasing x values */ -static find_spline_derivs(x, y, n, y2, u) - double *x, *y, *y2, *u; - int n; +static void +find_spline_derivs(double *x, double *y, int n, + double *y2, double *u) { int i, k; double p, sig; @@ -36,9 +36,9 @@ static find_spline_derivs(x, y, n, y2, u) } /* interpolate or extrapolate value at x using results of find_spline_derivs */ -static spline_interp(xa, ya, y2a, n, x, y) - double *xa, *ya, *y2a, x, *y; - int n; +static void +spline_interp(double *xa, double *ya, double *y2a, + int n, double x, double *y) { int klo, khi, k; double h, b, a; @@ -82,9 +82,9 @@ static spline_interp(xa, ya, y2a, n, x, y) } } -fit_spline(n, x, y, ns, xs, ys, work) - int n, ns; - double *x, *y, *xs, *ys, *work; +int +fit_spline(int n, double *x, double *y, + int ns, double *xs, double *ys, double *work) { int i; double *y2, *u; diff --git a/lib/studentbase.c b/lib/studentbase.c index 8fe7e41..907f6c6 100644 --- a/lib/studentbase.c +++ b/lib/studentbase.c @@ -4,12 +4,15 @@ #define HALF_PI 1.5707963268 #define EPSILON .000001 +extern void betabase(double *, double *, double *,int *, int *,double *); +extern void normbase(double *,double *); + extern double ppnd(), ppbeta(); /* CACM Algorithm 395, by G. W. Hill */ -studentbase(x, df, cdf) - double *x, *df, *cdf; +void +studentbase(double *x, double *df, double *cdf) { double t, y, b, a, z, j, n; @@ -28,10 +31,11 @@ studentbase(x, df, cdf) dx = db / (db + da * t); betabase(&dx, &db, &da, &ib, &ia, &dp); *cdf = (*x >= 0) ? 1.0 - .5 * dp : .5 * dp; - } - else { + } else { /* asymptotic series for large or non-integer df */ - if(y > EPSILON) y = log(b); + if(y > EPSILON) { + y = log(b); + } a = n - 0.5; b = 48.0 * a * a; y = a * y; @@ -42,14 +46,12 @@ studentbase(x, df, cdf) normbase(&y, cdf); if (*x > 0.0) *cdf = 1.0 - *cdf; } - } - else { + } else { /* nested summation of cosine series */ if (n < 20 && t < 4.0) { a = y = sqrt(y); if(n == 1.0) a = 0.0; - } - else { + } else { a = sqrt(b); y = a * n; for(j = 2; fabs(a - z) > EPSILON; j += 2.0) { @@ -64,8 +66,11 @@ studentbase(x, df, cdf) for(n = n - 2.0; n > 1.0; n -= 2.0) a = ((n - 1.0) / (b * n)) * a + y; a = (fabs(n) < EPSILON) ? a/sqrt(b) : TWOVRPI * (atan(y) + a / b); *cdf = z - a; - if(*x > 0.0) *cdf = 1.0 - 0.5 * *cdf; - else *cdf = 0.5 * *cdf; + if(*x > 0.0) { + *cdf = 1.0 - 0.5 * *cdf; + } else { + *cdf = 0.5 * *cdf; + } } } diff --git a/lib/svdecomp.c b/lib/svdecomp.c index b83f7ef..ed8aad2 100644 --- a/lib/svdecomp.c +++ b/lib/svdecomp.c @@ -20,10 +20,9 @@ static double PYTHAG(a, b) #define SWAPD(a, b) (temp = (a), (a) = (b), (b) = temp) -static sort_sv(m, n, k, a, w, v) - int m, n, k; - RMatrix a, v; - RVector w; +static void +sort_sv(int m, int n, int k, + double **a, double *w, double **v) { int i, j; double temp; @@ -39,10 +38,8 @@ static double maxarg1, maxarg2; #define Max(a, b) (maxarg1 = (a), maxarg2 = (b), (maxarg1) > (maxarg2) ? (maxarg1) : (maxarg2)) #define SIGN(a, b) ((b) >= 0.0 ? fabs(a) : -fabs(a)) -svdcmp(a, m, n, w, v) - RMatrix a, v; - RVector w; - int m, n; +int +svdcmp(double **a, int m, int n, double *w, double **v) { int flag, i, its, j, jj, k, l, nm; double c, f, h, s, x, y, z; @@ -173,8 +170,8 @@ svdcmp(a, m, n, w, v) w[i] = h; if (h == 0.0) { char s[100]; - sprintf(s, "h = %f, f = %f, g = %f\n", f, g); - stdputstr(s); + sprintf(s, "h = %f, f = %f, g = %f\n", h, f, g); + fprintf(stdout,"%s",s); } h = 1.0 / h; c = g * h; -- 2.11.4.GIT