Mercurial > octave-libgccjit
annotate scripts/general/accumarray.m @ 8934:c2099a4d12ea
partially optimize accumarray
| author | Jaroslav Hajek <highegg@gmail.com> |
|---|---|
| date | Mon, 09 Mar 2009 10:59:19 +0100 |
| parents | eb63fbe60fab |
| children | 1bf0ce0930be |
| rev | line source |
|---|---|
| 8920 | 1 ## Copyright (C) 2007, 2008, 2009 David Bateman |
|
8934
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
2 ## Copyright (C) 2009 VZLU Prague |
| 6770 | 3 ## |
| 7016 | 4 ## This file is part of Octave. |
| 6770 | 5 ## |
| 7016 | 6 ## Octave is free software; you can redistribute it and/or modify it |
| 7 ## under the terms of the GNU General Public License as published by | |
| 8 ## the Free Software Foundation; either version 3 of the License, or (at | |
| 9 ## your option) any later version. | |
| 10 ## | |
| 11 ## Octave is distributed in the hope that it will be useful, but | |
| 12 ## WITHOUT ANY WARRANTY; without even the implied warranty of | |
| 13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
| 14 ## General Public License for more details. | |
| 6770 | 15 ## |
| 16 ## You should have received a copy of the GNU General Public License | |
| 7016 | 17 ## along with Octave; see the file COPYING. If not, see |
| 18 ## <http://www.gnu.org/licenses/>. | |
| 6770 | 19 |
| 20 ## -*- texinfo -*- | |
| 21 ## @deftypefn {Function File} {} accumarray (@var{subs}, @var{vals}, @var{sz}, @var{fun}, @var{fillval}, @var{issparse}) | |
| 22 ## @deftypefnx {Function File} {} accumarray (@var{csubs}, @var{vals}, @dots{}) | |
| 23 ## | |
| 7186 | 24 ## Create an array by accumulating the elements of a vector into the |
| 25 ## positions defined by their subscripts. The subscripts are defined by | |
| 6770 | 26 ## the rows of the matrix @var{subs} and the values by @var{vals}. Each row |
| 27 ## of @var{subs} corresponds to one of the values in @var{vals}. | |
| 28 ## | |
| 29 ## The size of the matrix will be determined by the subscripts themselves. | |
| 30 ## However, if @var{sz} is defined it determines the matrix size. The length | |
| 31 ## of @var{sz} must correspond to the number of columns in @var{subs}. | |
| 32 ## | |
| 33 ## The default action of @code{accumarray} is to sum the elements with the | |
| 34 ## same subscripts. This behavior can be modified by defining the @var{fun} | |
| 35 ## function. This should be a function or function handle that accepts a | |
| 36 ## column vector and returns a scalar. The result of the function should not | |
| 37 ## depend on the order of the subscripts. | |
| 38 ## | |
|
8325
b93ac0586e4b
spelling corrections
Brian Gough<bjg@network-theory.co.uk>
parents:
7186
diff
changeset
|
39 ## The elements of the returned array that have no subscripts associated with |
| 6770 | 40 ## them are set to zero. Defining @var{fillval} to some other value allows |
| 41 ## these values to be defined. | |
| 42 ## | |
| 43 ## By default @code{accumarray} returns a full matrix. If @var{issparse} is | |
| 44 ## logically true, then a sparse matrix is returned instead. | |
| 45 ## | |
| 46 ## An example of the use of @code{accumarray} is: | |
| 47 ## | |
| 48 ## @example | |
| 49 ## @group | |
| 7186 | 50 ## accumarray ([1,1,1;2,1,2;2,3,2;2,1,2;2,3,2], 101:105) |
| 6770 | 51 ## @result{} ans(:,:,1) = [101, 0, 0; 0, 0, 0] |
| 52 ## ans(:,:,2) = [0, 0, 0; 206, 0, 208] | |
| 53 ## @end group | |
| 54 ## @end example | |
| 55 ## @end deftypefn | |
| 56 | |
| 57 function A = accumarray (subs, val, sz, fun, fillval, isspar) | |
| 58 | |
| 59 if (nargin < 2 || nargin > 6) | |
| 60 print_usage (); | |
| 61 endif | |
| 62 | |
| 7186 | 63 if (iscell (subs)) |
| 8507 | 64 subs = cell2mat (cellfun (@(x) x(:), subs, "UniformOutput", false)); |
| 6770 | 65 endif |
| 66 ndims = size (subs, 2); | |
| 67 | |
| 68 if (nargin < 5 || isempty (fillval)) | |
| 69 fillval = 0; | |
| 70 endif | |
| 71 | |
| 72 if (nargin < 6 || isempty (isspar)) | |
| 73 isspar = false; | |
| 74 endif | |
| 7186 | 75 |
| 76 if (isspar && ndims > 2) | |
| 77 error ("accumarray: sparse matrices limited to 2 dimensions"); | |
| 6770 | 78 endif |
| 79 | |
|
8934
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
80 if (nargin < 4 || isempty (fun)) |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
81 fun = @sum; |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
82 ## This is the fast summation case. Unlike the general case, |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
83 ## this case will be handled using an O(N) algorithm. |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
84 |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
85 if (isspar && fillval == 0) |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
86 ## The "sparse" function can handle this case. |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
87 |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
88 if ((nargin < 3 || isempty (sz))) |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
89 A = sparse (subs(:,1), subs(:,2), val); |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
90 elseif (length (sz) == 2) |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
91 A = sparse (subs(:,1), subs(:,2), val, sz(1), sz(2)); |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
92 else |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
93 error ("accumarray: dimensions mismatch") |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
94 endif |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
95 else |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
96 ## This case is handled by an internal function. |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
97 |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
98 if (ndims > 1) |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
99 if ((nargin < 3 || isempty (sz))) |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
100 sz = max (subs); |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
101 elseif (ndims != length (sz)) |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
102 error ("accumarray: dimensions mismatch") |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
103 elseif (any (max (subs) > sz)) |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
104 error ("accumarray: index out of range") |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
105 endif |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
106 |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
107 ## Convert multidimensional subscripts. |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
108 subs = sub2ind (sz, mat2cell (subs, rows (subs), ones (1, ndims)){:}); |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
109 elseif (nargin < 3) |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
110 ## In case of linear indexing, the fast built-in accumulator |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
111 ## will determine the extent for us. |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
112 sz = []; |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
113 endif |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
114 |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
115 ## Call the built-in accumulator. |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
116 if (isempty (sz)) |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
117 A = __accumarray_sum__ (subs, val); |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
118 else |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
119 A = __accumarray_sum__ (subs, val, prod (sz)); |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
120 ## set proper shape. |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
121 A = reshape (A, sz); |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
122 endif |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
123 |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
124 ## we fill in nonzero fill value. |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
125 if (fillval != 0) |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
126 mask = true (size (A)); |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
127 mask(subs) = false; |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
128 A(mask) = fillval; |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
129 endif |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
130 endif |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
131 |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
132 return |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
133 endif |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
134 |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
135 if (nargin < 3 || isempty (sz)) |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
136 sz = max (subs); |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
137 if (isscalar(sz)) |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
138 sz = [sz, 1]; |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
139 endif |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
140 elseif (length (sz) != ndims |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
141 && (ndims != 1 || length (sz) != 2 || sz(2) != 1)) |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
142 error ("accumarray: inconsistent dimensions"); |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
143 endif |
|
c2099a4d12ea
partially optimize accumarray
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
144 |
| 6770 | 145 [subs, idx] = sortrows (subs); |
|
8820
89b95972e178
fix previously introduced problem in octave_sort, improve design
Jaroslav Hajek <highegg@gmail.com>
parents:
8507
diff
changeset
|
146 |
| 6770 | 147 if (isscalar (val)) |
| 148 val = val * ones (size (idx)); | |
| 149 else | |
| 150 val = val(idx); | |
| 151 endif | |
| 7186 | 152 cidx = find ([true; (sum (abs (diff (subs)), 2) != 0)]); |
| 6770 | 153 idx = cell (1, ndims); |
| 154 for i = 1:ndims | |
| 155 idx{i} = subs (cidx, i); | |
| 156 endfor | |
| 157 x = cellfun (fun, mat2cell (val(:), diff ([cidx; length(val) + 1]))); | |
| 158 if (isspar && fillval == 0) | |
| 159 A = sparse (idx{1}, idx{2}, x, sz(1), sz(2)); | |
| 160 else | |
| 161 if (iscell (x)) | |
| 162 ## Why did matlab choose to reverse the order of the elements | |
| 8507 | 163 x = cellfun (@(x) flipud (x(:)), x, "UniformOutput", false); |
| 6770 | 164 A = cell (sz); |
| 165 elseif (fillval == 0) | |
| 7186 | 166 A = zeros (sz, class (x)); |
| 6770 | 167 else |
| 168 A = fillval .* ones (sz); | |
| 169 endif | |
| 7186 | 170 A(sub2ind (sz, idx{:})) = x; |
| 6770 | 171 endif |
| 172 endfunction | |
| 173 | |
| 174 %!error (accumarray (1:5)) | |
| 175 %!error (accumarray ([1,2,3],1:2)) | |
| 176 %!assert (accumarray ([1;2;4;2;4],101:105), [101;206;0;208]) | |
| 177 %!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])) | |
| 178 %!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]))) | |
| 179 %!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])) | |
| 180 %!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]) | |
| 181 %!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)) | |
| 182 %!assert (accumarray ([1 1; 2 1; 2 3; 2 1; 2 3],1,[2,4]), [1,0,0,0;2,0,2,0]) | |
| 183 %!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]) | |
| 184 %!test | |
| 185 %! A = accumarray ([1 1; 2 1; 2 3; 2 1; 2 3],101:105,[2,4],@(x){x}); | |
| 186 %! assert (A{2},[104;102]) |