/*
   igraph library.
   Copyright (C) 2007-2020 The igraph development team

   This program is free software; you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation; either version 2 of the License, or
   (at your option) any later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software
   Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
   02110-1301 USA

*/

#include "igraph_community.h"

#include "igraph_interface.h"
#include "igraph_structural.h"

#include "community/community_internal.h"

/**
 * \function igraph_modularity
 * \brief Calculates the modularity of a graph with respect to some clusters or vertex types.
 *
 * The modularity of a graph with respect to some clustering of the vertices
 * (or assignment of vertex types)
 * measures how strongly separated the different clusters are from each
 * other compared to a random null model. It is defined as
 *
 * </para><para>
 * <code>Q = 1/(2m) sum_ij (A_ij - γ k_i k_j / (2m)) δ(c_i,c_j)</code>,
 *
 * </para><para>
 * where \c m is the number of edges, <code>A_ij</code> is the adjacency matrix,
 * \c k_i is the degree of vertex \c i, \c c_i is the cluster that vertex \c i belongs to
 * (or its vertex type), <code>δ(i,j)=1</code> if <code>i=j</code> and 0 otherwise,
 * and the sum goes over all \c i, \c j pairs of vertices. Note that in this formula,
 * the diagonal of the adjacency matrix contains twice the number of self-loops.
 *
 * </para><para>
 * The resolution parameter \c γ allows weighting the random null model, which
 * might be useful when finding partitions with a high modularity. Maximizing modularity
 * with higher values of the resolution parameter typically results in more, smaller clusters
 * when finding partitions with a high modularity. Lower values typically results in
 * fewer, larger clusters. The original definition of modularity is retrieved
 * when setting <code>γ = 1</code>.
 *
 * </para><para>
 * Modularity can also be calculated on directed graphs. This only requires a relatively
 * modest change,
 *
 * </para><para>
 * <code>Q = 1/m sum_ij (A_ij - γ k^out_i k^in_j / m) δ(c_i,c_j)</code>,
 *
 * </para><para>
 * where \c k^out_i is the out-degree of node \c i and \c k^in_j is the in-degree of node \c j.
 *
 * </para><para>
 * Modularity on weighted graphs is also meaningful. When taking
 * edge weights into account, \c A_ij equals the weight of the corresponding edge
 * (or 0 if there is no edge), \c k_i is the strength (i.e. the weighted degree) of
 * vertex \c i, with similar counterparts for a directed graph, and \c m is the total
 * weight of all edges.
 *
 * </para><para>
 * Note that the modularity is not well-defined for graphs with no edges.
 * igraph returns \c NaN for graphs with no edges; see
 * https://github.com/igraph/igraph/issues/1539 for
 * a detailed discussion.
 *
 * </para><para>
 * For the original definition of modularity, see Newman, M. E. J., and Girvan, M.
 * (2004). Finding and evaluating community structure in networks.
 * Physical Review E 69, 026113. https://doi.org/10.1103/PhysRevE.69.026113
 *
 * </para><para>
 * For the directed definition of modularity, see Leicht, E. A., and Newman, M. E.
 * J. (2008). Community Structure in Directed Networks. Physical Review Letters 100,
 * 118703. https://doi.org/10.1103/PhysRevLett.100.118703
 *
 * </para><para>
 * For the introduction of the resolution parameter \c γ, see Reichardt, J., and
 * Bornholdt, S. (2006). Statistical mechanics of community detection. Physical
 * Review E 74, 016110. https://doi.org/10.1103/PhysRevE.74.016110
 *
 * \param graph      The input graph.
 * \param membership Numeric vector of integer values which gives the type of each
 *                   vertex, i.e. the cluster to which it belongs.
 *                   It does not have to be consecutive, i.e. empty communities
 *                   are allowed. For better performance, ensure that community
 *                   indices are nonnegative and smaller than the vertex count.
 *                   This can be ensured using \ref igraph_reindex_membership().
 * \param weights    Weight vector or \c NULL if no weights are specified.
 * \param resolution The resolution parameter \c γ. Must not be negative.
 *                   Set it to 1 to use the classical definition of modularity.
 * \param directed   Whether to use the directed or undirected version of modularity.
 *                   Ignored for undirected graphs.
 * \param modularity Pointer to a real number, the result will be
 *                   stored here.
 * \return Error code.
 *
 * \sa \ref igraph_modularity_matrix()
 *
 * Time complexity: O(|V|+|E|), the number of vertices plus the number
 * of edges, assuming that community indices are nonnegative and smaller
 * than the vertex count. Otherwise, O(|V| log |V| + |E|).
 */
igraph_error_t igraph_modularity(const igraph_t *graph,
                      const igraph_vector_int_t *membership,
                      const igraph_vector_t *weights,
                      const igraph_real_t resolution,
                      const igraph_bool_t directed,
                      igraph_real_t *modularity) {

    const igraph_int_t vcount = igraph_vcount(graph);
    const igraph_int_t ecount = igraph_ecount(graph);
    const igraph_vector_int_t *p_membership;
    igraph_vector_int_t i_membership;
    igraph_bool_t using_i_membership = false;
    igraph_vector_t k_out, k_in;
    igraph_int_t min_cluster_id, max_cluster_id, no_of_partitions;
    igraph_real_t e; /* count/fraction of edges/weights within partitions */
    igraph_real_t m; /* edge count / weight sum */
    igraph_int_t c1, c2;
    /* Only consider the graph as directed if it actually is directed */
    igraph_bool_t use_directed = directed && igraph_is_directed(graph);
    igraph_real_t directed_multiplier = (use_directed ? 1 : 2);

    if (igraph_vector_int_size(membership) != vcount) {
        IGRAPH_ERROR("Membership vector size differs from number of vertices.",
                     IGRAPH_EINVAL);
    }
    if (resolution < 0.0) {
      IGRAPH_ERROR("The resolution parameter must not be negative.", IGRAPH_EINVAL);
    }

    if (ecount == 0) {
        /* Special case: the modularity of graphs with no edges is not
         * well-defined */
        *modularity = IGRAPH_NAN;
        return IGRAPH_SUCCESS;
    }

    /* At this point, the 'membership' vector does not have length zero,
       thus it is safe to call igraph_vector_int_minmax(). */

    /* If community indices are outside of the standard range, automatically
     * reindex them. */
    igraph_vector_int_minmax(membership, &min_cluster_id, &max_cluster_id);
    if (min_cluster_id < 0 || max_cluster_id >= vcount) {
        IGRAPH_CHECK(igraph_vector_int_init_copy(&i_membership, membership));
        IGRAPH_FINALLY(igraph_vector_int_destroy, &i_membership);
        IGRAPH_CHECK(igraph_i_reindex_membership_large(&i_membership, NULL, &no_of_partitions));
        p_membership = &i_membership;
        using_i_membership = true;
    } else {
        no_of_partitions = max_cluster_id + 1;
        p_membership = membership;
    }

    IGRAPH_VECTOR_INIT_FINALLY(&k_out, no_of_partitions);
    IGRAPH_VECTOR_INIT_FINALLY(&k_in, no_of_partitions);

    e = 0.0;
    if (weights) {
        if (igraph_vector_size(weights) != ecount)
            IGRAPH_ERROR("Weight vector size differs from number of edges.",
                         IGRAPH_EINVAL);
        m = 0.0;
        for (igraph_int_t i = 0; i < ecount; i++) {
            igraph_real_t w = VECTOR(*weights)[i];
            if (w < 0) {
                IGRAPH_ERROR("Negative weight in weight vector.", IGRAPH_EINVAL);
            }
            c1 = VECTOR(*p_membership)[ IGRAPH_FROM(graph, i) ];
            c2 = VECTOR(*p_membership)[ IGRAPH_TO(graph, i) ];
            if (c1 == c2) {
                e += directed_multiplier * w;
            }
            VECTOR(k_out)[c1] += w;
            VECTOR(k_in)[c2]  += w;
            m += w;
        }
    } else {
        m = ecount;
        for (igraph_int_t i = 0; i < ecount; i++) {
            c1 = VECTOR(*p_membership)[ IGRAPH_FROM(graph, i) ];
            c2 = VECTOR(*p_membership)[ IGRAPH_TO(graph, i) ];
            if (c1 == c2) {
                e += directed_multiplier;
            }
            VECTOR(k_out)[c1] += 1;
            VECTOR(k_in)[c2]  += 1;
        }
    }

    if (!use_directed) {
        /* Graph is undirected, simply add vectors */
        igraph_vector_add(&k_out, &k_in);
        igraph_vector_update(&k_in, &k_out);
    }

    /* Divide all vectors by total weight. */
    igraph_vector_scale(&k_out, 1.0/( directed_multiplier * m ) );
    igraph_vector_scale(&k_in, 1.0/( directed_multiplier * m ) );
    e /= directed_multiplier * m;

    if (m > 0) {
        *modularity = e;
        for (igraph_int_t i = 0; i < no_of_partitions; i++) {
            *modularity -= resolution * VECTOR(k_out)[i] * VECTOR(k_in)[i];
        }
    } else {
        *modularity = IGRAPH_NAN;
    }

    igraph_vector_destroy(&k_out);
    igraph_vector_destroy(&k_in);
    IGRAPH_FINALLY_CLEAN(2);

    if (using_i_membership) {
        igraph_vector_int_destroy(&i_membership);
        IGRAPH_FINALLY_CLEAN(1);
    }

    return IGRAPH_SUCCESS;
}

static igraph_error_t igraph_i_modularity_matrix_get_adjacency(
           const igraph_t *graph, igraph_matrix_t *res,
           const igraph_vector_t *weights, igraph_bool_t directed) {

    /* Specifically used to handle weights and/or ignore direction */
    igraph_eit_t edgeit;
    igraph_int_t no_of_nodes = igraph_vcount(graph);
    igraph_int_t from, to;

    IGRAPH_CHECK(igraph_matrix_resize(res, no_of_nodes, no_of_nodes));
    igraph_matrix_null(res);
    IGRAPH_CHECK(igraph_eit_create(graph, igraph_ess_all(IGRAPH_EDGEORDER_ID), &edgeit));
    IGRAPH_FINALLY(igraph_eit_destroy, &edgeit);

    if (weights) {
        for (; !IGRAPH_EIT_END(edgeit); IGRAPH_EIT_NEXT(edgeit)) {
            igraph_int_t edge = IGRAPH_EIT_GET(edgeit);
            from = IGRAPH_FROM(graph, edge);
            to = IGRAPH_TO(graph, edge);
            MATRIX(*res, from, to) += VECTOR(*weights)[edge];
            if (!directed) {
                MATRIX(*res, to, from) += VECTOR(*weights)[edge];
            }
        }
    } else {
        for (; !IGRAPH_EIT_END(edgeit); IGRAPH_EIT_NEXT(edgeit)) {
            igraph_int_t edge = IGRAPH_EIT_GET(edgeit);
            from = IGRAPH_FROM(graph, edge);
            to = IGRAPH_TO(graph, edge);
            MATRIX(*res, from, to) += 1;
            if (!directed) {
                MATRIX(*res, to, from) += 1;
            }
        }
    }

    igraph_eit_destroy(&edgeit);
    IGRAPH_FINALLY_CLEAN(1);
    return IGRAPH_SUCCESS;
}

/**
 * \function igraph_modularity_matrix
 * \brief Calculates the modularity matrix.
 *
 * This function returns the modularity matrix, which is defined as
 *
 * </para><para>
 * <code>B_ij = A_ij - γ k_i k_j / (2m)</code>
 *
 * </para><para>
 * for undirected graphs, where \c A_ij is the adjacency matrix, \c γ is the
 * resolution parameter, \c k_i is the degree of vertex \c i, and \c m is the
 * number of edges in the graph. When there are no edges, or the weights add up
 * to zero, the result is undefined.
 *
 * </para><para>
 * For directed graphs the modularity matrix is changed to
 *
 * </para><para>
 * <code>B_ij = A_ij - γ k^out_i k^in_j / m</code>
 *
 * </para><para>
 * where <code>k^out_i</code> is the out-degree of node \c i and <code>k^in_j</code> is the
 * in-degree of node \c j.
 *
 * </para><para>
 * Note that self-loops in undirected graphs are multiplied by 2 in this
 * implementation. If weights are specified, the weighted counterparts of the adjacency
 * matrix and degrees are used.
 *
 * \param graph      The input graph.
 * \param weights    Edge weights, pointer to a vector. If this is a null pointer
 *                   then every edge is assumed to have a weight of 1.
 * \param resolution The resolution parameter \c γ. Must not be negative.
 *                   Default is 1. Lower values favor fewer, larger communities;
 *                   higher values favor more, smaller communities.
 * \param modmat     Pointer to an initialized matrix in which the modularity
 *                   matrix is stored.
 * \param directed   For directed graphs: if the edges should be treated as
 *                   undirected. For undirected graphs this is ignored.
 * \return Error code.
 *
 * \sa \ref igraph_modularity()
 */
igraph_error_t igraph_modularity_matrix(const igraph_t *graph,
                             const igraph_vector_t *weights,
                             const igraph_real_t resolution,
                             igraph_matrix_t *modmat,
                             igraph_bool_t directed) {

    const igraph_int_t no_of_nodes = igraph_vcount(graph);
    const igraph_int_t no_of_edges = igraph_ecount(graph);
    const igraph_real_t sw = weights ? igraph_vector_sum(weights) : no_of_edges;
    igraph_vector_t deg, in_deg, out_deg;
    igraph_real_t scaling_factor;

    if (weights && igraph_vector_size(weights) != no_of_edges) {
        IGRAPH_ERROR("Invalid weight vector length.", IGRAPH_EINVAL);
    }

    if (resolution < 0.0) {
        IGRAPH_ERROR("The resolution parameter must not be negative.", IGRAPH_EINVAL);
    }

    if (!igraph_is_directed(graph)) {
        directed = false;
    }
    IGRAPH_CHECK(igraph_i_modularity_matrix_get_adjacency(graph, modmat, weights, directed));

    /* Performance notes:
     *  - Iterating in column-major order makes a large difference.
     *  - Applying the scaling_factor to in_deg (or out_deg) first to reduce the
     *    number of multiplications does not make an appreciable performance
     *    difference. However, doing this in the undirected case causes the result
     *    matrix to sometimes not be strictly symmetric due to the non-associativity
     *    of floating point multiplication.
     */

    if (directed) {
        IGRAPH_VECTOR_INIT_FINALLY(&in_deg, no_of_nodes);
        IGRAPH_VECTOR_INIT_FINALLY(&out_deg, no_of_nodes);

        IGRAPH_CHECK(igraph_strength(graph, &in_deg, igraph_vss_all(), IGRAPH_IN,
                                     IGRAPH_LOOPS, weights));
        IGRAPH_CHECK(igraph_strength(graph, &out_deg, igraph_vss_all(), IGRAPH_OUT,
                                     IGRAPH_LOOPS, weights));

        scaling_factor = resolution / sw;

        for (igraph_int_t j = 0; j < no_of_nodes; j++) {
            for (igraph_int_t i = 0; i < no_of_nodes; i++) {
                MATRIX(*modmat, i, j) -= VECTOR(out_deg)[i] * VECTOR(in_deg)[j] * scaling_factor;
            }
        }
        igraph_vector_destroy(&in_deg);
        igraph_vector_destroy(&out_deg);
        IGRAPH_FINALLY_CLEAN(2);
    } else {
        IGRAPH_VECTOR_INIT_FINALLY(&deg, no_of_nodes);

        IGRAPH_CHECK(igraph_strength(graph, &deg, igraph_vss_all(), IGRAPH_ALL,
                                     IGRAPH_LOOPS, weights));

        scaling_factor = resolution / 2.0 / sw;

        for (igraph_int_t j = 0; j < no_of_nodes; j++) {
            for (igraph_int_t i = 0; i < no_of_nodes; i++) {
                MATRIX(*modmat, i, j) -= VECTOR(deg)[i] * VECTOR(deg)[j] * scaling_factor;
            }
        }
        igraph_vector_destroy(&deg);
        IGRAPH_FINALLY_CLEAN(1);
    }

    return IGRAPH_SUCCESS;
}