Chris@16: // Copyright (C) 2005-2006 The Trustees of Indiana University.
Chris@16: 
Chris@16: // Use, modification and distribution is subject to the Boost Software
Chris@16: // License, Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
Chris@16: // http://www.boost.org/LICENSE_1_0.txt)
Chris@16: 
Chris@16: //  Authors: Douglas Gregor
Chris@16: //           Andrew Lumsdaine
Chris@16: #ifndef BOOST_GRAPH_DISTRIBUTED_FRUCHTERMAN_REINGOLD_HPP
Chris@16: #define BOOST_GRAPH_DISTRIBUTED_FRUCHTERMAN_REINGOLD_HPP
Chris@16: 
Chris@16: #ifndef BOOST_GRAPH_USE_MPI
Chris@16: #error "Parallel BGL files should not be included unless <boost/graph/use_mpi.hpp> has been included"
Chris@16: #endif
Chris@16: 
Chris@16: #include <boost/graph/fruchterman_reingold.hpp>
Chris@16: 
Chris@16: namespace boost { namespace graph { namespace distributed {
Chris@16: 
Chris@16: class simple_tiling
Chris@16: {
Chris@16:  public:
Chris@16:   simple_tiling(int columns, int rows, bool flip = true) 
Chris@16:     : columns(columns), rows(rows), flip(flip)
Chris@16:   {
Chris@16:   }
Chris@16: 
Chris@16:   // Convert from a position (x, y) in the tiled display into a
Chris@16:   // processor ID number
Chris@16:   int operator()(int x, int y) const
Chris@16:   {
Chris@16:     return flip? (rows - y - 1) * columns + x : y * columns + x;
Chris@16:   }
Chris@16: 
Chris@16:   // Convert from a process ID to a position (x, y) in the tiled
Chris@16:   // display
Chris@16:   std::pair<int, int> operator()(int id)
Chris@16:   {
Chris@16:     int my_col = id % columns;
Chris@16:     int my_row = flip? rows - (id / columns) - 1 : id / columns;
Chris@16:     return std::make_pair(my_col, my_row);
Chris@16:   }
Chris@16: 
Chris@16:   int columns, rows;
Chris@16: 
Chris@16:  private:
Chris@16:   bool flip;
Chris@16: };
Chris@16: 
Chris@16: // Force pairs function object that does nothing
Chris@16: struct no_force_pairs
Chris@16: {
Chris@16:   template<typename Graph, typename ApplyForce>
Chris@16:   void operator()(const Graph&, const ApplyForce&)
Chris@16:   {
Chris@16:   }
Chris@16: };
Chris@16: 
Chris@16: // Computes force pairs in the distributed case.
Chris@16: template<typename PositionMap, typename DisplacementMap, typename LocalForces,
Chris@16:          typename NonLocalForces = no_force_pairs>
Chris@16: class distributed_force_pairs_proxy
Chris@16: {
Chris@16:  public:
Chris@16:   distributed_force_pairs_proxy(const PositionMap& position, 
Chris@16:                                 const DisplacementMap& displacement,
Chris@16:                                 const LocalForces& local_forces,
Chris@16:                                 const NonLocalForces& nonlocal_forces = NonLocalForces())
Chris@16:     : position(position), displacement(displacement), 
Chris@16:       local_forces(local_forces), nonlocal_forces(nonlocal_forces)
Chris@16:   {
Chris@16:   }
Chris@16: 
Chris@16:   template<typename Graph, typename ApplyForce>
Chris@16:   void operator()(const Graph& g, ApplyForce apply_force)
Chris@16:   {
Chris@16:     // Flush remote displacements
Chris@16:     displacement.flush();
Chris@16: 
Chris@16:     // Receive updated positions for all of our neighbors
Chris@16:     synchronize(position);
Chris@16: 
Chris@16:     // Reset remote displacements 
Chris@16:     displacement.reset();
Chris@16: 
Chris@16:     // Compute local repulsive forces
Chris@16:     local_forces(g, apply_force);
Chris@16: 
Chris@16:     // Compute neighbor repulsive forces
Chris@16:     nonlocal_forces(g, apply_force);
Chris@16:   }
Chris@16: 
Chris@16:  protected:
Chris@16:   PositionMap position;
Chris@16:   DisplacementMap displacement;
Chris@16:   LocalForces local_forces;
Chris@16:   NonLocalForces nonlocal_forces;
Chris@16: };
Chris@16: 
Chris@16: template<typename PositionMap, typename DisplacementMap, typename LocalForces>
Chris@16: inline 
Chris@16: distributed_force_pairs_proxy<PositionMap, DisplacementMap, LocalForces>
Chris@16: make_distributed_force_pairs(const PositionMap& position, 
Chris@16:                              const DisplacementMap& displacement,
Chris@16:                              const LocalForces& local_forces)
Chris@16: {
Chris@16:   typedef 
Chris@16:     distributed_force_pairs_proxy<PositionMap, DisplacementMap, LocalForces>
Chris@16:     result_type;
Chris@16:   return result_type(position, displacement, local_forces);
Chris@16: }
Chris@16: 
Chris@16: template<typename PositionMap, typename DisplacementMap, typename LocalForces,
Chris@16:          typename NonLocalForces>
Chris@16: inline 
Chris@16: distributed_force_pairs_proxy<PositionMap, DisplacementMap, LocalForces,
Chris@16:                               NonLocalForces>
Chris@16: make_distributed_force_pairs(const PositionMap& position, 
Chris@16:                              const DisplacementMap& displacement,
Chris@16:                              const LocalForces& local_forces,
Chris@16:                              const NonLocalForces& nonlocal_forces)
Chris@16: {
Chris@16:   typedef 
Chris@16:     distributed_force_pairs_proxy<PositionMap, DisplacementMap, LocalForces,
Chris@16:                                   NonLocalForces>
Chris@16:       result_type;
Chris@16:   return result_type(position, displacement, local_forces, nonlocal_forces);
Chris@16: }
Chris@16: 
Chris@16: // Compute nonlocal force pairs based on the shared borders with
Chris@16: // adjacent tiles.
Chris@16: template<typename PositionMap>
Chris@16: class neighboring_tiles_force_pairs
Chris@16: {
Chris@16:  public:
Chris@16:   typedef typename property_traits<PositionMap>::value_type Point;
Chris@16:   typedef typename point_traits<Point>::component_type Dim;
Chris@16: 
Chris@16:   enum bucket_position { left, top, right, bottom, end_position };
Chris@16:   
Chris@16:   neighboring_tiles_force_pairs(PositionMap position, Point origin,
Chris@16:                                 Point extent, simple_tiling tiling)
Chris@16:     : position(position), origin(origin), extent(extent), tiling(tiling)
Chris@16:   {
Chris@16:   }
Chris@16: 
Chris@16:   template<typename Graph, typename ApplyForce>
Chris@16:   void operator()(const Graph& g, ApplyForce apply_force)
Chris@16:   {
Chris@16:     // TBD: Do this some smarter way
Chris@16:     if (tiling.columns == 1 && tiling.rows == 1)
Chris@16:       return;
Chris@16: 
Chris@16:     typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;
Chris@16: #ifndef BOOST_NO_STDC_NAMESPACE
Chris@16:     using std::sqrt;
Chris@16: #endif // BOOST_NO_STDC_NAMESPACE
Chris@16: 
Chris@16:     // TBD: num_vertices(g) should be the global number of vertices?
Chris@16:     Dim two_k = Dim(2) * sqrt(extent[0] * extent[1] / num_vertices(g));
Chris@16: 
Chris@16:     std::vector<vertex_descriptor> my_vertices[4];
Chris@16:     std::vector<vertex_descriptor> neighbor_vertices[4];
Chris@16:    
Chris@16:     // Compute cutoff positions
Chris@16:     Dim cutoffs[4];
Chris@16:     cutoffs[left] = origin[0] + two_k;
Chris@16:     cutoffs[top] = origin[1] + two_k;
Chris@16:     cutoffs[right] = origin[0] + extent[0] - two_k;
Chris@16:     cutoffs[bottom] = origin[1] + extent[1] - two_k;
Chris@16: 
Chris@16:     // Compute neighbors
Chris@16:     typename PositionMap::process_group_type pg = position.process_group();
Chris@16:     std::pair<int, int> my_tile = tiling(process_id(pg));
Chris@16:     int neighbors[4] = { -1, -1, -1, -1 } ;
Chris@16:     if (my_tile.first > 0)
Chris@16:       neighbors[left] = tiling(my_tile.first - 1, my_tile.second);
Chris@16:     if (my_tile.second > 0)
Chris@16:       neighbors[top] = tiling(my_tile.first, my_tile.second - 1);
Chris@16:     if (my_tile.first < tiling.columns - 1)
Chris@16:       neighbors[right] = tiling(my_tile.first + 1, my_tile.second);
Chris@16:     if (my_tile.second < tiling.rows - 1)
Chris@16:       neighbors[bottom] = tiling(my_tile.first, my_tile.second + 1);
Chris@16: 
Chris@16:     // Sort vertices along the edges into buckets
Chris@16:     BGL_FORALL_VERTICES_T(v, g, Graph) {
Chris@16:       if (position[v][0] <= cutoffs[left]) my_vertices[left].push_back(v); 
Chris@16:       if (position[v][1] <= cutoffs[top]) my_vertices[top].push_back(v); 
Chris@16:       if (position[v][0] >= cutoffs[right]) my_vertices[right].push_back(v); 
Chris@16:       if (position[v][1] >= cutoffs[bottom]) my_vertices[bottom].push_back(v); 
Chris@16:     }
Chris@16: 
Chris@16:     // Send vertices to neighbors, and gather our neighbors' vertices
Chris@16:     bucket_position pos;
Chris@16:     for (pos = left; pos < end_position; pos = bucket_position(pos + 1)) {
Chris@16:       if (neighbors[pos] != -1) {
Chris@16:         send(pg, neighbors[pos], 0, my_vertices[pos].size());
Chris@16:         if (!my_vertices[pos].empty())
Chris@16:           send(pg, neighbors[pos], 1, 
Chris@16:                &my_vertices[pos].front(), my_vertices[pos].size());
Chris@16:       }
Chris@16:     }
Chris@16: 
Chris@16:     // Pass messages around
Chris@16:     synchronize(pg);
Chris@16:     
Chris@16:     // Receive neighboring vertices
Chris@16:     for (pos = left; pos < end_position; pos = bucket_position(pos + 1)) {
Chris@16:       if (neighbors[pos] != -1) {
Chris@16:         std::size_t incoming_vertices;
Chris@16:         receive(pg, neighbors[pos], 0, incoming_vertices);
Chris@16: 
Chris@16:         if (incoming_vertices) {
Chris@16:           neighbor_vertices[pos].resize(incoming_vertices);
Chris@16:           receive(pg, neighbors[pos], 1, &neighbor_vertices[pos].front(),
Chris@16:                   incoming_vertices);
Chris@16:         }
Chris@16:       }
Chris@16:     }
Chris@16: 
Chris@16:     // For each neighboring vertex, we need to get its current position
Chris@16:     for (pos = left; pos < end_position; pos = bucket_position(pos + 1))
Chris@16:       for (typename std::vector<vertex_descriptor>::iterator i = 
Chris@16:              neighbor_vertices[pos].begin(); 
Chris@16:            i != neighbor_vertices[pos].end();
Chris@16:            ++i)
Chris@16:         request(position, *i);
Chris@16:     synchronize(position);
Chris@16: 
Chris@16:     // Apply forces in adjacent bins. This is O(n^2) in the worst
Chris@16:     // case. Oh well.
Chris@16:     for (pos = left; pos < end_position; pos = bucket_position(pos + 1)) {
Chris@16:       for (typename std::vector<vertex_descriptor>::iterator i = 
Chris@16:              my_vertices[pos].begin(); 
Chris@16:            i != my_vertices[pos].end();
Chris@16:            ++i)
Chris@16:         for (typename std::vector<vertex_descriptor>::iterator j = 
Chris@16:                neighbor_vertices[pos].begin(); 
Chris@16:              j != neighbor_vertices[pos].end();
Chris@16:              ++j)
Chris@16:           apply_force(*i, *j);
Chris@16:     }
Chris@16:   }
Chris@16: 
Chris@16:  protected:
Chris@16:   PositionMap position;
Chris@16:   Point origin;
Chris@16:   Point extent;
Chris@16:   simple_tiling tiling;
Chris@16: };
Chris@16: 
Chris@16: template<typename PositionMap>
Chris@16: inline neighboring_tiles_force_pairs<PositionMap>
Chris@16: make_neighboring_tiles_force_pairs
Chris@16:  (PositionMap position, 
Chris@16:   typename property_traits<PositionMap>::value_type origin,
Chris@16:   typename property_traits<PositionMap>::value_type extent,
Chris@16:   simple_tiling tiling)
Chris@16: {
Chris@16:   return neighboring_tiles_force_pairs<PositionMap>(position, origin, extent,
Chris@16:                                                     tiling);
Chris@16: }
Chris@16: 
Chris@16: template<typename DisplacementMap, typename Cooling>
Chris@16: class distributed_cooling_proxy
Chris@16: {
Chris@16:  public:
Chris@16:   typedef typename Cooling::result_type result_type;
Chris@16: 
Chris@16:   distributed_cooling_proxy(const DisplacementMap& displacement,
Chris@16:                             const Cooling& cooling)
Chris@16:     : displacement(displacement), cooling(cooling)
Chris@16:   {
Chris@16:   }
Chris@16: 
Chris@16:   result_type operator()()
Chris@16:   {
Chris@16:     // Accumulate displacements computed on each processor
Chris@16:     synchronize(displacement.data->process_group);
Chris@16: 
Chris@16:     // Allow the underlying cooling to occur
Chris@16:     return cooling();
Chris@16:   }
Chris@16: 
Chris@16:  protected:
Chris@16:   DisplacementMap displacement;
Chris@16:   Cooling cooling;
Chris@16: };
Chris@16: 
Chris@16: template<typename DisplacementMap, typename Cooling>
Chris@16: inline distributed_cooling_proxy<DisplacementMap, Cooling>
Chris@16: make_distributed_cooling(const DisplacementMap& displacement,
Chris@16:                          const Cooling& cooling)
Chris@16: {
Chris@16:   typedef distributed_cooling_proxy<DisplacementMap, Cooling> result_type;
Chris@16:   return result_type(displacement, cooling);
Chris@16: }
Chris@16: 
Chris@16: template<typename Point>
Chris@16: struct point_accumulating_reducer {
Chris@16:   BOOST_STATIC_CONSTANT(bool, non_default_resolver = true);
Chris@16: 
Chris@16:   template<typename K>
Chris@16:   Point operator()(const K&) const { return Point(); }
Chris@16: 
Chris@16:   template<typename K>
Chris@16:   Point operator()(const K&, const Point& p1, const Point& p2) const 
Chris@16:   { return Point(p1[0] + p2[0], p1[1] + p2[1]); }
Chris@16: };
Chris@16: 
Chris@16: template<typename Graph, typename PositionMap, 
Chris@16:          typename AttractiveForce, typename RepulsiveForce,
Chris@16:          typename ForcePairs, typename Cooling, typename DisplacementMap>
Chris@16: void
Chris@16: fruchterman_reingold_force_directed_layout
Chris@16:  (const Graph&    g,
Chris@16:   PositionMap     position,
Chris@16:   typename property_traits<PositionMap>::value_type const& origin,
Chris@16:   typename property_traits<PositionMap>::value_type const& extent,
Chris@16:   AttractiveForce attractive_force,
Chris@16:   RepulsiveForce  repulsive_force,
Chris@16:   ForcePairs      force_pairs,
Chris@16:   Cooling         cool,
Chris@16:   DisplacementMap displacement)
Chris@16: {
Chris@16:   typedef typename property_traits<PositionMap>::value_type Point;
Chris@16: 
Chris@16:   // Reduction in the displacement map involves summing the forces
Chris@16:   displacement.set_reduce(point_accumulating_reducer<Point>());
Chris@16: 
Chris@16:   // We need to track the positions of all of our neighbors
Chris@16:   BGL_FORALL_VERTICES_T(u, g, Graph)
Chris@16:     BGL_FORALL_ADJ_T(u, v, g, Graph)
Chris@16:       request(position, v);
Chris@16: 
Chris@16:   // Invoke the "sequential" Fruchterman-Reingold implementation
Chris@16:   boost::fruchterman_reingold_force_directed_layout
Chris@16:     (g, position, origin, extent,
Chris@16:      attractive_force, repulsive_force,
Chris@16:      make_distributed_force_pairs(position, displacement, force_pairs),
Chris@16:      make_distributed_cooling(displacement, cool),
Chris@16:      displacement);
Chris@16: }
Chris@16: 
Chris@16: template<typename Graph, typename PositionMap, 
Chris@16:          typename AttractiveForce, typename RepulsiveForce,
Chris@16:          typename ForcePairs, typename Cooling, typename DisplacementMap>
Chris@16: void
Chris@16: fruchterman_reingold_force_directed_layout
Chris@16:  (const Graph&    g,
Chris@16:   PositionMap     position,
Chris@16:   typename property_traits<PositionMap>::value_type const& origin,
Chris@16:   typename property_traits<PositionMap>::value_type const& extent,
Chris@16:   AttractiveForce attractive_force,
Chris@16:   RepulsiveForce  repulsive_force,
Chris@16:   ForcePairs      force_pairs,
Chris@16:   Cooling         cool,
Chris@16:   DisplacementMap displacement,
Chris@16:   simple_tiling   tiling)
Chris@16: {
Chris@16:   typedef typename property_traits<PositionMap>::value_type Point;
Chris@16: 
Chris@16:   // Reduction in the displacement map involves summing the forces
Chris@16:   displacement.set_reduce(point_accumulating_reducer<Point>());
Chris@16: 
Chris@16:   // We need to track the positions of all of our neighbors
Chris@16:   BGL_FORALL_VERTICES_T(u, g, Graph)
Chris@16:     BGL_FORALL_ADJ_T(u, v, g, Graph)
Chris@16:       request(position, v);
Chris@16: 
Chris@16:   // Invoke the "sequential" Fruchterman-Reingold implementation
Chris@16:   boost::fruchterman_reingold_force_directed_layout
Chris@16:     (g, position, origin, extent,
Chris@16:      attractive_force, repulsive_force,
Chris@16:      make_distributed_force_pairs
Chris@16:       (position, displacement, force_pairs,
Chris@16:        make_neighboring_tiles_force_pairs(position, origin, extent, tiling)),
Chris@16:      make_distributed_cooling(displacement, cool),
Chris@16:      displacement);
Chris@16: }
Chris@16: 
Chris@16: } } } // end namespace boost::graph::distributed
Chris@16: 
Chris@16: #endif // BOOST_GRAPH_DISTRIBUTED_FRUCHTERMAN_REINGOLD_HPP