How to adapt a boost::property_map in order to use it like a ublas::vector?

Question

I am looking for a smart way to use the Boost Graph Library in conjunction with Boost uBLAS. More precisely, I need to update a given vertex property for each vertex by using the result of the scalar product between the graphs adjacency matrix and a vector containing some other vertex property for each vertex. Let me give you a (unfortunately lenghty) minimal example to illustrate the problem:

#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/iteration_macros.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>

using namespace boost;
namespace ublas = boost::numeric::ublas;


struct Vertex {             //Using bundled vertex properties

        double old_potential; 
        double new_potential; 
};


typedef adjacency_list< listS, vecS, directedS, Vertex > Graph;

int main(){ 

  //[Prepare a graph with two vertex property maps and initialize
  Graph graph; 
  add_edge (0, 1, graph); 
  add_edge (0, 3, graph); 
  add_edge (1, 2, graph); 

  auto  v_old_potential = boost::get( &Vertex::old_potential, graph );
  auto  v_new_potential = boost::get( &Vertex::new_potential, graph );

  unsigned int source_strength = 0;

  BGL_FORALL_VERTICES( v, graph, Graph )    {

    v_old_potential[v] = source_strength++;
    v_new_potential[v] = 0;
  }
  //]  


  //[ Extracting the adjacency matrix by iteration through all edges --> uBLAS matrix
  ublas::zero_matrix<int> zero_matrix( num_vertices(graph) , num_vertices(graph) ); 
  ublas::matrix<int> adjacency_matrix( zero_matrix );       //initialize properly

  BGL_FORALL_EDGES( e, graph, Graph )   {

     adjacency_matrix( source(e,graph), target(e,graph) ) = 1;
     adjacency_matrix( target(e,graph), source(e,graph) ) = 1;
  }
  //]


  //[ Extracting the old potentials by iteration through all vertices --> uBLAS vector
  ublas::zero_vector<double> zero_vector( num_vertices(graph) );
  ublas::vector<double> old_potential_vector( zero_vector );    //initialize properly
  ublas::vector<double> new_potential_vector( zero_vector );    //initialize properly

  BGL_FORALL_VERTICES(v, graph, Graph)    {

        old_potential_vector( vertex(v,graph) ) = v_old_potential[v];
  }
  //]


  //[ Compute new potentials = A . old potentials ! 
  new_potential_vector = ublas::prod ( adjacency_matrix, old_potential_vector );        //     new = A.old
  //]


  //[ Updating the property map for the new potentials with the newly computed values from above
  BGL_FORALL_VERTICES(v, graph, Graph)    {

        v_new_potential[v] = old_potential_vector( vertex(v,graph) );
  }
  //]

  std::cout << adjacency_matrix << std::endl;       //output = [4,4]((0,1,0,1),(1,0,1,0),(0,1,0,0),(1,0,0,0))
  std::cout << old_potential_vector << std::endl;   //output = [4](0,1,2,3)
  std::cout << new_potential_vector << std::endl;   //output = [4](4,2,1,0)

}

Now, while my suggestion is a possible solution, i am not quite satisfied with it. The problem is, (a) I copy the whole old_potential property map to the associated ublas::vector in order to compute the scalar product. And (b) I need to traverse the new_potential property map as well to get the newly computed values back into the graph. I suspect that these operations will be repeated a lot in my application that is why I want to get this part done as clean as possible from the start.

Ideally, I'd like to be done with all this copying back and forth and instead use some sort of adapter to make a boost::property_map work as a ublas::vector in the call of prod(). It would be awesome to use something like this:

adapter(new_potential) = ublas::prod( adjacency_matrix, adapter(old_potential) );

If anyone has an idea how to achieve such functionality or how to improve on my solution I would be very grateful.

Answer 1

#include <iostream>
#include <memory>
#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/iteration_macros.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>

using namespace boost;
namespace ublas = boost::numeric::ublas;

enum vertex_old_potential_t { vertex_old_potential };
enum vertex_new_potential_t { vertex_new_potential };

namespace boost
{
  BOOST_INSTALL_PROPERTY(vertex, new_potential);
  BOOST_INSTALL_PROPERTY(vertex, old_potential);
}

typedef property<vertex_new_potential_t, double, property<vertex_old_potential_t,double> > Vertex;


typedef adjacency_list< listS, vecS, directedS, Vertex > Graph;

struct ublas_vector_map;

namespace boost {
  template<>
  struct property_map< Graph, vertex_new_potential_t > {
    typedef ublas_vector_map type;
    typedef ublas_vector_map const_type;
  };

  template<>
  struct property_map< Graph, vertex_old_potential_t > {
    typedef ublas_vector_map type;
    typedef ublas_vector_map const_type;
  };
}

struct ublas_vector_map : put_get_helper<double&,ublas_vector_map> {
  typedef double value_type;
  typedef value_type& reference;
  typedef typename graph_traits<Graph>::vertex_descriptor key_type;
  typedef boost::lvalue_property_map_tag category;

  ublas_vector_map(Graph* g, vertex_old_potential_t&):vec(new ublas::vector<double>(num_vertices(*g),0.0)){}
  ublas_vector_map(Graph* g, vertex_new_potential_t&):vec(new ublas::vector<double>(num_vertices(*g),0.0)){}


  reference operator[](key_type v) const {
    return (*vec)(v);
  }

  ublas::vector<double>& vector() { return *vec; }

  private:
   std::unique_ptr<ublas::vector<double> > vec;
};

int main(){ 

  //[Prepare a graph with two vertex property maps and initialize
  Graph graph; 
  add_edge (0, 1, graph); 
  add_edge (0, 3, graph); 
  add_edge (1, 2, graph); 

  auto  v_old_potential = boost::get( vertex_old_potential, graph );
  auto  v_new_potential = boost::get( vertex_new_potential, graph );

  unsigned int source_strength = 0;

  BGL_FORALL_VERTICES( v, graph, Graph )    {

    v_old_potential[v] = source_strength++;
  }
  //]  


  //[ Extracting the adjacency matrix by iteration through all edges --> uBLAS matrix
  ublas::zero_matrix<int> zero_matrix( num_vertices(graph) , num_vertices(graph) ); 
  ublas::matrix<int> adjacency_matrix( zero_matrix );       //initialize properly

  BGL_FORALL_EDGES( e, graph, Graph )   {

     adjacency_matrix( source(e,graph), target(e,graph) ) = 1;
     adjacency_matrix( target(e,graph), source(e,graph) ) = 1;
  }
  //]



  //[ Compute new potentials = A . old potentials ! 
  v_new_potential.vector() = ublas::prod ( adjacency_matrix, v_old_potential.vector() );        //     new = A.old
  //]




  std::cout << adjacency_matrix << std::endl;       //output = [4,4]((0,1,0,1),(1,0,1,0),(0,1,0,0),(1,0,0,0))
  std::cout << v_old_potential.vector() << std::endl;   //output = [4](0,1,2,3)
  std::cout << v_new_potential.vector() << std::endl;   //output = [4](4,2,1,0)

  //You must access the properties via v_new_potential and v_old_potential, if you use get... again it resets
  std::cout << v_new_potential[0] << std::endl;
  std::cout << get(vertex_new_potential, graph)[0] << std::endl;

}