23
$\begingroup$

When we arrive at V13.1, we have a great LibraryFunctionDeclaration now, which can call C function in MMA as this article:

Needs["CCompilerDriver`"]
lib = CCompilerDriver`CreateLibrary["#include \"WolframLibrary.h\"
#include \"WolframCompileLibrary.h\"

template<typename T>
void selection_sort(T* a, int low, int high)
{
int min;
    for (int i = low; i <= high; i++)
    {
    min = i;
    for (int j = i + 1; j <= high; j++)
    {
    if (a[j] < a[min])
    min = j;
    }
    if (min != i)
    {
    auto temp = a[min];
    a[min] = a[i];
    a[i] = temp;
    }
    }
}

EXTERN_C DLLEXPORT MTensor _selectionSort_Cpp(MTensor in)
{
auto in_data = (mreal*)in->data;
(in->refcount)++;
selection_sort(in_data, 0, in->nelems - 1);
return in;
}", "selectionSort", Language -> "C++"]
dec = LibraryFunctionDeclaration[
   selectionSort -> "_selectionSort_Cpp", 
   lib, {"PackedArray"["Real64", 1]} -> "PackedArray"["Real64", 1]];
fc = FunctionCompile[dec, 
  Function[Typed[n, "PackedArray"["Real64", 1]], selectionSort[n]]]

Then we can use C function(selection_sort) directly in MMA now:

lis = RandomReal[{-10, 10}, 100000];
fc[lis]; // Timing

{1.54688, Null}

It is awesome. But since we can use the Language -> "C++". Why not use C++ grammar? Such as we have a function sort_add2 like in the following:

#include<iostream>
#include<vector>
#include <algorithm>

using namespace std;

vector<int> sort_add2(vector<int> vec) {
    sort(vec.begin(), vec.end());
    for (int &i : vec) {
        i=i+2;
    }
    return vec;
}

int main()
{
    vector<int> vec = {5,9,20,3};
    vector<int> new_vec = sort_add2(vec);
    for(int i : new_vec){
        cout << i << endl;
    }

    return 0;
}

How to call that sort_add2 directly in MMA 13.1? Is it enough to construct a structure with MMA's existing compiled type so that MMA can recognize the input and output of sort_add2? Is possible?

$\endgroup$
5
  • 1
    $\begingroup$ Are you editing the question just to bump it to the front page? I don't see any actual content changes in the 6 edits made so far. If you want to motivate a response, perhaps you can place a bounty on the question? $\endgroup$
    – Jason B.
    Jul 26, 2022 at 19:53
  • $\begingroup$ @JasonB. Hello. The first three edits were to fix minor issues in the post. The last three edits were really just to activate it wanting more potential answerers to see them. bounty can only be set after 3 days. But the question is important to me and it's popular judging by the votes, so I had to use this dumb method $\endgroup$
    – yode
    Jul 27, 2022 at 1:41
  • 2
    $\begingroup$ No worries, I was just “busting your chops”. The question is certainly relevant to my interests and an answer would be of great utility $\endgroup$
    – Jason B.
    Jul 27, 2022 at 3:00
  • 1
    $\begingroup$ At least be cleverer than this: mathematica.stackexchange.com/posts/228098/revisions $\endgroup$
    – xzczd
    Jul 27, 2022 at 8:59
  • 1
    $\begingroup$ @xzczd I noticed that you asked about how to compile I Ragged List when using the old version of compile, now my last update has obviously solved this problem well. :) $\endgroup$
    – yode
    Aug 2, 2022 at 6:02

2 Answers 2

14
+500
$\begingroup$

This post is written with Mathematica version 13.1 and is subject to change in the future.

Mathematica's C Types does not include vector type but includes CArray which can be converted to C vector type easily.

Compiling, in general, and in this case, will change symbol names that's why you used EXTERN_C DLLEXPORT which are extern "C" and __declspec(dllexport) respectively to make sure they're exported and their names won't change. So I added a little bit of interface to make the code work.

  1. Let's assume our code looks like this:
cppCode = "
#include <vector>
#include <algorithm>

using namespace std;

// accept a reference to change the original data
void sort_add2(vector<int>& vec) {
    
    sort(vec.begin(), vec.end());
    for (int &i : vec) {
        i=i+2;
    }
}

extern \"C\" {
    __declspec(dllexport) int* sort_add2_interface(int* arr,int len){
        
        // build a vector from our array
        vector<int> temp(arr, arr + len);

        // do the operation
        sort_add2(temp);

        // return the first element address
        return &(temp[0]);
    }
}
";
  1. Now build the code:
lib = CCompilerDriver`CreateLibrary[cppCode, "selectionSort", 
  Language -> "C++"]
  1. Declare the function:
declaration = 
 LibraryFunctionDeclaration["sort_add2_interface", 
  lib, {"CArray"::["CInt"], "CInt"} -> "CArray"::["CInt"]]
  1. Load the exported library:
LibraryLoad[lib];
  1. On the Mathematica side, we convert our list to CArray using CreateTypeInstance, then call the library and build the result from reading the output with FromRawPointer (it get casts to Integer64).
fn = FunctionCompile[declaration, 
  Function[Typed[arr, "ListVector"::["Integer32"]], 
   Module[{temp, len, result},

    len = Cast[Length[arr], "CInt", "BitCast"];

    temp = 
     LibraryFunction["sort_add2_interface"][
      CreateTypeInstance["Managed"::["CArray"::["CInt"]], arr], len];

    Table[Cast[FromRawPointer[temp, i - 1], "Integer64"], {i, len}]

    ]
   ]
  ]

Testing:

fn[{2, 3, 1}]

(* Out: {3, 4, 5} *)

Remember that if you have a library that contains sort_add2 and is exposed, you can build another library that acts as an interface. In this case, I combine the two but they can be separated.

If they are separated, you have to include their declaration to your code and use CCompilerDriver`CreateLibrary's "Libraries" option to make it work.

For example, if you have a lib1.dll (with a lib1.lib) which contains sort_add2, you code will look like this:

cppCode2 = "
#include <vector>
#include <algorithm>

using namespace std;

// declaration
extern \"C\" void sort_add2(vector<int>& vec);

extern \"C\" __declspec(dllexport) int* sort_add2_interface(int* \
arr,int len){
    vector<int> temp(arr,arr+len);
    sort_add2(temp);
    return &(temp[0]);
}
";

You would use this code to compile:

lib = CCompilerDriver`CreateLibrary[cppCode2, "library_interface", 
  Language -> "C++", "Libraries" -> {"C:\\PATH_TO_LIBRARY\\lib1.lib"}]

And the rest is similar (follow from step 3, in step 4 you need to load both libraries).


I was more curious about the performance, so here is the result (assuming data = RandomInteger[{0, 2^31 - 1}, 10000]):

Code Timing Memory Usage
Sort[data]+2 0.00025 160,456
fn[data] 0.0013 (5x slower) 281,160

I decide to go deeper and this is the result of different parts of calling fn1:

Code Timing (~)
building vector (CreateTypeInstance) 0.00083
calling sort_add2 0.00048
building the result (Table[ ... FromRawPointer, ... ]) 0.00001

The newly added feature has much more power and along with the LibraryLink which can give direct access to the data (and I assume is faster), Mathematica become more connected than ever.

There are many things that I don't and I'm sure the code above can improve further. (I'll welcome any suggestion/correction)

Update 1

Thanks to @yode's comment and this post using "PackedArray" instead of "ListVector" will boost the speed by a factor of 2, so the new code will look like this:

fn2 = FunctionCompile[declaration, 
 Function[Typed[x, "PackedArray"["MachineInteger", 1]],
  Module[{in, len, out},
   len = Cast[Length[x], "CInt", "BitCast"];
   in = CreateTypeInstance["CArray"["CInt"], len];
   Do[ToRawPointer[in, i - 1, Cast[x[[i]], "CInt", "BitCast"]], {i, 
     len}];
   out = LibraryFunction["sort_add2_interface"][in, len];
   Table[Cast[FromRawPointer[out, i - 1], "Integer64"], {i, len}]]]]

In pursuing performance I did try sort_add2 with the i64 type, so we might gain speed by removing all the Cast but it didn't have an effect.

My next attempt was building a library and hoping by using LibraryFunctionLoad, we could gain speed which didn't happen (both Automatic and "Shared" memory management types). It was around ~20% slower than fn2. Interestingly whether it was "Shared" (in-place operation) or it was a copy of the original data, the speed didn't have a difference.

$\endgroup$
5
  • 2
    $\begingroup$ As far as I know, the ListVector is a time killer. If you change your "ListVector"::["Integer32"] into "PackedArray"::["MachineInteger", 1], your code will get a great improvement $\endgroup$
    – yode
    Jul 28, 2022 at 3:08
  • 1
    $\begingroup$ I think my following self-answer solves my problem actually. It's much more concise and efficient. The only thing I'm confused about is when I wake up and all users vote you but not me. :) But I really feel your warm reply, so I give you all the bounty. :) Thank you very much. $\endgroup$
    – yode
    Jul 28, 2022 at 3:10
  • 2
    $\begingroup$ FunctionCompile[declaration,Typed[x,"PackedArray"["MachineInteger",1]]|->Module[{in=CreateTypeInstance["CArray"["CInt"],Length@x],len,out},len=Cast[Length[x],"CInt","BitCast"]; Do[ToRawPointer[in,i-1,Cast[x[[i]],"CInt","BitCast"]],{i,Length@x}]; out=LibraryFunction["sort_add2_interface"][in,len]; Table[Cast[FromRawPointer[out,i-1],"Integer64"],{i,Length@x}]]] more faster from the author in my post link... $\endgroup$
    – yode
    Jul 28, 2022 at 14:50
  • $\begingroup$ @yode Thanks for your kind words. I really appreciate both your suggestion and bounty ;) I should also point out that the difficulty in the interface was due to the Wolfram libraries that I didn't use. So I have to build the array interface. But the performance was really a suprised. $\endgroup$
    – Ben Izd
    Jul 29, 2022 at 16:46
  • $\begingroup$ I don't know how to remove those Cast...Can you add an example? $\endgroup$
    – yode
    Aug 7, 2022 at 9:55
8
$\begingroup$

This is my program, in fact most of it comes from an unknown author. And it exactly is what I'm after. For more details you need to see follow this article for updates


Needs["CCompilerDriver`"]
lib = CCompilerDriver`CreateLibrary["
#include <vector>
#include <algorithm>
#include <WolframCompileLibrary.h>

using namespace std;

vector<mint> sort_add2(vector<mint> vec)
{
    sort(vec.begin(), vec.end());
    for (mint &i : vec)
    {
        i = i + 2;
    }
    return vec;
}

EXTERN_C DLLEXPORT MTensor sort_add2_fc(MTensor inmma)
{
    auto pinmma = (mint*)inmma->data;//or MTensor_getIntegerDataMacro(inmma);
    vector<mint> inc(pinmma, pinmma + inmma->nelems);
    vector<mint> outc = sort_add2(inc);
    copy(outc.begin(),outc.end(), pinmma);
    (inmma->refcount)++;
    return inmma;
}", "sortAdd2File", Language -> "C++"];

dec = LibraryFunctionDeclaration[fun -> "sort_add2_fc", 
   lib, {"PackedArray"::["MachineInteger", 1]} -> 
    "PackedArray"::["MachineInteger", 1]];
fc = FunctionCompile[dec, 
  Function[Typed[n, "PackedArray"::["MachineInteger", 1]], fun[n]]]

Then we can use it now:

fc[{2, 3, 1}]

{3, 4, 5}


Performace with Ben Izd's fn:

Needs["GeneralUtilities`"]
GeneralUtilities`BenchmarkPlot[<|"Sort[lis]+2" -> (Sort[#] + 2 &), 
  "yode" -> fc, "Ben Izd" -> fn|>, RandomInteger[100, #] &, 
 Range[100, 1000000, 1000], "IncludeFits" -> True]

enter image description here


Add another case to calculate the mean for a ragged list here. As the direction of the author, we can greatly boost the speed by a factor of 5 than the uncompiled solution. It is a significant dig, so I share it here for the community.

The huge trouble is the c++ cannot receive the "ListVector"["PackedArray"["Real64", 1]], we can see the compile string with FunctionCompileExportString for the format:

FunctionCompileExportString[
 Typed[x, "ListVector"["PackedArray"["Real64", 1]]] |-> x]

line6~9 enter image description here

We have built a struct receive the parameter from MMA in c++ like:

struct ManagedArray {
    MTensor* data;
    mint length;
};

struct FixedArray{
    ManagedArray data;
    mint refcount;
};

struct RaggedTensor{
    FixedArray* data;
    mint refcount;
};

So the memory for both the parameters and the arguments are aligned. Here is the complete code:

Needs["CCompilerDriver`"]
lib=CCompilerDriver`CreateLibrary["
#include<vector>
#include<numeric>
#include<WolframCompileLibrary.h>

using namespace std;

struct ManagedArray {
    MTensor* data;
    mint length;
};

struct FixedArray {
    ManagedArray data;
    mint refcount;
};

struct RaggedTensor {
    FixedArray* data;
    mint refcount;
public:
    mint size() {
        return data->data.length;
    }

    MTensor* begin() {
        return data->data.data;
    }

    MTensor* end() {
        return this->begin() + this->size();
    }
    vector<vector<mreal>> toVecs() {
        vector<vector<mreal>> vecs;
        for (auto tensor : *this) {
            auto ptensor = (mreal*)tensor->data;
            vecs.push_back(vector<mreal>(ptensor, ptensor + tensor->nelems));
        }
        return vecs;
    }
};

mreal mean(vector<vector<mreal>> vecs) {
    vector<mreal> tmp;
    for (vector<mreal> vec : vecs) {
        tmp.insert(tmp.end(), vec.begin(), vec.end());
    }
    return accumulate(tmp.begin(), tmp.end(), 0.) / tmp.size();
}

EXTERN_C DLLEXPORT mreal ragged_mean(RaggedTensor in) {
    return mean(in.toVecs());
}","Ragged_Mean",Language->"C++"];

dec=LibraryFunctionDeclaration[fun->"ragged_mean",lib,{"ListVector"::["PackedArray"::["Real64",1]]}->"Real64"];
fc=FunctionCompile[dec,Function[{Typed[lis,"ListVector"::["PackedArray"::["Real64",1]]]},fun[lis]]]

We can use it now:

raglist = N[{{5, 9, 3}, {7, 9}, {9, 11, 8, 33}}];
fc[Flatten[raglist], Length /@ raglist]

10.4444

Performace:

lis=DeleteCases[Table[ResourceFunction["RandomSplit"][RandomReal[255, 1000], 10], 
    100000], {}, {2}];

RepeatedTiming[fc /@ lis;]
RepeatedTiming[Mean@*Flatten /@ lis;]

{0.673275, Null}

{2.80634, Null}

$\endgroup$
1
  • $\begingroup$ When compiling "dec = LibraryFunctionDeclaration[fun -> "sort_add2_fc", lib, {"PackedArray"::["MachineInteger", 1]} -> "PackedArray"::["MachineInteger", 1]];", I got this error message "Syntax::sntxf: "{" cannot be followed by "PackedArray::[MachineInteger,1]}".". Could you tell me why? Thank you. $\endgroup$
    – xpaul
    Mar 22 at 18:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.