Message Boards Message Boards

GROUPS:

Mathematica++ A C++ Library with Sugar

Posted 4 months ago
1232 Views
|
16 Replies
|
23 Total Likes
|

I was working on a library for my simulations. The difficult part was the shape of its input and output requires sophisticated UI. So after few days of manipulation with Manipulate I switched back to C++/Qt based UI and Mathematica for computation. So I was working on a C++ library to interact with Mathematica using WSTP API. It is usable but, not yet complete, however nothing is ever complete. I named it mathematica++ and released it under FreeBSD License. it is hosted in gitlab. gitlab project page.

I don't want to bloat this post, These are some of the examples copy pasted from the above links

symbol x("x");
value  res;
std::string method = "Newton";

shell << Values(FindRoot(ArcTan(1000 * Cos(x)), List(x, 1, 2),  Rule("Method") = method));
shell >> res;
std::vector<double> results = cast<std::vector<double>>(res);
std::cout << results[0] << std::endl; // Prints 10.9956

Example 2

mathematica::m mata = Table(Mod(i + j, 2), List(i, 1, 2), List(j, 1, 2));
mathematica::m matb = Table(Mod(i + j, 3), List(i, 1, 2), List(j, 1, 2)];
mathematica::m matc = Dot(mata, matb);
mathematica::m matd = Det(matc);

// Execute mathematica constructs and fetch the response
shell << matd;
shell >> determinant;

// determinant can be converted to C++ machine sized types
std::cout << determinant << std::endl; // Prints -2

cast<T> brings the output returned by mathematica to STL containers

value result;
typedef std::vector<std::vector<int>> ivv_type;
shell << FactorInteger(2434500);
shell >> result;
ivv_type prime_powers = mathematica::cast<ivv_type>(result);
for(auto pp: prime_powers){
    std::cout << pp[0] << " ^ " << pp[1] << std::endl;
}
// Prints the following output
// 2 ^ 2
// 3 ^ 2
// 5 ^ 3
// 541 ^ 1

Not limited to std::vector only

Thank You.

16 Replies

That's exactly what I was hoping/looking for, before I decided to also learn the Wolfram Language!! Absolutely amazing, thanks for your hard, fruit-bearing work which I know from gitlab. I've tested some other sample codes and so far so good, everything works like a charm, fantastic!!

Mathematica++

Simply brilliant!!!

Thank you for the very first review. I am glad to know that it worked nicely.

This looks quite interesting. Since I don't have the time to delve into it right now, would you mind answering a few questions?

Do I understand it right that this is syntactic sugar for putting expressions onto a MathLink/WSTP link?

Does it only work with wstp.h or also mathlink.h? I believe LibraryLink (which is usually the best way to extend Mma) still requires MathLink and won't work with the WSTP headers. I might be wrong though as I have not investigated this in 11.3.

Based on the example you posted, it would appear that you implemented an independent way to store Mma expressions. These expressions are built at runtime, then sent to the link. Is this correct? If yes, what is the overhead from building them at runtime?

I was always wondering if it was possible to use some C++ template magic to do the same at compile time, i.e. create convenient syntax that effectively translates directly to MathLink code. E.g. something like

link << expr("Plus", 42, expr("Times", 2, symbol("x"));

would directly translate to the equivalent of

MLPutFunction(link, "Plus", 2);
    MLPutInteger(link, 42);
    MLPutFunction(link, "Times", 2);
        MLPutInteger(link, 2);
        MLPutSymbol(link, "x");

at compile time without the need to build an intermediate representation at runtime.

It is not clear to me if it is technically possible to do something like this. Have you seen the mlstream header of my LTemplate tool? It would be a lovely extension. Currently, mlstream handles putting/getting basic types (including arrays), with a focus on performance. It does not handle generic expressions.

While we're at the topic of performance, does it handle large numerical array efficiently, e.g. MLPutReal64List or MLPutReal64Array?

Yes it uses template magic in the expected way, while overloading operators for symbols and composite constructs. Here the symbols are declared using

mathematica::symbol x("x");
mathematica::symbol x("y");

Expressions like x + y translates to Plus[x, y] Composite expressions like Sin[x] are represented using Sin(x) However this Sin function is previously declared using

MATHEMATICA_DECLARE(Sin)

Now this actually creates an object mathematica::m("Sin") that have operator() overloaded. and operator overloads works for symbols as well as m constructs. Thus things gets chained.

Actually I am storing chain of delayed callbacks that are fired in order. When they fire the corresponding WS function gets invoked.

In the code there is a debug macro that can be set to see to 1 to check what it is sending and receiving to/from mathematica. It is more or less similar to the example in your post.

It uses wstp.h only, and only WS* functions. However user code does not need to include that header. A good amount of WS functions and ML functions share their name so it is more or less the same (one macro away). It is trivial to keep a macro that switches between WS functions and ML functions.

I never knew about LTemplate. Checked just now. It looks like LTemplate makes it robust for mathematica module development. In mathematica++ I focuced on application <--> mathematica part. In Ltemplate you focused on mathematica <--> module part. Certainly an integration of ltemplate + mathematica++ will contribute to a complete library.

I have never checked performance. While implementing I avoided functions like PutReal64List Instead I have gone in a more generic pathway. creating a list first and then Put numbers one by one. Using such functions will be a case by case optimization that I left for latter.

Amy actual work was related to optimization. The output was several (about < 25) Export[] of images and 3d graphics. What I remember it used to take significant amount of time after all computations are over. I suspect that is mostly due to so many exports.

Yes it uses template magic in the expected way ...

... Actually I am storing chain of delayed callbacks that are fired in order.

What I was trying to say was that it would be nice to have zero runtime overhead---if at all possible. As I understand from your answer, this library will effectively translate the convenience syntax to WS... calls at runtime, not at compile time. There is extra code that runs when using this library, that would be avoided if using WS... functions directly.

My question was: do you have any idea of how great this overhead is?

While implementing I avoided functions like PutReal64List Instead I have gone in a more generic pathway. creating a list first and then Put numbers one by one.

This is however very inefficient when working with numerical arrays, a very common task in scientific programming ...

It uses wstp.h only, and only WS* functions. However user code does not need to include that header. A good amount of WS functions and ML functions share their name so it is more or less the same (one macro away). It is trivial to keep a macro that switches between WS functions and ML functions.

MathLink and WSTP are in fact exactly the same, not more or less the same. "Replacing" MathLink with WSTP was an entirely pointless marketing move that is causing real inconvenience and extra work for people who work with this technology.

The problem is that if one decides to write a program that runs in a separate process and communicates with the kernel through MathLink, one does have the option to choose either MathLink (ML prefix) or WSTP (WS prefix), or switch between the two with the preprocessor technique you mentioned. However, when using LibraryLink, which is usually a better option than installable MathLink programs, one must use MathLink (ML-prefix functions). Thus basing such a library on the wstp.h header will limit its usability (no LibraryLink). Basing it on mathlink.h wouldn't.

I should have spelt out that I am very interested in using this library if it can be made compatible with LibraryLink, and it there are workarounds for potential performance issues (e.g. use a PutReal64Array manually when advantageous, even as part of a larger expression).

With LibraryLink one is basically given an MLINK (WSLINK) object and then uses it to read off the arguments passed to the corresponding Mathematica function, then writes the return value to it.

The limitation is that LibraryLink still only works with mathlink.h. It is possible that it can be made to work with wstp.h, but this would not be straightforward and so far I haven't tried very hard (as I would not have gotten anything in return--it would only have made compilation more difficult).

Thanks for the suggestion. I am not sure about the first one but the second one is easy to accomplish. Although PutReal64Array is feasible I am doubtful about GetReal64Array

What I was trying to say was that it would be nice to have zero runtime overhead---if at all possible. As I understand from your answer, this library will effectively translate the convenience syntax to WS... calls at runtime, not at compile time. There is extra code that runs when using this library, that would be avoided if using WS... functions directly.

I am not sure about the feasibility of all compile time solution.I may work to some extent but there are scenarios that are difficult to tackle this way. Given a very simple use case

List[1, 2, 3, 4] // mathematica syntax

A similar C++ syntax List<1, 2, 3> which build the WS call chain on compile time, but it is not possible to pass string in this way.

A syntax like a = List(1, 2, 3, "hallo") can builds the call chain at compile time. But a scenario like the following

a = List(); // empty List initiated
for (auto x: l){
    a.arg(x);
}

expecting to build a List (not append) at the end of the loop requires a runtime storage for the call chain. Otherwise if it entirely built on templates some more magic like b = arg(a, x) can be used to abandon a and generate a new b. i.e. once an object is created it wont be mutable. Whether it is a good design or bad design is a matter of debate and depends on use case. I am not very much sure about such feasibility. And I am afraid that things may turn ugly.

I thought about this for quite a while and decided to store and populate it on runtime with delayed callbacks. Even if I generate the call chain on compile time there will be some wrapper around the WSTP functions. That part is unavoidable. The only extra overhead is the runtime storage of the call chain. Each delayed call ftor requires few words. So for the second part of your question, yes there is an extra code that is running but I don't consider it as a bottleneck.

This is however very inefficient when working with numerical arrays, a very common task in scientific programming

Yes then there needs to be case by case exception in the fetch call. However the problem is while fetching the library knows that it is fetching a List. But the library doesn't know that it is fetching a list of all 64 bit integers. While making Put calls the exceptional rules may still work but for Get calls I am not sure about its feasibility.

The problem is that if one decides to write a program that runs in a separate process and communicates with the kernel through MathLink, one does have the option to choose either MathLink (ML prefix) or WSTP (WS prefix), or switch between the two with the preprocessor technique you mentioned. However, when using LibraryLink, which is usually a better option than installable MathLink programs, one must use MathLink (ML-prefix functions). Thus basing such a library on the wstp.h header will limit its usability (no LibraryLink). Basing it on mathlink.h wouldn't.

I don't think there is any benefit criticizing wolfarm's decision in this regard this is what is happening for few years. It doesn't seem to go in the ML direction in future. So I targeted with WSTP only. However I kept all low level (WS) calls in few places so that it never turns hectic to change from WS to ML. Currently I am busy for an week. I'll see if I can push that change in few days.

While making Put calls the exceptional rules may still work ...

To me, it would seem that the library is most useful for the Put case.

but for Get calls I am not sure about its feasibility.

If you mean that you don't know how to detect whether the link has a packed array on it: this can be done, but only with not-well-documented functions. We had a long discussion about it the other day. See here.

But a more practical solution would be to design the library in such a way that users can (partially) fall back to direct use of the MathLink API when this is advantageous.

I can detect a packed array but I cannot detect whether all elements in that are are int64 or int32 or double or string because List is not uniform. So I wrap it in a variant type and fetch one by one.

I can detect a packed array but I cannot detect whether all elements in that are are int64 or int32 or double or string because List is not uniform.

Do you mean that you can detect a List expression (not a packed array)?

A packed array is always homogeneous. All elements will have the same numeric type (no strings allowed). See the link I posted (chatroom discussion) on how to detect not a generic List, but specifically a packed array, and how to determine the type of its elements.

Also take a look here: https://mathematica.stackexchange.com/questions/3496/what-is-a-mathematica-packed-array

A packed array is essentially what you might call a C array---a contiguous block of numerical data (e.g. doubles). It is the storage format that is the most efficient as well as the easiest to handle in C (or C++). In most scientific applications that I came across, it is the single most useful data structure that one might want to exchange with Mathematica. This is also why LibraryLink supports it directly (see MTensor).

Aha I see. Yes I was talking about detecting a List. I'll have a look into packed array

I was working on switching between ML and WS. Now -DWITH_ML=ON can be set on cmake to use MathLink instead of using WSTP. I just pushed the changes into develop branch. Things are working as expected.

I was working on switching between ML and WS.

That's great! I'll try it out on the weekend :-)

I just pushed some changes to use The PutArray functions is the input is a nesting of std::vector<std::vector<...>> and that vector forms a matrix or higher dimensional matrix or tensor of proper shape.

typedef std::vector<std::vector<std::vector<int>>> ivv_type; // 2 x 2 x 3
ivv_type matrix;
std::vector<int> row11 = {111, 112, 113};
std::vector<int> row12 = {121, 122, 123};
std::vector<int> row21 = {211, 212, 213};
std::vector<int> row22 = {221, 222, 223};

std::vector<std::vector<int>> row1;
std::vector<std::vector<int>> row2;

row1.push_back(row11);
row1.push_back(row12);

row2.push_back(row21);
row2.push_back(row22);

matrix.push_back(row1);
matrix.push_back(row2);

value res;

shell.enable(mathematica::bulk_io); // ENABLE this feature
shell << Total(Flatten(matrix));
shell >> res;

BOOST_CHECK(*res == 2004);

The above will use MLPutInteger32Array/ WSPutInteger32Array

These two changes are currently in develop branch. I'll have some more tests and then I'll put them in master branch.

After that I'll use Armadillo Eigen etc.. as optional dependencies.

I have recently added LibraryLink Support on mathematica++. It is still in develop branch. It'll take a few days before releasing it to the master branch.

It already have features to stream WSPut/MLPut calls streamed through the mlink or wslink.

shell << Values(FindRoot(ArcTan(1000 * Cos(x)), List(x, 1, 2),  Rule("Method") = method));

And it can cast tokens fetched from mathematica into generic stl container, like std::vector, std::pair, boost::tuple etc..

Exploiting these two features the aim of the Library Link support is to wrap only the difficult parts while letting the programmer focus on the implementation functions. I started the development process with the following 4 concerns in mind.

  • Unharmed Skeleton The boilerplate of 4 functions, initialize, uninitialize, version and the actual function to be exported is left intact with the exact signature, so that the developer can easily follow up after reading the LibraryLink documentation on Mathematica website. It also doesn't restrict the developer to code whatever he/she could with plain C API.
  • C++ Environment The developer should be able to develop and unit test the module only using C++ based environment without using Mathematica fronted. However the developer is not barred from using the frontend in the middle of the development process. The concern is a C++ developer should feel home like while developing a module.
  • Support Generics standard C++ generic containers should be supported.
  • Extensible The library should be extensible to support 3rd party C++ libraries. (Eigen Matrix is already supported)

Following is an example how link based Library development is abstracted. The native_link can still be used with the C API. The resolver redirects calls to different impl functions either based on heads or based on number of arguments.

EXTERN_C DLLEXPORT int SomeFunctionX(WolframLibraryData libData, WMK_LINK native_link){
   mathematica::wtransport shell(libData, native_link);

   try{
        mathematica::resolver resolver(shell);
        resolver, overload(&some_function_impl_geo, shell) = {"GeoPosition", "GeoPosition"}
                , overload(&some_function_impl_complex) = {"Complex", "Complex"}
                , overload(&some_function_impl_binary)
                , overload(&some_function_impl_unary);
        return resolver.resolve();
    }catch(...){
        return shell.pass();
    }
    return 0;
}

Following are the signatures of the impl functions in the above example. overload checks the function argument types and cast accordingly. This can can handle user defined types like the point_2d in the example through different specialized serialization policies.

double some_function_impl_geo(mathematica::transport& shell, point_2d<double> p1, point_2d<double> p2);
double some_function_impl_complex(std::complex<double> p1, std::complex<double> p2)
int some_function_impl_binary(double x, double y);
mathematica::m some_function_impl_unary(double x);

The shell (first argument in the first overload) comes from libData and can be used to do computations with mathematica. In the example some_function_impl_geo converts point_2d to GeoPosition and then passes to GeoDistance as arguments.

double some_function_impl_geo(mathematica::transport& shell, point_2d<double> p1, point_2d<double> p2){    
    std::string unit("Kilometers");
    double res;
    shell << QuantityMagnitude(GeoDistance(p1, p2, Rule("!UnitSystem") = unit));
    shell >> res;
    return res;
}

The returned output is sent back to the frontend. The last overload in the example returns a mathematica expression.


As previously mentioned no flavours are abandoned. The MArgument based style is abstracted in mathematica::mtransport. mtransport doesn't support overloads, type of its input is expected to be fixed with LibraryFunctionLoad. Input output with tensor of any type, rank and dimension is supported. Type of nesting level of std::vector determines the expected type and rank of the tensor. Currently only std::vector is supported, however there are plans to support eigen and other containers.

args = shell is used to capture the input arguments, whereas shell = ret is used to get the returned value.

EXTERN_C DLLEXPORT int SomeFunctionMXT(WolframLibraryData libData, mint argc, MArgument* argv, MArgument res){
    mathematica::mtransport shell(libData, argc, argv, res);
    typedef std::vector<std::vector<double>> matrix_type;
    try{
        boost::tuple<matrix_type, matrix_type> args = shell; // fetch input
        matrix_type matl, matr, mato;   // declare variables to hold the input and output
        boost::tie(matl, matr) = args;  // tie tuple to variables
        shell << Dot(matl, matr);       // To lazy to write matrix multiplication code in C++ so using mathematica instead
                                  // This converts matl and matr to the corresponding List[..] expression
        shell >> mato;            // retrieve multiplication result as List[...] and cast to matrix_type
        shell = mato;             // return the result
    }catch(...){
        return shell.pass();
    }
    return LIBRARY_NO_ERROR;
}

shell.pass() passes the exceptions to the frontend and aborts.

I'd be glad to hear reviews and comments from others.

Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract