Lots more information in the documentation that will likely prove useful to call C from Julia.
Can call C functions directly from shared C libraries using ccall
Shared C libraries called by (:function, "library")
:function is the exported C function name
"library" is the name of the shared libraryIf C library is available as a shared library only needs to be on system LOAD_PATH, Otherwise need to compile it and put the library on the LOAD_PATH. Will see in example later.
ccallccall is like an ordinary function call, can use in REPLccall( (:function, “library”),
Return Type (i.e. Int32),
Tuple of Input Types. (i.e. (Int32, Int64, ...) ),
Arguments to pass to C from Julia )Note: ccall automatically converts Julia type to C type. See Julia documentation for a list of type correspondences.
In practice, can be useful to wrap ccall in a Julia function to:
A Julia package that wraps for John Burharkdt's C code to generate nodes and weights for sparse grid integration.
Approximate a $D$-variate integral of the form
$$\int_{\Omega} g(x)f(x)dx$$where $\Omega\in D$ and $f(x)$ is a weighting function. A quadrature rule prescribes a $D\times R$ matrix of nodes $X$ and a $R$-vector of weights $w$ depending on $\Omega$ and $f$. The approximation is calculated as
$$g(X)w$$The built-in families of integration rules are identified by a 3-letter key:
As a general rule, the versions with the nested univariate rules are more efficient in the sense that they require fewer nodes to achieve the same level of polynomial exactness.
The C function I want to use, nwspgr, calls nested functions. See lines 3377-3428 in 'sparse_grids_hw.c'.
void nwspgr ( void rule ( int n, double x[], double w[] ),
int rule_order ( int l ),
int dim, int k, int r_size, int *s_size, double nodes[], double weights[] )Problems: Nested C functions are not compatible with ccall (as far as i know!)
$\rightarrow$ Added 3 functions in C so can use Julia's ccall to run nwspgr See lines 6061-6492 in 'sparse_grids_hw.c' and the header file: 'sparse_grids_hw.h'.
- int nwspgr_rule_size ( void rule ( int n, double x[], double w[] ),
int rule_order ( int l ), int dim, int k, int r_size);
- int nwspgr_rule_size_wrapper (int rule , int d, int k);
- void nwspgr_wrapper (int rule , int d, int k, double x[], double w[]);Two steps:
To use the code in Julia compile the C-code to create a dynamic library ready to be called by Julia. If using the 'sparse_grids_hw_jl.sh', it creates a dynamic library called 'sparselib.dylib' in a new folder, /lib.
To enable Julia to find the newly created library, add the path where the library is stored to Julia's DL_LOAD_PATH. One way to do this is to add the line
push!(Sys.DL_LOAD_PATH,"[Path to new library's folder]")to ~/.juliarc.jl.
To output nodes and weights specify which integration rule to use, the number of dimensions, and the level accuracy of the integration rule. The constructor function is:
nw = nwspgr( rule, dims, level )where rule is one of {"GQU","KPU","GQN","KPN"} and nw is an instance of
type nwspgr
rule :: ASCIIString
dims :: Int64
level :: Int64
nodes :: Array{Float64,2}
weights :: Vector{Float64}
endNote: The maximum dimension supported by the code is 20.
Julia code to where the Julia nwspgr type is defined is contained in the file sparse_grids_hw.jl.
Recall that nwspgr in C-code was nested $\rightarrow$ added separate C functions to make compatible with Julia ccall syntax.
So, in the constructor function in Julia there are ccall's to these `ccall' compatible C functions.
ccall for nwspgr_rule_size_wrapper to tell Julia the correct size of the Array needed to store the nodes and weights.
ccall for nwspgr_wrapper is passed a pointer to correctly sized arrays created in Julia and fills them in with nodes and weights.
Note: ccall on passes a reference location in memory of array, not the array itself.
Managing memory in non-standard ccall of C functions (i.e. nested functions)
General rule stated in Julia docs: The program that created the memory needs to clean it up once it is finished with it.
How to avoid problems? If passing Julia created array for C-code to fill in, make sure Julia array has correct size when array is initiated. Helps to avoid segmentation faults...
Note: Segmentation faults occur due to bad memory management. For example, can occur when a program tries to grab memory it thinks is free only to find there is something already there (or the reverse).
Integral to approximate:
$$\mathcal{I}=\int_{0}^{1}\cdots\int_{0}^{1}\underbrace{{\displaystyle \prod_{d=1}^{D}\frac{1}{\sqrt{2\pi\sigma^{2}}}e^{-\frac{1}{2}\left(\frac{x_{d}-\mu}{\sigma}\right)^{2}}}}_{g(x)}dx_{D}\cdots dx_{1}$$where $x_{d}\sim N\left(\mu,\sigma\right)$.
Let $\mu=0$ and $\sigma=2$, then the "truth":
$$\mathcal{I} = \frac{1}{2} \operatorname{erf} \left(\frac{1}{ \sigma\sqrt{2}}\right)$$because $\Phi(x) = \frac{1}{2}\left(1+\operatorname{erf}\left(\frac{x-\mu}{ \sigma\sqrt{2}}\right)\right)$
In [ ]:
include("./Git/SparseGridsHW/test/nwspgr_demo.jl")
Integral to approximate:
$$\mathcal{I}=\int_{-\infty}^{\infty}\cdots\int_{-\infty}^{\infty}\mathcal{P}(x)f\left(x\right)dx_{D}\cdots dx_{1}$$where $\mathcal{P}(x)=4x_{1}^{2}-3x_{2}^{3}+x_{3}^{9}-2x_{4}^{4}$ and $x\sim N\left(\mu,\Sigma\right)$.
In [ ]:
include("./Git/SparseGridsHW/test/nwspgr_mvrn.jl")