How to make a C library accessible to the Csnippet

[MarieAugerMethe]

Hi again,

I'm trying to sample from a truncated normal distribution and found a C library that does just that. The source code is available at:
https://people.sc.fsu.edu/~jburkardt/c_src/truncated_normal/truncated_normal.c
https://people.sc.fsu.edu/~jburkardt/c_src/truncated_normal/truncated_normal.h

But I'm not sure how to make that C library available to my Csnippet.

I've tried a few things. For example, giving the the name of the file in the pomp argument shlib.args and having a line in a global C function that mentions to include that library. Also tried compiling it with r cmd shlib. Nothing appears to work. I can provide a specific example if that's helpful, but I was wondering first whether there was a specific way in which I should try to include a C library. I feel like I'm so off base that my reproducible example might not be useful. If possible, it would be great if that way was done directly from R, so I can distribute the R code easily. Any tips would be great!

Many thanks,

Marie

[kingaa]

@MarieAugerMethe:

Sorry to be so slow in replying. Here is an example I cobbled together to demonstrate two different ways of accomplishing this.

The first compiles the auxiliary C code into an object file, "truncated_normal.o", and then links this with the pomp-generated C snippet codes. Note the use of the globals argument, which puts the #include line into the C snippet file, and the shlib.args argument, which tells pomp to tell R CMD SHLIB to link with "truncated_normal.o". Note also that I've here specified the C snippet file name (via the cfile parameter) and the directory in which the compilation will take place (via the cdir parameter).

library(pomp)
library(magrittr)
stopifnot(packageVersion("pomp")>="1.16.4.1")

pompExample(ricker)

## newly created object(s):
##  ricker

system("R CMD COMPILE truncated_normal.c")

set.seed(2300595)
pomp(ricker,rmeasure=Csnippet("
        y = rnorm(N,1);
        y = (y > 0) ? y : 0;
        "),
     dmeasure=Csnippet("
        lik = truncated_normal_a_pdf(y,N,1,0);"),
     statenames="N",
     globals="#include \"truncated_normal.h\"",
     shlib.args="truncated_normal.o",
     cdir=getwd(),cfile="tnrick") %>%
  simulate() -> rick1

# test the 'dmeasure':
replicate(10,pfilter(rick1,Np=1000) %>% logLik()) %>% logmeanexp(se=TRUE)

##                      se 
## -80.2869566   0.1048461

Note that it isn't necessary to first compile "truncated_normal.c", and it is possible to do everything in the default place in the temporary directory. In this case, since it uses make, R CMD SHLIB is smart enough to figure out that it needs to compile "truncated_normal.c" to get "truncated_normal.o". However, I do have to put the path to the header file "truncated_normal.h" explicitly into the C snippet file so that the preprocessor can find the header file.

set.seed(2300595)
pomp(ricker,rmeasure=Csnippet("
        y = rnorm(N,1);
        y = (y > 0) ? y : 0;
        "),
  dmeasure=Csnippet("
        lik = truncated_normal_a_pdf(y,N,1,0);"),
  statenames="N",
  globals=sprintf("#include \"%s/truncated_normal.h\"",getwd()),
  shlib.args="truncated_normal.o") %>%
  simulate() -> rick2

# test the 'dmeasure':
replicate(10,pfilter(rick2,Np=1000) %>% logLik()) %>% logmeanexp(se=TRUE)

##                      se 
## -80.2869566   0.1048461

The second method just inserts the entire contents of the header file and the C file into the pomp-generated C snippet file via the globals argument. Obviously, this works fine here but might be suboptimal with very large codes.

Note that the first two lines read the .c and .h files into the R session. The third line actually inserts the contents of the .h file into the appropriate place in the .c file (which is what the C preprocessor normally does). The fourth line just concatenates all the lines into one long string.

readLines("truncated_normal.h") %>% paste(collapse="\n") -> tn.h
readLines("truncated_normal.c") -> tn.c
tn.c[grepl("include \"truncated_normal.h\"",tn.c)] <- tn.h
extract(tn.c,!grepl("include \"truncated_normal.h\"",tn.c)) %>%
  paste(collapse="\n") -> tn.c

set.seed(2300595)
pomp(ricker,rmeasure=Csnippet("
        y = rnorm(N,1);
        y = (y > 0) ? y : 0;
        "),
  dmeasure=Csnippet("
        lik = truncated_normal_a_pdf(y,N,1,0);"),
  statenames="N",
  globals=tn.c) %>%
  simulate() -> rick3

# test the 'dmeasure':
replicate(10,pfilter(rick3,Np=1000) %>% logLik()) %>% logmeanexp(se=TRUE)

##                      se 
## -80.2869566   0.1048461

I hope this helps. Let me know if it does not, or if you have more questions.

I want to thank you for bringing the question up: it will help make a useful entry to the FAQ.

To make this issue self-contained, I've upload the .c and .h files here:

truncated_normal.c
truncated_normal.h

(Somewhat annoyingly, they're treated as .txt files, so after downloading, the .txt file extension will have to be removed manually.)

[MarieAugerMethe]

@kingaa Fantastic, thank you! All three of these work.

I'm not sure I fully understand why the first option has to create the .o file and then the specific .so .c files for pomp... but it works. I'm also not fully clear on the third option... but I think it's just because I don't fully understand C.

I've created two completely artificial examples based on your two first examples that add more to the global file (because in my real model I have pretty extensive global files). It's pretty straightforward, but it took me a few minutes to figure it out. So here it is for posterity.

Thank you!

Marie

library(pomp)
library(magrittr)
stopifnot(packageVersion("pomp")>="1.16.4.1")

pompExample(ricker)

system("R CMD COMPILE truncated_normal.c")

set.seed(2300595)

rmeas <- Csnippet("
                  y = rnorm(N,1);
                  y = (y > 0) ? y : 0;
                  testt();
                  ")

dmeas <- Csnippet("lik = truncated_normal_a_pdf(y,N,1,0);")

c_global_cp <- ('
#include \"truncated_normal.h\"
#include<stdio.h>
void testt(){
    printf("It is a test...\\n");
}

')

pomp(ricker,
     rmeasure=rmeas,
     dmeasure=dmeas,
     statenames="N",
     #globals="#include \"truncated_normal.h\"",
     globals=c_global_cp,
     shlib.args="truncated_normal.o",
     cdir=getwd(),cfile="tnrick") %>%
  simulate() -> rick1

# test the 'dmeasure':
replicate(10,pfilter(rick1,Np=1000) %>% logLik()) %>% logmeanexp(se=TRUE)

#c_global_cp_2 <- sprintf("#include \"%s/truncated_normal.h\"",getwd())
c_global_cp_2 <- paste( 
sprintf("#include \"%s/truncated_normal.h\"",getwd()),
('
#include<stdio.h>
void testt(){
 printf("It is a test...\\n");
 }
'))

set.seed(2300595)
pomp(ricker,rmeasure= rmeas,
     dmeasure=dmeas,
     statenames="N",
     globals=c_global_cp_2,
     shlib.args="truncated_normal.o") %>%
  simulate() -> rick2

# test the 'dmeasure':
replicate(10,pfilter(rick2,Np=1000) %>% logLik()) %>% logmeanexp(se=TRUE)

[kingaa]

Am re-opening this to remind me to use this material in an FAQ entry. Thanks again, @MarieAugerMethe!

[chvandorp]

Hi,
I have a similar question. I would like to use a c++ library instead of c. Is this at all possible?
Thanks!

[kingaa]

But of course!

[kingaa]

The following updates the codes above for pomp versions ≥ 4. The codes below link the truncated_normal library using C snippets. This corresponds to my option two above. Note that things have gotten a bit simpler by version 4 of pomp.

ricker() |>
  simulate(
    rmeasure=Csnippet("
        y = rnorm(N,1);
        y = (y > 0) ? y : 0;
        "),
    dmeasure=Csnippet("
        lik = truncated_normal_a_pdf(y,N,1,0);
        if (give_log) lik = log(lik);"),
    statenames="N",
    globals="#include \"truncated_normal.h\"",
    shlib.args="truncated_normal.o"
) -> rick2

replicate(n=10,
  rick2 |>
    pfilter(Np=1000) |>
    logLik()
) |>
  logmeanexp(se=TRUE)

The two files mentioned (truncated_normal.h, truncated_normal.c) are included in this zip file: truncated_normal.zip