octave_packages/m/general/accumdim.m

   1 ## Copyright (C) 2010-2012 VZLU Prague
   2 ##
   3 ## This file is part of Octave.
   4 ##
   5 ## Octave is free software; you can redistribute it and/or modify it
   6 ## under the terms of the GNU General Public License as published by
   7 ## the Free Software Foundation; either version 3 of the License, or (at
   8 ## your option) any later version.
   9 ##
  10 ## Octave is distributed in the hope that it will be useful, but
  11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 ## General Public License for more details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with Octave; see the file COPYING.  If not, see
  17 ## <http://www.gnu.org/licenses/>.
  18
  19 ## -*- texinfo -*-
  20 ## @deftypefn {Function File} {} accumdim (@var{subs}, @var{vals}, @var{dim}, @var{n}, @var{func}, @var{fillval})
  21 ## Create an array by accumulating the slices of an array into the
  22 ## positions defined by their subscripts along a specified dimension.
  23 ## The subscripts are defined by the index vector @var{subs}.
  24 ## The dimension is specified by @var{dim}.  If not given, it defaults
  25 ## to the first non-singleton dimension.  The length of @var{subs} must
  26 ## be equal to @code{size (@var{vals}, @var{dim})}.
  27 ##
  28 ## The extent of the result matrix in the working dimension will be
  29 ## determined by the subscripts themselves.  However, if @var{n} is
  30 ## defined it determines this extent.
  31 ##
  32 ## The default action of @code{accumdim} is to sum the subarrays with the
  33 ## same subscripts.  This behavior can be modified by defining the
  34 ## @var{func} function.  This should be a function or function handle
  35 ## that accepts an array and a dimension, and reduces the array along
  36 ## this dimension.  As a special exception, the built-in @code{min} and
  37 ## @code{max} functions can be used directly, and @code{accumdim}
  38 ## accounts for the middle empty argument that is used in their calling.
  39 ##
  40 ## The slices of the returned array that have no subscripts associated
  41 ## with them are set to zero.  Defining @var{fillval} to some other
  42 ## value allows these values to be defined.
  43 ##
  44 ## An example of the use of @code{accumdim} is:
  45 ##
  46 ## @example
  47 ## @group
  48 ## accumdim ([1, 2, 1, 2, 1], [ 7, -10,   4;
  49 ##                             -5, -12,   8;
  50 ##                            -12,   2,   8;
  51 ##                            -10,   9,  -3;
  52 ##                             -5,  -3, -13])
  53 ## @result{} [-10,-11,-1;-15,-3,5]
  54 ## @end group
  55 ## @end example
  56 ##
  57 ## @seealso{accumarray}
  58 ## @end deftypefn
  59
  60 function A = accumdim (subs, vals, dim, n = 0, func = [], fillval = 0)
  61
  62   if (nargin < 2 || nargin > 5)
  63     print_usage ();
  64   endif
  65
  66   if (isempty (fillval))
  67     fillval = 0;
  68   endif
  69
  70   if (! isvector (subs))
  71     error ("accumdim: SUBS must be a subscript vector");
  72   elseif (! isindex (subs)) # creates index cache
  73     error ("accumdim: indices must be positive integers");
  74   else
  75     m = max (subs);
  76     if (n == 0 || isempty (n))
  77       n = m;
  78     elseif (n < m)
  79       error ("accumdim: N index out of range");
  80     endif
  81   endif
  82
  83   sz = size (vals);
  84
  85   if (nargin < 3)
  86     [~, dim] = max (sz != 1); # first non-singleton dim
  87   elseif (! isindex (dim))
  88     error ("accumdim: DIM must be a valid dimension");
  89   elseif (dim > length (sz))
  90     sz(end+1:dim) = 1;
  91   endif
  92   sz(dim) = n;
  93
  94   if (length (subs) != size (vals, dim))
  95     error ("accumdim: dimension mismatch")
  96   endif
  97
  98   if (isempty (func) || func == @sum)
  99     ## Fast summation case.
 100     A = __accumdim_sum__ (subs, vals, dim, n);
 101
 102     ## Fill in nonzero fill value
 103     if (fillval != 0)
 104       mask = true (n, 1);
 105       mask(subs) = false;
 106       subsc = {':'}(ones (1, length (sz)));
 107       subsc{dim} = mask;
 108       A(subsc{:}) = fillval;
 109     endif
 110     return
 111   endif
 112
 113   ## The general case.
 114   ns = length (subs);
 115   ## Sort indices.
 116   [subs, idx] = sort (subs(:));
 117   ## Identify runs.
 118   jdx = find (subs(1:ns-1) != subs(2:ns));
 119   jdx = [jdx; ns];
 120   ## Collect common slices.
 121   szc = num2cell (sz);
 122   szc{dim} = diff ([0; jdx]);
 123   subsc = {':'}(ones (1, length (sz)));
 124   subsc{dim} = idx;
 125   vals = mat2cell (vals(subsc{:}), szc{:});
 126   ## Apply reductions. Special case min, max.
 127   if (func == @min || func == @max)
 128     vals = cellfun (func, vals, {[]}, {dim}, "uniformoutput", false);
 129   else
 130     vals = cellfun (func, vals, {dim}, "uniformoutput", false);
 131   endif
 132   subs = subs(jdx);
 133
 134   ## Concatenate reduced slices.
 135   vals = cat (dim, vals{:});
 136
 137   ## Construct matrix of fillvals.
 138   if (fillval == 0)
 139     A = zeros (sz, class (vals));
 140   else
 141     A = repmat (fillval, sz);
 142   endif
 143
 144   ## Set the reduced values.
 145   subsc{dim} = subs;
 146   A(subsc{:}) = vals;
 147
 148 endfunction
 149
 150 %%test accumdim vs. accumarray
 151
 152 %!shared a
 153 %! a = rand (5, 5, 5);
 154
 155 %!assert (accumdim ([1;3;1;3;3], a)(:,2,3), accumarray ([1;3;1;3;3], a(:,2,3)))
 156 %!assert (accumdim ([2;3;2;2;2], a, 2, 4)(4,:,2), accumarray ([2;3;2;2;2], a(4,:,2), [1,4]))
 157 %!assert (accumdim ([2;3;2;1;2], a, 3, 3, @min)(1,5,:), accumarray ([2;3;2;1;2], a(1,5,:), [1,1,3], @min))
 158 %!assert (accumdim ([1;3;2;2;1], a, 2, 3, @median)(4,:,5), accumarray ([1;3;2;2;1], a(4,:,5), [1,3], @median))