#ifndef _MODULARITY_HPP_
#define _MODULARITY_HPP_

#include <boost/graph/iteration_macros.hpp>



//Just an ugly hack to make eclipse belief that these macros exist.
//I don't know where it is looking for them, but it certainly isn't in iteration_macros.hpp
#ifndef BGL_CAT
#define BGL_CAT(x,y) x ## y
#define BGL_RANGE(linenum) BGL_CAT(bgl_range_,linenum)
#define BGL_FIRST(linenum) (BGL_RANGE(linenum).first)
#define BGL_LAST(linenum) (BGL_RANGE(linenum).second)
#endif

#include <sferes/dbg/dbg.hpp>

#include <boost/graph/graph_traits.hpp>
#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/graphviz.hpp>
#include <boost/foreach.hpp>

#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>

#include "igraph.hpp"

#define MOD_FINETUNING

// Community Structure in Directed Networks
// EA Leicht, MEJ Newman - Physical Review Letters, 2008 - APS
// http://arxiv.org/pdf/0709.4500
namespace mod {
namespace null_model {
// the null model used by Leicht
struct Directed {
    template<typename G>
    double operator()(const G& g,
            const typename G::vertex_descriptor& v1,
            const typename G::vertex_descriptor& v2) const {
        size_t m = num_edges(g);
        size_t k_i_in = in_degree(v1, g);
        size_t k_j_out = out_degree(v2, g);
        if(m==0){
            std::cout << "(Directed): The number of edges  is 0" << std::endl;
            std::cout << "The returned value is " << k_i_in * k_j_out / (double) m << std::endl;
        }
        dbg::assertion(DBG_ASSERTION(m > 0));
        return k_i_in * k_j_out / (double) m;
    }
};

// the null model used by Newman
struct Undirected {
    template<typename G>
    double operator()(const G& g,
            const typename G::vertex_descriptor& v1,
            const typename G::vertex_descriptor& v2) const {
        size_t m = num_edges(g);
        size_t k_i_in = degree(v1, g);
        size_t k_j_out = degree(v2, g);
        dbg::assertion(DBG_ASSERTION(m > 0));
        return k_i_in * k_j_out / (double) m;
    }
};

// for layered networks (e.g. multi-layer perceptron)
// the map contains the id of the layer for each node
// WARNING: the last layer is 0!
template<typename G>
struct Layered {
    typedef std::map<typename G::vertex_descriptor, int> map_t;
    typedef typename map_t::const_iterator map_it_t;
    Layered(const G& g, const map_t & m) : _m(m) {
        // find the max layer
        int cm = -1;
        for (typename map_t::const_iterator it = m.begin(); it != m.end(); ++it)
            if (it->second > cm)
                cm = it->second;
        _layers.resize(cm + 1);
        for (typename map_t::const_iterator it = m.begin(); it != m.end(); ++it)
            _layers[it->second].push_back(it->first);
        _num_edges.resize(_layers.size());
        for (size_t i = 0; i < _layers.size(); ++i) {
            int k = 0;
            for (size_t j = 0; j < _layers[i].size(); ++j)
                k += out_degree(_layers[i][j], g);
            _num_edges[i] = k;
        }
    }

    double operator()(const G& g,
            const typename G::vertex_descriptor& v1,
            const typename G::vertex_descriptor& v2) const {
//        for (map_it_t it=_m.begin(); it!=_m.end(); ++it){
//            std::cout << it->first << " => " << it->second << '\n';
//        }
        int k1 = _m.find(v1)->second;
        int k2 = _m.find(v2)->second;
        dbg::out(dbg::info, "null_model") << "v1: " << g[v1].get_id()  << " k1: " << k1  << " v2: " << g[v2].get_id()  << " k2: " << k2 << std::endl;
        dbg::assertion(DBG_ASSERTION(_m.find(v1) != _m.end()));
        dbg::assertion(DBG_ASSERTION(_m.find(v2) != _m.end()));

        if (k1 <= k2){
            return 0;
        }

        //JH: More strict version
        if ((k1 - k2) != 1){
            return 0;
        }

        size_t m = _num_edges[k1];
        size_t k_j_out = out_degree(v1, g);
        size_t k_i_in = in_degree(v2, g);

        dbg::out(dbg::info, "null_model") << "m: " << m << " k_j_out: " << k_j_out << " k_i_in: " << k_i_in << " P: " << k_i_in * k_j_out / (double) m << std::endl;
        dbg::assertion(DBG_ASSERTION(m > 0));
        return k_i_in * k_j_out / (double) m;
    }
    map_t _m;
    std::vector<std::vector<typename G::vertex_descriptor> > _layers;
    std::vector<int> _num_edges;
};
}

template<typename G1, typename G2>
void convert_graph(const G1& src, G2& g,
        std::map <
        typename boost::graph_traits<G2>::vertex_descriptor,
        typename boost::graph_traits<G1>::vertex_descriptor > & rmap) {
    using namespace boost;
    typedef typename graph_traits<G2>::vertex_descriptor vertex_desc_t;
    std::map<typename graph_traits<G1>::vertex_descriptor, vertex_desc_t> neuron_map;

    BGL_FORALL_VERTICES_T(v, src, G1) {
        vertex_desc_t d = add_vertex(g);
        neuron_map[v] = d;
        rmap[d] = v;
    }

    BGL_FORALL_EDGES_T(v, src, G1)
    add_edge(neuron_map[source(v, src)],
            neuron_map[target(v, src)],
            g);
}

template<typename G>
inline bool _linked(const G& g,
        const typename G::vertex_descriptor& v1,
        const typename G::vertex_descriptor& v2) {
    BGL_FORALL_OUTEDGES_T(v1, e, g, G)
            if (target(e, g) == v2)
                return true;
    return false;
}



template<typename G, typename NullModel>
inline void _b_matrix(const G& g, gsl_matrix *& b, const NullModel& null_model) {
    dbg::trace trace("modularity" , DBG_HERE);
    size_t nb_vertices = num_vertices(g);
    dbg::out(dbg::info, "modularity") << "Nb of vertices: " << nb_vertices << std::endl;
    b = gsl_matrix_alloc(nb_vertices, nb_vertices);

//    size_t m = num_edges(g);
    size_t i = 0;
    BGL_FORALL_VERTICES_T(v1, g, G) {
        size_t j = 0;
        BGL_FORALL_VERTICES_T(v2, g, G) {
            size_t k_i_in = in_degree(v1, g);
            size_t k_j_out = out_degree(v2, g);
            double b_ij = _linked(g, v1, v2) - null_model(g, v1, v2);
            gsl_matrix_set(b, i, j, b_ij);
            ++j;
        }
        ++i;
    }
}

// res = b^T + b, allocate res
static void _bbt(gsl_matrix*b, gsl_matrix *& res) {
    gsl_matrix*bt = gsl_matrix_alloc(b->size1, b->size2);
    res = gsl_matrix_alloc(b->size1, b->size2);
    gsl_matrix_memcpy(res, b);

    gsl_matrix_transpose_memcpy(bt, b);
    gsl_matrix_add(res, bt);
    gsl_matrix_free(bt);
}

// split in two according to b
static void _split_eigen(gsl_matrix*b, size_t n, std::vector<int>& s) {
    assert(n);
    s.resize(n);
    gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc(n);
    gsl_vector *eval = gsl_vector_alloc(n);
    gsl_matrix *evec = gsl_matrix_alloc(n, n);
    gsl_matrix *b_copy = gsl_matrix_alloc(n, n);
    gsl_matrix_memcpy(b_copy, b);
    gsl_eigen_symmv(b_copy, eval, evec, w);

    gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_VAL_DESC);
    for (size_t i = 0; i < n; ++i)
        s[i] = (gsl_matrix_get(evec, i, 0) > 0 ? 1 : -1);

    gsl_eigen_symmv_free(w);
    gsl_vector_free(eval);
    gsl_matrix_free(evec);
    gsl_matrix_free(b_copy);
    assert(s.size() == n);
}

static void _bg(gsl_matrix*b, const std::vector<size_t>& sg1,
        gsl_matrix *& bg) {
    size_t n1 = sg1.size();
    assert(n1);
    bg = gsl_matrix_alloc(n1, n1);
    for (size_t i = 0; i < n1; ++i)
        for (size_t j = 0; j < n1; ++j) {
            double b_ij = gsl_matrix_get(b, sg1[i], sg1[j]);
            if (i == j)
                for (size_t k = 0; k < n1; ++k)
                    b_ij -= gsl_matrix_get(b, sg1[i], sg1[k]);
            gsl_matrix_set(bg, i, j, b_ij);
        }
}


double _delta_q(gsl_matrix*bbt, const std::vector<int>& sg, size_t m) {
    // 1/(4 * m) * s^T * bbt * s
    assert(sg.size());
    gsl_matrix*sv = gsl_matrix_alloc(sg.size(), 1);
    gsl_matrix*svt = gsl_matrix_alloc(1, sg.size());
    gsl_matrix*tmp1 = gsl_matrix_alloc(sg.size(), 1);
    gsl_matrix*tmp2 = gsl_matrix_alloc(1, 1);

    for (size_t i = 0; i < sg.size(); ++i) {
        gsl_matrix_set(sv, i, 0, sg[i]);
        gsl_matrix_set(svt, 0, i, sg[i]);
    }
    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans,
            1.0, bbt, sv,
            0.0, tmp1);
    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans,
            1.0, svt, tmp1,
            0.0, tmp2);
    double delta_q = 1.0 / (4.0 * m) * gsl_matrix_get(tmp2, 0, 0);

    gsl_matrix_free(sv);
    gsl_matrix_free(svt);
    gsl_matrix_free(tmp1);
    gsl_matrix_free(tmp2);

    return delta_q;
}

template<typename G, typename NullModel>
float compute_modularity(const G& g,
        const std::vector<std::string>& modules,
        const NullModel& null_model)
{
    double m = boost::num_edges(g);
    if(m == 0) return 0;
    double q = 0.0;
    size_t v1_i = 0;
    BGL_FORALL_VERTICES_T(v1, g, G) {
        size_t v2_i = 0;
        BGL_FORALL_VERTICES_T(v2, g, G) {
            if (modules[v1_i] == modules[v2_i]) {
                double b_ij = _linked(g, v1, v2) - null_model(g, v1, v2);
                q += b_ij;
            }
            ++v2_i;
        }
        ++v1_i;
    }
    return q / m;
}

/**
 * Finds the best (or least detrimental) vertices to switch.
 *
 * @param g           (input)        The boost graph.
 * @param null_model  (input)        The null model used for calculating modularity.
 * @param modules     (input/output) The tentative modular decomposition.
 * @param to_be_moved (input/output) A set containing all indices that should be considered for moving.
 *                                   The index of the vertex that is moved is removed from this set when
 *                                   the function returns.
 */
template<typename G, typename NullModel>
inline float _move_best(const G& g, std::vector<std::string>& modules, const NullModel& null_model, std::set<size_t>& to_be_moved) {
    size_t depth = modules[size_t(*to_be_moved.begin())].size()-1;
    float max_mod = -2;
    float current_mod = 0;
    size_t max_mod_index = 0;
    char old_module;
    std::set<size_t>::iterator it;
    for (it = to_be_moved.begin(); it != to_be_moved.end(); ++it){
        size_t i = *it; // Note the "*" here

        //Store previous module
        old_module = modules[i][depth];

        //Move vertices from one module to the other module
        if(modules[i][depth] == '0'){
            modules[i][depth] = '1';
        } else {
            modules[i][depth] = '0';
        }

        //Calculate the modularity of this new split
        current_mod = compute_modularity(g, modules, null_model);

        //Record the node that has the maximum impact on modularity
        if(current_mod > max_mod){
            max_mod = current_mod;
            max_mod_index = i;
        }

        //Restore what we have done
        modules[i][depth] = old_module;
    }

    //Move the vertex that leads to maximum modularity
    if(modules[max_mod_index][depth] == '0'){
        modules[max_mod_index][depth] = '1';
    } else {
        modules[max_mod_index][depth] = '0';
    }

    //Remove the vertex index from the set
    to_be_moved.erase(max_mod_index);

    return max_mod;
}

/**
 * Apply the fine-tuning process to the provided modular split.
 *
 * @param g          (input)        The boost graph.
 * @param null_model (input)        The null model used for calculating modularity.
 * @param modules    (input/output) The tentative modular decomposition.
 * @param sg         (input)        A vector of indices for the sub-group currently being split.
 * @param sg1        (output)       A vector of indices for the first sub-sub-group after fine-tuning.
 * @param sg2        (output)       A vector of indices for the second sub-sub-group after fine-tuning.
 */
template<typename G, typename NullModel>
void _finetune(const G& g, const NullModel& null_model,
        std::vector<std::string>& modules, const std::vector<size_t>& sg,
        std::vector<size_t>& sg1, std::vector<size_t>& sg2){
    //JH: Horribly inefficient implementation for now,
    //but we just need to see if this works
    float best_mod;
    float prev_mod;
    float temp_mod;
    std::set<size_t> to_be_moved;
    std::vector<std::string> best_modules;
    do {

        //Store the current modularity of the network
        prev_mod = compute_modularity(g, modules, null_model);
        best_mod = prev_mod;
        best_modules = modules;

        //Create a new set including all vertex indices.
        to_be_moved.clear();
        for (size_t i = 0; i < sg.size(); ++i){
            to_be_moved.insert(sg[i]);
        }

        //Find the best split possible
        for (size_t i = 0; i < sg.size(); ++i){
            temp_mod = _move_best(g, modules, null_model, to_be_moved);
            dbg::out(dbg::info, "modularity") << "  sub-split; temp_mod: " << temp_mod << std::endl;
            dbg::out(dbg::info, "modularity") << "  sub-split; comp_mod: " << compute_modularity(g, modules, null_model) << std::endl;
            if(temp_mod > best_mod){
                best_modules = modules;
                best_mod = temp_mod;
            }
        }

        //If the best split is better than the previous split
        //make it the current split and repeat.
        dbg::out(dbg::info, "modularity") << "sub-split; best mod: " << best_mod <<" prev mod: " << prev_mod << std::endl;
//        if(best_mod > prev_mod){
            modules = best_modules;
//        }


    } while(best_mod > prev_mod);
    for(size_t i=0; i<modules.size(); ++i){
        dbg::out(dbg::info, "modularity") << "sub-split; Modules: " << modules[i] << std::endl;
    }
    //Copy the new modules into the current split
    sg1.clear();
    sg2.clear();
    for (size_t i = 0; i < sg.size(); ++i){
        if (modules[sg[i]][modules[sg[i]].size()-1] == '0') {
            sg1.push_back(sg[i]);
        } else {
            sg2.push_back(sg[i]);
        }
    }
}


/**
 * Forward declaration to allow for recursion.
 */
template<typename G, typename NullModel>
void _split_sub_group_fine(const G& g, const NullModel& null_model, gsl_matrix*b, const std::vector<size_t>& sg, size_t m, std::vector<std::string>& modules);

/**
 * Recursively split the network, includes fine-tuning.
 *
 * @param g          (input)        The boost graph.
 * @param null_model (input)        The null model used for calculating modularity.
 * @param b          (input)        The modularity matrix, calculated once.
 * @param sg         (input)        A vector of indices for the sub-group currently being split.
 * @param m          (input)        The number of vertices of the network.
 * @param modules    (input/output) The best modular decomposition before the previous step.
 * @param s          (input)        The best split found through the Eigenvector method in the previous step.
 */
template<typename G, typename NullModel>
void _split_fine_req(const G& g, const NullModel& null_model, gsl_matrix*b, const std::vector<size_t>& sg, size_t m,
        std::vector<std::string>& modules, const std::vector<int>& s) {
    //Define modules in a way that the compute modularity method understands.
    assert(sg.size() == s.size());
    for (size_t i = 0; i < s.size(); ++i){
        if (s[i] == -1) {
            modules[sg[i]] += '0';
        } else {
            modules[sg[i]] += '1';
        }
    }

    //Perform fine-tuning of the sub-groups
    std::vector<size_t> sg1, sg2;
    _finetune(g, null_model, modules, sg, sg1, sg2);

    //Recursively call split for sub-groups
    if (sg1.size() > 0 && sg2.size() > 0) {
        _split_sub_group_fine(g, null_model, b, sg1, m, modules);
        _split_sub_group_fine(g, null_model, b, sg2, m, modules);
    }
}

/**
 * Recursively split sub-group, including fine-tuning.
 *
 * @param g          (input)        The boost graph.
 * @param null_model (input)        The null model used for calculating modularity.
 * @param b          (input)        The modularity matrix, calculated once.
 * @param sg         (input)        A vector of indices for the sub-group currently being split.
 * @param m          (input)        The number of vertices of the network.
 * @param modules    (input/output) The best modular decomposition found thus far.
 */
template<typename G, typename NullModel>
void _split_sub_group_fine(const G& g, const NullModel& null_model, gsl_matrix*b, const std::vector<size_t>& sg, size_t m,
        std::vector<std::string>& modules) {
    assert(sg.size());

    //Declare necessary pointers
    gsl_matrix*bbt, *bg;

    //Calculate modularity matrix for the sub-network
    _bg(b, sg, bg);

    //Calculate the optimal split for the sub-network
    _bbt(bg, bbt);
    std::vector<int> s;
    _split_eigen(bbt, sg.size(), s);

    //Calculate the delta for the sub-network
    double delta_q = _delta_q(bbt, s, m);

    //Cleanup
    gsl_matrix_free(bbt);
    gsl_matrix_free(bg);

    //If our delta is below the threshold, return
    if (delta_q <= 1e-5) return;

    //Otherwise, continue with recursion
    _split_fine_req(g, null_model, b, sg, m, modules, s);
}


/**
 * Calculates the optimal split for a network, including a fine-tuning process.
 *
 * Calculates the optimal split for a network based on Newman's Eigenvector method
 * (also know as the spectral method). Also includes an optimization in the form
 * of a fine-tuning process. The fine-tuning process works as follows:
 *
 * First, a tentative split into two modules is determined through the Eigenvector method.
 * Second, we determine the vertice that, if it is moved to the other module, will lead
 * to the greatest increase (or lowest decrease) in modularity. Now, excluding this first
 * vertice, we determine the second vertice that has the most beneficial (or least
 * detrimental) effect on modularity. We do this until all nodes have moved exactly once.
 * Third, if, during the above process, the modularity of the network ever became greater
 * than the starting modularity, we assign vertices according to the best distribution
 * encountered, and we repeat step 2. Otherwise, we do not change the modular decomposition
 * of the network, and exit.
 *
 * @param g          (input)  The boost graph.
 * @param modules    (output) The best modular decomposition found thus far.
 * @param null_model (input)  The null model used for calculating modularity.
 */
template<typename G, typename NullModel>
inline void split_fine(const G& g, std::vector<std::string>& modules, const NullModel& null_model) {
    dbg::trace trace("modularity" , DBG_HERE);

    //Get number of edges and vertices and define the necessary pointers
    gsl_matrix *b, *bbt;
    size_t m = num_edges(g);
    size_t n = num_vertices(g);
    modules.resize(n);
    dbg::out(dbg::info, "modularity") << "Num vertices: " << n << " Num edges: " << m << std::endl;

    //If there are no vertices, nothing needs to be done.
    if(n == 0) return;

    //If there are no connections in the network, there does not exist a sensible split
    //Just put all nodes in the same module and call it a day.
    if(m==0){
        for (size_t i = 0; i < n; ++i){
            modules[i] += '0';
        }
        return;
    }

    //Initialize a sub-group vector. Initially, all vertices are in the same group.
    std::vector<size_t> sg(n);
    for (size_t i = 0; i < n; ++i){
        sg[i] = i;
    }

    //Compute the modularity matrix
    _b_matrix(g, b, null_model);

    //Compute the optimal split according to the Eigen-vector method.
    _bbt(b, bbt);
    std::vector<int> s;
    _split_eigen(bbt, n, s);

    //Recursively split the network
    _split_fine_req(g, null_model, b, sg, m, modules, s);

    //Cleanup
    gsl_matrix_free(bbt);
    gsl_matrix_free(b);
}

template<typename G, typename NullModel>
void _split(const G& g, const NullModel& null_model, gsl_matrix*b, const std::vector<size_t>& sg, size_t m,
        std::vector<std::string>& modules) {
    assert(sg.size());

    gsl_matrix*bbt, *bg;
    _bg(b, sg, bg);
    _bbt(bg, bbt);
    std::vector<int> s;
    _split_eigen(bbt, sg.size(), s);

    double delta_q = _delta_q(bbt, s, m);
    if (delta_q <= 1e-5) {
        gsl_matrix_free(bbt);
        gsl_matrix_free(bg);
        return;
    }

    // new sg
    std::vector<size_t> sg1, sg2;
    assert(sg.size() == s.size());
    for (size_t i = 0; i < s.size(); ++i){
        if (s[i] == -1) {
            modules[sg[i]] += '0';
            sg1.push_back(sg[i]);
        } else {
            modules[sg[i]] += '1';
            sg2.push_back(sg[i]);
        }
    }

    gsl_matrix_free(bbt);
    gsl_matrix_free(bg);

    if (sg1.size() > 0 && sg2.size() > 0) {
        _split(g, null_model, b, sg1, m, modules);
        _split(g, null_model, b, sg2, m, modules);
    }

}

template<typename G, typename NullModel>
inline void split(const G& g, std::vector<std::string>& modules, const NullModel& null_model) {
    dbg::trace trace("modularity" , DBG_HERE);
    dbg::out(dbg::info, "modularity") << std::endl;
    dbg::out(dbg::info, "modularity") << "**************************************" << std::endl;
    dbg::out(dbg::info, "modularity") << "********** SPLIT START      **********" << std::endl;
    dbg::out(dbg::info, "modularity") << "**************************************" << std::endl;

    size_t m = num_edges(g);
    size_t n = num_vertices(g);
    modules.resize(num_vertices(g));
    if(n == 0) return;
    //If there are no connections in the network, there does not exist a sensible split
    //Just put all nodes in the same module and call it a day.
    if(m==0){
        for (size_t i = 0; i < n; ++i){
            modules[i] += '0';
        }
        return;
    }

    dbg::out(dbg::info, "modularity") << "Num edges: " << m << std::endl;
    gsl_matrix*b, *bbt;

    _b_matrix(g, b, null_model);
    _bbt(b, bbt);


    std::vector<int> s;
    _split_eigen(bbt, n, s);

    size_t n1 = std::count(s.begin(), s.end(), -1);
    size_t n2 = n - n1;

    std::vector<size_t> sg1(n1), sg2(n2);
    size_t j1 = 0, j2 = 0;
    for (size_t i = 0; i < n; ++i){
        if (s[i] == -1) {
            modules[i] += '0';
            assert(j1 < sg1.size());
            sg1[j1] = i;
            ++j1;
        } else {
            modules[i] += '1';
            assert(j2 < sg2.size());
            sg2[j2] = i;
            ++j2;
        }
    }

    if (sg1.size() > 0 && sg2.size() > 0) {
        _split(g, null_model, b, sg1, m, modules);
        _split(g, null_model, b, sg2, m, modules);
    }

    gsl_matrix_free(b);
    gsl_matrix_free(bbt);
}

template<typename G>
void _h_modules(const std::set<typename boost::graph_traits<G>::vertex_descriptor>& nodes,
        const std::vector<std::string>& ids,
        std::vector<std::set<typename boost::graph_traits<G>::vertex_descriptor> >& mods,
        size_t n) {
    assert(nodes.size());
    typedef typename boost::graph_traits<G>::vertex_descriptor vertex_desc_t;
    std::set<vertex_desc_t> s1, s2;
    BOOST_FOREACH(vertex_desc_t v, nodes) {
        if (n >= ids[v].size()) {
            assert(s1.empty());
            assert(s2.empty());
            return;
        }

        if (ids[v][n] == '0')
            s1.insert(v);
        else
            s2.insert(v);
    }
    mods.push_back(s1);
    mods.push_back(s2);
    assert(s1.size());
    assert(s2.size());
    if (s1.size())
        _h_modules<G>(s1, ids, mods, n + 1);
    if (s2.size())
        _h_modules<G>(s2, ids, mods, n + 1);
}

// ALL the modules (high and low levels)
template<typename G, typename NullModel>
double h_modules(const G& g_origin,
        std::vector<std::set<typename boost::graph_traits<G>::vertex_descriptor> >& mods_r,
        const NullModel& null_model) {
    using namespace boost;
    typedef adjacency_list<vecS, vecS, bidirectionalS> graph_t;
    typedef graph_traits<graph_t>::vertex_descriptor vertex_desc_t;

    graph_t g;
    std::vector<std::set<vertex_desc_t> > mods;
    std::map<vertex_desc_t, typename boost::graph_traits<G>::vertex_descriptor> rmap;
    convert_graph(g_origin, g, rmap);

    assert(num_vertices(g) != 0);
    std::vector<std::string> ids;
    split(g, ids, null_model);
    typedef std::set<vertex_desc_t> vset_t;
    vset_t nodes;

    BGL_FORALL_VERTICES_T(v, g, graph_t)
    nodes.insert(v);

    assert(nodes.size());
    _h_modules<graph_t>(nodes, ids, mods, 0);

    float m = compute_modularity(g, ids, null_model);
    // create mod_r
    mods_r.resize(mods.size());
    for (size_t i = 0; i < mods.size(); ++i)
        BOOST_FOREACH(vertex_desc_t v, mods[i])
        mods_r[i].insert(rmap[v]);
    return m;
}

// true if set is included in Vm
template<typename S>
bool _includes(const S& s, const std::vector<S>& vm) {
    for (size_t i = 0; i < vm.size(); ++i)
        if (s != vm[i] &&
                std::includes(vm[i].begin(), vm[i].end(), s.begin(), s.end()))
            return true;
    return false;
}

template<typename G, typename Vm>
void _write_modules(const G& g, const Vm& vm, std::ostream& os) {
    Vm vv;
    // find the sets that are not included in another set
    for (size_t i = 0; i < vm.size(); ++i)
        if (!_includes(vm[i], vm))
            vv.push_back(vm[i]);
    // split the vm list
    std::vector<Vm> n_vm(vv.size());
    for (size_t k = 0; k < vv.size(); ++k)
        for (size_t i = 0; i < vm.size(); ++i)
            if (vm[i] != vv[k] &&
                    std::includes(vv[k].begin(), vv[k].end(),
                            vm[i].begin(), vm[i].end()))
                n_vm[k].push_back(vm[i]);

    // recurse
    for (size_t i = 0; i < n_vm.size(); ++i) {
        os << "subgraph cluster_" << rand() << " {" << std::endl;
        if (!n_vm[i].empty())
            _write_modules(g, n_vm[i], os);
        else
            for (typename Vm::value_type::const_iterator it = vv[i].begin();
                    it != vv[i].end(); ++it)
                os << g[*it]._id << std::endl;
        os << "}" << std::endl;
    }
}


// hierarchical dot file (slow version)
template<typename G, typename NullModel>
void write_modules(const G& g, std::ostream& os, const NullModel& null_model) {
    typedef std::vector<std::set<typename boost::graph_traits<G>::vertex_descriptor> >
    v_mod_t;
    v_mod_t mods_r;
    h_modules(g, mods_r, null_model);
    os << "digraph G {" << std::endl;
    _write_modules(g, mods_r, os);
    BGL_FORALL_EDGES_T(e, g, G) {
        os << g[source(e, g)]._id
                << " -> " << g[target(e, g)]._id
                << "[label=\"" << g[e].get_weight() << "\"]" << std::endl;
    }

    os << "}" << std::endl;
}


// post-processing : modules in a more convenient form (list of sets)
// ONLY leaf-modules (not the 'higher-level' modules')
// TODO : convert graph !
template<typename G, typename NullModel>
inline double modules(
        const G& g,
        std::vector<std::vector<typename boost::graph_traits<G>::vertex_descriptor> >& mods,
        const NullModel& null_model) {
    dbg::trace trace("modularity" , DBG_HERE);
    std::vector<std::string> modules;
    split(g, modules, null_model);
    dbg::out(dbg::info, "modularity") << "Nb of modules: " << modules.size() << std::endl;
    dbg::out(dbg::info, "modularity") << "Nb of nodes: " << num_vertices(g) << std::endl;
    dbg::assertion(DBG_ASSERTION(modules.size() == num_vertices(g)));

    typedef typename boost::graph_traits<G>::vertex_descriptor vertex_desc_t;
    typedef std::map<std::string, std::vector<vertex_desc_t> > m_map_t;
    typedef typename m_map_t::const_iterator it_t;
    m_map_t m_map;
    size_t i = 0;
    BGL_FORALL_VERTICES_T(v, g, G){
        m_map[modules[i++]].push_back(v);
    }

    // copy to res
    mods.resize(m_map.size());
    size_t k = 0;
    for (it_t it = m_map.begin(); it != m_map.end(); ++it) {
        BOOST_FOREACH(vertex_desc_t v, it->second){
              mods[k].push_back(v);
        }
        ++k;
    }

    float q = compute_modularity(g, modules, null_model);
    return q;
}


template<typename G, typename NullModel>
inline double modularity(const G& g, const NullModel& null_model) {
    dbg::trace trace("modularity" , DBG_HERE);
    std::vector<std::vector<typename boost::graph_traits<G>::vertex_descriptor> > mods;
    return modules(g, mods, null_model);
}



template<typename G, typename NullModel>
std::pair<double, double> random(const G& g, const NullModel& null_model,
        size_t k = 1000, size_t k2 = 50) {
    using namespace boost;
    typedef typename boost::graph_traits<G>::edge_descriptor edge_desc_t;
    typedef typename boost::graph_traits<G>::vertex_descriptor vertex_desc_t;

    double m = 0;
    double mm = 0;
    for (size_t i = 0; i < k; ++i) {
        if (i % 100 == 0) {
            std::cout << "[" << i << "] ";
            std::cout.flush();
        }
        G rg = g;
        for (size_t j = 0; j < num_vertices(rg) * k2; ++j) {
            edge_desc_t e1 = igraph::random_edge(rg);
            edge_desc_t e2 = igraph::random_edge(rg);
            if (e1 == e2)
                continue;
            vertex_desc_t src_1 = source(e1, rg);
            vertex_desc_t src_2 = source(e2, rg);
            vertex_desc_t tgt_1 = target(e1, rg);
            vertex_desc_t tgt_2 = target(e2, rg);
            remove_edge(e1, rg);
            remove_edge(e2, rg);
            add_edge(src_1, tgt_2, rg);
            add_edge(src_2, tgt_1, rg);
        }
        double q = modularity(rg, null_model);
        m += q;
        mm = std::max(mm, q);
    }
    return std::make_pair(m / k, mm);
}


// the source network (g) has to be layered!
template<typename G>
std::pair<double, double> random_feedforward(const G& g, size_t k = 1000, size_t k2 = 50) {
    using namespace boost;
    typedef typename boost::graph_traits<G>::edge_descriptor edge_desc_t;
    typedef typename boost::graph_traits<G>::vertex_descriptor vertex_desc_t;
    typedef std::map<typename G::vertex_descriptor, int> map_t;
    // reconstruct layers
    map_t layer_map_g = compute_layers(g);
    int cm = -1;
    for (typename map_t::const_iterator it = layer_map_g.begin(); it != layer_map_g.end(); ++it)
        if (it->second > cm)
            cm = it->second;
    std::vector<std::vector<edge_desc_t> > edge_list(cm);
    double m = 0;
    double mm = 0;

    for (size_t i = 0; i < k; ++i) {
        if (i % 100 == 0) {
            std::cout << "[" << i << "] ";
            std::cout.flush();
        }
        G rg = g;
        map_t layer_map = compute_layers(rg);

        for (size_t j = 0; j < num_vertices(rg) * k; ++j) {
            for (size_t t = 0; t < edge_list.size(); ++t)
                edge_list[t].clear();
            BGL_FORALL_EDGES_T(e, rg, G)
            edge_list[cm - layer_map[source(e, rg)]].push_back(e);
            int n = igraph::misc::rand<int>(0, edge_list.size());
            int r1 = igraph::misc::rand<int>(0, edge_list[n].size());
            int r2 = igraph::misc::rand<int>(0, edge_list[n].size());
            edge_desc_t e1 = edge_list[n][r1];
            edge_desc_t e2 = edge_list[n][r2];
            if (e1 == e2)
                continue;
            vertex_desc_t src_1 = source(e1, rg);
            vertex_desc_t src_2 = source(e2, rg);
            vertex_desc_t tgt_1 = target(e1, rg);
            vertex_desc_t tgt_2 = target(e2, rg);
            remove_edge(e1, rg);
            remove_edge(e2, rg);
            add_edge(src_1, tgt_2, rg);
            add_edge(src_2, tgt_1, rg);
        }
        double q = modularity(rg, null_model::Layered<igraph::graph_t>(rg, layer_map));
        m += q;
        mm = std::max(mm, q);
    }
    return std::make_pair(m / k, mm);
}

}

#endif
