pomp
Inference for partially observed Markov processes
All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
bspline.c File Reference
#include "internal.h"
Include dependency graph for bspline.c:

Go to the source code of this file.

Functions

static void bspline_eval (double *y, const double *x, int nx, int i, int degree, int deriv, const double *knots)
 
SEXP bspline_basis (SEXP range, SEXP x, SEXP nbasis, SEXP degree, SEXP deriv)
 
SEXP periodic_bspline_basis (SEXP x, SEXP nbasis, SEXP degree, SEXP period, SEXP deriv)
 
void bspline_basis_eval_deriv (double x, double *knots, int degree, int nbasis, int deriv, double *y)
 
void periodic_bspline_basis_eval_deriv (double x, double period, int degree, int nbasis, int deriv, double *y)
 

Function Documentation

◆ bspline_basis()

SEXP bspline_basis ( SEXP  range,
SEXP  x,
SEXP  nbasis,
SEXP  degree,
SEXP  deriv 
)

Definition at line 51 of file bspline.c.

51  {
52  SEXP y, xr;
53  int nx = LENGTH(x);
54  int nb = INTEGER_VALUE(nbasis);
55  int deg = INTEGER_VALUE(degree);
56  int d = INTEGER_VALUE(deriv);
57  int nk = nb+deg+1;
58  double dx, minx, maxx;
59  double *knots = NULL;
60  double *xdata, *ydata;
61  int i;
62  if (deg < 0) err("must have degree > 0");
63  if (nb <= deg) err("must have nbasis > degree");
64  if (d < 0) err("must have deriv >= 0");
65  knots = (double *) R_Calloc(nk,double);
66  PROTECT(xr = AS_NUMERIC(x));
67  PROTECT(y = allocMatrix(REALSXP,nx,nb));
68  xdata = REAL(xr);
69  ydata = REAL(y);
70  minx = REAL(range)[0];
71  maxx = REAL(range)[1];
72  if (minx >= maxx) err("improper range 'rg'");
73  dx = (maxx-minx)/((double) (nb-deg));
74  knots[0] = minx-deg*dx;
75  for (i = 1; i < nk; i++) knots[i] = knots[i-1]+dx;
76  for (i = 0; i < nb; i++) {
77  bspline_eval(ydata,xdata,nx,i,deg,d,&knots[0]);
78  ydata += nx;
79  }
80  R_Free(knots);
81  UNPROTECT(2);
82  return(y);
83 }
static void bspline_eval(double *y, const double *x, int nx, int i, int degree, int deriv, const double *knots)
Definition: bspline.c:10
#define err(...)
Definition: pomp.h:21
Here is the call graph for this function:

◆ bspline_basis_eval_deriv()

void bspline_basis_eval_deriv ( double  x,
double *  knots,
int  degree,
int  nbasis,
int  deriv,
double *  y 
)

Definition at line 115 of file bspline.c.

116  {
117  for (int i = 0; i < nbasis; i++) {
118  bspline_eval(&y[i],&x,1,i,degree,deriv,knots);
119  }
120 }
Here is the call graph for this function:
Here is the caller graph for this function:

◆ bspline_eval()

static void bspline_eval ( double *  y,
const double *  x,
int  nx,
int  i,
int  degree,
int  deriv,
const double *  knots 
)
static

Definition at line 9 of file bspline.c.

12 {
13  int j;
14  if (deriv > degree) {
15  for (j = 0; j < nx; j++) y[j] = 0.0;
16  } else if (deriv > 0) {
17  int i2 = i+1, p2 = degree-1, d2 = deriv-1;
18  double *y1 = (double *) R_Calloc(nx,double);
19  double *y2 = (double *) R_Calloc(nx,double);
20  double a, b;
21  bspline_eval(y1,x,nx,i,p2,d2,knots);
22  bspline_eval(y2,x,nx,i2,p2,d2,knots);
23  for (j = 0; j < nx; j++) {
24  a = degree / (knots[i+degree]-knots[i]);
25  b = degree / (knots[i2+degree]-knots[i2]);
26  y[j] = a * y1[j] - b * y2[j];
27  }
28  R_Free(y1); R_Free(y2);
29  } else { // deriv == 0
30  int i2 = i+1, p2 = degree-1;
31  if (degree > 0) { // case deriv < degree
32  double *y1 = (double *) R_Calloc(nx,double);
33  double *y2 = (double *) R_Calloc(nx,double);
34  double a, b;
35  bspline_eval(y1,x,nx,i,p2,deriv,knots);
36  bspline_eval(y2,x,nx,i2,p2,deriv,knots);
37  for (j = 0; j < nx; j++) {
38  a = (x[j]-knots[i]) / (knots[i+degree]-knots[i]);
39  b = (knots[i2+degree]-x[j]) / (knots[i2+degree]-knots[i2]);
40  y[j] = a * y1[j] + b * y2[j];
41  }
42  R_Free(y1); R_Free(y2);
43  } else if (degree == 0) { // case deriv == degree
44  for (j = 0; j < nx; j++)
45  y[j] = (double) ((knots[i] <= x[j]) && (x[j] < knots[i2]));
46  }
47  }
48 }
Here is the caller graph for this function:

◆ periodic_bspline_basis()

SEXP periodic_bspline_basis ( SEXP  x,
SEXP  nbasis,
SEXP  degree,
SEXP  period,
SEXP  deriv 
)

Definition at line 85 of file bspline.c.

86  {
87 
88  SEXP y, xr;
89  int nx = LENGTH(x);
90  int nb = INTEGER_VALUE(nbasis);
91  int deg = INTEGER_VALUE(degree);
92  int d = INTEGER_VALUE(deriv);
93  double pd = NUMERIC_VALUE(period);
94  int j, k;
95  double *xrd, *ydata, *val;
96  PROTECT(xr = AS_NUMERIC(x));
97  xrd = REAL(xr);
98  PROTECT(y = allocMatrix(REALSXP,nx,nb));
99  ydata = REAL(y);
100  val = (double *) R_alloc(nb,sizeof(double));
101  for (j = 0; j < nx; j++) {
102  periodic_bspline_basis_eval_deriv(xrd[j],pd,deg,nb,d,val);
103  for (k = 0; k < nb; k++) ydata[j+nx*k] = val[k];
104  }
105  UNPROTECT(2);
106  return y;
107 }
void periodic_bspline_basis_eval_deriv(double x, double period, int degree, int nbasis, int deriv, double *y)
Definition: bspline.c:122
Here is the call graph for this function:

◆ periodic_bspline_basis_eval_deriv()

void periodic_bspline_basis_eval_deriv ( double  x,
double  period,
int  degree,
int  nbasis,
int  deriv,
double *  y 
)

Definition at line 122 of file bspline.c.

124 {
125  int nknots = nbasis+2*degree+1;
126  int shift = (degree-1)/2;
127  double *knots = NULL, *yy = NULL;
128  double dx;
129  int j, k;
130 
131  if (period <= 0.0) err("must have period > 0");
132  if (nbasis <= 0) err("must have nbasis > 0");
133  if (degree < 0) err("must have degree >= 0");
134  if (nbasis < degree) err("must have nbasis >= degree");
135  if (deriv < 0) err("must have deriv >= 0");
136 
137  knots = (double *) R_Calloc(nknots+degree+1,double);
138  yy = (double *) R_Calloc(nknots,double);
139 
140  dx = period/((double) nbasis);
141  for (k = -degree; k <= nbasis+degree; k++) {
142  knots[degree+k] = k*dx;
143  }
144  x = fmod(x,period);
145  if (x < 0.0) x += period;
146  for (k = 0; k < nknots; k++) {
147  bspline_eval(&yy[k],&x,1,k,degree,deriv,knots);
148  }
149  for (k = 0; k < degree; k++) yy[k] += yy[nbasis+k];
150  for (k = 0; k < nbasis; k++) {
151  j = (shift+k)%nbasis;
152  y[k] = yy[j];
153  }
154  R_Free(yy); R_Free(knots);
155 }
Here is the call graph for this function:
Here is the caller graph for this function: