Skip to main content

development - How to simplify writing LibraryLink code?


LibraryLink is an API for extending Mathematica through C or C++. It is very fast because it gives direct access to Mathematica's packed array data structure, without even needing to make a copy of the arrays.


Unfortunately, working with LibraryLink involves writing a lot of tedious boilerplate code.


Mathematica 10 introduced managed library expressions, a way to have C-side data structures automatically destroyed when Mathematica no longer keeps a reference to them. I find this feature useful for almost all non-trivial LibraryLink projects I make, but unfortunately it requires even more boilerplate.


For me, writing all this repetitive code (and looking up how to do it in the documentation) has been the single biggest barrier to starting any LibraryLink project. It just takes too much time. Is there an easier way?



Answer



I wrote a package to automatically generate all the boilerplate needed for LibraryLink and for managed library expressions based on a template that describes a class interface:



Here's how it works




  • Write a template (a Mathematica) expression that describes a C++ class interface

  • This template is used to generate LibraryLink-compatible C functions that call the class's member functions, to compile them, and to finally load them

  • Instances of the class are created as managed library expressions; all the boilerplate code for this is auto-generated

  • The package is designed to be embedded into Mathematica applications. This avoids conflicts between different applications using different LTemplate versions. LTemplate can also be used standalone for experimentation and initial development.


Examples


Note: The package comes with a tutorial and many more examples in its Documentation subdirectory.


While LTemplate is designed to be embedded into other applications, for this small example I'll demonstrate interactive standalone use.


<< LTemplate`


Let us use the system's temporary directory as our working directory for this small example, so we don't leave a lot of mess behind. Skip this to use the current directory instead.


SetDirectory[$TemporaryDirectory]

To be able to demonstrate managed library expressions later, we turn off input/output history tracking:


$HistoryLength = 0;

Let's write a small C++ class for calculating the mean and variance of an array:


code = "
#include


class MeanVariance {
double m, v;

public:
MeanVariance() { mma::print(\"constructor called\"); }
~MeanVariance() { mma::print(\"destructor called\"); }

void compute(mma::RealTensorRef vec) {
double sum = 0.0, sum2 = 0.0;
for (mint i=0; i < vec.length(); ++i) {

sum += vec[i];
sum2 += vec[i]*vec[i];
}
m = sum / vec.length();
v = sum2 / vec.length() - m*m;
}

double mean() { return m; }
double variance() { return v; }
};

";

Export["MeanVariance.h", code, "String"]

LTemplate functions live in the mma C++ namespace. mma::RealTensorRef represents a reference to a Mathematica packed array of Reals. It is just a wrapper for the MTensor type. There may be more than one reference to the same array, so array creation and destruction must still be managed manually.


RealTensorRef (an alias for TensorRef) has a number of convenience features:




  • Linear indexing into the array using the [ ] operator. RealMatrixRef and RealCubeRef are subclasses of RealTensorRef, which additionally allow indexing into 2D or 3D arrays using the ( ) operator.





  • Iteration through elements using begin() and end(). Range-based for loops work, e.g. we could have used for (auto &x : vec) sum += x;.




We have three member functions: one for computing both the mean and variance in one go, and two others for retrieving each result.


The class must go in a header file with the same name, MeanVariance.h in this case.


Let's write the corresponding template now:


template =
LClass["MeanVariance",
{

LFun["compute", {{Real, _, "Constant"}}, "Void"],
LFun["mean", {}, Real],
LFun["variance", {}, Real]
}
];

We can optionally check that the template has no errors using ValidTemplateQ[template].


Compiling and loading is now as simple as


In[]:= CompileTemplate[template]


Current directory is: /private/var/folders/31/l_62jfs110lf0dh7k5n_y2th0000gn/T

Unloading library MeanVariance ...

Generating library code ...

Compiling library code ...

Out[]= "/Users/szhorvat/Library/Mathematica/SystemFiles/LibraryResources/MacOSX-x86-64/MeanVariance.dylib"


and then


LoadTemplate[template]

The compilation step created a file called LTemplate-MeanVariance.cpp in the same directory, the source code of which I attached at the end, for illustration.


We can now create a managed library expression corresponding to this class:


obj = Make["MeanVariance"]

During evaluation of In[]: constructor called

(* MeanVariance[1] *)


arr = RandomReal[1, 10];

obj@"compute"[arr]

{obj@"mean"[], obj@"variance"[]}
(* {0.482564, 0.104029} *)

We can check the list of live expressions of type MeanVariance using


LExpressionList["MeanVariance"]

(* {MeanVariance[1]} *)

As soon as Mathematica has no more references to this expression, it gets automatically destroyed


obj =.

During evaluation of In[]: destructor called

LExpressionList["MeanVariance"]
(* {} *)


The reason why we had to use $HistoryLength = 0 above is to prevent Out keeping a reference to the expression.


One practical way to write a Mathematica functions that exposes this functionality is


meanVariance[arr_] :=
Block[{obj = Make[MeanVariance]},
obj@"compute"[arr];
{obj@"mean"[], obj@"variance"[]}
]

As soon as the Block finishes, obj gets automatically destroyed:


meanVariance[arr]


During evaluation of In[]: constructor called

During evaluation of In[]: destructor called

(* {0.618649, 0.033828} *)

This is one of those special cases when using Block over Module may be worth it for performance reasons. (The usual caveats about Block apply though.)


Notice that the expression has the form MeanVariance[1]. The integer index 1 is the ManagedLibraryExpressionID. The symbol MeanVariance is created in the context


LClassContext[]

(* "LTemplate`Classes`" *)

This context is added to the $ContextPath when using LTemplate interactively as a standalone package, but not when it's loaded privately by another application. We can check the usage message of the symbol:


?MeanVariance

class MeanVariance:
Void compute(Constant List)
Real mean()
Real variance()




The package is too large to present fully in a StackExchange post, so if you are interested, download it and read the tutorial, and take a look at the many examples that come with the package!


LTemplate is continually under development, and breaking changes are possible (though unlikely). However, since the recommended way to deploy it is to embed it fully in the Mathematica application that uses it, this should not be a problem. I am using LTemplate in the IGraph/M package, which proves its feasibility for use in large projects.


There are additional features such as:



  • Multiple related classes in the same template

  • Support for additional data types, such as SparseArray, RawArray, Image/Image3D

  • Pass another managed library expression to a function, and receive it as object reference on the C++ side (LExpressionID)

  • Format templates to be human-readable (FormatTemplate)

  • User-friendly error messages


  • Error handling through exceptions (mma::LibraryError); unknown exceptions are also caught to prevent killing the kernel when possible.

  • Calling Print for debugging (also through a C++ streams interface), massert macro to replace C's standard assert and avoid killing the Mathematica kernel.

  • Calling Message, setting a symbol to associate standard messages with

  • Argument passing and return using MathLink (LinkObject passing)

  • mlstream.h auxiliary header for easier LinkObject-based passing


The documentation isn't complete, but if you have questions, feel free to comment here or email me.




Questions and limitations


Can I use plain C instead of C++? No, LTemplate requires the use of C++. However, the only C++ feature the end-user programmer must use is creating a basic class.



Why do I have to create a class? I only want a few functions. I didn't implement free functions due to lack of time and need. There's no reason why this shouldn't be added. However, you can always create a single instance of a class, and keep calling functions on it. The overhead will be minimal according to my measurements.


Why can't I use underscores in class or function names? LTemplate currently only supports names that are valid both in Mathematica and C++. This excludes underscores and \$ signs (even though some C++ compilers support $ in identifiers). This also helps avoid name conflicts with auxiliary functions LTemplate generates (which always have underscores).


How do I write library initialization and cleanup code? Currently LTemplate doesn't support injecting code into WolframLibrary_initialize and WolframLibrary_uninitialize. Initialization code can be called manually from Mathematica. Create a single instance of a special class, put the initialization code in one of its member functions, and call it from Mathematica right after loading the template. The uninitialization code can go in the destructor of the class. All objects are destroyed when the library is unloaded (e.g. when Mathematica quits). Warning: when using this method, there's no guarantee about which expression will be destroyed last! To fix this, initialization/uninitialization support is planned for later.


Can I create a function that takes an MTensor of unspecified type? No, LTemplate requires specifying the data type (but not the rank) of the MTensor. {Real, 2} is a valid type specifier and so is {Real, _}. {_, _} is not allowed in LTemplate, even though it's valid in standard LibraryLink. The same applies to MSparseArrays (mma::SparseArrayRef). However, RawArray and Image may be passed without an explicit element-type specification.


Which LibraryLink features are not supported?




  • The numerical type underlying tensors must be explicitly specified. Tensors without explicitly specified types are not supported. In the future tensors of unspecified types may be handled through C++ templates.





  • There's no explicit support for library callback functions yet, but they can be used by accessing the standard LibraryLink API (function pointers in mma::libData).




Feedback


Contributions and ideas for improvements are most welcome! Feel free to email me. If you find a bug, file a bug report.




Source of LTemplate-MeanVariance.cpp:


#define LTEMPLATE_MMA_VERSION  1120

#include "LTemplate.h"

#include "LTemplateHelpers.h"
#include "MeanVariance.h"


#define LTEMPLATE_MESSAGE_SYMBOL "LTemplate`LTemplate"

#include "LTemplate.inc"


std::map MeanVariance_collection;


DLLEXPORT void MeanVariance_manager_fun(WolframLibraryData libData, mbool mode, mint id)
{
if (mode == 0) { // create
MeanVariance_collection[id] = new MeanVariance();
} else { // destroy
if (MeanVariance_collection.find(id) == MeanVariance_collection.end()) {
libData->Message("noinst");
return;
}

delete MeanVariance_collection[id];
MeanVariance_collection.erase(id);
}
}

extern "C" DLLEXPORT int MeanVariance_get_collection(WolframLibraryData libData, mint Argc, MArgument * Args, MArgument Res)
{
mma::TensorRef res = mma::detail::get_collection(MeanVariance_collection);
mma::detail::setTensor(Res, res);
return LIBRARY_NO_ERROR;

}


extern "C" DLLEXPORT mint WolframLibrary_getVersion()
{
return 3;
}

extern "C" DLLEXPORT int WolframLibrary_initialize(WolframLibraryData libData)
{

mma::libData = libData;
{
int err;
err = (*libData->registerLibraryExpressionManager)("MeanVariance", MeanVariance_manager_fun);
if (err != LIBRARY_NO_ERROR) return err;
}
return LIBRARY_NO_ERROR;
}

extern "C" DLLEXPORT void WolframLibrary_uninitialize(WolframLibraryData libData)

{
(*libData->unregisterLibraryExpressionManager)("MeanVariance");
return;
}


extern "C" DLLEXPORT int MeanVariance_compute(WolframLibraryData libData, mint Argc, MArgument * Args, MArgument Res)
{
mma::detail::MOutFlushGuard flushguard;
const mint id = MArgument_getInteger(Args[0]);

if (MeanVariance_collection.find(id) == MeanVariance_collection.end()) { libData->Message("noinst"); return LIBRARY_FUNCTION_ERROR; }

try
{
mma::TensorRef var1 = mma::detail::getTensor(Args[1]);

(MeanVariance_collection[id])->compute(var1);
}
catch (const mma::LibraryError & libErr)
{

libErr.report();
return libErr.error_code();
}
catch (const std::exception & exc)
{
mma::detail::handleUnknownException(exc.what(), "MeanVariance::compute()");
return LIBRARY_FUNCTION_ERROR;
}
catch (...)
{

mma::detail::handleUnknownException(NULL, "MeanVariance::compute()");
return LIBRARY_FUNCTION_ERROR;
}

return LIBRARY_NO_ERROR;
}


extern "C" DLLEXPORT int MeanVariance_mean(WolframLibraryData libData, mint Argc, MArgument * Args, MArgument Res)
{

mma::detail::MOutFlushGuard flushguard;
const mint id = MArgument_getInteger(Args[0]);
if (MeanVariance_collection.find(id) == MeanVariance_collection.end()) { libData->Message("noinst"); return LIBRARY_FUNCTION_ERROR; }

try
{
double res = (MeanVariance_collection[id])->mean();
MArgument_setReal(Res, res);
}
catch (const mma::LibraryError & libErr)

{
libErr.report();
return libErr.error_code();
}
catch (const std::exception & exc)
{
mma::detail::handleUnknownException(exc.what(), "MeanVariance::mean()");
return LIBRARY_FUNCTION_ERROR;
}
catch (...)

{
mma::detail::handleUnknownException(NULL, "MeanVariance::mean()");
return LIBRARY_FUNCTION_ERROR;
}

return LIBRARY_NO_ERROR;
}


extern "C" DLLEXPORT int MeanVariance_variance(WolframLibraryData libData, mint Argc, MArgument * Args, MArgument Res)

{
mma::detail::MOutFlushGuard flushguard;
const mint id = MArgument_getInteger(Args[0]);
if (MeanVariance_collection.find(id) == MeanVariance_collection.end()) { libData->Message("noinst"); return LIBRARY_FUNCTION_ERROR; }

try
{
double res = (MeanVariance_collection[id])->variance();
MArgument_setReal(Res, res);
}

catch (const mma::LibraryError & libErr)
{
libErr.report();
return libErr.error_code();
}
catch (const std::exception & exc)
{
mma::detail::handleUnknownException(exc.what(), "MeanVariance::variance()");
return LIBRARY_FUNCTION_ERROR;
}

catch (...)
{
mma::detail::handleUnknownException(NULL, "MeanVariance::variance()");
return LIBRARY_FUNCTION_ERROR;
}

return LIBRARY_NO_ERROR;
}

Comments

Popular posts from this blog

plotting - Filling between two spheres in SphericalPlot3D

Manipulate[ SphericalPlot3D[{1, 2 - n}, {θ, 0, Pi}, {ϕ, 0, 1.5 Pi}, Mesh -> None, PlotPoints -> 15, PlotRange -> {-2.2, 2.2}], {n, 0, 1}] I cant' seem to be able to make a filling between two spheres. I've already tried the obvious Filling -> {1 -> {2}} but Mathematica doesn't seem to like that option. Is there any easy way around this or ... Answer There is no built-in filling in SphericalPlot3D . One option is to use ParametricPlot3D to draw the surfaces between the two shells: Manipulate[ Show[SphericalPlot3D[{1, 2 - n}, {θ, 0, Pi}, {ϕ, 0, 1.5 Pi}, PlotPoints -> 15, PlotRange -> {-2.2, 2.2}], ParametricPlot3D[{ r {Sin[t] Cos[1.5 Pi], Sin[t] Sin[1.5 Pi], Cos[t]}, r {Sin[t] Cos[0 Pi], Sin[t] Sin[0 Pi], Cos[t]}}, {r, 1, 2 - n}, {t, 0, Pi}, PlotStyle -> Yellow, Mesh -> {2, 15}]], {n, 0, 1}]

plotting - Plot 4D data with color as 4th dimension

I have a list of 4D data (x position, y position, amplitude, wavelength). I want to plot x, y, and amplitude on a 3D plot and have the color of the points correspond to the wavelength. I have seen many examples using functions to define color but my wavelength cannot be expressed by an analytic function. Is there a simple way to do this? Answer Here a another possible way to visualize 4D data: data = Flatten[Table[{x, y, x^2 + y^2, Sin[x - y]}, {x, -Pi, Pi,Pi/10}, {y,-Pi,Pi, Pi/10}], 1]; You can use the function Point along with VertexColors . Now the points are places using the first three elements and the color is determined by the fourth. In this case I used Hue, but you can use whatever you prefer. Graphics3D[ Point[data[[All, 1 ;; 3]], VertexColors -> Hue /@ data[[All, 4]]], Axes -> True, BoxRatios -> {1, 1, 1/GoldenRatio}]

plotting - Adding a thick curve to a regionplot

Suppose we have the following simple RegionPlot: f[x_] := 1 - x^2 g[x_] := 1 - 0.5 x^2 RegionPlot[{y < f[x], f[x] < y < g[x], y > g[x]}, {x, 0, 2}, {y, 0, 2}] Now I'm trying to change the curve defined by $y=g[x]$ into a thick black curve, while leaving all other boundaries in the plot unchanged. I've tried adding the region $y=g[x]$ and playing with the plotstyle, which didn't work, and I've tried BoundaryStyle, which changed all the boundaries in the plot. Now I'm kinda out of ideas... Any help would be appreciated! Answer With f[x_] := 1 - x^2 g[x_] := 1 - 0.5 x^2 You can use Epilog to add the thick line: RegionPlot[{y < f[x], f[x] < y < g[x], y > g[x]}, {x, 0, 2}, {y, 0, 2}, PlotPoints -> 50, Epilog -> (Plot[g[x], {x, 0, 2}, PlotStyle -> {Black, Thick}][[1]]), PlotStyle -> {Directive[Yellow, Opacity[0.4]], Directive[Pink, Opacity[0.4]],