Re: [ublas] sorting a ublas vector and get index

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