5
$\begingroup$

I have lots of data which looks like this example:

data = Sort@Flatten[{SeedRandom[42]; RandomReal[5, 2^8 - 2^2],     
                     RandomReal[25, 2^2] + 5}];

I need a binning function which is as fast as possible. In addition to the data, the binning function should have a binwidth argument and should output all frequencies up to a max number. The first bin interval is $0$ to binwidth. For the dataset data, I use binwidth=2^-1 and max=32. In total, the output should be a list of length max/binwidth. In short, the binning function should behave like

BinCounts[data, {0, 32, 2^-1}]

So I searched this site and the web and found the following:

ClearAll[myBinCounts, myBinCounts2, myBinCounts3]
myBinCounts[data_, binwidth_, max_] := 
  Module[{dat = Floor[1 + data/binwidth], res}, 
    System`SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> 1}];
    res = SparseArray[
            Flatten@{dat, max/binwidth} -> Flatten@{Table[1, {Length[dat]}], 0}];
    System`SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> 0}]; 
    Normal@res]
myBinCounts2[data_, binwidth_, max_] := 
  Module[{s = SortBy[Tally@Quotient[data, binwidth], First], num = Floor[max/binwidth], res},
    res = ConstantArray[0, num]; 
    Part[res, s[[All, 1]] + 1] = s[[All, 2]]; res]
myBinCounts3[data_, binwidth_, max_] := 
  Module[{s = Sort[Tally@Quotient[data, binwidth]], num = Floor[max/binwidth], res},
    res = ConstantArray[0, num]; 
    Part[res, s[[All, 1]] + 1] = s[[All, 2]]; res]

The idea of myBinCounts is from mathematica-fast-2d-binning-algorithm, the idea for myBinCounts2 from Szabolcs in this thread. The latter design is about 5 times faster then the former for this problem size. So I wrote compileable code and substitutet ConstantArray with Table (and SortBy by Sort from myBinCounts2 to myBinCounts3).

ClearAll[CmyBinCounts]
CmyBinCounts = 
  Compile[{{data, _Real, 1}, {binwidth, _Real, 0}, {max, _Integer, 0}},
    Module[{s = Sort[Tally@Quotient[data, binwidth]], num = Floor[max/binwidth], res},
      res = Table[0, {num}]; Part[res, s[[All, 1]] + 1] = s[[All, 2]]; res],
    CompilationTarget -> "C",
    Parallelization -> False,
    (*RuntimeAttributes -> {Listable},*)
    RuntimeOptions -> {"Speed", "EvaluateSymbolically" -> False}
  ]

The compiled function does not call MainEvaluate anymore:

StringFreeQ[CompiledFunctionTools`CompilePrint@CmyBinCounts, "MainEvaluate"]

The original Mathematica function BinCounts and my versions all give the same output:

1 == Length@DeleteDuplicates@FlattenAt[{BinCounts[data, {0, 32, 2^-1}], 
       Table[Thread[f[data, 2^-1, 32]], {f, 
         {myBinCounts, myBinCounts2, myBinCounts3, CmyBinCounts}
         }]}, 2]
(* True *)

Timing all versions, I get on my Windows 8 PC with CPU i7-2600 and MMA 10:

t = With[{k = 10(*adjust to your CPU*)}, 
      FlattenAt[{BinCounts[data, {0, 32, 2^-1}]~Do~{2^k}//AbsoluteTiming//First, 
        Table[Thread[f[data, 2^-1, 32]]~Do~{2^k}//AbsoluteTiming//First, {f,
          {myBinCounts, myBinCounts2, myBinCounts3, CmyBinCounts}
        }]}, 2]]
t/Min[t]
(* {0.207138, 0.173115, 0.041027, 0.029019, 0.011007} *)
(* {18.82, 15.73, 3.727, 2.636, 1.000} *)

PackedArrays are fine with me, I am after the fastest solution. It seems compiled code with no explicit SparseArray is fastest, but I am happy to learn. Changing SparseArrayOptions every time seems a waste of time. But I couldn't get a function to run with localized variables and the option changed globally (and my attempts were not much faster).

PS: I am relatively new to Mathematica, I am using it for about 1 month now. If there are some major drawbacks in the code or the way I program, please let me know. Still trying to understand all the different concepts, this site is a great learning resource.

$\endgroup$
  • $\begingroup$ I wonder if it would be worth writing this in C instead, possibly parallelized with OpenMP. $\endgroup$ – Szabolcs Oct 7 '15 at 12:05
  • 1
    $\begingroup$ Your sample list is very small, with only 256 entries. Do you need a binning function that will be run many times on small lists (thus overhead is important), or one that will be run only a few times on a long list (thus overhead is negligible)? $\endgroup$ – Szabolcs Oct 7 '15 at 12:07
  • $\begingroup$ Overhead is important, I have to run the binning function many times on small lists of about 100-500 entries. About 500 list entries is the longest which makes sense in my specific setting. $\endgroup$ – Marco Breitig Oct 7 '15 at 12:13
3
$\begingroup$

Here's a C++ implementation using LTemplate. I'm using LTemplate because it made it easy enough to write the code that I didn't give up before starting ;-)

<< LTemplate`

SetDirectory[$TemporaryDirectory]; (* currently LTemplate writes and reads files to/from the current directory *)
code = "
  #include <cmath>

  struct Binner {
    mma::IntTensorRef bin(mma::RealTensorRef t, double binwidth, double max) {
        mint n = std::ceil(max/binwidth);
        mma::IntTensorRef res = mma::makeVector<mint>(n);
        std::fill(res.begin(), res.end(), 0);
        for (double *i = t.begin(); i != t.end(); ++i) {
            mint b = std::floor((*i)/binwidth);
            if (0 <= b && b < n)
                res[b]++;   
        }
        return res;
    }
  };
  ";

Export["Binner.h", code, "String"];

template = 
  LClass["Binner", {LFun["bin", {{Real, 1, "Constant"}, Real, Real}, {Integer, 1}]}];

CompileTemplate[template]
LoadTemplate[template]

Here's the function to call:

binner = Make["Binner"]; (* create object once, and re-use it later, to reduce overhead *)
binCountsSz[data_, binwidth_, max_] := binner@"bin"[data, binwidth, max]

Let's test it:

data = Sort@Flatten[{SeedRandom[42]; RandomReal[5, 2^8 - 2^2], RandomReal[25, 2^2] + 5}];

Measure:

TimeIt@binCountsSz[data, 0.5, 32]
(* 2.61037*10^-6 *)

TimeIt@myBinCounts3[data, 0.5, 32]
(* 0.0000208091 *)

TimeIt@CmyBinCounts[data, 0.5, 32]
(* 9.53201*10^-6 *)

binCountsSz[data, 0.5, 32] == CmyBinCounts[data, 0.5, 32]
(* True *)

It does about 3.5 times better than the Compile version. To be fair, it's also some 3 times longer ... but still fairly short.

If you write it using pure LibraryLink instead of LTemplate, the overhead may be reduced further. I haven't tested this for this particular application.


TimeIt is something I use for benchmarking occasionally. It evaluates the expression a sufficient number of times that the timing is at least 1 second.

SetAttributes[TimeIt, HoldAll]
TimeIt[expr_, duration_ : 1.] :=
    Module[{t = 0., n = 1/2, d = duration},
      While[t < d,
        n *= 2;
        t = First@AbsoluteTiming@Do[expr, {n}]
      ];
      t/n
    ]
$\endgroup$
  • $\begingroup$ Thank you very much. I got it working, but now I have to understand LTemplates and how to use them in the future. On my PC, your version is about 1.8 to 1.9 times faster than my compiled code (which derived from your idea anyway), but I am happy about every speed-up I can get. :-) $\endgroup$ – Marco Breitig Oct 7 '15 at 13:49
  • 1
    $\begingroup$ @MarcoBreitig What compiler are you using? MSVC should be good on Windows. Let me know what you find difficult to understand about LTemplate and I'll improve the documentation. As far as I know, no one (but me) has used it before, so I wasn't terribly motivated to document it well. $\endgroup$ – Szabolcs Oct 7 '15 at 14:14
  • $\begingroup$ I used MinGWCompiler but reading your comment changed to MSVC 14. But compairing speed the difference seems to be minimal and negligible. Do you have some tips or a benchmark case on how to best compare different compilers? $\endgroup$ – Marco Breitig Oct 7 '15 at 15:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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