Fast 1D BinCounts Alternative

Question

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.

Answer 1

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
    ]