--- author: Stéphane Laurent date: '2017-09-26' highlighter: 'pandoc-solarized' output: html_document: highlight: haddock keep_md: no md_document: preserve_yaml: True toc: no variant: markdown rbloggers: yes tags: 'haskell, R' title: Passing a R function to Haskell --- Passing R objects to Haskell ---------------------------- In two previous posts I have shown some examples of calling Haskell from R. More precisely, the procedure consists in building a DLL with Haskell and using this DLL in R, with the help of the .C function. We can obviously pass an integer, a double or a character string in the .C function. Thanks to the inline-r Haskell library, we can do more: namely, it is possible to pass any R object, since this library implements the type SEXP. Let's give an example. In this example we pass a R vector of doubles to Haskell, we calculate the square of each component in Haskell and we send the result to R.  {.haskell .numberLines} {-# LANGUAGE DataKinds #-} {-# LANGUAGE ForeignFunctionInterface #-} module Lib where import qualified Data.Vector.SEXP as VS import Foreign import Foreign.C import qualified Foreign.R.Type as R foreign export ccall squaredDoubles1 :: Ptr (SEXP s 'R.Real) -> Ptr (SEXP s 'R.Real) -> IO () squaredDoubles1 :: Ptr (SEXP s 'R.Real) -> Ptr (SEXP s 'R.Real) -> IO () squaredDoubles1 input result = do input <- peek input let inputAsList = (VS.toList . VS.fromSEXP) input let outputAsList = map (\x -> x*x) inputAsList let output = (VS.toSEXP . VS.fromList) outputAsList :: SEXP s 'R.Real poke result output  To call in R with the .C function, the R objects must be encapsulated in list():  {.r} > .C("squaredDoubles1", input = list(c(1,2,3)), result=list(0))$result[[1]] [1] 1 4 9  Instead of using VS.toList . VS.fromSEXP to convert the R vector to a Haskell list, we could use the real function of the Foreign.R module (this is a port of the C function REAL):  {.haskell .numberLines} ... import qualified Foreign.R as FR foreign export ccall squaredDoubles2 :: Ptr (SEXP s 'R.Real) -> Ptr (SEXP s 'R.Real) -> IO () squaredDoubles2 :: Ptr (SEXP s 'R.Real) -> Ptr (SEXP s 'R.Real) -> IO () squaredDoubles2 input result = do input <- peek input inputAsListPtr <- FR.real input l <- FR.length input inputAsList <- peekArray l inputAsListPtr let outputAsList = map (\x -> x*x) inputAsList let output = (VS.toSEXP . VS.fromList) outputAsList :: SEXP s 'R.Real poke result output  The performance is a bit better:  {.r} > library(microbenchmark) > x <- rnorm(100000) > microbenchmark( + H1 = .C("squaredDoubles1", input = list(x), result=list(0))$result[[1]], + H2 = .C("squaredDoubles2", input = list(x), result=list(0))$result[[1]] + ) Unit: milliseconds expr min lq mean median uq max neval cld H1 26.96348 34.70504 44.02896 38.77741 42.31139 205.2244 100 b H2 24.33826 30.25337 34.39467 32.80317 35.54754 160.4622 100 a  Alternatively, we can avoid the pointers and use the .Call function instead of the .C function:  {.haskell} foreign export ccall squaredDoubles3 :: SEXP s 'R.Real -> SEXP s 'R.Real squaredDoubles3 :: SEXP s 'R.Real -> SEXP s 'R.Real squaredDoubles3 input = (VS.toSEXP . VS.fromList) (map (\x -> x*x) ((VS.toList . VS.fromSEXP) input))   {.r} > .Call("squaredDoubles3", c(1,2,3)) [1] 1 4 9  More advanced usage: resorting to the FFI ----------------------------------------- Now we will show how to evaluate a R function. The function below is written in C. It takes as arguments a R function f (that is, a SEXP object of class CLOSXP), a double x, and it evaluates f(x). I'm using the C language and not inline-r for two reasons: - there's no port of the C functions allocSExp and defineVar in inline-r; - even if these two functions were available in Haskell (we could import them with the FFI), the Haskell code would be similar to the C code.  {.c .numberLines} #include #include double myeval(SEXP f, double x) { // convert x to SEXP SEXP xR; PROTECT(xR = allocVector(REALSXP, 1)); REAL(xR)[0] = x; UNPROTECT(1); // put f in an environment SEXP envir = allocSExp(ENVSXP); SEXP f_symbol = install("f"); defineVar(f_symbol, f, envir); // evaluate f(x) - like eval(call("f", x), envir) in R SEXP call = Rf_lang2(f_symbol, xR); return(REAL(eval(call, envir))[0]); }  Now we need to import this function:  {.haskell .numberLines} {-# LANGUAGE DataKinds #-} {-# LANGUAGE ForeignFunctionInterface #-} module Lib where import Foreign.C.Types import Foreign.R (SEXP, SEXP0, unsexp) import qualified Foreign.R as R import qualified Foreign.R.Type as R foreign import ccall unsafe "myeval" c_myeval :: SEXP0 -> CDouble -> CDouble myeval :: SEXP s 'R.Closure -> Double -> Double myeval f x = realToFrac (c_myeval (unsexp f) (realToFrac x))  Let us try it. The numerous realToFrac's could seem silly but for a more serious application we prefer the signature SEXP s 'R.Closure -> Double -> Double rather than SEXP s 'R.Closure -> CDouble -> CDouble.  {.haskell} foreign export ccall myevalR :: Ptr (SEXP s 'R.Closure) -> Ptr CDouble -> Ptr CDouble -> IO () myevalR :: Ptr (SEXP s 'R.Closure) -> Ptr CDouble -> Ptr CDouble -> IO () myevalR f x result = do f <- peek f x <- peek x poke result$ realToFrac $myeval f (realToFrac x :: Double)   {.r} > .C("myevalR", f=list(function(x) x+1), x=3, result=0)$result [1] 4  Thus, myeval f is a Haskell function of signature Double -> Double, though the evaluation is not performed by Haskell. Let us see an example of application. Form R, we will call the function  {.haskell} chebyshevFit :: Int -> (Double -> Double) -> [Double]  of the polynomial library.  {.haskell .numberLines} ... import Math.Polynomial.Chebyshev foreign export ccall chebyshevFitR :: Ptr (SEXP s 'R.Closure) -> Ptr CInt -> Ptr (SEXP V 'R.Real) -> IO () chebyshevFitR :: Ptr (SEXP s 'R.Closure) -> Ptr CInt -> Ptr (SEXP V 'R.Real) -> IO () chebyshevFitR f n result = do n <- peek n f <- peek f let fit = chebyshevFit (fromIntegral n :: Int) (myeval f) poke result $(VS.toSEXP . VS.fromList) fit  We will apply it to the function$x \mapsto \cos(4\arccos(x))$, which is the Chebyshev polynomial of order$4$for$|x| \leq 1$. Therefore, for any$n \geq 5$, the result must theoretically be$0, 0, 0, 0, 1, 0, \ldots, 0$.  {.r} > f <- function(x) cos(4*acos(x)) > .C("chebyshevFitR", f=list(f), n=6L, result=list(0))$result[[1]] [1] -1.110223e-16 3.145632e-16 -1.480297e-16 4.255855e-16 1.000000e+00 2.775558e-16