// Copyright 2004, 2005 The Trustees of Indiana University.

// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)

//  Authors: Douglas Gregor
//           Andrew Lumsdaine
#ifndef BOOST_GRAPH_FRUCHTERMAN_REINGOLD_FORCE_DIRECTED_LAYOUT_HPP
#define BOOST_GRAPH_FRUCHTERMAN_REINGOLD_FORCE_DIRECTED_LAYOUT_HPP

#include <boost/config/no_tr1/cmath.hpp>
#include <boost/graph/graph_traits.hpp>
#include <boost/graph/named_function_params.hpp>
#include <boost/graph/iteration_macros.hpp>
#include <boost/graph/topology.hpp> // For topology concepts
#include <boost/graph/detail/mpi_include.hpp>
#include <vector>
#include <list>
#include <algorithm> // for std::min and std::max
#include <numeric> // for std::accumulate
#include <cmath> // for std::sqrt and std::fabs
#include <functional>

namespace boost
{

struct square_distance_attractive_force
{
    template < typename Graph, typename T >
    T operator()(typename graph_traits< Graph >::edge_descriptor, T k, T d,
        const Graph&) const
    {
        return d * d / k;
    }
};

struct square_distance_repulsive_force
{
    template < typename Graph, typename T >
    T operator()(typename graph_traits< Graph >::vertex_descriptor,
        typename graph_traits< Graph >::vertex_descriptor, T k, T d,
        const Graph&) const
    {
        return k * k / d;
    }
};

template < typename T > struct linear_cooling
{
    typedef T result_type;

    linear_cooling(std::size_t iterations)
    : temp(T(iterations) / T(10)), step(0.1)
    {
    }

    linear_cooling(std::size_t iterations, T temp)
    : temp(temp), step(temp / T(iterations))
    {
    }

    T operator()()
    {
        T old_temp = temp;
        temp -= step;
        if (temp < T(0))
            temp = T(0);
        return old_temp;
    }

private:
    T temp;
    T step;
};

struct all_force_pairs
{
    template < typename Graph, typename ApplyForce >
    void operator()(const Graph& g, ApplyForce apply_force)
    {
        typedef typename graph_traits< Graph >::vertex_iterator vertex_iterator;
        vertex_iterator v, end;
        for (boost::tie(v, end) = vertices(g); v != end; ++v)
        {
            vertex_iterator u = v;
            for (++u; u != end; ++u)
            {
                apply_force(*u, *v);
                apply_force(*v, *u);
            }
        }
    }
};

template < typename Topology, typename PositionMap > struct grid_force_pairs
{
    typedef typename property_traits< PositionMap >::value_type Point;
    BOOST_STATIC_ASSERT(Point::dimensions == 2);
    typedef typename Topology::point_difference_type point_difference_type;

    template < typename Graph >
    explicit grid_force_pairs(
        const Topology& topology, PositionMap position, const Graph& g)
    : topology(topology), position(position)
    {
        two_k = 2. * this->topology.volume(this->topology.extent())
            / std::sqrt((double)num_vertices(g));
    }

    template < typename Graph, typename ApplyForce >
    void operator()(const Graph& g, ApplyForce apply_force)
    {
        typedef typename graph_traits< Graph >::vertex_iterator vertex_iterator;
        typedef
            typename graph_traits< Graph >::vertex_descriptor vertex_descriptor;
        typedef std::list< vertex_descriptor > bucket_t;
        typedef std::vector< bucket_t > buckets_t;

        std::size_t columns = std::size_t(topology.extent()[0] / two_k + 1.);
        std::size_t rows = std::size_t(topology.extent()[1] / two_k + 1.);
        buckets_t buckets(rows * columns);
        vertex_iterator v, v_end;
        for (boost::tie(v, v_end) = vertices(g); v != v_end; ++v)
        {
            std::size_t column = std::size_t(
                (get(position, *v)[0] + topology.extent()[0] / 2) / two_k);
            std::size_t row = std::size_t(
                (get(position, *v)[1] + topology.extent()[1] / 2) / two_k);

            if (column >= columns)
                column = columns - 1;
            if (row >= rows)
                row = rows - 1;
            buckets[row * columns + column].push_back(*v);
        }

        for (std::size_t row = 0; row < rows; ++row)
            for (std::size_t column = 0; column < columns; ++column)
            {
                bucket_t& bucket = buckets[row * columns + column];
                typedef typename bucket_t::iterator bucket_iterator;
                for (bucket_iterator u = bucket.begin(); u != bucket.end(); ++u)
                {
                    // Repulse vertices in this bucket
                    bucket_iterator v = u;
                    for (++v; v != bucket.end(); ++v)
                    {
                        apply_force(*u, *v);
                        apply_force(*v, *u);
                    }

                    std::size_t adj_start_row = row == 0 ? 0 : row - 1;
                    std::size_t adj_end_row = row == rows - 1 ? row : row + 1;
                    std::size_t adj_start_column = column == 0 ? 0 : column - 1;
                    std::size_t adj_end_column
                        = column == columns - 1 ? column : column + 1;
                    for (std::size_t other_row = adj_start_row;
                         other_row <= adj_end_row; ++other_row)
                        for (std::size_t other_column = adj_start_column;
                             other_column <= adj_end_column; ++other_column)
                            if (other_row != row || other_column != column)
                            {
                                // Repulse vertices in this bucket
                                bucket_t& other_bucket
                                    = buckets[other_row * columns
                                        + other_column];
                                for (v = other_bucket.begin();
                                     v != other_bucket.end(); ++v)
                                {
                                    double dist = topology.distance(
                                        get(position, *u), get(position, *v));
                                    if (dist < two_k)
                                        apply_force(*u, *v);
                                }
                            }
                }
            }
    }

private:
    const Topology& topology;
    PositionMap position;
    double two_k;
};

template < typename PositionMap, typename Topology, typename Graph >
inline grid_force_pairs< Topology, PositionMap > make_grid_force_pairs(
    const Topology& topology, const PositionMap& position, const Graph& g)
{
    return grid_force_pairs< Topology, PositionMap >(topology, position, g);
}

template < typename Graph, typename PositionMap, typename Topology >
void scale_graph(const Graph& g, PositionMap position, const Topology& topology,
    typename Topology::point_type upper_left,
    typename Topology::point_type lower_right)
{
    if (num_vertices(g) == 0)
        return;

    typedef typename Topology::point_type Point;
    typedef typename Topology::point_difference_type point_difference_type;

    // Find min/max ranges
    Point min_point = get(position, *vertices(g).first), max_point = min_point;
    BGL_FORALL_VERTICES_T(v, g, Graph)
    {
        min_point = topology.pointwise_min(min_point, get(position, v));
        max_point = topology.pointwise_max(max_point, get(position, v));
    }

    Point old_origin = topology.move_position_toward(min_point, 0.5, max_point);
    Point new_origin
        = topology.move_position_toward(upper_left, 0.5, lower_right);
    point_difference_type old_size = topology.difference(max_point, min_point);
    point_difference_type new_size
        = topology.difference(lower_right, upper_left);

    // Scale to bounding box provided
    BGL_FORALL_VERTICES_T(v, g, Graph)
    {
        point_difference_type relative_loc
            = topology.difference(get(position, v), old_origin);
        relative_loc = (relative_loc / old_size) * new_size;
        put(position, v, topology.adjust(new_origin, relative_loc));
    }
}

namespace detail
{
    template < typename Topology, typename PropMap, typename Vertex >
    void maybe_jitter_point(const Topology& topology, const PropMap& pm,
        Vertex v, const typename Topology::point_type& p2)
    {
        double too_close = topology.norm(topology.extent()) / 10000.;
        if (topology.distance(get(pm, v), p2) < too_close)
        {
            put(pm, v,
                topology.move_position_toward(
                    get(pm, v), 1. / 200, topology.random_point()));
        }
    }

    template < typename Topology, typename PositionMap,
        typename DisplacementMap, typename RepulsiveForce, typename Graph >
    struct fr_apply_force
    {
        typedef
            typename graph_traits< Graph >::vertex_descriptor vertex_descriptor;
        typedef typename Topology::point_type Point;
        typedef typename Topology::point_difference_type PointDiff;

        fr_apply_force(const Topology& topology, const PositionMap& position,
            const DisplacementMap& displacement, RepulsiveForce repulsive_force,
            double k, const Graph& g)
        : topology(topology)
        , position(position)
        , displacement(displacement)
        , repulsive_force(repulsive_force)
        , k(k)
        , g(g)
        {
        }

        void operator()(vertex_descriptor u, vertex_descriptor v)
        {
            if (u != v)
            {
                // When the vertices land on top of each other, move the
                // first vertex away from the boundaries.
                maybe_jitter_point(topology, position, u, get(position, v));

                double dist
                    = topology.distance(get(position, u), get(position, v));
                typename Topology::point_difference_type dispv
                    = get(displacement, v);
                if (dist == 0.)
                {
                    for (std::size_t i = 0; i < Point::dimensions; ++i)
                    {
                        dispv[i] += 0.01;
                    }
                }
                else
                {
                    double fr = repulsive_force(u, v, k, dist, g);
                    dispv += (fr / dist)
                        * topology.difference(
                            get(position, v), get(position, u));
                }
                put(displacement, v, dispv);
            }
        }

    private:
        const Topology& topology;
        PositionMap position;
        DisplacementMap displacement;
        RepulsiveForce repulsive_force;
        double k;
        const Graph& g;
    };

} // end namespace detail

template < typename Topology, typename Graph, typename PositionMap,
    typename AttractiveForce, typename RepulsiveForce, typename ForcePairs,
    typename Cooling, typename DisplacementMap >
void fruchterman_reingold_force_directed_layout(const Graph& g,
    PositionMap position, const Topology& topology,
    AttractiveForce attractive_force, RepulsiveForce repulsive_force,
    ForcePairs force_pairs, Cooling cool, DisplacementMap displacement)
{
    typedef typename graph_traits< Graph >::vertex_iterator vertex_iterator;
    typedef typename graph_traits< Graph >::vertex_descriptor vertex_descriptor;
    typedef typename graph_traits< Graph >::edge_iterator edge_iterator;

    double volume = topology.volume(topology.extent());

    // assume positions are initialized randomly
    double k = pow(volume / num_vertices(g),
        1. / (double)(Topology::point_difference_type::dimensions));

    detail::fr_apply_force< Topology, PositionMap, DisplacementMap,
        RepulsiveForce, Graph >
        apply_force(topology, position, displacement, repulsive_force, k, g);

    do
    {
        // Calculate repulsive forces
        vertex_iterator v, v_end;
        for (boost::tie(v, v_end) = vertices(g); v != v_end; ++v)
            put(displacement, *v, typename Topology::point_difference_type());
        force_pairs(g, apply_force);

        // Calculate attractive forces
        edge_iterator e, e_end;
        for (boost::tie(e, e_end) = edges(g); e != e_end; ++e)
        {
            vertex_descriptor v = source(*e, g);
            vertex_descriptor u = target(*e, g);

            // When the vertices land on top of each other, move the
            // first vertex away from the boundaries.
            ::boost::detail::maybe_jitter_point(
                topology, position, u, get(position, v));

            typename Topology::point_difference_type delta
                = topology.difference(get(position, v), get(position, u));
            double dist = topology.distance(get(position, u), get(position, v));
            double fa = attractive_force(*e, k, dist, g);

            put(displacement, v, get(displacement, v) - (fa / dist) * delta);
            put(displacement, u, get(displacement, u) + (fa / dist) * delta);
        }

        if (double temp = cool())
        {
            // Update positions
            BGL_FORALL_VERTICES_T(v, g, Graph)
            {
                BOOST_USING_STD_MIN();
                BOOST_USING_STD_MAX();
                double disp_size = topology.norm(get(displacement, v));
                put(position, v,
                    topology.adjust(get(position, v),
                        get(displacement, v)
                            * (min BOOST_PREVENT_MACRO_SUBSTITUTION(
                                   disp_size, temp)
                                / disp_size)));
                put(position, v, topology.bound(get(position, v)));
            }
        }
        else
        {
            break;
        }
    } while (true);
}

namespace detail
{
    template < typename DisplacementMap > struct fr_force_directed_layout
    {
        template < typename Topology, typename Graph, typename PositionMap,
            typename AttractiveForce, typename RepulsiveForce,
            typename ForcePairs, typename Cooling, typename Param, typename Tag,
            typename Rest >
        static void run(const Graph& g, PositionMap position,
            const Topology& topology, AttractiveForce attractive_force,
            RepulsiveForce repulsive_force, ForcePairs force_pairs,
            Cooling cool, DisplacementMap displacement,
            const bgl_named_params< Param, Tag, Rest >&)
        {
            fruchterman_reingold_force_directed_layout(g, position, topology,
                attractive_force, repulsive_force, force_pairs, cool,
                displacement);
        }
    };

    template <> struct fr_force_directed_layout< param_not_found >
    {
        template < typename Topology, typename Graph, typename PositionMap,
            typename AttractiveForce, typename RepulsiveForce,
            typename ForcePairs, typename Cooling, typename Param, typename Tag,
            typename Rest >
        static void run(const Graph& g, PositionMap position,
            const Topology& topology, AttractiveForce attractive_force,
            RepulsiveForce repulsive_force, ForcePairs force_pairs,
            Cooling cool, param_not_found,
            const bgl_named_params< Param, Tag, Rest >& params)
        {
            typedef typename Topology::point_difference_type PointDiff;
            std::vector< PointDiff > displacements(num_vertices(g));
            fruchterman_reingold_force_directed_layout(g, position, topology,
                attractive_force, repulsive_force, force_pairs, cool,
                make_iterator_property_map(displacements.begin(),
                    choose_const_pmap(
                        get_param(params, vertex_index), g, vertex_index),
                    PointDiff()));
        }
    };

} // end namespace detail

template < typename Topology, typename Graph, typename PositionMap,
    typename Param, typename Tag, typename Rest >
void fruchterman_reingold_force_directed_layout(const Graph& g,
    PositionMap position, const Topology& topology,
    const bgl_named_params< Param, Tag, Rest >& params)
{
    typedef typename get_param_type< vertex_displacement_t,
        bgl_named_params< Param, Tag, Rest > >::type D;

    detail::fr_force_directed_layout< D >::run(g, position, topology,
        choose_param(get_param(params, attractive_force_t()),
            square_distance_attractive_force()),
        choose_param(get_param(params, repulsive_force_t()),
            square_distance_repulsive_force()),
        choose_param(get_param(params, force_pairs_t()),
            make_grid_force_pairs(topology, position, g)),
        choose_param(
            get_param(params, cooling_t()), linear_cooling< double >(100)),
        get_param(params, vertex_displacement_t()), params);
}

template < typename Topology, typename Graph, typename PositionMap >
void fruchterman_reingold_force_directed_layout(
    const Graph& g, PositionMap position, const Topology& topology)
{
    fruchterman_reingold_force_directed_layout(g, position, topology,
        attractive_force(square_distance_attractive_force()));
}

} // end namespace boost

#include BOOST_GRAPH_MPI_INCLUDE(<boost/graph/distributed/fruchterman_reingold.hpp>)

#endif // BOOST_GRAPH_FRUCHTERMAN_REINGOLD_FORCE_DIRECTED_LAYOUT_HPP