Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 09:37:27

0001 // Copyright 2004, 2005 The Trustees of Indiana University.
0002 
0003 // Distributed under the Boost Software License, Version 1.0.
0004 // (See accompanying file LICENSE_1_0.txt or copy at
0005 // http://www.boost.org/LICENSE_1_0.txt)
0006 
0007 //  Authors: Douglas Gregor
0008 //           Andrew Lumsdaine
0009 #ifndef BOOST_GRAPH_FRUCHTERMAN_REINGOLD_FORCE_DIRECTED_LAYOUT_HPP
0010 #define BOOST_GRAPH_FRUCHTERMAN_REINGOLD_FORCE_DIRECTED_LAYOUT_HPP
0011 
0012 #include <boost/config/no_tr1/cmath.hpp>
0013 #include <boost/graph/graph_traits.hpp>
0014 #include <boost/graph/named_function_params.hpp>
0015 #include <boost/graph/iteration_macros.hpp>
0016 #include <boost/graph/topology.hpp> // For topology concepts
0017 #include <boost/graph/detail/mpi_include.hpp>
0018 #include <vector>
0019 #include <list>
0020 #include <algorithm> // for std::min and std::max
0021 #include <numeric> // for std::accumulate
0022 #include <cmath> // for std::sqrt and std::fabs
0023 #include <functional>
0024 
0025 namespace boost
0026 {
0027 
0028 struct square_distance_attractive_force
0029 {
0030     template < typename Graph, typename T >
0031     T operator()(typename graph_traits< Graph >::edge_descriptor, T k, T d,
0032         const Graph&) const
0033     {
0034         return d * d / k;
0035     }
0036 };
0037 
0038 struct square_distance_repulsive_force
0039 {
0040     template < typename Graph, typename T >
0041     T operator()(typename graph_traits< Graph >::vertex_descriptor,
0042         typename graph_traits< Graph >::vertex_descriptor, T k, T d,
0043         const Graph&) const
0044     {
0045         return k * k / d;
0046     }
0047 };
0048 
0049 template < typename T > struct linear_cooling
0050 {
0051     typedef T result_type;
0052 
0053     linear_cooling(std::size_t iterations)
0054     : temp(T(iterations) / T(10)), step(0.1)
0055     {
0056     }
0057 
0058     linear_cooling(std::size_t iterations, T temp)
0059     : temp(temp), step(temp / T(iterations))
0060     {
0061     }
0062 
0063     T operator()()
0064     {
0065         T old_temp = temp;
0066         temp -= step;
0067         if (temp < T(0))
0068             temp = T(0);
0069         return old_temp;
0070     }
0071 
0072 private:
0073     T temp;
0074     T step;
0075 };
0076 
0077 struct all_force_pairs
0078 {
0079     template < typename Graph, typename ApplyForce >
0080     void operator()(const Graph& g, ApplyForce apply_force)
0081     {
0082         typedef typename graph_traits< Graph >::vertex_iterator vertex_iterator;
0083         vertex_iterator v, end;
0084         for (boost::tie(v, end) = vertices(g); v != end; ++v)
0085         {
0086             vertex_iterator u = v;
0087             for (++u; u != end; ++u)
0088             {
0089                 apply_force(*u, *v);
0090                 apply_force(*v, *u);
0091             }
0092         }
0093     }
0094 };
0095 
0096 template < typename Topology, typename PositionMap > struct grid_force_pairs
0097 {
0098     typedef typename property_traits< PositionMap >::value_type Point;
0099     BOOST_STATIC_ASSERT(Point::dimensions == 2);
0100     typedef typename Topology::point_difference_type point_difference_type;
0101 
0102     template < typename Graph >
0103     explicit grid_force_pairs(
0104         const Topology& topology, PositionMap position, const Graph& g)
0105     : topology(topology), position(position)
0106     {
0107         two_k = 2. * this->topology.volume(this->topology.extent())
0108             / std::sqrt((double)num_vertices(g));
0109     }
0110 
0111     template < typename Graph, typename ApplyForce >
0112     void operator()(const Graph& g, ApplyForce apply_force)
0113     {
0114         typedef typename graph_traits< Graph >::vertex_iterator vertex_iterator;
0115         typedef
0116             typename graph_traits< Graph >::vertex_descriptor vertex_descriptor;
0117         typedef std::list< vertex_descriptor > bucket_t;
0118         typedef std::vector< bucket_t > buckets_t;
0119 
0120         std::size_t columns = std::size_t(topology.extent()[0] / two_k + 1.);
0121         std::size_t rows = std::size_t(topology.extent()[1] / two_k + 1.);
0122         buckets_t buckets(rows * columns);
0123         vertex_iterator v, v_end;
0124         for (boost::tie(v, v_end) = vertices(g); v != v_end; ++v)
0125         {
0126             std::size_t column = std::size_t(
0127                 (get(position, *v)[0] + topology.extent()[0] / 2) / two_k);
0128             std::size_t row = std::size_t(
0129                 (get(position, *v)[1] + topology.extent()[1] / 2) / two_k);
0130 
0131             if (column >= columns)
0132                 column = columns - 1;
0133             if (row >= rows)
0134                 row = rows - 1;
0135             buckets[row * columns + column].push_back(*v);
0136         }
0137 
0138         for (std::size_t row = 0; row < rows; ++row)
0139             for (std::size_t column = 0; column < columns; ++column)
0140             {
0141                 bucket_t& bucket = buckets[row * columns + column];
0142                 typedef typename bucket_t::iterator bucket_iterator;
0143                 for (bucket_iterator u = bucket.begin(); u != bucket.end(); ++u)
0144                 {
0145                     // Repulse vertices in this bucket
0146                     bucket_iterator v = u;
0147                     for (++v; v != bucket.end(); ++v)
0148                     {
0149                         apply_force(*u, *v);
0150                         apply_force(*v, *u);
0151                     }
0152 
0153                     std::size_t adj_start_row = row == 0 ? 0 : row - 1;
0154                     std::size_t adj_end_row = row == rows - 1 ? row : row + 1;
0155                     std::size_t adj_start_column = column == 0 ? 0 : column - 1;
0156                     std::size_t adj_end_column
0157                         = column == columns - 1 ? column : column + 1;
0158                     for (std::size_t other_row = adj_start_row;
0159                          other_row <= adj_end_row; ++other_row)
0160                         for (std::size_t other_column = adj_start_column;
0161                              other_column <= adj_end_column; ++other_column)
0162                             if (other_row != row || other_column != column)
0163                             {
0164                                 // Repulse vertices in this bucket
0165                                 bucket_t& other_bucket
0166                                     = buckets[other_row * columns
0167                                         + other_column];
0168                                 for (v = other_bucket.begin();
0169                                      v != other_bucket.end(); ++v)
0170                                 {
0171                                     double dist = topology.distance(
0172                                         get(position, *u), get(position, *v));
0173                                     if (dist < two_k)
0174                                         apply_force(*u, *v);
0175                                 }
0176                             }
0177                 }
0178             }
0179     }
0180 
0181 private:
0182     const Topology& topology;
0183     PositionMap position;
0184     double two_k;
0185 };
0186 
0187 template < typename PositionMap, typename Topology, typename Graph >
0188 inline grid_force_pairs< Topology, PositionMap > make_grid_force_pairs(
0189     const Topology& topology, const PositionMap& position, const Graph& g)
0190 {
0191     return grid_force_pairs< Topology, PositionMap >(topology, position, g);
0192 }
0193 
0194 template < typename Graph, typename PositionMap, typename Topology >
0195 void scale_graph(const Graph& g, PositionMap position, const Topology& topology,
0196     typename Topology::point_type upper_left,
0197     typename Topology::point_type lower_right)
0198 {
0199     if (num_vertices(g) == 0)
0200         return;
0201 
0202     typedef typename Topology::point_type Point;
0203     typedef typename Topology::point_difference_type point_difference_type;
0204 
0205     // Find min/max ranges
0206     Point min_point = get(position, *vertices(g).first), max_point = min_point;
0207     BGL_FORALL_VERTICES_T(v, g, Graph)
0208     {
0209         min_point = topology.pointwise_min(min_point, get(position, v));
0210         max_point = topology.pointwise_max(max_point, get(position, v));
0211     }
0212 
0213     Point old_origin = topology.move_position_toward(min_point, 0.5, max_point);
0214     Point new_origin
0215         = topology.move_position_toward(upper_left, 0.5, lower_right);
0216     point_difference_type old_size = topology.difference(max_point, min_point);
0217     point_difference_type new_size
0218         = topology.difference(lower_right, upper_left);
0219 
0220     // Scale to bounding box provided
0221     BGL_FORALL_VERTICES_T(v, g, Graph)
0222     {
0223         point_difference_type relative_loc
0224             = topology.difference(get(position, v), old_origin);
0225         relative_loc = (relative_loc / old_size) * new_size;
0226         put(position, v, topology.adjust(new_origin, relative_loc));
0227     }
0228 }
0229 
0230 namespace detail
0231 {
0232     template < typename Topology, typename PropMap, typename Vertex >
0233     void maybe_jitter_point(const Topology& topology, const PropMap& pm,
0234         Vertex v, const typename Topology::point_type& p2)
0235     {
0236         double too_close = topology.norm(topology.extent()) / 10000.;
0237         if (topology.distance(get(pm, v), p2) < too_close)
0238         {
0239             put(pm, v,
0240                 topology.move_position_toward(
0241                     get(pm, v), 1. / 200, topology.random_point()));
0242         }
0243     }
0244 
0245     template < typename Topology, typename PositionMap,
0246         typename DisplacementMap, typename RepulsiveForce, typename Graph >
0247     struct fr_apply_force
0248     {
0249         typedef
0250             typename graph_traits< Graph >::vertex_descriptor vertex_descriptor;
0251         typedef typename Topology::point_type Point;
0252         typedef typename Topology::point_difference_type PointDiff;
0253 
0254         fr_apply_force(const Topology& topology, const PositionMap& position,
0255             const DisplacementMap& displacement, RepulsiveForce repulsive_force,
0256             double k, const Graph& g)
0257         : topology(topology)
0258         , position(position)
0259         , displacement(displacement)
0260         , repulsive_force(repulsive_force)
0261         , k(k)
0262         , g(g)
0263         {
0264         }
0265 
0266         void operator()(vertex_descriptor u, vertex_descriptor v)
0267         {
0268             if (u != v)
0269             {
0270                 // When the vertices land on top of each other, move the
0271                 // first vertex away from the boundaries.
0272                 maybe_jitter_point(topology, position, u, get(position, v));
0273 
0274                 double dist
0275                     = topology.distance(get(position, u), get(position, v));
0276                 typename Topology::point_difference_type dispv
0277                     = get(displacement, v);
0278                 if (dist == 0.)
0279                 {
0280                     for (std::size_t i = 0; i < Point::dimensions; ++i)
0281                     {
0282                         dispv[i] += 0.01;
0283                     }
0284                 }
0285                 else
0286                 {
0287                     double fr = repulsive_force(u, v, k, dist, g);
0288                     dispv += (fr / dist)
0289                         * topology.difference(
0290                             get(position, v), get(position, u));
0291                 }
0292                 put(displacement, v, dispv);
0293             }
0294         }
0295 
0296     private:
0297         const Topology& topology;
0298         PositionMap position;
0299         DisplacementMap displacement;
0300         RepulsiveForce repulsive_force;
0301         double k;
0302         const Graph& g;
0303     };
0304 
0305 } // end namespace detail
0306 
0307 template < typename Topology, typename Graph, typename PositionMap,
0308     typename AttractiveForce, typename RepulsiveForce, typename ForcePairs,
0309     typename Cooling, typename DisplacementMap >
0310 void fruchterman_reingold_force_directed_layout(const Graph& g,
0311     PositionMap position, const Topology& topology,
0312     AttractiveForce attractive_force, RepulsiveForce repulsive_force,
0313     ForcePairs force_pairs, Cooling cool, DisplacementMap displacement)
0314 {
0315     typedef typename graph_traits< Graph >::vertex_iterator vertex_iterator;
0316     typedef typename graph_traits< Graph >::vertex_descriptor vertex_descriptor;
0317     typedef typename graph_traits< Graph >::edge_iterator edge_iterator;
0318 
0319     double volume = topology.volume(topology.extent());
0320 
0321     // assume positions are initialized randomly
0322     double k = pow(volume / num_vertices(g),
0323         1. / (double)(Topology::point_difference_type::dimensions));
0324 
0325     detail::fr_apply_force< Topology, PositionMap, DisplacementMap,
0326         RepulsiveForce, Graph >
0327         apply_force(topology, position, displacement, repulsive_force, k, g);
0328 
0329     do
0330     {
0331         // Calculate repulsive forces
0332         vertex_iterator v, v_end;
0333         for (boost::tie(v, v_end) = vertices(g); v != v_end; ++v)
0334             put(displacement, *v, typename Topology::point_difference_type());
0335         force_pairs(g, apply_force);
0336 
0337         // Calculate attractive forces
0338         edge_iterator e, e_end;
0339         for (boost::tie(e, e_end) = edges(g); e != e_end; ++e)
0340         {
0341             vertex_descriptor v = source(*e, g);
0342             vertex_descriptor u = target(*e, g);
0343 
0344             // When the vertices land on top of each other, move the
0345             // first vertex away from the boundaries.
0346             ::boost::detail::maybe_jitter_point(
0347                 topology, position, u, get(position, v));
0348 
0349             typename Topology::point_difference_type delta
0350                 = topology.difference(get(position, v), get(position, u));
0351             double dist = topology.distance(get(position, u), get(position, v));
0352             double fa = attractive_force(*e, k, dist, g);
0353 
0354             put(displacement, v, get(displacement, v) - (fa / dist) * delta);
0355             put(displacement, u, get(displacement, u) + (fa / dist) * delta);
0356         }
0357 
0358         if (double temp = cool())
0359         {
0360             // Update positions
0361             BGL_FORALL_VERTICES_T(v, g, Graph)
0362             {
0363                 BOOST_USING_STD_MIN();
0364                 BOOST_USING_STD_MAX();
0365                 double disp_size = topology.norm(get(displacement, v));
0366                 put(position, v,
0367                     topology.adjust(get(position, v),
0368                         get(displacement, v)
0369                             * (min BOOST_PREVENT_MACRO_SUBSTITUTION(
0370                                    disp_size, temp)
0371                                 / disp_size)));
0372                 put(position, v, topology.bound(get(position, v)));
0373             }
0374         }
0375         else
0376         {
0377             break;
0378         }
0379     } while (true);
0380 }
0381 
0382 namespace detail
0383 {
0384     template < typename DisplacementMap > struct fr_force_directed_layout
0385     {
0386         template < typename Topology, typename Graph, typename PositionMap,
0387             typename AttractiveForce, typename RepulsiveForce,
0388             typename ForcePairs, typename Cooling, typename Param, typename Tag,
0389             typename Rest >
0390         static void run(const Graph& g, PositionMap position,
0391             const Topology& topology, AttractiveForce attractive_force,
0392             RepulsiveForce repulsive_force, ForcePairs force_pairs,
0393             Cooling cool, DisplacementMap displacement,
0394             const bgl_named_params< Param, Tag, Rest >&)
0395         {
0396             fruchterman_reingold_force_directed_layout(g, position, topology,
0397                 attractive_force, repulsive_force, force_pairs, cool,
0398                 displacement);
0399         }
0400     };
0401 
0402     template <> struct fr_force_directed_layout< param_not_found >
0403     {
0404         template < typename Topology, typename Graph, typename PositionMap,
0405             typename AttractiveForce, typename RepulsiveForce,
0406             typename ForcePairs, typename Cooling, typename Param, typename Tag,
0407             typename Rest >
0408         static void run(const Graph& g, PositionMap position,
0409             const Topology& topology, AttractiveForce attractive_force,
0410             RepulsiveForce repulsive_force, ForcePairs force_pairs,
0411             Cooling cool, param_not_found,
0412             const bgl_named_params< Param, Tag, Rest >& params)
0413         {
0414             typedef typename Topology::point_difference_type PointDiff;
0415             std::vector< PointDiff > displacements(num_vertices(g));
0416             fruchterman_reingold_force_directed_layout(g, position, topology,
0417                 attractive_force, repulsive_force, force_pairs, cool,
0418                 make_iterator_property_map(displacements.begin(),
0419                     choose_const_pmap(
0420                         get_param(params, vertex_index), g, vertex_index),
0421                     PointDiff()));
0422         }
0423     };
0424 
0425 } // end namespace detail
0426 
0427 template < typename Topology, typename Graph, typename PositionMap,
0428     typename Param, typename Tag, typename Rest >
0429 void fruchterman_reingold_force_directed_layout(const Graph& g,
0430     PositionMap position, const Topology& topology,
0431     const bgl_named_params< Param, Tag, Rest >& params)
0432 {
0433     typedef typename get_param_type< vertex_displacement_t,
0434         bgl_named_params< Param, Tag, Rest > >::type D;
0435 
0436     detail::fr_force_directed_layout< D >::run(g, position, topology,
0437         choose_param(get_param(params, attractive_force_t()),
0438             square_distance_attractive_force()),
0439         choose_param(get_param(params, repulsive_force_t()),
0440             square_distance_repulsive_force()),
0441         choose_param(get_param(params, force_pairs_t()),
0442             make_grid_force_pairs(topology, position, g)),
0443         choose_param(
0444             get_param(params, cooling_t()), linear_cooling< double >(100)),
0445         get_param(params, vertex_displacement_t()), params);
0446 }
0447 
0448 template < typename Topology, typename Graph, typename PositionMap >
0449 void fruchterman_reingold_force_directed_layout(
0450     const Graph& g, PositionMap position, const Topology& topology)
0451 {
0452     fruchterman_reingold_force_directed_layout(g, position, topology,
0453         attractive_force(square_distance_attractive_force()));
0454 }
0455 
0456 } // end namespace boost
0457 
0458 #include BOOST_GRAPH_MPI_INCLUDE(<boost/graph/distributed/fruchterman_reingold.hpp>)
0459 
0460 #endif // BOOST_GRAPH_FRUCHTERMAN_REINGOLD_FORCE_DIRECTED_LAYOUT_HPP