Returning multiple results from a LibraryLink function

Question

What is the best way to return multiple inhomogeneous results from a LibraryLink function? Imagine that the result of a single computation is several tensors of different dimensions and several numbers. They should all be returned at the same time.

I see two ways:

1. Get a MathLink connection and use MathLink API functions to return the results. Is there any advantage at all to using LibraryLink (instead of pure MathLink) in this case? Did anyone compare LibraryLink/MathLink data transfer performance so I can more easily assess what the performance hit would be?

2. Have several LibraryLink functions: one will send the input, do the computation and store the result; there would be separate functions for retrieving each result; finally there would be a function to free up the library memory used to store those results.

   I'm not particularly keen on doing this. It seems more trouble and it forces me to manually manage memory from Mathematica.

Are there any better options?

---

Update for version 10

Version 10 brings changes both to MathLink (now WSTP) and to LibraryLink. MathLink now has the "IntraProcess" protocol for fast communication within a single process. LibraryLink now has managed library expressions which allows Mathematica to auto-free library side data structures (at the cost of more boilerplate code). We have now new data types and there are examples of generic looking type specifications in the documentation such as {_,_} (same linked page).

Does the answer change for version 10.0 or 10.2? Are there less cumbersome solutions available now? Imagine having to write an interface to a library that will commonly return multiple results from the same function.

Answer 1

I think number 2 is the way to go. This was one of the reasons why LibraryLink was written. The TetGenLink interface is an example where multiple LibrayLink functions query a single tetgen c++ instance and this works very efficiently.

Unfortunately, there is no LibrayFree function that you could connect to.

Answer 2

I see a third way, which is a combination of your two choices, but its usefulness depends on the problem you are solving. I used this once to get the best from both worlds. What I needed was a similar thing to ComponentMeasurements in 3D. Therefore, I needed a component labeling algorithm (similar to MorphologicalComponents) which takes a volume, binarizes it and assigns each separate component a unique ID. For this I needed only read access to the volume and therefore, I wanted to use it inside the library without copying it.

For this part I wrote a function using your second point which stored the result, the labeled volume, statically inside the library instance.

For each component in the labeled volume I wanted to have certain measures. My result was therefore a list which consisted of a set of several measurements for each component. This was something like

{"ID"->1, "Volume"->123, "BoundingBox"->{{..},{..},{..}}}

for each component, so highly inhomogeneous. To retrieve this result I wrote another function using MathLink where it didn't matter that I couldn't use the memory advantages of LibraryLink because I didn't have to transfer the large volume data.

I personally like how you can build very complex expressions with MathLink and for me it would have been very inconvenient if I had to write a LibraryLink function to retrieve every single property of my components.

void MLPutComponentProperties(MLINK mlp, const Component3d<mint> &c) {
    MLPutFunction(mlp, "List", c.number_of_properties);
    MLPutFunction(mlp, "Rule", 2); MLPutString(mlp, "Id"); MLPutInteger(mlp, c.label);
    MLPutFunction(mlp, "Rule", 2); MLPutString(mlp, "Volume"); MLPutInteger(mlp, c.volume);
    MLPutFunction(mlp, "Rule", 2); MLPutString(mlp, "BoundingBox");
        MLPutFunction(mlp, "List", 3);
        MLPutFunction(mlp, "List", 2);
          MLPutInteger(mlp, c.bounding_box.x1+1);
          MLPutInteger(mlp, c.bounding_box.x2+1);
        MLPutFunction(mlp, "List", 2);
          MLPutInteger(mlp, c.bounding_box.y1+1);
          MLPutInteger(mlp, c.bounding_box.y2+1);
        MLPutFunction(mlp, "List", 2);
          MLPutInteger(mlp, c.bounding_box.z1+1);
          MLPutInteger(mlp, c.bounding_box.z2+1);
}

void MLPutComponentList(MLINK mlp, const Component3dList components) {
    MLPutFunction(mlp, "List", components->size());
    Component3dListIterator it;
    for (it = components->begin(); it != components->end(); ++it) {
        MLPutComponentProperties(mlp, **it);
    }
}

Answer 3

Summary: I vouch for option 2 as LibraryLink is great with memory (sharing as opposed to copying). Furthermore I did a little investigation into the possible types of tensors with a negative/conservative result. The main point of posting the analysis to prevent that other people do the same.

Memory

Addressing: "Did anyone compare LibraryLink/MathLink data transfer performance so I can more easily assess what the performance hit would be"?

I think sending data over a MathLink connection always causes a copy of the data. In LibraryLink this is not necessary.

Firstly, you can load a function with the instruction that a tensor argument has to be shared between Mathematica and LibraryLink. This prevents copying information when giving the information from the kernel to the C program. Secondly, I think any tensor "returned" from a LibraryLink function using MArgument_setMTensor most likely does not have to be copied either. The main evidence for this is that the docs say

"Arguments passed to and from a library function can share data, saving on memory consumption and the time to copy large amounts of data."

in LibraryLink/tutorial/Introduction#163333181. This doc page also lists other differences between LibraryLink and MathLink. In particular the overhead of calling a LibraryLink function is much lower.

Example

Please see my answer here for an example of returning multiple outputs. It involves calling a separate function to access a global variable. I hope to learn more about this so that I can improve this code. At the moment it crashes for large inputs.

Analysis

From the LibraryLink user guide we have: "You can exchange not only C-like data types such as integers, reals, packed arrays, and strings, but also arbitrary Mathematica expressions".

Of course, from LibraryLink we can call MathLink, so I guess in this sense that statement is true (I am now quite confident that MathLink is necessary for this). However, I was hoping we could maybe make ragged arrays.

I investigated a bit, and found in WolframLibrary.h

#define MType_Integer 2
#define MType_Real 3
#define MType_Complex 4

My hope was now that there would be another type corresponding to 1. The following code can be used to investigate this further

DLLEXPORT int raggedArray_T_I(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res) {

    MTensor T_arg = MArgument_getMTensor(Args[0]);
    int type = libData->MTensor_getType(T_arg);
    MArgument_setInteger(Res, type);
    return 0;
}

It doesn't even really matter what the code above is. The point was to see if we could call a library function with something other than the usual "packable" arrays as arguments. I loaded the function above as follows

rAFu = LibraryFunctionLoad[libraryName, "raggedArray_T_I", {{_, _}}, 
  Integer]

Of course for the usual arrays the output was expected. One interesting thing is that any combinations of empty lists apparently has type 3 as a tensor, corresponding to an array of reals.

Unfortunately, I was unable to find any other input that qualified as a tensor. For example raFu[{{1},1}] gives an error that {{1},1} is not a tensor (of rank -1 of Removed[$$Failure] oddly).

Far fetched

My suspicion is now that maybe the valid types are 2,3 and 4 instead of 1, 2 and 3 is because of the following reason. Possibly expressions contain an indicator (maybe in the sense that they are structs have the indicator as an attribute, I read somewhere that expressions "are" structs (and unions)) containing the information what kind of expression the expression is. Then 1 could stand for "normal expression" and 2, 3 and 4 for the PackedArrays. Other types of expressions could be Integer, Real or Complex. That is what I did when I programmed my own little prototype of Mathematica anyway. I would be glad if somebody could shine more light on this.

Answer 4

Motivation

I was inspired by the discussion in this thread to investigate using JSON to return multiple data types in a single LibraryLink call. I'm posting this here so others might have an example of how easy it is to do this.

By using a header-only JSON library, we can quickly create a list of json data types (integers, reals, strings, Objects, or arrays thereof), and serialize it to a string which is returned via as a "UTF8String" to Mathematica. This has the added benefit of being able to use a single exported LibraryLink function to return different data types depending on the input.

The general plan on the library side is

• Create a new json object
• Add the data to the object,
  • You can do this either by adding key/value pairs, or using push_back
  • Alternatively you can instantiate the json object directly using the assignment operator.
• Serialize this json object to a string, and return this through a LibraryLink call
• In the Wolfram Language, your LibraryFunction will return a string. Now feed this to either Developer`ReadRawJSONString or Developer`ReadExpressionJSONString

The overhead involved in creating a json object, serializing, and then deserializing is smaller than you'd think. The Developer`functions are blazingly fast. If you are passing large amounts of numeric data, I'd imagine MathLink would be faster. But timing tests by Szabolcs in the linked community discussion seem to say the opposite.

The reason I like using json is simply the speed of writing the code. So much better than MathLink. I use this in conjunction with the LTemplate package, where I've made some modifications that apply the Developer` automatically. But the example below is self-contained.

Example

The contents of "test_json.cpp" are pasted at the end. The library is compiled via

<<CCompilerDriver`
lib = CreateLibrary[
    {"~/test_json/test_json.cpp"},
    "test_json_lib","TargetDirectory"->"~/test_json/",
    "CompileOptions"->{"-std=c++11"}
]

There are three exported functions,

testRawJSON = LibraryFunctionLoad[lib, "returnRawJSON", {}, "UTF8String"];
testGetExpression =  LibraryFunctionLoad[lib, 
      "returnExprs", {}, "UTF8String"];
testGetKeyValue = LibraryFunctionLoad[lib, 
      "getKeyValue", {"UTF8String", "UTF8String"}, "UTF8String"];

The first example returns a simple json list containing an Integer, a boolean, and an Association (what do we call these in json? Objects?)

testRawJSON[]
(* "[1,true,{\"BoundingBox\":[[0.0,0.0,0.0],[1.0,1.0,1.0]],\"ID\":1,\"Volume\":123}]" *)

Notice that the Association returned has the keys sorted automatically on the library side, in contrast to an Association.

Developer`ReadRawJSONString@%
(* {1, True, <|"BoundingBox" -> {{0., 0., 0.}, {1., 1., 1.}}, 
  "ID" -> 1, "Volume" -> 123|>} *)

Now let's use the library to return an arbitrary WL expression,

testGetExpression[]
(* "[\"List\",[\"List\",\"'Hello'\",123.5,1,[\"Style\",\"'Test'\
\",\"Bold\"]],[\"List\",1,2,3],[\"List\",[\"DateObject\",\"Now\"]]]" *)

Convert it to an expression, and notice how the DateObject evaluates,

Developer`ReadExpressionJSONString@%

Here is an example of passing in a json object to the library, so that you can access the values by their keys.

str = 
  Developer`WriteRawJSONString@<|"BoundingBox" -> {{0.`, 0.`, 
       0.`}, {1.`, 1.`, 1.`}}, "ID" -> 1, "Volume" -> 123|>;
testGetKeyValue[str, "BoundingBox"] // Developer`ReadRawJSONString
testGetKeyValue[str, "Volume"] // Developer`ReadRawJSONString
(* {{0., 0., 0.}, {1., 1., 1.}} *)
(* 123 *)

Here is the source code for "test_json.cpp", you need to make sure you have the json library by Niels Lohmann, available on github.

#include "WolframLibrary.h"


// downloaded from
// https://github.com/nlohmann/json/blob/develop/src/json.hpp
#include "json.hpp"


using json = nlohmann::json;


char *outString = NULL;


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


    // creating a list of rules this way is
    // so much faster for development than using
    // mathlink, code can be written nearly as fast
    // as in WL.  Here is how we can create an
    // Association
    json assoc;
    assoc["ID"] = 1;
    assoc["Volume"] = 123;
    assoc["BoundingBox"] = {
        {0.0, 0.0, 0.0},
        {1.0, 1.0, 1.0}
    };

    // you don't need to initialize the data when declaring the
    // object, you can grow it piece by piece
    json list;
    list.push_back(1);
    list.push_back(true);
    list.push_back(assoc);


    auto rawJSONString = list.dump();


    // some boilerplate to return a string
    if (outString!=NULL) {
      delete outString;
    }
    outString = new char[rawJSONString.length() + 1];
    std::strcpy(outString , rawJSONString.c_str());
    MArgument_setUTF8String(Res, outString);
        return LIBRARY_NO_ERROR;
}

// In an ExpresssionJSON string, symbols are represented as strings,
// therefore Strings need to be "escaped" with the ' character
std::string stringify(std::string input) {
    return "'" + input + "'";
}

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

    json styleObject = {
        "Style",
        stringify("Test"),
        "Bold"
    };
    json arg1 = {
        "List",
        stringify("Hello"),
        123.5,
        1,
        styleObject
    };
    json arg2 = {"List",1, 2, 3};
    json dateObject = {
        "DateObject",
        "Now"
    };
    json arg3 = {"List", dateObject};

    // here's an example of growing a JSON list item by item,
    // this is a major advantage over mathlink, where you need
    // to declare the number of arguments when you declare the head
    json res;
    res.push_back("List");
    res.push_back(arg1);
    res.push_back(arg2);
    res.push_back(arg3);


    auto exprJSONString = res.dump();
    if (outString!=NULL) {
      delete outString;
    }

    outString = new char[exprJSONString.length() + 1];
    std::strcpy(outString , exprJSONString.c_str());
    MArgument_setUTF8String(Res, outString);
    return LIBRARY_NO_ERROR;
}

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

    const char* assocRaw;
    const char* key;

    assocRaw = MArgument_getUTF8String(Args[0]);
    key = MArgument_getUTF8String(Args[1]);

    // here I'm assuming the association that was passed in
    // is a valid json object - with string-valued keys and
    // values that are themselves valid json objects.  If this
    // weren't the case, then Developer`WriteRawJSONString would
    // have issued a message.
    json assoc = json::parse((std::string) assocRaw);

    // now access the value associated with the second key
    // given as the second argument.
    json res = assoc[ key];

        auto rawJSONString = res.dump();
    if (outString!=NULL) {
      delete outString;
    }

    outString = new char[rawJSONString.length() + 1];
    std::strcpy(outString , rawJSONString.c_str());
    MArgument_setUTF8String(Res, outString);

    libData->UTF8String_disown(const_cast<char *>(assocRaw));
    libData->UTF8String_disown(const_cast<char *>(key));
    return LIBRARY_NO_ERROR;
}

Please feel free to suggest improvements on the c-code, I am new to it.

Answer 5

Another potential solution is using "Shared" passing, and writing a function that modifies its arguments. In other words: use pass-by-reference.

I am not quite certain that this is among the intended uses of "Shared" argument passing, so I advise caution. I have not yet seriously used this method. Pass-by-reference isn't really compatible with how Mathematica works. To use this method successfully, one must have some familiarity of how Mathematica manages its expressions internally (copy-on-write paradigm).

---

Here's a demonstration using LTemplate. (LTemplate is not necessary at all for this. It just made it quicker for me to experiment.)

SetDirectory[$TemporaryDirectory];

<< LTemplate`

tem = LClass["Ret",
   {LFun["incdec", {{Integer, _, "Shared"}, {Integer, _, "Shared"}}, 
     "Void"]}
   ];

code = "
  struct Ret {
    void incdec(mma::IntTensorRef a, mma::IntTensorRef b) {
        for (mint i=0; i < a.size(); ++i)
            a[i] = i+1;
        for (mint i=0; i < b.size(); ++i)
            b[i] = b.size() - i;
        a.disown();
        b.disown();
    }
  };
  ";
Export["Ret.h", code, "String"];

CompileTemplate[tem]

LoadTemplate[tem]

obj = Make[Ret];

This function fills its first and second integer vector argument with increasing and decreasing numbers. It can be used like this:

a = ConstantArray[0, 10];
b = ConstantArray[0, 5];

{a, b}
(* {{0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0}} *)

obj@"incdec"[a, b]

{a, b}
(* {{1, 2, 3, 4, 5, 6, 7, 8, 9, 10}, {5, 4, 3, 2, 1}} *)

Notice that both a and b were modified.

One should be very careful when doing such things because they violate the basic principles of how Mathematica is expected to work. You must be aware, for example, that Mathematica uses copy-on-write. Thus when setting c=b, the contents of b are not internally copied, only referenced. The actual copying only happens when modifying a sub-part of c. However, modifications done by a LibraryLink library cannot be detected by the Mathematica kernel, so we may end up in a situation when more than one value appears to get modified.

a = ConstantArray[0, 10];
b = ConstantArray[0, 5];
c = b;

{a, b, c}
(* {{0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0}, {0, 0, 0, 0, 0}} *)

obj@"incdec"[a, b]

{a, b, c}
(* {{1, 2, 3, 4, 5, 6, 7, 8, 9, 10}, {5, 4, 3, 2, 1}, {5, 4, 3, 2, 1}} *)

Notice that we only modified b, but c changed as well.

We must carefully ensure that the packed-array arguments we pass to this LibraryLink function are referenced by only a single Mathematica value. We must create these packed arrays ourselves, and never use a value that comes from a source we don't directly control.

We can create a deep copy of a packed array like this:

a = ConstantArray[0, 10];
b = ConstantArray[0, 5];
c = b[[ ;; ]];

{a, b, c}
(* {{0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0}, {0, 0, 0, 0, 0}} *)

obj@"incdec"[a, b]

{a, b, c}
(* {{1, 2, 3, 4, 5, 6, 7, 8, 9, 10}, {5, 4, 3, 2, 1}, {0, 0, 0, 0, 0}} *)

There are other things to pay attention to. LibraryLink functions can only work with packed arrays. If we pass in a non-packed one, it will be packed, and therefore copied. The library function will modify the copy, not the original array.

a = ConstantArray[0, 10];
b = Developer`FromPackedArray@ConstantArray[0, 5];

{a, b}
(* {{0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0}} *)

obj@"incdec"[a, b]

LibraryFunction::shrnopack: The argument {0,0,0,0,0} at position 3 was not a PackedArray but was successfully converted to a PackedArray.

{a, b}
(* {{1, 2, 3, 4, 5, 6, 7, 8, 9, 10}, {0, 0, 0, 0, 0}} *)

Notice that b was not modified.

Yet another thing that can go wrong is that the input array is packed, but not of the type expected by the library. In this case, it will be silently converted to the correct type (and thus copied again).

a = ConstantArray[0, 10];
b = N@ConstantArray[0, 5];

{a, b}
(* {{0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0., 0., 0., 0., 0.}} *)

obj@"incdec"[a, b]

{a, b}
(* {{1, 2, 3, 4, 5, 6, 7, 8, 9, 10}, {0., 0., 0., 0., 0.}} *)

Again, b was not modified.

As you can see, there are plenty of things that can go wrong ... Be very cautious if you decide to use this method to return multiple results.

Answer 6

What is the best way to return multiple inhomogeneous results from a LibraryLink function?

There are two ways to pass arguments and return result from LibraryLink:

• Using the standard types such as MTensor, which is very fast. This allows for returning one result only per function call.

• Using LinkObject-based passing. This is done through MathLink/WSTP (a loopback link) and allows returning any number of results, but it is slower.

Thus there are two practical solutions to returning multiple results.

1. Use separate function calls for retrieving each result. The sequence of calls will like this:

• Compute all results and store them
• Retrieve 1st result
• Retrieve 2nd result
• ...
• Free up space used for storing results.

Each of these will be a separate call.

The LTemplate package makes it a bit easier to do this by automating freeing up results. Check the "MeanVariance" class example.

2. Use LinkObject-based passing. This is often slower than several plain LibraryLink function calls.

Did anyone compare LibraryLink/MathLink data transfer performance so I can more easily assess what the performance hit would be?

I did a partial comparison:

• MathLink (WSTP) vs LibraryLink performance

This can be used a guideline for deciding whether a MathLink or pure-LibraryLink based solution may be faster for a given use case. The best solution will depend on the number of results that need to be returned, as well as on the amount of data that must be transferred.