Skip to content

Commit

Permalink
use perlre for template trigger matching;
Browse files Browse the repository at this point in the history
Csnippet interface to ordinary B spline bases
  • Loading branch information
kingaa committed Nov 30, 2023
1 parent 2782c5c commit 658bb34
Show file tree
Hide file tree
Showing 22 changed files with 151 additions and 30 deletions.
4 changes: 2 additions & 2 deletions DESCRIPTION
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical Inference for Partially Observed Markov Processes
Version: 5.4.3.0
Date: 2023-11-26
Version: 5.4.3.1
Date: 2023-11-30
Authors@R: c(person(given=c("Aaron","A."),family="King",role=c("aut","cre"),email="[email protected]",comment=c(ORCID="0000-0001-6159-3207")),
person(given=c("Edward","L."),family="Ionides",role="aut",comment=c(ORCID="0000-0002-4190-0174")) ,
person(given="Carles",family="Bretó",role="aut",comment=c(ORCID="0000-0003-4695-4902")),
Expand Down
3 changes: 2 additions & 1 deletion R/builder.R
Original file line number Diff line number Diff line change
Expand Up @@ -197,7 +197,8 @@ Cbuilder <- function (..., templates, name = NULL, dir = NULL,
\(x) any(
grepl(
pomp_templates$utilities[[x]]$trigger,
c(snippets,on_load)
c(snippets,on_load),
perl=TRUE
)
),
logical(1L)
Expand Down
10 changes: 9 additions & 1 deletion R/templates.R
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,14 @@ static periodic_bspline_basis_eval_t *__pomp_periodic_bspline_basis_eval;
#define periodic_bspline_basis_eval(X,Y,M,N,Z) (__pomp_periodic_bspline_basis_eval((X),(Y),(M),(N),(Z)))}",
reg=r"{
__pomp_periodic_bspline_basis_eval = (periodic_bspline_basis_eval_t *) R_GetCCallable("pomp","periodic_bspline_basis_eval");}"
),
bspline_basis=list(
trigger=r"{(?<!periodic_)bspline_basis_eval}",
header=r"{
static bspline_basis_eval_t *__pomp_bspline_basis_eval;
#define bspline_basis_eval(X,K,D,N,Y) (__pomp_bspline_basis_eval((X),(K),(D),(N),(Y)))}",
reg=r"{
__pomp_bspline_basis_eval = (bspline_basis_eval_t *) R_GetCCallable("pomp","bspline_basis_eval");}"
),
get_userdata_int=list(
trigger="get_userdata_int",
Expand All @@ -43,7 +51,7 @@ static get_userdata_double_t *__pomp_get_userdata_double;
__pomp_get_userdata_double = (get_userdata_double_t *) R_GetCCallable("pomp","get_userdata_double");}"
),
get_userdata=list(
trigger=r"{get_userdata(\b|[^_])}",
trigger=r"{get_userdata(?!_)}",
header=r"{
static get_userdata_t *__pomp_get_userdata;
#define get_userdata(X) (__pomp_get_userdata(X))}",
Expand Down
12 changes: 10 additions & 2 deletions inst/NEWS
Original file line number Diff line number Diff line change
Expand Up @@ -2,8 +2,16 @@ _N_e_w_s _f_o_r _p_a_c_k_a_g_e '_p_o_m_p'

_C_h_a_n_g_e_s _i_n '_p_o_m_p' _v_e_r_s_i_o_n _5._4._3:

• Modify formatting of printed messages to forestall warnings
in latest R-devel.
• There is now access to ordinary B-spline basis functions at
the C snippet level. It will still almost always be
preferable to construct a spline basis once and pass it to
‘pomp’ functions as part of a ‘covariate_table’, but this
functionality may occasionally be useful.

• ‘bspline_eval’ is now deprecated as part of the ‘pomp’ C API.

• The formatting of some printed messages has been modified to
forestall warnings in the latest R-devel.

_C_h_a_n_g_e_s _i_n '_p_o_m_p' _v_e_r_s_i_o_n _5._4._2:

Expand Down
5 changes: 4 additions & 1 deletion inst/NEWS.Rd
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,10 @@
\title{News for package `pomp'}
\section{Changes in \pkg{pomp} version 5.4.3}{
\itemize{
\item Modify formatting of printed messages to forestall warnings in latest R-devel.
\item There is now access to ordinary B-spline basis functions at the C snippet level.
It will still almost always be preferable to construct a spline basis once and pass it to \pkg{pomp} functions as part of a \code{covariate_table}, but this functionality may occasionally be useful.
\item \code{bspline_eval} is now deprecated as part of the \pkg{pomp} C API.
\item The formatting of some printed messages has been modified to forestall warnings in the latest R-devel.
}
}
\section{Changes in \pkg{pomp} version 5.4.2}{
Expand Down
11 changes: 6 additions & 5 deletions inst/include/pomp.h
Original file line number Diff line number Diff line change
Expand Up @@ -21,13 +21,14 @@
#define err(...) errorcall(R_NilValue,__VA_ARGS__)
#define warn(...) warningcall(R_NilValue,__VA_ARGS__)

typedef void bspline_eval_t (double *y, const double *x, int nx, int i, int p,
int d, const double *knots);
typedef void bspline_basis_eval_t (double x, double *knots, int degree,
int nbasis, double *y);
typedef void periodic_bspline_basis_eval_t (double x, double period, int degree,
int nbasis, double *y);
typedef void periodic_bspline_basis_eval_deriv_t (double x, double period,
int degree, int nbasis,
int deriv, double *y);

// The following is deprecated and will be removed from the API.
typedef void bspline_eval_t (double *y, const double *x, int nx, int i, int p,
int d, const double *knots);

typedef const SEXP get_userdata_t (const char *name);
typedef const int *get_userdata_int_t (const char *name);
Expand Down
22 changes: 14 additions & 8 deletions src/bspline.c
Original file line number Diff line number Diff line change
Expand Up @@ -4,10 +4,10 @@

// The following function computes the derivative of order 'deriv' of the i-th
// B-spline of given degree with given knots at each of the nx points in x.
// knots must point to an array of length nbasis+p+1
// knots must point to an array of length nbasis+degree+1
// The results are stored in y.
void bspline_eval (double *y, const double *x, int nx, int i,
int degree, int deriv, const double *knots)
int degree, int deriv, const double *knots)
{
int j;
if (deriv > degree) {
Expand Down Expand Up @@ -105,15 +105,21 @@ SEXP periodic_bspline_basis (SEXP x, SEXP nbasis, SEXP degree, SEXP period,
return y;
}

/*
void bspline_basis_eval (double x, int degree, int nbasis, double *y) {
// The following function computes value at x of all the nbasis
// B-spline basis functions of the specified degree with the given knots.
// 'knots' must point to an array of length nbasis+degree+1
// The results are stored in y.
void bspline_basis_eval (double x, double *knots, int degree,
int nbasis, double *y) {
bspline_basis_eval_deriv(x,knots,degree,nbasis,0,y);
}

void bspline_basis_eval_deriv (double x, double period, int degree,
int nbasis, int deriv, double *y)
{}
*/
void bspline_basis_eval_deriv (double x, double *knots, int degree,
int nbasis, int deriv, double *y) {
for (int i = 0; i < nbasis; i++) {
bspline_eval(&y[i],&x,1,i,degree,deriv,knots);
}
}

void periodic_bspline_basis_eval (double x, double period, int degree,
int nbasis, double *y) {
Expand Down
2 changes: 2 additions & 0 deletions src/decls.h
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,8 @@
extern void bspline_eval(double *y, const double *x, int nx, int i, int degree, int deriv, const double *knots);
extern SEXP bspline_basis(SEXP range, SEXP x, SEXP nbasis, SEXP degree, SEXP deriv);
extern SEXP periodic_bspline_basis(SEXP x, SEXP nbasis, SEXP degree, SEXP period, SEXP deriv);
extern void bspline_basis_eval(double x, double *knots, int degree, int nbasis, double *y);
extern void bspline_basis_eval_deriv(double x, double *knots, int degree, int nbasis, int deriv, double *y);
extern void periodic_bspline_basis_eval(double x, double period, int degree, int nbasis, double *y);
extern void periodic_bspline_basis_eval_deriv(double x, double period, int degree, int nbasis, int deriv, double *y);
/* src/dinit.c */
Expand Down
1 change: 1 addition & 0 deletions src/init.c
Original file line number Diff line number Diff line change
Expand Up @@ -52,6 +52,7 @@ static const R_CallMethodDef callMethods[] = {

void R_init_pomp (DllInfo *info) {
// C functions provided for users
R_RegisterCCallable("pomp","bspline_basis_eval",(DL_FUNC) &bspline_basis_eval);
R_RegisterCCallable("pomp","bspline_eval",(DL_FUNC) &bspline_eval);
R_RegisterCCallable("pomp","periodic_bspline_basis_eval",(DL_FUNC) &periodic_bspline_basis_eval);
R_RegisterCCallable("pomp","get_userdata",(DL_FUNC) &get_userdata);
Expand Down
11 changes: 6 additions & 5 deletions src/pomp.h
Original file line number Diff line number Diff line change
Expand Up @@ -21,13 +21,14 @@
#define err(...) errorcall(R_NilValue,__VA_ARGS__)
#define warn(...) warningcall(R_NilValue,__VA_ARGS__)

typedef void bspline_eval_t (double *y, const double *x, int nx, int i, int p,
int d, const double *knots);
typedef void bspline_basis_eval_t (double x, double *knots, int degree,
int nbasis, double *y);
typedef void periodic_bspline_basis_eval_t (double x, double period, int degree,
int nbasis, double *y);
typedef void periodic_bspline_basis_eval_deriv_t (double x, double period,
int degree, int nbasis,
int deriv, double *y);

// The following is deprecated and will be removed from the API.
typedef void bspline_eval_t (double *y, const double *x, int nx, int i, int p,
int d, const double *knots);

typedef const SEXP get_userdata_t (const char *name);
typedef const int *get_userdata_int_t (const char *name);
Expand Down
File renamed without changes
File renamed without changes
File renamed without changes
File renamed without changes
File renamed without changes
File renamed without changes
File renamed without changes
2 changes: 1 addition & 1 deletion tests/bsplines.R → tests/bsplines1.R
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
options(digits=3)
png(filename="bsplines-%02d.png",res=100)
png(filename="bsplines1-%02d.png",res=100)
library(pomp)

x <- seq(0,2,by=0.01)
Expand Down
6 changes: 2 additions & 4 deletions tests/bsplines.Rout.save → tests/bsplines1.Rout.save
Original file line number Diff line number Diff line change
@@ -1,14 +1,12 @@

R version 4.3.1 (2023-06-16) -- "Beagle Scouts"
R version 4.3.2 (2023-10-31) -- "Eye Holes"
Copyright (C) 2023 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Expand All @@ -18,7 +16,7 @@ Type 'demo()' for some demos, 'help()' for on-line help, or
Type 'q()' to quit R.

> options(digits=3)
> png(filename="bsplines-%02d.png",res=100)
> png(filename="bsplines1-%02d.png",res=100)
> library(pomp)

Welcome to pomp!
Expand Down
Binary file added tests/bsplines2-01.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
32 changes: 32 additions & 0 deletions tests/bsplines2.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,32 @@
options(digits=3)
png(filename="bsplines2-%02d.png",res=100)
library(pomp)

trajectory(
t0=0, times=seq(0,2,by=0.02),
skeleton=vectorfield(Csnippet(r"{
const double *a = &a1;
double knots[] = {-3, -2, -1, 0, 1, 2, 3, 4, 5};
double s[5], p[5];
bspline_basis_eval(t,knots,3,5,s);
periodic_bspline_basis_eval(t,1,3,5,p);
Dx = dot_product(5,a,s);
Dy = dot_product(5,a,&trend_1);
Dz = dot_product(5,a,p);
Dw = dot_product(5,a,&seas_1);
}")),
params = c(
a1=-1,a2=1,a3=0,a4=-3,a5=18,
x_0=0,y_0=0,z_0=0,w_0=0
),
paramnames=c("a1","a2","a3"),
statenames=c("x","y","w","z"),
covar=covariate_table(
times=seq(0,2.1,by=0.001),
trend=bspline_basis(x=times,nbasis=5,degree=3,rg=c(0,2)),
seas=periodic_bspline_basis(x=times,period=1,nbasis=5,degree=3)
)
) |>
plot()

dev.off()
60 changes: 60 additions & 0 deletions tests/bsplines2.Rout.save
Original file line number Diff line number Diff line change
@@ -0,0 +1,60 @@

R version 4.3.2 (2023-10-31) -- "Eye Holes"
Copyright (C) 2023 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> options(digits=3)
> png(filename="bsplines2-%02d.png",res=100)
> library(pomp)

Welcome to pomp!

As of version 4.6, no user-visible pomp function has a name that
includes a dot ('.'). Function names have been changed to replace the
dot with an underscore ('_'). For more information, see the pomp blog:
https://kingaa.github.io/pomp/blog.html.

>
> trajectory(
+ t0=0, times=seq(0,2,by=0.02),
+ skeleton=vectorfield(Csnippet(r"{
+ const double *a = &a1;
+ double knots[] = {-3, -2, -1, 0, 1, 2, 3, 4, 5};
+ double s[5], p[5];
+ bspline_basis_eval(t,knots,3,5,s);
+ periodic_bspline_basis_eval(t,1,3,5,p);
+ Dx = dot_product(5,a,s);
+ Dy = dot_product(5,a,&trend_1);
+ Dz = dot_product(5,a,p);
+ Dw = dot_product(5,a,&seas_1);
+ }")),
+ params = c(
+ a1=-1,a2=1,a3=0,a4=-3,a5=18,
+ x_0=0,y_0=0,z_0=0,w_0=0
+ ),
+ paramnames=c("a1","a2","a3"),
+ statenames=c("x","y","w","z"),
+ covar=covariate_table(
+ times=seq(0,2.1,by=0.001),
+ trend=bspline_basis(x=times,nbasis=5,degree=3,rg=c(0,2)),
+ seas=periodic_bspline_basis(x=times,period=1,nbasis=5,degree=3)
+ )
+ ) |>
+ plot()
>
> dev.off()
null device
1
>

0 comments on commit 658bb34

Please sign in to comment.