# Ublas :

Subject: Re: [ublas] sorting a ublas vector and get index
From: Gunter Winkler (guwi17_at_[hidden])
Date: 2009-12-07 17:37:02

Am Wednesday 02 December 2009 schrieb Kaveh Kohan:
> Hi All,
>
> I have a simple question but it seems that the answer is not that
> trivial not at least to me. Could you please tell me how to sort a
> vector and get corresponding indexes? Let me give and example from
> MATLAB:
>
> [a,b]=sort([5, 8, 7])
>
> a = 5 7 8
> b = 1 3 2
>
> a is the sorted vector and b is corresponding indexes. It seems that
> the same functionality exist in gsl:
>

This feature is used in the coordinate_matrix::sort() method. There is a
(experimental) class "index_pair_array" which can be used to generate
proxy classes.

Alternatively, you create a vector holding the indices:

VECTOR v(...);
vector<int> idx(v.size());
// set all v[i] = i; // or use ublas::permutation_matrix
std::stable_sort(idx.begin(), idx.end(), vector_less< VECTOR >(v));

using the class

template <class VEC>
class vector_less
{
private:
typedef typename VEC::size_type size_type;
typedef typename VEC::value_type value_type;
vector_less();

const VEC& data;
public:
vector_less(const VEC& vec) : data(vec) { }
bool operator() (const size_type& left, const size_type& right)
const
{
return std::less<value_type> () ( data(left), data(right) );
}
};

to cite some more code, here is part of a routine sorting elements of a
mesh:

template <class VEC>
VEC inverse_permutation(VEC & p)
{
typedef typename VEC::size_type size_type;
VEC res(p.size());

for (size_type i=0; i<p.size(); ++i) {
res(p(i)) = i;
}

return res;
}

...

/// sort nodes by node level
void sort_node_level()
{
ublas::permutation_matrix<size_type> P( node_count() );

std::stable_sort(P.begin(), P.end(), vector_less< node_level_type >
( node_level_vector ) );

#ifndef NDEBUG
std::cout << P << std::endl;
#endif

ublas::permutation_matrix<size_type> IP = inverse_permutation(P);

#ifndef NDEBUG
std::cout << IP << std::endl;
#endif

for (size_type i=0; i<edge_count(); ++i) {
edgematrix(i,0) = IP(edgematrix(i,0));
edgematrix(i,1) = IP(edgematrix(i,1));
}

for (size_type i=0; i<face_count(); ++i) {
for (size_type j=0; j<facematrix.size2(); ++j) {
facematrix(i,j) = IP(facematrix(i,j));
}
}

#ifndef NDEBUG
ublas::permutation_matrix<size_type> P2(P);
coordmatrix_type C2(coordmatrix);
coordmatrix_type C3(coordmatrix.size());
for (size_type i=0; i<P2.size(); ++i) {
for (size_type j=0; j < dimension; ++j) {
C3[i][j] = C2[P2(i)][j];
}
}
#endif

permutation_to_swap(P);

#ifndef NDEBUG
std::cout << P << std::endl;
#endif

boost::numeric::ublas::swap_rows(P, node_level_vector);
swap_rows(P, coordmatrix);
boost::numeric::ublas::swap_rows(P, coorddata);

#ifndef NDEBUG
std::cout << "# check perm" << std::endl;
for (size_type i=0; i<C3.size(); ++i) {
if ((coordmatrix[i][0] != C3[i][0])
|| (coordmatrix[i][1] != C3[i][1])) {
std::cout << i << ": " << C3[i] << std::endl;
}
}
#endif
...

mfg
Gunter