HISTCOUNTS a true replacement of HISTC?

Since few latest releases we get the code warning "histc is not recommended. Use histcounts instead".

However how to replace HISTC with specified DIM argument, for example what is the command to get C1 and C2 in this example?

I hope you won't tell me I need a for-loop or some sort disguised loop.

>> A=reshape(linspace(0,1,20),[4 5])

A =

           0    0.2105    0.4211    0.6316    0.8421
      0.0526    0.2632    0.4737    0.6842    0.8947
      0.1053    0.3158    0.5263    0.7368    0.9474
      0.1579    0.3684    0.5789    0.7895    1.0000

>> b=linspace(0,1,3)

b =

           0    0.5000    1.0000

>> c1=histc(A,b,1)

c1 =

       4     4     2     0     0
       0     0     2     4     3
       0     0     0     0     1

>> c2=histc(A,b,2)

c2 =

       3     2     0
       3     2     0
       2     3     0
       2     2     1

>>

Answer 1

[Bruno Luong]

Third solution, may be the best in term of speed/memory

    % replacement of [n,i] = histc(x,edges,dim);
    [~,~,i] = histcounts(x, [edges,Inf]);
    s = uint32(size(x));
    nbins = uint32(length(edges)); 
    p1 = uint32(prod(s(1:dim-1)));
    p2 = uint32(prod(s(dim+1:end)));
    if p1 == 1 % happens for dim==1
        k = reshape(uint32(i),[s(dim),p2]);
        l = nbins*(0:p2-1);
        ilin = l + k;
    else
        j = reshape(1:p1,[p1,1,1]);
        k = reshape(uint32(i-1)*p1,[p1,s(dim),p2]);
        l = reshape((p1*nbins)*(0:p2-1),[1,1,p2]);
        ilin = (j + l) + k;
    end
    ilin = ilin(:);
    s(dim) = nbins;
    n = accumarray(ilin(i~=0),1,[prod(s),1]);
    n = reshape(n,s);

EDIT: slight improve code by (1) avoid reduce calculation for dim=1, (2) indexing cast to UINT32

[Bruno Luong]

Sorry, I must Accept my own answer but I think it's fair. Promise, I'll accept another solution if it's better.

[Bruno Luong]

Timing of the above algo for arrays of 100x100x100

Working along dimension #1
  Time HISTC = 15.549 [ms]
  Time HISTCOUNTS = 19.481 [ms]
  Penalty factor = 1.3
Working along dimension #2
  Time HISTC = 20.104 [ms]
  Time HISTCOUNTS = 16.454 [ms]
  Penalty factor = 0.8
Working along dimension #3
  Time HISTC = 21.446 [ms]
  Time HISTCOUNTS = 18.485 [ms]
  Penalty factor = 0.9

The test code (in function) is as following:

x = reshape(linspace(0,1,1e6),[100 100 100]);
edges = linspace(0,1,100);
t = zeros(2,ndims(x));
for dim=1:ndims(x)
    tic
    [nn,~] = histc(x,edges,dim);
    t(1,dim) = toc;   
    tic
    [~,~,i] = histcounts(x, [edges,Inf]);
    s = uint32(size(x));
    nbins = uint32(length(edges)); 
    p1 = uint32(prod(s(1:dim-1)));
    p2 = uint32(prod(s(dim+1:end)));
    if p1 == 1
        k = reshape(uint32(i),[s(dim),p2]);
        l = reshape(nbins*(0:p2-1),1,[]);
        ilin = l + k;
    else
        j = reshape(1:p1,[],1,1);
        k = reshape(uint32(i-1)*p1,[p1,s(dim),p2]);
        l = reshape((p1*nbins)*(0:p2-1),1,1,[]);
        ilin = (j + l) + k;
    end
    ilin = ilin(:);
    s(dim) = nbins;
    n = accumarray(ilin(i~=0),1,[prod(s),1]);
    n = reshape(n,s);
    t(2,dim) = toc;
end
t = t*1e3;
for dim=1:ndims(x)
    fprintf('Working along dimension #%d\n', dim);
    fprintf('\tTime HISTC = %1.3f [ms]\n', t(1,dim));
    fprintf('\tTime HISTCOUNTS = %1.3f [ms]\n', t(2,dim));
    fprintf('\tPenalty factor = %1.1f\n', t(2,dim)/t(1,dim));
end

Answer 2

[Bruno Luong]

OK, here is the answer of my own question.

The replacement of

[n,i] = histc(x, edges, dim);

is

[~,~,i] = histcounts(x, edges);
s = size(x);
nd = length(s);
idx = arrayfun(@(n) 1:n, s, 'unif', 0);
[idx{:}] = ndgrid(idx{:});
idx{dim} = i;
s(dim) = length(edges);
idx = cat(nd+1,idx{:});
idx = reshape(idx,[],nd);
n = accumarray(idx(i~=0,:),1,s);

That's a lot of code. If an authority of TMW can shed a light if there is a shorter, faster, cleaner solution, I'm in.

NOTE: there is still a subtle difference caused by data right at the outer edge

[Stephen23]

Nice. This definitely wins the Answer of the Day award.

Perhaps the next step would be to make an enhancement request.

Answer 3

[OCDER]

Interesting finding that histcounts doesn't have the Dim option. You should ask TMW to add that feature - they must have overlooked that.

"I hope you won't tell me I need a for-loop or some sort disguised loop."

But, there really is no need to avoid a simpler for-loop solution. Also, your answer has a "disguised loop" inside the arrayfun - arrayfun uses an internal for loop.

A =[       0    0.2105    0.4211    0.6316    0.8421;
      0.0526    0.2632    0.4737    0.6842    0.8947;
      0.1053    0.3158    0.5263    0.7368    0.9474;
      0.1579    0.3684    0.5789    0.7895    1.0000]; 
Bin = linspace(0, 1, 3);

C1 = histc(A, Bin, 1);
C2 = myHistC(A, Bin, 1);

D1 = histc(A, Bin, 2);
D2 = myHistC(A, Bin, 2);

assert(isequal(C1, C2) && isequal(D1, D2), 'Not the same thing');

function C = myHistC(A, Bin, Dim)
    if nargin < 3; Dim = 1; end

      if Dim == 1
          C = zeros(numel(Bin), size(A, 2));
          for j = 1:size(A, 2)
              C(:, j) = histcounts(A(:, j), [Bin(:); Inf]);
          end
      else
          C = zeros(size(A, 1), numel(Bin));
          for j = 1:size(A, 1)
              C(j, :) = histcounts(A(j, :), [Bin(:); Inf]);
          end
      end
  end

[Bruno Luong]

What I try to avoid is the loop on HISTCOUNTS, the small parasite loops I don't mind.

BTW your code can't handle nd array as it is.

[OCDER]

Oh, I see. Yeah, TMW should step in and suggest a better solution (if any). Plus, my solution also cannot handle the variable outputs generated by histc/histcounts. I'm curious as to what's the alternative to histc.

[Guillaume]

"they must have overlooked that"

I'm not speaking for Mathworks, but I'd say they're very much aware of it. You just have to look at Replace Discouraged Instances of hist and histc which clearly states under Code Updates for histc:

"histcounts treats the input matrix as a single tall vector and calculates the bin counts for the entire matrix [...] Use a for-loop to calculate bin counts over each column. [...] If performance is a problem due to a large number of columns in the matrix, then consider continuing to use histc for the column-wise bin counts."

[OCDER]

I see. But if the suggestion is to use histc, I'm wondering why histcounts doesn't implement this directly using a name-value input like

histcounts(..., 'dim', N) ?

Overall, I find it confusing that they acknowledged the shortcoming of histcounts, offered a workaround code in the example but do not implement this within histcounts, and then, made histc generate warning messages to use histcounts instead, which has no replacement for histc's column/row-wise binning. It's an interesting decision.

From the website:

Use a for-loop to calculate bin counts over each column.

A = randn(100,10);
nbins = 10;
N = zeros(nbins, size(A,2));
for k = 1:size(A,2)
   N(:,k) = histcounts(A(:,k),nbins);
end

[Bruno Luong]

+1 OCDER

And the example in the document merely deals with a workaround for 2D case. It would be a nightmare for someone who is not familiar with MATLAB and find a decend alternative for nd-array case.

Answer 4

[Bruno Luong]

Another solution of replacement, still long, but consumes less memory than my first solution at the cost of array permutations:

    % replacement of [n,i] = histc(x,edges,dim);
    s = size(x);
    nd = length(s);
    p = 1:nd;
    p([1 dim]) = p([dim 1]);
    q(1:nd) = p;
    xp = permute(x, p);
    sp = s(p);
    [~,~,i] = histcounts(xp, [edges, Inf]);
    m = prod(sp(2:end));
    sp1 = sp(1);
    nbins = length(edges);
    j = repmat(nbins*(0:m-1),[sp1,1]);
    keep = i~=0;
    linidx = i(keep) + j(keep);
    sp(1) = nbins;
    n = accumarray(linidx(:),1,[prod(sp) 1]);
    n = permute(reshape(n,sp),q);
    i = permute(i,q);