/* -*- mode: C -*- */
/*
IGraph library.
Copyright (C) 2005 Gabor Csardi <csardi@rmki.kfki.hu>
MTA RMKI, Konkoly-Thege Miklos st. 29-33, Budapest 1121, Hungary
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.h"
#include "memory.h"
#include "random.h"
#include <string.h>
#include <limits.h>
/**
* \section about_structural
*
* <para>These functions usually calculate some structural property
* of a graph, like its diameter, the degree of the nodes, etc.</para>
*/
/**
* \ingroup structural
* \function igraph_diameter
* \brief Calculates the diameter of a graph (longest geodesic).
*
* \param graph The graph object.
* \param pres Pointer to an integer, if not \c NULL then it will contain
* the diameter (the actual distance).
* \param pfrom Pointer to an integer, if not \c NULL it will be set to the
* source vertex of the diameter path.
* \param pto Pointer to an integer, if not \c NULL it will be set to the
* target vertex of the diameter path.
* \param path Pointer to an initialized vector. If not \c NULL the actual
* longest geodesic path will be stored here. The vector will be
* resized as needed.
* \param directed Boolean, whether to consider directed
* paths. Ignored for undirected graphs.
* \param unconn What to do if the graph is not connected. If
* \c TRUE the longest geodesic within a component
* will be returned, otherwise the number of vertices is
* returned. (The rationale behind the latter is that this is
* always longer than the longest possible diameter in a
* graph.)
* \return Error code:
* \c IGRAPH_ENOMEM, not enough memory for
* temporary data.
*
* Time complexity: O(|V||E|), the
* number of vertices times the number of edges.
*/
int igraph_diameter(const igraph_t *graph, igraph_integer_t *pres,
igraph_integer_t *pfrom, igraph_integer_t *pto,
igraph_vector_t *path,
igraph_bool_t directed, igraph_bool_t unconn) {
long int no_of_nodes=igraph_vcount(graph);
long int i, j, n;
long int *already_added;
long int nodes_reached;
long int from=0, to=0;
long int res=0;
igraph_dqueue_t q=IGRAPH_DQUEUE_NULL;
igraph_vector_t *neis;
igraph_integer_t dirmode;
igraph_i_adjlist_t allneis;
if (directed) { dirmode=IGRAPH_OUT; } else { dirmode=IGRAPH_ALL; }
already_added=Calloc(no_of_nodes, long int);
if (already_added==0) {
IGRAPH_ERROR("diameter failed", IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(igraph_free, already_added);
IGRAPH_DQUEUE_INIT_FINALLY(&q, 100);
igraph_i_adjlist_init(graph, &allneis, dirmode);
IGRAPH_FINALLY(igraph_i_adjlist_destroy, &allneis);
for (i=0; i<no_of_nodes; i++) {
nodes_reached=1;
IGRAPH_CHECK(igraph_dqueue_push(&q, i));
IGRAPH_CHECK(igraph_dqueue_push(&q, 0));
already_added[i]=i+1;
igraph_progress("Diameter: ", 100.0*i/no_of_nodes, NULL);
IGRAPH_ALLOW_INTERRUPTION();
while (!igraph_dqueue_empty(&q)) {
long int actnode=igraph_dqueue_pop(&q);
long int actdist=igraph_dqueue_pop(&q);
if (actdist>res) {
res=actdist;
from=i;
to=actnode;
}
neis=igraph_i_adjlist_get(&allneis, actnode);
n=igraph_vector_size(neis);
for (j=0; j<n; j++) {
long int neighbor=VECTOR(*neis)[j];
if (already_added[neighbor] == i+1) { continue; }
already_added[neighbor]=i+1;
nodes_reached++;
IGRAPH_CHECK(igraph_dqueue_push(&q, neighbor));
IGRAPH_CHECK(igraph_dqueue_push(&q, actdist+1));
}
} /* while !igraph_dqueue_empty */
/* not connected, return largest possible */
if (nodes_reached != no_of_nodes && !unconn) {
res=no_of_nodes;
from=-1;
to=-1;
break;
}
} /* for i<no_of_nodes */
igraph_progress("Diameter: ", 100.0, NULL);
/* return the requested info */
if (pres != 0) {
*pres=res;
}
if (pfrom != 0) {
*pfrom=from;
}
if (pto != 0) {
*pto=to;
}
if (path != 0) {
if (res==no_of_nodes) {
igraph_vector_clear(path);
} else {
igraph_vector_ptr_t tmpptr;
igraph_vector_ptr_init(&tmpptr, 1);
IGRAPH_FINALLY(igraph_vector_ptr_destroy, &tmpptr);
VECTOR(tmpptr)[0]=path;
IGRAPH_CHECK(igraph_get_shortest_paths(graph, &tmpptr, from,
igraph_vss_1(to), dirmode));
igraph_vector_ptr_destroy(&tmpptr);
IGRAPH_FINALLY_CLEAN(1);
}
}
/* clean */
Free(already_added);
igraph_dqueue_destroy(&q);
igraph_i_adjlist_destroy(&allneis);
IGRAPH_FINALLY_CLEAN(3);
return 0;
}
/**
* \ingroup structural
* \function igraph_average_path_length
* \brief Calculates the average geodesic length in a graph.
*
* \param graph The graph object.
* \param res Pointer to a real number, this will contain the result.
* \param directed Boolean, whether to consider directed
* paths. Ignored for undirected graphs.
* \param unconn What to do if the graph is not connected. If
* \c TRUE the average of thr geodesics
* within the components
* will be returned, otherwise the number of vertices is
* used for the length of non-existing geodesics. (The rationale
* behind this is that this is always longer than the longest
* possible geodesic in a graph.)
* \return Error code:
* \c IGRAPH_ENOMEM, not enough memory for
* data structures
*
* Time complexity: O(|V||E|), the
* number of vertices times the number of edges.
*/
int igraph_average_path_length(const igraph_t *graph, igraph_real_t *res,
igraph_bool_t directed, igraph_bool_t unconn) {
long int no_of_nodes=igraph_vcount(graph);
long int i, j, n;
long int *already_added;
long int nodes_reached=0;
igraph_real_t normfact=0.0;
igraph_dqueue_t q=IGRAPH_DQUEUE_NULL;
igraph_vector_t *neis;
igraph_integer_t dirmode;
igraph_i_adjlist_t allneis;
*res=0;
if (directed) { dirmode=IGRAPH_OUT; } else { dirmode=IGRAPH_ALL; }
already_added=Calloc(no_of_nodes, long int);
if (already_added==0) {
IGRAPH_ERROR("average path length failed", IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(free, already_added); /* TODO: hack */
IGRAPH_DQUEUE_INIT_FINALLY(&q, 100);
igraph_i_adjlist_init(graph, &allneis, dirmode);
IGRAPH_FINALLY(igraph_i_adjlist_destroy, &allneis);
for (i=0; i<no_of_nodes; i++) {
nodes_reached=0;
IGRAPH_CHECK(igraph_dqueue_push(&q, i));
IGRAPH_CHECK(igraph_dqueue_push(&q, 0));
already_added[i]=i+1;
IGRAPH_ALLOW_INTERRUPTION();
while (!igraph_dqueue_empty(&q)) {
long int actnode=igraph_dqueue_pop(&q);
long int actdist=igraph_dqueue_pop(&q);
neis=igraph_i_adjlist_get(&allneis, actnode);
n=igraph_vector_size(neis);
for (j=0; j<n; j++) {
long int neighbor=VECTOR(*neis)[j];
if (already_added[neighbor] == i+1) { continue; }
already_added[neighbor]=i+1;
nodes_reached++;
*res += actdist+1;
normfact+=1;
IGRAPH_CHECK(igraph_dqueue_push(&q, neighbor));
IGRAPH_CHECK(igraph_dqueue_push(&q, actdist+1));
}
} /* while !igraph_dqueue_empty */
/* not connected, return largest possible */
if (!unconn) {
*res += (no_of_nodes * (no_of_nodes-1-nodes_reached));
normfact += no_of_nodes-1-nodes_reached;
}
} /* for i<no_of_nodes */
*res /= normfact;
/* clean */
Free(already_added);
igraph_dqueue_destroy(&q);
igraph_i_adjlist_destroy(&allneis);
IGRAPH_FINALLY_CLEAN(3);
return 0;
}
/**
* \ingroup structural
* \function igraph_minimum_spanning_tree_unweighted
* \brief Calculates one minimum spanning tree of an unweighted graph.
*
* </para><para>
* If the graph has more minimum spanning trees (this is always the
* case, except if it is a forest) this implementation returns only
* the same one.
*
* </para><para>
* Directed graphs are considered as undirected for this computation.
*
* </para><para>
* If the graph is not connected then its minimum spanning forest is
* returned. This is the set of the minimum spanning trees of each
* component.
* \param graph The graph object.
* \param mst The minimum spanning tree, another graph object. Do
* \em not initialize this object before passing it to
* this function, but be sure to call \ref igraph_destroy() on it if
* you don't need it any more.
* \return Error code:
* \c IGRAPH_ENOMEM, not enough memory for
* temporary data.
*
* Time complexity: O(|V|+|E|),
* |V| is the
* number of vertices, |E| the number
* of edges in the graph.
*
* \sa \ref igraph_minimum_spanning_tree_prim() for weighted graphs.
*/
int igraph_minimum_spanning_tree_unweighted(const igraph_t *graph,
igraph_t *mst) {
long int no_of_nodes=igraph_vcount(graph);
long int no_of_edges=igraph_ecount(graph);
char *already_added;
char *added_edges;
igraph_dqueue_t q=IGRAPH_DQUEUE_NULL;
igraph_vector_t edges=IGRAPH_VECTOR_NULL;
igraph_vector_t tmp=IGRAPH_VECTOR_NULL;
long int i, j;
added_edges=Calloc(no_of_edges, char);
if (added_edges==0) {
IGRAPH_ERROR("unweighted spanning tree failed", IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(igraph_free, added_edges);
IGRAPH_VECTOR_INIT_FINALLY(&edges, 0);
already_added=Calloc(no_of_nodes, char);
if (already_added==0) {
IGRAPH_ERROR("unweighted spanning tree failed", IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(igraph_free, already_added);
IGRAPH_VECTOR_INIT_FINALLY(&tmp, 0);
IGRAPH_DQUEUE_INIT_FINALLY(&q, 100);
for (i=0; i<no_of_nodes; i++) {
if (already_added[i]>0) { continue; }
IGRAPH_ALLOW_INTERRUPTION();
already_added[i]=1;
IGRAPH_CHECK(igraph_dqueue_push(&q, i));
while (! igraph_dqueue_empty(&q)) {
long int act_node=igraph_dqueue_pop(&q);
IGRAPH_CHECK(igraph_adjacent(graph, &tmp, act_node, IGRAPH_ALL));
for (j=0; j<igraph_vector_size(&tmp); j++) {
long int edge=VECTOR(tmp)[j];
if (added_edges[edge]==0) {
igraph_integer_t from, to;
igraph_edge(graph, edge, &from, &to);
if (act_node==to) { to=from; }
if (already_added[(long int) to]==0) {
already_added[(long int) to]=1;
added_edges[edge]=1;
IGRAPH_CHECK(igraph_dqueue_push(&q, to));
}
}
}
}
}
igraph_dqueue_destroy(&q);
Free(already_added);
igraph_vector_destroy(&tmp);
IGRAPH_FINALLY_CLEAN(3);
/* summarize the edges to delete */
j=0;
for (i=0; i<no_of_edges; i++) {
if (added_edges[i]==0) {
j++;
}
}
IGRAPH_CHECK(igraph_vector_resize(&edges, j));
j=0;
for (i=0; i<no_of_edges; i++) {
if (added_edges[i]==0) {
VECTOR(edges)[j++]=i;
}
}
IGRAPH_CHECK(igraph_copy(mst, graph));
IGRAPH_FINALLY(igraph_destroy, mst);
IGRAPH_CHECK(igraph_delete_edges(mst, igraph_ess_vector(&edges)));
igraph_vector_destroy(&edges);
Free(added_edges);
IGRAPH_FINALLY_CLEAN(3);
return 0;
}
/**
* \ingroup structural
* \function igraph_minimum_spanning_tree_prim
* \brief Calculates one minimum spanning tree of a weighted graph.
*
* </para><para>
* This function uses Prim's method for carrying out the computation,
* see Prim, R.C.: Shortest connection networks and some
* generalizations, Bell System Technical
* Journal, Vol. 36,
* 1957, 1389--1401.
*
* </para><para>
* If the graph has more than one minimum spanning tree, the current
* implementation returns always the same one.
*
* </para><para>
* Directed graphs are considered as undirected for this computation.
*
* </para><para>
* If the graph is not connected then its minimum spanning forest is
* returned. This is the set of the minimum spanning trees of each
* component.
*
* \param graph The graph object.
* \param mst The result of the computation, a graph object containing
* the minimum spanning tree of the graph.
* Do \em not initialize this object before passing it to
* this function, but be sure to call \ref igraph_destroy() on it if
* you don't need it any more.
* \param weights A vector containing the weights of the the edges.
* in the same order as the simple edge iterator visits them.
* \return Error code:
* \c IGRAPH_ENOMEM, not enough memory.
* \c IGRAPH_EINVAL, length of weight vector does not
* match number of edges.
*
* Time complexity: O(|V|+|E|),
* |V| is the number of vertices,
* |E| the number of edges in the
* graph.
*
* \sa \ref igraph_minimum_spanning_tree_unweighted() for unweighted graphs.
*/
int igraph_minimum_spanning_tree_prim(const igraph_t *graph, igraph_t *mst,
const igraph_vector_t *weights) {
long int no_of_nodes=igraph_vcount(graph);
long int no_of_edges=igraph_ecount(graph);
char *already_added;
char *added_edges;
igraph_d_indheap_t heap=IGRAPH_D_INDHEAP_NULL;
igraph_vector_t edges=IGRAPH_VECTOR_NULL;
igraph_integer_t mode=IGRAPH_ALL;
igraph_vector_t adj;
long int i, j;
if (igraph_vector_size(weights) != igraph_ecount(graph)) {
IGRAPH_ERROR("Invalid weights length", IGRAPH_EINVAL);
}
added_edges=Calloc(no_of_edges, char);
if (added_edges==0) {
IGRAPH_ERROR("prim spanning tree failed", IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(igraph_destroy, added_edges);
IGRAPH_VECTOR_INIT_FINALLY(&edges, 0);
already_added=Calloc(no_of_nodes, char);
if (already_added == 0) {
IGRAPH_ERROR("prim spanning tree failed", IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(igraph_free, already_added);
IGRAPH_CHECK(igraph_d_indheap_init(&heap, 0));
IGRAPH_FINALLY(igraph_d_indheap_destroy, &heap);
IGRAPH_VECTOR_INIT_FINALLY(&adj, 0);
for (i=0; i<no_of_nodes; i++) {
if (already_added[i]>0) { continue; }
IGRAPH_ALLOW_INTERRUPTION();
already_added[i]=1;
/* add all edges of the first vertex */
igraph_adjacent(graph, &adj, i, mode);
for (j=0; j<igraph_vector_size(&adj); j++) {
long int edgeno=VECTOR(adj)[j];
igraph_integer_t edgefrom, edgeto;
long int neighbor;
igraph_edge(graph, edgeno, &edgefrom, &edgeto);
neighbor= edgefrom != i ? edgefrom : edgeto;
if (already_added[neighbor] == 0) {
IGRAPH_CHECK(igraph_d_indheap_push(&heap, -VECTOR(*weights)[edgeno], i,
edgeno));
}
}
while(! igraph_d_indheap_empty(&heap)) {
/* Get minimal edge */
long int from, edge;
igraph_integer_t tmp, to;
igraph_d_indheap_max_index(&heap, &from, &edge);
igraph_edge(graph, edge, &tmp, &to);
/* Erase it */
igraph_d_indheap_delete_max(&heap);
/* Is this edge already included? */
if (added_edges[edge]==0) {
if (from==to) { to=tmp; }
/* Does it point to a visited node? */
if (already_added[(long int)to]==0) {
already_added[(long int)to]=1;
added_edges[edge]=1;
/* add all outgoing edges */
igraph_adjacent(graph, &adj, to, mode);
for (j=0; j<igraph_vector_size(&adj); j++) {
long int edgeno=VECTOR(adj)[j];
igraph_integer_t edgefrom, edgeto;
long int neighbor;
igraph_edge(graph, edgeno, &edgefrom, &edgeto);
neighbor= edgefrom != to ? edgefrom : edgeto;
if (already_added[neighbor] == 0) {
IGRAPH_CHECK(igraph_d_indheap_push(&heap, -VECTOR(*weights)[edgeno], to,
edgeno));
}
}
} /* for */
} /* if !already_added */
} /* while in the same component */
} /* for all nodes */
igraph_d_indheap_destroy(&heap);
Free(already_added);
igraph_vector_destroy(&adj);
IGRAPH_FINALLY_CLEAN(3);
/* Ok, collect the edges to delete */
j=0;
for (i=0; i<no_of_edges; i++) {
if (added_edges[i]==0) {
j++;
}
}
IGRAPH_CHECK(igraph_vector_resize(&edges, j));
j=0;
for (i=0; i<no_of_edges; i++) {
if (added_edges[i]==0) {
VECTOR(edges)[j++]=i;
}
}
IGRAPH_CHECK(igraph_copy(mst, graph));
IGRAPH_FINALLY(igraph_destroy, mst);
IGRAPH_CHECK(igraph_delete_edges(mst, igraph_ess_vector(&edges)));
igraph_vector_destroy(&edges);
Free(added_edges);
IGRAPH_FINALLY_CLEAN(3);
return 0;
}
/**
* \ingroup structural
* \function igraph_closeness
* \brief Closeness centrality calculations for some vertices.
*
* </para><para>
* The closeness centrality of a vertex measures how easily other
* vertices can be reached from it (or the other way: how easily it
* can be reached from the other vertices). It is defined as the
* number of the number of vertices minus one divided by the sum of the
* lengths of all geodesics from/to the given vertex.
*
* </para><para>
* If the graph is not connected, and there is no path between two
* vertices, the number of vertices is used instead the length of the
* geodesic. This is always longer than the longest possible geodesic.
*
* \param graph The graph object.
* \param res The result of the computation, a vector containing the
* closeness centrality scores for the given vertices.
* \param vids Vector giving the vertices for which the closeness
* centrality scores will be computed.
* \param mode The type of shortest paths to be used for the
* calculation in directed graphs. Possible values:
* \clist
* \cli IGRAPH_OUT
* the lengths of the outgoing paths are calculated.
* \cli IGRAPH_IN
* the lengths of the incoming paths are calculated.
* \cli IGRAPH_ALL
* the directed graph is considered as an
* undirected one for the computation.
* \endclist
* \return Error code:
* \clist
* \cli IGRAPH_ENOMEM
* not enough memory for temporary data.
* \cli IGRAPH_EINVVID
* invalid vertex id passed.
* \cli IGRAPH_EINVMODE
* invalid mode argument.
* \endclist
*
* Time complexity: O(n|E|),
* n is the number
* of vertices for which the calculation is done and
* |E| is the number
* of edges in the graph.
*
* \sa Other centrality types: \ref igraph_degree(), \ref igraph_betweenness().
*/
int igraph_closeness(const igraph_t *graph, igraph_vector_t *res,
igraph_vs_t vids,
igraph_neimode_t mode) {
long int no_of_nodes=igraph_vcount(graph);
igraph_vector_t already_counted;
long int i, j;
long int nodes_reached;
igraph_dqueue_t q;
long int nodes_to_calc;
igraph_vector_t tmp;
igraph_vit_t vit;
IGRAPH_CHECK(igraph_vit_create(graph, vids, &vit));
IGRAPH_FINALLY(igraph_vit_destroy, &vit);
nodes_to_calc=IGRAPH_VIT_SIZE(vit);
if (mode != IGRAPH_OUT && mode != IGRAPH_IN &&
mode != IGRAPH_ALL) {
IGRAPH_ERROR("calculating closeness", IGRAPH_EINVMODE);
}
IGRAPH_VECTOR_INIT_FINALLY(&already_counted, no_of_nodes);
IGRAPH_VECTOR_INIT_FINALLY(&tmp, 0);
IGRAPH_DQUEUE_INIT_FINALLY(&q, 100);
IGRAPH_CHECK(igraph_vector_resize(res, nodes_to_calc));
igraph_vector_null(res);
for (IGRAPH_VIT_RESET(vit), i=0;
!IGRAPH_VIT_END(vit);
IGRAPH_VIT_NEXT(vit), i++) {
IGRAPH_CHECK(igraph_dqueue_push(&q, IGRAPH_VIT_GET(vit)));
IGRAPH_CHECK(igraph_dqueue_push(&q, 0));
nodes_reached=1;
VECTOR(already_counted)[(long int)IGRAPH_VIT_GET(vit)]=i+1;
IGRAPH_ALLOW_INTERRUPTION();
while (!igraph_dqueue_empty(&q)) {
long int act=igraph_dqueue_pop(&q);
long int actdist=igraph_dqueue_pop(&q);
VECTOR(*res)[i] += actdist;
IGRAPH_CHECK(igraph_neighbors(graph, &tmp, act, mode));
for (j=0; j<igraph_vector_size(&tmp); j++) {
long int neighbor=VECTOR(tmp)[j];
if (VECTOR(already_counted)[neighbor] == i+1) { continue; }
VECTOR(already_counted)[neighbor] = i+1;
nodes_reached++;
IGRAPH_CHECK(igraph_dqueue_push(&q, neighbor));
IGRAPH_CHECK(igraph_dqueue_push(&q, actdist+1));
}
}
VECTOR(*res)[i] += (no_of_nodes * (no_of_nodes-nodes_reached));
VECTOR(*res)[i] = (no_of_nodes-1) / VECTOR(*res)[i];
}
/* Clean */
igraph_dqueue_destroy(&q);
igraph_vector_destroy(&tmp);
igraph_vector_destroy(&already_counted);
igraph_vit_destroy(&vit);
IGRAPH_FINALLY_CLEAN(4);
return 0;
}
/**
* \ingroup structural
* \function igraph_shortest_paths
* \brief The length of the shortest paths between vertices.
*
* \param graph The graph object.
* \param res The result of the calculation, a matrix. It has the same
* number of rows as the length of the \c from
* argument, and its number of columns is the number of
* vertices in the graph. One row of the matrix shows the
* distances from/to a given vertex to all the others in the
* graph, the order is fixed by the vertex ids.
* \param from Vector of the vertex ids for which the path length
* calculations are done.
* \param mode The type of shortest paths to be use for the
* calculation in directed graphs. Possible values:
* \clist
* \cli IGRAPH_OUT
* the lengths of the outgoing paths are calculated.
* \cli IGRAPH_IN
* the lengths of the incoming paths are calculated.
* \cli IGRAPH_ALL
* the directed graph is considered as an undirected one for
* the computation.
* \endclist
* \return Error code:
* \clist
* \cli IGRAPH_ENOMEM
* not enough memory for temporary
* data.
* \cli IGRAPH_EINVVID
* invalid vertex id passed.
* \cli IGRAPH_EINVMODE
* invalid mode argument.
* \endclist
*
* Time complexity: O(n(|V|+|E|)),
* n is the
* number of vertices to calculate, |V| and
* |E| are the number of vertices and
* edges in the graph.
*
* \sa \ref igraph_get_shortest_paths() to get the paths themselves.
*/
int igraph_shortest_paths(const igraph_t *graph, igraph_matrix_t *res,
igraph_vs_t from, igraph_neimode_t mode) {
long int no_of_nodes=igraph_vcount(graph);
long int no_of_from;
long int *already_counted;
igraph_dqueue_t q=IGRAPH_DQUEUE_NULL;
long int i, j;
igraph_vector_t tmp=IGRAPH_VECTOR_NULL;
igraph_vit_t fromvit;
IGRAPH_CHECK(igraph_vit_create(graph, from, &fromvit));
IGRAPH_FINALLY(igraph_vit_destroy, &fromvit);
no_of_from=IGRAPH_VIT_SIZE(fromvit);
if (mode != IGRAPH_OUT && mode != IGRAPH_IN &&
mode != IGRAPH_ALL) {
IGRAPH_ERROR("Invalid mode argument", IGRAPH_EINVMODE);
}
already_counted=Calloc(no_of_nodes, long int);
if (already_counted==0) {
IGRAPH_ERROR("shortest paths failed", IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(free, already_counted);
IGRAPH_VECTOR_INIT_FINALLY(&tmp, 0);
IGRAPH_DQUEUE_INIT_FINALLY(&q, 100);
IGRAPH_CHECK(igraph_matrix_resize(res, no_of_from, no_of_nodes));
igraph_matrix_null(res);
for (IGRAPH_VIT_RESET(fromvit), i=0;
!IGRAPH_VIT_END(fromvit);
IGRAPH_VIT_NEXT(fromvit), i++) {
long int reached=1;
IGRAPH_CHECK(igraph_dqueue_push(&q, IGRAPH_VIT_GET(fromvit)));
IGRAPH_CHECK(igraph_dqueue_push(&q, 0));
already_counted[ (long int) IGRAPH_VIT_GET(fromvit) ] = i+1;
IGRAPH_ALLOW_INTERRUPTION();
while (!igraph_dqueue_empty(&q)) {
long int act=igraph_dqueue_pop(&q);
long int actdist=igraph_dqueue_pop(&q);
MATRIX(*res, i, act)=actdist;
IGRAPH_CHECK(igraph_neighbors(graph, &tmp, act, mode));
for (j=0; j<igraph_vector_size(&tmp); j++) {
long int neighbor=VECTOR(tmp)[j];
if (already_counted[neighbor] == i+1) { continue; }
already_counted[neighbor] = i+1;
reached++;
IGRAPH_CHECK(igraph_dqueue_push(&q, neighbor));
IGRAPH_CHECK(igraph_dqueue_push(&q, actdist+1));
}
}
/* Plus the unreachable nodes */
j=0;
while (reached < no_of_nodes) {
if (MATRIX(*res, i, j) == 0 && j != IGRAPH_VIT_GET(fromvit)) {
MATRIX(*res, i, j)=no_of_nodes;
reached++;
}
j++;
}
}
/* Clean */
igraph_vector_destroy(&tmp);
Free(already_counted);
igraph_dqueue_destroy(&q);
igraph_vit_destroy(&fromvit);
IGRAPH_FINALLY_CLEAN(4);
return 0;
}
/**
* \ingroup structural
* \function igraph_get_shortest_paths
* \brief Calculates the shortest paths from/to one vertex.
*
* </para><para>
* If there is more than one geodesic between two vertices, this
* function gives only one of them.
* \param graph The graph object.
* \param res The result, this is a pointer vector, each element points
* to a vector
* object. These should be initialized before passing them to
* the function, which will properly clear and/or resize them
* and fill the ids of the vertices along the geodesics from/to
* the vertices.
* \param from The id of the vertex from/to which the geodesics are
* calculated.
* \param to Vertex sequence with the ids of the vertices to/from which the
* shortest paths will be calculated. A vertex might be given multiple
* times.
* \param mode The type of shortest paths to be use for the
* calculation in directed graphs. Possible values:
* \clist
* \cli IGRAPH_OUT
* the outgoing paths are calculated.
* \cli IGRAPH_IN
* the incoming paths are calculated.
* \cli IGRAPH_ALL
* the directed graph is considered as an
* undirected one for the computation.
* \endclist
* \return Error code:
* \clist
* \cli IGRAPH_ENOMEM
* not enough memory for temporary data.
* \cli IGRAPH_EINVVID
* \p from is invalid vertex id, or the length of \p to is
* not the same as the length of \p res.
* \cli IGRAPH_EINVMODE
* invalid mode argument.
* \endclist
*
* Time complexity: O(|V|+|E|),
* |V| is the number of vertices,
* |E| the number of edges in the
* graph.
*
* \sa \ref igraph_shortest_paths() if you only need the path length but
* not the paths themselves.
*/
int igraph_get_shortest_paths(const igraph_t *graph, igraph_vector_ptr_t *res,
igraph_integer_t from, igraph_vs_t to,
igraph_neimode_t mode) {
/* TODO: use adjlist_t if to is long (longer than 1?) */
long int no_of_nodes=igraph_vcount(graph);
long int *father;
igraph_dqueue_t q=IGRAPH_DQUEUE_NULL;
long int j;
igraph_vector_t tmp=IGRAPH_VECTOR_NULL;
igraph_vit_t vit;
long int to_reach;
long int reached=0;
if (from<0 || from>=no_of_nodes) {
IGRAPH_ERROR("cannot get shortest paths", IGRAPH_EINVVID);
}
if (mode != IGRAPH_OUT && mode != IGRAPH_IN &&
mode != IGRAPH_ALL) {
IGRAPH_ERROR("Invalid mode argument", IGRAPH_EINVMODE);
}
IGRAPH_CHECK(igraph_vit_create(graph, to, &vit));
IGRAPH_FINALLY(igraph_vit_destroy, &vit);
if (IGRAPH_VIT_SIZE(vit) != igraph_vector_ptr_size(res)) {
IGRAPH_ERROR("Size of the `res' and the `to' should match", IGRAPH_EINVAL);
}
father=Calloc(no_of_nodes, long int);
if (father==0) {
IGRAPH_ERROR("cannot get shortest paths", IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(free, father); /* TODO: hack */
IGRAPH_VECTOR_INIT_FINALLY(&tmp, 0);
IGRAPH_DQUEUE_INIT_FINALLY(&q, 100);
/* Mark the vertices we need to reach */
to_reach=IGRAPH_VIT_SIZE(vit);
for (IGRAPH_VIT_RESET(vit); !IGRAPH_VIT_END(vit); IGRAPH_VIT_NEXT(vit)) {
if (father[ (long int) IGRAPH_VIT_GET(vit) ] == 0) {
father[ (long int) IGRAPH_VIT_GET(vit) ] = -1;
} else {
to_reach--; /* this node was given multiple times */
}
}
IGRAPH_CHECK(igraph_dqueue_push(&q, from+1));
if (father[ (long int) from ] < 0) { reached++; }
father[ (long int)from ] = from+1;
while (!igraph_dqueue_empty(&q) && reached < to_reach) {
long int act=igraph_dqueue_pop(&q);
IGRAPH_CHECK(igraph_neighbors(graph, &tmp, act-1, mode));
for (j=0; j<igraph_vector_size(&tmp); j++) {
long int neighbor=VECTOR(tmp)[j]+1;
if (father[neighbor-1] > 0) {
continue;
} else if (father[neighbor-1] < 0) {
reached++;
}
father[neighbor-1] = act;
IGRAPH_CHECK(igraph_dqueue_push(&q, neighbor));
}
}
if (reached < to_reach) {
IGRAPH_WARNING("Couldn't reach some vertices");
}
for (IGRAPH_VIT_RESET(vit), j=0;
!IGRAPH_VIT_END(vit);
IGRAPH_VIT_NEXT(vit), j++) {
long int node=IGRAPH_VIT_GET(vit);
igraph_vector_t *vec=VECTOR(*res)[j];
igraph_vector_clear(vec);
IGRAPH_ALLOW_INTERRUPTION();
if (father[node]>0) {
long int act=node+1;
long int size=0;
while (father[act-1] != act) {
size++;
act=father[act-1];
}
size++;
IGRAPH_CHECK(igraph_vector_resize(vec, size));
VECTOR(*vec)[--size]=node;
act=node+1;
while (father[act-1] != act) {
VECTOR(*vec)[--size]=father[act-1]-1;
act=father[act-1];
}
}
}
/* Clean */
Free(father);
igraph_dqueue_destroy(&q);
igraph_vector_destroy(&tmp);
igraph_vit_destroy(&vit);
IGRAPH_FINALLY_CLEAN(4);
return 0;
}
void igraph_i_gasp_paths_destroy(igraph_vector_ptr_t *v) {
long int i;
for (i=0; i<igraph_vector_ptr_size(v); i++) {
if (VECTOR(*v)[i] != 0) {
igraph_vector_destroy(VECTOR(*v)[i]);
Free(VECTOR(*v)[i]);
}
}
igraph_vector_ptr_destroy(v);
}
/**
* \function igraph_get_all_shortest_paths
* \brief Finds all shortest paths (geodesics) from a vertex to all
* other vertices
*
* \param graph The graph object.
* \param res Pointer to an initialized pointer vector, the result
* will be stored here in igraph_vector_t objects. Each vector
* object contains the vertices along a shortest path from \p from
* to another vertex. The vectors are ordered according to their
* target vertex: first the shortest paths to vertex 0, then to
* vertex 1, etc. No data is included for unreachable vertices.
* \param nrgeo Pointer to an initialized igraph_vector_t object or
* NULL. If not NULL the number of shortest paths from \p from are
* is stored here for every vertex in the graph.
* \param from The id of the vertex from/to which the geodesics are
* calculated.
* \param mode The type of shortest paths to be use for the
* calculation in directed graphs. Possible values:
* \clist
* \cli IGRAPH_OUT
* the lengths of the outgoing paths are calculated.
* \cli IGRAPH_IN
* the lengths of the incoming paths are calculated.
* \cli IGRAPH_ALL
* the directed graph is considered as an
* undirected one for the computation.
* \endclist
* \return Error code:
* \clist
* \cli IGRAPH_ENOMEM
* not enough memory for temporary data.
* \cli IGRAPH_EINVVID
* \p from is invalid vertex id.
* \cli IGRAPH_EINVMODE
* invalid mode argument.
* \endclist
*
* Added in version 0.2.</para><para>
*
* Time complexity: O(|V|+|E|) for most graphs, O(|V|^2) in the worst
* case.
*/
int igraph_get_all_shortest_paths(const igraph_t *graph,
igraph_vector_ptr_t *res,
igraph_vector_t *nrgeo,
igraph_integer_t from, igraph_vs_t to,
igraph_neimode_t mode) {
long int no_of_nodes=igraph_vcount(graph);
long int *geodist;
igraph_vector_ptr_t paths;
igraph_dqueue_t q;
igraph_vector_t *vptr;
igraph_vector_t neis;
igraph_vector_t ptrlist;
igraph_vector_t ptrhead;
long int n, j, i;
igraph_vit_t vit;
if (from<0 || from>=no_of_nodes) {
IGRAPH_ERROR("cannot get shortest paths", IGRAPH_EINVVID);
}
if (mode != IGRAPH_OUT && mode != IGRAPH_IN &&
mode != IGRAPH_ALL) {
IGRAPH_ERROR("Invalid mode argument", IGRAPH_EINVMODE);
}
IGRAPH_CHECK(igraph_vit_create(graph, to, &vit));
IGRAPH_FINALLY(igraph_vit_destroy, &vit);
IGRAPH_CHECK(igraph_vector_ptr_init(&paths, 0));
IGRAPH_FINALLY(igraph_i_gasp_paths_destroy, &paths);
IGRAPH_VECTOR_INIT_FINALLY(&neis, 0);
IGRAPH_VECTOR_INIT_FINALLY(&ptrlist, 0);
IGRAPH_VECTOR_INIT_FINALLY(&ptrhead, no_of_nodes);
geodist=Calloc(no_of_nodes, long int);
if (geodist==0) {
IGRAPH_ERROR("Cannot calculate shortest paths", IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(igraph_free, geodist);
IGRAPH_CHECK(igraph_dqueue_init(&q, 100));
IGRAPH_FINALLY(igraph_dqueue_destroy, &q);
if (nrgeo) {
IGRAPH_CHECK(igraph_vector_resize(nrgeo, no_of_nodes));
igraph_vector_null(nrgeo);
}
/* from -> from */
vptr=Calloc(1, igraph_vector_t); /* TODO: dirty */
IGRAPH_CHECK(igraph_vector_ptr_push_back(&paths, vptr));
IGRAPH_CHECK(igraph_vector_init(vptr, 1));
VECTOR(*vptr)[0]=from;
geodist[(long int)from]=1;
VECTOR(ptrhead)[(long int)from]=1;
IGRAPH_CHECK(igraph_vector_push_back(&ptrlist, 0));
if (nrgeo) { VECTOR(*nrgeo)[(long int)from]=1; }
/* Init queue */
IGRAPH_CHECK(igraph_dqueue_push(&q, from));
IGRAPH_CHECK(igraph_dqueue_push(&q, 0.0));
while (!igraph_dqueue_empty(&q)) {
long int actnode=igraph_dqueue_pop(&q);
long int actdist=igraph_dqueue_pop(&q);
IGRAPH_ALLOW_INTERRUPTION();
IGRAPH_CHECK(igraph_neighbors(graph, &neis, actnode, mode));
n=igraph_vector_size(&neis);
for (j=0; j<n; j++) {
long int neighbor=VECTOR(neis)[j];
long int fatherptr=VECTOR(ptrhead)[actnode];
if (geodist[neighbor] != 0 &&
geodist[neighbor]-1 < actdist+1) { continue; }
if (nrgeo) { VECTOR(*nrgeo)[neighbor] += VECTOR(*nrgeo)[actnode]; }
if (geodist[neighbor] == 0) {
IGRAPH_CHECK(igraph_dqueue_push(&q, neighbor));
IGRAPH_CHECK(igraph_dqueue_push(&q, actdist+1));
}
geodist[neighbor]=actdist+2;
/* copy all existing paths to the parent */
while (fatherptr != 0) {
vptr=Calloc(1, igraph_vector_t);
IGRAPH_CHECK(igraph_vector_ptr_push_back(&paths, vptr));
IGRAPH_CHECK(igraph_vector_copy(vptr, VECTOR(paths)[fatherptr-1]));
IGRAPH_CHECK(igraph_vector_reserve(vptr, actdist+2));
igraph_vector_push_back(vptr, neighbor);
IGRAPH_CHECK(igraph_vector_push_back(&ptrlist,
VECTOR(ptrhead)[neighbor]));
VECTOR(ptrhead)[neighbor]=igraph_vector_size(&ptrlist);
fatherptr=VECTOR(ptrlist)[fatherptr-1];
}
}
}
igraph_dqueue_destroy(&q);
IGRAPH_FINALLY_CLEAN(1);
/* Copy to the result */
memset(geodist, 0, sizeof(long int)*no_of_nodes);
for (IGRAPH_VIT_RESET(vit); !IGRAPH_VIT_END(vit); IGRAPH_VIT_NEXT(vit)) {
geodist[ (long int) IGRAPH_VIT_GET(vit) ] = 1;
}
n=0;
for (i=0; i<no_of_nodes; i++) {
long int fatherptr=VECTOR(ptrhead)[i];
if (geodist[i] > 0) {
while (fatherptr != 0) {
n++;
fatherptr=VECTOR(ptrlist)[fatherptr-1];
}
}
}
IGRAPH_CHECK(igraph_vector_ptr_resize(res, n));
j=0;
for (i=0; i<no_of_nodes; i++) {
long int fatherptr=VECTOR(ptrhead)[i];
IGRAPH_ALLOW_INTERRUPTION();
if (geodist[i] > 0) {
while (fatherptr != 0) {
VECTOR(*res)[j++]=VECTOR(paths)[fatherptr-1];
fatherptr=VECTOR(ptrlist)[fatherptr-1];
}
} else {
while (fatherptr != 0) {
igraph_vector_destroy(VECTOR(paths)[fatherptr-1]);
Free(VECTOR(paths)[fatherptr-1]);
fatherptr=VECTOR(ptrlist)[fatherptr-1];
}
}
}
Free(geodist);
igraph_vector_destroy(&ptrlist);
igraph_vector_destroy(&ptrhead);
igraph_vector_destroy(&neis);
igraph_vector_ptr_destroy(&paths);
igraph_vit_destroy(&vit);
IGRAPH_FINALLY_CLEAN(6);
return 0;
}
/**
* \ingroup structural
* \function igraph_subcomponent
* \brief The vertices in the same component as a given vertex.
*
* \param graph The graph object.
* \param res The result, vector with the ids of the vertices in the
* same component.
* \param vertex The id of the vertex of which the component is
* searched.
* \param mode Type of the component for directed graphs, possible
* values:
* \clist
* \cli IGRAPH_OUT
* the set of vertices reachable \em from the
* \p vertex,
* \cli IGRAPH_IN
* the set of vertices from which the
* \p vertex is reachable.
* \cli IGRAPH_ALL
* the graph is considered as an
* undirected graph. Note that this is \em not the same
* as the union of the previous two.
* \endclist
* \return Error code:
* \clist
* \cli IGRAPH_ENOMEM
* not enough memory for temporary data.
* \cli IGRAPH_EINVVID
* \p vertex is an invalid vertex id
* \cli IGRAPH_EINVMODE
* invalid mode argument passed.
* \endclist
*
* Time complexity: O(|V|+|E|),
* |V| and
* |E| are the number of vertices and
* edges in the graph.
*
* \sa \ref igraph_subgraph() if you want a graph object consisting only
* a given set of vertices and the edges between them.
*/
int igraph_subcomponent(const igraph_t *graph, igraph_vector_t *res, igraph_real_t vertex,
igraph_neimode_t mode) {
long int no_of_nodes=igraph_vcount(graph);
igraph_dqueue_t q=IGRAPH_DQUEUE_NULL;
char *already_added;
long int i;
igraph_vector_t tmp=IGRAPH_VECTOR_NULL;
if (vertex<0 || vertex>=no_of_nodes) {
IGRAPH_ERROR("subcomponent failed", IGRAPH_EINVVID);
}
if (mode != IGRAPH_OUT && mode != IGRAPH_IN &&
mode != IGRAPH_ALL) {
IGRAPH_ERROR("invalid mode argument", IGRAPH_EINVMODE);
}
already_added=Calloc(no_of_nodes, char);
if (already_added==0) {
IGRAPH_ERROR("subcomponent failed",IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(free, already_added); /* TODO: hack */
IGRAPH_VECTOR_INIT_FINALLY(&tmp, 0);
IGRAPH_DQUEUE_INIT_FINALLY(&q, 100);
IGRAPH_CHECK(igraph_dqueue_push(&q, vertex));
IGRAPH_CHECK(igraph_vector_push_back(res, vertex));
already_added[(long int)vertex]=1;
while (!igraph_dqueue_empty(&q)) {
long int actnode=igraph_dqueue_pop(&q);
IGRAPH_ALLOW_INTERRUPTION();
IGRAPH_CHECK(igraph_neighbors(graph, &tmp, actnode, mode));
for (i=0; i<igraph_vector_size(&tmp); i++) {
long int neighbor=VECTOR(tmp)[i];
if (already_added[neighbor]) { continue; }
already_added[neighbor]=1;
IGRAPH_CHECK(igraph_vector_push_back(res, neighbor));
IGRAPH_CHECK(igraph_dqueue_push(&q, neighbor));
}
}
igraph_dqueue_destroy(&q);
igraph_vector_destroy(&tmp);
Free(already_added);
IGRAPH_FINALLY_CLEAN(3);
return 0;
}
/**
* \ingroup structural
* \function igraph_betweenness
* \brief Betweenness centrality of some vertices.
*
* </para><para>
* The betweenness centrality of a vertex is the number of geodesics
* going through it. If there are more than one geodesic between two
* vertices, the value of these geodesics are weighted by one over the
* number of geodesics.
* \param graph The graph object.
* \param res The result of the computation, a vector containing the
* betweenness scores for the specified vertices.
* \param vids The vertices of which the betweenness centrality scores
* will be calculated.
* \param directed Logical, if true directed paths will be considered
* for directed graphs. It is ignored for undirected graphs.
* \return Error code:
* \c IGRAPH_ENOMEM, not enough memory for
* temporary data.
* \c IGRAPH_EINVVID, invalid vertex id passed in
* \p vids.
*
* Time complexity: O(|V||E|),
* |V| and
* |E| are the number of vertices and
* edges in the graph.
* Note that the time complexity is independent of the number of
* vertices for which the score is calculated.
*
* \sa Other centrality types: \ref igraph_degree(), \ref igraph_closeness().
* See \ref igraph_edge_betweenness() for calculating the betweenness score
* of the edges in a graph.
*/
int igraph_betweenness (const igraph_t *graph, igraph_vector_t *res,
igraph_vs_t vids,
igraph_bool_t directed) {
long int no_of_nodes=igraph_vcount(graph);
igraph_dqueue_t q=IGRAPH_DQUEUE_NULL;
long int *distance;
long int *nrgeo;
double *tmpscore;
igraph_stack_t stack=IGRAPH_STACK_NULL;
long int source;
long int j;
igraph_vector_t tmp=IGRAPH_VECTOR_NULL;
igraph_integer_t modein, modeout;
igraph_vit_t vit;
IGRAPH_CHECK(igraph_vit_create(graph, vids, &vit));
IGRAPH_FINALLY(igraph_vit_destroy, &vit);
if (directed)
{ modeout=IGRAPH_OUT; modein=IGRAPH_IN; }
else
{ modeout=modein=IGRAPH_ALL; }
distance=Calloc(no_of_nodes, long int);
if (distance==0) {
IGRAPH_ERROR("betweenness failed", IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(free, distance); /* TODO: hack */
nrgeo=Calloc(no_of_nodes, long int);
if (nrgeo==0) {
IGRAPH_ERROR("betweenness failed", IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(free, nrgeo); /* TODO: hack */
tmpscore=Calloc(no_of_nodes, double);
if (tmpscore==0) {
IGRAPH_ERROR("betweenness failed", IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(free, tmpscore); /* TODO: hack */
IGRAPH_VECTOR_INIT_FINALLY(&tmp, 0);
IGRAPH_DQUEUE_INIT_FINALLY(&q, 100);
igraph_stack_init(&stack, no_of_nodes);
IGRAPH_FINALLY(igraph_stack_destroy, &stack);
IGRAPH_CHECK(igraph_vector_resize(res, IGRAPH_VIT_SIZE(vit)));
igraph_vector_null(res);
/* here we go */
for (source=0; source<no_of_nodes; source++) {
IGRAPH_ALLOW_INTERRUPTION();
memset(distance, 0, no_of_nodes*sizeof(long int));
memset(nrgeo, 0, no_of_nodes*sizeof(long int));
memset(tmpscore, 0, no_of_nodes*sizeof(double));
igraph_stack_clear(&stack); /* it should be empty anyway... */
IGRAPH_CHECK(igraph_dqueue_push(&q, source));
nrgeo[source]=1;
distance[source]=0;
while (!igraph_dqueue_empty(&q)) {
long int actnode=igraph_dqueue_pop(&q);
IGRAPH_CHECK(igraph_neighbors(graph, &tmp, actnode, modeout));
for (j=0; j<igraph_vector_size(&tmp); j++) {
long int neighbor=VECTOR(tmp)[j];
if (nrgeo[neighbor] != 0) {
/* we've already seen this node, another shortest path? */
if (distance[neighbor]==distance[actnode]+1) {
nrgeo[neighbor]+=nrgeo[actnode];
}
} else {
/* we haven't seen this node yet */
nrgeo[neighbor]+=nrgeo[actnode];
distance[neighbor]=distance[actnode]+1;
IGRAPH_CHECK(igraph_dqueue_push(&q, neighbor));
IGRAPH_CHECK(igraph_stack_push(&stack, neighbor));
}
}
} /* while !igraph_dqueue_empty */
/* Ok, we've the distance of each node and also the number of
shortest paths to them. Now we do an inverse search, starting
with the farthest nodes. */
while (!igraph_stack_empty(&stack)) {
long int actnode=igraph_stack_pop(&stack);
if (distance[actnode]<=1) { continue; } /* skip source node */
/* set the temporary score of the friends */
IGRAPH_CHECK(igraph_neighbors(graph, &tmp, actnode, modein));
for (j=0; j<igraph_vector_size(&tmp); j++) {
long int neighbor=VECTOR(tmp)[j];
if (distance[neighbor]==distance[actnode]-1 &&
nrgeo[neighbor] != 0) {
tmpscore[neighbor] +=
(tmpscore[actnode]+1)*nrgeo[neighbor]/nrgeo[actnode];
}
}
}
/* Ok, we've the scores for this source */
for (IGRAPH_VIT_RESET(vit); !IGRAPH_VIT_END(vit); IGRAPH_VIT_NEXT(vit)) {
long int node=IGRAPH_VIT_GET(vit);
VECTOR(*res)[node] += tmpscore[node];
tmpscore[node] = 0.0; /* in case a node is in vids multiple times */
}
} /* for source < no_of_nodes */
/* divide by 2 for undirected graph */
if (!directed || !igraph_is_directed(graph)) {
for (j=0; j<igraph_vector_size(res); j++) {
VECTOR(*res)[j] /= 2.0;
}
}
/* clean */
Free(distance);
Free(nrgeo);
Free(tmpscore);
igraph_dqueue_destroy(&q);
igraph_stack_destroy(&stack);
igraph_vector_destroy(&tmp);
igraph_vit_destroy(&vit);
IGRAPH_FINALLY_CLEAN(7);
return 0;
}
/**
* \ingroup structural
* \function igraph_edge_betweenness
* \brief Betweenness centrality of the edges.
*
* </para><para>
* The betweenness centrality of an edge is the number of geodesics
* going through it. If there are more than one geodesics between two
* vertices, the value of these geodesics are weighted by one over the
* number of geodesics.
* \param graph The graph object.
* \param result The result of the computation, vector containing the
* betweenness scores for the edges.
* \param directed Logical, if true directed paths will be considered
* for directed graphs. It is ignored for undirected graphs.
* \return Error code:
* \c IGRAPH_ENOMEM, not enough memory for
* temporary data.
*
* Time complexity: O(|V||E|),
* |V| and
* |E| are the number of vertices and
* edges in the graph.
*
* \sa Other centrality types: \ref igraph_degree(), \ref igraph_closeness().
* See \ref igraph_edge_betweenness() for calculating the betweenness score
* of the edges in a graph.
*/
int igraph_edge_betweenness (const igraph_t *graph, igraph_vector_t *result,
igraph_bool_t directed) {
long int no_of_nodes=igraph_vcount(graph);
long int no_of_edges=igraph_ecount(graph);
igraph_dqueue_t q=IGRAPH_DQUEUE_NULL;
long int *distance;
long int *nrgeo;
double *tmpscore;
igraph_stack_t stack=IGRAPH_STACK_NULL;
long int source;
long int j;
igraph_i_adjedgelist_t elist_out, elist_in;
igraph_i_adjedgelist_t *elist_out_p, *elist_in_p;
igraph_vector_t *neip;
long int neino;
long int i;
igraph_integer_t modein, modeout;
directed=directed && igraph_is_directed(graph);
if (directed) {
modeout=IGRAPH_OUT;
modein=IGRAPH_IN;
IGRAPH_CHECK(igraph_i_adjedgelist_init(graph, &elist_out, IGRAPH_OUT));
IGRAPH_FINALLY(igraph_i_adjedgelist_destroy, &elist_out);
IGRAPH_CHECK(igraph_i_adjedgelist_init(graph, &elist_in, IGRAPH_IN));
IGRAPH_FINALLY(igraph_i_adjedgelist_destroy, &elist_in);
elist_out_p=&elist_out;
elist_in_p=&elist_in;
} else {
modeout=modein=IGRAPH_ALL;
IGRAPH_CHECK(igraph_i_adjedgelist_init(graph,&elist_out, IGRAPH_ALL));
IGRAPH_FINALLY(igraph_i_adjedgelist_destroy, &elist_out);
elist_out_p=elist_in_p=&elist_out;
}
distance=Calloc(no_of_nodes, long int);
if (distance==0) {
IGRAPH_ERROR("edge betweenness failed", IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(free, distance); /* TODO: hack */
nrgeo=Calloc(no_of_nodes, long int);
if (nrgeo==0) {
IGRAPH_ERROR("edge betweenness failed", IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(free, nrgeo); /* TODO: hack */
tmpscore=Calloc(no_of_nodes, double);
if (tmpscore==0) {
IGRAPH_ERROR("edge betweenness failed", IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(free, tmpscore); /* TODO: hack */
IGRAPH_DQUEUE_INIT_FINALLY(&q, 100);
IGRAPH_CHECK(igraph_stack_init(&stack, no_of_nodes));
IGRAPH_FINALLY(igraph_stack_destroy, &stack);
IGRAPH_CHECK(igraph_vector_resize(result, no_of_edges));
igraph_vector_null(result);
/* here we go */
for (source=0; source<no_of_nodes; source++) {
IGRAPH_ALLOW_INTERRUPTION();
memset(distance, 0, no_of_nodes*sizeof(long int));
memset(nrgeo, 0, no_of_nodes*sizeof(long int));
memset(tmpscore, 0, no_of_nodes*sizeof(double));
igraph_stack_clear(&stack); /* it should be empty anyway... */
IGRAPH_CHECK(igraph_dqueue_push(&q, source));
nrgeo[source]=1;
distance[source]=0;
while (!igraph_dqueue_empty(&q)) {
long int actnode=igraph_dqueue_pop(&q);
neip=igraph_i_adjedgelist_get(elist_out_p, actnode);
neino=igraph_vector_size(neip);
for (i=0; i<neino; i++) {
igraph_integer_t edge=VECTOR(*neip)[i], from, to;
long int neighbor;
igraph_edge(graph, edge, &from, &to);
neighbor = actnode!=from ? from : to;
if (nrgeo[neighbor] != 0) {
/* we've already seen this node, another shortest path? */
if (distance[neighbor]==distance[actnode]+1) {
nrgeo[neighbor]+=nrgeo[actnode];
}
} else {
/* we haven't seen this node yet */
nrgeo[neighbor]+=nrgeo[actnode];
distance[neighbor]=distance[actnode]+1;
IGRAPH_CHECK(igraph_dqueue_push(&q, neighbor));
IGRAPH_CHECK(igraph_stack_push(&stack, neighbor));
}
}
} /* while !igraph_dqueue_empty */
/* Ok, we've the distance of each node and also the number of
shortest paths to them. Now we do an inverse search, starting
with the farthest nodes. */
while (!igraph_stack_empty(&stack)) {
long int actnode=igraph_stack_pop(&stack);
if (distance[actnode]<1) { continue; } /* skip source node */
/* set the temporary score of the friends */
neip=igraph_i_adjedgelist_get(elist_in_p, actnode);
neino=igraph_vector_size(neip);
for (i=0; i<neino; i++) {
igraph_integer_t from, to;
long int neighbor;
long int edgeno=VECTOR(*neip)[i];
igraph_edge(graph, edgeno, &from, &to);
neighbor= actnode != from ? from : to;
if (distance[neighbor]==distance[actnode]-1 &&
nrgeo[neighbor] != 0) {
tmpscore[neighbor] +=
(tmpscore[actnode]+1)*nrgeo[neighbor]/nrgeo[actnode];
VECTOR(*result)[edgeno] +=
(tmpscore[actnode]+1)*nrgeo[neighbor]/nrgeo[actnode];
}
}
}
/* Ok, we've the scores for this source */
} /* for source <= no_of_nodes */
/* clean and return */
Free(distance);
Free(nrgeo);
Free(tmpscore);
igraph_dqueue_destroy(&q);
igraph_stack_destroy(&stack);
IGRAPH_FINALLY_CLEAN(5);
if (directed) {
igraph_i_adjedgelist_destroy(&elist_out);
igraph_i_adjedgelist_destroy(&elist_in);
IGRAPH_FINALLY_CLEAN(2);
} else {
igraph_i_adjedgelist_destroy(&elist_out);
IGRAPH_FINALLY_CLEAN(1);
}
/* divide by 2 for undirected graph */
if (!directed || !igraph_is_directed(graph)) {
for (j=0; j<igraph_vector_size(result); j++) {
VECTOR(*result)[j] /= 2.0;
}
}
return 0;
}
/**
* \ingroup structural
* \function igraph_pagerank
* \brief Calculates the Google PageRank for the specified vertices.
*
* </para><para>
* Please note that the PageRank of a given vertex depends on the PageRank
* of all other vertices, so even if you want to calculate the PageRank for
* only some of the vertices, all of them must be calculated. Requesting
* the PageRank for only some of the vertices does not result in any
* performance increase at all.
* </para>
* <para>
* Since the calculation is an iterative
* process, the algorithm is stopped after a given count of iterations
* or if the PageRank value differences between iterations are less than
* a predefined value.
* </para>
*
* <para>
* For the explanation of the PageRank algorithm, see the following
* webpage:
* http://www-db.stanford.edu/~backrub/google.html, or the
* following reference:
* </para>
*
* <para>
* Sergey Brin and Larry Page: The Anatomy of a Large-Scale Hypertextual
* Web Search Engine. Proceedings of the 7th World-Wide Web Conference,
* Brisbane, Australia, April 1998.
* </para>
* <para>
* \param graph The graph object.
* \param res The result vector containing the PageRank values for the
* given nodes.
* \param vids Vector with the vertex ids
* \param directed Logical, if true directed paths will be considered
* for directed graphs. It is ignored for undirected graphs.
* \param niter The maximum number of iterations to perform
* \param eps The algorithm will consider the calculation as complete
* if the difference of PageRank values between iterations change
* less than this value for every node
* \param damping The damping factor ("d" in the original paper)
* \return Error code:
* \c IGRAPH_ENOMEM, not enough memory for
* temporary data.
* \c IGRAPH_EINVVID, invalid vertex id in
* \p vids.
*
* Time complexity: TODO.
*/
int igraph_pagerank(const igraph_t *graph, igraph_vector_t *res,
igraph_vs_t vids, igraph_bool_t directed, igraph_integer_t niter,
igraph_real_t eps, igraph_real_t damping) {
long int no_of_nodes=igraph_vcount(graph);
long int i, j, n, nodes_to_calc;
igraph_real_t *prvec, *prvec_new, *prvec_aux, *prvec_scaled;
igraph_vector_t *neis, outdegree;
igraph_integer_t dirmode;
igraph_i_adjlist_t allneis;
igraph_real_t maxdiff=eps;
igraph_vit_t vit;
if (niter<=0) IGRAPH_ERROR("Invalid iteration count", IGRAPH_EINVAL);
if (eps<=0) IGRAPH_ERROR("Invalid epsilon value", IGRAPH_EINVAL);
if (damping<=0 || damping>=1) IGRAPH_ERROR("Invalid damping factor", IGRAPH_EINVAL);
IGRAPH_CHECK(igraph_vit_create(graph, vids, &vit));
IGRAPH_FINALLY(igraph_vit_destroy, &vit);
nodes_to_calc=IGRAPH_VIT_SIZE(vit);
IGRAPH_CHECK(igraph_vector_resize(res, nodes_to_calc));
igraph_vector_null(res);
IGRAPH_VECTOR_INIT_FINALLY(&outdegree, no_of_nodes);
prvec=Calloc(no_of_nodes, igraph_real_t);
if (prvec==0) {
IGRAPH_ERROR("pagerank failed", IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(igraph_free, prvec);
prvec_new=Calloc(no_of_nodes, igraph_real_t);
if (prvec_new==0) {
IGRAPH_ERROR("pagerank failed", IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(igraph_free, prvec_new);
prvec_scaled=Calloc(no_of_nodes, igraph_real_t);
if (prvec_scaled==0) {
IGRAPH_ERROR("pagerank failed", IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(igraph_free, prvec_scaled);
if (directed) { dirmode=IGRAPH_IN; } else { dirmode=IGRAPH_ALL; }
igraph_i_adjlist_init(graph, &allneis, dirmode);
IGRAPH_FINALLY(igraph_i_adjlist_destroy, &allneis);
/* Calculate outdegrees for every node */
igraph_degree(graph, &outdegree, igraph_vss_all(),
directed?IGRAPH_OUT:IGRAPH_ALL, 0);
/* Initialize PageRank values */
for (i=0; i<no_of_nodes; i++) {
prvec[i]=1-damping;
/* The next line is necessary to avoid division by zero in the
* calculation of prvec_scaled. This won't cause any problem,
* since if a node doesn't have any outgoing links, its
* prvec_scaled value won't be used anywhere */
if (VECTOR(outdegree)[i]==0) VECTOR(outdegree)[i]=1;
}
/* We will always calculate the new PageRank values into prvec_new
* based on the existing values from prvec. To avoid unnecessary
* copying from prvec_new to prvec at the end of every iteration,
* the pointers are swapped after every iteration */
while (niter>0 && maxdiff >= eps) {
niter--;
maxdiff=0;
/* Calculate the quotient of the actual PageRank value and the
* outdegree for every node */
for (i=0; i<no_of_nodes; i++) {
prvec_scaled[i]=prvec[i]/VECTOR(outdegree)[i];
}
/* Calculate new PageRank values based on the old ones */
for (i=0; i<no_of_nodes; i++) {
IGRAPH_ALLOW_INTERRUPTION();
prvec_new[i]=0;
neis=igraph_i_adjlist_get(&allneis, i);
n=igraph_vector_size(neis);
for (j=0; j<n; j++) {
long int neighbor=VECTOR(*neis)[j];
prvec_new[i]+=prvec_scaled[neighbor];
}
prvec_new[i]*=damping;
prvec_new[i]+=(1-damping);
if (prvec_new[i]-prvec[i]>maxdiff)
maxdiff=prvec_new[i]-prvec[i];
else if (prvec[i]-prvec_new[i]>maxdiff)
maxdiff=prvec[i]-prvec_new[i];
}
/* Swap the vectors */
prvec_aux=prvec_new;
prvec_new=prvec;
prvec=prvec_aux;
}
/* Copy results from prvec to res */
for (IGRAPH_VIT_RESET(vit), i=0;
!IGRAPH_VIT_END(vit);
IGRAPH_VIT_NEXT(vit), i++) {
long int vid=IGRAPH_VIT_GET(vit);
VECTOR(*res)[i]=prvec[vid];
}
igraph_i_adjlist_destroy(&allneis);
igraph_vit_destroy(&vit);
igraph_vector_destroy(&outdegree);
Free(prvec);
Free(prvec_new);
Free(prvec_scaled);
IGRAPH_FINALLY_CLEAN(6);
return 0;
}
/**
* \ingroup structural
* \function igraph_rewire
* \brief Randomly rewires a graph while preserving the degree distribution.
*
* </para><para>
* This function generates a new graph based on the original one by randomly
* rewiring edges while preserving the original graph's degree distribution.
* Please note that the rewiring is done "in place", so no new graph will
* be generated. If you would like to keep the original graph intact, use
* \ref igraph_copy() before.
*
* \param graph The graph object to be rewired.
* \param n Number of rewiring trials to perform.
* \param mode The rewiring algorithm to be used. It can be one of the following:
* \c IGRAPH_REWIRING_SIMPLE: simple rewiring algorithm which
* chooses two arbitrary edges in each step (namely (a,b) and (c,d))
* and substitutes them with (a,d) and (c,b) if they don't exist.
* Time complexity: TODO.
* \return Error code:
* \clist
* \cli IGRAPH_EINVMODE
* Invalid rewiring mode.
* \cli IGRAPH_EINVAL
* Graph unsuitable for rewiring (e.g. it has
* less than 4 nodes in case of \c IGRAPH_REWIRING_SIMPLE)
* \cli IGRAPH_ENOMEM
* Not enough memory for temporary data.
* \endclist
*
* Time complexity: TODO.
*/
int igraph_rewire(igraph_t *graph, igraph_integer_t n, igraph_rewiring_t mode) {
long int no_of_nodes=igraph_vcount(graph);
long int i, a, b, c, d;
igraph_i_adjlist_t allneis;
igraph_vector_t *neis[2], edgevec;
igraph_es_t es;
if (mode == IGRAPH_REWIRING_SIMPLE && no_of_nodes<4)
IGRAPH_ERROR("graph unsuitable for rewiring", IGRAPH_EINVAL);
RNG_BEGIN();
igraph_i_adjlist_init(graph, &allneis, IGRAPH_OUT);
IGRAPH_FINALLY(igraph_i_adjlist_destroy, &allneis);
igraph_vector_init(&edgevec, 4);
IGRAPH_FINALLY(igraph_vector_destroy, &edgevec);
while (n>0) {
IGRAPH_ALLOW_INTERRUPTION();
switch (mode) {
case IGRAPH_REWIRING_SIMPLE:
a=RNG_INTEGER(0, no_of_nodes-1);
do { c=RNG_INTEGER(0, no_of_nodes-1); } while (c==a);
/* We don't want the algorithm to get stuck in an infinite loop when
* it can't choose two edges satisfying the conditions. Instead of
* this, we choose two arbitrary edges and if they have endpoints
* in common, we just decrease the number of trials left and continue
* (so unsuccessful rewirings still count as a trial)
*/
neis[0]=igraph_i_adjlist_get(&allneis, a);
i=igraph_vector_size(neis[0]);
if (i==0) b=c;
else b=VECTOR(*neis[0])[RNG_INTEGER(0, i-1)];
neis[1]=igraph_i_adjlist_get(&allneis, c);
i=igraph_vector_size(neis[1]);
if (i==0) d=a;
else d=VECTOR(*neis[1])[RNG_INTEGER(0, i-1)];
/* Okay, we have two edges. Can they be rewired?
* neis[0] mustn't contain d and neis[1] mustn't contain b
*/
if (!igraph_vector_search(neis[0], 0, d, NULL) &&
!igraph_vector_search(neis[1], 0, b, NULL) &&
b!=c && a!=d && a!=b && c!=d) {
/* TODO: maybe it would be more efficient to implement an
* igraph_replace_edges function?
*/
VECTOR(edgevec)[0]=a; VECTOR(edgevec)[1]=b;
VECTOR(edgevec)[2]=c; VECTOR(edgevec)[3]=d;
/*printf("Deleting: %d -> %d, %d -> %d\n", a, b, c, d);*/
IGRAPH_CHECK(igraph_es_pairs(&es, &edgevec, IGRAPH_DIRECTED));
IGRAPH_FINALLY(igraph_es_destroy, &es);
IGRAPH_CHECK(igraph_delete_edges(graph, es));
igraph_es_destroy(&es);
IGRAPH_FINALLY_CLEAN(1);
VECTOR(edgevec)[0]=a; VECTOR(edgevec)[1]=d;
VECTOR(edgevec)[2]=c; VECTOR(edgevec)[3]=b;
/*printf("Adding: %d -> %d, %d -> %d\n", a, d, c, b);*/
igraph_add_edges(graph, &edgevec, 0);
/* We have to adjust the adjacency list view as well.
It's a luck that we have the pointers in neis[0] and neis[1] */
for (i=igraph_vector_size(neis[0])-1; i>=0; i--)
if (VECTOR(*neis[0])[i]==b) { VECTOR(*neis[0])[i]=d; break; }
for (i=igraph_vector_size(neis[1])-1; i>=0; i--)
if (VECTOR(*neis[1])[i]==d) { VECTOR(*neis[1])[i]=b; break; }
}
break;
default:
RNG_END();
IGRAPH_ERROR("unknown rewiring mode", IGRAPH_EINVMODE);
}
n--;
}
igraph_i_adjlist_destroy(&allneis);
igraph_vector_destroy(&edgevec);
IGRAPH_FINALLY_CLEAN(2);
RNG_END();
return 0;
}
/**
* \ingroup structural
* \function igraph_subgraph
* \brief Creates a subgraph with the specified vertices.
*
* </para><para>
* This function collects the specified vertices and all edges between
* them to a new graph.
* As the vertex ids in a graph always start with one, this function
* very likely needs to reassign ids to the vertices.
* \param graph The graph object.
* \param res The subgraph, another graph object will be stored here,
* do \em not initialize this object before calling this
* function, and call \ref igraph_destroy() on it if you don't need
* it any more.
* \param vids Vector with the vertex ids to put in the subgraph.
* \return Error code:
* \c IGRAPH_ENOMEM, not enough memory for
* temporary data.
* \c IGRAPH_EINVVID, invalid vertex id in
* \p vids.
*
* Time complexity: O(|V|+|E|),
* |V| and
* |E| are the number of vertices and
* edges in the original graph.
*
* \sa \ref igraph_delete_vertices() to delete the specified set of
* vertices from a graph, the opposite of this function.
*/
int igraph_subgraph(const igraph_t *graph, igraph_t *res,
igraph_vs_t vids) {
long int no_of_nodes=igraph_vcount(graph);
igraph_vector_t delete=IGRAPH_VECTOR_NULL;
char *remain;
long int i;
igraph_vit_t vit;
IGRAPH_CHECK(igraph_vit_create(graph, vids, &vit));
IGRAPH_FINALLY(igraph_vit_destroy, &vit);
IGRAPH_VECTOR_INIT_FINALLY(&delete, 0);
remain=Calloc(no_of_nodes, char);
if (remain==0) {
IGRAPH_ERROR("subgraph failed", IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(free, remain); /* TODO: hack */
IGRAPH_CHECK(igraph_vector_reserve(&delete, no_of_nodes-IGRAPH_VIT_SIZE(vit)));
for (IGRAPH_VIT_RESET(vit); !IGRAPH_VIT_END(vit); IGRAPH_VIT_NEXT(vit)) {
remain[ (long int) IGRAPH_VIT_GET(vit) ] = 1;
}
for (i=0; i<no_of_nodes; i++) {
IGRAPH_ALLOW_INTERRUPTION();
if (remain[i] == 0) {
IGRAPH_CHECK(igraph_vector_push_back(&delete, i));
}
}
Free(remain);
IGRAPH_FINALLY_CLEAN(1);
/* must set res->attr to 0 before calling igraph_copy */
res->attr=0; /* Why is this needed? TODO */
IGRAPH_CHECK(igraph_copy(res, graph));
IGRAPH_FINALLY(igraph_destroy, res);
IGRAPH_CHECK(igraph_delete_vertices(res, igraph_vss_vector(&delete)));
igraph_vector_destroy(&delete);
igraph_vit_destroy(&vit);
IGRAPH_FINALLY_CLEAN(3);
return 0;
}
/**
* \ingroup structural
* \function igraph_simplify
* \brief Removes loop and/or multiple edges from the graph.
*
* \param graph The graph object.
* \param multiple Logical, if true, multiple edges will be removed.
* \param loops Logical, if true, loops (self edges) will be removed.
* \return Error code:
* \c IGRAPH_ENOMEM if we are out of memory.
*
* Time complexity: O(|V|+|E|).
*/
int igraph_simplify(igraph_t *graph, igraph_bool_t multiple, igraph_bool_t loops) {
igraph_vector_t edges=IGRAPH_VECTOR_NULL;
igraph_vector_t neis=IGRAPH_VECTOR_NULL;
long int no_of_nodes=igraph_vcount(graph);
long int i, j;
igraph_es_t es;
igraph_bool_t directed=igraph_is_directed(graph);
IGRAPH_VECTOR_INIT_FINALLY(&edges, 0);
IGRAPH_VECTOR_INIT_FINALLY(&neis, 0);
if (directed) {
for (i=0; i<no_of_nodes; i++) {
IGRAPH_CHECK(igraph_neighbors(graph, &neis, i, IGRAPH_OUT));
IGRAPH_ALLOW_INTERRUPTION();
if (loops) {
for (j=0; j<igraph_vector_size(&neis); j++) {
if (VECTOR(neis)[j]==i) {
IGRAPH_CHECK(igraph_vector_push_back(&edges, i));
IGRAPH_CHECK(igraph_vector_push_back(&edges, i));
}
}
} /* if loops */
if (multiple) {
for (j=1; j<igraph_vector_size(&neis); j++) {
if (VECTOR(neis)[j]==VECTOR(neis)[j-1] &&
(!loops || VECTOR(neis)[j] != i) ) {
IGRAPH_CHECK(igraph_vector_push_back(&edges, i));
IGRAPH_CHECK(igraph_vector_push_back(&edges, VECTOR(neis)[j]));
}
}
}
}
} else { /* not directed */
for (i=0; i<no_of_nodes; i++) {
int flip=0;
IGRAPH_CHECK(igraph_neighbors(graph, &neis, i, IGRAPH_OUT));
IGRAPH_ALLOW_INTERRUPTION();
if (loops) {
for (j=0; j<igraph_vector_size(&neis); j++) {
if (VECTOR(neis)[j]==i) {
flip=1-flip;
if (flip==0) {
IGRAPH_CHECK(igraph_vector_push_back(&edges, i));
IGRAPH_CHECK(igraph_vector_push_back(&edges, i));
}
}
}
} /* if loops */
if (multiple) {
for (j=1; j<igraph_vector_size(&neis); j++) {
if (VECTOR(neis)[j] > i && VECTOR(neis)[j]==VECTOR(neis)[j-1]) {
IGRAPH_CHECK(igraph_vector_push_back(&edges, i));
IGRAPH_CHECK(igraph_vector_push_back(&edges, VECTOR(neis)[j]));
}
if (VECTOR(neis)[j]==i && VECTOR(neis)[j-1]==i && !loops) {
flip=1-flip;
if (flip==0) {
IGRAPH_CHECK(igraph_vector_push_back(&edges, i));
IGRAPH_CHECK(igraph_vector_push_back(&edges, i));
}
}
}
}
}
}
igraph_vector_destroy(&neis);
IGRAPH_FINALLY_CLEAN(1);
IGRAPH_CHECK(igraph_es_multipairs(&es, &edges, IGRAPH_DIRECTED));
IGRAPH_FINALLY(igraph_es_destroy, &es);
IGRAPH_CHECK(igraph_delete_edges(graph, es));
igraph_es_destroy(&es);
IGRAPH_FINALLY_CLEAN(1);
igraph_vector_destroy(&edges);
IGRAPH_FINALLY_CLEAN(1);
return 0;
}
/**
* \function igraph_transitivity_avglocal_undirected
* \brief Average local transitivity (clustering coefficient)
*
* The transitivity measures the probability that two neighbors of a
* vertex are connected. In case of the average local transitivity
* this probability if calculated for each vertex and then the average
* is taken for those vertices which have at least two neighbors. If
* there are no such vertices then \c NaN is returned.
* \param graph The input graph, directed graphs are considered as
* undirected ones.
* \param res Pointer to a real variable, the result will be stored here.
* \return Error code.
*
* \sa \ref igraph_transitivity_undirected(), \ref
* igraph_transitivity_local_undirected().
*
* Time complexity: O(|V|*d^2), |V| is the number of vertices in the
* graph and d is the average degree.
*/
int igraph_transitivity_avglocal_undirected(const igraph_t *graph,
igraph_real_t *res) {
long int no_of_nodes=igraph_vcount(graph);
igraph_real_t sum=0.0;
igraph_integer_t count=0;
long int node, i, j, nn;
igraph_i_adjlist_t allneis;
igraph_vector_t *neis1, *neis2;
long int neilen1, neilen2;
igraph_integer_t triples;
long int *neis;
long int maxdegree;
igraph_vector_t order;
igraph_vector_t rank;
igraph_vector_t degree;
igraph_vector_t triangles;
IGRAPH_VECTOR_INIT_FINALLY(&order, no_of_nodes);
IGRAPH_VECTOR_INIT_FINALLY(°ree, no_of_nodes);
IGRAPH_CHECK(igraph_degree(graph, °ree, igraph_vss_all(), IGRAPH_ALL,
IGRAPH_LOOPS));
maxdegree=igraph_vector_max(°ree)+1;
igraph_vector_order1(°ree, &order, maxdegree);
igraph_vector_destroy(°ree);
IGRAPH_FINALLY_CLEAN(1);
IGRAPH_VECTOR_INIT_FINALLY(&rank, no_of_nodes);
for (i=0; i<no_of_nodes; i++) {
VECTOR(rank)[ (long int) VECTOR(order)[i] ] = no_of_nodes-i-1;
}
IGRAPH_CHECK(igraph_i_adjlist_init(graph, &allneis, IGRAPH_ALL));
IGRAPH_FINALLY(igraph_i_adjlist_destroy, &allneis);
IGRAPH_CHECK(igraph_i_adjlist_simplify(&allneis));
neis=Calloc(no_of_nodes, long int);
if (neis==0) {
IGRAPH_ERROR("undirected average local transitivity failed",
IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(igraph_free, neis);
IGRAPH_VECTOR_INIT_FINALLY(&triangles, no_of_nodes);
for (nn=no_of_nodes-1; nn >= 0; nn--) {
node=VECTOR(order)[nn];
IGRAPH_ALLOW_INTERRUPTION();
neis1=igraph_i_adjlist_get(&allneis, node);
neilen1=igraph_vector_size(neis1);
triples = (double)neilen1 * (neilen1-1) / 2;
/* Mark the neighbors of 'node' */
for (i=0; i<neilen1; i++) {
neis[ (long int)VECTOR(*neis1)[i] ] = node+1;
}
for (i=0; i<neilen1; i++) {
long int nei=VECTOR(*neis1)[i];
if (VECTOR(rank)[nei] > VECTOR(rank)[node]) {
neis2=igraph_i_adjlist_get(&allneis, nei);
neilen2=igraph_vector_size(neis2);
for (j=0; j<neilen2; j++) {
long int nei2=VECTOR(*neis2)[j];
if (VECTOR(rank)[nei2] < VECTOR(rank)[nei]) {
continue;
}
if (neis[nei2] == node+1) {
VECTOR(triangles)[nei2] += 1;
VECTOR(triangles)[nei] += 1;
VECTOR(triangles)[node] += 1;
}
}
}
}
if (triples != 0) {
sum += VECTOR(triangles)[node] / triples;
count++;
}
}
*res = sum/count;
igraph_vector_destroy(&triangles);
Free(neis);
igraph_i_adjlist_destroy(&allneis);
igraph_vector_destroy(&rank);
igraph_vector_destroy(&order);
IGRAPH_FINALLY_CLEAN(5);
return 0;
}
int igraph_transitivity_local_undirected1(const igraph_t *graph,
igraph_vector_t *res,
const igraph_vs_t vids) {
long int no_of_nodes=igraph_vcount(graph);
igraph_vit_t vit;
long int nodes_to_calc;
igraph_vector_t *neis1, *neis2;
igraph_real_t triples, triangles;
long int i, j, k;
long int neilen1, neilen2;
long int *neis;
igraph_i_lazy_adjlist_t adjlist;
IGRAPH_CHECK(igraph_vit_create(graph, vids, &vit));
IGRAPH_FINALLY(igraph_vit_destroy, &vit);
nodes_to_calc=IGRAPH_VIT_SIZE(vit);
neis=Calloc(no_of_nodes, long int);
if (neis==0) {
IGRAPH_ERROR("local undirected transitivity failed", IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(igraph_free, neis);
IGRAPH_CHECK(igraph_vector_resize(res, nodes_to_calc));
igraph_i_lazy_adjlist_init(graph, &adjlist, IGRAPH_ALL, IGRAPH_I_SIMPLIFY);
IGRAPH_FINALLY(igraph_i_lazy_adjlist_destroy, &adjlist);
for (i=0; !IGRAPH_VIT_END(vit); IGRAPH_VIT_NEXT(vit), i++) {
long int node=IGRAPH_VIT_GET(vit);
IGRAPH_ALLOW_INTERRUPTION();
neis1=igraph_i_lazy_adjlist_get(&adjlist, node);
neilen1=igraph_vector_size(neis1);
for (j=0; j<neilen1; j++) {
neis[ (long int)VECTOR(*neis1)[j] ] = i+1;
}
triples = (double)neilen1*(neilen1-1);
triangles = 0;
for (j=0; j<neilen1; j++) {
long int v=VECTOR(*neis1)[j];
neis2=igraph_i_lazy_adjlist_get(&adjlist, v);
neilen2=igraph_vector_size(neis2);
for (k=0; k<neilen2; k++) {
long int v2=VECTOR(*neis2)[k];
if (neis[v2] == i+1) {
triangles += 1.0;
}
}
}
VECTOR(*res)[i] = triangles/triples;
/* fprintf(stderr, "%f %f\n", triangles, triples); */
}
igraph_i_lazy_adjlist_destroy(&adjlist);
Free(neis);
igraph_vit_destroy(&vit);
IGRAPH_FINALLY_CLEAN(3);
return 0;
}
int igraph_transitivity_local_undirected2(const igraph_t *graph,
igraph_vector_t *res,
const igraph_vs_t vids) {
long int no_of_nodes=igraph_vcount(graph);
igraph_vit_t vit;
long int nodes_to_calc, affected_nodes;
long int maxdegree=0;
long int i, j, k, nn;
igraph_i_lazy_adjlist_t adjlist;
igraph_vector_t index, avids, rank, order, triangles, degree;
long int *neis;
IGRAPH_CHECK(igraph_vit_create(graph, vids, &vit));
IGRAPH_FINALLY(igraph_vit_destroy, &vit);
nodes_to_calc=IGRAPH_VIT_SIZE(vit);
IGRAPH_CHECK(igraph_i_lazy_adjlist_init(graph, &adjlist, IGRAPH_ALL,
IGRAPH_I_SIMPLIFY));
IGRAPH_FINALLY(igraph_i_lazy_adjlist_destroy, &adjlist);
IGRAPH_VECTOR_INIT_FINALLY(&index, no_of_nodes);
IGRAPH_VECTOR_INIT_FINALLY(&avids, 0);
IGRAPH_CHECK(igraph_vector_reserve(&avids, nodes_to_calc));
k=0;
for (i=0; i<nodes_to_calc; IGRAPH_VIT_NEXT(vit), i++) {
long int v=IGRAPH_VIT_GET(vit);
igraph_vector_t *neis;
long int neilen;
if (VECTOR(index)[v]==0) {
VECTOR(index)[v]=k+1; k++;
IGRAPH_CHECK(igraph_vector_push_back(&avids, v));
}
neis=igraph_i_lazy_adjlist_get(&adjlist, v);
neilen=igraph_vector_size(neis);
for (j=0; j<neilen; j++) {
long int nei=VECTOR(*neis)[j];
if (VECTOR(index)[nei]==0) {
VECTOR(index)[nei]=k+1; k++;
IGRAPH_CHECK(igraph_vector_push_back(&avids, nei));
}
}
}
/* Degree, ordering, ranking */
affected_nodes=igraph_vector_size(&avids);
IGRAPH_VECTOR_INIT_FINALLY(&order, 0);
IGRAPH_VECTOR_INIT_FINALLY(°ree, affected_nodes);
for (i=0; i<affected_nodes; i++) {
long int v=VECTOR(avids)[i];
igraph_vector_t *neis;
long int deg;
neis=igraph_i_lazy_adjlist_get(&adjlist, v);
VECTOR(degree)[i]=deg=igraph_vector_size(neis);
if (deg > maxdegree) { maxdegree = deg; }
}
igraph_vector_order1(°ree, &order, maxdegree+1);
igraph_vector_destroy(°ree);
IGRAPH_FINALLY_CLEAN(1);
IGRAPH_VECTOR_INIT_FINALLY(&rank, affected_nodes);
for (i=0; i<affected_nodes; i++) {
VECTOR(rank)[ (long int) VECTOR(order)[i] ] = affected_nodes-i-1;
}
neis=Calloc(no_of_nodes, long int);
if (neis==0) {
IGRAPH_ERROR("local transitivity calculation failed", IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(igraph_free, neis);
IGRAPH_VECTOR_INIT_FINALLY(&triangles, affected_nodes);
for (nn=affected_nodes-1; nn>=0; nn--) {
long int node=VECTOR(avids) [ (long int) VECTOR(order)[nn] ];
igraph_vector_t *neis1, *neis2;
long int neilen1, neilen2;
long int nodeindex=VECTOR(index)[node];
long int noderank=VECTOR(rank) [nodeindex-1];
/* fprintf(stderr, "node %li (index %li, rank %li)\n", node, */
/* (long int)VECTOR(index)[node]-1, noderank); */
IGRAPH_ALLOW_INTERRUPTION();
neis1=igraph_i_lazy_adjlist_get(&adjlist, node);
neilen1=igraph_vector_size(neis1);
for (i=0; i<neilen1; i++) {
long int nei=VECTOR(*neis1)[i];
neis[nei] = node+1;
}
for (i=0; i<neilen1; i++) {
long int nei=VECTOR(*neis1)[i];
long int neiindex=VECTOR(index)[nei];
long int neirank=VECTOR(rank)[neiindex-1];
/* fprintf(stderr, " nei %li (index %li, rank %li)\n", nei, */
/* neiindex, neirank); */
if (neirank > noderank) {
neis2=igraph_i_lazy_adjlist_get(&adjlist, nei);
neilen2=igraph_vector_size(neis2);
for (j=0; j<neilen2; j++) {
long int nei2=VECTOR(*neis2)[j];
long int nei2index=VECTOR(index)[nei2];
long int nei2rank=VECTOR(rank)[nei2index-1];
/* fprintf(stderr, " triple %li %li %li\n", node, nei, nei2); */
if (nei2rank < neirank) {
continue;
}
if (neis[nei2] == node+1) {
/* fprintf(stderr, " triangle\n"); */
VECTOR(triangles) [ nei2index-1 ] += 1;
VECTOR(triangles) [ neiindex-1 ] += 1;
VECTOR(triangles) [ nodeindex-1 ] += 1;
}
}
}
}
}
/* Ok, for all affected vertices the number of triangles were counted */
IGRAPH_CHECK(igraph_vector_resize(res, nodes_to_calc));
IGRAPH_VIT_RESET(vit);
for (i=0; i<nodes_to_calc; i++, IGRAPH_VIT_NEXT(vit)) {
long int node=IGRAPH_VIT_GET(vit);
long int idx=VECTOR(index)[node]-1;
igraph_vector_t *neis=igraph_i_lazy_adjlist_get(&adjlist, node);
long int deg=igraph_vector_size(neis);
igraph_real_t triples=(double) deg * (deg-1) / 2;
VECTOR(*res)[i] = VECTOR(triangles)[idx] / triples;
/* fprintf(stderr, "%f %f\n", VECTOR(triangles)[idx], triples); */
}
igraph_vector_destroy(&triangles);
igraph_free(neis);
igraph_vector_destroy(&rank);
igraph_vector_destroy(&order);
igraph_vector_destroy(&avids);
igraph_vector_destroy(&index);
igraph_i_lazy_adjlist_destroy(&adjlist);
igraph_vit_destroy(&vit);
IGRAPH_FINALLY_CLEAN(8);
return 0;
}
/* We don't use this, it is theoretically good, but practically not.
*/
/* int igraph_transitivity_local_undirected3(const igraph_t *graph, */
/* igraph_vector_t *res, */
/* const igraph_vs_t vids) { */
/* igraph_vit_t vit; */
/* long int nodes_to_calc; */
/* igraph_i_lazy_adjlist_t adjlist; */
/* long int i, j; */
/* IGRAPH_CHECK(igraph_vit_create(graph, vids, &vit)); */
/* IGRAPH_FINALLY(igraph_vit_destroy, &vit); */
/* nodes_to_calc=IGRAPH_VIT_SIZE(vit); */
/* IGRAPH_CHECK(igraph_i_lazy_adjlist_init(graph, &adjlist, IGRAPH_ALL, */
/* IGRAPH_I_SIMPLIFY)); */
/* IGRAPH_FINALLY(igraph_i_lazy_adjlist_destroy, &adjlist); */
/* IGRAPH_CHECK(igraph_vector_resize(res, nodes_to_calc)); */
/* for (i=0, IGRAPH_VIT_RESET(vit); !IGRAPH_VIT_END(vit); */
/* i++, IGRAPH_VIT_NEXT(vit)) { */
/* long int node=IGRAPH_VIT_GET(vit); */
/* igraph_vector_t *neis=igraph_i_lazy_adjlist_get(&adjlist, node); */
/* long int n1=igraph_vector_size(neis); */
/* igraph_real_t triangles=0; */
/* igraph_real_t triples=(double)n1*(n1-1); */
/* IGRAPH_ALLOW_INTERRUPTION(); */
/* for (j=0; j<n1; j++) { */
/* long int node2=VECTOR(*neis)[j]; */
/* igraph_vector_t *neis2=igraph_i_lazy_adjlist_get(&adjlist, node2); */
/* long int n2=igraph_vector_size(neis2); */
/* long int l1=0, l2=0; */
/* while (l1 < n1 && l2 < n2) { */
/* long int nei1=VECTOR(*neis)[l1]; */
/* long int nei2=VECTOR(*neis2)[l2]; */
/* if (nei1 < nei2) { */
/* l1++; */
/* } else if (nei1 > nei2) { */
/* l2++; */
/* } else { */
/* triangles+=1; */
/* l1++; l2++; */
/* } */
/* } */
/* } */
/* /\* We're done with 'node' *\/ */
/* VECTOR(*res)[i] = triangles / triples; */
/* } */
/* igraph_i_lazy_adjlist_destroy(&adjlist); */
/* igraph_vit_destroy(&vit); */
/* IGRAPH_FINALLY_CLEAN(2); */
/* return 0; */
/* } */
int igraph_transitivity_local_undirected4(const igraph_t *graph,
igraph_vector_t *res,
const igraph_vs_t vids) {
long int no_of_nodes=igraph_vcount(graph);
long int node, i, j, nn;
igraph_i_adjlist_t allneis;
igraph_vector_t *neis1, *neis2;
long int neilen1, neilen2;
igraph_integer_t triples;
long int *neis;
long int maxdegree;
igraph_vector_t order;
igraph_vector_t rank;
igraph_vector_t degree;
if (!igraph_vs_is_all((igraph_vs_t*)&vids)) {
IGRAPH_ERROR("Internal error, wrong transitivity function called",
IGRAPH_EINVAL);
}
IGRAPH_VECTOR_INIT_FINALLY(&order, no_of_nodes);
IGRAPH_VECTOR_INIT_FINALLY(°ree, no_of_nodes);
IGRAPH_CHECK(igraph_degree(graph, °ree, igraph_vss_all(), IGRAPH_ALL,
IGRAPH_LOOPS));
maxdegree=igraph_vector_max(°ree)+1;
igraph_vector_order1(°ree, &order, maxdegree);
igraph_vector_destroy(°ree);
IGRAPH_FINALLY_CLEAN(1);
IGRAPH_VECTOR_INIT_FINALLY(&rank, no_of_nodes);
for (i=0; i<no_of_nodes; i++) {
VECTOR(rank)[ (long int)VECTOR(order)[i] ] = no_of_nodes-i-1;
}
IGRAPH_CHECK(igraph_i_adjlist_init(graph, &allneis, IGRAPH_ALL));
IGRAPH_FINALLY(igraph_i_adjlist_destroy, &allneis);
IGRAPH_CHECK(igraph_i_adjlist_simplify(&allneis));
neis=Calloc(no_of_nodes, long int);
if (neis==0) {
IGRAPH_ERROR("undirected local transitivity failed", IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(igraph_free, neis);
IGRAPH_CHECK(igraph_vector_resize(res, no_of_nodes));
igraph_vector_null(res);
for (nn=no_of_nodes-1; nn>=0; nn--) {
node=VECTOR(order)[nn];
IGRAPH_ALLOW_INTERRUPTION();
neis1=igraph_i_adjlist_get(&allneis, node);
neilen1=igraph_vector_size(neis1);
triples=(double)neilen1*(neilen1-1)/2;
/* Mark the neighbors of the node */
for (i=0; i<neilen1; i++) {
neis[ (long int) VECTOR(*neis1)[i] ] = node+1;
}
for (i=0; i<neilen1; i++) {
long int nei=VECTOR(*neis1)[i];
if (VECTOR(rank)[nei] > VECTOR(rank)[node]) {
neis2=igraph_i_adjlist_get(&allneis, nei);
neilen2=igraph_vector_size(neis2);
for (j=0; j<neilen2; j++) {
long int nei2=VECTOR(*neis2)[j];
if (VECTOR(rank)[nei2] < VECTOR(rank)[nei]) {
continue;
}
if (neis[nei2] == node+1) {
VECTOR(*res)[nei2] += 1;
VECTOR(*res)[nei] += 1;
VECTOR(*res)[node] += 1;
}
}
}
}
VECTOR(*res)[node] /= triples;
}
igraph_free(neis);
igraph_i_adjlist_destroy(&allneis);
igraph_vector_destroy(&rank);
igraph_vector_destroy(&order);
IGRAPH_FINALLY_CLEAN(4);
return 0;
}
/**
* \function igraph_transitivity_local_undirected
* \brief Calculates the local transitivity (clustering coefficient)
* of a graph
*
* The transitivity measures the probability that two neighbors of a
* vertex are connected. In case of the local transitivity, this
* probability is calculated separately for each vertex.
* \param graph The input graph, it can be directed but direction of
* the edges will be ignored.
* \param res Pointer to an initialized vector, the result will be
* stored here. It will be resized as needed.
* \param vids Vertex set, the vertices for which the local
* transitivity will be calculated.
* \return Error code.
*
* \sa \ref igraph_transitivity_undirected(), \ref
* igraph_transitivity_avglocal_undirected().
*
* Time complexity: O(n*d^2), n is the number of vertices for which
* the transitivity is calculated, d is the average vertex degree.
*/
int igraph_transitivity_local_undirected(const igraph_t *graph,
igraph_vector_t *res,
const igraph_vs_t vids) {
if (igraph_vs_is_all((igraph_vs_t*)&vids)) {
return igraph_transitivity_local_undirected4(graph, res, vids);
} else {
igraph_vit_t vit;
long int size;
IGRAPH_CHECK(igraph_vit_create(graph, vids, &vit));
IGRAPH_FINALLY(igraph_vit_destroy, &vit);
size=IGRAPH_VIT_SIZE(vit);
igraph_vit_destroy(&vit);
if (size < 100) {
return igraph_transitivity_local_undirected1(graph, res, vids);
} else {
return igraph_transitivity_local_undirected2(graph, res, vids);
}
}
return 0;
}
/**
* \ingroup structural
* \function igraph_transitivity_undirected
* \brief Calculates the transitivity (clustering coefficient) of a graph.
*
* </para><para>
* The transitivity measures the probability that two neighbors of a
* vertex are connected. More precisely this is the ratio of the
* triangles and connected triples in the graph, the result is a
* single real number or NaN (0/0) if there are no connected triples
* in the graph. Directed graphs are considered as undirected ones.
* \param graph The graph object.
* \param res Pointer to a real variable, the result will be stored here.
* \return Error code:
* \c IGRAPH_ENOMEM: not enough memory for
* temporary data.
*
* \sa \ref igraph_transitivity_local_undirected(),
* \ref igraph_transitivity_avglocal_undirected().
*
* Time complexity: O(|V|*d^2), |V| is the number of vertices in
* the graph, d is the average node degree.
*/
int igraph_transitivity_undirected(const igraph_t *graph,
igraph_real_t *res) {
long int no_of_nodes=igraph_vcount(graph);
igraph_real_t triples=0, triangles=0;
long int node, nn;
long int maxdegree;
long int *neis;
igraph_vector_t order;
igraph_vector_t rank;
igraph_vector_t degree;
igraph_i_adjlist_t allneis;
igraph_vector_t *neis1, *neis2;
long int i, j, neilen1, neilen2;
IGRAPH_VECTOR_INIT_FINALLY(&order, no_of_nodes);
IGRAPH_VECTOR_INIT_FINALLY(°ree, no_of_nodes);
IGRAPH_CHECK(igraph_degree(graph, °ree, igraph_vss_all(), IGRAPH_ALL,
IGRAPH_LOOPS));
maxdegree=igraph_vector_max(°ree)+1;
igraph_vector_order1(°ree, &order, maxdegree);
igraph_vector_destroy(°ree);
IGRAPH_FINALLY_CLEAN(1);
IGRAPH_VECTOR_INIT_FINALLY(&rank, no_of_nodes);
for (i=0; i<no_of_nodes; i++) {
VECTOR(rank)[ (long int) VECTOR(order)[i] ]=no_of_nodes-i-1;
}
IGRAPH_CHECK(igraph_i_adjlist_init(graph, &allneis, IGRAPH_ALL));
IGRAPH_FINALLY(igraph_i_adjlist_destroy, &allneis);
IGRAPH_CHECK(igraph_i_adjlist_simplify(&allneis));
neis=Calloc(no_of_nodes, long int);
if (neis==0) {
IGRAPH_ERROR("undirected transitivity failed", IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(igraph_free, neis);
for (nn=no_of_nodes-1; nn >=0; nn--) {
node=VECTOR(order)[nn];
IGRAPH_ALLOW_INTERRUPTION();
neis1=igraph_i_adjlist_get(&allneis, node);
neilen1=igraph_vector_size(neis1);
triples += (double)neilen1 * (neilen1-1);
/* Mark the neighbors of 'node' */
for (i=0; i<neilen1; i++) {
long int nei=VECTOR(*neis1)[i];
neis[nei] = node+1;
}
for (i=0; i<neilen1; i++) {
long int nei=VECTOR(*neis1)[i];
/* If 'nei' is not ready yet */
if (VECTOR(rank)[nei] > VECTOR(rank)[node]) {
neis2=igraph_i_adjlist_get(&allneis, nei);
neilen2=igraph_vector_size(neis2);
for (j=0; j<neilen2; j++) {
long int nei2=VECTOR(*neis2)[j];
if (neis[nei2] == node+1) {
triangles += 1.0;
}
}
}
}
}
Free(neis);
igraph_i_adjlist_destroy(&allneis);
igraph_vector_destroy(&rank);
igraph_vector_destroy(&order);
IGRAPH_FINALLY_CLEAN(4);
*res = triangles / triples * 2.0;
return 0;
}
/**
* \ingroup structural
* \function igraph_reciprocity
* \brief Calculates the reciprocity of a directed graph.
*
* </para><para>
* A vertex pair (A, B) is said to be reciprocal if there are edges
* between them in both directions. The reciprocity of a directed graph
* is the proportion of all possible (A, B) pairs which are reciprocal,
* provided there is at least one edge between A and B. The reciprocity
* of an empty graph is undefined (results in an error code). Undirected
* graphs always have a reciprocity of 1.0 unless they are empty.
*
* \param graph The graph object.
* \param res Pointer to an \c igraph_real_t which will contain the result.
* \param ignore_loops Whether to ignore loop edges.
* \return Error code:
* \c IGRAPH_EINVAL: graph has no edges
* \c IGRAPH_ENOMEM: not enough memory for
* temporary data.
*
* Time complexity: O(|V|+|E|), |V| is the number of vertices,
* |E| is the number of edges.
*/
int igraph_reciprocity(const igraph_t *graph, igraph_real_t *res,
igraph_bool_t ignore_loops) {
igraph_integer_t nonrec=0, rec=0;
igraph_vector_t inneis, outneis;
long int i;
long int no_of_nodes=igraph_vcount(graph);
/* THIS IS AN EXIT HERE !!!!!!!!!!!!!! */
if (!igraph_is_directed(graph)) {
*res=1.0;
return 0;
}
IGRAPH_VECTOR_INIT_FINALLY(&inneis, 0);
IGRAPH_VECTOR_INIT_FINALLY(&outneis, 0);
for (i=0; i<no_of_nodes; i++) {
long int ip, op;
igraph_neighbors(graph, &inneis, i, IGRAPH_IN);
igraph_neighbors(graph, &outneis, i, IGRAPH_OUT);
ip=op=0;
while (ip < igraph_vector_size(&inneis) &&
op < igraph_vector_size(&outneis)) {
if (VECTOR(inneis)[ip] < VECTOR(outneis)[op]) {
nonrec += 1;
ip++;
} else if (VECTOR(inneis)[ip] > VECTOR(outneis)[op]) {
nonrec += 1;
op++;
} else {
/* loop edge? */
if (!ignore_loops || VECTOR(inneis)[ip] != i) {
rec += 1;
}
ip++;
op++;
}
}
nonrec += (igraph_vector_size(&inneis)-ip) +
(igraph_vector_size(&outneis)-op);
}
*res= rec/(rec+nonrec);
igraph_vector_destroy(&inneis);
igraph_vector_destroy(&outneis);
IGRAPH_FINALLY_CLEAN(2);
return 0;
}
/**
* \function igraph_constraint
* \brief Burt's constraint scores
*
* </para><para>
* This function calculates Burt's constraint scores for the given
* vertices, also known as structural holes.
*
* </para><para>
* Burt's constraint is higher if ego has less, or mutually stronger
* related (i.e. more redundant) contacts. Burt's measure of
* constraint, C[i], of vertex i's ego network V[i], is defined for
* directed and valued graphs,
* <blockquote><para>
* C[i] = sum( sum( (p[i,q] p[q,j])^2, q in V[i], q != i,j ), j in
* V[], j != i)
* </para></blockquote>
* for a graph of order (ie. number od vertices) N, where proportional
* tie strengths are defined as
* <blockquote><para>
* p[i,j]=(a[i,j]+a[j,i]) / sum(a[i,k]+a[k,i], k in V[i], k != i),
* </para></blockquote>
* a[i,j] are elements of A and
* the latter being the graph adjacency matrix. For isolated vertices,
* constraint is undefined.
*
* </para><para>
* Burt, R.S. (2004). Structural holes and good ideas. American
* Journal of Sociology 110, 349-399.
*
* </para><para>
* The first R version of this function was contributed by Jeroen
* Bruggeman.
* \param graph A graph object.
* \param res Pointer to an initialized vector, the result will be
* stored here. The vector will be resized to have the
* appropriate size for holding the result.
* \param vids Vertex selector containing the vertices for which the
* constraint should be calculated.
* \param weights Vector giving the weights of the edges. If it is
* \c NULL then each edge is supposed to have the same weight.
* \return Error code.
*
* Time complexity: O(|V|+E|+n*d^2), n is the number of vertices for
* which the constraint is calculated and d is the average degree, |V|
* is the number of vertices, |E| the number of edges in the
* graph. If the weights argument is \c NULL then the time complexity
* is O(|V|+n*d^2).
*/
int igraph_constraint(const igraph_t *graph, igraph_vector_t *res,
igraph_vs_t vids, const igraph_vector_t *weights) {
long int no_of_nodes=igraph_vcount(graph);
long int no_of_edges=igraph_ecount(graph);
igraph_vit_t vit;
long int nodes_to_calc;
long int a, b, c, i, j, q;
igraph_integer_t edge, from, to, edge2, from2, to2;
igraph_vector_t contrib;
igraph_vector_t degree;
igraph_vector_t ineis_in, ineis_out, jneis_in, jneis_out;
if (weights != 0 && igraph_vector_size(weights) != no_of_edges) {
IGRAPH_ERROR("Invalid length of weight vector", IGRAPH_EINVAL);
}
IGRAPH_VECTOR_INIT_FINALLY(&contrib, no_of_nodes);
IGRAPH_VECTOR_INIT_FINALLY(°ree, no_of_nodes);
IGRAPH_VECTOR_INIT_FINALLY(&ineis_in, 0);
IGRAPH_VECTOR_INIT_FINALLY(&ineis_out, 0);
IGRAPH_VECTOR_INIT_FINALLY(&jneis_in, 0);
IGRAPH_VECTOR_INIT_FINALLY(&jneis_out, 0);
IGRAPH_CHECK(igraph_vit_create(graph, vids, &vit));
IGRAPH_FINALLY(igraph_vit_destroy, &vit);
nodes_to_calc=IGRAPH_VIT_SIZE(vit);
if (weights==0) {
IGRAPH_CHECK(igraph_degree(graph, °ree, igraph_vss_all(),
IGRAPH_ALL, IGRAPH_NO_LOOPS));
} else {
for (a=0; a<no_of_edges; a++) {
igraph_edge(graph, a, &from, &to);
if (from != to) {
VECTOR(degree)[(long int) from] += VECTOR(*weights)[a];
VECTOR(degree)[(long int) to ] += VECTOR(*weights)[a];
}
}
}
IGRAPH_CHECK(igraph_vector_resize(res, nodes_to_calc));
igraph_vector_null(res);
for (a=0; a<nodes_to_calc; a++, IGRAPH_VIT_NEXT(vit)) {
i=IGRAPH_VIT_GET(vit);
/* get neighbors of i */
IGRAPH_CHECK(igraph_adjacent(graph, &ineis_in, i, IGRAPH_IN));
IGRAPH_CHECK(igraph_adjacent(graph, &ineis_out, i, IGRAPH_OUT));
/* NaN for isolates */
if (igraph_vector_size(&ineis_in) == 0 &&
igraph_vector_size(&ineis_out) == 0) {
VECTOR(*res)[a] = 0.0 / 0.0;
}
/* zero their contribution */
for (b=0; b<igraph_vector_size(&ineis_in); b++) {
edge=VECTOR(ineis_in)[b];
igraph_edge(graph, edge, &from, &to);
if (to==i) { to=from; }
j=to;
VECTOR(contrib)[j]=0.0;
}
for (b=0; b<igraph_vector_size(&ineis_out); b++) {
edge=VECTOR(ineis_out)[b];
igraph_edge(graph, edge, &from, &to);
if (to==i) { to=from; }
j=to;
VECTOR(contrib)[j]=0.0;
}
/* add the direct contibutions, in-neighbors and out-neighbors */
for (b=0; b<igraph_vector_size(&ineis_in); b++) {
edge=VECTOR(ineis_in)[b];
igraph_edge(graph, edge, &from, &to);
if (to==i) { to=from; }
j=to;
if (i != j) { /* excluding loops */
if (weights) {
VECTOR(contrib)[j] +=
VECTOR(*weights)[(long int)edge]/VECTOR(degree)[i];
} else {
VECTOR(contrib)[j] += 1.0/VECTOR(degree)[i];
}
}
}
if (igraph_is_directed(graph)) {
for (b=0; b<igraph_vector_size(&ineis_out); b++) {
edge=VECTOR(ineis_out)[b];
igraph_edge(graph, edge, &from, &to);
if (to==i) { to=from; }
j=to;
if (i != j) {
if (weights) {
VECTOR(contrib)[j] +=
VECTOR(*weights)[(long int)edge]/VECTOR(degree)[i];
} else {
VECTOR(contrib)[j] += 1.0/VECTOR(degree)[i];
}
}
}
}
/* add the indirect contributions, in-in, in-out, out-in, out-out */
for (b=0; b<igraph_vector_size(&ineis_in); b++) {
edge=VECTOR(ineis_in)[b];
igraph_edge(graph, edge, &from, &to);
if (to==i) { to=from; }
j=to;
if (i == j) { continue; }
IGRAPH_CHECK(igraph_adjacent(graph, &jneis_in, j, IGRAPH_IN));
IGRAPH_CHECK(igraph_adjacent(graph, &jneis_out, j, IGRAPH_OUT));
for (c=0; c<igraph_vector_size(&jneis_in); c++) {
edge2=VECTOR(jneis_in)[c];
igraph_edge(graph, edge2, &from2, &to2);
if (to2==j) { to2=from2; }
q=to2;
if (j != q) {
if (weights) {
VECTOR(contrib)[q] +=
VECTOR(*weights)[(long int)edge]*
VECTOR(*weights)[(long int)edge2]/
VECTOR(degree)[i]/VECTOR(degree)[j];
} else {
VECTOR(contrib)[q] += 1/VECTOR(degree)[i]/VECTOR(degree)[j];
}
}
}
if (igraph_is_directed(graph)) {
for (c=0; c<igraph_vector_size(&jneis_out); c++) {
edge2=VECTOR(jneis_out)[c];
igraph_edge(graph, edge2, &from2, &to2);
if (to2==j) { to2=from2; }
q=to2;
if (j != q) {
if (weights) {
VECTOR(contrib)[q] +=
VECTOR(*weights)[(long int)edge]*
VECTOR(*weights)[(long int)edge2]/
VECTOR(degree)[i]/VECTOR(degree)[j];
} else {
VECTOR(contrib)[q] += 1/VECTOR(degree)[i]/VECTOR(degree)[j];
}
}
}
}
}
if (igraph_is_directed(graph)) {
for (b=0; b<igraph_vector_size(&ineis_out); b++) {
edge=VECTOR(ineis_out)[b];
igraph_edge(graph, edge, &from, &to);
if (to==i) { to=from; }
j=to;
if (i == j) { continue; }
IGRAPH_CHECK(igraph_adjacent(graph, &jneis_in, j, IGRAPH_IN));
IGRAPH_CHECK(igraph_adjacent(graph, &jneis_out, j, IGRAPH_OUT));
for (c=0; c<igraph_vector_size(&jneis_in); c++) {
edge2=VECTOR(jneis_in)[c];
igraph_edge(graph, edge2, &from2, &to2);
if (to2==j) { to2=from2; }
q=to2;
if (j != q) {
if (weights) {
VECTOR(contrib)[q] +=
VECTOR(*weights)[(long int)edge]*
VECTOR(*weights)[(long int)edge2]/
VECTOR(degree)[i]/VECTOR(degree)[j];
} else {
VECTOR(contrib)[q] += 1/VECTOR(degree)[i]/VECTOR(degree)[j];
}
}
}
for (c=0; c<igraph_vector_size(&jneis_out); c++) {
edge2=VECTOR(jneis_out)[c];
igraph_edge(graph, edge2, &from2, &to2);
if (to2==j) { to2=from2; }
q=to2;
if (j != q) {
if (weights) {
VECTOR(contrib)[q] +=
VECTOR(*weights)[(long int)edge]*
VECTOR(*weights)[(long int)edge2]/
VECTOR(degree)[i]/VECTOR(degree)[j];
} else {
VECTOR(contrib)[q] += 1/VECTOR(degree)[i]/VECTOR(degree)[j];
}
}
}
}
}
/* squared sum of the contributions */
for (b=0; b<igraph_vector_size(&ineis_in); b++) {
edge=VECTOR(ineis_in)[b];
igraph_edge(graph, edge, &from, &to);
if (to==i) { to=from; }
j=to;
if (i == j) { continue; }
VECTOR(*res)[a] += VECTOR(contrib)[j] * VECTOR(contrib)[j];
VECTOR(contrib)[j]=0.0;
}
if (igraph_is_directed(graph)) {
for (b=0; b<igraph_vector_size(&ineis_out); b++) {
edge=VECTOR(ineis_out)[b];
igraph_edge(graph, edge, &from, &to);
if (to==i) { to=from; }
j=to;
if (i == j) { continue; }
VECTOR(*res)[a] += VECTOR(contrib)[j] * VECTOR(contrib)[j];
VECTOR(contrib)[j]=0.0;
}
}
}
igraph_vit_destroy(&vit);
igraph_vector_destroy(&jneis_out);
igraph_vector_destroy(&jneis_in);
igraph_vector_destroy(&ineis_out);
igraph_vector_destroy(&ineis_in);
igraph_vector_destroy(°ree);
igraph_vector_destroy(&contrib);
IGRAPH_FINALLY_CLEAN(7);
return 0;
}
/**
* \function igraph_maxdegree
* \brief Calculate the maximum degree in a graph (or set of vertices).
*
* </para><para>
* The largest in-, out- or total degree of the specified vertices is
* calculated.
* \param graph The input graph.
* \param res Pointer to an integer (\c igraph_integer_t), the result
* will be stored here.
* \param mode Defines the type of the degree.
* \c IGRAPH_OUT, out-degree,
* \c IGRAPH_IN, in-degree,
* \c IGRAPH_ALL, total degree (sum of the
* in- and out-degree).
* This parameter is ignored for undirected graphs.
* \param loops Boolean, gives whether the self-loops should be
* counted.
* \return Error code:
* \c IGRAPH_EINVVID: invalid vertex id.
* \c IGRAPH_EINVMODE: invalid mode argument.
*
* Time complexity: O(v) if
* loops is
* TRUE, and
* O(v*d)
* otherwise. v is the number
* vertices for which the degree will be calculated, and
* d is their (average) degree.
*/
int igraph_maxdegree(const igraph_t *graph, igraph_integer_t *res,
igraph_vs_t vids, igraph_neimode_t mode,
igraph_bool_t loops) {
igraph_vector_t tmp;
IGRAPH_VECTOR_INIT_FINALLY(&tmp, 0);
igraph_degree(graph, &tmp, vids, mode, loops);
*res=igraph_vector_max(&tmp);
igraph_vector_destroy(&tmp);
IGRAPH_FINALLY_CLEAN(1);
return 0;
}
/**
* \function igraph_density
* Calculate the density of a graph.
*
* </para><para>The density of a graph is simply the ratio number of
* edges and the number of possible edges. Note that density is
* ill-defined for graphs with multiple and/or loop edges, so consider
* calling \ref igraph_simplify() on the graph if you know that it
* contains multiple or loop edges.
* \param graph The input graph object.
* \param res Pointer to a real number, the result will be stored
* here.
* \param loops Logical constant, whether to include loops in the
* calculation. If this constant is TRUE then
* loop edges are thought to be possible in the graph (this does not
* neccessary means that the graph really contains any loops). If
* this FALSE then the result is only correct if the graph does not
* contain loops.
* \return Error code.
*
* Time complexity: O(1).
*/
int igraph_density(const igraph_t *graph, igraph_real_t *res,
igraph_bool_t loops) {
igraph_integer_t no_of_nodes=igraph_vcount(graph);
igraph_integer_t no_of_edges=igraph_ecount(graph);
igraph_bool_t directed=igraph_is_directed(graph);
if (!loops) {
if (directed) {
*res = no_of_edges / (no_of_nodes*(no_of_nodes-1));
} else {
*res = no_of_edges / (no_of_nodes*(no_of_nodes-1)/2);
}
} else {
if (directed) {
*res = no_of_edges / (no_of_nodes*no_of_nodes);
} else {
*res = no_of_edges / (no_of_nodes*no_of_nodes/2);
}
}
return 0;
}
/**
* \function igraph_neighborhood_size
* \brief Calculates the size of the neighborhood of a given vertex
*
* The neighborhood of a given order of a vertex includes all vertices
* which are closer to the vertex than the order. Ie. order 0 is
* always the vertex itself, order 1 is the vertex plus its immediate
* neighbors, order 2 is order 1 plus the immediate neighbors of the
* vertices in order 1, etc.
*
* </para><para> This function calculates the size of the neighborhood
* of the given order for the given vertices.
* \param graph The input graph.
* \param res Pointer to an initialized vector, the result will be
* stored here. It will be resized as needed.
* \param vids The vertices for which the calculation is performed.
* \param order Integer giving the order of the neighborhood.
* \param mode Specifies how to use the direction of the edges if a
* directed graph is analyzed. For \c IGRAPH_OUT only the outgoing
* edges are followed, so all vertices reachable from the source
* vertex in at most \c order steps are counted. For \c IGRAPH_IN
* all vertices from which the source vertex is reachable in at most
* \c order steps are counted. \c IGRAPH_ALL ignores the direction
* of the edges. This argument is ignored for undirected graphs.
* \return Error code.
*
* \sa \ref igraph_neighborhood() for calculating the actual neighborhood,
* \ref igraph_neighborhood_graphs() for creating separate graphs from
* the neighborhoods.
*
* Time complexity: O(n*d*o), where n is the number vertices for which
* the calculation is performed, d is the average degree, o is the order.
*/
int igraph_neighborhood_size(const igraph_t *graph, igraph_vector_t *res,
igraph_vs_t vids, igraph_integer_t order,
igraph_neimode_t mode) {
long int no_of_nodes=igraph_vcount(graph);
igraph_dqueue_t q;
igraph_vit_t vit;
long int i, j;
long int *added;
igraph_vector_t neis;
if (order < 0) {
IGRAPH_ERROR("Negative order in neighborhood size", IGRAPH_EINVAL);
}
added=Calloc(no_of_nodes, long int);
if (added==0) {
IGRAPH_ERROR("Cannot calculate neighborhood size", IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(igraph_free, added);
IGRAPH_DQUEUE_INIT_FINALLY(&q, 100);
IGRAPH_CHECK(igraph_vit_create(graph, vids, &vit));
IGRAPH_FINALLY(igraph_vit_destroy, &vit);
IGRAPH_VECTOR_INIT_FINALLY(&neis, 0);
IGRAPH_CHECK(igraph_vector_resize(res, IGRAPH_VIT_SIZE(vit)));
for (i=0; !IGRAPH_VIT_END(vit); IGRAPH_VIT_NEXT(vit), i++) {
long int node=IGRAPH_VIT_GET(vit);
long int size=1;
added[node]=i+1;
igraph_dqueue_clear(&q);
if (order > 0) {
igraph_dqueue_push(&q, node);
igraph_dqueue_push(&q, 0);
}
while (!igraph_dqueue_empty(&q)) {
long int actnode=igraph_dqueue_pop(&q);
long int actdist=igraph_dqueue_pop(&q);
long int n;
igraph_neighbors(graph, &neis, actnode, mode);
n=igraph_vector_size(&neis);
if (actdist<order-1) {
/* we add them to the q */
for (j=0; j<n; j++) {
long int nei=VECTOR(neis)[j];
if (added[nei] != i+1) {
added[nei]=i+1;
IGRAPH_CHECK(igraph_dqueue_push(&q, nei));
IGRAPH_CHECK(igraph_dqueue_push(&q, actdist+1));
size++;
}
}
} else {
/* we just count them, but don't add them */
for (j=0; j<n; j++) {
long int nei=VECTOR(neis)[j];
if (added[nei] != i+1) {
added[nei]=i+1;
size++;
}
}
}
} /* while q not empty */
VECTOR(*res)[i]=size;
} /* for VIT, i */
igraph_vector_destroy(&neis);
igraph_vit_destroy(&vit);
igraph_dqueue_destroy(&q);
Free(added);
IGRAPH_FINALLY_CLEAN(4);
return 0;
}
/**
* \function igraph_neighborhood
* Calculate the neighborhood of vertices
*
* The neighborhood of a given order of a vertex includes all vertices
* which are closer to the vertex than the order. Ie. order 0 is
* always the vertex itself, order 1 is the vertex plus its immediate
* neighbors, order 2 is order 1 plus the immediate neighbors of the
* vertices in order 1, etc.
*
* </para><para> This function calculates the vertices within the
* neighborhood of the specified vertices.
* \param graph The input graph.
* \param res An initialized pointer vector. Note that the objects
* (pointers) in the vector will \em not be freed, but the pointer
* vector will be resized as needed. The result of the calculation
* will be stored here in \c vector_t objects.
* \param vids The vertices for which the calculation is performed.
* \param order Integer giving the order of the neighborhood.
* \param mode Specifies how to use the direction of the edges if a
* directed graph is analyzed. For \c IGRAPH_OUT only the outgoing
* edges are followed, so all vertices reachable from the source
* vertex in at most \c order steps are included. For \c IGRAPH_IN
* all vertices from which the source vertex is reachable in at most
* \c order steps are included. \c IGRAPH_ALL ignores the direction
* of the edges. This argument is ignored for undirected graphs.
* \return Error code.
*
* \sa \ref igraph_neighborhood_size() to calculate the size of the
* neighborhood, \ref igraph_neighborhood_graphs() for creating
* graphs from the neighborhoods.
*
* Time complexity: O(n*d*o), n is the number of vertices for which
* the calculation is performed, d is the average degree, o is the
* order.
*/
int igraph_neighborhood(const igraph_t *graph, igraph_vector_ptr_t *res,
igraph_vs_t vids, igraph_integer_t order,
igraph_neimode_t mode) {
long int no_of_nodes=igraph_vcount(graph);
igraph_dqueue_t q;
igraph_vit_t vit;
long int i, j;
long int *added;
igraph_vector_t neis;
igraph_vector_t tmp;
igraph_vector_t *newv;
if (order < 0) {
IGRAPH_ERROR("Negative order in neighborhood size", IGRAPH_EINVAL);
}
added=Calloc(no_of_nodes, long int);
if (added==0) {
IGRAPH_ERROR("Cannot calculate neighborhood size", IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(igraph_free, added);
IGRAPH_DQUEUE_INIT_FINALLY(&q, 100);
IGRAPH_CHECK(igraph_vit_create(graph, vids, &vit));
IGRAPH_FINALLY(igraph_vit_destroy, &vit);
IGRAPH_VECTOR_INIT_FINALLY(&neis, 0);
IGRAPH_VECTOR_INIT_FINALLY(&tmp, 0);
IGRAPH_CHECK(igraph_vector_ptr_resize(res, IGRAPH_VIT_SIZE(vit)));
for (i=0; !IGRAPH_VIT_END(vit); IGRAPH_VIT_NEXT(vit), i++) {
long int node=IGRAPH_VIT_GET(vit);
added[node]=i+1;
igraph_vector_clear(&tmp);
IGRAPH_CHECK(igraph_vector_push_back(&tmp, node));
if (order > 0) {
igraph_dqueue_push(&q, node);
igraph_dqueue_push(&q, 0);
}
while (!igraph_dqueue_empty(&q)) {
long int actnode=igraph_dqueue_pop(&q);
long int actdist=igraph_dqueue_pop(&q);
long int n;
igraph_neighbors(graph, &neis, actnode, mode);
n=igraph_vector_size(&neis);
if (actdist<order-1) {
/* we add them to the q */
for (j=0; j<n; j++) {
long int nei=VECTOR(neis)[j];
if (added[nei] != i+1) {
added[nei]=i+1;
IGRAPH_CHECK(igraph_dqueue_push(&q, nei));
IGRAPH_CHECK(igraph_dqueue_push(&q, actdist+1));
IGRAPH_CHECK(igraph_vector_push_back(&tmp, nei));
}
}
} else {
/* we just count them but don't add them to q */
for (j=0; j<n; j++) {
long int nei=VECTOR(neis)[j];
if (added[nei] != i+1) {
added[nei]=i+1;
IGRAPH_CHECK(igraph_vector_push_back(&tmp, nei));
}
}
}
} /* while q not empty */
newv=Calloc(1, igraph_vector_t);
if (newv==0) {
IGRAPH_ERROR("Cannot calculate neighborhood", IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(igraph_free, newv);
IGRAPH_CHECK(igraph_vector_copy(newv, &tmp));
VECTOR(*res)[i]=newv;
IGRAPH_FINALLY_CLEAN(1);
}
igraph_vector_destroy(&tmp);
igraph_vector_destroy(&neis);
igraph_vit_destroy(&vit);
igraph_dqueue_destroy(&q);
Free(added);
IGRAPH_FINALLY_CLEAN(5);
return 0;
}
/**
* \function igraph_neighborhood_graphs
* Create graphs from the neighborhood(s) of some vertex/vertices
*
* The neighborhood of a given order of a vertex includes all vertices
* which are closer to the vertex than the order. Ie. order 0 is
* always the vertex itself, order 1 is the vertex plus its immediate
* neighbors, order 2 is order 1 plus the immediate neighbors of the
* vertices in order 1, etc.
*
* </para><para> This function finds every vertex in the neighborhood
* of a given parameter vertex and creates a graph from these
* vertices.
*
* </para><para> The first version of this function was written by
* Vincent Matossian, thanks Vincent.
* \param graph The input graph.
* \param res Pointer to a pointer vector, the result will be stored
* here, ie. \c res will contain pointers to \c igraph_t
* objects. It will be resized if needed but note that the
* objects in the pointer vector will not be freed.
* \param vids The vertices for which the calculation is performed.
* \param order Integer giving the order of the neighborhood.
* \param mode Specifies how to use the direction of the edges if a
* directed graph is analyzed. For \c IGRAPH_OUT only the outgoing
* edges are followed, so all vertices reachable from the source
* vertex in at most \c order steps are counted. For \c IGRAPH_IN
* all vertices from which the source vertex is reachable in at most
* \c order steps are counted. \c IGRAPH_ALL ignores the direction
* of the edges. This argument is ignored for undirected graphs.
* \return Error code.
*
* \sa \ref igraph_neighborhood_size() for calculating the neighborhood
* sizes only, \ref igraph_neighborhood() for calculating the
* neighborhoods (but not creating graphs).
*
* Time complexity: O(n*(|V|+|E|)), where n is the number vertices for
* which the calculation is performed, |V| and |E| are the number of
* vertices and edges in the original input graph.
*/
int igraph_neighborhood_graphs(const igraph_t *graph, igraph_vector_ptr_t *res,
igraph_vs_t vids, igraph_integer_t order,
igraph_neimode_t mode) {
long int no_of_nodes=igraph_vcount(graph);
igraph_dqueue_t q;
igraph_vit_t vit;
long int i, j;
long int *added;
igraph_vector_t neis;
igraph_vector_t tmp;
igraph_t *newg;
if (order < 0) {
IGRAPH_ERROR("Negative order in neighborhood size", IGRAPH_EINVAL);
}
added=Calloc(no_of_nodes, long int);
if (added==0) {
IGRAPH_ERROR("Cannot calculate neighborhood size", IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(igraph_free, added);
IGRAPH_DQUEUE_INIT_FINALLY(&q, 100);
IGRAPH_CHECK(igraph_vit_create(graph, vids, &vit));
IGRAPH_FINALLY(igraph_vit_destroy, &vit);
IGRAPH_VECTOR_INIT_FINALLY(&neis, 0);
IGRAPH_VECTOR_INIT_FINALLY(&tmp, 0);
IGRAPH_CHECK(igraph_vector_ptr_resize(res, IGRAPH_VIT_SIZE(vit)));
for (i=0; !IGRAPH_VIT_END(vit); IGRAPH_VIT_NEXT(vit), i++) {
long int node=IGRAPH_VIT_GET(vit);
added[node]=i+1;
igraph_vector_clear(&tmp);
IGRAPH_CHECK(igraph_vector_push_back(&tmp, node));
if (order > 0) {
igraph_dqueue_push(&q, node);
igraph_dqueue_push(&q, 0);
}
while (!igraph_dqueue_empty(&q)) {
long int actnode=igraph_dqueue_pop(&q);
long int actdist=igraph_dqueue_pop(&q);
long int n;
igraph_neighbors(graph, &neis, actnode, mode);
n=igraph_vector_size(&neis);
if (actdist<order-1) {
/* we add them to the q */
for (j=0; j<n; j++) {
long int nei=VECTOR(neis)[j];
if (added[nei] != i+1) {
added[nei]=i+1;
IGRAPH_CHECK(igraph_dqueue_push(&q, nei));
IGRAPH_CHECK(igraph_dqueue_push(&q, actdist+1));
IGRAPH_CHECK(igraph_vector_push_back(&tmp, nei));
}
}
} else {
/* we just count them but don't add them to q */
for (j=0; j<n; j++) {
long int nei=VECTOR(neis)[j];
if (added[nei] != i+1) {
added[nei]=i+1;
IGRAPH_CHECK(igraph_vector_push_back(&tmp, nei));
}
}
}
} /* while q not empty */
newg=Calloc(1, igraph_t);
if (newg==0) {
IGRAPH_ERROR("Cannot create neighborhood graph", IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(igraph_free, newg);
if (igraph_vector_size(&tmp) < no_of_nodes) {
IGRAPH_CHECK(igraph_subgraph(graph, newg, igraph_vss_vector(&tmp)));
} else {
IGRAPH_CHECK(igraph_copy(newg, graph));
}
VECTOR(*res)[i]=newg;
IGRAPH_FINALLY_CLEAN(1);
}
igraph_vector_destroy(&tmp);
igraph_vector_destroy(&neis);
igraph_vit_destroy(&vit);
igraph_dqueue_destroy(&q);
Free(added);
IGRAPH_FINALLY_CLEAN(5);
return 0;
}
/**
* \function igraph_topological_sorting
* Calculate a possible topological sorting of the graph
*
* </para><para>
* A topological sorting of a directed acyclic graph is a linear ordering
* of its nodes where each node comes before all nodes to which it has
* edges. Every DAG has at least one topological sort, and may have many.
* This function returns a possible topological sort among them. If the
* graph is not acyclic (it has at least one cycle), a partial topological
* sort is returned and a warning is issued.
*
* \param graph The input graph.
* \param res Pointer to a vector, the result will be stored here.
* It will be resized if needed.
* \param mode Specifies how to use the direction of the edges.
* For \c IGRAPH_OUT, the sorting order ensures that each node comes
* before all nodes to which it has edges, so nodes with no incoming
* edges go first. For \c IGRAPH_IN, it is quite the opposite: each
* node comes before all nodes from which it receives edges. Nodes
* with no outgoing edges go first.
* \return Error code.
*
* Time complexity: O(|V|+|E|), where |V| and |E| are the number of
* vertices and edges in the original input graph.
*/
int igraph_topological_sorting(const igraph_t* graph, igraph_vector_t *res,
igraph_neimode_t mode) {
long int no_of_nodes=igraph_vcount(graph);
igraph_vector_t degrees, neis;
igraph_dqueue_t sources;
igraph_neimode_t deg_mode;
long int node, i, j;
if (mode == IGRAPH_ALL || !igraph_is_directed(graph)) {
IGRAPH_ERROR("topological sorting does not make sense for undirected graphs", IGRAPH_EINVAL);
} else if (mode == IGRAPH_OUT) {
deg_mode = IGRAPH_IN;
} else if (mode == IGRAPH_IN) {
deg_mode = IGRAPH_OUT;
} else {
IGRAPH_ERROR("invalid mode", IGRAPH_EINVAL);
}
IGRAPH_VECTOR_INIT_FINALLY(°rees, no_of_nodes);
IGRAPH_VECTOR_INIT_FINALLY(&neis, 0);
IGRAPH_CHECK(igraph_dqueue_init(&sources, 0));
IGRAPH_FINALLY(igraph_dqueue_destroy, &sources);
IGRAPH_CHECK(igraph_degree(graph, °rees, igraph_vss_all(), deg_mode, 0));
igraph_vector_clear(res);
/* Do we have nodes with no incoming vertices? */
for (i=0; i<no_of_nodes; i++) {
if (VECTOR(degrees)[i] == 0)
IGRAPH_CHECK(igraph_dqueue_push(&sources, i));
}
/* Take all nodes with no incoming vertices and remove them */
while (!igraph_dqueue_empty(&sources)) {
node=(long)igraph_dqueue_pop(&sources);
/* Add the node to the result vector */
igraph_vector_push_back(res, node);
/* Exclude the node from further source searches */
VECTOR(degrees)[node]=-1;
/* Get the neighbors and decrease their degrees by one */
IGRAPH_CHECK(igraph_neighbors(graph, &neis, node, mode));
j=igraph_vector_size(&neis);
for (i=0; i<j; i++) {
VECTOR(degrees)[(long)VECTOR(neis)[i]]--;
if (VECTOR(degrees)[(long)VECTOR(neis)[i]] == 0)
IGRAPH_CHECK(igraph_dqueue_push(&sources, VECTOR(neis)[i]));
}
}
if (igraph_vector_size(res)<no_of_nodes)
IGRAPH_WARNING("graph contains a cycle, partial result is returned");
igraph_vector_destroy(°rees);
igraph_vector_destroy(&neis);
igraph_dqueue_destroy(&sources);
IGRAPH_FINALLY_CLEAN(3);
return 0;
}
int igraph_is_loop(const igraph_t *graph, igraph_vector_t *res, igraph_es_t es) {
igraph_eit_t eit;
long int i;
IGRAPH_CHECK(igraph_eit_create(graph, es, &eit));
IGRAPH_FINALLY(igraph_eit_destroy, &eit);
IGRAPH_CHECK(igraph_vector_resize(res, IGRAPH_EIT_SIZE(eit)));
for (i=0; !IGRAPH_EIT_END(eit); i++, IGRAPH_EIT_NEXT(eit)) {
long int e=IGRAPH_EIT_GET(eit);
VECTOR(*res)[i] = (IGRAPH_FROM(graph, e)==IGRAPH_TO(graph, e)) ? 1 : 0;
}
igraph_eit_destroy(&eit);
IGRAPH_FINALLY_CLEAN(1);
return 0;
}
int igraph_is_multiple(const igraph_t *graph, igraph_vector_t *res, igraph_es_t es) {
igraph_eit_t eit;
long int i;
igraph_i_lazy_adjedgelist_t adjlist;
IGRAPH_CHECK(igraph_eit_create(graph, es, &eit));
IGRAPH_FINALLY(igraph_eit_destroy, &eit);
IGRAPH_CHECK(igraph_i_lazy_adjedgelist_init(graph, &adjlist, IGRAPH_OUT));
IGRAPH_FINALLY(igraph_i_lazy_adjedgelist_destroy, &adjlist);
IGRAPH_CHECK(igraph_vector_resize(res, IGRAPH_EIT_SIZE(eit)));
for (i=0; !IGRAPH_EIT_END(eit); i++, IGRAPH_EIT_NEXT(eit)) {
long int e=IGRAPH_EIT_GET(eit);
long int from=IGRAPH_FROM(graph, e);
long int to=IGRAPH_TO(graph, e);
igraph_vector_t *neis=igraph_i_lazy_adjedgelist_get(&adjlist, from);
long int j, n=igraph_vector_size(neis);
VECTOR(*res)[i]=0;
for (j=0; j<n; j++) {
long int e2=VECTOR(*neis)[j];
long int to2=IGRAPH_OTHER(graph,e2,from);
if (to2==to && e2<e) {
VECTOR(*res)[i]=1;
}
}
}
igraph_i_lazy_adjedgelist_destroy(&adjlist);
igraph_eit_destroy(&eit);
IGRAPH_FINALLY_CLEAN(2);
return 0;
}
/**
* \function igraph_girth
* \brief The girth of a graph is the length of the shortest circle in it.
*
* </para><para>
* The current implementation works for undirected graphs only,
* directed graphs are treated as undirected graphs. Loop edges and
* multiple edges are ignored.
* </para><para>
* If the graph is a forest (ie. acyclic), then zero is returned.
* </para><para>
* This implementation is based on Alon Itai and Michael Rodeh:
* Finding a minimum circuit in a graph
* \emb Proceedings of the ninth annual ACM symposium on Theory of
* computing \eme, 1-10, 1977. The first implementation of this
* function was done by Keith Briggs, thanks Keith.
* \param graph The input graph.
* \param girth Pointer to an integer, if not \c NULL then the result
* will be stored here.
* \param circle Pointer to an initialized vector, the vertex ids in
* the shortest circle will be stored here. If \c NULL then it is
* ignored.
* \return Error code.
*
* Time complexity: O((|V|+|E|)^2), |V| is the number of vertices, |E|
* is the number of edges in the general case. If the graph has no
* circles at all then the function needs O(|V|+|E|) time to realize
* this and then it stops.
*/
int igraph_girth(const igraph_t *graph, igraph_integer_t *girth,
igraph_vector_t *circle) {
long int no_of_nodes=igraph_vcount(graph);
igraph_dqueue_t q;
igraph_i_lazy_adjlist_t adjlist;
long int mincirc=LONG_MAX, minvertex=0;
long int node;
igraph_bool_t triangle=0;
igraph_vector_t *neis;
igraph_vector_long_t level;
long int stoplevel=no_of_nodes+1;
igraph_bool_t anycircle=0;
long int t1=0, t2=0;
IGRAPH_CHECK(igraph_i_lazy_adjlist_init(graph, &adjlist, IGRAPH_ALL,
IGRAPH_I_SIMPLIFY));
IGRAPH_FINALLY(igraph_i_lazy_adjlist_destroy, &adjlist);
IGRAPH_DQUEUE_INIT_FINALLY(&q, 100);
IGRAPH_CHECK(igraph_vector_long_init(&level, no_of_nodes));
IGRAPH_FINALLY(igraph_vector_long_destroy, &level);
for (node=0; !triangle && node<no_of_nodes; node++) {
/* Are there circles in this graph at all? */
if (node==1 && anycircle==0) {
igraph_bool_t conn;
IGRAPH_CHECK(igraph_is_connected(graph, &conn, IGRAPH_WEAK));
if (conn) {
/* No, there are none */
break;
}
}
anycircle=0;
igraph_dqueue_clear(&q);
igraph_vector_long_null(&level);
IGRAPH_CHECK(igraph_dqueue_push(&q, node));
VECTOR(level)[node]=1;
IGRAPH_ALLOW_INTERRUPTION();
while (!igraph_dqueue_empty(&q)) {
long int actnode=igraph_dqueue_pop(&q);
long int actlevel=VECTOR(level)[actnode];
long int i, n;
if (actlevel>=stoplevel) { break; }
neis=igraph_i_lazy_adjlist_get(&adjlist, actnode);
n=igraph_vector_size(neis);
for (i=0; i<n; i++) {
long int nei=VECTOR(*neis)[i];
long int neilevel=VECTOR(level)[nei];
if (neilevel != 0) {
if (neilevel==actlevel-1) {
continue;
} else {
/* found circle */
stoplevel=neilevel;
anycircle=1;
if (actlevel<mincirc) {
/* Is it a minimum circle? */
mincirc=actlevel+neilevel-1;
minvertex=node;
t1=actnode; t2=nei;
if (neilevel==2) {
/* Is it a triangle? */
triangle=1;
}
}
if (neilevel==actlevel) {
break;
}
}
} else {
igraph_dqueue_push(&q, nei);
VECTOR(level)[nei]=actlevel+1;
}
}
} /* while q !empty */
} /* node */
if (girth) {
if (mincirc==LONG_MAX) {
*girth=mincirc=0;
} else {
*girth=mincirc;
}
}
/* Store the actual circle, if needed */
if (circle) {
IGRAPH_CHECK(igraph_vector_resize(circle, mincirc));
if (mincirc != 0) {
long int i, n, idx=0;
igraph_dqueue_clear(&q);
igraph_vector_long_null(&level); /* used for father pointers */
#define FATHER(x) (VECTOR(level)[(x)])
IGRAPH_CHECK(igraph_dqueue_push(&q, minvertex));
FATHER(minvertex)=minvertex;
while (FATHER(t1)==0 || FATHER(t2)==0) {
long int actnode=igraph_dqueue_pop(&q);
neis=igraph_i_lazy_adjlist_get(&adjlist, actnode);
n=igraph_vector_size(neis);
for (i=0; i<n; i++) {
long int nei=VECTOR(*neis)[i];
if (FATHER(nei) == 0) {
FATHER(nei)=actnode+1;
igraph_dqueue_push(&q, nei);
}
}
} /* while q !empty */
/* Ok, now use FATHER to create the path */
while (t1!=minvertex) {
VECTOR(*circle)[idx++]=t1;
t1=FATHER(t1)-1;
}
VECTOR(*circle)[idx]=minvertex;
idx=mincirc-1;
while (t2!=minvertex) {
VECTOR(*circle)[idx--]=t2;
t2=FATHER(t2)-1;
}
} /* anycircle */
} /* circle */
#undef FATHER
igraph_vector_long_destroy(&level);
igraph_dqueue_destroy(&q);
igraph_i_lazy_adjlist_destroy(&adjlist);
IGRAPH_FINALLY_CLEAN(3);
return 0;
}
syntax highlighted by Code2HTML, v. 0.9.1