accumarray

[David Bateman]

I have no idea if the attached function is really useful, but it is a
Matlab core function. I saw it as a request on the Wiki and thought it
would be fairly trivial to implement a version of it.

Matlab has this as a builtin function, so I'm sure that Matlab is
significantly faster for its version of this function. Thought I'd tried
to keep the code vectorized as much as I could. Include this function in
the CVS Jon if you want it..

Cheers
David

## Copyright (C) 2007 David Bateman
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; if not, write to the Free Software
## Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
## 02110-1301  USA

## -*- texinfo -*-
## @deftypefn {Function File} {} accumarray (@var{subs}, @var{vals}, @var{sz}, 
@var{fun}, @var{fillval}, @var{issparse})
## @deftypefnx {Function File} {} accumarray (@var{csubs}, @var{vals}, @dots{})
##
## Creates an array by accumulating the elements of a vector into the
## positions of defined by their subscripts. The subscripts are defined by
## the rows of the matrix @var{subs} and the values by @var{vals}. Each row
## of @var{subs} corresponds to one of the values in @var{vals}.
##
## The size of the matrix will be determined by the subscripts themselves.
## However, if @var{sz} is defined it determines the matrix size. The length
## of @var{sz} must correspond to the number of columns in @var{subs}.
##
## The default action of @code{accumarray} is to sum the elements with the
## same subscripts. This behavior can be modified by defining the @var{fun}
## function. This should be a function or function handle that accepts a 
## column vector and returns a scalar. The result of the function should not
## depend on the order of the subscripts.
##
## The elements of the returned array that have no subscripts assoicated with
## them are set to zero. Defining @var{fillval} to some other value allows
## these values to be defined.
##
## By default @code{accumarray} returns a full matrix. If @var{issparse} is
## logically true, then a sparse matrix is returned instead.
##
## An example of the use of @code{accumarray} is:
##
## @example
## @group
## accumarray ([1,1,1;2,1,2;2,3,2;2,1,2;2,3,2],101:105)
## @result{} ans(:,:,1) = [101, 0, 0; 0, 0, 0]
##    ans(:,:,2) = [0, 0, 0; 206, 0, 208]
## @end group
## @end example
## @end deftypefn

function A = accumarray (subs, val, sz, fun, fillval, isspar)  
  if (nargin < 2 || nargin > 6)
    print_usage ();
  endif

  if (iscell(subs))
    subs = cell2mat (cellfun (@(x) x(:), subs, 'UniformOutput', false));
  endif
  ndims = size (subs, 2);

  if (nargin < 3 || isempty (sz))
    sz = max (subs);
    if (isscalar(sz))
      sz = [sz, 1];
    endif
  else
    if (length (sz) != ndims)
      error ("accumarray: inconsistent dimensions");
    endif
  endif
  
  if (nargin < 4 || isempty (fun))
    fun = @sum;
  endif

  if (nargin < 5 || isempty (fillval))
    fillval = 0;
  endif

  if (nargin < 6 || isempty (isspar))
    isspar = false;
  endif
  if (isspar)
    if (ndims > 2)
      error ("Can not have more than 2 dimensions in a sparse matrix");
    endif
  endif

  [subs, idx] = sortrows (subs);
  if (isscalar (val))
    val = val * ones (size (idx));
  else
    val = val(idx);
  endif
  cidx = find([true; (sum (abs (diff (subs)), 2) != 0)]);
  idx = cell (1, ndims);
  for i = 1:ndims
    idx{i} = subs (cidx, i);
  endfor
  for i = 1 : length(cidx)-1
    x(i) = feval (fun, val(cidx(i) : cidx(i + 1) - 1));
  endfor
  x(length(cidx)) = feval (fun, val(cidx(end) : end));
  x = x(:);
  if (isspar && fillval == 0)
    A = sparse (idx{1}, idx{2}, x, sz(1), sz(2));
  else
    if (iscell (x))
      ## Why did matlab choose to reverse the order of the elements
      x = cellfun (@(x) flipud(x(:)), x, 'UniformOutput', false);
      A = cell (sz);
    elseif (fillval == 0)
      A = zeros (sz, class(x));
    else 
      A = fillval .* ones (sz);
    endif
    A (sub2ind (sz, idx{:})) = x;
  endif
endfunction

%!error (accumarray (1:5))
%!assert (accumarray ([1;2;4;2;4],101:105), [101;206;0;208])
%!assert (accumarray ([1,1,1;2,1,2;2,3,2;2,1,2;2,3,2],101:105),cat(3, 
[101,0,0;0,0,0],[0,0,0;206,0,208]))
%!assert (accumarray 
([1,1,1;2,1,2;2,3,2;2,1,2;2,3,2],101:105,[],@(x)sin(sum(x))),sin(cat(3, 
[101,0,0;0,0,0],[0,0,0;206,0,208])))
%!assert (accumarray ({[1 3 3 2 3 1 2 2 3 3 1 2],[3 4 2 1 4 3 4 2 2 4 3 4],[1 1 
2 2 1 1 2 1 1 1 2 
2]},101:112),cat(3,[0,0,207,0;0,108,0,0;0,109,0,317],[0,0,111,0;104,0,0,219;0,103,0,0]))
%!assert (accumarray 
([1,1;2,1;2,3;2,1;2,3],101:105,[2,4],@max,NaN),[101,NaN,NaN,NaN;104,NaN,105,NaN])
%!assert (accumarray ([1 1; 2 1; 2 3; 2 1; 2 3],101:105,[2 
4],@prod,0,true),sparse([1,2,2],[1,1,3],[101,10608,10815],2,4))
%!assert (accumarray ([1 1; 2 1; 2 3; 2 1; 2 3],1,[2,4]), [1,0,0,0;2,0,2,0])
%!assert (accumarray ([1 1; 2 1; 2 3; 2 1; 2 
3],101:105,[2,4],@(x)length(x)>1),[false,false,false,false;true,false,true,false])
%!test
%! A = accumarray ([1 1; 2 1; 2 3; 2 1; 2 3],101:105,[2,4],@(x){x});
%! assert (A{2},[104;102])