Examples - How To Extend UBLAS

Example 1 - "Vectorize" standard math functions

Q: Did you ever want to apply a simple function to all elements of a vector or a matrix? Usually this means iteration over all elements. But is there a better way? Can we hide the loops and even use lazy computation?

A: Yes, we can.

Example: We want to apply log() to all elements of a vector.

The logarithm is a unary function, so we have to look for vector_unary. This is already implemented for the unary minus and the complex conjugate of a vector. In order to have a flexible notation we want to create a new function similar to conj() that returns a vector expression applying a (at compile time) known functor to any element of a vector. Therefore we copy the unary minus from vector_expression.hpp and modify some types.

  // (op v) [i] = op( v [i] )
  template<class OP, class E> 
  BOOST_UBLAS_INLINE
  typename vector_unary_traits<E, OP>::result_type
  apply_to_all (const vector_expression<E> &e, const OP& op = OP() ) {
    typedef typename vector_unary_traits<E, OP>::expression_type expression_type;
    return expression_type (e ());
  }

Now we are able to write apply_to_all(v, op) where op can be any functor with following structure:

namespace functor {

  template <class T>
  class log
  {
  public:
    typedef T value_type;
    typedef T result_type;
    log() { }
    
    static
    result_type apply(const value_type& x)
    {
      return std::log(x); // insert any function you want
    }

  };

}

All done. The usage is very simple:

  vector<double> v( 100 );
  // ... fill v ...                                     
  std::cout << apply_to_all( v, functor::log<double>() ) << std::endl; 
  std::cout << apply_to_all<functor::log<double> >( v ) << std::endl; 

Example 2 - Compute a scalar from a vector

Example: You need a special vector to scalar mapping, which should work optimal for dense, packed and sparse vectors. Here you want to compute the minimal logarithm of all vector elements exclusive the zeros (which would give minus infinity) and negative values (for which the logarithm is not defined).

The requested function is similar to sum(). Therefore we use the vector_sum as a template and generalize the interface. We have to provide three versions of apply, one general, one for dense iterators and one for general (forward-)iterators.

  // pitfalls: 
  //   the result is OP::initial_value() for "empty" sparse vectors,
  //   zero elements that are not stored may be ignored
  template<class T, class OP>
  struct general_vector_scalar_unary: 
    public vector_scalar_unary_functor<T> {
    typedef typename vector_scalar_unary_functor<T>::size_type size_type;
    typedef typename vector_scalar_unary_functor<T>::difference_type difference_type;
    typedef typename vector_scalar_unary_functor<T>::value_type value_type;
    typedef typename vector_scalar_unary_functor<T>::result_type result_type;

    // general case, access by index
    template<class E>
    static BOOST_UBLAS_INLINE
    result_type apply (const vector_expression<E> &e) { 
      result_type t = OP::initial_value();
      size_type size (e ().size ());
      for (size_type i = 0; i < size; ++ i)
        OP::update(t, e () (i));
      return t;
    }
    // Dense case
    template<class I>
    static BOOST_UBLAS_INLINE
    result_type apply (difference_type size, I it) { 
      result_type t = OP::initial_value();
      while (-- size >= 0)
        OP::update(t, *it), ++ it;
      return t; 
    }
    // Sparse case
    template<class I>
    static BOOST_UBLAS_INLINE
    result_type apply (I it, const I &it_end) {
      result_type t = OP::initial_value();
      while (it != it_end) 
        OP::update(t, *it), ++ it;
      return t; 
    }
  };

Now we have the general framework that uses some functor OP. The implementation of this functor is easy, because we only have to provide two static functions.

namespace functor {

  template <class T>
  struct log_minimum
  {
  public:
    typedef T value_type;
    typedef T result_type;
    
    static
    result_type initial_value()
    {
      return std::log(std::numeric_limits<result_type>::max());
    }
    static
    void update(result_type& t, const value_type& x)
    {
      if ( (x>0) && (std::log(x)<t) ) t = std::log(x);
    }
  };

}

In order to have a convenient syntax we write a free function similar to sum that chooses the correct types.

  // log_min v = min ( log(v [i]), i=1..N if v[i]>0 )
  template<class E>
  BOOST_UBLAS_INLINE
  typename vector_scalar_unary_traits<E, 
      general_vector_scalar_unary< typename E::value_type, 
                                   functor::log_minimum<typename E::value_type> > 
    >::result_type
  log_min (const vector_expression<E> &e) {
    typedef typename vector_scalar_unary_traits<E, 
        general_vector_scalar_unary< typename E::value_type, 
                                     functor::log_minimum<typename E::value_type> > 
      >::expression_type expression_type;
    return expression_type (e ());
  }

Finally, we can use the new function:

  vector<double> v( 100 );
  // ... fill v ...                                     
  std::cout << log_min( v ) << std::endl;