igraph/src/clustertool.cpp

/* -*- mode: C -*-  */
/*
   IGraph library.
   Copyright (C) 2006-2012  Gabor Csardi <csardi.gabor@gmail.com>
   334 Harvard street, Cambridge, MA 02139 USA

   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

*/

/* The original version of this file was written by Joerg Reichardt
   The original copyright notice follows here */

/***************************************************************************
                          main.cpp  -  description
                             -------------------
    begin                : Tue Jul 13 11:26:47 CEST 2004
    copyright            : (C) 2004 by
    email                :
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   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.                                   *
 *                                                                         *
 ***************************************************************************/

#ifdef HAVE_CONFIG_H
    #include <config.h>
#endif

#include "NetDataTypes.h"
#include "NetRoutines.h"
#include "pottsmodel_2.h"

#include "igraph_community.h"
#include "igraph_error.h"
#include "igraph_random.h"
#include "igraph_math.h"
#include "igraph_interface.h"
#include "igraph_components.h"
#include "igraph_interrupt_internal.h"
#include "igraph_handle_exceptions.h"

static int igraph_i_community_spinglass_orig(
        const igraph_t *graph,
        const igraph_vector_t *weights,
        igraph_real_t *modularity,
        igraph_real_t *temperature,
        igraph_vector_t *membership,
        igraph_vector_t *csize,
        igraph_integer_t spins,
        igraph_bool_t parupdate,
        igraph_real_t starttemp,
        igraph_real_t stoptemp,
        igraph_real_t coolfact,
        igraph_spincomm_update_t update_rule,
        igraph_real_t gamma);

static int igraph_i_community_spinglass_negative(
        const igraph_t *graph,
        const igraph_vector_t *weights,
        igraph_real_t *modularity,
        igraph_real_t *temperature,
        igraph_vector_t *membership,
        igraph_vector_t *csize,
        igraph_integer_t spins,
        igraph_bool_t parupdate,
        igraph_real_t starttemp,
        igraph_real_t stoptemp,
        igraph_real_t coolfact,
        igraph_spincomm_update_t update_rule,
        igraph_real_t gamma,
        /* igraph_matrix_t *adhesion, */
        /* igraph_matrix_t *normalised_adhesion, */
        /* igraph_real_t *polarization, */
        igraph_real_t gamma_minus);

/**
 * \function igraph_community_spinglass
 * \brief Community detection based on statistical mechanics
 *
 * This function implements the community structure detection
 * algorithm proposed by Joerg Reichardt and Stefan Bornholdt.
 * The algorithm is described in their paper: Statistical Mechanics of
 * Community Detection, http://arxiv.org/abs/cond-mat/0603718 .
 *
 * </para><para> From version 0.6 igraph also supports an extension to
 * the algorithm that allows negative edge weights. This is described
 * in  V.A. Traag and Jeroen Bruggeman: Community detection in networks
 * with positive and negative links, http://arxiv.org/abs/0811.2329 .
 * \param graph The input graph, it may be directed but the direction
 *     of the edge is not used in the algorithm.
 * \param weights The vector giving the edge weights, it may be \c NULL,
 *     in which case all edges are weighted equally. Edge weights
 *     should be positive, altough this is not tested.
 * \param modularity Pointer to a real number, if not \c NULL then the
 *     modularity score of the solution will be stored here. This is the
 *     gereralized modularity that simplifies to the one defined in
 *     M. E. J. Newman and M. Girvan, Phys. Rev. E 69, 026113 (2004),
 *     if the gamma parameter is one.
 * \param temperature Pointer to a real number, if not \c NULL then
 *     the temperature at the end of the algorithm will be stored
 *     here.
 * \param membership Pointer to an initialized vector or \c NULL. If
 *     not \c NULL then the result of the clustering will be stored
 *     here, for each vertex the number of its cluster is given, the
 *     first cluster is numbered zero. The vector will be resized as
 *     needed.
 * \param csize Pointer to an initialized vector or \c NULL. If not \c
 *     NULL then the sizes of the clusters will stored here in cluster
 *     number order. The vector will be resized as needed.
 * \param spins Integer giving the number of spins, ie. the maximum
 *     number of clusters. Usually it is not a program to give a high
 *     number here, the default was 25 in the original code. Even if
 *     the number of spins is high the number of clusters in the
 *     result might small.
 * \param parupdate A logical constant, whether to update all spins in
 *     parallel. The default for this argument was \c FALSE (ie. 0) in
 *     the original code. It is not implemented in the \c
 *     IGRAPH_SPINCOMM_INP_NEG implementation.
 * \param starttemp Real number, the temperature at the start. The
 *     value of this argument was 1.0 in the original code.
 * \param stoptemp Real number, the algorithm stops at this
 *     temperature. The default was 0.01 in the original code.
 * \param coolfact Real number, the coolinf factor for the simulated
 *     annealing. The default was 0.99 in the original code.
 * \param update_rule The type of the update rule. Possible values: \c
 *     IGRAPH_SPINCOMM_UPDATE_SIMPLE and \c
 *     IGRAPH_SPINCOMM_UPDATE_CONFIG. Basically this parameter defined
 *     the null model based on which the actual clustering is done. If
 *     this is \c IGRAPH_SPINCOMM_UPDATE_SIMPLE then the random graph
 *     (ie. G(n,p)), if it is \c IGRAPH_SPINCOMM_UPDATE then the
 *     configuration model is used. The configuration means that the
 *     baseline for the clustering is a random graph with the same
 *     degree distribution as the input graph.
 * \param gamma Real number. The gamma parameter of the
 *     algorithm. This defined the weight of the missing and existing
 *     links in the quality function for the clustering. The default
 *     value in the original code was 1.0, which is equal weight to
 *     missing and existing edges. Smaller values make the existing
 *     links contibute more to the energy function which is minimized
 *     in the algorithm. Bigger values make the missing links more
 *     important. (If my understanding is correct.)
 * \param implementation Constant, chooses between the two
 *     implementations of the spin-glass algorithm that are included
 *     in igraph. \c IGRAPH_SPINCOMM_IMP_ORIG selects the original
 *     implementation, this is faster, \c IGRAPH_SPINCOMM_INP_NEG selects
 *     a new implementation by Vincent Traag that allows negative edge
 *     weights.
 * \param gamma_minus Real number. Parameter for the \c
 *     IGRAPH_SPINCOMM_IMP_NEG implementation. This
 *     specifies the balance between the importance of present and
 *     non-present negative weighted edges in a community. Smaller values of
 *     \p gamma_minus lead to communities with lesser
 *     negative intra-connectivity.
 *     If this argument is set to zero, the algorithm reduces to a graph
 *     coloring algorithm, using the number of spins as the number of
 *     colors.
 * \return Error code.
 *
 * \sa igraph_community_spinglass_single() for calculating the community
 * of a single vertex.
 *
 * Time complexity: TODO.
 *
 * \example examples/simple/spinglass.c
 */

int igraph_community_spinglass(const igraph_t *graph,
                               const igraph_vector_t *weights,
                               igraph_real_t *modularity,
                               igraph_real_t *temperature,
                               igraph_vector_t *membership,
                               igraph_vector_t *csize,
                               igraph_integer_t spins,
                               igraph_bool_t parupdate,
                               igraph_real_t starttemp,
                               igraph_real_t stoptemp,
                               igraph_real_t coolfact,
                               igraph_spincomm_update_t update_rule,
                               igraph_real_t gamma,
                               /* the rest is for the NegSpin implementation */
                               igraph_spinglass_implementation_t implementation,
                               /*                 igraph_matrix_t *adhesion, */
                               /*                 igraph_matrix_t *normalised_adhesion, */
                               /*                 igraph_real_t *polarization, */
                               igraph_real_t gamma_minus) {

    IGRAPH_HANDLE_EXCEPTIONS(
        switch (implementation) {
        case IGRAPH_SPINCOMM_IMP_ORIG:
            return igraph_i_community_spinglass_orig(graph, weights, modularity,
                    temperature, membership, csize,
                    spins, parupdate, starttemp,
                    stoptemp, coolfact, update_rule,
                    gamma);
            break;
        case IGRAPH_SPINCOMM_IMP_NEG:
            return igraph_i_community_spinglass_negative(graph, weights, modularity,
                    temperature, membership, csize,
                    spins, parupdate, starttemp,
                    stoptemp, coolfact,
                    update_rule, gamma,
                    /*                       adhesion, normalised_adhesion, */
                    /*                       polarization, */
                    gamma_minus);
            break;
        default:
            IGRAPH_ERROR("Unknown `implementation' in spinglass community finding",
                         IGRAPH_EINVAL);
        }
    );
    return 0;
}

static int igraph_i_community_spinglass_orig(
        const igraph_t *graph,
        const igraph_vector_t *weights,
        igraph_real_t *modularity,
        igraph_real_t *temperature,
        igraph_vector_t *membership,
        igraph_vector_t *csize,
        igraph_integer_t spins,
        igraph_bool_t parupdate,
        igraph_real_t starttemp,
        igraph_real_t stoptemp,
        igraph_real_t coolfact,
        igraph_spincomm_update_t update_rule,
        igraph_real_t gamma) {

    unsigned long changes, runs;
    igraph_bool_t use_weights = 0;
    bool zeroT;
    double kT, acc, prob;
    ClusterList<NNode*> *cl_cur;
    network *net;
    PottsModel *pm;

    /* Check arguments */

    if (spins < 2 || spins > 500) {
        IGRAPH_ERROR("Invalid number of spins", IGRAPH_EINVAL);
    }
    if (update_rule != IGRAPH_SPINCOMM_UPDATE_SIMPLE &&
        update_rule != IGRAPH_SPINCOMM_UPDATE_CONFIG) {
        IGRAPH_ERROR("Invalid update rule", IGRAPH_EINVAL);
    }
    if (weights) {
        if (igraph_vector_size(weights) != igraph_ecount(graph)) {
            IGRAPH_ERROR("Invalid weight vector length", IGRAPH_EINVAL);
        }
        use_weights = 1;
    }
    if (coolfact < 0 || coolfact >= 1.0) {
        IGRAPH_ERROR("Invalid cooling factor", IGRAPH_EINVAL);
    }
    if (gamma < 0.0) {
        IGRAPH_ERROR("Invalid gamma value", IGRAPH_EINVAL);
    }
    if (starttemp / stoptemp < 1.0) {
        IGRAPH_ERROR("starttemp should be larger in absolute value than stoptemp",
                     IGRAPH_EINVAL);
    }

    /* Check whether we have a single component */
    igraph_bool_t conn;
    IGRAPH_CHECK(igraph_is_connected(graph, &conn, IGRAPH_WEAK));
    if (!conn) {
        IGRAPH_ERROR("Cannot work with unconnected graph", IGRAPH_EINVAL);
    }

    net = new network;
    net->node_list   = new DL_Indexed_List<NNode*>();
    net->link_list   = new DL_Indexed_List<NLink*>();
    net->cluster_list = new DL_Indexed_List<ClusterList<NNode*>*>();

    /* Transform the igraph_t */
    IGRAPH_CHECK(igraph_i_read_network(graph, weights,
                                       net, use_weights, 0));

    prob = 2.0 * net->sum_weights / double(net->node_list->Size())
           / double(net->node_list->Size() - 1);

    pm = new PottsModel(net, (unsigned int)spins, update_rule);

    /* initialize the random number generator */
    RNG_BEGIN();

    if ((stoptemp == 0.0) && (starttemp == 0.0)) {
        zeroT = true;
    } else {
        zeroT = false;
    }
    if (!zeroT) {
        kT = pm->FindStartTemp(gamma, prob, starttemp);
    } else {
        kT = stoptemp;
    }
    /* assign random initial configuration */
    pm->assign_initial_conf(-1);
    runs = 0;
    changes = 1;

    while (changes > 0 && (kT / stoptemp > 1.0 || (zeroT && runs < 150))) {

        IGRAPH_ALLOW_INTERRUPTION(); /* This is not clean.... */

        runs++;
        if (!zeroT) {
            kT *= coolfact;
            if (parupdate) {
                changes = pm->HeatBathParallelLookup(gamma, prob, kT, 50);
            } else {
                acc = pm->HeatBathLookup(gamma, prob, kT, 50);
                if (acc < (1.0 - 1.0 / double(spins)) * 0.01) {
                    changes = 0;
                } else {
                    changes = 1;
                }
            }
        } else {
            if (parupdate) {
                changes = pm->HeatBathParallelLookupZeroTemp(gamma, prob, 50);
            } else {
                acc = pm->HeatBathLookupZeroTemp(gamma, prob, 50);
                /* less than 1 percent acceptance ratio */
                if (acc < (1.0 - 1.0 / double(spins)) * 0.01) {
                    changes = 0;
                } else {
                    changes = 1;
                }
            }
        }
    } /* while loop */

    pm->WriteClusters(modularity, temperature, csize, membership, kT, gamma);

    while (net->link_list->Size()) {
        delete net->link_list->Pop();
    }
    while (net->node_list->Size()) {
        delete net->node_list->Pop();
    }
    while (net->cluster_list->Size()) {
        cl_cur = net->cluster_list->Pop();
        while (cl_cur->Size()) {
            cl_cur->Pop();
        }
        delete cl_cur;
    }
    delete net->link_list;
    delete net->node_list;
    delete net->cluster_list;

    RNG_END();

    delete net;
    delete pm;

    return 0;
}

/**
 * \function igraph_community_spinglass_single
 * \brief Community of a single node based on statistical mechanics
 *
 * This function implements the community structure detection
 * algorithm proposed by Joerg Reichardt and Stefan Bornholdt. It is
 * described in their paper: Statistical Mechanics of
 * Community Detection, http://arxiv.org/abs/cond-mat/0603718 .
 *
 * </para><para>
 * This function calculates the community of a single vertex without
 * calculating all the communities in the graph.
 *
 * \param graph The input graph, it may be directed but the direction
 *    of the edges is not used in the algorithm.
 * \param weights Pointer to a vector with the weights of the edges.
 *    Alternatively \c NULL can be supplied to have the same weight
 *    for every edge.
 * \param vertex The vertex id of the vertex of which ths community is
 *    calculated.
 * \param community Pointer to an initialized vector, the result, the
 *    ids of the vertices in the community of the input vertex will be
 *    stored here. The vector will be resized as needed.
 * \param cohesion Pointer to a real variable, if not \c NULL the
 *     cohesion index of the community will be stored here.
 * \param adhesion Pointer to a real variable, if not \c NULL the
 *     adhesion index of the community will be stored here.
 * \param inner_links Pointer to an integer, if not \c NULL the
 *     number of edges within the community is stored here.
 * \param outer_links Pointer to an integer, if not \c NULL the
 *     number of edges between the community and the rest of the graph
 *     will be stored here.
 * \param spins The number of spins to use, this can be higher than
 *    the actual number of clusters in the network, in which case some
 *    clusters will contain zero vertices.
 * \param update_rule The type of the update rule. Possible values: \c
 *     IGRAPH_SPINCOMM_UPDATE_SIMPLE and \c
 *     IGRAPH_SPINCOMM_UPDATE_CONFIG. Basically this parameter defined
 *     the null model based on which the actual clustering is done. If
 *     this is \c IGRAPH_SPINCOMM_UPDATE_SIMPLE then the random graph
 *     (ie. G(n,p)), if it is \c IGRAPH_SPINCOMM_UPDATE then the
 *     configuration model is used. The configuration means that the
 *     baseline for the clustering is a random graph with the same
 *     degree distribution as the input graph.
 * \param gamma Real number. The gamma parameter of the
 *     algorithm. This defined the weight of the missing and existing
 *     links in the quality function for the clustering. The default
 *     value in the original code was 1.0, which is equal weight to
 *     missing and existing edges. Smaller values make the existing
 *     links contibute more to the energy function which is minimized
 *     in the algorithm. Bigger values make the missing links more
 *     important. (If my understanding is correct.)
 * \return Error code.
 *
 * \sa igraph_community_spinglass() for the traditional version of the
 * algorithm.
 *
 * Time complexity: TODO.
 */

int igraph_community_spinglass_single(const igraph_t *graph,
                                      const igraph_vector_t *weights,
                                      igraph_integer_t vertex,
                                      igraph_vector_t *community,
                                      igraph_real_t *cohesion,
                                      igraph_real_t *adhesion,
                                      igraph_integer_t *inner_links,
                                      igraph_integer_t *outer_links,
                                      igraph_integer_t spins,
                                      igraph_spincomm_update_t update_rule,
                                      igraph_real_t gamma) {
    IGRAPH_HANDLE_EXCEPTIONS(
        igraph_bool_t use_weights = 0;
        double prob;
        ClusterList<NNode*> *cl_cur;
        network *net;
        PottsModel *pm;
        char startnode[255];

        /* Check arguments */

        if (spins < 2 || spins > 500) {
            IGRAPH_ERROR("Invalid number of spins", IGRAPH_EINVAL);
        }
        if (update_rule != IGRAPH_SPINCOMM_UPDATE_SIMPLE &&
            update_rule != IGRAPH_SPINCOMM_UPDATE_CONFIG) {
            IGRAPH_ERROR("Invalid update rule", IGRAPH_EINVAL);
        }
        if (weights) {
            if (igraph_vector_size(weights) != igraph_ecount(graph)) {
                IGRAPH_ERROR("Invalid weight vector length", IGRAPH_EINVAL);
            }
            use_weights = 1;
        }
        if (gamma < 0.0) {
            IGRAPH_ERROR("Invalid gamme value", IGRAPH_EINVAL);
        }
        if (vertex < 0 || vertex > igraph_vcount(graph)) {
            IGRAPH_ERROR("Invalid vertex id", IGRAPH_EINVAL);
        }

        /* Check whether we have a single component */
        igraph_bool_t conn;
        IGRAPH_CHECK(igraph_is_connected(graph, &conn, IGRAPH_WEAK));
        if (!conn) {
            IGRAPH_ERROR("Cannot work with unconnected graph", IGRAPH_EINVAL);
        }

        net = new network;
        net->node_list   = new DL_Indexed_List<NNode*>();
        net->link_list   = new DL_Indexed_List<NLink*>();
        net->cluster_list = new DL_Indexed_List<ClusterList<NNode*>*>();

        /* Transform the igraph_t */
        IGRAPH_CHECK(igraph_i_read_network(graph, weights,
                                           net, use_weights, 0));

        prob = 2.0 * net->sum_weights / double(net->node_list->Size())
               / double(net->node_list->Size() - 1);

        pm = new PottsModel(net, (unsigned int)spins, update_rule);

        /* initialize the random number generator */
        RNG_BEGIN();

        /* to be exected, if we want to find the community around a particular node*/
        /* the initial conf is needed, because otherwise,
           the degree of the nodes is not in the weight property, stupid!!! */
        pm->assign_initial_conf(-1);
        snprintf(startnode, 255, "%li", (long int)vertex + 1);
        pm->FindCommunityFromStart(gamma, prob, startnode, community,
                                   cohesion, adhesion, inner_links, outer_links);

        while (net->link_list->Size()) {
            delete net->link_list->Pop();
        }
        while (net->node_list->Size()) {
            delete net->node_list->Pop();
        }
        while (net->cluster_list->Size()) {
            cl_cur = net->cluster_list->Pop();
            while (cl_cur->Size()) {
                cl_cur->Pop();
            }
            delete cl_cur;
        }
        delete net->link_list;
        delete net->node_list;
        delete net->cluster_list;

        RNG_END();

        delete net;
        delete pm;
    );

    return 0;
}

static int igraph_i_community_spinglass_negative(
        const igraph_t *graph,
        const igraph_vector_t *weights,
        igraph_real_t *modularity,
        igraph_real_t *temperature,
        igraph_vector_t *membership,
        igraph_vector_t *csize,
        igraph_integer_t spins,
        igraph_bool_t parupdate,
        igraph_real_t starttemp,
        igraph_real_t stoptemp,
        igraph_real_t coolfact,
        igraph_spincomm_update_t update_rule,
        igraph_real_t gamma,
        /* igraph_matrix_t *adhesion, */
        /* igraph_matrix_t *normalised_adhesion, */
        /* igraph_real_t *polarization, */
        igraph_real_t gamma_minus) {

    unsigned long changes, runs;
    igraph_bool_t use_weights = 0;
    bool zeroT;
    double kT, acc;
    ClusterList<NNode*> *cl_cur;
    network *net;
    PottsModelN *pm;
    igraph_real_t d_n;
    igraph_real_t d_p;

    /* Check arguments */

    if (parupdate) {
        IGRAPH_ERROR("Parallel spin update not implemented with "
                     "negative gamma", IGRAPH_UNIMPLEMENTED);
    }

    if (spins < 2 || spins > 500) {
        IGRAPH_ERROR("Invalid number of spins", IGRAPH_EINVAL);
    }
    if (update_rule != IGRAPH_SPINCOMM_UPDATE_SIMPLE &&
        update_rule != IGRAPH_SPINCOMM_UPDATE_CONFIG) {
        IGRAPH_ERROR("Invalid update rule", IGRAPH_EINVAL);
    }
    if (weights) {
        if (igraph_vector_size(weights) != igraph_ecount(graph)) {
            IGRAPH_ERROR("Invalid weight vector length", IGRAPH_EINVAL);
        }
        use_weights = 1;
    }
    if (coolfact < 0 || coolfact >= 1.0) {
        IGRAPH_ERROR("Invalid cooling factor", IGRAPH_EINVAL);
    }
    if (gamma < 0.0) {
        IGRAPH_ERROR("Invalid gamma value", IGRAPH_EINVAL);
    }
    if (starttemp / stoptemp < 1.0) {
        IGRAPH_ERROR("starttemp should be larger in absolute value than stoptemp",
                     IGRAPH_EINVAL);
    }

    /* Check whether we have a single component */
    igraph_bool_t conn;
    IGRAPH_CHECK(igraph_is_connected(graph, &conn, IGRAPH_WEAK));
    if (!conn) {
        IGRAPH_ERROR("Cannot work with unconnected graph", IGRAPH_EINVAL);
    }

    if (weights) {
        igraph_vector_minmax(weights, &d_n, &d_p);
    } else {
        d_n = d_p = 1;
    }

    if (d_n > 0) {
        d_n = 0;
    }
    if (d_p < 0) {
        d_p = 0;
    }
    d_n = -d_n;

    net = new network;
    net->node_list   = new DL_Indexed_List<NNode*>();
    net->link_list   = new DL_Indexed_List<NLink*>();
    net->cluster_list = new DL_Indexed_List<ClusterList<NNode*>*>();

    /* Transform the igraph_t */
    IGRAPH_CHECK(igraph_i_read_network(graph, weights,
                                       net, use_weights, 0));

    bool directed = igraph_is_directed(graph);

    pm = new PottsModelN(net, (unsigned int)spins, directed);

    /* initialize the random number generator */
    RNG_BEGIN();

    if ((stoptemp == 0.0) && (starttemp == 0.0)) {
        zeroT = true;
    } else {
        zeroT = false;
    }

    //Begin at a high enough temperature
    kT = pm->FindStartTemp(gamma, gamma_minus, starttemp);

    /* assign random initial configuration */
    pm->assign_initial_conf(true);

    runs = 0;
    changes = 1;
    acc = 0;
    while (changes > 0 && (kT / stoptemp > 1.0 || (zeroT && runs < 150))) {

        IGRAPH_ALLOW_INTERRUPTION(); /* This is not clean.... */

        runs++;
        kT = kT * coolfact;
        acc = pm->HeatBathLookup(gamma, gamma_minus, kT, 50);
        if (acc < (1.0 - 1.0 / double(spins)) * 0.001) {
            changes = 0;
        } else {
            changes = 1;
        }

    } /* while loop */

    /* These are needed, otherwise 'modularity' is not calculated */
    igraph_matrix_t adhesion, normalized_adhesion;
    igraph_real_t polarization;
    IGRAPH_MATRIX_INIT_FINALLY(&adhesion, 0, 0);
    IGRAPH_MATRIX_INIT_FINALLY(&normalized_adhesion, 0, 0);
    pm->WriteClusters(modularity, temperature, csize, membership,
                      &adhesion, &normalized_adhesion, &polarization,
                      kT, d_p, d_n, gamma, gamma_minus);
    igraph_matrix_destroy(&normalized_adhesion);
    igraph_matrix_destroy(&adhesion);
    IGRAPH_FINALLY_CLEAN(2);

    while (net->link_list->Size()) {
        delete net->link_list->Pop();
    }
    while (net->node_list->Size()) {
        delete net->node_list->Pop();
    }
    while (net->cluster_list->Size()) {
        cl_cur = net->cluster_list->Pop();
        while (cl_cur->Size()) {
            cl_cur->Pop();
        }
        delete cl_cur;
    }

    RNG_END();

    return 0;
}