|
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 |
|
4853
|
243 #if defined (HAVE_IEEE754_DATA_FORMAT) && defined (EIGHT_BYTE_INT) |
|
4850
|
244 |
|
4999
|
245 static inline unsigned EIGHT_BYTE_INT |
|
|
246 FloatFlip (unsigned EIGHT_BYTE_INT f) |
|
2928
|
247 { |
|
5000
|
248 unsigned EIGHT_BYTE_INT mask |
|
|
249 = -(EIGHT_BYTE_INT)(f >> 63) | 0x8000000000000000ULL; |
|
4999
|
250 |
|
4850
|
251 return f ^ mask; |
|
|
252 } |
|
2928
|
253 |
|
5009
|
254 static inline unsigned EIGHT_BYTE_INT |
|
5000
|
255 IFloatFlip (unsigned EIGHT_BYTE_INT f) |
|
4850
|
256 { |
|
|
257 unsigned EIGHT_BYTE_INT mask = ((f >> 63) - 1) | 0x8000000000000000ULL; |
|
4999
|
258 |
|
4850
|
259 return f ^ mask; |
|
|
260 } |
|
2928
|
261 |
|
5009
|
262 template <> |
|
4996
|
263 bool |
|
|
264 ascending_compare (unsigned EIGHT_BYTE_INT a, |
|
4998
|
265 unsigned EIGHT_BYTE_INT b) |
|
4850
|
266 { |
|
4996
|
267 return (a < b); |
|
|
268 } |
|
2928
|
269 |
|
5009
|
270 template <> |
|
4850
|
271 bool |
|
4996
|
272 ascending_compare (vec_index<unsigned EIGHT_BYTE_INT> *a, |
|
|
273 vec_index<unsigned EIGHT_BYTE_INT> *b) |
|
4850
|
274 { |
|
|
275 return (a->vec < b->vec); |
|
2928
|
276 } |
|
|
277 |
|
5009
|
278 template <> |
|
4996
|
279 bool |
|
|
280 descending_compare (unsigned EIGHT_BYTE_INT a, |
|
|
281 unsigned EIGHT_BYTE_INT b) |
|
4850
|
282 { |
|
4996
|
283 return (a > b); |
|
4850
|
284 } |
|
2928
|
285 |
|
5009
|
286 template <> |
|
4850
|
287 bool |
|
4996
|
288 descending_compare (vec_index<unsigned EIGHT_BYTE_INT> *a, |
|
|
289 vec_index<unsigned EIGHT_BYTE_INT> *b) |
|
4850
|
290 { |
|
4996
|
291 return (a->vec > b->vec); |
|
4850
|
292 } |
|
|
293 |
|
4996
|
294 template class octave_sort<unsigned EIGHT_BYTE_INT>; |
|
|
295 template class vec_index<unsigned EIGHT_BYTE_INT>; |
|
|
296 template class octave_sort<vec_index<unsigned EIGHT_BYTE_INT> *>; |
|
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 |
|
5010
|
317 unsigned EIGHT_BYTE_INT *p = reinterpret_cast<unsigned EIGHT_BYTE_INT *> (v); |
|
4850
|
318 |
|
4996
|
319 octave_sort<unsigned EIGHT_BYTE_INT> sort; |
|
|
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)) |
|
|
350 if (mode == UNDEFINED || mode == ASCENDING) |
|
|
351 { |
|
5275
|
352 octave_idx_type i = 0; |
|
5009
|
353 double *vtmp = (double *)p; |
|
|
354 while (xisnan (vtmp[i++]) && i < ns); |
|
5275
|
355 for (octave_idx_type l = 0; l < ns - i + 1; l++) |
|
5009
|
356 vtmp[l] = vtmp[l+i-1]; |
|
5275
|
357 for (octave_idx_type l = ns - i + 1; l < ns; l++) |
|
5009
|
358 vtmp[l] = octave_NaN; |
|
|
359 } |
|
|
360 else |
|
|
361 { |
|
5275
|
362 octave_idx_type i = ns; |
|
5009
|
363 double *vtmp = (double *)p; |
|
|
364 while (xisnan (vtmp[--i]) && i > 0); |
|
5275
|
365 for (octave_idx_type l = i; l >= 0; l--) |
|
5009
|
366 vtmp[l-i+ns-1] = vtmp[l]; |
|
5275
|
367 for (octave_idx_type l = 0; l < ns - i - 1; l++) |
|
5009
|
368 vtmp[l] = octave_NaN; |
|
|
369 } |
|
4996
|
370 |
|
|
371 p += ns; |
|
|
372 } |
|
|
373 } |
|
|
374 else |
|
|
375 { |
|
|
376 OCTAVE_LOCAL_BUFFER (unsigned EIGHT_BYTE_INT, vi, ns); |
|
|
377 |
|
5275
|
378 for (octave_idx_type j = 0; j < iter; j++) |
|
4850
|
379 { |
|
5275
|
380 octave_idx_type offset = j; |
|
|
381 octave_idx_type offset2 = 0; |
|
4850
|
382 while (offset >= stride) |
|
|
383 { |
|
|
384 offset -= stride; |
|
|
385 offset2++; |
|
|
386 } |
|
|
387 offset += offset2 * stride * ns; |
|
2928
|
388 |
|
4853
|
389 // Flip the data in the vector so that int compares on |
|
|
390 // IEEE754 give the correct ordering. |
|
|
391 |
|
5275
|
392 for (octave_idx_type i = 0; i < ns; i++) |
|
4996
|
393 vi[i] = FloatFlip (p[i*stride + offset]); |
|
4850
|
394 |
|
4996
|
395 sort.sort (vi, ns); |
|
4850
|
396 |
|
4996
|
397 // Flip the data out of the vector so that int compares |
|
|
398 // on IEEE754 give the correct ordering. |
|
4853
|
399 |
|
5275
|
400 for (octave_idx_type i = 0; i < ns; i++) |
|
4996
|
401 p[i*stride + offset] = IFloatFlip (vi[i]); |
|
|
402 |
|
|
403 // There are two representations of NaN. One will be |
|
|
404 // sorted to the beginning of the vector and the other |
|
|
405 // to the end. If it will be sorted to the beginning, |
|
|
406 // fix things up. |
|
2928
|
407 |
|
5009
|
408 if (lo_ieee_signbit (octave_NaN)) |
|
|
409 if (mode == UNDEFINED || mode == ASCENDING) |
|
|
410 { |
|
5275
|
411 octave_idx_type i = 0; |
|
5009
|
412 while (xisnan (v[i++*stride + offset]) && i < ns); |
|
5275
|
413 for (octave_idx_type l = 0; l < ns - i + 1; l++) |
|
5009
|
414 v[l*stride + offset] = v[(l+i-1)*stride + offset]; |
|
5275
|
415 for (octave_idx_type l = ns - i + 1; l < ns; l++) |
|
5009
|
416 v[l*stride + offset] = octave_NaN; |
|
|
417 } |
|
|
418 else |
|
|
419 { |
|
5275
|
420 octave_idx_type i = ns; |
|
5009
|
421 while (xisnan (v[--i*stride + offset]) && i > 0); |
|
5275
|
422 for (octave_idx_type l = i; l >= 0; l--) |
|
5009
|
423 v[(l-i+ns-1)*stride + offset] = v[l*stride + offset]; |
|
5275
|
424 for (octave_idx_type l = 0; l < ns - i - 1; l++) |
|
5009
|
425 v[l*stride + offset] = octave_NaN; |
|
|
426 } |
|
2928
|
427 } |
|
|
428 } |
|
|
429 |
|
4998
|
430 retval = m; |
|
|
431 |
|
2928
|
432 return retval; |
|
|
433 } |
|
|
434 |
|
5275
|
435 // Should other overloaded functions have their static keywords removed? |
|
4996
|
436 template <> |
|
5275
|
437 octave_value_list |
|
4996
|
438 mx_sort_indexed (ArrayN<double> &m, int dim, sortmode mode) |
|
2928
|
439 { |
|
|
440 octave_value_list retval; |
|
|
441 |
|
5366
|
442 dim_vector dv = m.dims (); |
|
|
443 |
|
4850
|
444 if (m.length () < 1) |
|
5360
|
445 { |
|
5419
|
446 retval(1) = NDArray (dv); |
|
|
447 retval(0) = ArrayN<double> (dv); |
|
5360
|
448 return retval; |
|
|
449 } |
|
4850
|
450 |
|
5275
|
451 octave_idx_type ns = dv(dim); |
|
|
452 octave_idx_type iter = dv.numel () / ns; |
|
|
453 octave_idx_type stride = 1; |
|
|
454 for (int i = 0; i < dim; i++) |
|
4850
|
455 stride *= dv(i); |
|
2928
|
456 |
|
4996
|
457 double *v = m.fortran_vec (); |
|
4850
|
458 |
|
5010
|
459 unsigned EIGHT_BYTE_INT *p = reinterpret_cast<unsigned EIGHT_BYTE_INT *> (v); |
|
4996
|
460 |
|
|
461 octave_sort<vec_index<unsigned EIGHT_BYTE_INT> *> indexed_sort; |
|
2928
|
462 |
|
4996
|
463 if (mode == ASCENDING) |
|
|
464 indexed_sort.set_compare (ascending_compare); |
|
|
465 else if (mode == DESCENDING) |
|
|
466 indexed_sort.set_compare (descending_compare); |
|
4850
|
467 |
|
4996
|
468 OCTAVE_LOCAL_BUFFER (vec_index<unsigned EIGHT_BYTE_INT> *, vi, ns); |
|
|
469 OCTAVE_LOCAL_BUFFER (vec_index<unsigned EIGHT_BYTE_INT>, vix, ns); |
|
|
470 |
|
5275
|
471 for (octave_idx_type i = 0; i < ns; i++) |
|
4850
|
472 vi[i] = &vix[i]; |
|
2928
|
473 |
|
4850
|
474 NDArray idx (dv); |
|
4996
|
475 |
|
5275
|
476 for (octave_idx_type j = 0; j < iter; j++) |
|
2928
|
477 { |
|
5275
|
478 octave_idx_type offset = j; |
|
|
479 octave_idx_type offset2 = 0; |
|
4996
|
480 while (offset >= stride) |
|
2928
|
481 { |
|
4996
|
482 offset -= stride; |
|
|
483 offset2++; |
|
|
484 } |
|
|
485 offset += offset2 * stride * ns; |
|
2928
|
486 |
|
4996
|
487 // Flip the data in the vector so that int compares on |
|
|
488 // IEEE754 give the correct ordering. |
|
|
489 |
|
5275
|
490 for (octave_idx_type i = 0; i < ns; i++) |
|
4996
|
491 { |
|
|
492 vi[i]->vec = FloatFlip (p[i*stride + offset]); |
|
|
493 vi[i]->indx = i + 1; |
|
4850
|
494 } |
|
4996
|
495 |
|
|
496 indexed_sort.sort (vi, ns); |
|
|
497 |
|
|
498 // Flip the data out of the vector so that int compares on |
|
|
499 // IEEE754 give the correct ordering |
|
|
500 |
|
5275
|
501 for (octave_idx_type i = 0; i < ns; i++) |
|
4850
|
502 { |
|
4996
|
503 p[i*stride + offset] = IFloatFlip (vi[i]->vec); |
|
|
504 idx(i*stride + offset) = vi[i]->indx; |
|
|
505 } |
|
|
506 |
|
|
507 // There are two representations of NaN. One will be sorted |
|
|
508 // to the beginning of the vector and the other to the end. |
|
|
509 // If it will be sorted to the beginning, fix things up. |
|
2928
|
510 |
|
5009
|
511 if (lo_ieee_signbit (octave_NaN)) |
|
|
512 if (mode == UNDEFINED || mode == ASCENDING) |
|
|
513 { |
|
5275
|
514 octave_idx_type i = 0; |
|
5009
|
515 while (xisnan (v[i++*stride+offset]) && i < ns); |
|
|
516 OCTAVE_LOCAL_BUFFER (double, itmp, i - 1); |
|
5275
|
517 for (octave_idx_type l = 0; l < i -1; l++) |
|
5009
|
518 itmp[l] = idx(l*stride + offset); |
|
5275
|
519 for (octave_idx_type l = 0; l < ns - i + 1; l++) |
|
5009
|
520 { |
|
|
521 v[l*stride + offset] = v[(l+i-1)*stride + offset]; |
|
|
522 idx(l*stride + offset) = idx((l+i-1)*stride + offset); |
|
|
523 } |
|
5275
|
524 for (octave_idx_type k = 0, l = ns - i + 1; l < ns; l++, k++) |
|
5009
|
525 { |
|
|
526 v[l*stride + offset] = octave_NaN; |
|
|
527 idx(l*stride + offset) = itmp[k]; |
|
|
528 } |
|
|
529 } |
|
|
530 else |
|
|
531 { |
|
5275
|
532 octave_idx_type i = ns; |
|
5009
|
533 while (xisnan (v[--i*stride+offset]) && i > 0); |
|
|
534 OCTAVE_LOCAL_BUFFER (double, itmp, ns - i - 1); |
|
5275
|
535 for (octave_idx_type l = 0; l < ns - i -1; l++) |
|
5009
|
536 itmp[l] = idx((l+i+1)*stride + offset); |
|
5275
|
537 for (octave_idx_type l = i; l >= 0; l--) |
|
5009
|
538 { |
|
|
539 v[(l-i+ns-1)*stride + offset] = v[l*stride + offset]; |
|
|
540 idx((l-i+ns-1)*stride + offset) = idx(l*stride + offset); |
|
|
541 } |
|
5275
|
542 for (octave_idx_type k = 0, l = 0; l < ns - i - 1; l++, k++) |
|
5009
|
543 { |
|
|
544 v[l*stride + offset] = octave_NaN; |
|
|
545 idx(l*stride + offset) = itmp[k]; |
|
|
546 } |
|
|
547 } |
|
2928
|
548 } |
|
|
549 |
|
4997
|
550 retval(1) = idx; |
|
|
551 retval(0) = m; |
|
4998
|
552 |
|
2928
|
553 return retval; |
|
|
554 } |
|
|
555 |
|
4996
|
556 #else |
|
|
557 |
|
5009
|
558 template <> |
|
4996
|
559 bool |
|
|
560 ascending_compare (double a, double b) |
|
4991
|
561 { |
|
4997
|
562 return (xisnan (b) || (a < b)); |
|
4996
|
563 } |
|
|
564 |
|
5009
|
565 template <> |
|
4996
|
566 bool |
|
|
567 ascending_compare (vec_index<double> *a, vec_index<double> *b) |
|
|
568 { |
|
4997
|
569 return (xisnan (b->vec) || (a->vec < b->vec)); |
|
4996
|
570 } |
|
|
571 |
|
5009
|
572 template <> |
|
4996
|
573 bool |
|
|
574 descending_compare (double a, double b) |
|
|
575 { |
|
4997
|
576 return (xisnan (a) || (a > b)); |
|
4996
|
577 } |
|
4991
|
578 |
|
5009
|
579 template <> |
|
4991
|
580 bool |
|
4996
|
581 descending_compare (vec_index<double> *a, vec_index<double> *b) |
|
|
582 { |
|
4997
|
583 return (xisnan (a->vec) || (a->vec > b->vec)); |
|
4996
|
584 } |
|
|
585 |
|
|
586 template class octave_sort<double>; |
|
|
587 template class vec_index<double>; |
|
|
588 template class octave_sort<vec_index<double> *>; |
|
|
589 |
|
|
590 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL) |
|
|
591 static octave_value_list |
|
|
592 mx_sort (ArrayN<double> &m, int dim, sortmode mode); |
|
|
593 |
|
|
594 static octave_value_list |
|
|
595 mx_sort_indexed (ArrayN<double> &m, int dim, sortmode mode); |
|
|
596 #endif |
|
|
597 #endif |
|
|
598 |
|
|
599 // std::abs(Inf) returns NaN!! |
|
|
600 static inline double |
|
|
601 xabs (const Complex& x) |
|
|
602 { |
|
|
603 return (xisinf (x.real ()) || xisinf (x.imag ())) ? octave_Inf : abs (x); |
|
|
604 } |
|
|
605 |
|
5009
|
606 template <> |
|
4996
|
607 bool |
|
|
608 ascending_compare (vec_index<Complex> *a, vec_index<Complex> *b) |
|
|
609 { |
|
4997
|
610 return (xisnan (b->vec) |
|
|
611 || (xabs (a->vec) < xabs (b->vec)) |
|
|
612 || ((xabs (a->vec) == xabs (b->vec)) |
|
|
613 && (arg (a->vec) < arg (b->vec)))); |
|
4996
|
614 } |
|
|
615 |
|
5009
|
616 template <> |
|
4996
|
617 bool |
|
|
618 descending_compare (vec_index<Complex> *a, vec_index<Complex> *b) |
|
|
619 { |
|
4997
|
620 return (xisnan (a->vec) |
|
|
621 || (xabs (a->vec) > xabs (b->vec)) |
|
|
622 || ((xabs (a->vec) == xabs (b->vec)) |
|
|
623 && (arg (a->vec) > arg (b->vec)))); |
|
4996
|
624 } |
|
|
625 |
|
|
626 template class vec_index<Complex>; |
|
|
627 template class octave_sort<vec_index<Complex> *>; |
|
|
628 |
|
|
629 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL) |
|
|
630 static octave_value_list |
|
|
631 mx_sort_indexed (ArrayN<Complex> &m, int dim, sortmode mode); |
|
|
632 #endif |
|
|
633 |
|
4991
|
634 template class octave_sort<char>; |
|
4996
|
635 template class vec_index<char>; |
|
|
636 template class octave_sort<vec_index<char> *>; |
|
|
637 |
|
|
638 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL) |
|
5009
|
639 bool |
|
|
640 ascending_compare (char a, char b); |
|
|
641 |
|
|
642 bool |
|
|
643 ascending_compare (vec_index<char> *a, vec_index<char> *b); |
|
|
644 |
|
|
645 bool |
|
|
646 descending_compare (char a, char b); |
|
|
647 |
|
|
648 bool |
|
|
649 descending_compare (vec_index<char> *a, vec_index<char> *b); |
|
|
650 |
|
4996
|
651 static octave_value_list |
|
|
652 mx_sort (ArrayN<char> &m, int dim, sortmode mode); |
|
4991
|
653 |
|
|
654 static octave_value_list |
|
4996
|
655 mx_sort_indexed (ArrayN<char> &m, int dim, sortmode mode); |
|
|
656 #endif |
|
4991
|
657 |
|
5009
|
658 template <> |
|
4996
|
659 bool |
|
|
660 ascending_compare (vec_index<octave_value> *a, vec_index<octave_value> *b) |
|
|
661 { |
|
|
662 return (a->vec.string_value () < b->vec.string_value ()); |
|
|
663 } |
|
4991
|
664 |
|
5009
|
665 template <> |
|
4996
|
666 bool |
|
|
667 descending_compare (vec_index<octave_value> *a, vec_index<octave_value> *b) |
|
|
668 { |
|
|
669 return (a->vec.string_value () > b->vec.string_value ()); |
|
|
670 } |
|
4991
|
671 |
|
4996
|
672 template class vec_index<octave_value>; |
|
|
673 template class octave_sort<vec_index<octave_value> *>; |
|
4991
|
674 |
|
4996
|
675 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL) |
|
|
676 static octave_value_list |
|
|
677 mx_sort_indexed (ArrayN<octave_value> &m, int dim, sortmode mode); |
|
|
678 #endif |
|
4991
|
679 |
|
2928
|
680 DEFUN_DLD (sort, args, nargout, |
|
3369
|
681 "-*- texinfo -*-\n\ |
|
|
682 @deftypefn {Loadable Function} {[@var{s}, @var{i}] =} sort (@var{x})\n\ |
|
4850
|
683 @deftypefnx {Loadable Function} {[@var{s}, @var{i}] =} sort (@var{x}, @var{dim})\n\ |
|
5000
|
684 @deftypefnx {Loadable Function} {[@var{s}, @var{i}] =} sort (@var{x}, @var{mode})\n\ |
|
|
685 @deftypefnx {Loadable Function} {[@var{s}, @var{i}] =} sort (@var{x}, @var{dim}, @var{mode})\n\ |
|
5499
|
686 Return a copy of @var{x} with the elements arranged in increasing\n\ |
|
|
687 order. For matrices, @code{sort} orders the elements in each column.\n\ |
|
3369
|
688 \n\ |
|
|
689 For example,\n\ |
|
|
690 \n\ |
|
|
691 @example\n\ |
|
|
692 @group\n\ |
|
|
693 sort ([1, 2; 2, 3; 3, 1])\n\ |
|
|
694 @result{} 1 1\n\ |
|
|
695 2 2\n\ |
|
|
696 3 3\n\ |
|
|
697 @end group\n\ |
|
|
698 @end example\n\ |
|
2928
|
699 \n\ |
|
3369
|
700 The @code{sort} function may also be used to produce a matrix\n\ |
|
|
701 containing the original row indices of the elements in the sorted\n\ |
|
|
702 matrix. For example,\n\ |
|
|
703 \n\ |
|
|
704 @example\n\ |
|
|
705 @group\n\ |
|
|
706 [s, i] = sort ([1, 2; 2, 3; 3, 1])\n\ |
|
|
707 @result{} s = 1 1\n\ |
|
|
708 2 2\n\ |
|
|
709 3 3\n\ |
|
|
710 @result{} i = 1 3\n\ |
|
|
711 2 1\n\ |
|
|
712 3 2\n\ |
|
|
713 @end group\n\ |
|
|
714 @end example\n\ |
|
4850
|
715 \n\ |
|
|
716 If the optional argument @var{dim} is given, then the matrix is sorted\n\ |
|
4996
|
717 along the dimension defined by @var{dim}. The optional argument @code{mode}\n\ |
|
|
718 defines the order in which the values will be sorted. Valid values of\n\ |
|
|
719 @code{mode} are `ascend' or `descend'.\n\ |
|
4850
|
720 \n\ |
|
|
721 For equal elements, the indices are such that the equal elements are listed\n\ |
|
|
722 in the order that appeared in the original list.\n\ |
|
|
723 \n\ |
|
4996
|
724 The @code{sort} function may also be used to sort strings and cell arrays\n\ |
|
5499
|
725 of strings, in which case the dictionary order of the strings is used.\n\ |
|
4996
|
726 \n\ |
|
4850
|
727 The algorithm used in @code{sort} is optimized for the sorting of partially\n\ |
|
|
728 ordered lists.\n\ |
|
3369
|
729 @end deftypefn") |
|
2928
|
730 { |
|
|
731 octave_value_list retval; |
|
|
732 |
|
|
733 int nargin = args.length (); |
|
4996
|
734 sortmode smode = ASCENDING; |
|
2928
|
735 |
|
5001
|
736 if (nargin < 1 || nargin > 3) |
|
2928
|
737 { |
|
|
738 print_usage ("sort"); |
|
|
739 return retval; |
|
|
740 } |
|
|
741 |
|
4850
|
742 bool return_idx = nargout > 1; |
|
2928
|
743 |
|
|
744 octave_value arg = args(0); |
|
|
745 |
|
4850
|
746 int dim = 0; |
|
4996
|
747 if (nargin > 1) |
|
|
748 { |
|
|
749 if (args(1).is_string ()) |
|
|
750 { |
|
|
751 std::string mode = args(1).string_value(); |
|
|
752 if (mode == "ascend") |
|
|
753 smode = ASCENDING; |
|
|
754 else if (mode == "descend") |
|
|
755 smode = DESCENDING; |
|
|
756 else |
|
|
757 { |
|
5000
|
758 error ("sort: mode must be either \"ascend\" or \"descend\""); |
|
4996
|
759 return retval; |
|
|
760 } |
|
|
761 } |
|
|
762 else |
|
|
763 dim = args(1).nint_value () - 1; |
|
|
764 } |
|
|
765 |
|
|
766 if (nargin > 2) |
|
|
767 { |
|
|
768 if (args(1).is_string ()) |
|
|
769 { |
|
|
770 print_usage ("sort"); |
|
|
771 return retval; |
|
|
772 } |
|
|
773 |
|
|
774 if (! args(2).is_string ()) |
|
|
775 { |
|
|
776 error ("sort: mode must be a string"); |
|
|
777 return retval; |
|
|
778 } |
|
|
779 std::string mode = args(2).string_value(); |
|
|
780 if (mode == "ascend") |
|
|
781 smode = ASCENDING; |
|
|
782 else if (mode == "descend") |
|
|
783 smode = DESCENDING; |
|
|
784 else |
|
|
785 { |
|
5000
|
786 error ("sort: mode must be either \"ascend\" or \"descend\""); |
|
4996
|
787 return retval; |
|
|
788 } |
|
|
789 } |
|
4850
|
790 |
|
5275
|
791 dim_vector dv = arg.dims (); |
|
4850
|
792 if (error_state) |
|
|
793 { |
|
|
794 gripe_wrong_type_arg ("sort", arg); |
|
|
795 return retval; |
|
|
796 } |
|
4996
|
797 if (nargin == 1 || args(1).is_string ()) |
|
4850
|
798 { |
|
|
799 // Find first non singleton dimension |
|
|
800 for (int i = 0; i < dv.length (); i++) |
|
|
801 if (dv(i) > 1) |
|
|
802 { |
|
|
803 dim = i; |
|
|
804 break; |
|
|
805 } |
|
|
806 } |
|
|
807 else |
|
|
808 { |
|
|
809 if (dim < 0 || dim > dv.length () - 1) |
|
|
810 { |
|
|
811 error ("sort: dim must be a valid dimension"); |
|
|
812 return retval; |
|
|
813 } |
|
|
814 } |
|
|
815 |
|
2928
|
816 if (arg.is_real_type ()) |
|
|
817 { |
|
4850
|
818 NDArray m = arg.array_value (); |
|
2928
|
819 |
|
|
820 if (! error_state) |
|
4996
|
821 { |
|
|
822 #ifdef HAVE_IEEE754_DATA_FORMAT |
|
|
823 // As operator > gives the right result, can special case here |
|
|
824 if (! return_idx && smode == ASCENDING) |
|
|
825 retval = mx_sort (m, dim); |
|
|
826 else |
|
|
827 #endif |
|
|
828 { |
|
|
829 if (return_idx) |
|
|
830 retval = mx_sort_indexed (m, dim, smode); |
|
|
831 else |
|
|
832 retval = mx_sort (m, dim, smode); |
|
|
833 } |
|
|
834 } |
|
2928
|
835 } |
|
|
836 else if (arg.is_complex_type ()) |
|
|
837 { |
|
4850
|
838 ComplexNDArray cm = arg.complex_array_value (); |
|
2928
|
839 |
|
4996
|
840 // Don't have unindexed version as no ">" operator |
|
2928
|
841 if (! error_state) |
|
4996
|
842 retval = mx_sort_indexed (cm, dim, smode); |
|
2928
|
843 } |
|
4991
|
844 else if (arg.is_string ()) |
|
|
845 { |
|
|
846 charNDArray chm = arg.char_array_value (); |
|
|
847 |
|
|
848 if (! error_state) |
|
4996
|
849 { |
|
|
850 // As operator > gives the right result, can special case here |
|
|
851 if (! return_idx && smode == ASCENDING) |
|
|
852 retval = mx_sort (chm, dim); |
|
|
853 else |
|
|
854 { |
|
|
855 if (return_idx) |
|
|
856 retval = mx_sort_indexed (chm, dim, smode); |
|
|
857 else |
|
|
858 retval = mx_sort (chm, dim, smode); |
|
|
859 } |
|
|
860 |
|
|
861 // XXX FIXME XXX It would have been better to call |
|
|
862 // "octave_value(m, true)" but how can that be done |
|
|
863 // within the template |
|
|
864 retval(0) = retval(0).convert_to_str (false, true); |
|
|
865 } |
|
|
866 } |
|
|
867 else if (arg.is_cell ()) |
|
|
868 { |
|
|
869 Cell cellm = arg.cell_value (); |
|
|
870 |
|
|
871 // Need to check that all elements are strings |
|
5275
|
872 for (octave_idx_type i = 0; i < cellm.numel (); i++) |
|
4996
|
873 if (! cellm(i).is_string ()) |
|
|
874 { |
|
|
875 gripe_wrong_type_arg ("sort", arg); |
|
|
876 break; |
|
|
877 } |
|
|
878 |
|
|
879 // Don't have unindexed version as ">" operator doesn't return bool |
|
|
880 if (!error_state) |
|
|
881 retval = mx_sort_indexed (cellm, dim, smode); |
|
4991
|
882 } |
|
2928
|
883 else |
|
|
884 gripe_wrong_type_arg ("sort", arg); |
|
|
885 |
|
|
886 return retval; |
|
|
887 } |
|
|
888 |
|
|
889 /* |
|
|
890 ;;; Local Variables: *** |
|
|
891 ;;; mode: C++ *** |
|
|
892 ;;; End: *** |
|
|
893 */ |
