I recently ran into the problem that Julias’s (GNU Scientific Library) GSL.jl library is too new, i.e. some portions don’t work as of early 2016. However, to solve my problem I needed access to adaptive Monte Carlo integration methods. This means it was time to go in depth in Julia’s C interface. I will go step by step on how I created and compiled the C code, and called it from Julia.

What I wished to use was the GSL Vegas functions. Luckily, there is a good example on this page for how to use the library. I started with this code. I then modified it a bit. To start, I changed the main function header to:

int monte_carlo_integrate (double (*integrand)(double *,size_t,void *),double* retRes,double* retErr,int mode,int dim,double* xl,double* xu,size_t calls){

Before it took in no values, but I will want to be able to do most things from Julia. First of all, I changed the function name from main to avoid namespace issues. This is just good practice when working with shared libraries. then I added a bunch of arguments. The first argument is a function pointer, where the function is defined just as g is in the example. Then I pass in pointers which we become the return values. I add mode (1,2,3 are the different integration types), dim, xl, xu, and calls as passed in values to easily extend the functionality. The rest of the C-code stays mostly the same: I change the parts which reference g to integrand (I like the function name better) and comment out the print functions (they tend to no print correctly). So I end up with a C-function as follows:

#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_plain.h>
#include <gsl/gsl_monte_miser.h>
#include <gsl/gsl_monte_vegas.h>
int monte_carlo_integrate (double (*integrand)(double *,size_t,void *),double* retRes,double* retErr,int mode,int dim,double* xl,double* xu,size_t calls){
  double res, err;
  const gsl_rng_type *T;
  gsl_rng *r;
  gsl_monte_function G = { *integrand, dim, 0 };
  gsl_rng_env_setup ();
  T = gsl_rng_default;
  r = gsl_rng_alloc (T);
  if (mode==1){
    gsl_monte_plain_state *s = gsl_monte_plain_alloc (dim);
    gsl_monte_plain_integrate (&G, xl, xu, dim, calls, r, s,
                               &res, &err);
    gsl_monte_plain_free (s);
    /* display_results ("plain", res, err); */
  }
  if (mode==2){
    gsl_monte_miser_state *s = gsl_monte_miser_alloc (dim);
    gsl_monte_miser_integrate (&G, xl, xu, dim, calls, r, s,
                               &res, &err);
    gsl_monte_miser_free (s);
    /* display_results ("miser", res, err); */
  }
  if (mode==3){
    gsl_monte_vegas_state *s = gsl_monte_vegas_alloc (dim);
    gsl_monte_vegas_integrate (&G, xl, xu, dim, 10000, r, s,
                               &res, &err);
    /* display_results ("vegas warm-up", res, err); */
    /* printf ("converging...\n"); */
    do
      {
        gsl_monte_vegas_integrate (&G, xl, xu, dim, calls/5, r, s,
                                   &res, &err);
        /* printf ("result = % .6f sigma = % .6f "
                "chisq/dof = %.1f\n", res, err, gsl_monte_vegas_chisq (s)); */
      }
    while (fabs (gsl_monte_vegas_chisq (s) - 1.0) > 0.5);
    /* display_results ("vegas final", res, err); */
    gsl_monte_vegas_free (s);
  }
  gsl_rng_free (r);
  *retRes = res;
  *retErr = err;
  return 0;
}

Notice how nothing was specifically changed to be “Julia-style”, this is all normal C-code. So how do I use it? First I need to compile it to a shared library. This means I use the -shared tag and save it to a .so. I will need the GSL libraries loaded when compiling, so I add the -lgsl, -lgslcblas, and -lm (math) tags to link those libraries. Then I add the -fPIC tag since Julia’s documentation says so, and tell it the file to compile. The complete code is:

gcc -Wall -shared -o libMonte.so -lgsl -lgslcblas -lm -fPIC monteCarloExample.c

This will give you a .so file. Now we need to use it. I will describe passing in the function last. Here’s the setup. The first of the other variables are the two pointers retRes and retErr where the results will be saved. In Julia, to get a pointer, I make two 1-element arrays as follows:

x = Array{Float64}(1)
y = Array{Float64}(1)

The only other peculiar thing I had to do was to change pi from Julia’s irrational type to a Float64 (Cdouble) and use that to make the array. So for the other variables I used:

fpi = convert(Float64,pi)
xl = [0.; 0.; 0.]
xu = [fpi; fpi; fpi]
calls = 500000
mode = 3
dim = 3

Now I pass this all to C as follows:

ccall((:monte_carlo_integrate,"/path/to/library/libMonte.so"),Int32,(Ptr{Void},Ptr{Cdouble},Ptr{Cdouble},Int32,Int32,Ptr{Cdouble},Ptr{Cdouble},Csize_t),integrand_c,x,y,mode,dim,xl,xu,calls)

The first argument a tuple where the first place is the symbol which says which function in the library to use and the second place is the library name or the path to the library. The second argument is the return type (here it’s Int32 because our function returns int 0 when complete). The next portion is a tuple which defined the types of the arguments we are passing. The first is Ptr{Void} which is used for things like function pointers. Next there are two Cdouble pointers for retRes and retErr. Then two integers, two more pointers, and a Csize_t. Lastly we put in all of the variables we want to pass in.

Recall that the result was stored into the pointers for the second and third variable passed in. These are the pointers x and y. So to get the values of res and err, I de-reference them:

res=x[1]
err=y[1]

Now res and err hold the values for the solution and the error.

I am not done yet since I didn’t talk about the function! To do this, we define the function in Julia. We have to use parametric types or Julia will yet at us, and we have to ensure that the returned value is a Cdouble (in our case). So we re-write the g function in Julia as:

function integrand{T,T2,T3}(x::T,dim::T2,params::T3)
  A = 1.0 / (pi * pi * pi)
  return A / (1.0 - cos(unsafe_load(x,1))*cos(unsafe_load(x,2))*cos(unsafe_load(x,3)))::Cdouble
end

Notice that instead of x[1], we have to use the Julia function unsafe_load(x,1) to de-reference a pointer. However, since this part is in Julia, other things are much safer, like how we can use pi directly without having to convert it to a float. Also notice that we can add print statements in this function, and they will print directly to Julia. You can use this to modify the display_results function to be a Julia function which prints. However, this itself is still not able to be passed into C. To do that, you have to translate it to a C function:

integrand_c = cfunction(integrand,Cdouble,(Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}))

Here we used the Julia cfunction function where the first argument is the function, the second argument is what the function returns, and the third argument is a tuple of C-types that the function will take on. If you look back at the ccall, this integrand_c is what we passed as the first argument to the actual C-function.

Check to see that this all works. It worked for me. For completeness I will put the full Julia code below. Happy programming!

x = Array{Float64}(1)
y = Array{Float64}(1)
fpi = convert(Float64,pi)
xl = [0.; 0.; 0.]
xu = [fpi; fpi; fpi]
 
function integrand{T,T2,T3}(x::T,dim::T2,params::T3)
  A = 1.0 / (pi * pi * pi)
  return 3*A / (1.0 - cos(unsafe_load(x,1))*cos(unsafe_load(x,2))*cos(unsafe_load(x,3)))::Cdouble
end
 
integrand_c = cfunction(integrand,Cdouble,(Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}))
 
calls = 500000
mode = 3
dim = 3
ccall((:monte_carlo_integrate,"/home/crackauc/Public/libMonte.so"),Int32,(Ptr{Void},Ptr{Cdouble},Ptr{Cdouble},Int32,Int32,Ptr{Cdouble},Ptr{Cdouble},Csize_t),integrand_c,x,y,mode,dim,xl,xu,calls)
res=x[1]
err=y[1]