|
5000
|
1 /* |
|
2928
|
2 |
|
|
3 Copyright (C) 1996, 1997 John W. Eaton |
|
4850
|
4 Copyright (C) 2004 David Bateman |
|
2928
|
5 |
|
|
6 This file is part of Octave. |
|
|
7 |
|
|
8 Octave is free software; you can redistribute it and/or modify it |
|
|
9 under the terms of the GNU General Public License as published by the |
|
|
10 Free Software Foundation; either version 2, or (at your option) any |
|
|
11 later version. |
|
|
12 |
|
|
13 Octave is distributed in the hope that it will be useful, but WITHOUT |
|
|
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
|
|
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
|
|
16 for more details. |
|
|
17 |
|
|
18 You should have received a copy of the GNU General Public License |
|
|
19 along with Octave; see the file COPYING. If not, write to the Free |
|
5307
|
20 Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA |
|
|
21 02110-1301, USA. |
|
2928
|
22 |
|
|
23 */ |
|
|
24 |
|
|
25 #ifdef HAVE_CONFIG_H |
|
|
26 #include <config.h> |
|
|
27 #endif |
|
|
28 |
|
5164
|
29 #include <vector> |
|
|
30 |
|
3826
|
31 #include "lo-mappers.h" |
|
4153
|
32 #include "quit.h" |
|
3826
|
33 |
|
2928
|
34 #include "defun-dld.h" |
|
|
35 #include "error.h" |
|
|
36 #include "gripes.h" |
|
|
37 #include "oct-obj.h" |
|
4850
|
38 #include "lo-ieee.h" |
|
|
39 #include "data-conv.h" |
|
|
40 #include "ov-cx-mat.h" |
|
4996
|
41 #include "ov-cell.h" |
|
4850
|
42 #include "oct-sort.cc" |
|
2928
|
43 |
|
4997
|
44 enum sortmode { UNDEFINED, ASCENDING, DESCENDING }; |
|
4996
|
45 |
|
|
46 template <class T> |
|
|
47 class |
|
|
48 vec_index |
|
|
49 { |
|
4997
|
50 public: |
|
4996
|
51 T vec; |
|
5275
|
52 octave_idx_type indx; |
|
4996
|
53 }; |
|
|
54 |
|
|
55 template <class T> |
|
5009
|
56 bool |
|
|
57 ascending_compare (T a, T b) |
|
|
58 { |
|
|
59 return (a < b); |
|
|
60 } |
|
|
61 |
|
|
62 template <class T> |
|
|
63 bool |
|
|
64 descending_compare (T a, T b) |
|
|
65 { |
|
|
66 return (a > b); |
|
|
67 } |
|
|
68 |
|
|
69 template <class T> |
|
|
70 bool |
|
|
71 ascending_compare (vec_index<T> *a, vec_index<T> *b) |
|
|
72 { |
|
|
73 return (a->vec < b->vec); |
|
|
74 } |
|
|
75 |
|
|
76 template <class T> |
|
|
77 bool |
|
|
78 descending_compare (vec_index<T> *a, vec_index<T> *b) |
|
|
79 { |
|
|
80 return (a->vec > b->vec); |
|
|
81 } |
|
|
82 |
|
|
83 template <class T> |
|
4998
|
84 static octave_value |
|
4996
|
85 mx_sort (ArrayN<T> &m, int dim, sortmode mode = UNDEFINED) |
|
|
86 { |
|
4998
|
87 octave_value retval; |
|
4996
|
88 |
|
5366
|
89 dim_vector dv = m.dims (); |
|
|
90 |
|
4996
|
91 if (m.length () < 1) |
|
5366
|
92 return ArrayN<T> (dv); |
|
4996
|
93 |
|
5275
|
94 octave_idx_type ns = dv(dim); |
|
|
95 octave_idx_type iter = dv.numel () / ns; |
|
|
96 octave_idx_type stride = 1; |
|
|
97 for (int i = 0; i < dim; i++) |
|
4996
|
98 stride *= dv(i); |
|
|
99 |
|
|
100 T *v = m.fortran_vec (); |
|
|
101 octave_sort<T> sort; |
|
|
102 |
|
|
103 if (mode == ASCENDING) |
|
|
104 sort.set_compare (ascending_compare); |
|
|
105 else if (mode == DESCENDING) |
|
|
106 sort.set_compare (descending_compare); |
|
|
107 |
|
|
108 if (stride == 1) |
|
4997
|
109 { |
|
5275
|
110 for (octave_idx_type j = 0; j < iter; j++) |
|
4997
|
111 { |
|
|
112 sort.sort (v, ns); |
|
|
113 v += ns; |
|
|
114 } |
|
|
115 } |
|
4996
|
116 else |
|
|
117 { |
|
|
118 OCTAVE_LOCAL_BUFFER (T, vi, ns); |
|
5275
|
119 for (octave_idx_type j = 0; j < iter; j++) |
|
4996
|
120 { |
|
5275
|
121 octave_idx_type offset = j; |
|
|
122 octave_idx_type offset2 = 0; |
|
4996
|
123 while (offset >= stride) |
|
|
124 { |
|
|
125 offset -= stride; |
|
|
126 offset2++; |
|
|
127 } |
|
|
128 offset += offset2 * stride * ns; |
|
|
129 |
|
5275
|
130 for (octave_idx_type i = 0; i < ns; i++) |
|
4996
|
131 vi[i] = v[i*stride + offset]; |
|
|
132 |
|
|
133 sort.sort (vi, ns); |
|
|
134 |
|
5275
|
135 for (octave_idx_type i = 0; i < ns; i++) |
|
4996
|
136 v[i*stride + offset] = vi[i]; |
|
|
137 } |
|
|
138 } |
|
|
139 |
|
4998
|
140 retval = m; |
|
4997
|
141 |
|
4996
|
142 return retval; |
|
|
143 } |
|
|
144 |
|
|
145 template <class T> |
|
|
146 static octave_value_list |
|
|
147 mx_sort_indexed (ArrayN<T> &m, int dim, sortmode mode = UNDEFINED) |
|
|
148 { |
|
|
149 octave_value_list retval; |
|
|
150 |
|
5366
|
151 dim_vector dv = m.dims (); |
|
|
152 |
|
4996
|
153 if (m.length () < 1) |
|
5360
|
154 { |
|
5366
|
155 retval(1) = NDArray (dv); |
|
|
156 retval(0) = ArrayN<T> (dv); |
|
5360
|
157 return retval; |
|
|
158 } |
|
4996
|
159 |
|
5275
|
160 octave_idx_type ns = dv(dim); |
|
|
161 octave_idx_type iter = dv.numel () / ns; |
|
|
162 octave_idx_type stride = 1; |
|
|
163 for (int i = 0; i < dim; i++) |
|
4996
|
164 stride *= dv(i); |
|
|
165 |
|
|
166 T *v = m.fortran_vec (); |
|
|
167 octave_sort<vec_index<T> *> indexed_sort; |
|
|
168 |
|
|
169 if (mode == ASCENDING) |
|
|
170 indexed_sort.set_compare (ascending_compare); |
|
|
171 else if (mode == DESCENDING) |
|
|
172 indexed_sort.set_compare (descending_compare); |
|
|
173 |
|
|
174 OCTAVE_LOCAL_BUFFER (vec_index<T> *, vi, ns); |
|
|
175 OCTAVE_LOCAL_BUFFER (vec_index<T>, vix, ns); |
|
|
176 |
|
5275
|
177 for (octave_idx_type i = 0; i < ns; i++) |
|
4996
|
178 vi[i] = &vix[i]; |
|
|
179 |
|
|
180 NDArray idx (dv); |
|
|
181 |
|
|
182 if (stride == 1) |
|
|
183 { |
|
5275
|
184 for (octave_idx_type j = 0; j < iter; j++) |
|
4996
|
185 { |
|
5275
|
186 octave_idx_type offset = j * ns; |
|
4996
|
187 |
|
5275
|
188 for (octave_idx_type i = 0; i < ns; i++) |
|
4996
|
189 { |
|
|
190 vi[i]->vec = v[i]; |
|
|
191 vi[i]->indx = i + 1; |
|
|
192 } |
|
|
193 |
|
|
194 indexed_sort.sort (vi, ns); |
|
|
195 |
|
5275
|
196 for (octave_idx_type i = 0; i < ns; i++) |
|
4996
|
197 { |
|
|
198 v[i] = vi[i]->vec; |
|
|
199 idx(i + offset) = vi[i]->indx; |
|
|
200 } |
|
|
201 v += ns; |
|
|
202 } |
|
|
203 } |
|
|
204 else |
|
|
205 { |
|
5275
|
206 for (octave_idx_type j = 0; j < iter; j++) |
|
4996
|
207 { |
|
5275
|
208 octave_idx_type offset = j; |
|
|
209 octave_idx_type offset2 = 0; |
|
4996
|
210 while (offset >= stride) |
|
|
211 { |
|
|
212 offset -= stride; |
|
|
213 offset2++; |
|
|
214 } |
|
|
215 offset += offset2 * stride * ns; |
|
|
216 |
|
5275
|
217 for (octave_idx_type i = 0; i < ns; i++) |
|
4996
|
218 { |
|
|
219 vi[i]->vec = v[i*stride + offset]; |
|
|
220 vi[i]->indx = i + 1; |
|
|
221 } |
|
|
222 |
|
|
223 indexed_sort.sort (vi, ns); |
|
|
224 |
|
5275
|
225 for (octave_idx_type i = 0; i < ns; i++) |
|
4996
|
226 { |
|
|
227 v[i*stride+offset] = vi[i]->vec; |
|
|
228 idx(i*stride+offset) = vi[i]->indx; |
|
|
229 } |
|
|
230 } |
|
|
231 } |
|
|
232 |
|
4997
|
233 retval(1) = idx; |
|
|
234 retval(0) = octave_value (m); |
|
4999
|
235 |
|
4996
|
236 return retval; |
|
|
237 } |
|
|
238 |
|
4853
|
239 // If we have IEEE 754 data format, then we can use the trick of |
|
|
240 // casting doubles as unsigned eight byte integers, and with a little |
|
|
241 // bit of magic we can automatically sort the NaN's correctly. |
|
2928
|
242 |
|
5828
|
243 #if defined (HAVE_IEEE754_DATA_FORMAT) |
|
4850
|
244 |
|
5828
|
245 static inline uint64_t |
|
|
246 FloatFlip (uint64_t f) |
|
2928
|
247 { |
|
5828
|
248 uint64_t mask |
|
|
249 = -static_cast<int64_t>(f >> 63) | 0x8000000000000000ULL; |
|
4999
|
250 |
|
4850
|
251 return f ^ mask; |
|
|
252 } |
|
2928
|
253 |
|
5828
|
254 static inline uint64_t |
|
|
255 IFloatFlip (uint64_t f) |
|
4850
|
256 { |
|
5828
|
257 uint64_t mask = ((f >> 63) - 1) | 0x8000000000000000ULL; |
|
4999
|
258 |
|
4850
|
259 return f ^ mask; |
|
|
260 } |
|
2928
|
261 |
|
5009
|
262 template <> |
|
4996
|
263 bool |
|
5828
|
264 ascending_compare (uint64_t a, |
|
|
265 uint64_t b) |
|
4850
|
266 { |
|
4996
|
267 return (a < b); |
|
|
268 } |
|
2928
|
269 |
|
5009
|
270 template <> |
|
4850
|
271 bool |
|
5828
|
272 ascending_compare (vec_index<uint64_t> *a, |
|
|
273 vec_index<uint64_t> *b) |
|
4850
|
274 { |
|
|
275 return (a->vec < b->vec); |
|
2928
|
276 } |
|
|
277 |
|
5009
|
278 template <> |
|
4996
|
279 bool |
|
5828
|
280 descending_compare (uint64_t a, |
|
|
281 uint64_t b) |
|
4850
|
282 { |
|
4996
|
283 return (a > b); |
|
4850
|
284 } |
|
2928
|
285 |
|
5009
|
286 template <> |
|
4850
|
287 bool |
|
5828
|
288 descending_compare (vec_index<uint64_t> *a, |
|
|
289 vec_index<uint64_t> *b) |
|
4850
|
290 { |
|
4996
|
291 return (a->vec > b->vec); |
|
4850
|
292 } |
|
|
293 |
|
5828
|
294 template class octave_sort<uint64_t>; |
|
|
295 template class vec_index<uint64_t>; |
|
|
296 template class octave_sort<vec_index<uint64_t> *>; |
|
2928
|
297 |
|
4996
|
298 template <> |
|
5275
|
299 octave_value |
|
4996
|
300 mx_sort (ArrayN<double> &m, int dim, sortmode mode) |
|
2928
|
301 { |
|
4998
|
302 octave_value retval; |
|
2928
|
303 |
|
5366
|
304 dim_vector dv = m.dims (); |
|
|
305 |
|
4850
|
306 if (m.length () < 1) |
|
5366
|
307 return ArrayN<double> (dv); |
|
4850
|
308 |
|
5275
|
309 octave_idx_type ns = dv(dim); |
|
|
310 octave_idx_type iter = dv.numel () / ns; |
|
|
311 octave_idx_type stride = 1; |
|
|
312 for (int i = 0; i < dim; i++) |
|
4850
|
313 stride *= dv(i); |
|
|
314 |
|
|
315 double *v = m.fortran_vec (); |
|
|
316 |
|
5828
|
317 uint64_t *p = reinterpret_cast<uint64_t *> (v); |
|
4850
|
318 |
|
5828
|
319 octave_sort<uint64_t> sort; |
|
4996
|
320 |
|
|
321 if (mode == ASCENDING) |
|
|
322 sort.set_compare (ascending_compare); |
|
|
323 else if (mode == DESCENDING) |
|
|
324 sort.set_compare (descending_compare); |
|
|
325 |
|
|
326 if (stride == 1) |
|
4850
|
327 { |
|
5275
|
328 for (octave_idx_type j = 0; j < iter; j++) |
|
4996
|
329 { |
|
|
330 // Flip the data in the vector so that int compares on |
|
|
331 // IEEE754 give the correct ordering. |
|
|
332 |
|
5275
|
333 for (octave_idx_type i = 0; i < ns; i++) |
|
4996
|
334 p[i] = FloatFlip (p[i]); |
|
|
335 |
|
|
336 sort.sort (p, ns); |
|
|
337 |
|
|
338 // Flip the data out of the vector so that int compares |
|
|
339 // on IEEE754 give the correct ordering. |
|
|
340 |
|
5275
|
341 for (octave_idx_type i = 0; i < ns; i++) |
|
4996
|
342 p[i] = IFloatFlip (p[i]); |
|
4850
|
343 |
|
4996
|
344 // There are two representations of NaN. One will be |
|
|
345 // sorted to the beginning of the vector and the other |
|
5009
|
346 // to the end. If it will be sorted incorrectly, fix |
|
|
347 // things up. |
|
4850
|
348 |
|
5009
|
349 if (lo_ieee_signbit (octave_NaN)) |
|
6484
|
350 { |
|
|
351 if (mode == UNDEFINED || mode == ASCENDING) |
|
|
352 { |
|
|
353 octave_idx_type i = 0; |
|
|
354 double *vtmp = reinterpret_cast<double *> (p); |
|
|
355 while (xisnan (vtmp[i++]) && i < ns); |
|
|
356 for (octave_idx_type l = 0; l < ns - i + 1; l++) |
|
|
357 vtmp[l] = vtmp[l+i-1]; |
|
|
358 for (octave_idx_type l = ns - i + 1; l < ns; l++) |
|
|
359 vtmp[l] = octave_NaN; |
|
|
360 } |
|
|
361 else |
|
|
362 { |
|
|
363 octave_idx_type i = ns; |
|
|
364 double *vtmp = reinterpret_cast<double *> (p); |
|
|
365 while (xisnan (vtmp[--i]) && i > 0); |
|
|
366 for (octave_idx_type l = i; l >= 0; l--) |
|
|
367 vtmp[l-i+ns-1] = vtmp[l]; |
|
|
368 for (octave_idx_type l = 0; l < ns - i - 1; l++) |
|
|
369 vtmp[l] = octave_NaN; |
|
|
370 } |
|
|
371 } |
|
4996
|
372 |
|
|
373 p += ns; |
|
|
374 } |
|
|
375 } |
|
|
376 else |
|
|
377 { |
|
5828
|
378 OCTAVE_LOCAL_BUFFER (uint64_t, vi, ns); |
|
4996
|
379 |
|
5275
|
380 for (octave_idx_type j = 0; j < iter; j++) |
|
4850
|
381 { |
|
5275
|
382 octave_idx_type offset = j; |
|
|
383 octave_idx_type offset2 = 0; |
|
4850
|
384 while (offset >= stride) |
|
|
385 { |
|
|
386 offset -= stride; |
|
|
387 offset2++; |
|
|
388 } |
|
|
389 offset += offset2 * stride * ns; |
|
2928
|
390 |
|
4853
|
391 // Flip the data in the vector so that int compares on |
|
|
392 // IEEE754 give the correct ordering. |
|
|
393 |
|
5275
|
394 for (octave_idx_type i = 0; i < ns; i++) |
|
4996
|
395 vi[i] = FloatFlip (p[i*stride + offset]); |
|
4850
|
396 |
|
4996
|
397 sort.sort (vi, ns); |
|
4850
|
398 |
|
4996
|
399 // Flip the data out of the vector so that int compares |
|
|
400 // on IEEE754 give the correct ordering. |
|
4853
|
401 |
|
5275
|
402 for (octave_idx_type i = 0; i < ns; i++) |
|
4996
|
403 p[i*stride + offset] = IFloatFlip (vi[i]); |
|
|
404 |
|
|
405 // There are two representations of NaN. One will be |
|
|
406 // sorted to the beginning of the vector and the other |
|
|
407 // to the end. If it will be sorted to the beginning, |
|
|
408 // fix things up. |
|
2928
|
409 |
|
5009
|
410 if (lo_ieee_signbit (octave_NaN)) |
|
6484
|
411 { |
|
|
412 if (mode == UNDEFINED || mode == ASCENDING) |
|
|
413 { |
|
|
414 octave_idx_type i = 0; |
|
|
415 while (xisnan (v[i++*stride + offset]) && i < ns); |
|
|
416 for (octave_idx_type l = 0; l < ns - i + 1; l++) |
|
|
417 v[l*stride + offset] = v[(l+i-1)*stride + offset]; |
|
|
418 for (octave_idx_type l = ns - i + 1; l < ns; l++) |
|
|
419 v[l*stride + offset] = octave_NaN; |
|
|
420 } |
|
|
421 else |
|
|
422 { |
|
|
423 octave_idx_type i = ns; |
|
|
424 while (xisnan (v[--i*stride + offset]) && i > 0); |
|
|
425 for (octave_idx_type l = i; l >= 0; l--) |
|
|
426 v[(l-i+ns-1)*stride + offset] = v[l*stride + offset]; |
|
|
427 for (octave_idx_type l = 0; l < ns - i - 1; l++) |
|
|
428 v[l*stride + offset] = octave_NaN; |
|
|
429 } |
|
|
430 } |
|
2928
|
431 } |
|
|
432 } |
|
|
433 |
|
4998
|
434 retval = m; |
|
|
435 |
|
2928
|
436 return retval; |
|
|
437 } |
|
|
438 |
|
5275
|
439 // Should other overloaded functions have their static keywords removed? |
|
4996
|
440 template <> |
|
5275
|
441 octave_value_list |
|
4996
|
442 mx_sort_indexed (ArrayN<double> &m, int dim, sortmode mode) |
|
2928
|
443 { |
|
|
444 octave_value_list retval; |
|
|
445 |
|
5366
|
446 dim_vector dv = m.dims (); |
|
|
447 |
|
4850
|
448 if (m.length () < 1) |
|
5360
|
449 { |
|
5419
|
450 retval(1) = NDArray (dv); |
|
|
451 retval(0) = ArrayN<double> (dv); |
|
5360
|
452 return retval; |
|
|
453 } |
|
4850
|
454 |
|
5275
|
455 octave_idx_type ns = dv(dim); |
|
|
456 octave_idx_type iter = dv.numel () / ns; |
|
|
457 octave_idx_type stride = 1; |
|
|
458 for (int i = 0; i < dim; i++) |
|
4850
|
459 stride *= dv(i); |
|
2928
|
460 |
|
4996
|
461 double *v = m.fortran_vec (); |
|
4850
|
462 |
|
5828
|
463 uint64_t *p = reinterpret_cast<uint64_t *> (v); |
|
4996
|
464 |
|
5828
|
465 octave_sort<vec_index<uint64_t> *> indexed_sort; |
|
2928
|
466 |
|
4996
|
467 if (mode == ASCENDING) |
|
|
468 indexed_sort.set_compare (ascending_compare); |
|
|
469 else if (mode == DESCENDING) |
|
|
470 indexed_sort.set_compare (descending_compare); |
|
4850
|
471 |
|
5828
|
472 OCTAVE_LOCAL_BUFFER (vec_index<uint64_t> *, vi, ns); |
|
|
473 OCTAVE_LOCAL_BUFFER (vec_index<uint64_t>, vix, ns); |
|
4996
|
474 |
|
5275
|
475 for (octave_idx_type i = 0; i < ns; i++) |
|
4850
|
476 vi[i] = &vix[i]; |
|
2928
|
477 |
|
4850
|
478 NDArray idx (dv); |
|
4996
|
479 |
|
5275
|
480 for (octave_idx_type j = 0; j < iter; j++) |
|
2928
|
481 { |
|
5275
|
482 octave_idx_type offset = j; |
|
|
483 octave_idx_type offset2 = 0; |
|
4996
|
484 while (offset >= stride) |
|
2928
|
485 { |
|
4996
|
486 offset -= stride; |
|
|
487 offset2++; |
|
|
488 } |
|
|
489 offset += offset2 * stride * ns; |
|
2928
|
490 |
|
4996
|
491 // Flip the data in the vector so that int compares on |
|
|
492 // IEEE754 give the correct ordering. |
|
|
493 |
|
5275
|
494 for (octave_idx_type i = 0; i < ns; i++) |
|
4996
|
495 { |
|
|
496 vi[i]->vec = FloatFlip (p[i*stride + offset]); |
|
|
497 vi[i]->indx = i + 1; |
|
4850
|
498 } |
|
4996
|
499 |
|
|
500 indexed_sort.sort (vi, ns); |
|
|
501 |
|
|
502 // Flip the data out of the vector so that int compares on |
|
|
503 // IEEE754 give the correct ordering |
|
|
504 |
|
5275
|
505 for (octave_idx_type i = 0; i < ns; i++) |
|
4850
|
506 { |
|
4996
|
507 p[i*stride + offset] = IFloatFlip (vi[i]->vec); |
|
|
508 idx(i*stride + offset) = vi[i]->indx; |
|
|
509 } |
|
|
510 |
|
|
511 // There are two representations of NaN. One will be sorted |
|
|
512 // to the beginning of the vector and the other to the end. |
|
|
513 // If it will be sorted to the beginning, fix things up. |
|
2928
|
514 |
|
5009
|
515 if (lo_ieee_signbit (octave_NaN)) |
|
6484
|
516 { |
|
|
517 if (mode == UNDEFINED || mode == ASCENDING) |
|
|
518 { |
|
|
519 octave_idx_type i = 0; |
|
|
520 while (xisnan (v[i++*stride+offset]) && i < ns); |
|
|
521 OCTAVE_LOCAL_BUFFER (double, itmp, i - 1); |
|
|
522 for (octave_idx_type l = 0; l < i -1; l++) |
|
|
523 itmp[l] = idx(l*stride + offset); |
|
|
524 for (octave_idx_type l = 0; l < ns - i + 1; l++) |
|
|
525 { |
|
|
526 v[l*stride + offset] = v[(l+i-1)*stride + offset]; |
|
|
527 idx(l*stride + offset) = idx((l+i-1)*stride + offset); |
|
|
528 } |
|
|
529 for (octave_idx_type k = 0, l = ns - i + 1; l < ns; l++, k++) |
|
|
530 { |
|
|
531 v[l*stride + offset] = octave_NaN; |
|
|
532 idx(l*stride + offset) = itmp[k]; |
|
|
533 } |
|
|
534 } |
|
|
535 else |
|
|
536 { |
|
|
537 octave_idx_type i = ns; |
|
|
538 while (xisnan (v[--i*stride+offset]) && i > 0); |
|
|
539 OCTAVE_LOCAL_BUFFER (double, itmp, ns - i - 1); |
|
|
540 for (octave_idx_type l = 0; l < ns - i -1; l++) |
|
|
541 itmp[l] = idx((l+i+1)*stride + offset); |
|
|
542 for (octave_idx_type l = i; l >= 0; l--) |
|
|
543 { |
|
|
544 v[(l-i+ns-1)*stride + offset] = v[l*stride + offset]; |
|
|
545 idx((l-i+ns-1)*stride + offset) = idx(l*stride + offset); |
|
|
546 } |
|
|
547 for (octave_idx_type k = 0, l = 0; l < ns - i - 1; l++, k++) |
|
|
548 { |
|
|
549 v[l*stride + offset] = octave_NaN; |
|
|
550 idx(l*stride + offset) = itmp[k]; |
|
|
551 } |
|
|
552 } |
|
|
553 } |
|
2928
|
554 } |
|
|
555 |
|
4997
|
556 retval(1) = idx; |
|
|
557 retval(0) = m; |
|
4998
|
558 |
|
2928
|
559 return retval; |
|
|
560 } |
|
|
561 |
|
4996
|
562 #else |
|
|
563 |
|
5009
|
564 template <> |
|
4996
|
565 bool |
|
|
566 ascending_compare (double a, double b) |
|
4991
|
567 { |
|
4997
|
568 return (xisnan (b) || (a < b)); |
|
4996
|
569 } |
|
|
570 |
|
5009
|
571 template <> |
|
4996
|
572 bool |
|
|
573 ascending_compare (vec_index<double> *a, vec_index<double> *b) |
|
|
574 { |
|
4997
|
575 return (xisnan (b->vec) || (a->vec < b->vec)); |
|
4996
|
576 } |
|
|
577 |
|
5009
|
578 template <> |
|
4996
|
579 bool |
|
|
580 descending_compare (double a, double b) |
|
|
581 { |
|
4997
|
582 return (xisnan (a) || (a > b)); |
|
4996
|
583 } |
|
4991
|
584 |
|
5009
|
585 template <> |
|
4991
|
586 bool |
|
4996
|
587 descending_compare (vec_index<double> *a, vec_index<double> *b) |
|
|
588 { |
|
4997
|
589 return (xisnan (a->vec) || (a->vec > b->vec)); |
|
4996
|
590 } |
|
|
591 |
|
|
592 template class octave_sort<double>; |
|
|
593 template class vec_index<double>; |
|
|
594 template class octave_sort<vec_index<double> *>; |
|
|
595 |
|
|
596 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL) |
|
|
597 static octave_value_list |
|
|
598 mx_sort (ArrayN<double> &m, int dim, sortmode mode); |
|
|
599 |
|
|
600 static octave_value_list |
|
|
601 mx_sort_indexed (ArrayN<double> &m, int dim, sortmode mode); |
|
|
602 #endif |
|
|
603 #endif |
|
|
604 |
|
|
605 // std::abs(Inf) returns NaN!! |
|
|
606 static inline double |
|
|
607 xabs (const Complex& x) |
|
|
608 { |
|
|
609 return (xisinf (x.real ()) || xisinf (x.imag ())) ? octave_Inf : abs (x); |
|
|
610 } |
|
|
611 |
|
5009
|
612 template <> |
|
4996
|
613 bool |
|
|
614 ascending_compare (vec_index<Complex> *a, vec_index<Complex> *b) |
|
|
615 { |
|
4997
|
616 return (xisnan (b->vec) |
|
|
617 || (xabs (a->vec) < xabs (b->vec)) |
|
|
618 || ((xabs (a->vec) == xabs (b->vec)) |
|
|
619 && (arg (a->vec) < arg (b->vec)))); |
|
4996
|
620 } |
|
|
621 |
|
5009
|
622 template <> |
|
4996
|
623 bool |
|
|
624 descending_compare (vec_index<Complex> *a, vec_index<Complex> *b) |
|
|
625 { |
|
4997
|
626 return (xisnan (a->vec) |
|
|
627 || (xabs (a->vec) > xabs (b->vec)) |
|
|
628 || ((xabs (a->vec) == xabs (b->vec)) |
|
|
629 && (arg (a->vec) > arg (b->vec)))); |
|
4996
|
630 } |
|
|
631 |
|
|
632 template class vec_index<Complex>; |
|
|
633 template class octave_sort<vec_index<Complex> *>; |
|
|
634 |
|
|
635 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL) |
|
|
636 static octave_value_list |
|
|
637 mx_sort_indexed (ArrayN<Complex> &m, int dim, sortmode mode); |
|
|
638 #endif |
|
|
639 |
|
4991
|
640 template class octave_sort<char>; |
|
4996
|
641 template class vec_index<char>; |
|
|
642 template class octave_sort<vec_index<char> *>; |
|
|
643 |
|
|
644 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL) |
|
5009
|
645 bool |
|
|
646 ascending_compare (char a, char b); |
|
|
647 |
|
|
648 bool |
|
|
649 ascending_compare (vec_index<char> *a, vec_index<char> *b); |
|
|
650 |
|
|
651 bool |
|
|
652 descending_compare (char a, char b); |
|
|
653 |
|
|
654 bool |
|
|
655 descending_compare (vec_index<char> *a, vec_index<char> *b); |
|
|
656 |
|
4996
|
657 static octave_value_list |
|
|
658 mx_sort (ArrayN<char> &m, int dim, sortmode mode); |
|
4991
|
659 |
|
|
660 static octave_value_list |
|
4996
|
661 mx_sort_indexed (ArrayN<char> &m, int dim, sortmode mode); |
|
|
662 #endif |
|
4991
|
663 |
|
5009
|
664 template <> |
|
4996
|
665 bool |
|
|
666 ascending_compare (vec_index<octave_value> *a, vec_index<octave_value> *b) |
|
|
667 { |
|
|
668 return (a->vec.string_value () < b->vec.string_value ()); |
|
|
669 } |
|
4991
|
670 |
|
5009
|
671 template <> |
|
4996
|
672 bool |
|
|
673 descending_compare (vec_index<octave_value> *a, vec_index<octave_value> *b) |
|
|
674 { |
|
|
675 return (a->vec.string_value () > b->vec.string_value ()); |
|
|
676 } |
|
4991
|
677 |
|
4996
|
678 template class vec_index<octave_value>; |
|
|
679 template class octave_sort<vec_index<octave_value> *>; |
|
4991
|
680 |
|
4996
|
681 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL) |
|
|
682 static octave_value_list |
|
|
683 mx_sort_indexed (ArrayN<octave_value> &m, int dim, sortmode mode); |
|
|
684 #endif |
|
4991
|
685 |
|
2928
|
686 DEFUN_DLD (sort, args, nargout, |
|
3369
|
687 "-*- texinfo -*-\n\ |
|
|
688 @deftypefn {Loadable Function} {[@var{s}, @var{i}] =} sort (@var{x})\n\ |
|
4850
|
689 @deftypefnx {Loadable Function} {[@var{s}, @var{i}] =} sort (@var{x}, @var{dim})\n\ |
|
5000
|
690 @deftypefnx {Loadable Function} {[@var{s}, @var{i}] =} sort (@var{x}, @var{mode})\n\ |
|
|
691 @deftypefnx {Loadable Function} {[@var{s}, @var{i}] =} sort (@var{x}, @var{dim}, @var{mode})\n\ |
|
5499
|
692 Return a copy of @var{x} with the elements arranged in increasing\n\ |
|
|
693 order. For matrices, @code{sort} orders the elements in each column.\n\ |
|
3369
|
694 \n\ |
|
|
695 For example,\n\ |
|
|
696 \n\ |
|
|
697 @example\n\ |
|
|
698 @group\n\ |
|
|
699 sort ([1, 2; 2, 3; 3, 1])\n\ |
|
|
700 @result{} 1 1\n\ |
|
|
701 2 2\n\ |
|
|
702 3 3\n\ |
|
|
703 @end group\n\ |
|
|
704 @end example\n\ |
|
2928
|
705 \n\ |
|
3369
|
706 The @code{sort} function may also be used to produce a matrix\n\ |
|
|
707 containing the original row indices of the elements in the sorted\n\ |
|
|
708 matrix. For example,\n\ |
|
|
709 \n\ |
|
|
710 @example\n\ |
|
|
711 @group\n\ |
|
|
712 [s, i] = sort ([1, 2; 2, 3; 3, 1])\n\ |
|
|
713 @result{} s = 1 1\n\ |
|
|
714 2 2\n\ |
|
|
715 3 3\n\ |
|
|
716 @result{} i = 1 3\n\ |
|
|
717 2 1\n\ |
|
|
718 3 2\n\ |
|
|
719 @end group\n\ |
|
|
720 @end example\n\ |
|
4850
|
721 \n\ |
|
|
722 If the optional argument @var{dim} is given, then the matrix is sorted\n\ |
|
4996
|
723 along the dimension defined by @var{dim}. The optional argument @code{mode}\n\ |
|
|
724 defines the order in which the values will be sorted. Valid values of\n\ |
|
|
725 @code{mode} are `ascend' or `descend'.\n\ |
|
4850
|
726 \n\ |
|
|
727 For equal elements, the indices are such that the equal elements are listed\n\ |
|
|
728 in the order that appeared in the original list.\n\ |
|
|
729 \n\ |
|
4996
|
730 The @code{sort} function may also be used to sort strings and cell arrays\n\ |
|
5499
|
731 of strings, in which case the dictionary order of the strings is used.\n\ |
|
4996
|
732 \n\ |
|
4850
|
733 The algorithm used in @code{sort} is optimized for the sorting of partially\n\ |
|
|
734 ordered lists.\n\ |
|
3369
|
735 @end deftypefn") |
|
2928
|
736 { |
|
|
737 octave_value_list retval; |
|
|
738 |
|
|
739 int nargin = args.length (); |
|
4996
|
740 sortmode smode = ASCENDING; |
|
2928
|
741 |
|
5001
|
742 if (nargin < 1 || nargin > 3) |
|
2928
|
743 { |
|
5823
|
744 print_usage (); |
|
2928
|
745 return retval; |
|
|
746 } |
|
|
747 |
|
4850
|
748 bool return_idx = nargout > 1; |
|
2928
|
749 |
|
|
750 octave_value arg = args(0); |
|
|
751 |
|
4850
|
752 int dim = 0; |
|
4996
|
753 if (nargin > 1) |
|
|
754 { |
|
|
755 if (args(1).is_string ()) |
|
|
756 { |
|
|
757 std::string mode = args(1).string_value(); |
|
|
758 if (mode == "ascend") |
|
|
759 smode = ASCENDING; |
|
|
760 else if (mode == "descend") |
|
|
761 smode = DESCENDING; |
|
|
762 else |
|
|
763 { |
|
5000
|
764 error ("sort: mode must be either \"ascend\" or \"descend\""); |
|
4996
|
765 return retval; |
|
|
766 } |
|
|
767 } |
|
|
768 else |
|
|
769 dim = args(1).nint_value () - 1; |
|
|
770 } |
|
|
771 |
|
|
772 if (nargin > 2) |
|
|
773 { |
|
|
774 if (args(1).is_string ()) |
|
|
775 { |
|
5823
|
776 print_usage (); |
|
4996
|
777 return retval; |
|
|
778 } |
|
|
779 |
|
|
780 if (! args(2).is_string ()) |
|
|
781 { |
|
|
782 error ("sort: mode must be a string"); |
|
|
783 return retval; |
|
|
784 } |
|
|
785 std::string mode = args(2).string_value(); |
|
|
786 if (mode == "ascend") |
|
|
787 smode = ASCENDING; |
|
|
788 else if (mode == "descend") |
|
|
789 smode = DESCENDING; |
|
|
790 else |
|
|
791 { |
|
5000
|
792 error ("sort: mode must be either \"ascend\" or \"descend\""); |
|
4996
|
793 return retval; |
|
|
794 } |
|
|
795 } |
|
4850
|
796 |
|
5275
|
797 dim_vector dv = arg.dims (); |
|
4850
|
798 if (error_state) |
|
|
799 { |
|
|
800 gripe_wrong_type_arg ("sort", arg); |
|
|
801 return retval; |
|
|
802 } |
|
4996
|
803 if (nargin == 1 || args(1).is_string ()) |
|
4850
|
804 { |
|
|
805 // Find first non singleton dimension |
|
|
806 for (int i = 0; i < dv.length (); i++) |
|
|
807 if (dv(i) > 1) |
|
|
808 { |
|
|
809 dim = i; |
|
|
810 break; |
|
|
811 } |
|
|
812 } |
|
|
813 else |
|
|
814 { |
|
|
815 if (dim < 0 || dim > dv.length () - 1) |
|
|
816 { |
|
|
817 error ("sort: dim must be a valid dimension"); |
|
|
818 return retval; |
|
|
819 } |
|
|
820 } |
|
|
821 |
|
2928
|
822 if (arg.is_real_type ()) |
|
|
823 { |
|
4850
|
824 NDArray m = arg.array_value (); |
|
2928
|
825 |
|
|
826 if (! error_state) |
|
4996
|
827 { |
|
|
828 #ifdef HAVE_IEEE754_DATA_FORMAT |
|
|
829 // As operator > gives the right result, can special case here |
|
|
830 if (! return_idx && smode == ASCENDING) |
|
|
831 retval = mx_sort (m, dim); |
|
|
832 else |
|
|
833 #endif |
|
|
834 { |
|
|
835 if (return_idx) |
|
|
836 retval = mx_sort_indexed (m, dim, smode); |
|
|
837 else |
|
|
838 retval = mx_sort (m, dim, smode); |
|
|
839 } |
|
|
840 } |
|
2928
|
841 } |
|
|
842 else if (arg.is_complex_type ()) |
|
|
843 { |
|
4850
|
844 ComplexNDArray cm = arg.complex_array_value (); |
|
2928
|
845 |
|
4996
|
846 // Don't have unindexed version as no ">" operator |
|
2928
|
847 if (! error_state) |
|
4996
|
848 retval = mx_sort_indexed (cm, dim, smode); |
|
2928
|
849 } |
|
4991
|
850 else if (arg.is_string ()) |
|
|
851 { |
|
|
852 charNDArray chm = arg.char_array_value (); |
|
|
853 |
|
|
854 if (! error_state) |
|
4996
|
855 { |
|
|
856 // As operator > gives the right result, can special case here |
|
|
857 if (! return_idx && smode == ASCENDING) |
|
|
858 retval = mx_sort (chm, dim); |
|
|
859 else |
|
|
860 { |
|
|
861 if (return_idx) |
|
|
862 retval = mx_sort_indexed (chm, dim, smode); |
|
|
863 else |
|
|
864 retval = mx_sort (chm, dim, smode); |
|
|
865 } |
|
|
866 |
|
5775
|
867 // FIXME It would have been better to call |
|
4996
|
868 // "octave_value(m, true)" but how can that be done |
|
|
869 // within the template |
|
|
870 retval(0) = retval(0).convert_to_str (false, true); |
|
|
871 } |
|
|
872 } |
|
|
873 else if (arg.is_cell ()) |
|
|
874 { |
|
|
875 Cell cellm = arg.cell_value (); |
|
|
876 |
|
|
877 // Need to check that all elements are strings |
|
5275
|
878 for (octave_idx_type i = 0; i < cellm.numel (); i++) |
|
4996
|
879 if (! cellm(i).is_string ()) |
|
|
880 { |
|
|
881 gripe_wrong_type_arg ("sort", arg); |
|
|
882 break; |
|
|
883 } |
|
|
884 |
|
|
885 // Don't have unindexed version as ">" operator doesn't return bool |
|
|
886 if (!error_state) |
|
|
887 retval = mx_sort_indexed (cellm, dim, smode); |
|
4991
|
888 } |
|
2928
|
889 else |
|
|
890 gripe_wrong_type_arg ("sort", arg); |
|
|
891 |
|
|
892 return retval; |
|
|
893 } |
|
|
894 |
|
|
895 /* |
|
|
896 ;;; Local Variables: *** |
|
|
897 ;;; mode: C++ *** |
|
|
898 ;;; End: *** |
|
|
899 */ |
