mgcv.h 11.5 KB
Newer Older
1
/* main method routines */
2
#include <Rinternals.h>
3
#include <Rconfig.h>
4 5 6 7 8
/* Most compilers with openMP support supply 
   a pre-defined compiler macro _OPENMP. Following 
   facilitates selective turning off (by testing value 
   or defining multiple versions OPENMP_ON1, OPENMP_ON2...)  */

9
#if defined _OPENMP
10
#define OPENMP_ON 1 
11 12
#endif
/* ... note also that there is no actual *need* to protect #pragmas with 
13
  #ifdef OPENMP_ON, since C ignores undefined pragmas, but failing 
14 15 16 17 18 19
  to do so may produce alot of compilation warnings if openMP is not supported. 
  In contrast functions from omp.h must be protected, and there is 
  non-avoidable use of these in the mgcv code. */

//#define OMP_REPORT // define to have all routines using omp report on start and end.

20 21 22 23 24 25 26 27 28
/* sed -i 's/old-text/new-text/g' *.c
   is quite useful!!
*/
// For safe memory handling from R...
#define CALLOC R_chk_calloc
#define FREE R_chk_free
// Can reset to check for memory errors...
//#define CALLOC calloc
//#define FREE free
29 30
void magic(double *y,double *X,double *sp0,double *def_sp,double *S,double *H,double *L,
	   double *lsp0,double *gamma,double *scale, int *control,int *cS,double *rank_tol,
31
	   double *tol,double *b,double *rV,double *norm_const,int *n_score,int *nt);
32

33 34 35
void gdi1(double *X,double *E,double *Es,double *rS,double *U1,
	  double *sp,double *z,double *w,double *wf,double *alpha,double *mu,double *eta, double *y,
	 double *p_weights,double *g1,double *g2,double *g3,double *g4,double *V0,
36
	  double *V1,double *V2,double *V3,double *beta,double *b1,double *w1,double *D1,double *D2,
37
         double *P0, double *P1,double *P2,double *trA,
38
         double *trA1,double *trA2,double *rV,double *rank_tol,double *conv_tol, int *rank_est,
39
	 int *n,int *q, int *M,int *Mp,int *Enrow,int *rSncol,int *deriv,
40
	  int *REML,int *fisher,int *fixed_penalty,int *nthreads,double *dVkk);     
41

42
void gdi2(double *X,double *E,double *Es,double *rS,double *U1,
43
	  double *sp,double *theta,double *z,double *w,double *wz,double *wf,
44 45 46
          double *Dth,double *Det,double *Det2,double *Dth2,double *Det_th,
          double *Det2_th,double *Det3,double *Det_th2,
          double *Det4, double *Det3_th, double *Det2_th2,
47
          double *beta,double *b1,double *w1,double *D1,double *D2,double *P,double *P1,double *P2,
48
          double *ldet, double *ldet1,double *ldet2,double *rV,
49 50
          double *rank_tol,int *rank_est,
	  int *n,int *q, int *M,int *n_theta, int *Mp,int *Enrow,int *rSncol,int *deriv,
51
	  int *fixed_penalty,int *nt,int *type,double *dVkk);
52

53 54
void pls_fit1(double *y,double *X,double *w,double *wy,double *E,double *Es,int *n,int *q,int *rE,double *eta,
	      double *penalty,double *rank_tol,int *nt,int *use_wy);
55 56 57 58 59 60 61 62

void get_detS2(double *sp,double *sqrtS, int *rSncol, int *q,int *M, int * deriv, 
               double *det, double *det1, double *det2, double *d_tol,
               double *r_tol,int *fixed_penalty); /* stable determinant of sum evaluation */

void get_stableS(double *S,double *Qf,double *sp,double *sqrtS, int *rSncol, int *q,int *M, int * deriv, 
               double *det, double *det1, double *det2, double *d_tol,
		 double *r_tol,int *fixed_penalty);
63

64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79
/* cox model routines */

void coxpred(double *X,double *t,double *beta,double *Vb,double *a,double *h,double *q,
             double *tr,int *n,int *p, int *nt,double *s,double *se);
void coxpp(double *eta,double *X,int *r, int *d,double *h,double *q,double *km,
	   int *n,int *p, int *nt);
void coxlpl(double *eta,double *X,int *r, int *d,double *tr, 
            int *n,int *p, int *nt,double *lp,double *g,double *H,
            double *d1beta,double *d1H,double *d2beta,
            double *d2H,int *n_sp,int *deriv);

/* MVN smooth additive */
void mvn_ll(double *y,double *X,double *XX,double *beta,int *n,int *lpi,
            int *m,double *ll,double *lb,double *lbb,double *dbeta,
            double *dH,int *deriv,int *nsp,int *nt);

80
/* discretized covariate methods */
81
void XWXd(double *XWX,double *X,double *w,int *k,int *ks, int *m,int *p, int *n, int *nx, 
82 83
          int *ts, int *dt, int *nt,double *v,int *qc,int *nthreads,int *ar_stop,
          int *ar_row,double *ar_weights);
84
void XWyd(double *XWy,double *y,double *X,double *w,int *k, int *ks, int *m,int *p, int *n, 
85 86
	  int *nx, int *ts, int *dt, int *nt,double *v,int *qc,
          int *ar_stop,int *ar_row,double *ar_weights);
87
void Xbd(double *f,double *beta,double *X,int *k, int *ks, int *m,int *p, int *n, 
88
	 int *nx, int *ts, int *dt, int *nt,double *v,int *qc,int *bc);
89
void diagXVXt(double *diag,double *V,double *X,int *k,int *ks,int *m,int *p, int *n, 
90
	      int *nx, int *ts, int *dt, int *nt,double *v,int *qc,int *pv,int *nthreads);
91

92
/* various service routines */
93

94 95 96
void tweedious(double *w,double *w1,double *w2, double *w1p,double *w2p,double *w2pp, 
	       double *y,double *eps,int *n,
               double *th,double *rho,double *a, double *b);
97 98 99
void tweedious2(double *w,double *w1,double *w2, double *w1p,double *w2p,double *w2pp, 
	       double *y,double *eps,int *n,
               double *th,double *rho,double *a, double *b);
100
void psum(double *y, double *x,int *index,int *n);
101
void rwMatrix(int *stop,int *row,double *w,double *X,int *n,int *p,int *trans,double *work);
102
void in_out(double *bx, double *by, double *break_code, double *x,double *y,int *in, int *nb, int *n);
103
void Rlanczos(double *A,double *U,double *D,int *n, int *m, int *lm,double *tol,int *nt);
104
void RuniqueCombs(double *X,int *ind,int *r, int *c);
105
void  RPCLS(double *Xd,double *pd,double *yd, double *wd,double *Aind,double *bd,double *Afd,double *Sd,int *off,int *dim,double *theta, int *m,int *nar);
106
void RMonoCon(double *Ad,double *bd,double *xd,int *control,double *lower,double *upper,int *n);
107 108
/*void MinimumSeparation(double *gx,double *gy,int *gn,double *dx,double *dy, int *dn,double *dist);*/
void MinimumSeparation(double *x,int *n, int *d,double *t,int *m,double *dist);
109
void rksos(double *x,int *n,double *eps);
110
void pivoter(double *x,int *r,int *c,int *pivot, int *col, int *reverse);
111

112 113 114
/* Routines for linear algebra with direct access to linpack and lapack */ 
void band_chol(double *B,int *n,int *k,int *info);
void tri_chol(double *ld,double *sd,int *n,int *info);
115
void mgcv_omp(int *a);
116 117 118
void mgcv_chol(double *a,int *pivot,int *n,int *rank);
void mgcv_svd(double *x,double *u, double *d,int *r,int *c);
void mgcv_qrqy(double *b,double *a,double *tau,int *r,int *c,int *k,int *left,int *tp);
119
void mgcv_qrqy0(double *b,double *a,double *tau,int *r,int *c,int *k,int *left,int *tp);
120 121
void mgcv_backsolve(double *R,int *r,int *c,double *B,double *C, int *bc, int *right);
void mgcv_forwardsolve(double *R,int *r,int *c,double *B,double *C, int *bc, int *right);
122
void mgcv_qr(double *x, int *r, int *c,int *pivot,double *tau);
123
void mgcv_qr2(double *x, int *r, int *c,int *pivot,double *tau);
124
void update_qr(double *Q,double *R,int *n, int *q,double *lam, int *k);
125
extern void mgcv_mmult(double *A,double *B,double *C,int *bt,int *ct,int *r,int *c,int *n);
126
void mgcv_pmmult(double *A,double *B,double *C,int *bt,int *ct,int *r,int *c,int *n,int *nt);
127
SEXP mgcv_pmmult2(SEXP b, SEXP c,SEXP bt,SEXP ct, SEXP nthreads);
128
void mgcv_mmult0(double *A,double *B,double *C,int *bt,int *ct,int *r,int *c,int *n);
129
void mgcv_svd_full(double *x,double *vt,double *d,int *r,int *c);
130
void mgcv_symeig(double *A,double *ev,int *n,int *use_dsyevd, int *get_vectors,int *descending);
131
void mroot(double *A,int *rank,int *n);
132
void R_cond(double *R,int *r,int *c,double *work,double *Rcondition);
133 134
void mgcv_td_qy(double *S,double *tau,int *m,int *n, double *B,int *left,int *transpose);
void mgcv_tri_diag(double *S,int *n,double *tau);
135
void mgcv_trisymeig(double *d,double *g,double *v,int *n,int getvec,int descending); 
136 137 138 139
void getXtWX(double *XtWX, double *X,double *w,int *r,int *c,double *work);
void getXtX(double *XtX,double *X,int *r,int *c);
void getXtMX(double *XtMX,double *X,double *M,int *r,int *c,double *work);
void getXXt(double *XXt,double *X,int *r,int *c);
140
void read_mat(double *M,int *r,int*c, char *path);
141 142 143 144
void row_block_reorder(double *x,int *r,int *c,int *nb,int *reverse);
void mgcv_pqr(double *x,int *r, int *c,int *pivot, double *tau, int *nt);
void getRpqr(double *R,double *x,int *r, int *c,int *rr,int *nt);
void mgcv_pqrqy(double *b,double *a,double *tau,int *r,int *c,int *cb,int *tp,int *nt);
145
SEXP mgcv_Rpiqr(SEXP X, SEXP BETA,SEXP PIV,SEXP NT,SEXP NB);
146
void mgcv_tmm(SEXP x,SEXP t,SEXP D,SEXP M, SEXP N);
147 148 149
void mgcv_Rpbsi(SEXP A, SEXP NT);
void mgcv_RPPt(SEXP a,SEXP r, SEXP NT);
SEXP mgcv_Rpchol(SEXP Amat,SEXP PIV,SEXP NT,SEXP NB);
150 151 152
void dchol(double *dA, double *R, double *dR,int *p);
void vcorr(double *dR,double *Vr,double *Vb,int *p,int *M);
SEXP mgcv_Rpforwardsolve(SEXP R, SEXP B,SEXP NT);
153
SEXP mgcv_Rpbacksolve(SEXP R, SEXP B,SEXP NT);
154
SEXP mgcv_Rpcross(SEXP A, SEXP NT,SEXP NB);
155 156
SEXP mgcv_madi(SEXP a, SEXP b,SEXP ind,SEXP diag);

157

158 159
/* basis constructor/prediction routines*/

160
void crspl(double *x,int *n,double *xk, int *nk,double *X,double *S, double *F,int *Fsupplied);
161 162 163
void predict_tprs(double *x, int *d,int *n,int *m,int *k,int *M,double *Xu,int *nXu,
                  double *UZ,double *by,int *by_exists,double *X);
void construct_tprs(double *x,int *d,int *n,double *knt,int *nk,int *m,int *k,double *X,double *S,
164
                    double *UZ,double *Xu,int *nXu,double *C);
165
void gen_tps_poly_powers(int *pi,int *M,int *m, int *d);
166 167 168 169 170
void boundary(int *G, double *d, double *dto, double *x0, double *y0, double *dx, double *dy,
              int *nx, int *ny, double *x, double *y,double *break_code, int *n, int *nb);
void gridder(double *z,double *x,double *y,int *n,double *g, int *G,int *nx, int *ny,double *x0, 
             double *y0,double *dx,double *dy,double NA_code);
void pde_coeffs(int *G,double *x,int *ii,int *jj,int *n,int *nx,int *ny,double *dx,double *dy);
171 172

/* sparse smooth related routines */
173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191
typedef struct { /* defines structure for kd-tree box */
  double *lo,*hi;    /* box defining co-ordinates */
  int parent,child1,child2, /* indices of parent and 2 offspring */
      p0,p1;         /* indices of first and last point in box */
} box_type; 



typedef struct {
  box_type *box;
  int *ind, /* index of points in coordinate matrix which tree relates to */
      *rind, /* where is ith row of X in ind? */
      n_box, /* number of boxes */
      d, /* dimension */
    n; /* number of points that tree relates to */
  double huge; /* number indicating an open boundary */
} kdtree_type;

void k_newn_work(double *Xm,kdtree_type kd,double *X,double *dist,int *ni,int*m,int *n,int *d,int *k);
192
void k_nn(double *X,double *dist,double *a,int *ni,int *n,int *d,int *k,int *get_a);
193 194 195 196 197 198
//void Rkdtree(double *X,int *n, int *d,int *idat,double *ddat);
SEXP Rkdtree(SEXP x);
//void Rkdnearest(double *X,int *idat,double *ddat,int *n,double *x, int *m, int *ni, double *dist,int *k);
SEXP Rkdnearest(SEXP Xr, SEXP xr,SEXP k);
//void Rkradius(double *r,int *idat,double *ddat,double *X,double *x,int *m,int *off,int *ni,int *op);
SEXP Rkradius(SEXP Xr, SEXP xr,SEXP rr,SEXP offr);
199 200 201 202
double xidist(double *x,double *X,int i,int d, int n);
int closest(kdtree_type *kd, double *X,double *x,int n,int *ex,int nex);
void kd_tree(double *X,int *n, int *d,kdtree_type *kd);
void free_kdtree(kdtree_type kd);
203

204 205 206 207 208 209 210 211
void tri2nei(int *t,int *nt,int *n,int *d,int *off);
void nei_penalty(double *X,int *n,int *d,double *D,int *ni,int *ii,int *off,
		 int *m,int *a_weight,double *kappa);
void sspl_construct(double *lambda,double *x,double *w,double *U,double *V,
             double *diagA,double *lb,int *n,double *tol);
void sspl_mapply(double *y,double *x,double *w,double *U,double *V,int *n,int *nf,double *tol,int *m);