---
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 <R.h>
#include <Rinternals.h>

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
```