|
19007
|
1 /* |
|
|
2 |
|
|
3 Copyright (C) 2014 David Spies |
|
|
4 |
|
|
5 This file is part of Octave. |
|
|
6 |
|
|
7 Octave is free software; you can redistribute it and/or modify it |
|
|
8 under the terms of the GNU General Public License as published by the |
|
|
9 Free Software Foundation; either version 3 of the License, or (at your |
|
|
10 option) any later version. |
|
|
11 |
|
|
12 Octave is distributed in the hope that it will be useful, but WITHOUT |
|
|
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
|
|
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
|
|
15 for more details. |
|
|
16 |
|
|
17 You should have received a copy of the GNU General Public License |
|
|
18 along with Octave; see the file COPYING. If not, see |
|
|
19 <http://www.gnu.org/licenses/>. |
|
|
20 |
|
|
21 */ |
|
|
22 |
|
|
23 // A template for functions which accumulate all the values along |
|
|
24 // the rows or columns of a matrix into a single row or column vector |
|
|
25 // (more generally removing one dimension of a many-dimensional matrix) |
|
|
26 // Functions which should call this one: |
|
|
27 // min, max, sum, product, any, all |
|
|
28 // |
|
|
29 // TODO: Allow possible second return value for min and max |
|
|
30 #include "dim-vector.h" |
|
|
31 #include "Array.h" |
|
|
32 |
|
|
33 // TODO Support any/all with short-circuiting |
|
|
34 enum fun_type |
|
|
35 { |
|
|
36 min, max, sum, prod |
|
|
37 }; |
|
|
38 |
|
|
39 const bool MIN = false; |
|
|
40 const bool MAX = true; |
|
|
41 |
|
|
42 const int SUM = 0; |
|
|
43 const int PROD = 1; |
|
|
44 |
|
|
45 const int COL_ACCUM = 0; |
|
|
46 const int ROW_ACCUM = 1; |
|
|
47 |
|
|
48 // Pairs a value together with its index along the |
|
|
49 // accumulated dimension |
|
|
50 template<typename T> |
|
|
51 struct Elem |
|
|
52 { |
|
|
53 const octave_idx_type accum_index; |
|
|
54 const T value; |
|
|
55 Elem (octave_idx_type accum_index, T value) : |
|
|
56 accum_index (accum_index), value (value) { } |
|
|
57 }; |
|
|
58 |
|
|
59 //Chooses < or > based on template argument |
|
|
60 //(but only uses < operation) |
|
|
61 template<typename E, bool MinOrMax> |
|
|
62 bool |
|
|
63 exceeds (E v1, E v2) { return (MinOrMax == MIN) ? v1 < v2 : v2 < v1; } |
|
|
64 |
|
|
65 template<typename E, bool MinOrMax> |
|
|
66 struct FullMMAccum |
|
|
67 { |
|
|
68 public: |
|
|
69 |
|
|
70 Array<Elem<E> > res; |
|
|
71 Array<double> ind; |
|
|
72 |
|
|
73 FullMMAccum (const dim_vector& dims) : res (dims), ind (dims, 0) { } |
|
|
74 |
|
|
75 void |
|
|
76 add (octave_idx_type to_index, Elem<E> rhs) |
|
|
77 { |
|
|
78 if (ind (to_index) == 0 || exceeds<E, MinOrMax> (rhs.value, res (to_index))) |
|
|
79 { |
|
|
80 ind (to_index) = rhs.index + 1; |
|
|
81 res (to_index) = rhs.value; |
|
|
82 } |
|
|
83 } |
|
|
84 |
|
|
85 octave_value_list |
|
|
86 post_process (void) |
|
|
87 { |
|
|
88 return octave_value_list (res, ind); |
|
|
89 } |
|
|
90 }; |
|
|
91 |
|
|
92 // TODO Separate Sum and Prod |
|
|
93 // for short-circuiting Prod |
|
|
94 template<typename E, int SorP> |
|
|
95 struct FullSPAccum |
|
|
96 { |
|
|
97 |
|
|
98 Array<E> res; |
|
|
99 |
|
|
100 FullSPAccum (const dim_vector& dims) : res (dims, SorP) { } |
|
|
101 |
|
|
102 void |
|
|
103 add (octave_idx_type to_index, Elem<E> rhs) |
|
|
104 { |
|
|
105 if (SorP == PROD) |
|
|
106 res (to_index) *= rhs.value; |
|
|
107 else |
|
|
108 res (to_index) += rhs.value; |
|
|
109 } |
|
|
110 |
|
|
111 octave_value_list |
|
|
112 post_process () |
|
|
113 { |
|
|
114 return octave_value_list (res); |
|
|
115 } |
|
|
116 }; |
|
|
117 |
|
|
118 template<typename E, bool MinOrMax> |
|
|
119 struct SparseMMAccum |
|
|
120 { |
|
|
121 private: |
|
|
122 FullMMAccum<E, MinOrMax> my_accum; |
|
|
123 public: |
|
|
124 SparseMMAccum (dim_vector dims) : my_accum (dims) { } |
|
|
125 |
|
|
126 void |
|
|
127 add (octave_idx_type ind, Elem<E> value) |
|
|
128 { |
|
|
129 my_accum.add (ind, value); |
|
|
130 } |
|
|
131 |
|
|
132 octave_value_list |
|
|
133 post_process (void) |
|
|
134 { |
|
|
135 return octave_value_list (Sparse<E> (my_accum.res), my_accum.ind); |
|
|
136 } |
|
|
137 }; |
|
|
138 |
|
|
139 // TODO Separate Sum and Prod |
|
|
140 // for short-circuiting Prod |
|
|
141 template<typename E, int SorP> |
|
|
142 struct SparseSPAccum |
|
|
143 { |
|
|
144 private: |
|
|
145 FullSPAccum<E, SorP> my_accum; |
|
|
146 public: |
|
|
147 SparseSPAccum (dim_vector dims) : my_accum (dims) { } |
|
|
148 |
|
|
149 void |
|
|
150 add (octave_idx_type ind, Elem<E> value) |
|
|
151 { |
|
|
152 my_accum.add (ind, value); |
|
|
153 } |
|
|
154 |
|
|
155 octave_value_list |
|
|
156 post_process (void) |
|
|
157 { |
|
|
158 return octave_value_list (Sparse<E> (my_accum.res)); |
|
|
159 } |
|
|
160 }; |
|
|
161 |
|
|
162 template<typename E> |
|
|
163 octave_value_list |
|
|
164 sparseColumn (octave_idx_type height, |
|
|
165 std::vector<std::pair<octave_idx_type, E> > vec) |
|
|
166 { |
|
|
167 Sparse<E> res = Sparse<E> (height, 1, vec.size ()); |
|
|
168 res.xcidx (0) = 0; |
|
|
169 res.xcidx (1) = vec.size (); |
|
|
170 for (size_t i = 0; i < vec.size; ++i) |
|
|
171 { |
|
|
172 res.xridx (i) = vec[i].first; |
|
|
173 res.xdata (i) = vec[i].second; |
|
|
174 } |
|
|
175 return octave_value_list (res); |
|
|
176 } |
|
|
177 |
|
|
178 template<typename E, bool MinOrMax> |
|
|
179 struct TallMMAccum |
|
|
180 { |
|
|
181 private: |
|
|
182 const octave_idx_type height; |
|
|
183 std::vector<std::pair<octave_idx_type, E> > res_vec; |
|
|
184 |
|
|
185 public: |
|
|
186 TallMMAccum (octave_idx_type height) : height (height) { } |
|
|
187 |
|
|
188 void |
|
|
189 add_term (octave_idx_type row, Elem<E> value) |
|
|
190 { |
|
|
191 res_vec.push_back (std::pair (row, value.value)); |
|
|
192 } |
|
|
193 |
|
|
194 void |
|
|
195 join_term (Elem<E> value) |
|
|
196 { |
|
|
197 if (exceeds<E, MinOrMax> (value.value, res_vec.back ().second)) |
|
|
198 { |
|
|
199 res_vec.back ().second = value.value; |
|
|
200 } |
|
|
201 } |
|
|
202 |
|
|
203 octave_value_list |
|
|
204 post_process (void) |
|
|
205 { |
|
|
206 return sparseColumn (height, res_vec); |
|
|
207 } |
|
|
208 }; |
|
|
209 |
|
|
210 template<typename E, int SorP> |
|
|
211 struct TallSPAccum |
|
|
212 { |
|
|
213 private: |
|
|
214 const octave_idx_type height; |
|
|
215 std::vector<std::pair<octave_idx_type, E> > res_vec; |
|
|
216 |
|
|
217 public: |
|
|
218 TallSPAccum (octave_idx_type height) : height (height) { } |
|
|
219 |
|
|
220 void |
|
|
221 add_term (octave_idx_type row, Elem<E> value) |
|
|
222 { |
|
|
223 res_vec.push_back (std::pair (row, value.value)); |
|
|
224 } |
|
|
225 |
|
|
226 void |
|
|
227 join_term (Elem<E> value) |
|
|
228 { |
|
|
229 res_vec.back ().second += value.value; |
|
|
230 } |
|
|
231 |
|
|
232 octave_value_list |
|
|
233 post_process () |
|
|
234 { |
|
|
235 return sparseColumn (height, res_vec); |
|
|
236 } |
|
|
237 }; |
|
|
238 |
|
|
239 // M is the type of matrix. Right now I'm only handling |
|
|
240 // Array or Sparse. TODO What other types must I handle? |
|
|
241 template<typename M> |
|
|
242 octave_value_list |
|
|
243 accumulate (const M& matrix, int dim, int nargout, fun_type t) |
|
|
244 { |
|
|
245 if (dim < 0) |
|
|
246 { |
|
|
247 // TODO Raise exception: how? |
|
|
248 assert(false); |
|
|
249 } |
|
|
250 else if (dim >= matrix.ndims ()) |
|
|
251 return octave_value_list (matrix); |
|
|
252 // TODO Add is_sparse to SparseMatrix |
|
|
253 // Is there another way to check? |
|
|
254 else if (M::is_sparse) |
|
|
255 return select_sparse_accumulate<M::element_type> (matrix, dim, t, nargout); |
|
|
256 else |
|
|
257 return select_full_accumulate<M::element_type> (matrix, dim, t); |
|
|
258 } |
|
|
259 |
|
|
260 template<typename E> |
|
|
261 octave_value_list |
|
|
262 select_full_accumulate (const Array<E>& matrix, int dim, fun_type t) |
|
|
263 { |
|
|
264 // TODO Deal with h = 0 and w = 0 |
|
|
265 switch (t) |
|
|
266 { |
|
|
267 case min: |
|
|
268 return full_accumulate<E, FullMMAccum<E, MIN> > (matrix, dim); |
|
|
269 case max: |
|
|
270 return full_accumulate<E, FullMMAccum<E, MAX> > (matrix, dim); |
|
|
271 case sum: |
|
|
272 return full_accumulate<E, FullSPAccum<E, SUM> > (matrix, dim); |
|
|
273 case prod: |
|
|
274 return full_accumulate<E, FullSPAccum<E, PROD> > (matrix, dim); |
|
|
275 default: |
|
|
276 assert(false); |
|
|
277 // TODO Raise exception instead? |
|
|
278 } |
|
|
279 } |
|
|
280 |
|
|
281 // Change dim to a template parameter |
|
|
282 template<typename E> |
|
|
283 octave_value_list |
|
|
284 select_sparse_accumulate (const Sparse<E>& matrix, int dim, fun_type t, |
|
|
285 int nargout) |
|
|
286 { |
|
|
287 switch (dim) |
|
|
288 { |
|
|
289 case ROW_ACCUM: |
|
|
290 return select_sparse_accumulate<E, ROW_ACCUM> (matrix, t, nargout); |
|
|
291 case COL_ACCUM: |
|
|
292 return select_sparse_accumulate<E, COL_ACCUM> (matrix, t, nargout); |
|
|
293 default: |
|
|
294 // This was already checked by accumulate |
|
|
295 assert(false); |
|
|
296 } |
|
|
297 } |
|
|
298 |
|
|
299 // Accumulating over the columns of a sparse matrix is simple. |
|
|
300 // The rows are a little more complicated |
|
|
301 // There are two possible algorithms for accumulating over |
|
|
302 // the rows of a sparse matrix |
|
|
303 // the simple one runs in O(h+nnz) and just |
|
|
304 // constructs a full result vector (which can |
|
|
305 // be converted to a sparse vector). |
|
|
306 // The tricky one uses a heap to construct a sparse |
|
|
307 // result vector directly and runs in O(w+nnz*log(w)). |
|
|
308 // We need to look at matrix's dimensions to see |
|
|
309 // which one will work faster. If the matrix |
|
|
310 // is tall and sparse enough, we use the second one. |
|
|
311 |
|
|
312 template<typename E, int dim> |
|
|
313 octave_value_list |
|
|
314 select_sparse_accumulate (const Sparse<E>& matrix, fun_type t, int nargout) |
|
|
315 { |
|
|
316 octave_idx_type w = matrix.cols (); |
|
|
317 octave_idx_type h = matrix.rows (); |
|
|
318 octave_idx_type nnz = matrix.nnz (); |
|
|
319 //TODO Deal with h = 0 and w = 0 |
|
|
320 // For perfect optimality there should probably be some constants |
|
|
321 // before each of these terms (determined by extensive profiling). |
|
|
322 // But for big-O optimality, just making this comparison is good |
|
|
323 // enough. |
|
|
324 // Using float log because being precise isn't important. |
|
|
325 if (nargout == 1 && dim == ROW_ACCUM && w + nnz * log2f ((float) w) < h + nnz) |
|
|
326 { |
|
|
327 switch (t) |
|
|
328 { |
|
|
329 case min: |
|
|
330 return sparse_tall_row_accum<E, TallMMAccum<E, MIN> > (matrix); |
|
|
331 case max: |
|
|
332 return sparse_tall_row_accum<E, TallMMAccum<E, MAX> > (matrix); |
|
|
333 case sum: |
|
|
334 return sparse_tall_row_accum<E, TallSPAccum<E, SUM> > (matrix); |
|
|
335 case prod: |
|
|
336 return sparse_tall_row_accum<E, TallSPAccum<E, PROD> > (matrix); |
|
|
337 } |
|
|
338 } |
|
|
339 else |
|
|
340 { |
|
|
341 switch (t) |
|
|
342 { |
|
|
343 case min: |
|
|
344 return sparse_normal_accum<E, SparseMMAccum<E, MIN>, dim> (matrix); |
|
|
345 case max: |
|
|
346 return sparse_normal_accum<E, SparseMMAccum<E, MAX>, dim> (matrix); |
|
|
347 case sum: |
|
|
348 return sparse_normal_accum<E, SparseSPAccum<E, SUM>, dim> (matrix); |
|
|
349 case prod: |
|
|
350 return sparse_normal_accum<E, SparseSPAccum<E, PROD>, dim> (matrix); |
|
|
351 } |
|
|
352 } |
|
|
353 } |
|
|
354 |
|
|
355 // excluded_indexer takes a dim_vector v and an "excluded dimension" d |
|
|
356 // and iterates over the indices 0 to n (where n = v.numel()) while |
|
|
357 // simultaneously maintaining an iterator over (v excluding d) meaning |
|
|
358 // the index resulting from removing the dth element from v. |
|
|
359 // ex. let v = [3,4,5] and d = 1 |
|
|
360 // then we can use an excluded_indexer to iterate over the pairs |
|
|
361 // from_idx to_idx |
|
|
362 // |
|
|
363 // 0 0 (= [0,0,0], [0,0]) |
|
|
364 // 1 1 (= [1,0,0], [1,0]) |
|
|
365 // 2 2 (= [2,0,0], [2,0]) |
|
|
366 // 3 0 (= [0,1,0], [0,0]) |
|
|
367 // 4 1 (= [1,1,0], [1,0]) |
|
|
368 // 5 2 (= [2,1,0], [2,0]) |
|
|
369 // 6 0 (= [0,2,0], [0,0]) |
|
|
370 // 7 1 (= [1,2,0], [1,0]) |
|
|
371 // 8 2 (= [2,2,0], [2,0]) |
|
|
372 // 9 0 (= [0,3,0], [0,0]) |
|
|
373 // 10 1 (= [1,3,0], [1,0]) |
|
|
374 // 11 2 (= [2,3,0], [2,0]) |
|
|
375 // 12 3 (= [0,0,1], [0,1]) |
|
|
376 // 13 4 (= [1,0,1], [1,1]) |
|
|
377 // 14 5 (= [2,0,1], [2,1]) |
|
|
378 // 15 3 (= [0,1,1], [0,1]) |
|
|
379 // ...etc... |
|
|
380 // 58 13 (= [1,3,4], [1,3]) |
|
|
381 // 59 14 (= [2,3,4], [2,4]) |
|
|
382 // |
|
|
383 // This is done extremely efficiently using only incrementing and subtraction |
|
|
384 // (no multiplication,division modulo etc.) and is the most spatial-locality |
|
|
385 // friendly way to iterate over the two matrices. |
|
|
386 |
|
|
387 struct excluded_indexer |
|
|
388 { |
|
|
389 private: |
|
|
390 const octave_idx_type before_size; |
|
|
391 const octave_idx_type between_size; |
|
|
392 |
|
|
393 octave_idx_type from_idx; |
|
|
394 octave_idx_type to_idx; |
|
|
395 |
|
|
396 octave_idx_type before_ctr; |
|
|
397 octave_idx_type between_ctr; |
|
|
398 |
|
|
399 public: |
|
|
400 |
|
|
401 excluded_indexer (const dim_vector& sizes, int dim) : |
|
|
402 before_size (product_prefix (sizes, dim)), between_size (sizes (dim)), |
|
|
403 from_idx (0), to_idx (0), before_ctr (0), between_ctr (0) { } |
|
|
404 |
|
|
405 // Returns the index for the input (larger) matrix |
|
|
406 octave_idx_type |
|
|
407 get_from_idx (void) const |
|
|
408 { |
|
|
409 return from_idx; |
|
|
410 } |
|
|
411 |
|
|
412 // Returns the index for the output (smaller) matrix |
|
|
413 octave_idx_type |
|
|
414 get_to_idx (void) const |
|
|
415 { |
|
|
416 return to_idx; |
|
|
417 } |
|
|
418 |
|
|
419 octave_idx_type |
|
|
420 get_accum_idx (void) const |
|
|
421 { |
|
|
422 return between_ctr; |
|
|
423 } |
|
|
424 |
|
|
425 // Increments both indices and then subtracts from |
|
|
426 // the output index if necessary |
|
|
427 void |
|
|
428 operator++ (void) |
|
|
429 { |
|
|
430 ++from_idx; |
|
|
431 ++to_idx; |
|
|
432 ++before_ctr; |
|
|
433 if (before_ctr == before_size) |
|
|
434 { |
|
|
435 before_ctr = 0; |
|
|
436 ++between_ctr; |
|
|
437 if (between_ctr == between_size) |
|
|
438 between_ctr = 0; |
|
|
439 else |
|
|
440 to_idx -= before_size; |
|
|
441 } |
|
|
442 } |
|
|
443 |
|
|
444 // Compares the input index with from_size |
|
|
445 bool |
|
|
446 operator< (octave_idx_type from_size) const |
|
|
447 { |
|
|
448 return from_idx < from_size; |
|
|
449 } |
|
|
450 }; |
|
|
451 |
|
|
452 // Accumulate method for full matrices |
|
|
453 // Accumulator takes a dimension in its constructor. It supports methods |
|
|
454 // void add(octave_idx_type, Elem<E>) and octave_value_list post_process() |
|
|
455 template<typename E, typename Accumulator> |
|
|
456 octave_value_list |
|
|
457 full_accumulate (const Array<E>& matrix, int dim) |
|
|
458 { |
|
|
459 const dim_vector sizes = matrix.get_dims (); |
|
|
460 dim_vector res_dims = sizes; |
|
|
461 res_dims (dim) = 1; |
|
|
462 // TODO Can I do this? Is it safe to mutate a dim_vector |
|
|
463 // Does (res_dims = sizes) make a copy or copy the pointer? |
|
|
464 Accumulator res (res_dims); |
|
|
465 octave_idx_type from_numel = sizes.numel (); |
|
|
466 octave_idx_type to_numel = res_dims.numel (); |
|
|
467 for (excluded_indexer i (sizes, dim); i < from_numel; ++i) |
|
|
468 { |
|
|
469 res.add (i.get_to_idx (), |
|
|
470 Elem<E> (i.get_accum_idx (), matrix.xelem (i.get_from_idx ()))); |
|
|
471 } |
|
|
472 return res.post_process (); |
|
|
473 } |
|
|
474 |
|
|
475 // Takes the product of the first dim dimensions of sizes |
|
|
476 octave_idx_type |
|
|
477 product_prefix (const dim_vector& sizes, int dim) |
|
|
478 { |
|
|
479 assert(dim <= sizes.ndims ()); |
|
|
480 octave_idx_type res = 1; |
|
|
481 for (int i = 0; i < dim; ++i) |
|
|
482 { |
|
|
483 res *= sizes (i); |
|
|
484 } |
|
|
485 return res; |
|
|
486 } |
|
|
487 |
|
|
488 template<typename E, typename Accumulator, int accum_dir> |
|
|
489 octave_value_list |
|
|
490 sparse_normal_accum (const Sparse<E>& matrix) |
|
|
491 { |
|
|
492 assert(accum_dir == 1 || accum_dir == 0); |
|
|
493 dim_vector accum_dims = matrix.dims (); |
|
|
494 octave_idx_type idx_size = accum_dims (1 - accum_dir); |
|
|
495 octave_idx_type accum_size = accum_dims (accum_dir); |
|
|
496 accum_dims (accum_dir) = 1; |
|
|
497 Accumulator res (accum_dims); |
|
|
498 std::vector<octave_idx_type> known_mex (idx_size, 0); |
|
|
499 for (octave_idx_type col = 0; col < matrix.cols (); ++col) |
|
|
500 { |
|
|
501 for (octave_idx_type i = matrix.cidx (col); i < matrix.cidx (col + 1); |
|
|
502 ++i) |
|
|
503 { |
|
|
504 octave_idx_type res_idx; |
|
|
505 octave_idx_type accum_idx; |
|
|
506 if (accum_dir == COL_ACCUM) |
|
|
507 { |
|
|
508 res_idx = col; |
|
|
509 accum_idx = matrix.ridx (i); |
|
|
510 } |
|
|
511 else if (accum_dir == ROW_ACCUM) |
|
|
512 { |
|
|
513 res_idx = matrix.ridx (i); |
|
|
514 accum_idx = col; |
|
|
515 } |
|
|
516 else |
|
|
517 { |
|
|
518 assert(false); |
|
|
519 } |
|
|
520 if (known_mex[res_idx] != (octave_idx_type) (-1)) |
|
|
521 { |
|
|
522 if (accum_idx == known_mex[res_idx]) |
|
|
523 ++known_mex[res_idx]; |
|
|
524 else |
|
|
525 { |
|
|
526 assert(accum_idx > known_mex[res_idx]); |
|
|
527 res.add (res_idx, Elem<E> (known_mex[res_idx], 0)); |
|
|
528 known_mex[res_idx] = (octave_idx_type) (-1); |
|
|
529 } |
|
|
530 } |
|
|
531 res.add (res_idx, Elem<E> (accum_idx, matrix.data[i])); |
|
|
532 } |
|
|
533 } |
|
|
534 for (octave_idx_type i = 0; i < idx_size; ++i) |
|
|
535 { |
|
|
536 if (known_mex[i] != (octave_idx_type) (-1) && known_mex[i] < accum_size) |
|
|
537 { |
|
|
538 res.add (i, Elem<E> (known_mex[i], 0)); |
|
|
539 } |
|
|
540 } |
|
|
541 return res.post_process (); |
|
|
542 } |
|
|
543 |
|
|
544 // Knows everything about an index-value pair |
|
|
545 // in a sparse matrix (row, col, nnz_idx, value) |
|
|
546 template<typename E> |
|
|
547 struct SparseElem |
|
|
548 { |
|
|
549 const octave_idx_type row; |
|
|
550 const octave_idx_type col; |
|
|
551 const E value; |
|
|
552 const octave_idx_type idx; |
|
|
553 |
|
|
554 SparseElem (octave_idx_type row, octave_idx_type col, E value, |
|
|
555 octave_idx_type idx) : |
|
|
556 row (row), col (col), value (value), idx (idx) { } |
|
|
557 |
|
|
558 bool |
|
|
559 operator< (const SparseElem<E> rhs) const |
|
|
560 { |
|
|
561 return row < rhs.row || (row == rhs.row && col < rhs.col); |
|
|
562 } |
|
|
563 }; |
|
|
564 |
|
|
565 template<typename E, typename TallAccumulator> |
|
|
566 octave_value_list |
|
|
567 tall_sparse_row_accum (const Sparse<E>& matrix) |
|
|
568 { |
|
|
569 TallAccumulator res; |
|
|
570 |
|
|
571 std::priority_queue<SparseElem<E> > nextItem; |
|
|
572 for (octave_idx_type col = 0; col < matrix.cols (); ++col) |
|
|
573 { |
|
|
574 octave_idx_type idx = matrix.cidx[col]; |
|
|
575 nextItem.push ( |
|
|
576 SparseElem<E> (matrix.ridx[idx], col, matrix.data[idx], idx)); |
|
|
577 } |
|
|
578 |
|
|
579 octave_idx_type current_row = -1; |
|
|
580 |
|
|
581 // mex is short for minimum excluded value. |
|
|
582 // (meaning the smallest nonnegative integer which is not |
|
|
583 // in a set. accum_mex is the mex of the columns containing |
|
|
584 // nonzero elements in the current row. |
|
|
585 // It's uncommon to find a use for the function |
|
|
586 // outside of combinatoric game theory. In this case, |
|
|
587 // a row contains a zero if accum_mex < width. |
|
|
588 octave_idx_type accum_mex; |
|
|
589 |
|
|
590 // It's either macro or copy-paste. I chose the lesser |
|
|
591 // of two evils. If you don't like it, learn Lisp. |
|
|
592 // *** |
|
|
593 #define finish_row \ |
|
|
594 if (current_row >= 0 \ |
|
|
595 && accum_mex >= 0 \ |
|
|
596 && accum_mex < matrix.cols ()) \ |
|
|
597 res.join_term (Elem<E> (accum_mex, 0)); \ |
|
|
598 |
|
|
599 while (!nextItem.empty ()) |
|
|
600 { |
|
|
601 SparseElem<E> item = nextItem.front (); |
|
|
602 nextItem.pop (); |
|
|
603 if (item.row != current_row) |
|
|
604 { |
|
|
605 // macro; see *** above |
|
|
606 finish_row; |
|
|
607 |
|
|
608 current_row = item.row; |
|
|
609 if (item.col > 0) |
|
|
610 { |
|
|
611 res.add_term (item.row, Elem<E> (0, 0)); |
|
|
612 accum_mex = -1; |
|
|
613 } |
|
|
614 else |
|
|
615 accum_mex = 1; |
|
|
616 res.add_term (item.row, Elem<E> (item.col, item.value)); |
|
|
617 } |
|
|
618 else |
|
|
619 { |
|
|
620 if (accum_mex != (octave_idx_type) (-1)) |
|
|
621 { |
|
|
622 if (item.col == accum_mex) |
|
|
623 ++accum_mex; |
|
|
624 else |
|
|
625 { |
|
|
626 res.join_term (Elem<E> (accum_mex, 0)); |
|
|
627 accum_mex = -1; |
|
|
628 } |
|
|
629 } |
|
|
630 res.join_term (Elem<E> (item.col, item.value)); |
|
|
631 } |
|
|
632 |
|
|
633 octave_idx_type idx = item.idx + 1; |
|
|
634 octave_idx_type col = item.col; |
|
|
635 if (idx < matrix.cidx[col + 1]) |
|
|
636 { |
|
|
637 nextItem.push ( |
|
|
638 SparseElem<E> (matrix.ridx[idx], col, matrix.data[idx], idx)); |
|
|
639 } |
|
|
640 } |
|
|
641 |
|
|
642 // macro; see *** above |
|
|
643 finish_row; |
|
|
644 |
|
|
645 #undef finish_row |
|
|
646 |
|
|
647 return res.post_process (); |
|
|
648 } |
