/*
Copyright (C) 2003 David Bateman
This program is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2, or (at your option) any
later version.
This program is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with this program; see the file COPYING. If not, write to the Free
Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
In addition to the terms of the GPL, you are permitted to link
this program with any Open Source program, as defined by the
Open Source Initiative (www.opensource.org)
*/
#include <octave/oct.h>
#include <octave/lo-mappers.h>
#include <octave/config.h>
#include <octave/quit.h>
#include <octave/lo-ieee.h>
#include <octave/data-conv.h>
#include <octave/ov-cx-mat.h>
#include "oct-sort.cc"
// ======= Cruft to support ancient versions of Octave =========
#ifndef OCTAVE_LOCAL_BUFFER
#include <vector>
#define OCTAVE_LOCAL_BUFFER(T, buf, size) \
std::vector<T> buf ## _vector (size); \
T *buf = &(buf ## _vector[0])
#endif
#ifndef lo_ieee_signbit
#if defined (signbit)
#define lo_ieee_signbit(x) signbit (x)
#elif defined (HAVE_SIGNBIT)
#if defined (__MINGW32__)
extern int signbit (double);
#endif
#define lo_ieee_signbit(x) (x < 0 || signbit (x))
#elif defined (copysign)
#define lo_ieee_signbit(x) (copysign (1.0, x) < 0)
#elif defined (HAVE_COPYSIGN)
#define lo_ieee_signbit(x) (x < 0 || copysign (1.0, x) < 0)
#else
#define lo_ieee_signbit(x) 0
#endif
#endif
// ========= End Cruft ===================================
/* If we are IEEE 754 or IEEE 854 compliant, then we can use the trick of
* casting doubles as unsigned eight byte integers, and with a little
* bit of magic we can automatically sort the NaN's correctly.
*/
#if defined(HAVE_IEEE754_COMPLIANCE) && defined(EIGHT_BYTE_INT)
static inline unsigned EIGHT_BYTE_INT FloatFlip(unsigned EIGHT_BYTE_INT f)
{
unsigned EIGHT_BYTE_INT mask = -(EIGHT_BYTE_INT)(f >> 63) |
0x8000000000000000ULL;
return f ^ mask;
}
inline unsigned EIGHT_BYTE_INT IFloatFlip(unsigned EIGHT_BYTE_INT f)
{
unsigned EIGHT_BYTE_INT mask = ((f >> 63) - 1) | 0x8000000000000000ULL;
return f ^ mask;
}
struct vec_index
{
unsigned EIGHT_BYTE_INT vec;
int indx;
};
bool
ieee754_compare (vec_index *a, vec_index *b)
{
return (a->vec < b->vec);
}
template class octave_sort<unsigned EIGHT_BYTE_INT>;
template class octave_sort<vec_index *>;
#else
struct vec_index
{
double vec;
int indx;
};
bool
double_compare (double a, double b)
{
return (xisnan(b) || (a < b));
}
bool
double_compare (vec_index *a, vec_index *b)
{
return (xisnan(b->vec) || (a->vec < b->vec));
}
template class octave_sort<double>;
template class octave_sort<vec_index *>;
#endif
struct complex_vec_index
{
Complex vec;
int indx;
};
bool
complex_compare (complex_vec_index *a, complex_vec_index *b)
{
return (xisnan(b->vec) || (abs(a->vec) < abs(b->vec)));
}
template class octave_sort<complex_vec_index *>;
#ifdef HAVE_ND_ARRAYS
static octave_value_list
mx_sort (NDArray &m, bool return_idx, int dim)
{
octave_value_list retval;
if (m.length () < 1)
return retval;
dim_vector dv = m.dims ();
unsigned int ns = dv (dim);
unsigned int iter = dv.numel () / ns;
unsigned int stride = 1;
for (unsigned int i = 0; i < (unsigned int)dim; i++)
stride *= dv(i);
#if defined(HAVE_IEEE754_COMPLIANCE) && defined(EIGHT_BYTE_INT)
double *v = m.fortran_vec ();
unsigned EIGHT_BYTE_INT *p = (unsigned EIGHT_BYTE_INT *)v;
if (return_idx)
{
octave_sort<vec_index *> indexed_ieee754_sort (ieee754_compare);
OCTAVE_LOCAL_BUFFER (vec_index *, vi, ns);
OCTAVE_LOCAL_BUFFER (vec_index, vix, ns);
for (unsigned int i = 0; i < ns; i++)
vi[i] = &vix[i];
NDArray idx (dv);
for (unsigned int j = 0; j < iter; j++)
{
unsigned int offset = j;
unsigned int offset2 = 0;
while (offset >= stride)
{
offset -= stride;
offset2++;
}
offset += offset2 * stride * ns;
/* Flip the data in the vector so that int compares on
* IEEE754 give the correct ordering
*/
for (unsigned int i = 0; i < ns; i++)
{
vi[i]->vec = FloatFlip (p[i*stride + offset]);
vi[i]->indx = i + 1;
}
indexed_ieee754_sort.sort (vi, ns);
/* Flip the data out of the vector so that int compares on
* IEEE754 give the correct ordering
*/
for (unsigned int i = 0; i < ns; i++)
{
p[i*stride + offset] = IFloatFlip (vi[i]->vec);
idx(i*stride + offset) = vi[i]->indx;
}
/* There are two representations of NaN. One will be sorted to
* the beginning of the vector and the other to the end. If it
* will be sorted to the beginning, fix things up.
*/
if (lo_ieee_signbit (octave_NaN))
{
unsigned int i = 0;
while (xisnan(v[i++*stride+offset]) && i < ns);
OCTAVE_LOCAL_BUFFER (double, itmp, i - 1);
for (unsigned int l = 0; l < i -1; l++)
itmp[l] = idx(l*stride + offset);
for (unsigned int l = 0; l < ns - i + 1; l++)
{
v[l*stride + offset] = v[(l+i-1)*stride + offset];
idx(l*stride + offset) = idx((l+i-1)*stride + offset);
}
for (unsigned int k = 0, l = ns - i + 1; l < ns; l++, k++)
{
v[l*stride + offset] = octave_NaN;
idx(l*stride + offset) = itmp[k];
}
}
}
retval (1) = idx;
}
else
{
octave_sort<unsigned EIGHT_BYTE_INT> ieee754_sort;
if (stride == 1)
{
for (unsigned int j = 0; j < iter; j++)
{
/* Flip the data in the vector so that int compares on
* IEEE754 give the correct ordering
*/
for (unsigned int i = 0; i < ns; i++)
p[i] = FloatFlip (p[i]);
ieee754_sort.sort (p, ns);
/* Flip the data out of the vector so that int compares on
* IEEE754 give the correct ordering
*/
for (unsigned int i = 0; i < ns; i++)
p[i] = IFloatFlip (p[i]);
/* There are two representations of NaN. One will be sorted to
* the beginning of the vector and the other to the end. If it
* will be sorted to the beginning, fix things up.
*/
if (lo_ieee_signbit (octave_NaN))
{
unsigned int i = 0;
double *vtmp = (double *)p;
while (xisnan(vtmp[i++]) && i < ns);
for (unsigned int l = 0; l < ns - i + 1; l++)
vtmp[l] = vtmp[l+i-1];
for (unsigned int l = ns - i + 1; l < ns; l++)
vtmp[l] = octave_NaN;
}
p += ns;
}
}
else
{
OCTAVE_LOCAL_BUFFER (unsigned EIGHT_BYTE_INT, vi, ns);
for (unsigned int j = 0; j < iter; j++)
{
unsigned int offset = j;
unsigned int offset2 = 0;
while (offset >= stride)
{
offset -= stride;
offset2++;
}
offset += offset2 * stride * ns;
/* Flip the data in the vector so that int compares on
* IEEE754 give the correct ordering
*/
for (unsigned int i = 0; i < ns; i++)
vi[i] = FloatFlip (p[i*stride + offset]);
ieee754_sort.sort (vi, ns);
/* Flip the data out of the vector so that int compares on
* IEEE754 give the correct ordering
*/
for (unsigned int i = 0; i < ns; i++)
p[i*stride + offset] = IFloatFlip (vi[i]);
/* There are two representations of NaN. One will be sorted to
* the beginning of the vector and the other to the end. If it
* will be sorted to the beginning, fix things up.
*/
if (lo_ieee_signbit (octave_NaN))
{
unsigned int i = 0;
while (xisnan(v[i++*stride + offset]) && i < ns);
for (unsigned int l = 0; l < ns - i + 1; l++)
v[l*stride + offset] = v[(l+i-1)*stride + offset];
for (unsigned int l = ns - i + 1; l < ns; l++)
v[l*stride + offset] = octave_NaN;
}
}
}
}
#else
if (return_idx)
{
double *v = m.fortran_vec ();
octave_sort<vec_index *> indexed_double_sort (double_compare);
OCTAVE_LOCAL_BUFFER (vec_index *, vi, ns);
OCTAVE_LOCAL_BUFFER (vec_index, vix, ns);
for (unsigned int i = 0; i < ns; i++)
vi[i] = &vix[i];
NDArray idx (dv);
if (stride == 1)
{
for (unsigned int j = 0; j < iter; j++)
{
unsigned int offset = j * ns;
for (unsigned int i = 0; i < ns; i++)
{
vi[i]->vec = v[i];
vi[i]->indx = i + 1;
}
indexed_double_sort.sort (vi, ns);
for (unsigned int i = 0; i < ns; i++)
{
v[i] = vi[i]->vec;
idx(i + offset) = vi[i]->indx;
}
v += ns;
}
}
else
{
for (unsigned int j = 0; j < iter; j++)
{
unsigned int offset = j;
unsigned int offset2 = 0;
while (offset >= stride)
{
offset -= stride;
offset2++;
}
offset += offset2 * stride * ns;
for (unsigned int i = 0; i < ns; i++)
{
vi[i]->vec = v[i*stride + offset];
vi[i]->indx = i + 1;
}
indexed_double_sort.sort (vi, ns);
for (unsigned int i = 0; i < ns; i++)
{
v[i*stride+offset] = vi[i]->vec;
idx(i*stride+offset) = vi[i]->indx;
}
}
}
retval (1) = idx;
}
else
{
double *v = m.fortran_vec ();
octave_sort<double> double_sort (double_compare);
if (stride == 1)
for (unsigned int j = 0; j < iter; j++)
{
double_sort.sort (v, ns);
v += ns;
}
else
{
OCTAVE_LOCAL_BUFFER (double, vi, ns);
for (unsigned int j = 0; j < iter; j++)
{
unsigned int offset = j;
unsigned int offset2 = 0;
while (offset >= stride)
{
offset -= stride;
offset2++;
}
offset += offset2 * stride * ns;
for (unsigned int i = 0; i < ns; i++)
vi[i] = v[i*stride + offset];
double_sort.sort (vi, ns);
for (unsigned int i = 0; i < ns; i++)
v[i*stride + offset] = vi[i];
}
}
}
#endif
retval(0) = m;
return retval;
}
static octave_value_list
mx_sort (ComplexNDArray &m, bool return_idx, int dim)
{
octave_value_list retval;
if (m.length () < 1)
return retval;
dim_vector dv = m.dims ();
unsigned int ns = dv (dim);
unsigned int iter = dv.numel () / ns;
unsigned int stride = 1;
for (unsigned int i = 0; i < (unsigned int)dim; i++)
stride *= dv(i);
octave_sort<complex_vec_index *> indexed_double_sort (complex_compare);
Complex *v = m.fortran_vec ();
OCTAVE_LOCAL_BUFFER (complex_vec_index *, vi, ns);
OCTAVE_LOCAL_BUFFER (complex_vec_index, vix, ns);
for (unsigned int i = 0; i < ns; i++)
vi[i] = &vix[i];
NDArray idx (dv);
if (stride == 1)
{
for (unsigned int j = 0; j < iter; j++)
{
unsigned int offset = j * ns;
for (unsigned int i = 0; i < ns; i++)
{
vi[i]->vec = v[i];
vi[i]->indx = i + 1;
}
indexed_double_sort.sort (vi, ns);
if (return_idx)
{
for (unsigned int i = 0; i < ns; i++)
{
v[i] = vi[i]->vec;
idx(i + offset) = vi[i]->indx;
}
}
else
{
for (unsigned int i = 0; i < ns; i++)
v[i] = vi[i]->vec;
}
v += ns;
}
}
else
{
for (unsigned int j = 0; j < iter; j++)
{
unsigned int offset = j;
unsigned int offset2 = 0;
while (offset >= stride)
{
offset -= stride;
offset2++;
}
offset += offset2 * stride * ns;
for (unsigned int i = 0; i < ns; i++)
{
vi[i]->vec = v[i*stride + offset];
vi[i]->indx = i + 1;
}
indexed_double_sort.sort (vi, ns);
if (return_idx)
{
for (unsigned int i = 0; i < ns; i++)
{
v[i*stride + offset] = vi[i]->vec;
idx(i*stride + offset) = vi[i]->indx;
}
}
else
{
for (unsigned int i = 0; i < ns; i++)
v[i*stride + offset] = vi[i]->vec;
}
}
}
if (return_idx)
retval (1) = idx;
retval(0) = m;
return retval;
}
#else
static octave_value_list
mx_sort (Matrix &m, bool return_idx, int dim)
{
octave_value_list retval;
int nr = m.rows ();
int nc = m.columns ();
bool trans = false;
int nc2, nr2;
if (m.length () < 1)
return retval;
if (dim)
{
nr2 = nc;
nc2 = nr;
if (nr != 1)
{
trans = true;
m = m.transpose ();
}
}
else
{
nr2 = nr;
nc2 = nc;
}
#if defined(HAVE_IEEE754_COMPLIANCE) && defined(EIGHT_BYTE_INT)
double *v = m.fortran_vec ();
unsigned EIGHT_BYTE_INT *p = (unsigned EIGHT_BYTE_INT *)v;
if (return_idx)
{
octave_sort<vec_index *> indexed_ieee754_sort (ieee754_compare);
OCTAVE_LOCAL_BUFFER (vec_index *, vi, nr2);
OCTAVE_LOCAL_BUFFER (vec_index, vix, nr2);
for (int i = 0; i < nr2; i++)
vi[i] = &vix[i];
Matrix idx (nr2, nc2);
for (int j = 0; j < nc2; j++)
{
/* Flip the data in the vector so that int compares on
* IEEE754 give the correct ordering
*/
for (int i = 0; i < nr2; i++)
{
vi[i]->vec = FloatFlip (p[i]);
vi[i]->indx = i + 1;
}
indexed_ieee754_sort.sort (vi, nr2);
/* Flip the data out of the vector so that int compares on
* IEEE754 give the correct ordering
*/
for (int i = 0; i < nr2; i++)
{
p[i] = IFloatFlip (vi[i]->vec);
idx(i,j) = vi[i]->indx;
}
/* There are two representations of NaN. One will be sorted to
* the beginning of the vector and the other to the end. If it
* will be sorted to the beginning, fix things up.
*/
if (lo_ieee_signbit (octave_NaN))
{
unsigned int i = 0;
double *vtmp = (double *)p;
while (xisnan(vtmp[i++]));
OCTAVE_LOCAL_BUFFER (double, itmp, i - 1);
for (unsigned int l = 0; l < i -1; l++)
itmp[l] = idx(l,j);
for (unsigned int l = 0; l < nr2 - i + 1; l++)
{
vtmp[l] = vtmp[l+i-1];
idx(l,j) = idx(l+i-1,j);
}
for (unsigned int k = 0, l = nr2 - i + 1; l < nr2;
l++, k++)
{
vtmp[l] = octave_NaN;
idx(l,j) = itmp[k];
}
}
p += nr2;
}
if (dim)
retval (1) = idx.transpose ();
else
retval (1) = idx;
}
else
{
octave_sort<unsigned EIGHT_BYTE_INT> ieee754_sort;
for (int j = 0; j < nc2; j++)
{
/* Flip the data in the vector so that int compares on
* IEEE754 give the correct ordering
*/
for (int i = 0; i < nr2; i++)
p[i] = FloatFlip (p[i]);
ieee754_sort.sort (p, nr2);
/* Flip the data out of the vector so that int compares on
* IEEE754 give the correct ordering
*/
for (int i = 0; i < nr2; i++)
p[i] = IFloatFlip (p[i]);
/* There are two representations of NaN. One will be sorted to
* the beginning of the vector and the other to the end. If it
* will be sorted to the beginning, fix things up.
*/
if (lo_ieee_signbit (octave_NaN))
{
unsigned int i = 0;
double *vtmp = (double *)p;
while (xisnan(vtmp[i++]));
for (unsigned int l = 0; l < nr2 - i + 1; l++)
vtmp[l] = vtmp[l+i-1];
for (unsigned int l = nr2 - i + 1; l < nr2; l++)
vtmp[l] = octave_NaN;
}
p += nr2;
}
}
#else
if (return_idx)
{
double *v = m.fortran_vec ();
octave_sort<vec_index *> indexed_double_sort (double_compare);
OCTAVE_LOCAL_BUFFER (vec_index *, vi, nr2);
OCTAVE_LOCAL_BUFFER (vec_index, vix, nr2);
for (int i = 0; i < nr2; i++)
vi[i] = &vix[i];
Matrix idx (nr2, nc2);
for (int j = 0; j < nc2; j++)
{
for (int i = 0; i < nr2; i++)
{
vi[i]->vec = v[i];
vi[i]->indx = i + 1;
}
indexed_double_sort.sort (vi, nr2);
for (int i = 0; i < nr2; i++)
{
v[i] = vi[i]->vec;
idx(i,j) = vi[i]->indx;
}
v += nr2;
}
if (dim)
retval (1) = idx.transpose ();
else
retval (1) = idx;
}
else
{
double *v = m.fortran_vec ();
octave_sort<double> double_sort (double_compare);
for (int j = 0; j < nc2; j++)
{
double_sort.sort (v, nr2);
v += nr2;
}
}
#endif
if (trans)
retval(0) = m.transpose ();
else
retval(0) = m;
return retval;
}
// I don't think there is any sane way of doing this with IEEE754 integer
// compares. As the operation is fundamentally without sense, I don't
// care that its not optimized....
static octave_value_list
mx_sort (ComplexMatrix &m, bool return_idx, int dim)
{
octave_value_list retval;
int nr = m.rows ();
int nc = m.columns ();
bool trans = false;
int nc2, nr2;
if (m.length () < 1)
return retval;
if (dim)
{
nr2 = nc;
nc2 = nr;
if (nr != 1)
{
trans = true;
m = m.transpose ();
}
}
else
{
nr2 = nr;
nc2 = nc;
}
octave_sort<complex_vec_index *> indexed_double_sort (complex_compare);
Complex *v = m.fortran_vec ();
OCTAVE_LOCAL_BUFFER (complex_vec_index *, vi, nr2);
OCTAVE_LOCAL_BUFFER (complex_vec_index, vix, nr2);
for (int i = 0; i < nr2; i++)
vi[i] = &vix[i];
Matrix idx (nr2, nc2);
for (int j = 0; j < nc2; j++)
{
for (int i = 0; i < nr2; i++)
{
vi[i]->vec = v[i];
vi[i]->indx = i + 1;
}
indexed_double_sort.sort (vi, nr2);
if (return_idx)
{
for (int i = 0; i < nr2; i++)
{
v[i] = vi[i]->vec;
idx(i,j) = vi[i]->indx;
}
}
else
{
for (int i = 0; i < nr2; i++)
v[i] = vi[i]->vec;
}
v += nr2;
}
if (return_idx)
if (dim)
retval (1) = idx.transpose ();
else
retval (1) = idx;
if (trans)
retval(0) = m.transpose ();
else
retval(0) = m;
return retval;
}
#endif
DEFUN_DLD (sort, args, nargout,
"-*- texinfo -*-\n\
@deftypefn {Loadable Function} {[@var{s}, @var{i}] =} sort (@var{x})\n\
@deftypefnx {Loadable Function} {[@var{s}, @var{i}] =} sort (@var{x}, @var{dim})\n\
Return a copy of @var{x} with the elements elements arranged in\n\
increasing order. For matrices, @code{sort} orders the elements in each\n\
column.\n\
\n\
For example,\n\
\n\
@example\n\
@group\n\
sort ([1, 2; 2, 3; 3, 1])\n\
@result{} 1 1\n\
2 2\n\
3 3\n\
@end group\n\
@end example\n\
\n\
The @code{sort} function may also be used to produce a matrix\n\
containing the original row indices of the elements in the sorted\n\
matrix. For example,\n\
\n\
@example\n\
@group\n\
[s, i] = sort ([1, 2; 2, 3; 3, 1])\n\
@result{} s = 1 1\n\
2 2\n\
3 3\n\
@result{} i = 1 3\n\
2 1\n\
3 2\n\
@end group\n\
@end example\n\
\n\
If the optional argument @var{dim} is given, then the matrix is sorted\n\
along the dimension defined by @var{dim}.\n\
\n\
For equal elements, the indices are such that the equal elements\n\
are listed in the order that appeared in the original list.\n\
\n\
The algorithm used in @code{sort} is optimized for the sorting of partially\n\
ordered lists.\n\
@end deftypefn")
{
octave_value_list retval;
int nargin = args.length ();
if (nargin != 1 && nargin != 2)
{
print_usage ("sort");
return retval;
}
bool return_idx = nargout > 1;
octave_value arg = args(0);
int dim = 0;
if (nargin == 2)
dim = args(1).nint_value () - 1;
#ifdef HAVE_ND_ARRAYS
dim_vector dv = ((const octave_complex_matrix&) arg) .dims ();
if (error_state)
{
gripe_wrong_type_arg ("sort", arg);
return retval;
}
if (nargin != 2)
{
// Find first non singleton dimension
for (int i = 0; i < dv.length (); i++)
if (dv(i) > 1)
{
dim = i;
break;
}
}
else
{
if (dim < 0 || dim > dv.length () - 1)
{
error ("sort: dim must be a valid dimension");
return retval;
}
}
#else
int nr = ((const octave_complex_matrix&) arg) .rows ();
if (error_state)
{
gripe_wrong_type_arg ("sort", arg);
return retval;
}
if (nargin != 2)
dim = (nr == 1 ? 1 : 0);
else
{
if (dim < 0 || dim > 1)
{
error ("sort: dim must be a valid dimension");
return retval;
}
}
#endif
if (arg.is_real_type ())
{
#ifdef HAVE_ND_ARRAYS
NDArray m = arg.array_value ();
#else
Matrix m = arg.matrix_value ();
#endif
if (! error_state)
retval = mx_sort (m, return_idx, dim);
}
else if (arg.is_complex_type ())
{
#ifdef HAVE_ND_ARRAYS
ComplexNDArray cm = arg.complex_array_value ();
#else
ComplexMatrix cm = arg.complex_matrix_value ();
#endif
if (! error_state)
retval = mx_sort (cm, return_idx, dim);
}
else
gripe_wrong_type_arg ("sort", arg);
return retval;
}
/*
;;; Local Variables: ***
;;; mode: C++ ***
;;; End: ***
*/
syntax highlighted by Code2HTML, v. 0.9.1