Mercurial > hg > octave-nkf
annotate scripts/general/accumdim.m @ 10396:a0b51ac0f88a
optimize accumdim with summation
| author | Jaroslav Hajek <highegg@gmail.com> |
|---|---|
| date | Fri, 05 Mar 2010 12:31:30 +0100 |
| parents | aeb5b1e47978 |
| children | 180931276a52 |
| rev | line source |
|---|---|
| 10395 | 1 ## Copyright (C) 2010 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{sz}, @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. | |
| 26 ## | |
| 27 ## The extent of the result matrix in the working dimension will be determined | |
| 28 ## by the subscripts themselves. | |
| 29 ## However, if @var{n} is defined it determines this extent. | |
| 30 ## | |
| 31 ## The default action of @code{accumdim} is to sum the subarrays with the | |
| 32 ## same subscripts. This behavior can be modified by defining the @var{func} | |
| 33 ## function. This should be a function or function handle that accepts an | |
| 34 ## array and a dimension, and reduces the array along this dimension. | |
| 35 ## As a special exception, the built-in @code{min} and @code{max} functions | |
| 36 ## can be used directly, and @code{accumdim} accounts for the middle empty | |
| 37 ## argument that is used in their calling. | |
| 38 ## | |
| 39 ## The slices of the returned array that have no subscripts associated with | |
| 40 ## them are set to zero. Defining @var{fillval} to some other value allows | |
| 41 ## these values to be defined. | |
| 42 ## | |
| 43 ## An example of the use of @code{accumarray} is: | |
| 44 ## | |
| 45 ## @example | |
| 46 ## @group | |
| 47 ## accumarray ([1,1,1;2,1,2;2,3,2;2,1,2;2,3,2], 101:105) | |
| 48 ## @result{} ans(:,:,1) = [101, 0, 0; 0, 0, 0] | |
| 49 ## ans(:,:,2) = [0, 0, 0; 206, 0, 208] | |
| 50 ## @end group | |
| 51 ## @end example | |
| 52 ## | |
| 53 ## @seealso{accumarray} | |
| 54 ## @end deftypefn | |
| 55 | |
| 56 function A = accumdim (subs, val, dim, n = 0, func = [], fillval = 0) | |
| 57 | |
| 58 if (nargin < 2 || nargin > 5) | |
| 59 print_usage (); | |
| 60 endif | |
| 61 | |
| 62 if (isempty (fillval)) | |
| 63 fillval = 0; | |
| 64 endif | |
| 65 | |
| 66 if (! isvector (subs)) | |
| 67 error ("accumdim: subs must be a subscript vector"); | |
| 68 elseif (! isindex (subs)) # creates index cache | |
| 69 error ("accumdim: indices must be positive integers"); | |
| 70 else | |
| 71 m = max (subs); | |
| 72 if (n == 0) | |
| 73 n = m; | |
| 74 elseif (n < m) | |
| 75 error ("accumdim: index out of range") | |
| 76 endif | |
| 77 endif | |
| 78 | |
| 79 sz = size (val); | |
| 80 | |
| 81 if (nargin < 3) | |
| 82 [~, dim] = max (sz != 1); # first non-singleton dim | |
|
10396
a0b51ac0f88a
optimize accumdim with summation
Jaroslav Hajek <highegg@gmail.com>
parents:
10395
diff
changeset
|
83 elseif (! isindex (dim)) |
| 10395 | 84 error ("accumdim: dim must be a valid dimension"); |
|
10396
a0b51ac0f88a
optimize accumdim with summation
Jaroslav Hajek <highegg@gmail.com>
parents:
10395
diff
changeset
|
85 elseif (dim > length (sz)) |
|
a0b51ac0f88a
optimize accumdim with summation
Jaroslav Hajek <highegg@gmail.com>
parents:
10395
diff
changeset
|
86 sz(end+1:dim) = 1; |
| 10395 | 87 endif |
| 88 sz(dim) = n; | |
| 89 | |
|
10396
a0b51ac0f88a
optimize accumdim with summation
Jaroslav Hajek <highegg@gmail.com>
parents:
10395
diff
changeset
|
90 if (isempty (func) || func == @sum) |
|
a0b51ac0f88a
optimize accumdim with summation
Jaroslav Hajek <highegg@gmail.com>
parents:
10395
diff
changeset
|
91 ## Fast summation case. |
|
a0b51ac0f88a
optimize accumdim with summation
Jaroslav Hajek <highegg@gmail.com>
parents:
10395
diff
changeset
|
92 A = __accumdim_sum__ (subs, val, dim, n); |
|
a0b51ac0f88a
optimize accumdim with summation
Jaroslav Hajek <highegg@gmail.com>
parents:
10395
diff
changeset
|
93 |
|
a0b51ac0f88a
optimize accumdim with summation
Jaroslav Hajek <highegg@gmail.com>
parents:
10395
diff
changeset
|
94 ## Fill in nonzero fill value |
|
a0b51ac0f88a
optimize accumdim with summation
Jaroslav Hajek <highegg@gmail.com>
parents:
10395
diff
changeset
|
95 if (fillval != 0) |
|
a0b51ac0f88a
optimize accumdim with summation
Jaroslav Hajek <highegg@gmail.com>
parents:
10395
diff
changeset
|
96 mask = true (n, 1); |
|
a0b51ac0f88a
optimize accumdim with summation
Jaroslav Hajek <highegg@gmail.com>
parents:
10395
diff
changeset
|
97 mask(subs) = false; |
|
a0b51ac0f88a
optimize accumdim with summation
Jaroslav Hajek <highegg@gmail.com>
parents:
10395
diff
changeset
|
98 subsc = {':'}(ones (1, length (sz))); |
|
a0b51ac0f88a
optimize accumdim with summation
Jaroslav Hajek <highegg@gmail.com>
parents:
10395
diff
changeset
|
99 subsc{dim} = mask; |
|
a0b51ac0f88a
optimize accumdim with summation
Jaroslav Hajek <highegg@gmail.com>
parents:
10395
diff
changeset
|
100 A(subsc{:}) = fillval; |
|
a0b51ac0f88a
optimize accumdim with summation
Jaroslav Hajek <highegg@gmail.com>
parents:
10395
diff
changeset
|
101 endif |
|
a0b51ac0f88a
optimize accumdim with summation
Jaroslav Hajek <highegg@gmail.com>
parents:
10395
diff
changeset
|
102 return |
|
a0b51ac0f88a
optimize accumdim with summation
Jaroslav Hajek <highegg@gmail.com>
parents:
10395
diff
changeset
|
103 endif |
|
a0b51ac0f88a
optimize accumdim with summation
Jaroslav Hajek <highegg@gmail.com>
parents:
10395
diff
changeset
|
104 |
|
a0b51ac0f88a
optimize accumdim with summation
Jaroslav Hajek <highegg@gmail.com>
parents:
10395
diff
changeset
|
105 ## The general case. |
| 10395 | 106 ns = length (subs); |
| 107 ## Sort indices. | |
| 108 [subs, idx] = sort (subs(:)); | |
| 109 ## Identify runs. | |
| 110 jdx = find (subs(1:ns-1) != subs(2:ns)); | |
| 111 jdx = [jdx; ns]; | |
| 112 ## Collect common slices. | |
| 113 szc = num2cell (sz); | |
| 114 szc{dim} = diff ([0; jdx]); | |
| 115 subsc = {':'}(ones (1, length (sz))); | |
| 116 subsc{dim} = idx; | |
| 117 val = mat2cell (val(subsc{:}), szc{:}); | |
| 118 ## Apply reductions. Special case min, max. | |
| 119 if (func == @min || func == @max) | |
| 120 val = cellfun (func, val, {[]}, {dim}, "uniformoutput", false); | |
| 121 else | |
| 122 val = cellfun (func, val, {dim}, "uniformoutput", false); | |
| 123 endif | |
| 124 subs = subs(jdx); | |
| 125 | |
| 126 ## Concatenate reduced slices. | |
| 127 val = cat (dim, val{:}); | |
| 128 | |
| 129 ## Construct matrix of fillvals. | |
| 130 if (fillval == 0) | |
| 131 A = zeros (sz, class (val)); | |
| 132 else | |
| 133 A = repmat (fillval, sz); | |
| 134 endif | |
| 135 | |
| 136 ## Set the reduced values. | |
| 137 subsc{dim} = subs; | |
| 138 A(subsc{:}) = val; | |
| 139 | |
| 140 endfunction | |
| 141 | |
| 142 %%test accumdim vs. accumarray | |
| 143 | |
| 144 %!shared a | |
| 145 %! a = rand (5, 5, 5); | |
| 146 | |
| 147 %!assert (accumdim ([1;3;1;3;3], a)(:,2,3), accumarray ([1;3;1;3;3], a(:,2,3))) | |
| 148 %!assert (accumdim ([2;3;2;2;2], a, 2, 4)(4,:,2), accumarray ([2;3;2;2;2], a(4,:,2), [1,4])) | |
| 149 %!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)) | |
| 150 %!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)) |