haskell-igraph-0.8.0: igraph/src/structural_properties.c
/* -*- mode: C -*- */
/* vim:set ts=4 sw=4 sts=4 et: */
/*
IGraph library.
Copyright (C) 2005-2012 Gabor Csardi <csardi.gabor@gmail.com>
334 Harvard street, Cambridge, MA 02139 USA
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
02110-1301 USA
*/
#include "igraph_structural.h"
#include "igraph_transitivity.h"
#include "igraph_paths.h"
#include "igraph_math.h"
#include "igraph_memory.h"
#include "igraph_random.h"
#include "igraph_adjlist.h"
#include "igraph_interface.h"
#include "igraph_progress.h"
#include "igraph_interrupt_internal.h"
#include "igraph_centrality.h"
#include "igraph_components.h"
#include "igraph_constructors.h"
#include "igraph_conversion.h"
#include "igraph_types_internal.h"
#include "igraph_dqueue.h"
#include "igraph_attributes.h"
#include "igraph_neighborhood.h"
#include "igraph_topology.h"
#include "igraph_qsort.h"
#include "config.h"
#include "structural_properties_internal.h"
#include <assert.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.
*
* \example examples/simple/igraph_diameter.c
*/
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_int_t *neis;
igraph_neimode_t dirmode;
igraph_adjlist_t allneis;
if (directed) {
dirmode = IGRAPH_OUT;
} else {
dirmode = IGRAPH_ALL;
}
already_added = igraph_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_CHECK(igraph_adjlist_init(graph, &allneis, dirmode));
IGRAPH_FINALLY(igraph_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 = (long int) igraph_dqueue_pop(&q);
long int actdist = (long int) igraph_dqueue_pop(&q);
if (actdist > res) {
res = actdist;
from = i;
to = actnode;
}
neis = igraph_adjlist_get(&allneis, actnode);
n = igraph_vector_int_size(neis);
for (j = 0; j < n; j++) {
long int neighbor = (long int) 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 = (igraph_integer_t) res;
}
if (pfrom != 0) {
*pfrom = (igraph_integer_t) from;
}
if (pto != 0) {
*pto = (igraph_integer_t) 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, 0,
(igraph_integer_t) from,
igraph_vss_1((igraph_integer_t)to),
dirmode, 0, 0));
igraph_vector_ptr_destroy(&tmpptr);
IGRAPH_FINALLY_CLEAN(1);
}
}
/* clean */
igraph_Free(already_added);
igraph_dqueue_destroy(&q);
igraph_adjlist_destroy(&allneis);
IGRAPH_FINALLY_CLEAN(3);
return 0;
}
/**
* \ingroup structural
* \function igraph_average_path_length
* \brief Calculates the average shortest path length between all vertex pairs.
*
* \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, only those vertex pairs will be included in the calculation
* between which there is a path. If \c FALSE, the number of vertices is
* used as the distance between vertices unreachable from each other.
* 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.
*
* \example examples/simple/igraph_average_path_length.c
*/
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_int_t *neis;
igraph_neimode_t dirmode;
igraph_adjlist_t allneis;
*res = 0;
if (directed) {
dirmode = IGRAPH_OUT;
} else {
dirmode = IGRAPH_ALL;
}
already_added = igraph_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_adjlist_init(graph, &allneis, dirmode);
IGRAPH_FINALLY(igraph_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 = (long int) igraph_dqueue_pop(&q);
long int actdist = (long int) igraph_dqueue_pop(&q);
neis = igraph_adjlist_get(&allneis, actnode);
n = igraph_vector_int_size(neis);
for (j = 0; j < n; j++) {
long int neighbor = (long int) 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 */
if (normfact > 0) {
*res /= normfact;
} else {
*res = IGRAPH_NAN;
}
/* clean */
igraph_Free(already_added);
igraph_dqueue_destroy(&q);
igraph_adjlist_destroy(&allneis);
IGRAPH_FINALLY_CLEAN(3);
return 0;
}
/**
* \function igraph_path_length_hist
* Create a histogram of all shortest path lengths.
*
* This function calculates a histogram, by calculating the
* shortest path length between each pair of vertices. For directed
* graphs both directions might be considered and then every pair of vertices
* appears twice in the histogram.
* \param graph The input graph.
* \param res Pointer to an initialized vector, the result is stored
* here. The first (i.e. zeroth) element contains the number of
* shortest paths of length 1, etc. The supplied vector is resized
* as needed.
* \param unconnected Pointer to a real number, the number of
* pairs for which the second vertex is not reachable from the
* first is stored here.
* \param directed Whether to consider directed paths in a directed
* graph (if not zero). This argument is ignored for undirected
* graphs.
* \return Error code.
*
* Time complexity: O(|V||E|), the number of vertices times the number
* of edges.
*
* \sa \ref igraph_average_path_length() and \ref igraph_shortest_paths()
*/
int igraph_path_length_hist(const igraph_t *graph, igraph_vector_t *res,
igraph_real_t *unconnected, igraph_bool_t directed) {
long int no_of_nodes = igraph_vcount(graph);
long int i, j, n;
igraph_vector_long_t already_added;
long int nodes_reached;
igraph_dqueue_t q = IGRAPH_DQUEUE_NULL;
igraph_vector_int_t *neis;
igraph_neimode_t dirmode;
igraph_adjlist_t allneis;
igraph_real_t unconn = 0;
long int ressize;
if (directed) {
dirmode = IGRAPH_OUT;
} else {
dirmode = IGRAPH_ALL;
}
IGRAPH_CHECK(igraph_vector_long_init(&already_added, no_of_nodes));
IGRAPH_FINALLY(igraph_vector_long_destroy, &already_added);
IGRAPH_DQUEUE_INIT_FINALLY(&q, 100);
IGRAPH_CHECK(igraph_adjlist_init(graph, &allneis, dirmode));
IGRAPH_FINALLY(igraph_adjlist_destroy, &allneis);
IGRAPH_CHECK(igraph_vector_resize(res, 0));
ressize = 0;
for (i = 0; i < no_of_nodes; i++) {
nodes_reached = 1; /* itself */
IGRAPH_CHECK(igraph_dqueue_push(&q, i));
IGRAPH_CHECK(igraph_dqueue_push(&q, 0));
VECTOR(already_added)[i] = i + 1;
IGRAPH_PROGRESS("Path-hist: ", 100.0 * i / no_of_nodes, NULL);
IGRAPH_ALLOW_INTERRUPTION();
while (!igraph_dqueue_empty(&q)) {
long int actnode = (long int) igraph_dqueue_pop(&q);
long int actdist = (long int) igraph_dqueue_pop(&q);
neis = igraph_adjlist_get(&allneis, actnode);
n = igraph_vector_int_size(neis);
for (j = 0; j < n; j++) {
long int neighbor = (long int) VECTOR(*neis)[j];
if (VECTOR(already_added)[neighbor] == i + 1) {
continue;
}
VECTOR(already_added)[neighbor] = i + 1;
nodes_reached++;
if (actdist + 1 > ressize) {
IGRAPH_CHECK(igraph_vector_resize(res, actdist + 1));
for (; ressize < actdist + 1; ressize++) {
VECTOR(*res)[ressize] = 0;
}
}
VECTOR(*res)[actdist] += 1;
IGRAPH_CHECK(igraph_dqueue_push(&q, neighbor));
IGRAPH_CHECK(igraph_dqueue_push(&q, actdist + 1));
}
} /* while !igraph_dqueue_empty */
unconn += (no_of_nodes - nodes_reached);
} /* for i<no_of_nodes */
IGRAPH_PROGRESS("Path-hist: ", 100.0, NULL);
/* count every pair only once for an undirected graph */
if (!directed || !igraph_is_directed(graph)) {
for (i = 0; i < ressize; i++) {
VECTOR(*res)[i] /= 2;
}
unconn /= 2;
}
igraph_vector_long_destroy(&already_added);
igraph_dqueue_destroy(&q);
igraph_adjlist_destroy(&allneis);
IGRAPH_FINALLY_CLEAN(3);
if (unconnected) {
*unconnected = unconn;
}
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. A pointer to an
* initialized matrix, to be more precise. The matrix will be
* resized if needed. It will have 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 \c to argument. One row of the matrix shows the
* distances from/to a given vertex to the ones in \c to.
* For the unreachable vertices IGRAPH_INFINITY is returned.
* \param from Vector of the vertex ids for which the path length
* calculations are done.
* \param to Vector of the vertex ids to which the path length
* calculations are done. It is not allowed to have duplicated
* vertex ids here.
* \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(|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,
* \ref igraph_shortest_paths_dijkstra() for the weighted version.
*/
int igraph_shortest_paths(const igraph_t *graph, igraph_matrix_t *res,
const igraph_vs_t from, const igraph_vs_t to,
igraph_neimode_t mode) {
long int no_of_nodes = igraph_vcount(graph);
long int no_of_from, no_of_to;
long int *already_counted;
igraph_adjlist_t adjlist;
igraph_dqueue_t q = IGRAPH_DQUEUE_NULL;
igraph_vector_int_t *neis;
igraph_bool_t all_to;
long int i, j;
igraph_vit_t fromvit, tovit;
igraph_real_t my_infinity = IGRAPH_INFINITY;
igraph_vector_t indexv;
if (mode != IGRAPH_OUT && mode != IGRAPH_IN &&
mode != IGRAPH_ALL) {
IGRAPH_ERROR("Invalid mode argument", IGRAPH_EINVMODE);
}
IGRAPH_CHECK(igraph_vit_create(graph, from, &fromvit));
IGRAPH_FINALLY(igraph_vit_destroy, &fromvit);
no_of_from = IGRAPH_VIT_SIZE(fromvit);
IGRAPH_CHECK(igraph_adjlist_init(graph, &adjlist, mode));
IGRAPH_FINALLY(igraph_adjlist_destroy, &adjlist);
already_counted = igraph_Calloc(no_of_nodes, long int);
if (already_counted == 0) {
IGRAPH_ERROR("shortest paths failed", IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(free, already_counted);
IGRAPH_DQUEUE_INIT_FINALLY(&q, 100);
if ( (all_to = igraph_vs_is_all(&to)) ) {
no_of_to = no_of_nodes;
} else {
IGRAPH_VECTOR_INIT_FINALLY(&indexv, no_of_nodes);
IGRAPH_CHECK(igraph_vit_create(graph, to, &tovit));
IGRAPH_FINALLY(igraph_vit_destroy, &tovit);
no_of_to = IGRAPH_VIT_SIZE(tovit);
for (i = 0; !IGRAPH_VIT_END(tovit); IGRAPH_VIT_NEXT(tovit)) {
long int v = IGRAPH_VIT_GET(tovit);
if (VECTOR(indexv)[v]) {
IGRAPH_ERROR("Duplicate vertices in `to', this is not allowed",
IGRAPH_EINVAL);
}
VECTOR(indexv)[v] = ++i;
}
}
IGRAPH_CHECK(igraph_matrix_resize(res, no_of_from, no_of_to));
igraph_matrix_fill(res, my_infinity);
for (IGRAPH_VIT_RESET(fromvit), i = 0;
!IGRAPH_VIT_END(fromvit);
IGRAPH_VIT_NEXT(fromvit), i++) {
long int reached = 0;
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 = (long int) igraph_dqueue_pop(&q);
long int actdist = (long int) igraph_dqueue_pop(&q);
if (all_to) {
MATRIX(*res, i, act) = actdist;
} else {
if (VECTOR(indexv)[act]) {
MATRIX(*res, i, (long int)(VECTOR(indexv)[act] - 1)) = actdist;
reached++;
if (reached == no_of_to) {
igraph_dqueue_clear(&q);
break;
}
}
}
neis = igraph_adjlist_get(&adjlist, act);
for (j = 0; j < igraph_vector_int_size(neis); j++) {
long int neighbor = (long int) VECTOR(*neis)[j];
if (already_counted[neighbor] == i + 1) {
continue;
}
already_counted[neighbor] = i + 1;
IGRAPH_CHECK(igraph_dqueue_push(&q, neighbor));
IGRAPH_CHECK(igraph_dqueue_push(&q, actdist + 1));
}
}
}
/* Clean */
if (!all_to) {
igraph_vit_destroy(&tovit);
igraph_vector_destroy(&indexv);
IGRAPH_FINALLY_CLEAN(2);
}
igraph_Free(already_counted);
igraph_dqueue_destroy(&q);
igraph_vit_destroy(&fromvit);
igraph_adjlist_destroy(&adjlist);
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 vertices The result, the ids of the vertices along the paths.
* 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. Supply a null pointer here if you don't need
* these vectors.
* \param edges The result, the ids of the edges along the paths.
* 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. Supply a null pointer here if you don't need
* these vectors.
* \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 used 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
* \param predecessors A pointer to an initialized igraph vector or null.
* If not null, a vector containing the predecessor of each vertex in
* the single source shortest path tree is returned here. The
* predecessor of vertex i in the tree is the vertex from which vertex i
* was reached. The predecessor of the start vertex (in the \c from
* argument) is itself by definition. If the predecessor is -1, it means
* that the given vertex was not reached from the source during the
* search. Note that the search terminates if all the vertices in
* \c to are reached.
* \param inbound_edges A pointer to an initialized igraph vector or null.
* If not null, a vector containing the inbound edge of each vertex in
* the single source shortest path tree is returned here. The
* inbound edge of vertex i in the tree is the edge via which vertex i
* was reached. The start vertex and vertices that were not reached
* during the search will have -1 in the corresponding entry of the
* vector. Note that the search terminates if all the vertices in
* \c to are reached.
*
* \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.
*
* \example examples/simple/igraph_get_shortest_paths.c
*/
int igraph_get_shortest_paths(const igraph_t *graph,
igraph_vector_ptr_t *vertices,
igraph_vector_ptr_t *edges,
igraph_integer_t from, const igraph_vs_t to,
igraph_neimode_t mode,
igraph_vector_long_t *predecessors,
igraph_vector_long_t *inbound_edges) {
/* TODO: use inclist_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 i, 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 (vertices && IGRAPH_VIT_SIZE(vit) != igraph_vector_ptr_size(vertices)) {
IGRAPH_ERROR("Size of the `vertices' and the `to' should match", IGRAPH_EINVAL);
}
if (edges && IGRAPH_VIT_SIZE(vit) != igraph_vector_ptr_size(edges)) {
IGRAPH_ERROR("Size of the `edges' and the `to' should match", IGRAPH_EINVAL);
}
father = igraph_Calloc(no_of_nodes, long int);
if (father == 0) {
IGRAPH_ERROR("cannot get shortest paths", IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(igraph_free, father);
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 */
}
}
/* Meaning of father[i]:
*
* - If father[i] < 0, it means that vertex i has to be reached and has not
* been reached yet.
*
* - If father[i] = 0, it means that vertex i does not have to be reached and
* it has not been reached yet.
*
* - If father[i] = 1, it means that vertex i is the start vertex.
*
* - Otherwise, father[i] is the ID of the edge from which vertex i was
* reached plus 2.
*/
IGRAPH_CHECK(igraph_dqueue_push(&q, from + 1));
if (father[ (long int) from ] < 0) {
reached++;
}
father[ (long int)from ] = 1;
while (!igraph_dqueue_empty(&q) && reached < to_reach) {
long int act = (long int) igraph_dqueue_pop(&q) - 1;
IGRAPH_CHECK(igraph_incident(graph, &tmp, (igraph_integer_t) act, mode));
for (j = 0; j < igraph_vector_size(&tmp); j++) {
long int edge = (long int) VECTOR(tmp)[j];
long int neighbor = IGRAPH_OTHER(graph, edge, act);
if (father[neighbor] > 0) {
continue;
} else if (father[neighbor] < 0) {
reached++;
}
father[neighbor] = edge + 2;
IGRAPH_CHECK(igraph_dqueue_push(&q, neighbor + 1));
}
}
if (reached < to_reach) {
IGRAPH_WARNING("Couldn't reach some vertices");
}
/* Create `predecessors' if needed */
if (predecessors) {
IGRAPH_CHECK(igraph_vector_long_resize(predecessors, no_of_nodes));
for (i = 0; i < no_of_nodes; i++) {
if (father[i] <= 0) {
/* i was not reached */
VECTOR(*predecessors)[i] = -1;
} else if (father[i] == 1) {
/* i is the start vertex */
VECTOR(*predecessors)[i] = i;
} else {
/* i was reached via the edge with ID = father[i] - 2 */
VECTOR(*predecessors)[i] = IGRAPH_OTHER(graph, father[i] - 2, i);
}
}
}
/* Create `inbound_edges' if needed */
if (inbound_edges) {
IGRAPH_CHECK(igraph_vector_long_resize(inbound_edges, no_of_nodes));
for (i = 0; i < no_of_nodes; i++) {
if (father[i] <= 1) {
/* i was not reached or i is the start vertex */
VECTOR(*inbound_edges)[i] = -1;
} else {
/* i was reached via the edge with ID = father[i] - 2 */
VECTOR(*inbound_edges)[i] = father[i] - 2;
}
}
}
/* Create `vertices' and `edges' if needed */
if (vertices || edges) {
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 *vvec = 0, *evec = 0;
if (vertices) {
vvec = VECTOR(*vertices)[j];
igraph_vector_clear(vvec);
}
if (edges) {
evec = VECTOR(*edges)[j];
igraph_vector_clear(evec);
}
IGRAPH_ALLOW_INTERRUPTION();
if (father[node] > 0) {
long int act = node;
long int size = 0;
long int edge;
while (father[act] > 1) {
size++;
edge = father[act] - 2;
act = IGRAPH_OTHER(graph, edge, act);
}
if (vvec) {
IGRAPH_CHECK(igraph_vector_resize(vvec, size + 1));
VECTOR(*vvec)[size] = node;
}
if (evec) {
IGRAPH_CHECK(igraph_vector_resize(evec, size));
}
act = node;
while (father[act] > 1) {
size--;
edge = father[act] - 2;
act = IGRAPH_OTHER(graph, edge, act);
if (vvec) {
VECTOR(*vvec)[size] = act;
}
if (evec) {
VECTOR(*evec)[size] = edge;
}
}
}
}
}
/* Clean */
igraph_Free(father);
igraph_dqueue_destroy(&q);
igraph_vector_destroy(&tmp);
igraph_vit_destroy(&vit);
IGRAPH_FINALLY_CLEAN(4);
return 0;
}
/**
* \function igraph_get_shortest_path
* Shortest path from one vertex to another one.
*
* Calculates and returns a single unweighted shortest path from a
* given vertex to another one. If there are more than one shortest
* paths between the two vertices, then an arbitrary one is returned.
*
* </para><para>This function is a wrapper to \ref
* igraph_get_shortest_paths(), for the special case when only one
* target vertex is considered.
* \param graph The input graph, it can be directed or
* undirected. Directed paths are considered in directed
* graphs.
* \param vertices Pointer to an initialized vector or a null
* pointer. If not a null pointer, then the vertex ids along
* the path are stored here, including the source and target
* vertices.
* \param edges Pointer to an uninitialized vector or a null
* pointer. If not a null pointer, then the edge ids along the
* path are stored here.
* \param from The id of the source vertex.
* \param to The id of the target vertex.
* \param mode A constant specifying how edge directions are
* considered in directed graphs. Valid modes are:
* \c IGRAPH_OUT, follows edge directions;
* \c IGRAPH_IN, follows the opposite directions; and
* \c IGRAPH_ALL, ignores edge directions. This argument is
* ignored for undirected graphs.
* \return Error code.
*
* Time complexity: O(|V|+|E|), linear in the number of vertices and
* edges in the graph.
*
* \sa \ref igraph_get_shortest_paths() for the version with more target
* vertices.
*/
int igraph_get_shortest_path(const igraph_t *graph,
igraph_vector_t *vertices,
igraph_vector_t *edges,
igraph_integer_t from,
igraph_integer_t to,
igraph_neimode_t mode) {
igraph_vector_ptr_t vertices2, *vp = &vertices2;
igraph_vector_ptr_t edges2, *ep = &edges2;
if (vertices) {
IGRAPH_CHECK(igraph_vector_ptr_init(&vertices2, 1));
IGRAPH_FINALLY(igraph_vector_ptr_destroy, &vertices2);
VECTOR(vertices2)[0] = vertices;
} else {
vp = 0;
}
if (edges) {
IGRAPH_CHECK(igraph_vector_ptr_init(&edges2, 1));
IGRAPH_FINALLY(igraph_vector_ptr_destroy, &edges2);
VECTOR(edges2)[0] = edges;
} else {
ep = 0;
}
IGRAPH_CHECK(igraph_get_shortest_paths(graph, vp, ep, from,
igraph_vss_1(to), mode, 0, 0));
if (edges) {
igraph_vector_ptr_destroy(&edges2);
IGRAPH_FINALLY_CLEAN(1);
}
if (vertices) {
igraph_vector_ptr_destroy(&vertices2);
IGRAPH_FINALLY_CLEAN(1);
}
return 0;
}
void igraph_i_gasp_paths_destroy(igraph_vector_ptr_t *v);
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]);
igraph_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
* stored here for every vertex in the graph. Note that the values
* will be accurate only for those vertices that are in the target
* vertex sequence (see \p to), since the search terminates as soon
* as all the target vertices have been found.
* \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 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, const 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;
long int to_reach, reached = 0, maxdist = 0;
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);
/* paths will store the shortest paths during the search */
IGRAPH_CHECK(igraph_vector_ptr_init(&paths, 0));
IGRAPH_FINALLY(igraph_i_gasp_paths_destroy, &paths);
/* neis is a temporary vector holding the neighbors of the
* node being examined */
IGRAPH_VECTOR_INIT_FINALLY(&neis, 0);
/* ptrlist stores indices into the paths vector, in the order
* of how they were found. ptrhead is a second-level index that
* will be used to find paths that terminate in a given vertex */
IGRAPH_VECTOR_INIT_FINALLY(&ptrlist, 0);
/* ptrhead contains indices into ptrlist.
* ptrhead[i] = j means that element #j-1 in ptrlist contains
* the shortest path from the root to node i. ptrhead[i] = 0
* means that node i was not reached so far */
IGRAPH_VECTOR_INIT_FINALLY(&ptrhead, no_of_nodes);
/* geodist[i] == 0 if i was not reached yet and it is not in the
* target vertex sequence, or -1 if i was not reached yet and it
* is in the target vertex sequence. Otherwise it is
* one larger than the length of the shortest path from the
* source */
geodist = igraph_Calloc(no_of_nodes, long int);
if (geodist == 0) {
IGRAPH_ERROR("Cannot calculate shortest paths", IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(igraph_free, geodist);
/* dequeue to store the BFS queue -- odd elements are the vertex indices,
* even elements are the distances from the root */
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);
}
/* use geodist to count how many vertices we have to reach */
to_reach = IGRAPH_VIT_SIZE(vit);
for (IGRAPH_VIT_RESET(vit); !IGRAPH_VIT_END(vit); IGRAPH_VIT_NEXT(vit)) {
if (geodist[ (long int) IGRAPH_VIT_GET(vit) ] == 0) {
geodist[ (long int) IGRAPH_VIT_GET(vit) ] = -1;
} else {
to_reach--; /* this node was given multiple times */
}
}
if (geodist[ (long int) from ] < 0) {
reached++;
}
/* from -> from */
vptr = igraph_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 = (long int) igraph_dqueue_pop(&q);
long int actdist = (long int) igraph_dqueue_pop(&q);
IGRAPH_ALLOW_INTERRUPTION();
if (reached >= to_reach) {
/* all nodes were reached. Since we need all the shortest paths
* to all these nodes, we can stop the search only if the distance
* of the current node to the root is larger than the distance of
* any of the nodes we wanted to reach */
if (actdist > maxdist) {
/* safety check, maxdist should have been set when we reached the last node */
if (maxdist < 0) {
IGRAPH_ERROR("possible bug in igraph_get_all_shortest_paths, "
"maxdist is negative", IGRAPH_EINVAL);
}
break;
}
}
IGRAPH_CHECK(igraph_neighbors(graph, &neis, (igraph_integer_t) actnode,
mode));
n = igraph_vector_size(&neis);
for (j = 0; j < n; j++) {
long int neighbor = (long int) VECTOR(neis)[j];
long int fatherptr;
if (geodist[neighbor] > 0 &&
geodist[neighbor] - 1 < actdist + 1) {
/* this node was reached via a shorter path before */
continue;
}
/* yay, found another shortest path to neighbor */
if (nrgeo) {
/* the number of geodesics leading to neighbor must be
* increased by the number of geodesics leading to actnode */
VECTOR(*nrgeo)[neighbor] += VECTOR(*nrgeo)[actnode];
}
if (geodist[neighbor] <= 0) {
/* this node was not reached yet, push it into the queue */
IGRAPH_CHECK(igraph_dqueue_push(&q, neighbor));
IGRAPH_CHECK(igraph_dqueue_push(&q, actdist + 1));
if (geodist[neighbor] < 0) {
reached++;
}
if (reached == to_reach) {
maxdist = actdist;
}
}
geodist[neighbor] = actdist + 2;
/* copy all existing paths to the parent */
fatherptr = (long int) VECTOR(ptrhead)[actnode];
while (fatherptr != 0) {
/* allocate a new igraph_vector_t at the end of paths */
vptr = igraph_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_CHECK(igraph_vector_push_back(vptr, neighbor));
IGRAPH_CHECK(igraph_vector_push_back(&ptrlist,
VECTOR(ptrhead)[neighbor]));
VECTOR(ptrhead)[neighbor] = igraph_vector_size(&ptrlist);
fatherptr = (long int) VECTOR(ptrlist)[fatherptr - 1];
}
}
}
igraph_dqueue_destroy(&q);
IGRAPH_FINALLY_CLEAN(1);
/* mark the nodes for which we need the result */
memset(geodist, 0, sizeof(long int) * (size_t) no_of_nodes);
for (IGRAPH_VIT_RESET(vit); !IGRAPH_VIT_END(vit); IGRAPH_VIT_NEXT(vit)) {
geodist[ (long int) IGRAPH_VIT_GET(vit) ] = 1;
}
/* count the number of paths in the result */
n = 0;
for (i = 0; i < no_of_nodes; i++) {
long int fatherptr = (long int) VECTOR(ptrhead)[i];
if (geodist[i] > 0) {
while (fatherptr != 0) {
n++;
fatherptr = (long int) 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 = (long int) VECTOR(ptrhead)[i];
IGRAPH_ALLOW_INTERRUPTION();
/* do we need the paths leading to vertex i? */
if (geodist[i] > 0) {
/* yes, copy them to the result vector */
while (fatherptr != 0) {
VECTOR(*res)[j++] = VECTOR(paths)[fatherptr - 1];
fatherptr = (long int) VECTOR(ptrlist)[fatherptr - 1];
}
} else {
/* no, free them */
while (fatherptr != 0) {
igraph_vector_destroy(VECTOR(paths)[fatherptr - 1]);
igraph_Free(VECTOR(paths)[fatherptr - 1]);
fatherptr = (long int) VECTOR(ptrlist)[fatherptr - 1];
}
}
}
igraph_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 (!IGRAPH_FINITE(vertex) || 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 = igraph_Calloc(no_of_nodes, char);
if (already_added == 0) {
IGRAPH_ERROR("subcomponent failed", IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(free, already_added); /* TODO: hack */
igraph_vector_clear(res);
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 = (long int) igraph_dqueue_pop(&q);
IGRAPH_ALLOW_INTERRUPTION();
IGRAPH_CHECK(igraph_neighbors(graph, &tmp, (igraph_integer_t) actnode,
mode));
for (i = 0; i < igraph_vector_size(&tmp); i++) {
long int neighbor = (long int) 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);
igraph_Free(already_added);
IGRAPH_FINALLY_CLEAN(3);
return 0;
}
/**
* \ingroup structural
* \function igraph_pagerank_old
* \brief Calculates the Google PageRank for the specified vertices.
*
* </para><para>This is an old implementation,
* it is provided for compatibility with igraph versions earlier than
* 0.5. Please use the new implementation \ref igraph_pagerank() in
* new projects.
*
* </para><para>
* From version 0.7 this function is deprecated and its use gives a
* warning message.
*
* </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://infolab.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)
* \param old Boolean, whether to use the pre-igraph 0.5 way to
* calculate page rank. Not recommended for new applications,
* only included for compatibility. If this is non-zero then the damping
* factor is not divided by the number of vertices before adding it
* to the weighted page rank scores to calculate the
* new scores. I.e. the formula in the original PageRank paper
* is used. Furthermore, if this is non-zero then the PageRank
* vector is renormalized after each iteration.
* \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|) per iteration. A handful iterations
* should be enough. Note that if the old-style dumping is used then
* the iteration might not converge at all.
*
* \sa \ref igraph_pagerank() for the new implementation.
*/
int igraph_pagerank_old(const igraph_t *graph, igraph_vector_t *res,
const igraph_vs_t vids, igraph_bool_t directed,
igraph_integer_t niter, igraph_real_t eps,
igraph_real_t damping, igraph_bool_t old) {
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_int_t *neis;
igraph_vector_t outdegree;
igraph_neimode_t dirmode;
igraph_adjlist_t allneis;
igraph_real_t maxdiff = eps;
igraph_vit_t vit;
IGRAPH_WARNING("igraph_pagerank_old is deprecated from igraph 0.7, "
"use igraph_pagerank instead");
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 = igraph_Calloc(no_of_nodes, igraph_real_t);
if (prvec == 0) {
IGRAPH_ERROR("pagerank failed", IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(igraph_free, prvec);
prvec_new = igraph_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 = igraph_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_adjlist_init(graph, &allneis, dirmode);
IGRAPH_FINALLY(igraph_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) {
igraph_real_t sumfrom = 0, sum = 0;
niter--;
maxdiff = 0;
/* Calculate the quotient of the actual PageRank value and the
* outdegree for every node */
sumfrom = 0.0; sum = 0.0;
for (i = 0; i < no_of_nodes; i++) {
sumfrom += prvec[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_adjlist_get(&allneis, i);
n = igraph_vector_int_size(neis);
for (j = 0; j < n; j++) {
long int neighbor = (long int) VECTOR(*neis)[j];
prvec_new[i] += prvec_scaled[neighbor];
}
prvec_new[i] *= damping;
if (!old) {
prvec_new[i] += (1 - damping) / no_of_nodes;
} else {
prvec_new[i] += (1 - damping);
}
sum += prvec_new[i];
}
for (i = 0; i < no_of_nodes; i++) {
if (!old) {
prvec_new[i] /= sum;
}
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_adjlist_destroy(&allneis);
igraph_vit_destroy(&vit);
igraph_vector_destroy(&outdegree);
igraph_Free(prvec);
igraph_Free(prvec_new);
igraph_Free(prvec_scaled);
IGRAPH_FINALLY_CLEAN(6);
return 0;
}
/* Not declared static so that the testsuite can use it, but not part of the public API. */
int igraph_rewire_core(igraph_t *graph, igraph_integer_t n, igraph_rewiring_t mode, igraph_bool_t use_adjlist) {
long int no_of_nodes = igraph_vcount(graph);
long int no_of_edges = igraph_ecount(graph);
char message[256];
igraph_integer_t a, b, c, d, dummy, num_swaps, num_successful_swaps;
igraph_vector_t eids, edgevec, alledges;
igraph_bool_t directed, loops, ok;
igraph_es_t es;
igraph_adjlist_t al;
if (no_of_nodes < 4) {
IGRAPH_ERROR("graph unsuitable for rewiring", IGRAPH_EINVAL);
}
directed = igraph_is_directed(graph);
loops = (mode & IGRAPH_REWIRING_SIMPLE_LOOPS);
RNG_BEGIN();
IGRAPH_VECTOR_INIT_FINALLY(&eids, 2);
if (use_adjlist) {
/* As well as the sorted adjacency list, we maintain an unordered
* list of edges for picking a random edge in constant time.
*/
IGRAPH_CHECK(igraph_adjlist_init(graph, &al, IGRAPH_OUT));
IGRAPH_FINALLY(igraph_adjlist_destroy, &al);
IGRAPH_VECTOR_INIT_FINALLY(&alledges, no_of_edges * 2);
igraph_get_edgelist(graph, &alledges, /*bycol=*/ 0);
} else {
IGRAPH_VECTOR_INIT_FINALLY(&edgevec, 4);
es = igraph_ess_vector(&eids);
}
/* 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)
*/
num_swaps = num_successful_swaps = 0;
while (num_swaps < n) {
IGRAPH_ALLOW_INTERRUPTION();
if (num_swaps % 1000 == 0) {
snprintf(message, sizeof(message),
"Random rewiring (%.2f%% of the trials were successful)",
num_swaps > 0 ? ((100.0 * num_successful_swaps) / num_swaps) : 0.0);
IGRAPH_PROGRESS(message, (100.0 * num_swaps) / n, 0);
}
switch (mode) {
case IGRAPH_REWIRING_SIMPLE:
case IGRAPH_REWIRING_SIMPLE_LOOPS:
ok = 1;
/* Choose two edges randomly */
VECTOR(eids)[0] = RNG_INTEGER(0, no_of_edges - 1);
do {
VECTOR(eids)[1] = RNG_INTEGER(0, no_of_edges - 1);
} while (VECTOR(eids)[0] == VECTOR(eids)[1]);
/* Get the endpoints */
if (use_adjlist) {
a = VECTOR(alledges)[((igraph_integer_t)VECTOR(eids)[0]) * 2];
b = VECTOR(alledges)[(((igraph_integer_t)VECTOR(eids)[0]) * 2) + 1];
c = VECTOR(alledges)[((igraph_integer_t)VECTOR(eids)[1]) * 2];
d = VECTOR(alledges)[(((igraph_integer_t)VECTOR(eids)[1]) * 2) + 1];
} else {
IGRAPH_CHECK(igraph_edge(graph, (igraph_integer_t) VECTOR(eids)[0],
&a, &b));
IGRAPH_CHECK(igraph_edge(graph, (igraph_integer_t) VECTOR(eids)[1],
&c, &d));
}
/* For an undirected graph, we have two "variants" of each edge, i.e.
* a -- b and b -- a. Since some rewirings can be performed only when we
* "swap" the endpoints, we do it now with probability 0.5 */
if (!directed && RNG_UNIF01() < 0.5) {
dummy = c; c = d; d = dummy;
if (use_adjlist) {
/* Flip the edge in the unordered edge-list, so the update later on
* hits the correct end. */
VECTOR(alledges)[((igraph_integer_t)VECTOR(eids)[1]) * 2] = c;
VECTOR(alledges)[(((igraph_integer_t)VECTOR(eids)[1]) * 2) + 1] = d;
}
}
/* If we do not touch loops, check whether a == b or c == d and disallow
* the swap if needed */
if (!loops && (a == b || c == d)) {
ok = 0;
} else {
/* Check whether they are suitable for rewiring */
if (a == c || b == d) {
/* Swapping would have no effect */
ok = 0;
} else {
/* a != c && b != d */
/* If a == d or b == c, the swap would generate at least one loop, so
* we disallow them unless we want to have loops */
ok = loops || (a != d && b != c);
/* Also, if a == b and c == d and we allow loops, doing the swap
* would result in a multiple edge if the graph is undirected */
ok = ok && (directed || a != b || c != d);
}
}
/* All good so far. Now check for the existence of a --> d and c --> b to
* disallow the creation of multiple edges */
if (ok) {
if (use_adjlist) {
if (igraph_adjlist_has_edge(&al, a, d, directed)) {
ok = 0;
}
} else {
IGRAPH_CHECK(igraph_are_connected(graph, a, d, &ok));
ok = !ok;
}
}
if (ok) {
if (use_adjlist) {
if (igraph_adjlist_has_edge(&al, c, b, directed)) {
ok = 0;
}
} else {
IGRAPH_CHECK(igraph_are_connected(graph, c, b, &ok));
ok = !ok;
}
}
/* If we are still okay, we can perform the rewiring */
if (ok) {
/* printf("Deleting: %ld -> %ld, %ld -> %ld\n",
(long)a, (long)b, (long)c, (long)d); */
if (use_adjlist) {
// Replace entry in sorted adjlist:
IGRAPH_CHECK(igraph_adjlist_replace_edge(&al, a, b, d, directed));
IGRAPH_CHECK(igraph_adjlist_replace_edge(&al, c, d, b, directed));
// Also replace in unsorted edgelist:
VECTOR(alledges)[(((igraph_integer_t)VECTOR(eids)[0]) * 2) + 1] = d;
VECTOR(alledges)[(((igraph_integer_t)VECTOR(eids)[1]) * 2) + 1] = b;
} else {
IGRAPH_CHECK(igraph_delete_edges(graph, es));
VECTOR(edgevec)[0] = a; VECTOR(edgevec)[1] = d;
VECTOR(edgevec)[2] = c; VECTOR(edgevec)[3] = b;
/* printf("Adding: %ld -> %ld, %ld -> %ld\n",
(long)a, (long)d, (long)c, (long)b); */
igraph_add_edges(graph, &edgevec, 0);
}
num_successful_swaps++;
}
break;
default:
RNG_END();
IGRAPH_ERROR("unknown rewiring mode", IGRAPH_EINVMODE);
}
num_swaps++;
}
if (use_adjlist) {
/* Replace graph edges with the adjlist current state */
IGRAPH_CHECK(igraph_delete_edges(graph, igraph_ess_all(IGRAPH_EDGEORDER_ID)));
IGRAPH_CHECK(igraph_add_edges(graph, &alledges, 0));
}
IGRAPH_PROGRESS("Random rewiring: ", 100.0, 0);
if (use_adjlist) {
igraph_vector_destroy(&alledges);
igraph_adjlist_destroy(&al);
} else {
igraph_vector_destroy(&edgevec);
}
igraph_vector_destroy(&eids);
IGRAPH_FINALLY_CLEAN(use_adjlist ? 3 : 2);
RNG_END();
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 allocated. If you would like to keep the original graph intact, use
* \ref igraph_copy() beforehand.
*
* \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 flags:
* \clist
* \cli 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. The method will
* neither destroy nor create self-loops.
* \cli IGRAPH_REWIRING_SIMPLE_LOOPS
* Same as \c IGRAPH_REWIRING_SIMPLE but allows the creation or
* destruction of self-loops.
* \endclist
*
* \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.
*
* \example examples/simple/igraph_rewire.c
*/
#define REWIRE_ADJLIST_THRESHOLD 10
int igraph_rewire(igraph_t *graph, igraph_integer_t n, igraph_rewiring_t mode) {
igraph_bool_t use_adjlist = n >= REWIRE_ADJLIST_THRESHOLD;
return igraph_rewire_core(graph, n, mode, use_adjlist);
}
/**
* Subgraph creation, old version: it copies the graph and then deletes
* unneeded vertices.
*/
int igraph_i_subgraph_copy_and_delete(const igraph_t *graph, igraph_t *res,
const igraph_vs_t vids,
igraph_vector_t *map,
igraph_vector_t *invmap) {
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 = igraph_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));
}
}
igraph_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_idx(res, igraph_vss_vector(&delete),
map, invmap));
igraph_vector_destroy(&delete);
igraph_vit_destroy(&vit);
IGRAPH_FINALLY_CLEAN(3);
return 0;
}
/**
* Subgraph creation, new version: creates the new graph instead of
* copying the old one.
*/
int igraph_i_subgraph_create_from_scratch(const igraph_t *graph,
igraph_t *res,
const igraph_vs_t vids,
igraph_vector_t *map,
igraph_vector_t *invmap) {
igraph_bool_t directed = igraph_is_directed(graph);
long int no_of_nodes = igraph_vcount(graph);
long int no_of_new_nodes = 0;
long int i, j, n;
long int to;
igraph_integer_t eid;
igraph_vector_t vids_old2new, vids_new2old;
igraph_vector_t eids_new2old;
igraph_vector_t nei_edges;
igraph_vector_t new_edges;
igraph_vit_t vit;
igraph_vector_t *my_vids_old2new = &vids_old2new,
*my_vids_new2old = &vids_new2old;
/* The order of initialization is important here, they will be destroyed in the
* opposite order */
IGRAPH_VECTOR_INIT_FINALLY(&eids_new2old, 0);
if (invmap) {
my_vids_new2old = invmap;
igraph_vector_clear(my_vids_new2old);
} else {
IGRAPH_VECTOR_INIT_FINALLY(&vids_new2old, 0);
}
IGRAPH_VECTOR_INIT_FINALLY(&new_edges, 0);
IGRAPH_VECTOR_INIT_FINALLY(&nei_edges, 0);
if (map) {
my_vids_old2new = map;
IGRAPH_CHECK(igraph_vector_resize(map, no_of_nodes));
igraph_vector_null(map);
} else {
IGRAPH_VECTOR_INIT_FINALLY(&vids_old2new, no_of_nodes);
}
IGRAPH_CHECK(igraph_vit_create(graph, vids, &vit));
IGRAPH_FINALLY(igraph_vit_destroy, &vit);
/* Calculate the mapping from the old node IDs to the new ones. The other
* igraph_simplify implementation in igraph_i_simplify_copy_and_delete
* ensures that the order of vertex IDs is kept during remapping (i.e.
* if the old ID of vertex A is less than the old ID of vertex B, then
* the same will also be true for the new IDs). To ensure compatibility
* with the other implementation, we have to fetch the vertex IDs into
* a vector first and then sort it. We temporarily use new_edges for that.
*/
IGRAPH_CHECK(igraph_vit_as_vector(&vit, &nei_edges));
igraph_vit_destroy(&vit);
IGRAPH_FINALLY_CLEAN(1);
igraph_vector_sort(&nei_edges);
n = igraph_vector_size(&nei_edges);
for (i = 0; i < n; i++) {
long int vid = (long int) VECTOR(nei_edges)[i];
if (VECTOR(*my_vids_old2new)[vid] == 0) {
IGRAPH_CHECK(igraph_vector_push_back(my_vids_new2old, vid));
no_of_new_nodes++;
VECTOR(*my_vids_old2new)[vid] = no_of_new_nodes;
}
}
/* Create the new edge list */
for (i = 0; i < no_of_new_nodes; i++) {
long int old_vid = (long int) VECTOR(*my_vids_new2old)[i];
long int new_vid = i;
IGRAPH_CHECK(igraph_incident(graph, &nei_edges, old_vid, IGRAPH_OUT));
n = igraph_vector_size(&nei_edges);
if (directed) {
for (j = 0; j < n; j++) {
eid = (igraph_integer_t) VECTOR(nei_edges)[j];
to = (long int) VECTOR(*my_vids_old2new)[ (long int)IGRAPH_TO(graph, eid) ];
if (!to) {
continue;
}
IGRAPH_CHECK(igraph_vector_push_back(&new_edges, new_vid));
IGRAPH_CHECK(igraph_vector_push_back(&new_edges, to - 1));
IGRAPH_CHECK(igraph_vector_push_back(&eids_new2old, eid));
}
} else {
for (j = 0; j < n; j++) {
eid = (igraph_integer_t) VECTOR(nei_edges)[j];
if (IGRAPH_FROM(graph, eid) != old_vid) { /* avoid processing edges twice */
continue;
}
to = (long int) VECTOR(*my_vids_old2new)[ (long int)IGRAPH_TO(graph, eid) ];
if (!to) {
continue;
}
IGRAPH_CHECK(igraph_vector_push_back(&new_edges, new_vid));
IGRAPH_CHECK(igraph_vector_push_back(&new_edges, to - 1));
IGRAPH_CHECK(igraph_vector_push_back(&eids_new2old, eid));
}
}
}
/* Get rid of some vectors that are not needed anymore */
if (!map) {
igraph_vector_destroy(&vids_old2new);
IGRAPH_FINALLY_CLEAN(1);
}
igraph_vector_destroy(&nei_edges);
IGRAPH_FINALLY_CLEAN(1);
/* Create the new graph */
IGRAPH_CHECK(igraph_create(res, &new_edges, (igraph_integer_t)
no_of_new_nodes, directed));
IGRAPH_I_ATTRIBUTE_DESTROY(res);
/* Now we can also get rid of the new_edges vector */
igraph_vector_destroy(&new_edges);
IGRAPH_FINALLY_CLEAN(1);
/* Make sure that the newly created graph is destroyed if something happens from
* now on */
IGRAPH_FINALLY(igraph_destroy, res);
/* Copy the graph attributes */
IGRAPH_CHECK(igraph_i_attribute_copy(res, graph,
/* ga = */ 1, /* va = */ 0, /* ea = */ 0));
/* Copy the vertex attributes */
IGRAPH_CHECK(igraph_i_attribute_permute_vertices(graph, res,
my_vids_new2old));
/* Copy the edge attributes */
IGRAPH_CHECK(igraph_i_attribute_permute_edges(graph, res, &eids_new2old));
if (!invmap) {
igraph_vector_destroy(my_vids_new2old);
IGRAPH_FINALLY_CLEAN(1);
}
igraph_vector_destroy(&eids_new2old);
IGRAPH_FINALLY_CLEAN(2); /* 1 + 1 since we don't need to destroy res */
return 0;
}
/**
* \ingroup structural
* \function igraph_subgraph
* \brief Creates a subgraph induced by the specified vertices.
*
* </para><para>
* This function is an alias to \ref igraph_induced_subgraph(), it is
* left here to ensure API compatibility with igraph versions prior to 0.6.
*
* </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 zero, 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 A vertex selector describing which vertices to keep.
* \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,
const igraph_vs_t vids) {
IGRAPH_WARNING("igraph_subgraph is deprecated from igraph 0.6, "
"use igraph_induced_subgraph instead");
return igraph_induced_subgraph(graph, res, vids, IGRAPH_SUBGRAPH_AUTO);
}
/**
* \ingroup structural
* \function igraph_induced_subgraph
* \brief Creates a subgraph induced by 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 zero, 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 A vertex selector describing which vertices to keep.
* \param impl This parameter selects which implementation should we
* use when constructing the new graph. Basically there are two
* possibilities: \c IGRAPH_SUBGRAPH_COPY_AND_DELETE copies the
* existing graph and deletes the vertices that are not needed
* in the new graph, while \c IGRAPH_SUBGRAPH_CREATE_FROM_SCRATCH
* constructs the new graph from scratch without copying the old
* one. The latter is more efficient if you are extracting a
* relatively small subpart of a very large graph, while the
* former is better if you want to extract a subgraph whose size
* is comparable to the size of the whole graph. There is a third
* possibility: \c IGRAPH_SUBGRAPH_AUTO will select one of the
* two methods automatically based on the ratio of the number
* of vertices in the new and the old graph.
*
* \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_induced_subgraph(const igraph_t *graph, igraph_t *res,
const igraph_vs_t vids, igraph_subgraph_implementation_t impl) {
return igraph_induced_subgraph_map(graph, res, vids, impl, /* map= */ 0,
/* invmap= */ 0);
}
int igraph_i_induced_subgraph_suggest_implementation(
const igraph_t *graph, const igraph_vs_t vids,
igraph_subgraph_implementation_t *result) {
double ratio;
igraph_integer_t num_vs;
if (igraph_vs_is_all(&vids)) {
ratio = 1.0;
} else {
IGRAPH_CHECK(igraph_vs_size(graph, &vids, &num_vs));
ratio = (igraph_real_t) num_vs / igraph_vcount(graph);
}
/* TODO: needs benchmarking; threshold was chosen totally arbitrarily */
if (ratio > 0.5) {
*result = IGRAPH_SUBGRAPH_COPY_AND_DELETE;
} else {
*result = IGRAPH_SUBGRAPH_CREATE_FROM_SCRATCH;
}
return 0;
}
int igraph_induced_subgraph_map(const igraph_t *graph, igraph_t *res,
const igraph_vs_t vids,
igraph_subgraph_implementation_t impl,
igraph_vector_t *map,
igraph_vector_t *invmap) {
if (impl == IGRAPH_SUBGRAPH_AUTO) {
IGRAPH_CHECK(igraph_i_induced_subgraph_suggest_implementation(graph, vids, &impl));
}
switch (impl) {
case IGRAPH_SUBGRAPH_COPY_AND_DELETE:
return igraph_i_subgraph_copy_and_delete(graph, res, vids, map, invmap);
case IGRAPH_SUBGRAPH_CREATE_FROM_SCRATCH:
return igraph_i_subgraph_create_from_scratch(graph, res, vids, map,
invmap);
default:
IGRAPH_ERROR("unknown subgraph implementation type", IGRAPH_EINVAL);
}
return 0;
}
/**
* \ingroup structural
* \function igraph_subgraph_edges
* \brief Creates a subgraph with the specified edges and their endpoints.
*
* </para><para>
* This function collects the specified edges and their endpoints to a new
* graph.
* As the vertex ids in a graph always start with zero, 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 eids An edge selector describing which edges to keep.
* \param delete_vertices Whether to delete the vertices not incident on any
* of the specified edges as well. If \c FALSE, the number of vertices
* in the result graph will always be equal to the number of vertices
* in the input graph.
* \return Error code:
* \c IGRAPH_ENOMEM, not enough memory for
* temporary data.
* \c IGRAPH_EINVEID, invalid edge id in
* \p eids.
*
* Time complexity: O(|V|+|E|),
* |V| and
* |E| are the number of vertices and
* edges in the original graph.
*
* \sa \ref igraph_delete_edges() to delete the specified set of
* edges from a graph, the opposite of this function.
*/
int igraph_subgraph_edges(const igraph_t *graph, igraph_t *res,
const igraph_es_t eids, igraph_bool_t delete_vertices) {
long int no_of_nodes = igraph_vcount(graph);
long int no_of_edges = igraph_ecount(graph);
igraph_vector_t delete = IGRAPH_VECTOR_NULL;
char *vremain, *eremain;
long int i;
igraph_eit_t eit;
IGRAPH_CHECK(igraph_eit_create(graph, eids, &eit));
IGRAPH_FINALLY(igraph_eit_destroy, &eit);
IGRAPH_VECTOR_INIT_FINALLY(&delete, 0);
vremain = igraph_Calloc(no_of_nodes, char);
if (vremain == 0) {
IGRAPH_ERROR("subgraph_edges failed", IGRAPH_ENOMEM);
}
eremain = igraph_Calloc(no_of_edges, char);
if (eremain == 0) {
IGRAPH_ERROR("subgraph_edges failed", IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(free, vremain); /* TODO: hack */
IGRAPH_FINALLY(free, eremain); /* TODO: hack */
IGRAPH_CHECK(igraph_vector_reserve(&delete, no_of_edges - IGRAPH_EIT_SIZE(eit)));
/* Collect the vertex and edge IDs that will remain */
for (IGRAPH_EIT_RESET(eit); !IGRAPH_EIT_END(eit); IGRAPH_EIT_NEXT(eit)) {
igraph_integer_t from, to;
long int eid = (long int) IGRAPH_EIT_GET(eit);
IGRAPH_CHECK(igraph_edge(graph, (igraph_integer_t) eid, &from, &to));
eremain[eid] = vremain[(long int)from] = vremain[(long int)to] = 1;
}
/* Collect the edge IDs to be deleted */
for (i = 0; i < no_of_edges; i++) {
IGRAPH_ALLOW_INTERRUPTION();
if (eremain[i] == 0) {
IGRAPH_CHECK(igraph_vector_push_back(&delete, i));
}
}
igraph_Free(eremain);
IGRAPH_FINALLY_CLEAN(1);
/* Delete the unnecessary edges */
/* 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_edges(res, igraph_ess_vector(&delete)));
if (delete_vertices) {
/* Collect the vertex IDs to be deleted */
igraph_vector_clear(&delete);
for (i = 0; i < no_of_nodes; i++) {
IGRAPH_ALLOW_INTERRUPTION();
if (vremain[i] == 0) {
IGRAPH_CHECK(igraph_vector_push_back(&delete, i));
}
}
}
igraph_Free(vremain);
IGRAPH_FINALLY_CLEAN(1);
/* Delete the unnecessary vertices */
if (delete_vertices) {
IGRAPH_CHECK(igraph_delete_vertices(res, igraph_vss_vector(&delete)));
}
igraph_vector_destroy(&delete);
igraph_eit_destroy(&eit);
IGRAPH_FINALLY_CLEAN(3);
return 0;
}
void igraph_i_simplify_free(igraph_vector_ptr_t *p);
void igraph_i_simplify_free(igraph_vector_ptr_t *p) {
long int i, n = igraph_vector_ptr_size(p);
for (i = 0; i < n; i++) {
igraph_vector_t *v = VECTOR(*p)[i];
if (v) {
igraph_vector_destroy(v);
}
}
igraph_vector_ptr_destroy(p);
}
/**
* \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.
* \param edge_comb What to do with the edge attributes. See the igraph
* manual section about attributes for details.
* \return Error code:
* \c IGRAPH_ENOMEM if we are out of memory.
*
* Time complexity: O(|V|+|E|).
*
* \example examples/simple/igraph_simplify.c
*/
int igraph_simplify(igraph_t *graph, igraph_bool_t multiple,
igraph_bool_t loops,
const igraph_attribute_combination_t *edge_comb) {
igraph_vector_t edges = IGRAPH_VECTOR_NULL;
long int no_of_nodes = igraph_vcount(graph);
long int no_of_edges = igraph_ecount(graph);
long int edge;
igraph_bool_t attr = edge_comb && igraph_has_attribute_table();
long int from, to, pfrom = -1, pto = -2;
igraph_t res;
igraph_es_t es;
igraph_eit_t eit;
igraph_vector_t mergeinto;
long int actedge;
if (!multiple && !loops)
/* nothing to do */
{
return IGRAPH_SUCCESS;
}
if (!multiple) {
/* removing loop edges only, this is simple. No need to combine anything
* and the whole process can be done in-place */
IGRAPH_VECTOR_INIT_FINALLY(&edges, 0);
IGRAPH_CHECK(igraph_es_all(&es, IGRAPH_EDGEORDER_ID));
IGRAPH_FINALLY(igraph_es_destroy, &es);
IGRAPH_CHECK(igraph_eit_create(graph, es, &eit));
IGRAPH_FINALLY(igraph_eit_destroy, &eit);
while (!IGRAPH_EIT_END(eit)) {
edge = IGRAPH_EIT_GET(eit);
from = IGRAPH_FROM(graph, edge);
to = IGRAPH_TO(graph, edge);
if (from == to) {
IGRAPH_CHECK(igraph_vector_push_back(&edges, edge));
}
IGRAPH_EIT_NEXT(eit);
}
igraph_eit_destroy(&eit);
igraph_es_destroy(&es);
IGRAPH_FINALLY_CLEAN(2);
if (igraph_vector_size(&edges) > 0) {
IGRAPH_CHECK(igraph_delete_edges(graph, igraph_ess_vector(&edges)));
}
igraph_vector_destroy(&edges);
IGRAPH_FINALLY_CLEAN(1);
return IGRAPH_SUCCESS;
}
if (attr) {
IGRAPH_VECTOR_INIT_FINALLY(&mergeinto, no_of_edges);
}
IGRAPH_VECTOR_INIT_FINALLY(&edges, 0);
IGRAPH_CHECK(igraph_vector_reserve(&edges, no_of_edges * 2));
IGRAPH_CHECK(igraph_es_all(&es, IGRAPH_EDGEORDER_FROM));
IGRAPH_FINALLY(igraph_es_destroy, &es);
IGRAPH_CHECK(igraph_eit_create(graph, es, &eit));
IGRAPH_FINALLY(igraph_eit_destroy, &eit);
for (actedge = -1; !IGRAPH_EIT_END(eit); IGRAPH_EIT_NEXT(eit)) {
edge = IGRAPH_EIT_GET(eit);
from = IGRAPH_FROM(graph, edge);
to = IGRAPH_TO(graph, edge);
if (loops && from == to) {
/* Loop edge to be removed */
if (attr) {
VECTOR(mergeinto)[edge] = -1;
}
} else if (multiple && from == pfrom && to == pto) {
/* Multiple edge to be contracted */
if (attr) {
VECTOR(mergeinto)[edge] = actedge;
}
} else {
/* Edge to be kept */
igraph_vector_push_back(&edges, from);
igraph_vector_push_back(&edges, to);
if (attr) {
actedge++;
VECTOR(mergeinto)[edge] = actedge;
}
}
pfrom = from; pto = to;
}
igraph_eit_destroy(&eit);
igraph_es_destroy(&es);
IGRAPH_FINALLY_CLEAN(2);
IGRAPH_CHECK(igraph_create(&res, &edges, (igraph_integer_t) no_of_nodes,
igraph_is_directed(graph)));
igraph_vector_destroy(&edges);
IGRAPH_FINALLY_CLEAN(1);
IGRAPH_FINALLY(igraph_destroy, &res);
IGRAPH_I_ATTRIBUTE_DESTROY(&res);
IGRAPH_I_ATTRIBUTE_COPY(&res, graph, /*graph=*/ 1,
/*vertex=*/ 1, /*edge=*/ 0);
if (attr) {
igraph_fixed_vectorlist_t vl;
IGRAPH_CHECK(igraph_fixed_vectorlist_convert(&vl, &mergeinto,
actedge + 1));
IGRAPH_FINALLY(igraph_fixed_vectorlist_destroy, &vl);
IGRAPH_CHECK(igraph_i_attribute_combine_edges(graph, &res, &vl.v,
edge_comb));
igraph_fixed_vectorlist_destroy(&vl);
igraph_vector_destroy(&mergeinto);
IGRAPH_FINALLY_CLEAN(2);
}
IGRAPH_FINALLY_CLEAN(1);
igraph_destroy(graph);
*graph = res;
return 0;
}
/**
* \ingroup structural
* \function igraph_reciprocity
* \brief Calculates the reciprocity of a directed graph.
*
* </para><para>
* The measure of reciprocity defines the proportion of mutual
* connections, in a directed graph. It is most commonly defined as
* the probability that the opposite counterpart of a directed edge is
* also included in the graph. In adjacency matrix notation:
* <code>sum(i, j, (A.*A')ij) / sum(i, j, Aij)</code>, where
* <code>A.*A'</code> is the element-wise product of matrix
* <code>A</code> and its transpose. This measure is
* calculated if the \p mode argument is \c
* IGRAPH_RECIPROCITY_DEFAULT.
*
* </para><para>
* Prior to igraph version 0.6, another measure was implemented,
* defined as the probability of mutual connection between a vertex
* pair if we know that there is a (possibly non-mutual) connection
* between them. In other words, (unordered) vertex pairs are
* classified into three groups: (1) disconnected, (2)
* non-reciprocally connected, (3) reciprocally connected.
* The result is the size of group (3), divided by the sum of group
* sizes (2)+(3). This measure is calculated if \p mode is \c
* IGRAPH_RECIPROCITY_RATIO.
*
* \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.
* \param mode Type of reciprocity to calculate, possible values are
* \c IGRAPH_RECIPROCITY_DEFAULT and \c IGRAPH_RECIPROCITY_RATIO,
* please see their description above.
* \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.
*
* \example examples/simple/igraph_reciprocity.c
*/
int igraph_reciprocity(const igraph_t *graph, igraph_real_t *res,
igraph_bool_t ignore_loops,
igraph_reciprocity_t mode) {
igraph_integer_t nonrec = 0, rec = 0, loops = 0;
igraph_vector_t inneis, outneis;
long int i;
long int no_of_nodes = igraph_vcount(graph);
if (mode != IGRAPH_RECIPROCITY_DEFAULT &&
mode != IGRAPH_RECIPROCITY_RATIO) {
IGRAPH_ERROR("Invalid reciprocity type", IGRAPH_EINVAL);
}
/* 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, (igraph_integer_t) i, IGRAPH_IN);
igraph_neighbors(graph, &outneis, (igraph_integer_t) 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 (VECTOR(inneis)[ip] == i) {
loops += 1;
if (!ignore_loops) {
rec += 1;
}
} else {
rec += 1;
}
ip++;
op++;
}
}
nonrec += (igraph_vector_size(&inneis) - ip) +
(igraph_vector_size(&outneis) - op);
}
if (mode == IGRAPH_RECIPROCITY_DEFAULT) {
if (ignore_loops) {
*res = (igraph_real_t) rec / (igraph_ecount(graph) - loops);
} else {
*res = (igraph_real_t) rec / (igraph_ecount(graph));
}
} else if (mode == IGRAPH_RECIPROCITY_RATIO) {
*res = (igraph_real_t) 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 of 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, (igraph_integer_t) 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_incident(graph, &ineis_in, (igraph_integer_t) i,
IGRAPH_IN));
IGRAPH_CHECK(igraph_incident(graph, &ineis_out, (igraph_integer_t) i,
IGRAPH_OUT));
/* NaN for isolates */
if (igraph_vector_size(&ineis_in) == 0 &&
igraph_vector_size(&ineis_out) == 0) {
VECTOR(*res)[a] = IGRAPH_NAN;
}
/* zero their contribution */
for (b = 0; b < igraph_vector_size(&ineis_in); b++) {
edge = (igraph_integer_t) 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 = (igraph_integer_t) 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 contributions, in-neighbors and out-neighbors */
for (b = 0; b < igraph_vector_size(&ineis_in); b++) {
edge = (igraph_integer_t) 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 = (igraph_integer_t) 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 = (igraph_integer_t) VECTOR(ineis_in)[b];
igraph_edge(graph, edge, &from, &to);
if (to == i) {
to = from;
}
j = to;
if (i == j) {
continue;
}
IGRAPH_CHECK(igraph_incident(graph, &jneis_in, (igraph_integer_t) j,
IGRAPH_IN));
IGRAPH_CHECK(igraph_incident(graph, &jneis_out, (igraph_integer_t) j,
IGRAPH_OUT));
for (c = 0; c < igraph_vector_size(&jneis_in); c++) {
edge2 = (igraph_integer_t) 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 = (igraph_integer_t) 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 = (igraph_integer_t) VECTOR(ineis_out)[b];
igraph_edge(graph, edge, &from, &to);
if (to == i) {
to = from;
}
j = to;
if (i == j) {
continue;
}
IGRAPH_CHECK(igraph_incident(graph, &jneis_in, (igraph_integer_t) j,
IGRAPH_IN));
IGRAPH_CHECK(igraph_incident(graph, &jneis_out, (igraph_integer_t) j,
IGRAPH_OUT));
for (c = 0; c < igraph_vector_size(&jneis_in); c++) {
edge2 = (igraph_integer_t) 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 = (igraph_integer_t) 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 = (igraph_integer_t) 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 = (igraph_integer_t) 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 vids Vector giving the vertex IDs for which the maximum degree will
* be calculated.
* \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_integer_t) 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
* necessarily mean that the graph really contains any loops). If
* this is 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_real_t no_of_edges = igraph_ecount(graph);
igraph_bool_t directed = igraph_is_directed(graph);
if (no_of_nodes == 0) {
*res = IGRAPH_NAN;
return 0;
}
if (!loops) {
if (no_of_nodes == 1) {
*res = IGRAPH_NAN;
} else if (directed) {
*res = no_of_edges / no_of_nodes / (no_of_nodes - 1);
} else {
*res = no_of_edges / no_of_nodes * 2.0 / (no_of_nodes - 1);
}
} else {
if (directed) {
*res = no_of_edges / no_of_nodes / no_of_nodes;
} else {
*res = no_of_edges / no_of_nodes * 2.0 / (no_of_nodes + 1);
}
}
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.
* \param mindist The minimum distance to include a vertex in the counting.
* If this is one, then the starting vertex is not counted. If this is
* two, then its neighbors are not counted, either, etc.
* \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,
igraph_integer_t mindist) {
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);
}
if (mindist < 0 || mindist > order) {
IGRAPH_ERROR("Minimum distance should be between zero and order",
IGRAPH_EINVAL);
}
added = igraph_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 = mindist == 0 ? 1 : 0;
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 = (long int) igraph_dqueue_pop(&q);
long int actdist = (long int) igraph_dqueue_pop(&q);
long int n;
igraph_neighbors(graph, &neis, (igraph_integer_t) 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 = (long int) 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));
if (actdist + 1 >= mindist) {
size++;
}
}
}
} else {
/* we just count them, but don't add them */
for (j = 0; j < n; j++) {
long int nei = (long int) VECTOR(neis)[j];
if (added[nei] != i + 1) {
added[nei] = i + 1;
if (actdist + 1 >= mindist) {
size++;
}
}
}
}
} /* while q not empty */
VECTOR(*res)[i] = size;
} /* for VIT, i */
igraph_vector_destroy(&neis);
igraph_vit_destroy(&vit);
igraph_dqueue_destroy(&q);
igraph_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.
* \param mindist The minimum distance to include a vertex in the counting.
* If this is one, then the starting vertex is not counted. If this is
* two, then its neighbors are not counted, either, etc.
* \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, igraph_integer_t mindist) {
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);
}
if (mindist < 0 || mindist > order) {
IGRAPH_ERROR("Minimum distance should be between zero and order",
IGRAPH_EINVAL);
}
added = igraph_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);
if (mindist == 0) {
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 = (long int) igraph_dqueue_pop(&q);
long int actdist = (long int) igraph_dqueue_pop(&q);
long int n;
igraph_neighbors(graph, &neis, (igraph_integer_t) 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 = (long int) 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));
if (actdist + 1 >= mindist) {
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 = (long int) VECTOR(neis)[j];
if (added[nei] != i + 1) {
added[nei] = i + 1;
if (actdist + 1 >= mindist) {
IGRAPH_CHECK(igraph_vector_push_back(&tmp, nei));
}
}
}
}
} /* while q not empty */
newv = igraph_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);
igraph_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.
* \param mindist The minimum distance to include a vertex in the counting.
* If this is one, then the starting vertex is not counted. If this is
* two, then its neighbors are not counted, either, etc.
* \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,
igraph_integer_t mindist) {
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);
}
if (mindist < 0 || mindist > order) {
IGRAPH_ERROR("Minimum distance should be between zero and order",
IGRAPH_EINVAL);
}
added = igraph_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);
if (mindist == 0) {
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 = (long int) igraph_dqueue_pop(&q);
long int actdist = (long int) igraph_dqueue_pop(&q);
long int n;
igraph_neighbors(graph, &neis, (igraph_integer_t) 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 = (long int) 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));
if (actdist + 1 >= mindist) {
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 = (long int) VECTOR(neis)[j];
if (added[nei] != i + 1) {
added[nei] = i + 1;
if (actdist + 1 >= mindist) {
IGRAPH_CHECK(igraph_vector_push_back(&tmp, nei));
}
}
}
}
} /* while q not empty */
newg = igraph_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_induced_subgraph(graph, newg,
igraph_vss_vector(&tmp),
IGRAPH_SUBGRAPH_AUTO));
} 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);
igraph_Free(added);
IGRAPH_FINALLY_CLEAN(5);
return 0;
}
/**
* \function igraph_topological_sorting
* \brief 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.
*
* \sa \ref igraph_is_dag() if you are only interested in whether a given
* graph is a DAG or not, or \ref igraph_feedback_arc_set() to find a
* set of edges whose removal makes the graph a DAG.
*
* \example examples/simple/igraph_topological_sorting.c
*/
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)) {
igraph_real_t tmp = igraph_dqueue_pop(&sources); node = (long) tmp;
/* 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, (igraph_integer_t) 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;
}
/**
* \function igraph_is_dag
* Checks whether a graph is a directed acyclic graph (DAG) or not.
*
* </para><para>
* A directed acyclic graph (DAG) is a directed graph with no cycles.
*
* \param graph The input graph.
* \param res Pointer to a boolean constant, the result
* is stored here.
* \return Error code.
*
* Time complexity: O(|V|+|E|), where |V| and |E| are the number of
* vertices and edges in the original input graph.
*
* \sa \ref igraph_topological_sorting() to get a possible topological
* sorting of a DAG.
*/
int igraph_is_dag(const igraph_t* graph, igraph_bool_t *res) {
long int no_of_nodes = igraph_vcount(graph);
igraph_vector_t degrees, neis;
igraph_dqueue_t sources;
long int node, i, j, nei, vertices_left;
if (!igraph_is_directed(graph)) {
*res = 0;
return IGRAPH_SUCCESS;
}
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(), IGRAPH_OUT, 1));
vertices_left = no_of_nodes;
/* Do we have nodes with no incoming edges? */
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 edges and remove them */
while (!igraph_dqueue_empty(&sources)) {
igraph_real_t tmp = igraph_dqueue_pop(&sources); node = (long) tmp;
/* Exclude the node from further source searches */
VECTOR(degrees)[node] = -1;
vertices_left--;
/* Get the neighbors and decrease their degrees by one */
IGRAPH_CHECK(igraph_neighbors(graph, &neis, (igraph_integer_t) node,
IGRAPH_IN));
j = igraph_vector_size(&neis);
for (i = 0; i < j; i++) {
nei = (long)VECTOR(neis)[i];
if (nei == node) {
continue;
}
VECTOR(degrees)[nei]--;
if (VECTOR(degrees)[nei] == 0) {
IGRAPH_CHECK(igraph_dqueue_push(&sources, nei));
}
}
}
*res = (vertices_left == 0);
if (vertices_left < 0) {
IGRAPH_WARNING("vertices_left < 0 in igraph_is_dag, possible bug");
}
igraph_vector_destroy(°rees);
igraph_vector_destroy(&neis);
igraph_dqueue_destroy(&sources);
IGRAPH_FINALLY_CLEAN(3);
return IGRAPH_SUCCESS;
}
/**
* \function igraph_is_simple
* \brief Decides whether the input graph is a simple graph.
*
* </para><para>
* A graph is a simple graph if it does not contain loop edges and
* multiple edges.
*
* \param graph The input graph.
* \param res Pointer to a boolean constant, the result
* is stored here.
* \return Error code.
*
* \sa \ref igraph_is_loop() and \ref igraph_is_multiple() to
* find the loops and multiple edges, \ref igraph_simplify() to
* get rid of them, or \ref igraph_has_multiple() to decide whether
* there is at least one multiple edge.
*
* Time complexity: O(|V|+|E|).
*/
int igraph_is_simple(const igraph_t *graph, igraph_bool_t *res) {
long int vc = igraph_vcount(graph);
long int ec = igraph_ecount(graph);
if (vc == 0 || ec == 0) {
*res = 1;
} else {
igraph_vector_t neis;
long int i, j, n;
igraph_bool_t found = 0;
IGRAPH_VECTOR_INIT_FINALLY(&neis, 0);
for (i = 0; i < vc; i++) {
igraph_neighbors(graph, &neis, (igraph_integer_t) i, IGRAPH_OUT);
n = igraph_vector_size(&neis);
for (j = 0; j < n; j++) {
if (VECTOR(neis)[j] == i) {
found = 1; break;
}
if (j > 0 && VECTOR(neis)[j - 1] == VECTOR(neis)[j]) {
found = 1; break;
}
}
}
*res = !found;
igraph_vector_destroy(&neis);
IGRAPH_FINALLY_CLEAN(1);
}
return 0;
}
/**
* \function igraph_has_loop
* \brief Returns whether the graph has at least one loop edge.
*
* </para><para>
* A loop edge is an edge from a vertex to itself.
* \param graph The input graph.
* \param res Pointer to an initialized boolean vector for storing the result.
*
* \sa \ref igraph_simplify() to get rid of loop edges.
*
* Time complexity: O(e), the number of edges to check.
*
* \example examples/simple/igraph_has_loop.c
*/
int igraph_has_loop(const igraph_t *graph, igraph_bool_t *res) {
long int i, m = igraph_ecount(graph);
*res = 0;
for (i = 0; i < m; i++) {
if (IGRAPH_FROM(graph, i) == IGRAPH_TO(graph, i)) {
*res = 1;
break;
}
}
return 0;
}
/**
* \function igraph_is_loop
* \brief Find the loop edges in a graph.
*
* </para><para>
* A loop edge is an edge from a vertex to itself.
* \param graph The input graph.
* \param res Pointer to an initialized boolean vector for storing the result,
* it will be resized as needed.
* \param es The edges to check, for all edges supply \ref igraph_ess_all() here.
* \return Error code.
*
* \sa \ref igraph_simplify() to get rid of loop edges.
*
* Time complexity: O(e), the number of edges to check.
*
* \example examples/simple/igraph_is_loop.c
*/
int igraph_is_loop(const igraph_t *graph, igraph_vector_bool_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_bool_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;
}
/**
* \function igraph_has_multiple
* \brief Check whether the graph has at least one multiple edge.
*
* </para><para>
* An edge is a multiple edge if there is another
* edge with the same head and tail vertices in the graph.
*
* \param graph The input graph.
* \param res Pointer to a boolean variable, the result will be stored here.
* \return Error code.
*
* \sa \ref igraph_count_multiple(), \ref igraph_is_multiple() and \ref igraph_simplify().
*
* Time complexity: O(e*d), e is the number of edges to check and d is the
* average degree (out-degree in directed graphs) of the vertices at the
* tail of the edges.
*
* \example examples/simple/igraph_has_multiple.c
*/
int igraph_has_multiple(const igraph_t *graph, igraph_bool_t *res) {
long int vc = igraph_vcount(graph);
long int ec = igraph_ecount(graph);
igraph_bool_t directed = igraph_is_directed(graph);
if (vc == 0 || ec == 0) {
*res = 0;
} else {
igraph_vector_t neis;
long int i, j, n;
igraph_bool_t found = 0;
IGRAPH_VECTOR_INIT_FINALLY(&neis, 0);
for (i = 0; i < vc && !found; i++) {
IGRAPH_CHECK(igraph_neighbors(graph, &neis, (igraph_integer_t) i,
IGRAPH_OUT));
n = igraph_vector_size(&neis);
for (j = 1; j < n; j++) {
if (VECTOR(neis)[j - 1] == VECTOR(neis)[j]) {
/* If the graph is undirected, loop edges appear twice in the neighbor
* list, so check the next item as well */
if (directed) {
/* Directed, so this is a real multiple edge */
found = 1; break;
} else if (VECTOR(neis)[j - 1] != i) {
/* Undirected, but not a loop edge */
found = 1; break;
} else if (j < n - 1 && VECTOR(neis)[j] == VECTOR(neis)[j + 1]) {
/* Undirected, loop edge, multiple times */
found = 1; break;
}
}
}
}
*res = found;
igraph_vector_destroy(&neis);
IGRAPH_FINALLY_CLEAN(1);
}
return 0;
}
/**
* \function igraph_is_multiple
* \brief Find the multiple edges in a graph.
*
* </para><para>
* An edge is a multiple edge if there is another
* edge with the same head and tail vertices in the graph.
*
* </para><para>
* Note that this function returns true only for the second or more
* appearances of the multiple edges.
* \param graph The input graph.
* \param res Pointer to a boolean vector, the result will be stored
* here. It will be resized as needed.
* \param es The edges to check. Supply \ref igraph_ess_all() if you want
* to check all edges.
* \return Error code.
*
* \sa \ref igraph_count_multiple(), \ref igraph_has_multiple() and \ref igraph_simplify().
*
* Time complexity: O(e*d), e is the number of edges to check and d is the
* average degree (out-degree in directed graphs) of the vertices at the
* tail of the edges.
*
* \example examples/simple/igraph_is_multiple.c
*/
int igraph_is_multiple(const igraph_t *graph, igraph_vector_bool_t *res,
igraph_es_t es) {
igraph_eit_t eit;
long int i;
igraph_lazy_inclist_t inclist;
IGRAPH_CHECK(igraph_eit_create(graph, es, &eit));
IGRAPH_FINALLY(igraph_eit_destroy, &eit);
IGRAPH_CHECK(igraph_lazy_inclist_init(graph, &inclist, IGRAPH_OUT));
IGRAPH_FINALLY(igraph_lazy_inclist_destroy, &inclist);
IGRAPH_CHECK(igraph_vector_bool_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_lazy_inclist_get(&inclist,
(igraph_integer_t) from);
long int j, n = igraph_vector_size(neis);
VECTOR(*res)[i] = 0;
for (j = 0; j < n; j++) {
long int e2 = (long int) VECTOR(*neis)[j];
long int to2 = IGRAPH_OTHER(graph, e2, from);
if (to2 == to && e2 < e) {
VECTOR(*res)[i] = 1;
}
}
}
igraph_lazy_inclist_destroy(&inclist);
igraph_eit_destroy(&eit);
IGRAPH_FINALLY_CLEAN(2);
return 0;
}
/**
* \function igraph_count_multiple
* \brief Count the number of appearances of the edges in a graph.
*
* </para><para>
* If the graph has no multiple edges then the result vector will be
* filled with ones.
* (An edge is a multiple edge if there is another
* edge with the same head and tail vertices in the graph.)
*
* </para><para>
* \param graph The input graph.
* \param res Pointer to a vector, the result will be stored
* here. It will be resized as needed.
* \param es The edges to check. Supply \ref igraph_ess_all() if you want
* to check all edges.
* \return Error code.
*
* \sa \ref igraph_is_multiple() and \ref igraph_simplify().
*
* Time complexity: O(e*d), e is the number of edges to check and d is the
* average degree (out-degree in directed graphs) of the vertices at the
* tail of the edges.
*/
int igraph_count_multiple(const igraph_t *graph, igraph_vector_t *res, igraph_es_t es) {
igraph_eit_t eit;
long int i;
igraph_lazy_inclist_t inclist;
IGRAPH_CHECK(igraph_eit_create(graph, es, &eit));
IGRAPH_FINALLY(igraph_eit_destroy, &eit);
IGRAPH_CHECK(igraph_lazy_inclist_init(graph, &inclist, IGRAPH_OUT));
IGRAPH_FINALLY(igraph_lazy_inclist_destroy, &inclist);
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_lazy_inclist_get(&inclist,
(igraph_integer_t) from);
long int j, n = igraph_vector_size(neis);
VECTOR(*res)[i] = 0;
for (j = 0; j < n; j++) {
long int e2 = (long int) VECTOR(*neis)[j];
long int to2 = IGRAPH_OTHER(graph, e2, from);
if (to2 == to) {
VECTOR(*res)[i] += 1;
}
}
/* for loop edges, divide the result by two */
if (to == from) {
VECTOR(*res)[i] /= 2;
}
}
igraph_lazy_inclist_destroy(&inclist);
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.
*
* \example examples/simple/igraph_girth.c
*/
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_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_lazy_adjlist_init(graph, &adjlist, IGRAPH_ALL,
IGRAPH_SIMPLIFY));
IGRAPH_FINALLY(igraph_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 = (long int) igraph_dqueue_pop(&q);
long int actlevel = VECTOR(level)[actnode];
long int i, n;
if (actlevel >= stoplevel) {
break;
}
neis = igraph_lazy_adjlist_get(&adjlist, (igraph_integer_t) actnode);
n = igraph_vector_size(neis);
for (i = 0; i < n; i++) {
long int nei = (long int) 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 = (igraph_integer_t) 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 = (long int) igraph_dqueue_pop(&q);
neis = igraph_lazy_adjlist_get(&adjlist, (igraph_integer_t) actnode);
n = igraph_vector_size(neis);
for (i = 0; i < n; i++) {
long int nei = (long int) 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_lazy_adjlist_destroy(&adjlist);
IGRAPH_FINALLY_CLEAN(3);
return 0;
}
int igraph_i_linegraph_undirected(const igraph_t *graph, igraph_t *linegraph);
int igraph_i_linegraph_directed(const igraph_t *graph, igraph_t *linegraph);
/* Note to self: tried using adjacency lists instead of igraph_incident queries,
* with minimal performance improvements on a graph with 70K vertices and 360K
* edges. (1.09s instead of 1.10s). I think it's not worth the fuss. */
int igraph_i_linegraph_undirected(const igraph_t *graph, igraph_t *linegraph) {
long int no_of_edges = igraph_ecount(graph);
long int i, j, n;
igraph_vector_t adjedges, adjedges2;
igraph_vector_t edges;
long int prev = -1;
IGRAPH_VECTOR_INIT_FINALLY(&edges, 0);
IGRAPH_VECTOR_INIT_FINALLY(&adjedges, 0);
IGRAPH_VECTOR_INIT_FINALLY(&adjedges2, 0);
for (i = 0; i < no_of_edges; i++) {
long int from = IGRAPH_FROM(graph, i);
long int to = IGRAPH_TO(graph, i);
IGRAPH_ALLOW_INTERRUPTION();
if (from != prev) {
IGRAPH_CHECK(igraph_incident(graph, &adjedges, (igraph_integer_t) from,
IGRAPH_ALL));
}
n = igraph_vector_size(&adjedges);
for (j = 0; j < n; j++) {
long int e = (long int) VECTOR(adjedges)[j];
if (e < i) {
IGRAPH_CHECK(igraph_vector_push_back(&edges, i));
IGRAPH_CHECK(igraph_vector_push_back(&edges, e));
}
}
IGRAPH_CHECK(igraph_incident(graph, &adjedges2, (igraph_integer_t) to,
IGRAPH_ALL));
n = igraph_vector_size(&adjedges2);
for (j = 0; j < n; j++) {
long int e = (long int) VECTOR(adjedges2)[j];
if (e < i) {
IGRAPH_CHECK(igraph_vector_push_back(&edges, i));
IGRAPH_CHECK(igraph_vector_push_back(&edges, e));
}
}
prev = from;
}
igraph_vector_destroy(&adjedges);
igraph_vector_destroy(&adjedges2);
IGRAPH_FINALLY_CLEAN(2);
igraph_create(linegraph, &edges, (igraph_integer_t) no_of_edges,
igraph_is_directed(graph));
igraph_vector_destroy(&edges);
IGRAPH_FINALLY_CLEAN(1);
return 0;
}
int igraph_i_linegraph_directed(const igraph_t *graph, igraph_t *linegraph) {
long int no_of_edges = igraph_ecount(graph);
long int i, j, n;
igraph_vector_t adjedges;
igraph_vector_t edges;
long int prev = -1;
IGRAPH_VECTOR_INIT_FINALLY(&edges, 0);
IGRAPH_VECTOR_INIT_FINALLY(&adjedges, 0);
for (i = 0; i < no_of_edges; i++) {
long int from = IGRAPH_FROM(graph, i);
IGRAPH_ALLOW_INTERRUPTION();
if (from != prev) {
IGRAPH_CHECK(igraph_incident(graph, &adjedges, (igraph_integer_t) from,
IGRAPH_IN));
}
n = igraph_vector_size(&adjedges);
for (j = 0; j < n; j++) {
long int e = (long int) VECTOR(adjedges)[j];
IGRAPH_CHECK(igraph_vector_push_back(&edges, e));
IGRAPH_CHECK(igraph_vector_push_back(&edges, i));
}
prev = from;
}
igraph_vector_destroy(&adjedges);
IGRAPH_FINALLY_CLEAN(1);
igraph_create(linegraph, &edges, (igraph_integer_t) no_of_edges, igraph_is_directed(graph));
igraph_vector_destroy(&edges);
IGRAPH_FINALLY_CLEAN(1);
return 0;
}
/**
* \function igraph_linegraph
* \brief Create the line graph of a graph.
*
* The line graph L(G) of a G undirected graph is defined as follows.
* L(G) has one vertex for each edge in G and two vertices in L(G) are connected
* by an edge if their corresponding edges share an end point.
*
* </para><para>
* The line graph L(G) of a G directed graph is slightly different,
* L(G) has one vertex for each edge in G and two vertices in L(G) are connected
* by a directed edge if the target of the first vertex's corresponding edge
* is the same as the source of the second vertex's corresponding edge.
*
* </para><para>
* Edge \em i in the original graph will correspond to vertex \em i
* in the line graph.
*
* </para><para>
* The first version of this function was contributed by Vincent Matossian,
* thanks.
* \param graph The input graph, may be directed or undirected.
* \param linegraph Pointer to an uninitialized graph object, the
* result is stored here.
* \return Error code.
*
* Time complexity: O(|V|+|E|), the number of edges plus the number of vertices.
*/
int igraph_linegraph(const igraph_t *graph, igraph_t *linegraph) {
if (igraph_is_directed(graph)) {
return igraph_i_linegraph_directed(graph, linegraph);
} else {
return igraph_i_linegraph_undirected(graph, linegraph);
}
}
/**
* \function igraph_add_edge
* \brief Adds a single edge to a graph.
*
* </para><para>
* For directed graphs the edge points from \p from to \p to.
*
* </para><para>
* Note that if you want to add many edges to a big graph, then it is
* inefficient to add them one by one, it is better to collect them into
* a vector and add all of them via a single \ref igraph_add_edges() call.
* \param igraph The graph.
* \param from The id of the first vertex of the edge.
* \param to The id of the second vertex of the edge.
* \return Error code.
*
* \sa \ref igraph_add_edges() to add many edges, \ref
* igraph_delete_edges() to remove edges and \ref
* igraph_add_vertices() to add vertices.
*
* Time complexity: O(|V|+|E|), the number of edges plus the number of
* vertices.
*/
int igraph_add_edge(igraph_t *graph, igraph_integer_t from, igraph_integer_t to) {
igraph_vector_t edges;
int ret;
IGRAPH_VECTOR_INIT_FINALLY(&edges, 2);
VECTOR(edges)[0] = from;
VECTOR(edges)[1] = to;
IGRAPH_CHECK(ret = igraph_add_edges(graph, &edges, 0));
igraph_vector_destroy(&edges);
IGRAPH_FINALLY_CLEAN(1);
return ret;
}
/*
* \example examples/simple/graph_convergence_degree.c
*/
int igraph_convergence_degree(const igraph_t *graph, igraph_vector_t *result,
igraph_vector_t *ins, igraph_vector_t *outs) {
long int no_of_nodes = igraph_vcount(graph);
long int no_of_edges = igraph_ecount(graph);
long int i, j, k, n;
long int *geodist;
igraph_vector_int_t *eids;
igraph_vector_t *ins_p, *outs_p, ins_v, outs_v;
igraph_dqueue_t q;
igraph_inclist_t inclist;
igraph_bool_t directed = igraph_is_directed(graph);
if (result != 0) {
IGRAPH_CHECK(igraph_vector_resize(result, no_of_edges));
}
IGRAPH_CHECK(igraph_dqueue_init(&q, 100));
IGRAPH_FINALLY(igraph_dqueue_destroy, &q);
if (ins == 0) {
ins_p = &ins_v;
IGRAPH_VECTOR_INIT_FINALLY(ins_p, no_of_edges);
} else {
ins_p = ins;
IGRAPH_CHECK(igraph_vector_resize(ins_p, no_of_edges));
igraph_vector_null(ins_p);
}
if (outs == 0) {
outs_p = &outs_v;
IGRAPH_VECTOR_INIT_FINALLY(outs_p, no_of_edges);
} else {
outs_p = outs;
IGRAPH_CHECK(igraph_vector_resize(outs_p, no_of_edges));
igraph_vector_null(outs_p);
}
geodist = igraph_Calloc(no_of_nodes, long int);
if (geodist == 0) {
IGRAPH_ERROR("Cannot calculate convergence degrees", IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(igraph_free, geodist);
/* Collect shortest paths originating from/to every node to correctly
* determine input field sizes */
for (k = 0; k < (directed ? 2 : 1); k++) {
igraph_neimode_t neimode = (k == 0) ? IGRAPH_OUT : IGRAPH_IN;
igraph_real_t *vec;
IGRAPH_CHECK(igraph_inclist_init(graph, &inclist, neimode));
IGRAPH_FINALLY(igraph_inclist_destroy, &inclist);
vec = (k == 0) ? VECTOR(*ins_p) : VECTOR(*outs_p);
for (i = 0; i < no_of_nodes; i++) {
igraph_dqueue_clear(&q);
memset(geodist, 0, sizeof(long int) * (size_t) no_of_nodes);
geodist[i] = 1;
IGRAPH_CHECK(igraph_dqueue_push(&q, i));
IGRAPH_CHECK(igraph_dqueue_push(&q, 0.0));
while (!igraph_dqueue_empty(&q)) {
long int actnode = (long int) igraph_dqueue_pop(&q);
long int actdist = (long int) igraph_dqueue_pop(&q);
IGRAPH_ALLOW_INTERRUPTION();
eids = igraph_inclist_get(&inclist, actnode);
n = igraph_vector_int_size(eids);
for (j = 0; j < n; j++) {
long int neighbor = IGRAPH_OTHER(graph, VECTOR(*eids)[j], actnode);
if (geodist[neighbor] != 0) {
/* we've already seen this node, another shortest path? */
if (geodist[neighbor] - 1 == actdist + 1) {
/* Since this edge is in the BFS tree rooted at i, we must
* increase either the size of the infield or the outfield */
if (!directed) {
if (actnode < neighbor) {
VECTOR(*ins_p)[(long int)VECTOR(*eids)[j]] += 1;
} else {
VECTOR(*outs_p)[(long int)VECTOR(*eids)[j]] += 1;
}
} else {
vec[(long int)VECTOR(*eids)[j]] += 1;
}
} else if (geodist[neighbor] - 1 < actdist + 1) {
continue;
}
} else {
/* we haven't seen this node yet */
IGRAPH_CHECK(igraph_dqueue_push(&q, neighbor));
IGRAPH_CHECK(igraph_dqueue_push(&q, actdist + 1));
/* Since this edge is in the BFS tree rooted at i, we must
* increase either the size of the infield or the outfield */
if (!directed) {
if (actnode < neighbor) {
VECTOR(*ins_p)[(long int)VECTOR(*eids)[j]] += 1;
} else {
VECTOR(*outs_p)[(long int)VECTOR(*eids)[j]] += 1;
}
} else {
vec[(long int)VECTOR(*eids)[j]] += 1;
}
geodist[neighbor] = actdist + 2;
}
}
}
}
igraph_inclist_destroy(&inclist);
IGRAPH_FINALLY_CLEAN(1);
}
if (result != 0) {
for (i = 0; i < no_of_edges; i++)
VECTOR(*result)[i] = (VECTOR(*ins_p)[i] - VECTOR(*outs_p)[i]) /
(VECTOR(*ins_p)[i] + VECTOR(*outs_p)[i]);
if (!directed) {
for (i = 0; i < no_of_edges; i++)
if (VECTOR(*result)[i] < 0) {
VECTOR(*result)[i] = -VECTOR(*result)[i];
}
}
}
if (ins == 0) {
igraph_vector_destroy(ins_p);
IGRAPH_FINALLY_CLEAN(1);
}
if (outs == 0) {
igraph_vector_destroy(outs_p);
IGRAPH_FINALLY_CLEAN(1);
}
igraph_free(geodist);
igraph_dqueue_destroy(&q);
IGRAPH_FINALLY_CLEAN(2);
return 0;
}
/**
* \function igraph_shortest_paths_dijkstra
* Weighted shortest paths from some sources.
*
* This function is Dijkstra's algorithm to find the weighted
* shortest paths to all vertices from a single source. (It is run
* independently for the given sources.) It uses a binary heap for
* efficient implementation.
*
* \param graph The input graph, can be directed.
* \param res The result, a matrix. A pointer to an initialized matrix
* should be passed here. The matrix will be resized as needed.
* Each row contains the distances from a single source, to the
* vertices given in the \c to argument.
* Unreachable vertices has distance
* \c IGRAPH_INFINITY.
* \param from The source vertices.
* \param to The target vertices. It is not allowed to include a
* vertex twice or more.
* \param weights The edge weights. They must be all non-negative for
* Dijkstra's algorithm to work. An error code is returned if there
* is a negative edge weight in the weight vector. If this is a null
* pointer, then the
* unweighted version, \ref igraph_shortest_paths() is called.
* \param mode For directed graphs; whether to follow paths along edge
* directions (\c IGRAPH_OUT), or the opposite (\c IGRAPH_IN), or
* ignore edge directions completely (\c IGRAPH_ALL). It is ignored
* for undirected graphs.
* \return Error code.
*
* Time complexity: O(s*|E|log|E|+|V|), where |V| is the number of
* vertices, |E| the number of edges and s the number of sources.
*
* \sa \ref igraph_shortest_paths() for a (slightly) faster unweighted
* version or \ref igraph_shortest_paths_bellman_ford() for a weighted
* variant that works in the presence of negative edge weights (but no
* negative loops).
*
* \example examples/simple/dijkstra.c
*/
int igraph_shortest_paths_dijkstra(const igraph_t *graph,
igraph_matrix_t *res,
const igraph_vs_t from,
const igraph_vs_t to,
const igraph_vector_t *weights,
igraph_neimode_t mode) {
/* Implementation details. This is the basic Dijkstra algorithm,
with a binary heap. The heap is indexed, i.e. it stores not only
the distances, but also which vertex they belong to.
From now on we use a 2-way heap, so the distances can be queried
directly from the heap.
Dirty tricks:
- the opposite of the distance is stored in the heap, as it is a
maximum heap and we need a minimum heap.
- we don't use IGRAPH_INFINITY in the res matrix during the
computation, as IGRAPH_FINITE() might involve a function call
and we want to spare that. -1 will denote infinity instead.
*/
long int no_of_nodes = igraph_vcount(graph);
long int no_of_edges = igraph_ecount(graph);
igraph_2wheap_t Q;
igraph_vit_t fromvit, tovit;
long int no_of_from, no_of_to;
igraph_lazy_inclist_t inclist;
long int i, j;
igraph_real_t my_infinity = IGRAPH_INFINITY;
igraph_bool_t all_to;
igraph_vector_t indexv;
if (!weights) {
return igraph_shortest_paths(graph, res, from, to, mode);
}
if (igraph_vector_size(weights) != no_of_edges) {
IGRAPH_ERROR("Weight vector length does not match", IGRAPH_EINVAL);
}
if (igraph_vector_min(weights) < 0) {
IGRAPH_ERROR("Weight vector must be non-negative", IGRAPH_EINVAL);
}
IGRAPH_CHECK(igraph_vit_create(graph, from, &fromvit));
IGRAPH_FINALLY(igraph_vit_destroy, &fromvit);
no_of_from = IGRAPH_VIT_SIZE(fromvit);
IGRAPH_CHECK(igraph_2wheap_init(&Q, no_of_nodes));
IGRAPH_FINALLY(igraph_2wheap_destroy, &Q);
IGRAPH_CHECK(igraph_lazy_inclist_init(graph, &inclist, mode));
IGRAPH_FINALLY(igraph_lazy_inclist_destroy, &inclist);
if ( (all_to = igraph_vs_is_all(&to)) ) {
no_of_to = no_of_nodes;
} else {
IGRAPH_VECTOR_INIT_FINALLY(&indexv, no_of_nodes);
IGRAPH_CHECK(igraph_vit_create(graph, to, &tovit));
IGRAPH_FINALLY(igraph_vit_destroy, &tovit);
no_of_to = IGRAPH_VIT_SIZE(tovit);
for (i = 0; !IGRAPH_VIT_END(tovit); IGRAPH_VIT_NEXT(tovit)) {
long int v = IGRAPH_VIT_GET(tovit);
if (VECTOR(indexv)[v]) {
IGRAPH_ERROR("Duplicate vertices in `to', this is not allowed",
IGRAPH_EINVAL);
}
VECTOR(indexv)[v] = ++i;
}
}
IGRAPH_CHECK(igraph_matrix_resize(res, no_of_from, no_of_to));
igraph_matrix_fill(res, my_infinity);
for (IGRAPH_VIT_RESET(fromvit), i = 0;
!IGRAPH_VIT_END(fromvit);
IGRAPH_VIT_NEXT(fromvit), i++) {
long int reached = 0;
long int source = IGRAPH_VIT_GET(fromvit);
igraph_2wheap_clear(&Q);
igraph_2wheap_push_with_index(&Q, source, -1.0);
while (!igraph_2wheap_empty(&Q)) {
long int minnei = igraph_2wheap_max_index(&Q);
igraph_real_t mindist = -igraph_2wheap_deactivate_max(&Q);
igraph_vector_t *neis;
long int nlen;
if (all_to) {
MATRIX(*res, i, minnei) = mindist - 1.0;
} else {
if (VECTOR(indexv)[minnei]) {
MATRIX(*res, i, (long int)(VECTOR(indexv)[minnei] - 1)) = mindist - 1.0;
reached++;
if (reached == no_of_to) {
igraph_2wheap_clear(&Q);
break;
}
}
}
/* Now check all neighbors of 'minnei' for a shorter path */
neis = igraph_lazy_inclist_get(&inclist, (igraph_integer_t) minnei);
nlen = igraph_vector_size(neis);
for (j = 0; j < nlen; j++) {
long int edge = (long int) VECTOR(*neis)[j];
long int tto = IGRAPH_OTHER(graph, edge, minnei);
igraph_real_t altdist = mindist + VECTOR(*weights)[edge];
igraph_bool_t active = igraph_2wheap_has_active(&Q, tto);
igraph_bool_t has = igraph_2wheap_has_elem(&Q, tto);
igraph_real_t curdist = active ? -igraph_2wheap_get(&Q, tto) : 0.0;
if (!has) {
/* This is the first non-infinite distance */
IGRAPH_CHECK(igraph_2wheap_push_with_index(&Q, tto, -altdist));
} else if (altdist < curdist) {
/* This is a shorter path */
IGRAPH_CHECK(igraph_2wheap_modify(&Q, tto, -altdist));
}
}
} /* !igraph_2wheap_empty(&Q) */
} /* !IGRAPH_VIT_END(fromvit) */
if (!all_to) {
igraph_vit_destroy(&tovit);
igraph_vector_destroy(&indexv);
IGRAPH_FINALLY_CLEAN(2);
}
igraph_lazy_inclist_destroy(&inclist);
igraph_2wheap_destroy(&Q);
igraph_vit_destroy(&fromvit);
IGRAPH_FINALLY_CLEAN(3);
return 0;
}
/**
* \ingroup structural
* \function igraph_get_shortest_paths_dijkstra
* \brief Calculates the weighted shortest paths from/to one vertex.
*
* </para><para>
* If there is more than one path with the smallest weight between two vertices, this
* function gives only one of them.
* \param graph The graph object.
* \param vertices The result, the ids of the vertices along the paths.
* 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. Supply a null pointer here if you don't need
* these vectors. Normally, either this argument, or the \c
* edges should be non-null, but no error or warning is given
* if they are both null pointers.
* \param edges The result, the ids of the edges along the paths.
* 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. Supply a null pointer here if you don't need
* these vectors. Normally, either this argument, or the \c
* vertices should be non-null, but no error or warning is given
* if they are both null pointers.
* \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 weights a vector holding the edge weights. All weights must be
* positive.
* \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
* \param predecessors A pointer to an initialized igraph vector or null.
* If not null, a vector containing the predecessor of each vertex in
* the single source shortest path tree is returned here. The
* predecessor of vertex i in the tree is the vertex from which vertex i
* was reached. The predecessor of the start vertex (in the \c from
* argument) is itself by definition. If the predecessor is -1, it means
* that the given vertex was not reached from the source during the
* search. Note that the search terminates if all the vertices in
* \c to are reached.
* \param inbound_edges A pointer to an initialized igraph vector or null.
* If not null, a vector containing the inbound edge of each vertex in
* the single source shortest path tree is returned here. The
* inbound edge of vertex i in the tree is the edge via which vertex i
* was reached. The start vertex and vertices that were not reached
* during the search will have -1 in the corresponding entry of the
* vector. Note that the search terminates if all the vertices in
* \c to are reached.
* \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(|E|log|E|+|V|), where |V| is the number of
* vertices and |E| is the number of edges
*
* \sa \ref igraph_shortest_paths_dijkstra() if you only need the path length but
* not the paths themselves, \ref igraph_get_shortest_paths() if all edge
* weights are equal.
*
* \example examples/simple/igraph_get_shortest_paths_dijkstra.c
*/
int igraph_get_shortest_paths_dijkstra(const igraph_t *graph,
igraph_vector_ptr_t *vertices,
igraph_vector_ptr_t *edges,
igraph_integer_t from,
igraph_vs_t to,
const igraph_vector_t *weights,
igraph_neimode_t mode,
igraph_vector_long_t *predecessors,
igraph_vector_long_t *inbound_edges) {
/* Implementation details. This is the basic Dijkstra algorithm,
with a binary heap. The heap is indexed, i.e. it stores not only
the distances, but also which vertex they belong to. The other
mapping, i.e. getting the distance for a vertex is not in the
heap (that would by the double-indexed heap), but in the result
matrix.
Dirty tricks:
- the opposite of the distance is stored in the heap, as it is a
maximum heap and we need a minimum heap.
- we don't use IGRAPH_INFINITY in the distance vector during the
computation, as IGRAPH_FINITE() might involve a function call
and we want to spare that. So we store distance+1.0 instead of
distance, and zero denotes infinity.
- `parents' assigns the inbound edge IDs of all vertices in the
shortest path tree to the vertices. In this implementation, the
edge ID + 1 is stored, zero means unreachable vertices.
*/
long int no_of_nodes = igraph_vcount(graph);
long int no_of_edges = igraph_ecount(graph);
igraph_vit_t vit;
igraph_2wheap_t Q;
igraph_lazy_inclist_t inclist;
igraph_vector_t dists;
long int *parents;
igraph_bool_t *is_target;
long int i, to_reach;
if (!weights) {
return igraph_get_shortest_paths(graph, vertices, edges, from, to, mode,
predecessors, inbound_edges);
}
if (igraph_vector_size(weights) != no_of_edges) {
IGRAPH_ERROR("Weight vector length does not match", IGRAPH_EINVAL);
}
if (igraph_vector_min(weights) < 0) {
IGRAPH_ERROR("Weight vector must be non-negative", IGRAPH_EINVAL);
}
IGRAPH_CHECK(igraph_vit_create(graph, to, &vit));
IGRAPH_FINALLY(igraph_vit_destroy, &vit);
if (vertices && IGRAPH_VIT_SIZE(vit) != igraph_vector_ptr_size(vertices)) {
IGRAPH_ERROR("Size of `vertices' and `to' should match", IGRAPH_EINVAL);
}
if (edges && IGRAPH_VIT_SIZE(vit) != igraph_vector_ptr_size(edges)) {
IGRAPH_ERROR("Size of `edges' and `to' should match", IGRAPH_EINVAL);
}
IGRAPH_CHECK(igraph_2wheap_init(&Q, no_of_nodes));
IGRAPH_FINALLY(igraph_2wheap_destroy, &Q);
IGRAPH_CHECK(igraph_lazy_inclist_init(graph, &inclist, mode));
IGRAPH_FINALLY(igraph_lazy_inclist_destroy, &inclist);
IGRAPH_VECTOR_INIT_FINALLY(&dists, no_of_nodes);
igraph_vector_fill(&dists, -1.0);
parents = igraph_Calloc(no_of_nodes, long int);
if (parents == 0) {
IGRAPH_ERROR("Can't calculate shortest paths", IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(igraph_free, parents);
is_target = igraph_Calloc(no_of_nodes, igraph_bool_t);
if (is_target == 0) {
IGRAPH_ERROR("Can't calculate shortest paths", IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(igraph_free, is_target);
/* 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 (!is_target[ (long int) IGRAPH_VIT_GET(vit) ]) {
is_target[ (long int) IGRAPH_VIT_GET(vit) ] = 1;
} else {
to_reach--; /* this node was given multiple times */
}
}
VECTOR(dists)[(long int)from] = 0.0; /* zero distance */
parents[(long int)from] = 0;
igraph_2wheap_push_with_index(&Q, from, 0);
while (!igraph_2wheap_empty(&Q) && to_reach > 0) {
long int nlen, minnei = igraph_2wheap_max_index(&Q);
igraph_real_t mindist = -igraph_2wheap_delete_max(&Q);
igraph_vector_t *neis;
IGRAPH_ALLOW_INTERRUPTION();
if (is_target[minnei]) {
is_target[minnei] = 0;
to_reach--;
}
/* Now check all neighbors of 'minnei' for a shorter path */
neis = igraph_lazy_inclist_get(&inclist, (igraph_integer_t) minnei);
nlen = igraph_vector_size(neis);
for (i = 0; i < nlen; i++) {
long int edge = (long int) VECTOR(*neis)[i];
long int tto = IGRAPH_OTHER(graph, edge, minnei);
igraph_real_t altdist = mindist + VECTOR(*weights)[edge];
igraph_real_t curdist = VECTOR(dists)[tto];
if (curdist < 0) {
/* This is the first finite distance */
VECTOR(dists)[tto] = altdist;
parents[tto] = edge + 1;
IGRAPH_CHECK(igraph_2wheap_push_with_index(&Q, tto, -altdist));
} else if (altdist < curdist) {
/* This is a shorter path */
VECTOR(dists)[tto] = altdist;
parents[tto] = edge + 1;
IGRAPH_CHECK(igraph_2wheap_modify(&Q, tto, -altdist));
}
}
} /* !igraph_2wheap_empty(&Q) */
if (to_reach > 0) {
IGRAPH_WARNING("Couldn't reach some vertices");
}
/* Create `predecessors' if needed */
if (predecessors) {
IGRAPH_CHECK(igraph_vector_long_resize(predecessors, no_of_nodes));
for (i = 0; i < no_of_nodes; i++) {
if (i == from) {
/* i is the start vertex */
VECTOR(*predecessors)[i] = i;
} else if (parents[i] <= 0) {
/* i was not reached */
VECTOR(*predecessors)[i] = -1;
} else {
/* i was reached via the edge with ID = parents[i] - 1 */
VECTOR(*predecessors)[i] = IGRAPH_OTHER(graph, parents[i] - 1, i);
}
}
}
/* Create `inbound_edges' if needed */
if (inbound_edges) {
IGRAPH_CHECK(igraph_vector_long_resize(inbound_edges, no_of_nodes));
for (i = 0; i < no_of_nodes; i++) {
if (parents[i] <= 0) {
/* i was not reached */
VECTOR(*inbound_edges)[i] = -1;
} else {
/* i was reached via the edge with ID = parents[i] - 1 */
VECTOR(*inbound_edges)[i] = parents[i] - 1;
}
}
}
/* Reconstruct the shortest paths based on vertex and/or edge IDs */
if (vertices || edges) {
for (IGRAPH_VIT_RESET(vit), i = 0; !IGRAPH_VIT_END(vit); IGRAPH_VIT_NEXT(vit), i++) {
long int node = IGRAPH_VIT_GET(vit);
long int size, act, edge;
igraph_vector_t *vvec = 0, *evec = 0;
if (vertices) {
vvec = VECTOR(*vertices)[i];
igraph_vector_clear(vvec);
}
if (edges) {
evec = VECTOR(*edges)[i];
igraph_vector_clear(evec);
}
IGRAPH_ALLOW_INTERRUPTION();
size = 0;
act = node;
while (parents[act]) {
size++;
edge = parents[act] - 1;
act = IGRAPH_OTHER(graph, edge, act);
}
if (vvec) {
IGRAPH_CHECK(igraph_vector_resize(vvec, size + 1));
VECTOR(*vvec)[size] = node;
}
if (evec) {
IGRAPH_CHECK(igraph_vector_resize(evec, size));
}
act = node;
while (parents[act]) {
edge = parents[act] - 1;
act = IGRAPH_OTHER(graph, edge, act);
size--;
if (vvec) {
VECTOR(*vvec)[size] = act;
}
if (evec) {
VECTOR(*evec)[size] = edge;
}
}
}
}
igraph_lazy_inclist_destroy(&inclist);
igraph_2wheap_destroy(&Q);
igraph_vector_destroy(&dists);
igraph_Free(is_target);
igraph_Free(parents);
igraph_vit_destroy(&vit);
IGRAPH_FINALLY_CLEAN(6);
return 0;
}
/**
* \function igraph_get_shortest_path_dijkstra
* Weighted shortest path from one vertex to another one.
*
* Calculates a single (positively) weighted shortest path from
* a single vertex to another one, using Dijkstra's algorithm.
*
* </para><para>This function is a special case (and a wrapper) to
* \ref igraph_get_shortest_paths_dijkstra().
*
* \param graph The input graph, it can be directed or undirected.
* \param vertices Pointer to an initialized vector or a null
* pointer. If not a null pointer, then the vertex ids along
* the path are stored here, including the source and target
* vertices.
* \param edges Pointer to an uninitialized vector or a null
* pointer. If not a null pointer, then the edge ids along the
* path are stored here.
* \param from The id of the source vertex.
* \param to The id of the target vertex.
* \param weights Vector of edge weights, in the order of edge
* ids. They must be non-negative, otherwise the algorithm does
* not work.
* \param mode A constant specifying how edge directions are
* considered in directed graphs. \c IGRAPH_OUT follows edge
* directions, \c IGRAPH_IN follows the opposite directions,
* and \c IGRAPH_ALL ignores edge directions. This argument is
* ignored for undirected graphs.
* \return Error code.
*
* Time complexity: O(|E|log|E|+|V|), |V| is the number of vertices,
* |E| is the number of edges in the graph.
*
* \sa \ref igraph_get_shortest_paths_dijkstra() for the version with
* more target vertices.
*/
int igraph_get_shortest_path_dijkstra(const igraph_t *graph,
igraph_vector_t *vertices,
igraph_vector_t *edges,
igraph_integer_t from,
igraph_integer_t to,
const igraph_vector_t *weights,
igraph_neimode_t mode) {
igraph_vector_ptr_t vertices2, *vp = &vertices2;
igraph_vector_ptr_t edges2, *ep = &edges2;
if (vertices) {
IGRAPH_CHECK(igraph_vector_ptr_init(&vertices2, 1));
IGRAPH_FINALLY(igraph_vector_ptr_destroy, &vertices2);
VECTOR(vertices2)[0] = vertices;
} else {
vp = 0;
}
if (edges) {
IGRAPH_CHECK(igraph_vector_ptr_init(&edges2, 1));
IGRAPH_FINALLY(igraph_vector_ptr_destroy, &edges2);
VECTOR(edges2)[0] = edges;
} else {
ep = 0;
}
IGRAPH_CHECK(igraph_get_shortest_paths_dijkstra(graph, vp, ep,
from, igraph_vss_1(to),
weights, mode, 0, 0));
if (edges) {
igraph_vector_ptr_destroy(&edges2);
IGRAPH_FINALLY_CLEAN(1);
}
if (vertices) {
igraph_vector_ptr_destroy(&vertices2);
IGRAPH_FINALLY_CLEAN(1);
}
return 0;
}
int igraph_i_vector_tail_cmp(const void* path1, const void* path2);
/* Compares two paths based on their last elements. Required by
* igraph_get_all_shortest_paths_dijkstra to put the final result
* in order. Assumes that both paths are pointers to igraph_vector_t
* objects and that they are not empty
*/
int igraph_i_vector_tail_cmp(const void* path1, const void* path2) {
return (int) (igraph_vector_tail(*(const igraph_vector_t**)path1) -
igraph_vector_tail(*(const igraph_vector_t**)path2));
}
/**
* \ingroup structural
* \function igraph_get_all_shortest_paths_dijkstra
* \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
* stored here for every vertex in the graph. Note that the values
* will be accurate only for those vertices that are in the target
* vertex sequence (see \p to), since the search terminates as soon
* as all the target vertices have been found.
* \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 weights a vector holding the edge weights. All weights must be
* non-negative.
* \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(|E|log|E|+|V|), where |V| is the number of
* vertices and |E| is the number of edges
*
* \sa \ref igraph_shortest_paths_dijkstra() if you only need the path
* length but not the paths themselves, \ref igraph_get_all_shortest_paths()
* if all edge weights are equal.
*
* \example examples/simple/igraph_get_all_shortest_paths_dijkstra.c
*/
int igraph_get_all_shortest_paths_dijkstra(const igraph_t *graph,
igraph_vector_ptr_t *res,
igraph_vector_t *nrgeo,
igraph_integer_t from, igraph_vs_t to,
const igraph_vector_t *weights,
igraph_neimode_t mode) {
/* Implementation details: see igraph_get_shortest_paths_dijkstra,
it's basically the same.
*/
long int no_of_nodes = igraph_vcount(graph);
long int no_of_edges = igraph_ecount(graph);
igraph_vit_t vit;
igraph_2wheap_t Q;
igraph_lazy_inclist_t inclist;
igraph_vector_t dists, order;
igraph_vector_ptr_t parents;
unsigned char *is_target;
long int i, n, to_reach;
if (!weights) {
return igraph_get_all_shortest_paths(graph, res, nrgeo, from, to, mode);
}
if (res == 0 && nrgeo == 0) {
return IGRAPH_SUCCESS;
}
if (igraph_vector_size(weights) != no_of_edges) {
IGRAPH_ERROR("Weight vector length does not match", IGRAPH_EINVAL);
}
if (igraph_vector_min(weights) < 0) {
IGRAPH_ERROR("Weight vector must be non-negative", IGRAPH_EINVAL);
}
/* parents stores a vector for each vertex, listing the parent vertices
* of each vertex in the traversal */
IGRAPH_CHECK(igraph_vector_ptr_init(&parents, no_of_nodes));
IGRAPH_FINALLY(igraph_vector_ptr_destroy_all, &parents);
igraph_vector_ptr_set_item_destructor(&parents, (igraph_finally_func_t*)igraph_vector_destroy);
for (i = 0; i < no_of_nodes; i++) {
igraph_vector_t* parent_vec;
parent_vec = igraph_Calloc(1, igraph_vector_t);
if (parent_vec == 0) {
IGRAPH_ERROR("cannot run igraph_get_all_shortest_paths", IGRAPH_ENOMEM);
}
IGRAPH_CHECK(igraph_vector_init(parent_vec, 0));
VECTOR(parents)[i] = parent_vec;
}
/* distance of each vertex from the root */
IGRAPH_VECTOR_INIT_FINALLY(&dists, no_of_nodes);
igraph_vector_fill(&dists, -1.0);
/* order lists the order of vertices in which they were found during
* the traversal */
IGRAPH_VECTOR_INIT_FINALLY(&order, 0);
/* boolean array to mark whether a given vertex is a target or not */
is_target = igraph_Calloc(no_of_nodes, unsigned char);
if (is_target == 0) {
IGRAPH_ERROR("Can't calculate shortest paths", IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(igraph_free, is_target);
/* two-way heap storing vertices and distances */
IGRAPH_CHECK(igraph_2wheap_init(&Q, no_of_nodes));
IGRAPH_FINALLY(igraph_2wheap_destroy, &Q);
/* lazy adjacency edge list to query neighbours efficiently */
IGRAPH_CHECK(igraph_lazy_inclist_init(graph, &inclist, mode));
IGRAPH_FINALLY(igraph_lazy_inclist_destroy, &inclist);
/* Mark the vertices we need to reach */
IGRAPH_CHECK(igraph_vit_create(graph, to, &vit));
IGRAPH_FINALLY(igraph_vit_destroy, &vit);
to_reach = IGRAPH_VIT_SIZE(vit);
for (IGRAPH_VIT_RESET(vit); !IGRAPH_VIT_END(vit); IGRAPH_VIT_NEXT(vit)) {
if (!is_target[ (long int) IGRAPH_VIT_GET(vit) ]) {
is_target[ (long int) IGRAPH_VIT_GET(vit) ] = 1;
} else {
to_reach--; /* this node was given multiple times */
}
}
igraph_vit_destroy(&vit);
IGRAPH_FINALLY_CLEAN(1);
VECTOR(dists)[(long int)from] = 0.0; /* zero distance */
igraph_2wheap_push_with_index(&Q, from, 0);
while (!igraph_2wheap_empty(&Q) && to_reach > 0) {
long int nlen, minnei = igraph_2wheap_max_index(&Q);
igraph_real_t mindist = -igraph_2wheap_delete_max(&Q);
igraph_vector_t *neis;
IGRAPH_ALLOW_INTERRUPTION();
/*
printf("Reached vertex %ld, is_target[%ld] = %d, %ld to go\n",
minnei, minnei, (int)is_target[minnei], to_reach - is_target[minnei]);
*/
if (is_target[minnei]) {
is_target[minnei] = 0;
to_reach--;
}
/* Mark that we have reached this vertex */
IGRAPH_CHECK(igraph_vector_push_back(&order, minnei));
/* Now check all neighbors of 'minnei' for a shorter path */
neis = igraph_lazy_inclist_get(&inclist, (igraph_integer_t) minnei);
nlen = igraph_vector_size(neis);
for (i = 0; i < nlen; i++) {
long int edge = (long int) VECTOR(*neis)[i];
long int tto = IGRAPH_OTHER(graph, edge, minnei);
igraph_real_t altdist = mindist + VECTOR(*weights)[edge];
igraph_real_t curdist = VECTOR(dists)[tto];
igraph_vector_t *parent_vec;
if (curdist < 0) {
/* This is the first non-infinite distance */
VECTOR(dists)[tto] = altdist;
parent_vec = (igraph_vector_t*)VECTOR(parents)[tto];
IGRAPH_CHECK(igraph_vector_push_back(parent_vec, minnei));
IGRAPH_CHECK(igraph_2wheap_push_with_index(&Q, tto, -altdist));
} else if (altdist == curdist && VECTOR(*weights)[edge] > 0) {
/* This is an alternative path with exactly the same length.
* Note that we consider this case only if the edge via which we
* reached the node has a nonzero weight; otherwise we could create
* infinite loops in undirected graphs by traversing zero-weight edges
* back-and-forth */
parent_vec = (igraph_vector_t*)VECTOR(parents)[tto];
IGRAPH_CHECK(igraph_vector_push_back(parent_vec, minnei));
} else if (altdist < curdist) {
/* This is a shorter path */
VECTOR(dists)[tto] = altdist;
parent_vec = (igraph_vector_t*)VECTOR(parents)[tto];
igraph_vector_clear(parent_vec);
IGRAPH_CHECK(igraph_vector_push_back(parent_vec, minnei));
IGRAPH_CHECK(igraph_2wheap_modify(&Q, tto, -altdist));
}
}
} /* !igraph_2wheap_empty(&Q) */
if (to_reach > 0) {
IGRAPH_WARNING("Couldn't reach some vertices");
}
/* we don't need these anymore */
igraph_lazy_inclist_destroy(&inclist);
igraph_2wheap_destroy(&Q);
IGRAPH_FINALLY_CLEAN(2);
/*
printf("Order:\n");
igraph_vector_print(&order);
printf("Parent vertices:\n");
for (i = 0; i < no_of_nodes; i++) {
if (igraph_vector_size(VECTOR(parents)[i]) > 0) {
printf("[%ld]: ", (long int)i);
igraph_vector_print(VECTOR(parents)[i]);
}
}
*/
if (nrgeo) {
IGRAPH_CHECK(igraph_vector_resize(nrgeo, no_of_nodes));
igraph_vector_null(nrgeo);
/* Theoretically, we could calculate nrgeo in parallel with the traversal.
* However, that way we would have to check whether nrgeo is null or not
* every time we want to update some element in nrgeo. Since we need the
* order vector anyway for building the final result, we could just as well
* build nrgeo here.
*/
VECTOR(*nrgeo)[(long int)from] = 1;
n = igraph_vector_size(&order);
for (i = 1; i < n; i++) {
long int node, j, k;
igraph_vector_t *parent_vec;
node = (long int)VECTOR(order)[i];
/* now, take the parent vertices */
parent_vec = (igraph_vector_t*)VECTOR(parents)[node];
k = igraph_vector_size(parent_vec);
for (j = 0; j < k; j++) {
VECTOR(*nrgeo)[node] += VECTOR(*nrgeo)[(long int)VECTOR(*parent_vec)[j]];
}
}
}
if (res) {
igraph_vector_t *path, *paths_index, *parent_vec;
igraph_stack_t stack;
long int j, node;
/* a shortest path from the starting vertex to vertex i can be
* obtained by calculating the shortest paths from the "parents"
* of vertex i in the traversal. Knowing which of the vertices
* are "targets" (see is_target), we can collect for which other
* vertices do we need to calculate the shortest paths. We reuse
* is_target for that; is_target = 0 means that we don't need the
* vertex, is_target = 1 means that the vertex is a target (hence
* we need it), is_target = 2 means that the vertex is not a target
* but it stands between a shortest path between the root and one
* of the targets
*/
if (igraph_vs_is_all(&to)) {
memset(is_target, 1, sizeof(unsigned char) * (size_t) no_of_nodes);
} else {
memset(is_target, 0, sizeof(unsigned char) * (size_t) no_of_nodes);
IGRAPH_CHECK(igraph_stack_init(&stack, 0));
IGRAPH_FINALLY(igraph_stack_destroy, &stack);
/* Add the target vertices to the queue */
IGRAPH_CHECK(igraph_vit_create(graph, to, &vit));
IGRAPH_FINALLY(igraph_vit_destroy, &vit);
for (IGRAPH_VIT_RESET(vit); !IGRAPH_VIT_END(vit); IGRAPH_VIT_NEXT(vit)) {
i = (long int) IGRAPH_VIT_GET(vit);
if (!is_target[i]) {
is_target[i] = 1;
IGRAPH_CHECK(igraph_stack_push(&stack, i));
}
}
igraph_vit_destroy(&vit);
IGRAPH_FINALLY_CLEAN(1);
while (!igraph_stack_empty(&stack)) {
/* For each parent of node i, get its parents */
igraph_real_t el = igraph_stack_pop(&stack);
parent_vec = (igraph_vector_t*)VECTOR(parents)[(long int) el];
i = igraph_vector_size(parent_vec);
for (j = 0; j < i; j++) {
/* For each parent, check if it's already in the stack.
* If not, push it and mark it in is_target */
n = (long int) VECTOR(*parent_vec)[j];
if (!is_target[n]) {
is_target[n] = 2;
IGRAPH_CHECK(igraph_stack_push(&stack, n));
}
}
}
igraph_stack_destroy(&stack);
IGRAPH_FINALLY_CLEAN(1);
}
/* now, reconstruct the shortest paths from the parent list in the
* order we've found the nodes during the traversal.
* dists is being re-used as a vector where element i tells the
* index in res where the shortest paths leading to vertex i
* start, plus one (so that zero means that there are no paths
* for a given vertex).
*/
paths_index = &dists;
n = igraph_vector_size(&order);
igraph_vector_null(paths_index);
/* clear the paths vector */
igraph_vector_ptr_clear(res);
igraph_vector_ptr_set_item_destructor(res,
(igraph_finally_func_t*)igraph_vector_destroy);
/* by definition, the shortest path leading to the starting vertex
* consists of the vertex itself only */
path = igraph_Calloc(1, igraph_vector_t);
if (path == 0)
IGRAPH_ERROR("cannot run igraph_get_all_shortest_paths_dijkstra",
IGRAPH_ENOMEM);
IGRAPH_FINALLY(igraph_free, path);
IGRAPH_CHECK(igraph_vector_init(path, 1));
IGRAPH_CHECK(igraph_vector_ptr_push_back(res, path));
IGRAPH_FINALLY_CLEAN(1); /* ownership of path passed to res */
VECTOR(*path)[0] = from;
VECTOR(*paths_index)[(long int)from] = 1;
for (i = 1; i < n; i++) {
long int m, path_count;
igraph_vector_t *parent_path;
node = (long int) VECTOR(order)[i];
/* if we don't need the shortest paths for this node (because
* it is not standing in a shortest path between the source
* node and any of the target nodes), skip it */
if (!is_target[node]) {
continue;
}
IGRAPH_ALLOW_INTERRUPTION();
/* we are calculating the shortest paths of node now. */
/* first, we update the paths_index */
path_count = igraph_vector_ptr_size(res);
VECTOR(*paths_index)[node] = path_count + 1;
/* res_end = (igraph_vector_t*)&(VECTOR(*res)[path_count]); */
/* now, take the parent vertices */
parent_vec = (igraph_vector_t*)VECTOR(parents)[node];
m = igraph_vector_size(parent_vec);
/*
printf("Calculating shortest paths to vertex %ld\n", node);
printf("Parents are: ");
igraph_vector_print(parent_vec);
*/
for (j = 0; j < m; j++) {
/* for each parent, copy the shortest paths leading to that parent
* and add the current vertex in the end */
long int parent_node = (long int) VECTOR(*parent_vec)[j];
long int parent_path_idx = (long int) VECTOR(*paths_index)[parent_node] - 1;
/*
printf(" Considering parent: %ld\n", parent_node);
printf(" Paths to parent start at index %ld in res\n", parent_path_idx);
*/
assert(parent_path_idx >= 0);
for (; parent_path_idx < path_count; parent_path_idx++) {
parent_path = (igraph_vector_t*)VECTOR(*res)[parent_path_idx];
if (igraph_vector_tail(parent_path) != parent_node) {
break;
}
path = igraph_Calloc(1, igraph_vector_t);
if (path == 0)
IGRAPH_ERROR("cannot run igraph_get_all_shortest_paths_dijkstra",
IGRAPH_ENOMEM);
IGRAPH_FINALLY(igraph_free, path);
IGRAPH_CHECK(igraph_vector_copy(path, parent_path));
IGRAPH_CHECK(igraph_vector_ptr_push_back(res, path));
IGRAPH_FINALLY_CLEAN(1); /* ownership of path passed to res */
IGRAPH_CHECK(igraph_vector_push_back(path, node));
}
}
}
/* remove the destructor from the path vector */
igraph_vector_ptr_set_item_destructor(res, 0);
/* free those paths from the result vector which we won't need */
n = igraph_vector_ptr_size(res);
j = 0;
for (i = 0; i < n; i++) {
igraph_real_t tmp;
path = (igraph_vector_t*)VECTOR(*res)[i];
tmp = igraph_vector_tail(path);
if (is_target[(long int)tmp] == 1) {
/* we need this path, keep it */
VECTOR(*res)[j] = path;
j++;
} else {
/* we don't need this path, free it */
igraph_vector_destroy(path); free(path);
}
}
IGRAPH_CHECK(igraph_vector_ptr_resize(res, j));
/* sort the paths by the target vertices */
igraph_vector_ptr_sort(res, igraph_i_vector_tail_cmp);
}
/* free the allocated memory */
igraph_vector_destroy(&order);
igraph_Free(is_target);
igraph_vector_destroy(&dists);
igraph_vector_ptr_destroy_all(&parents);
IGRAPH_FINALLY_CLEAN(4);
return 0;
}
/**
* \function igraph_shortest_paths_bellman_ford
* Weighted shortest paths from some sources allowing negative weights.
*
* This function is the Bellman-Ford algorithm to find the weighted
* shortest paths to all vertices from a single source. (It is run
* independently for the given sources.). If there are no negative
* weights, you are better off with \ref igraph_shortest_paths_dijkstra() .
*
* \param graph The input graph, can be directed.
* \param res The result, a matrix. A pointer to an initialized matrix
* should be passed here, the matrix will be resized if needed.
* Each row contains the distances from a single source, to all
* vertices in the graph, in the order of vertex ids. For unreachable
* vertices the matrix contains \c IGRAPH_INFINITY.
* \param from The source vertices.
* \param weights The edge weights. There mustn't be any closed loop in
* the graph that has a negative total weight (since this would allow
* us to decrease the weight of any path containing at least a single
* vertex of this loop infinitely). If this is a null pointer, then the
* unweighted version, \ref igraph_shortest_paths() is called.
* \param mode For directed graphs; whether to follow paths along edge
* directions (\c IGRAPH_OUT), or the opposite (\c IGRAPH_IN), or
* ignore edge directions completely (\c IGRAPH_ALL). It is ignored
* for undirected graphs.
* \return Error code.
*
* Time complexity: O(s*|E|*|V|), where |V| is the number of
* vertices, |E| the number of edges and s the number of sources.
*
* \sa \ref igraph_shortest_paths() for a faster unweighted version
* or \ref igraph_shortest_paths_dijkstra() if you do not have negative
* edge weights.
*
* \example examples/simple/bellman_ford.c
*/
int igraph_shortest_paths_bellman_ford(const igraph_t *graph,
igraph_matrix_t *res,
const igraph_vs_t from,
const igraph_vs_t to,
const igraph_vector_t *weights,
igraph_neimode_t mode) {
long int no_of_nodes = igraph_vcount(graph);
long int no_of_edges = igraph_ecount(graph);
igraph_lazy_inclist_t inclist;
long int i, j, k;
long int no_of_from, no_of_to;
igraph_dqueue_t Q;
igraph_vector_t clean_vertices;
igraph_vector_t num_queued;
igraph_vit_t fromvit, tovit;
igraph_real_t my_infinity = IGRAPH_INFINITY;
igraph_bool_t all_to;
igraph_vector_t dist;
/*
- speedup: a vertex is marked clean if its distance from the source
did not change during the last phase. Neighbors of a clean vertex
are not relaxed again, since it would mean no change in the
shortest path values. Dirty vertices are queued. Negative loops can
be detected by checking whether a vertex has been queued at least
n times.
*/
if (!weights) {
return igraph_shortest_paths(graph, res, from, to, mode);
}
if (igraph_vector_size(weights) != no_of_edges) {
IGRAPH_ERROR("Weight vector length does not match", IGRAPH_EINVAL);
}
IGRAPH_CHECK(igraph_vit_create(graph, from, &fromvit));
IGRAPH_FINALLY(igraph_vit_destroy, &fromvit);
no_of_from = IGRAPH_VIT_SIZE(fromvit);
IGRAPH_DQUEUE_INIT_FINALLY(&Q, no_of_nodes);
IGRAPH_VECTOR_INIT_FINALLY(&clean_vertices, no_of_nodes);
IGRAPH_VECTOR_INIT_FINALLY(&num_queued, no_of_nodes);
IGRAPH_CHECK(igraph_lazy_inclist_init(graph, &inclist, mode));
IGRAPH_FINALLY(igraph_lazy_inclist_destroy, &inclist);
if ( (all_to = igraph_vs_is_all(&to)) ) {
no_of_to = no_of_nodes;
} else {
IGRAPH_CHECK(igraph_vit_create(graph, to, &tovit));
IGRAPH_FINALLY(igraph_vit_destroy, &tovit);
no_of_to = IGRAPH_VIT_SIZE(tovit);
}
IGRAPH_VECTOR_INIT_FINALLY(&dist, no_of_nodes);
IGRAPH_CHECK(igraph_matrix_resize(res, no_of_from, no_of_to));
for (IGRAPH_VIT_RESET(fromvit), i = 0;
!IGRAPH_VIT_END(fromvit);
IGRAPH_VIT_NEXT(fromvit), i++) {
long int source = IGRAPH_VIT_GET(fromvit);
igraph_vector_fill(&dist, my_infinity);
VECTOR(dist)[source] = 0;
igraph_vector_null(&clean_vertices);
igraph_vector_null(&num_queued);
/* Fill the queue with vertices to be checked */
for (j = 0; j < no_of_nodes; j++) {
IGRAPH_CHECK(igraph_dqueue_push(&Q, j));
}
while (!igraph_dqueue_empty(&Q)) {
igraph_vector_t *neis;
long int nlen;
j = (long int) igraph_dqueue_pop(&Q);
VECTOR(clean_vertices)[j] = 1;
VECTOR(num_queued)[j] += 1;
if (VECTOR(num_queued)[j] > no_of_nodes) {
IGRAPH_ERROR("cannot run Bellman-Ford algorithm", IGRAPH_ENEGLOOP);
}
/* If we cannot get to j in finite time yet, there is no need to relax
* its edges */
if (!IGRAPH_FINITE(VECTOR(dist)[j])) {
continue;
}
neis = igraph_lazy_inclist_get(&inclist, (igraph_integer_t) j);
nlen = igraph_vector_size(neis);
for (k = 0; k < nlen; k++) {
long int nei = (long int) VECTOR(*neis)[k];
long int target = IGRAPH_OTHER(graph, nei, j);
if (VECTOR(dist)[target] > VECTOR(dist)[j] + VECTOR(*weights)[nei]) {
/* relax the edge */
VECTOR(dist)[target] = VECTOR(dist)[j] + VECTOR(*weights)[nei];
if (VECTOR(clean_vertices)[target]) {
VECTOR(clean_vertices)[target] = 0;
IGRAPH_CHECK(igraph_dqueue_push(&Q, target));
}
}
}
}
/* Copy it to the result */
if (all_to) {
igraph_matrix_set_row(res, &dist, i);
} else {
for (IGRAPH_VIT_RESET(tovit), j = 0; !IGRAPH_VIT_END(tovit);
IGRAPH_VIT_NEXT(tovit), j++) {
long int v = IGRAPH_VIT_GET(tovit);
MATRIX(*res, i, j) = VECTOR(dist)[v];
}
}
}
igraph_vector_destroy(&dist);
IGRAPH_FINALLY_CLEAN(1);
if (!all_to) {
igraph_vit_destroy(&tovit);
IGRAPH_FINALLY_CLEAN(1);
}
igraph_vit_destroy(&fromvit);
igraph_dqueue_destroy(&Q);
igraph_vector_destroy(&clean_vertices);
igraph_vector_destroy(&num_queued);
igraph_lazy_inclist_destroy(&inclist);
IGRAPH_FINALLY_CLEAN(5);
return 0;
}
/**
* \function igraph_shortest_paths_johnson
* Calculate shortest paths from some sources using Johnson's algorithm.
*
* See Wikipedia at http://en.wikipedia.org/wiki/Johnson's_algorithm
* for Johnson's algorithm. This algorithm works even if the graph
* contains negative edge weights, and it is worth using it if we
* calculate the shortest paths from many sources.
*
* </para><para> If no edge weights are supplied, then the unweighted
* version, \ref igraph_shortest_paths() is called.
*
* </para><para> If all the supplied edge weights are non-negative,
* then Dijkstra's algorithm is used by calling
* \ref igraph_shortest_paths_dijkstra().
*
* \param graph The input graph, typically it is directed.
* \param res Pointer to an initialized matrix, the result will be
* stored here, one line for each source vertex, one column for each
* target vertex.
* \param from The source vertices.
* \param to The target vertices. It is not allowed to include a
* vertex twice or more.
* \param weights Optional edge weights. If it is a null-pointer, then
* the unweighted breadth-first search based \ref
* igraph_shortest_paths() will be called.
* \return Error code.
*
* Time complexity: O(s|V|log|V|+|V||E|), |V| and |E| are the number
* of vertices and edges, s is the number of source vertices.
*
* \sa \ref igraph_shortest_paths() for a faster unweighted version
* or \ref igraph_shortest_paths_dijkstra() if you do not have negative
* edge weights, \ref igraph_shortest_paths_bellman_ford() if you only
* need to calculate shortest paths from a couple of sources.
*/
int igraph_shortest_paths_johnson(const igraph_t *graph,
igraph_matrix_t *res,
const igraph_vs_t from,
const igraph_vs_t to,
const igraph_vector_t *weights) {
long int no_of_nodes = igraph_vcount(graph);
long int no_of_edges = igraph_ecount(graph);
igraph_t newgraph;
igraph_vector_t edges, newweights;
igraph_matrix_t bfres;
long int i, ptr;
long int nr, nc;
igraph_vit_t fromvit;
/* If no weights, then we can just run the unweighted version */
if (!weights) {
return igraph_shortest_paths(graph, res, from, to, IGRAPH_OUT);
}
if (igraph_vector_size(weights) != no_of_edges) {
IGRAPH_ERROR("Weight vector length does not match", IGRAPH_EINVAL);
}
/* If no negative weights, then we can run Dijkstra's algorithm */
if (igraph_vector_min(weights) >= 0) {
return igraph_shortest_paths_dijkstra(graph, res, from, to,
weights, IGRAPH_OUT);
}
if (!igraph_is_directed(graph)) {
IGRAPH_ERROR("Johnson's shortest path: undirected graph and negative weight",
IGRAPH_EINVAL);
}
/* ------------------------------------------------------------ */
/* -------------------- Otherwise proceed --------------------- */
IGRAPH_MATRIX_INIT_FINALLY(&bfres, 0, 0);
IGRAPH_VECTOR_INIT_FINALLY(&newweights, 0);
IGRAPH_CHECK(igraph_empty(&newgraph, (igraph_integer_t) no_of_nodes + 1,
igraph_is_directed(graph)));
IGRAPH_FINALLY(igraph_destroy, &newgraph);
/* Add a new node to the graph, plus edges from it to all the others. */
IGRAPH_VECTOR_INIT_FINALLY(&edges, no_of_edges * 2 + no_of_nodes * 2);
igraph_get_edgelist(graph, &edges, /*bycol=*/ 0);
igraph_vector_resize(&edges, no_of_edges * 2 + no_of_nodes * 2);
for (i = 0, ptr = no_of_edges * 2; i < no_of_nodes; i++) {
VECTOR(edges)[ptr++] = no_of_nodes;
VECTOR(edges)[ptr++] = i;
}
IGRAPH_CHECK(igraph_add_edges(&newgraph, &edges, 0));
igraph_vector_destroy(&edges);
IGRAPH_FINALLY_CLEAN(1);
IGRAPH_CHECK(igraph_vector_reserve(&newweights, no_of_edges + no_of_nodes));
igraph_vector_update(&newweights, weights);
igraph_vector_resize(&newweights, no_of_edges + no_of_nodes);
for (i = no_of_edges; i < no_of_edges + no_of_nodes; i++) {
VECTOR(newweights)[i] = 0;
}
/* Run Bellmann-Ford algorithm on the new graph, starting from the
new vertex. */
IGRAPH_CHECK(igraph_shortest_paths_bellman_ford(&newgraph, &bfres,
igraph_vss_1((igraph_integer_t) no_of_nodes),
igraph_vss_all(), &newweights, IGRAPH_OUT));
igraph_destroy(&newgraph);
IGRAPH_FINALLY_CLEAN(1);
/* Now the edges of the original graph are reweighted, using the
values from the BF algorithm. Instead of w(u,v) we will have
w(u,v) + h(u) - h(v) */
igraph_vector_resize(&newweights, no_of_edges);
for (i = 0; i < no_of_edges; i++) {
long int ffrom = IGRAPH_FROM(graph, i);
long int tto = IGRAPH_TO(graph, i);
VECTOR(newweights)[i] += MATRIX(bfres, 0, ffrom) - MATRIX(bfres, 0, tto);
}
/* Run Dijkstra's algorithm on the new weights */
IGRAPH_CHECK(igraph_shortest_paths_dijkstra(graph, res, from,
to, &newweights,
IGRAPH_OUT));
igraph_vector_destroy(&newweights);
IGRAPH_FINALLY_CLEAN(1);
/* Reweight the shortest paths */
nr = igraph_matrix_nrow(res);
nc = igraph_matrix_ncol(res);
IGRAPH_CHECK(igraph_vit_create(graph, from, &fromvit));
IGRAPH_FINALLY(igraph_vit_destroy, &fromvit);
for (i = 0; i < nr; i++, IGRAPH_VIT_NEXT(fromvit)) {
long int v1 = IGRAPH_VIT_GET(fromvit);
if (igraph_vs_is_all(&to)) {
long int v2;
for (v2 = 0; v2 < nc; v2++) {
igraph_real_t sub = MATRIX(bfres, 0, v1) - MATRIX(bfres, 0, v2);
MATRIX(*res, i, v2) -= sub;
}
} else {
long int j;
igraph_vit_t tovit;
IGRAPH_CHECK(igraph_vit_create(graph, to, &tovit));
IGRAPH_FINALLY(igraph_vit_destroy, &tovit);
for (j = 0, IGRAPH_VIT_RESET(tovit); j < nc; j++, IGRAPH_VIT_NEXT(tovit)) {
long int v2 = IGRAPH_VIT_GET(tovit);
igraph_real_t sub = MATRIX(bfres, 0, v1) - MATRIX(bfres, 0, v2);
MATRIX(*res, i, v2) -= sub;
}
igraph_vit_destroy(&tovit);
IGRAPH_FINALLY_CLEAN(1);
}
}
igraph_vit_destroy(&fromvit);
igraph_matrix_destroy(&bfres);
IGRAPH_FINALLY_CLEAN(2);
return 0;
}
/**
* \function igraph_unfold_tree
* Unfolding a graph into a tree, by possibly multiplicating its vertices.
*
* A graph is converted into a tree (or forest, if it is unconnected),
* by performing a breadth-first search on it, and replicating
* vertices that were found a second, third, etc. time.
* \param graph The input graph, it can be either directed or
* undirected.
* \param tree Pointer to an uninitialized graph object, the result is
* stored here.
* \param mode For directed graphs; whether to follow paths along edge
* directions (\c IGRAPH_OUT), or the opposite (\c IGRAPH_IN), or
* ignore edge directions completely (\c IGRAPH_ALL). It is ignored
* for undirected graphs.
* \param roots A numeric vector giving the root vertex, or vertices
* (if the graph is not connected), to start from.
* \param vertex_index Pointer to an initialized vector, or a null
* pointer. If not a null pointer, then a mapping from the vertices
* in the new graph to the ones in the original is created here.
* \return Error code.
*
* Time complexity: O(n+m), linear in the number vertices and edges.
*
*/
int igraph_unfold_tree(const igraph_t *graph, igraph_t *tree,
igraph_neimode_t mode, const igraph_vector_t *roots,
igraph_vector_t *vertex_index) {
long int no_of_nodes = igraph_vcount(graph);
long int no_of_edges = igraph_ecount(graph);
long int no_of_roots = igraph_vector_size(roots);
long int tree_vertex_count = no_of_nodes;
igraph_vector_t edges;
igraph_vector_bool_t seen_vertices;
igraph_vector_bool_t seen_edges;
igraph_dqueue_t Q;
igraph_vector_t neis;
long int i, n, r, v_ptr = no_of_nodes;
/* TODO: handle not-connected graphs, multiple root vertices */
IGRAPH_VECTOR_INIT_FINALLY(&edges, 0);
igraph_vector_reserve(&edges, no_of_edges * 2);
IGRAPH_DQUEUE_INIT_FINALLY(&Q, 100);
IGRAPH_VECTOR_INIT_FINALLY(&neis, 0);
IGRAPH_VECTOR_BOOL_INIT_FINALLY(&seen_vertices, no_of_nodes);
IGRAPH_VECTOR_BOOL_INIT_FINALLY(&seen_edges, no_of_edges);
if (vertex_index) {
IGRAPH_CHECK(igraph_vector_resize(vertex_index, no_of_nodes));
for (i = 0; i < no_of_nodes; i++) {
VECTOR(*vertex_index)[i] = i;
}
}
for (r = 0; r < no_of_roots; r++) {
long int root = (long int) VECTOR(*roots)[r];
VECTOR(seen_vertices)[root] = 1;
igraph_dqueue_push(&Q, root);
while (!igraph_dqueue_empty(&Q)) {
long int actnode = (long int) igraph_dqueue_pop(&Q);
IGRAPH_CHECK(igraph_incident(graph, &neis, (igraph_integer_t) actnode, mode));
n = igraph_vector_size(&neis);
for (i = 0; i < n; i++) {
long int edge = (long int) VECTOR(neis)[i];
long int from = IGRAPH_FROM(graph, edge);
long int to = IGRAPH_TO(graph, edge);
long int nei = IGRAPH_OTHER(graph, edge, actnode);
if (! VECTOR(seen_edges)[edge]) {
VECTOR(seen_edges)[edge] = 1;
if (! VECTOR(seen_vertices)[nei]) {
igraph_vector_push_back(&edges, from);
igraph_vector_push_back(&edges, to);
VECTOR(seen_vertices)[nei] = 1;
IGRAPH_CHECK(igraph_dqueue_push(&Q, nei));
} else {
tree_vertex_count++;
if (vertex_index) {
IGRAPH_CHECK(igraph_vector_push_back(vertex_index, nei));
}
if (from == nei) {
igraph_vector_push_back(&edges, v_ptr++);
igraph_vector_push_back(&edges, to);
} else {
igraph_vector_push_back(&edges, from);
igraph_vector_push_back(&edges, v_ptr++);
}
}
}
} /* for i<n */
} /* ! igraph_dqueue_empty(&Q) */
} /* r < igraph_vector_size(roots) */
igraph_vector_bool_destroy(&seen_edges);
igraph_vector_bool_destroy(&seen_vertices);
igraph_vector_destroy(&neis);
igraph_dqueue_destroy(&Q);
IGRAPH_FINALLY_CLEAN(4);
IGRAPH_CHECK(igraph_create(tree, &edges, tree_vertex_count, igraph_is_directed(graph)));
igraph_vector_destroy(&edges);
IGRAPH_FINALLY_CLEAN(1);
return 0;
}
/**
* \function igraph_is_mutual
* Check whether the edges of a directed graph are mutual.
*
* An (A,B) edge is mutual if the graph contains the (B,A) edge, too.
* </para>
*
* <para>An undirected graph only has mutual edges, by definition.
* </para>
*
* <para>Edge multiplicity is not considered here, e.g. if there are two
* (A,B) edges and one (B,A) edge, then all three are considered to be
* mutual.
*
* \param graph The input graph.
* \param res Pointer to an initialized vector, the result is stored
* here.
* \param es The sequence of edges to check. Supply
* <code>igraph_ess_all()</code> for all edges, see \ref
* igraph_ess_all().
* \return Error code.
*
* Time complexity: O(n log(d)), n is the number of edges supplied, d
* is the maximum in-degree of the vertices that are targets of the
* supplied edges. An upper limit of the time complexity is O(n log(|E|)),
* |E| is the number of edges in the graph.
*/
int igraph_is_mutual(igraph_t *graph, igraph_vector_bool_t *res, igraph_es_t es) {
igraph_eit_t eit;
igraph_lazy_adjlist_t adjlist;
long int i;
/* How many edges do we have? */
IGRAPH_CHECK(igraph_eit_create(graph, es, &eit));
IGRAPH_FINALLY(igraph_eit_destroy, &eit);
IGRAPH_CHECK(igraph_vector_bool_resize(res, IGRAPH_EIT_SIZE(eit)));
/* An undirected graph has mutual edges by definition,
res is already properly resized */
if (! igraph_is_directed(graph)) {
igraph_vector_bool_fill(res, 1);
igraph_eit_destroy(&eit);
IGRAPH_FINALLY_CLEAN(1);
return 0;
}
IGRAPH_CHECK(igraph_lazy_adjlist_init(graph, &adjlist, IGRAPH_OUT, IGRAPH_DONT_SIMPLIFY));
IGRAPH_FINALLY(igraph_lazy_adjlist_destroy, &adjlist);
for (i = 0; ! IGRAPH_EIT_END(eit); i++, IGRAPH_EIT_NEXT(eit)) {
long int edge = IGRAPH_EIT_GET(eit);
long int from = IGRAPH_FROM(graph, edge);
long int to = IGRAPH_TO(graph, edge);
/* Check whether there is a to->from edge, search for from in the
out-list of to. We don't search an empty vector, because
vector_binsearch seems to have a bug with this. */
igraph_vector_t *neis = igraph_lazy_adjlist_get(&adjlist,
(igraph_integer_t) to);
if (igraph_vector_empty(neis)) {
VECTOR(*res)[i] = 0;
} else {
VECTOR(*res)[i] = igraph_vector_binsearch2(neis, from);
}
}
igraph_lazy_adjlist_destroy(&adjlist);
igraph_eit_destroy(&eit);
IGRAPH_FINALLY_CLEAN(2);
return 0;
}
int igraph_i_avg_nearest_neighbor_degree_weighted(const igraph_t *graph,
igraph_vs_t vids,
igraph_neimode_t mode,
igraph_neimode_t neighbor_degree_mode,
igraph_vector_t *knn,
igraph_vector_t *knnk,
const igraph_vector_t *weights);
int igraph_i_avg_nearest_neighbor_degree_weighted(const igraph_t *graph,
igraph_vs_t vids,
igraph_neimode_t mode,
igraph_neimode_t neighbor_degree_mode,
igraph_vector_t *knn,
igraph_vector_t *knnk,
const igraph_vector_t *weights) {
long int no_of_nodes = igraph_vcount(graph);
igraph_vector_t neis, edge_neis;
long int i, j, no_vids;
igraph_vit_t vit;
igraph_vector_t my_knn_v, *my_knn = knn;
igraph_vector_t strength, deg;
igraph_integer_t maxdeg;
igraph_vector_t deghist;
igraph_real_t mynan = IGRAPH_NAN;
if (igraph_vector_size(weights) != igraph_ecount(graph)) {
IGRAPH_ERROR("Invalid weight vector size", IGRAPH_EINVAL);
}
IGRAPH_CHECK(igraph_vit_create(graph, vids, &vit));
IGRAPH_FINALLY(igraph_vit_destroy, &vit);
no_vids = IGRAPH_VIT_SIZE(vit);
if (!knn) {
IGRAPH_VECTOR_INIT_FINALLY(&my_knn_v, no_vids);
my_knn = &my_knn_v;
} else {
IGRAPH_CHECK(igraph_vector_resize(knn, no_vids));
}
// Get degree of neighbours
IGRAPH_VECTOR_INIT_FINALLY(°, no_of_nodes);
IGRAPH_CHECK(igraph_degree(graph, °, igraph_vss_all(),
neighbor_degree_mode, IGRAPH_LOOPS));
IGRAPH_VECTOR_INIT_FINALLY(&strength, no_of_nodes);
// Get strength of all nodes
IGRAPH_CHECK(igraph_strength(graph, &strength, igraph_vss_all(),
mode, IGRAPH_LOOPS, weights));
// Get maximum degree for initialization
IGRAPH_CHECK(igraph_maxdegree(graph, &maxdeg, igraph_vss_all(),
mode, IGRAPH_LOOPS));
IGRAPH_VECTOR_INIT_FINALLY(&neis, (long int)maxdeg);
IGRAPH_VECTOR_INIT_FINALLY(&edge_neis, (long int)maxdeg);
igraph_vector_resize(&neis, 0);
igraph_vector_resize(&edge_neis, 0);
if (knnk) {
IGRAPH_CHECK(igraph_vector_resize(knnk, (long int)maxdeg));
igraph_vector_null(knnk);
IGRAPH_VECTOR_INIT_FINALLY(°hist, (long int)maxdeg);
}
for (i = 0; !IGRAPH_VIT_END(vit); IGRAPH_VIT_NEXT(vit), i++) {
igraph_real_t sum = 0.0;
long int v = IGRAPH_VIT_GET(vit);
long int nv;
igraph_real_t str = VECTOR(strength)[v];
// Get neighbours and incident edges
IGRAPH_CHECK(igraph_neighbors(graph, &neis, (igraph_integer_t) v, mode));
IGRAPH_CHECK(igraph_incident(graph, &edge_neis, (igraph_integer_t) v, mode));
nv = igraph_vector_size(&neis);
for (j = 0; j < nv; j++) {
long int nei = (long int) VECTOR(neis)[j];
long int e = (long int) VECTOR(edge_neis)[j];
double w = VECTOR(*weights)[e];
sum += w * VECTOR(deg)[nei];
}
if (str != 0.0) {
VECTOR(*my_knn)[i] = sum / str;
} else {
VECTOR(*my_knn)[i] = mynan;
}
if (knnk && nv > 0) {
VECTOR(*knnk)[nv - 1] += VECTOR(*my_knn)[i];
VECTOR(deghist)[nv - 1] += 1;
}
}
if (knnk) {
for (i = 0; i < maxdeg; i++) {
igraph_real_t dh = VECTOR(deghist)[i];
if (dh != 0) {
VECTOR(*knnk)[i] /= dh;
} else {
VECTOR(*knnk)[i] = mynan;
}
}
igraph_vector_destroy(°hist);
IGRAPH_FINALLY_CLEAN(1);
}
igraph_vector_destroy(&neis);
igraph_vector_destroy(°);
IGRAPH_FINALLY_CLEAN(2);
if (!knn) {
igraph_vector_destroy(&my_knn_v);
IGRAPH_FINALLY_CLEAN(1);
}
igraph_vit_destroy(&vit);
IGRAPH_FINALLY_CLEAN(1);
return 0;
}
/**
* \function igraph_avg_nearest_neighbor_degree
* Average nearest neighbor degree.
*
* Calculates the average degree of the neighbors for each vertex, and
* optionally, the same quantity in the function of vertex degree.
*
* </para><para>For isolate vertices \p knn is set to \c
* IGRAPH_NAN. The same is done in \p knnk for vertex degrees that
* don't appear in the graph.
*
* \param graph The input graph, it can be directed but the
* directedness of the edges is ignored.
* \param vids The vertices for which the calculation is performed.
* \param mode The neighbors over which is averaged.
* \param neighbor_degree_mode The degree of the neighbors which is
* averaged.
* \param vids The vertices for which the calculation is performed.
* \param knn Pointer to an initialized vector, the result will be
* stored here. It will be resized as needed. Supply a NULL pointer
* here, if you only want to calculate \c knnk.
* \param knnk Pointer to an initialized vector, the average nearest
* neighbor degree in the function of vertex degree is stored
* here. The first (zeroth) element is for degree one vertices,
* etc. Supply a NULL pointer here if you don't want to calculate
* this.
* \param weights Optional edge weights. Supply a null pointer here
* for the non-weighted version. The weighted version computes
* a weighted average of the neighbor degrees, i.e.
*
* k_nn_i = 1/s_i sum_j w_ij k_j
*
* where s_i is the sum of the weights, the sum runs over
* the neighbors as indicated by \c mode (with appropriate weights)
* and k_j is the degree, specified by \c neighbor_degree_mode.
* \return Error code.
*
* Time complexity: O(|V|+|E|), linear in the number of vertices and
* edges.
*
* \example examples/simple/igraph_knn.c
*/
int igraph_avg_nearest_neighbor_degree(const igraph_t *graph,
igraph_vs_t vids,
igraph_neimode_t mode,
igraph_neimode_t neighbor_degree_mode,
igraph_vector_t *knn,
igraph_vector_t *knnk,
const igraph_vector_t *weights) {
long int no_of_nodes = igraph_vcount(graph);
igraph_vector_t neis;
long int i, j, no_vids;
igraph_vit_t vit;
igraph_vector_t my_knn_v, *my_knn = knn;
igraph_vector_t deg;
igraph_integer_t maxdeg;
igraph_vector_t deghist;
igraph_real_t mynan = IGRAPH_NAN;
igraph_bool_t simple;
IGRAPH_CHECK(igraph_is_simple(graph, &simple));
if (!simple) {
IGRAPH_ERROR("Average nearest neighbor degree works only with "
"simple graphs", IGRAPH_EINVAL);
}
if (weights) {
return igraph_i_avg_nearest_neighbor_degree_weighted(graph, vids,
mode, neighbor_degree_mode, knn, knnk, weights);
}
IGRAPH_CHECK(igraph_vit_create(graph, vids, &vit));
IGRAPH_FINALLY(igraph_vit_destroy, &vit);
no_vids = IGRAPH_VIT_SIZE(vit);
if (!knn) {
IGRAPH_VECTOR_INIT_FINALLY(&my_knn_v, no_vids);
my_knn = &my_knn_v;
} else {
IGRAPH_CHECK(igraph_vector_resize(knn, no_vids));
}
IGRAPH_VECTOR_INIT_FINALLY(°, no_of_nodes);
IGRAPH_CHECK(igraph_degree(graph, °, igraph_vss_all(),
neighbor_degree_mode, IGRAPH_LOOPS));
igraph_maxdegree(graph, &maxdeg, igraph_vss_all(), mode, IGRAPH_LOOPS);
IGRAPH_VECTOR_INIT_FINALLY(&neis, maxdeg);
igraph_vector_resize(&neis, 0);
if (knnk) {
IGRAPH_CHECK(igraph_vector_resize(knnk, (long int)maxdeg));
igraph_vector_null(knnk);
IGRAPH_VECTOR_INIT_FINALLY(°hist, (long int)maxdeg);
}
for (i = 0; !IGRAPH_VIT_END(vit); IGRAPH_VIT_NEXT(vit), i++) {
igraph_real_t sum = 0.0;
long int v = IGRAPH_VIT_GET(vit);
long int nv;
IGRAPH_CHECK(igraph_neighbors(graph, &neis, (igraph_integer_t) v, mode));
nv = igraph_vector_size(&neis);
for (j = 0; j < nv; j++) {
long int nei = (long int) VECTOR(neis)[j];
sum += VECTOR(deg)[nei];
}
if (nv != 0) {
VECTOR(*my_knn)[i] = sum / nv;
} else {
VECTOR(*my_knn)[i] = mynan;
}
if (knnk && nv > 0) {
VECTOR(*knnk)[nv - 1] += VECTOR(*my_knn)[i];
VECTOR(deghist)[nv - 1] += 1;
}
}
if (knnk) {
for (i = 0; i < maxdeg; i++) {
long int dh = (long int) VECTOR(deghist)[i];
if (dh != 0) {
VECTOR(*knnk)[i] /= dh;
} else {
VECTOR(*knnk)[i] = mynan;
}
}
igraph_vector_destroy(°hist);
IGRAPH_FINALLY_CLEAN(1);
}
igraph_vector_destroy(&neis);
igraph_vector_destroy(°);
igraph_vit_destroy(&vit);
IGRAPH_FINALLY_CLEAN(3);
if (!knn) {
igraph_vector_destroy(&my_knn_v);
IGRAPH_FINALLY_CLEAN(1);
}
return 0;
}
/**
* \function igraph_strength
* Strength of the vertices, weighted vertex degree in other words.
*
* In a weighted network the strength of a vertex is the sum of the
* weights of all incident edges. In a non-weighted network this is
* exactly the vertex degree.
* \param graph The input graph.
* \param res Pointer to an initialized vector, the result is stored
* here. It will be resized as needed.
* \param vids The vertices for which the calculation is performed.
* \param mode Gives whether to count only outgoing (\c IGRAPH_OUT),
* incoming (\c IGRAPH_IN) edges or both (\c IGRAPH_ALL).
* \param loops A logical scalar, whether to count loop edges as well.
* \param weights A vector giving the edge weights. If this is a NULL
* pointer, then \ref igraph_degree() is called to perform the
* calculation.
* \return Error code.
*
* Time complexity: O(|V|+|E|), linear in the number vertices and
* edges.
*
* \sa \ref igraph_degree() for the traditional, non-weighted version.
*/
int igraph_strength(const igraph_t *graph, igraph_vector_t *res,
const igraph_vs_t vids, igraph_neimode_t mode,
igraph_bool_t loops, const igraph_vector_t *weights) {
long int no_of_nodes = igraph_vcount(graph);
igraph_vit_t vit;
long int no_vids;
igraph_vector_t neis;
long int i;
if (!weights) {
return igraph_degree(graph, res, vids, mode, loops);
}
if (igraph_vector_size(weights) != igraph_ecount(graph)) {
IGRAPH_ERROR("Invalid weight vector length", IGRAPH_EINVAL);
}
IGRAPH_CHECK(igraph_vit_create(graph, vids, &vit));
IGRAPH_FINALLY(igraph_vit_destroy, &vit);
no_vids = IGRAPH_VIT_SIZE(vit);
IGRAPH_VECTOR_INIT_FINALLY(&neis, 0);
IGRAPH_CHECK(igraph_vector_reserve(&neis, no_of_nodes));
IGRAPH_CHECK(igraph_vector_resize(res, no_vids));
igraph_vector_null(res);
if (loops) {
for (i = 0; !IGRAPH_VIT_END(vit); IGRAPH_VIT_NEXT(vit), i++) {
long int vid = IGRAPH_VIT_GET(vit);
long int j, n;
IGRAPH_CHECK(igraph_incident(graph, &neis, (igraph_integer_t) vid, mode));
n = igraph_vector_size(&neis);
for (j = 0; j < n; j++) {
long int edge = (long int) VECTOR(neis)[j];
VECTOR(*res)[i] += VECTOR(*weights)[edge];
}
}
} else {
for (i = 0; !IGRAPH_VIT_END(vit); IGRAPH_VIT_NEXT(vit), i++) {
long int vid = IGRAPH_VIT_GET(vit);
long int j, n;
IGRAPH_CHECK(igraph_incident(graph, &neis, (igraph_integer_t) vid, mode));
n = igraph_vector_size(&neis);
for (j = 0; j < n; j++) {
long int edge = (long int) VECTOR(neis)[j];
long int from = IGRAPH_FROM(graph, edge);
long int to = IGRAPH_TO(graph, edge);
if (from != to) {
VECTOR(*res)[i] += VECTOR(*weights)[edge];
}
}
}
}
igraph_vit_destroy(&vit);
igraph_vector_destroy(&neis);
IGRAPH_FINALLY_CLEAN(2);
return 0;
}
/**
* \function igraph_diameter_dijkstra
* Weighted diameter using Dijkstra's algorithm, non-negative weights only.
*
* The diameter of a graph is its longest geodesic. I.e. the
* (weighted) shortest path is calculated for all pairs of vertices
* and the longest one is the diameter.
* \param graph The input graph, can be directed or undirected.
* \param pres Pointer to a real number, 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 \c IGRAPH_INFINITY is
* returned.
* \return Error code.
*
* Time complexity: O(|V||E|*log|E|), |V| is the number of vertices,
* |E| is the number of edges.
*/
int igraph_diameter_dijkstra(const igraph_t *graph,
const igraph_vector_t *weights,
igraph_real_t *pres,
igraph_integer_t *pfrom,
igraph_integer_t *pto,
igraph_vector_t *path,
igraph_bool_t directed,
igraph_bool_t unconn) {
/* Implementation details. This is the basic Dijkstra algorithm,
with a binary heap. The heap is indexed, i.e. it stores not only
the distances, but also which vertex they belong to.
From now on we use a 2-way heap, so the distances can be queried
directly from the heap.
Dirty tricks:
- the opposite of the distance is stored in the heap, as it is a
maximum heap and we need a minimum heap.
- we don't use IGRAPH_INFINITY during the computation, as IGRAPH_FINITE()
might involve a function call and we want to spare that. -1 will denote
infinity instead.
*/
long int no_of_nodes = igraph_vcount(graph);
long int no_of_edges = igraph_ecount(graph);
igraph_2wheap_t Q;
igraph_inclist_t inclist;
long int source, j;
igraph_neimode_t dirmode = directed ? IGRAPH_OUT : IGRAPH_ALL;
long int from = -1, to = -1;
igraph_real_t res = 0;
long int nodes_reached = 0;
if (!weights) {
igraph_integer_t diameter;
IGRAPH_CHECK(igraph_diameter(graph, &diameter, pfrom, pto, path, directed, unconn));
if (pres) {
*pres = diameter;
}
return IGRAPH_SUCCESS;
}
if (weights && igraph_vector_size(weights) != no_of_edges) {
IGRAPH_ERROR("Invalid weight vector length", IGRAPH_EINVAL);
}
if (igraph_vector_min(weights) < 0) {
IGRAPH_ERROR("Weight vector must be non-negative", IGRAPH_EINVAL);
}
IGRAPH_CHECK(igraph_2wheap_init(&Q, no_of_nodes));
IGRAPH_FINALLY(igraph_2wheap_destroy, &Q);
IGRAPH_CHECK(igraph_inclist_init(graph, &inclist, dirmode));
IGRAPH_FINALLY(igraph_inclist_destroy, &inclist);
for (source = 0; source < no_of_nodes; source++) {
IGRAPH_PROGRESS("Weighted diameter: ", source * 100.0 / no_of_nodes, NULL);
IGRAPH_ALLOW_INTERRUPTION();
igraph_2wheap_clear(&Q);
igraph_2wheap_push_with_index(&Q, source, -1.0);
nodes_reached = 0.0;
while (!igraph_2wheap_empty(&Q)) {
long int minnei = igraph_2wheap_max_index(&Q);
igraph_real_t mindist = -igraph_2wheap_deactivate_max(&Q);
igraph_vector_int_t *neis;
long int nlen;
if (mindist > res) {
res = mindist; from = source; to = minnei;
}
nodes_reached++;
/* Now check all neighbors of 'minnei' for a shorter path */
neis = igraph_inclist_get(&inclist, minnei);
nlen = igraph_vector_int_size(neis);
for (j = 0; j < nlen; j++) {
long int edge = (long int) VECTOR(*neis)[j];
long int tto = IGRAPH_OTHER(graph, edge, minnei);
igraph_real_t altdist = mindist + VECTOR(*weights)[edge];
igraph_bool_t active = igraph_2wheap_has_active(&Q, tto);
igraph_bool_t has = igraph_2wheap_has_elem(&Q, tto);
igraph_real_t curdist = active ? -igraph_2wheap_get(&Q, tto) : 0.0;
if (!has) {
/* First finite distance */
IGRAPH_CHECK(igraph_2wheap_push_with_index(&Q, tto, -altdist));
} else if (altdist < curdist) {
/* A shorter path */
IGRAPH_CHECK(igraph_2wheap_modify(&Q, tto, -altdist));
}
}
} /* !igraph_2wheap_empty(&Q) */
/* not connected, return infinity */
if (nodes_reached != no_of_nodes && !unconn) {
res = IGRAPH_INFINITY;
from = to = -1;
break;
}
} /* source < no_of_nodes */
/* Compensate for the +1 that we have added to distances */
res -= 1;
igraph_inclist_destroy(&inclist);
igraph_2wheap_destroy(&Q);
IGRAPH_FINALLY_CLEAN(2);
IGRAPH_PROGRESS("Weighted diameter: ", 100.0, NULL);
if (pres) {
*pres = res;
}
if (pfrom) {
*pfrom = (igraph_integer_t) from;
}
if (pto) {
*pto = (igraph_integer_t) to;
}
if (path) {
if (!igraph_finite(res)) {
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_dijkstra(graph,
/*vertices=*/ &tmpptr, /*edges=*/ 0,
(igraph_integer_t) from,
igraph_vss_1((igraph_integer_t) to),
weights, dirmode, /*predecessors=*/ 0,
/*inbound_edges=*/ 0));
igraph_vector_ptr_destroy(&tmpptr);
IGRAPH_FINALLY_CLEAN(1);
}
}
return 0;
}
/**
* \function igraph_sort_vertex_ids_by_degree
* \brief Calculate a list of vertex ids sorted by degree of the corresponding vertex.
*
* The list of vertex ids is returned in a vector that is sorted
* in ascending or descending order of vertex degree.
*
* \param graph The input graph.
* \param outvids Pointer to an initialized vector that will be
* resized and will contain the ordered vertex ids.
* \param vids Input vertex selector of vertex ids to include in
* calculation.
* \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.
* \param order Specifies whether the ordering should be ascending
* (\c IGRAPH_ASCENDING) or descending (\c IGRAPH_DESCENDING).
* \param only_indices If true, then return a sorted list of indices
* into a vector corresponding to \c vids, rather than a list
* of vertex ids. This parameter is ignored if \c vids is set
* to all vertices via igraph_vs_all() or igraph_vss_all(),
* because in this case the indices and vertex ids are the
* same.
* \return Error code:
* \c IGRAPH_EINVVID: invalid vertex id.
* \c IGRAPH_EINVMODE: invalid mode argument.
*
*/
int igraph_sort_vertex_ids_by_degree(const igraph_t *graph,
igraph_vector_t *outvids,
igraph_vs_t vids,
igraph_neimode_t mode,
igraph_bool_t loops,
igraph_order_t order,
igraph_bool_t only_indices) {
long int i;
igraph_vector_t degrees, vs_vec;
IGRAPH_VECTOR_INIT_FINALLY(°rees, 0);
IGRAPH_CHECK(igraph_degree(graph, °rees, vids, mode, loops));
IGRAPH_CHECK((int) igraph_vector_qsort_ind(°rees, outvids,
order == IGRAPH_DESCENDING));
if (only_indices || igraph_vs_is_all(&vids) ) {
igraph_vector_destroy(°rees);
IGRAPH_FINALLY_CLEAN(1);
} else {
IGRAPH_VECTOR_INIT_FINALLY(&vs_vec, 0);
IGRAPH_CHECK(igraph_vs_as_vector(graph, vids, &vs_vec));
for (i = 0; i < igraph_vector_size(outvids); i++) {
VECTOR(*outvids)[i] = VECTOR(vs_vec)[(long int)VECTOR(*outvids)[i]];
}
igraph_vector_destroy(&vs_vec);
igraph_vector_destroy(°rees);
IGRAPH_FINALLY_CLEAN(2);
}
return 0;
}
/**
* \function igraph_contract_vertices
* Replace multiple vertices with a single one.
*
* This function creates a new graph, by merging several
* vertices into one. The vertices in the new graph correspond
* to sets of vertices in the input graph.
* \param graph The input graph, it can be directed or
* undirected.
* \param mapping A vector giving the mapping. For each
* vertex in the original graph, it should contain
* its id in the new graph.
* \param vertex_comb What to do with the vertex attributes.
* See the igraph manual section about attributes for
* details.
* \return Error code.
*
* Time complexity: O(|V|+|E|), linear in the number
* or vertices plus edges.
*/
int igraph_contract_vertices(igraph_t *graph,
const igraph_vector_t *mapping,
const igraph_attribute_combination_t
*vertex_comb) {
igraph_vector_t edges;
long int no_of_nodes = igraph_vcount(graph);
long int no_of_edges = igraph_ecount(graph);
igraph_bool_t vattr = vertex_comb && igraph_has_attribute_table();
igraph_t res;
long int e, last = -1;
long int no_new_vertices;
if (igraph_vector_size(mapping) != no_of_nodes) {
IGRAPH_ERROR("Invalid mapping vector length",
IGRAPH_EINVAL);
}
IGRAPH_VECTOR_INIT_FINALLY(&edges, 0);
IGRAPH_CHECK(igraph_vector_reserve(&edges, no_of_edges * 2));
if (no_of_nodes > 0) {
last = (long int) igraph_vector_max(mapping);
}
for (e = 0; e < no_of_edges; e++) {
long int from = IGRAPH_FROM(graph, e);
long int to = IGRAPH_TO(graph, e);
long int nfrom = (long int) VECTOR(*mapping)[from];
long int nto = (long int) VECTOR(*mapping)[to];
igraph_vector_push_back(&edges, nfrom);
igraph_vector_push_back(&edges, nto);
if (nfrom > last) {
last = nfrom;
}
if (nto > last) {
last = nto;
}
}
no_new_vertices = last + 1;
IGRAPH_CHECK(igraph_create(&res, &edges, (igraph_integer_t) no_new_vertices,
igraph_is_directed(graph)));
igraph_vector_destroy(&edges);
IGRAPH_FINALLY_CLEAN(1);
IGRAPH_FINALLY(igraph_destroy, &res);
IGRAPH_I_ATTRIBUTE_DESTROY(&res);
IGRAPH_I_ATTRIBUTE_COPY(&res, graph, /*graph=*/ 1,
/*vertex=*/ 0, /*edge=*/ 1);
if (vattr) {
long int i;
igraph_vector_ptr_t merges;
igraph_vector_t sizes;
igraph_vector_t *vecs;
vecs = igraph_Calloc(no_new_vertices, igraph_vector_t);
if (!vecs) {
IGRAPH_ERROR("Cannot combine attributes while contracting"
" vertices", IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(igraph_free, vecs);
IGRAPH_CHECK(igraph_vector_ptr_init(&merges, no_new_vertices));
IGRAPH_FINALLY(igraph_i_simplify_free, &merges);
IGRAPH_VECTOR_INIT_FINALLY(&sizes, no_new_vertices);
for (i = 0; i < no_of_nodes; i++) {
long int to = (long int) VECTOR(*mapping)[i];
VECTOR(sizes)[to] += 1;
}
for (i = 0; i < no_new_vertices; i++) {
igraph_vector_t *v = &vecs[i];
IGRAPH_CHECK(igraph_vector_init(v, (long int) VECTOR(sizes)[i]));
igraph_vector_clear(v);
VECTOR(merges)[i] = v;
}
for (i = 0; i < no_of_nodes; i++) {
long int to = (long int) VECTOR(*mapping)[i];
igraph_vector_t *v = &vecs[to];
igraph_vector_push_back(v, i);
}
IGRAPH_CHECK(igraph_i_attribute_combine_vertices(graph, &res,
&merges,
vertex_comb));
igraph_vector_destroy(&sizes);
igraph_i_simplify_free(&merges);
igraph_free(vecs);
IGRAPH_FINALLY_CLEAN(3);
}
IGRAPH_FINALLY_CLEAN(1);
igraph_destroy(graph);
*graph = res;
return 0;
}
/* Create the transitive closure of a tree graph.
This is fairly simple, we just collect all ancestors of a vertex
using a depth-first search.
*/
int igraph_transitive_closure_dag(const igraph_t *graph,
igraph_t *closure) {
long int no_of_nodes = igraph_vcount(graph);
igraph_vector_t deg;
igraph_vector_t new_edges;
igraph_vector_t ancestors;
long int root;
igraph_vector_t neighbors;
igraph_stack_t path;
igraph_vector_bool_t done;
if (!igraph_is_directed(graph)) {
IGRAPH_ERROR("Tree transitive closure of a directed graph",
IGRAPH_EINVAL);
}
IGRAPH_VECTOR_INIT_FINALLY(&new_edges, 0);
IGRAPH_VECTOR_INIT_FINALLY(°, no_of_nodes);
IGRAPH_VECTOR_INIT_FINALLY(&ancestors, 0);
IGRAPH_VECTOR_INIT_FINALLY(&neighbors, 0);
IGRAPH_CHECK(igraph_stack_init(&path, 0));
IGRAPH_FINALLY(igraph_stack_destroy, &path);
IGRAPH_CHECK(igraph_vector_bool_init(&done, no_of_nodes));
IGRAPH_FINALLY(igraph_vector_bool_destroy, &done);
IGRAPH_CHECK(igraph_degree(graph, °, igraph_vss_all(),
IGRAPH_OUT, IGRAPH_LOOPS));
#define STAR (-1)
for (root = 0; root < no_of_nodes; root++) {
if (VECTOR(deg)[root] != 0) {
continue;
}
IGRAPH_CHECK(igraph_stack_push(&path, root));
while (!igraph_stack_empty(&path)) {
long int node = (long int) igraph_stack_top(&path);
if (node == STAR) {
/* Leaving a node */
long int j, n;
igraph_stack_pop(&path);
node = (long int) igraph_stack_pop(&path);
if (!VECTOR(done)[node]) {
igraph_vector_pop_back(&ancestors);
VECTOR(done)[node] = 1;
}
n = igraph_vector_size(&ancestors);
for (j = 0; j < n; j++) {
IGRAPH_CHECK(igraph_vector_push_back(&new_edges, node));
IGRAPH_CHECK(igraph_vector_push_back(&new_edges,
VECTOR(ancestors)[j]));
}
} else {
/* Getting into a node */
long int n, j;
if (!VECTOR(done)[node]) {
IGRAPH_CHECK(igraph_vector_push_back(&ancestors, node));
}
IGRAPH_CHECK(igraph_neighbors(graph, &neighbors,
(igraph_integer_t) node, IGRAPH_IN));
n = igraph_vector_size(&neighbors);
IGRAPH_CHECK(igraph_stack_push(&path, STAR));
for (j = 0; j < n; j++) {
long int nei = (long int) VECTOR(neighbors)[j];
IGRAPH_CHECK(igraph_stack_push(&path, nei));
}
}
}
}
#undef STAR
igraph_vector_bool_destroy(&done);
igraph_stack_destroy(&path);
igraph_vector_destroy(&neighbors);
igraph_vector_destroy(&ancestors);
igraph_vector_destroy(°);
IGRAPH_FINALLY_CLEAN(5);
IGRAPH_CHECK(igraph_create(closure, &new_edges, (igraph_integer_t)no_of_nodes,
IGRAPH_DIRECTED));
igraph_vector_destroy(&new_edges);
IGRAPH_FINALLY_CLEAN(1);
return 0;
}
/**
* \function igraph_diversity
* Structural diversity index of the vertices
*
* This measure was defined in Nathan Eagle, Michael Macy and Rob
* Claxton: Network Diversity and Economic Development, Science 328,
* 1029--1031, 2010.
*
* </para><para>
* It is simply the (normalized) Shannon entropy of the
* incident edges' weights. D(i)=H(i)/log(k[i]), and
* H(i) = -sum(p[i,j] log(p[i,j]), j=1..k[i]),
* where p[i,j]=w[i,j]/sum(w[i,l], l=1..k[i]), k[i] is the (total)
* degree of vertex i, and w[i,j] is the weight of the edge(s) between
* vertex i and j.
* \param graph The input graph, edge directions are ignored.
* \param weights The edge weights, in the order of the edge ids, must
* have appropriate length.
* \param res An initialized vector, the results are stored here.
* \param vids Vector with the vertex ids for which to calculate the
* measure.
* \return Error code.
*
* Time complexity: O(|V|+|E|), linear.
*
*/
int igraph_diversity(igraph_t *graph, const igraph_vector_t *weights,
igraph_vector_t *res, const igraph_vs_t vids) {
int no_of_nodes = igraph_vcount(graph);
int no_of_edges = igraph_ecount(graph);
igraph_vector_t incident;
igraph_vit_t vit;
igraph_real_t s, ent, w;
int i, j, k;
if (!weights) {
IGRAPH_ERROR("Edge weights must be given", IGRAPH_EINVAL);
}
if (igraph_vector_size(weights) != no_of_edges) {
IGRAPH_ERROR("Invalid edge weight vector length", IGRAPH_EINVAL);
}
IGRAPH_VECTOR_INIT_FINALLY(&incident, 10);
if (igraph_vs_is_all(&vids)) {
IGRAPH_CHECK(igraph_vector_resize(res, no_of_nodes));
for (i = 0; i < no_of_nodes; i++) {
s = ent = 0.0;
IGRAPH_CHECK(igraph_incident(graph, &incident, i, /*mode=*/ IGRAPH_ALL));
for (j = 0, k = (int) igraph_vector_size(&incident); j < k; j++) {
w = VECTOR(*weights)[(long int)VECTOR(incident)[j]];
s += w;
ent += (w * log(w));
}
VECTOR(*res)[i] = (log(s) - ent / s) / log(k);
}
} else {
IGRAPH_CHECK(igraph_vector_resize(res, 0));
IGRAPH_CHECK(igraph_vit_create(graph, vids, &vit));
IGRAPH_FINALLY(igraph_vit_destroy, &vit);
for (IGRAPH_VIT_RESET(vit), i = 0;
!IGRAPH_VIT_END(vit);
IGRAPH_VIT_NEXT(vit), i++) {
long int v = IGRAPH_VIT_GET(vit);
s = ent = 0.0;
IGRAPH_CHECK(igraph_incident(graph, &incident, (igraph_integer_t) v,
/*mode=*/ IGRAPH_ALL));
for (j = 0, k = (int) igraph_vector_size(&incident); j < k; j++) {
w = VECTOR(*weights)[(long int)VECTOR(incident)[j]];
s += w;
ent += (w * log(w));
}
IGRAPH_CHECK(igraph_vector_push_back(res, (log(s) - ent / s) / log(k)));
}
igraph_vit_destroy(&vit);
IGRAPH_FINALLY_CLEAN(1);
}
igraph_vector_destroy(&incident);
IGRAPH_FINALLY_CLEAN(1);
return 0;
}
#define SUCCEED { \
if (res) { \
*res = 1; \
} \
return IGRAPH_SUCCESS; \
}
#define FAIL { \
if (res) { \
*res = 0; \
} \
return IGRAPH_SUCCESS; \
}
/**
* \function igraph_is_degree_sequence
* Determines whether a degree sequence is valid.
*
* A sequence of n integers is a valid degree sequence if there exists some
* graph where the degree of the i-th vertex is equal to the i-th element of the
* sequence. Note that the graph may contain multiple or loop edges; if you are
* interested in whether the degrees of some \em simple graph may realize the
* given sequence, use \ref igraph_is_graphical_degree_sequence.
*
* </para><para>
* In particular, the function checks whether all the degrees are non-negative.
* For undirected graphs, it also checks whether the sum of degrees is even.
* For directed graphs, the function checks whether the lengths of the two
* degree vectors are equal and whether their sums are also equal. These are
* known sufficient and necessary conditions for a degree sequence to be
* valid.
*
* \param out_degrees an integer vector specifying the degree sequence for
* undirected graphs or the out-degree sequence for directed graphs.
* \param in_degrees an integer vector specifying the in-degrees of the
* vertices for directed graphs. For undirected graphs, this must be null.
* \param res pointer to a boolean variable, the result will be stored here
* \return Error code.
*
* Time complexity: O(n), where n is the length of the degree sequence.
*/
int igraph_is_degree_sequence(const igraph_vector_t *out_degrees,
const igraph_vector_t *in_degrees, igraph_bool_t *res) {
/* degrees must be non-negative */
if (igraph_vector_any_smaller(out_degrees, 0)) {
FAIL;
}
if (in_degrees && igraph_vector_any_smaller(in_degrees, 0)) {
FAIL;
}
if (in_degrees == 0) {
/* sum of degrees must be even */
if (((long int)igraph_vector_sum(out_degrees) % 2) != 0) {
FAIL;
}
} else {
/* length of the two degree vectors must be equal */
if (igraph_vector_size(out_degrees) != igraph_vector_size(in_degrees)) {
FAIL;
}
/* sum of in-degrees must be equal to sum of out-degrees */
if (igraph_vector_sum(out_degrees) != igraph_vector_sum(in_degrees)) {
FAIL;
}
}
SUCCEED;
return 0;
}
int igraph_i_is_graphical_degree_sequence_undirected(
const igraph_vector_t *degrees, igraph_bool_t *res);
int igraph_i_is_graphical_degree_sequence_directed(
const igraph_vector_t *out_degrees, const igraph_vector_t *in_degrees,
igraph_bool_t *res);
/**
* \function igraph_is_graphical_degree_sequence
* Determines whether a sequence of integers can be a degree sequence of some
* simple graph.
*
* </para><para>
* References:
*
* </para><para>
* Hakimi SL: On the realizability of a set of integers as degrees of the
* vertices of a simple graph. J SIAM Appl Math 10:496-506, 1962.
*
* </para><para>
* PL Erdos, I Miklos and Z Toroczkai: A simple Havel-Hakimi type algorithm
* to realize graphical degree sequences of directed graphs. The Electronic
* Journal of Combinatorics 17(1):R66, 2010.
*
* </para><para>
* Z Kiraly: Recognizing graphic degree sequences and generating all
* realizations. TR-2011-11, Egervary Research Group, H-1117, Budapest,
* Hungary. ISSN 1587-4451, 2012.
*
* \param out_degrees an integer vector specifying the degree sequence for
* undirected graphs or the out-degree sequence for directed graphs.
* \param in_degrees an integer vector specifying the in-degrees of the
* vertices for directed graphs. For undirected graphs, this must be null.
* \param res pointer to a boolean variable, the result will be stored here
* \return Error code.
*
* Time complexity: O(n log n) for undirected graphs, O(n^2) for directed
* graphs, where n is the length of the degree sequence.
*/
int igraph_is_graphical_degree_sequence(const igraph_vector_t *out_degrees,
const igraph_vector_t *in_degrees, igraph_bool_t *res) {
IGRAPH_CHECK(igraph_is_degree_sequence(out_degrees, in_degrees, res));
if (!*res) {
FAIL;
}
if (igraph_vector_size(out_degrees) == 0) {
SUCCEED;
}
if (in_degrees == 0) {
return igraph_i_is_graphical_degree_sequence_undirected(out_degrees, res);
} else {
return igraph_i_is_graphical_degree_sequence_directed(out_degrees, in_degrees, res);
}
}
int igraph_i_is_graphical_degree_sequence_undirected(
const igraph_vector_t *degrees, igraph_bool_t *res) {
igraph_vector_t work;
long int w, b, s, c, n, k;
IGRAPH_CHECK(igraph_vector_copy(&work, degrees));
IGRAPH_FINALLY(igraph_vector_destroy, &work);
igraph_vector_sort(&work);
/* This algorithm is outlined in TR-2011-11 of the Egervary Research Group,
* ISSN 1587-4451. The main loop of the algorithm is O(n) but it is dominated
* by an O(n log n) quicksort; this could in theory be brought down to
* O(n) with binsort but it's probably not worth the fuss.
*
* Variables names are mostly according to the technical report, apart from
* the degrees themselves. w and k are zero-based here; in the technical
* report they are 1-based */
*res = 1;
n = igraph_vector_size(&work);
w = n - 1; b = 0; s = 0; c = 0;
for (k = 0; k < n; k++) {
b += VECTOR(*degrees)[k];
c += w;
while (w > k && VECTOR(*degrees)[w] <= k + 1) {
s += VECTOR(*degrees)[w];
c -= (k + 1);
w--;
}
if (b > c + s) {
*res = 0;
break;
}
if (w == k) {
break;
}
}
igraph_vector_destroy(&work);
IGRAPH_FINALLY_CLEAN(1);
return 0;
}
typedef struct {
const igraph_vector_t* first;
const igraph_vector_t* second;
} igraph_i_qsort_dual_vector_cmp_data_t;
int igraph_i_qsort_dual_vector_cmp_desc(void* data, const void *p1, const void *p2) {
igraph_i_qsort_dual_vector_cmp_data_t* sort_data =
(igraph_i_qsort_dual_vector_cmp_data_t*)data;
long int index1 = *((long int*)p1);
long int index2 = *((long int*)p2);
if (VECTOR(*sort_data->first)[index1] < VECTOR(*sort_data->first)[index2]) {
return 1;
}
if (VECTOR(*sort_data->first)[index1] > VECTOR(*sort_data->first)[index2]) {
return -1;
}
if (VECTOR(*sort_data->second)[index1] < VECTOR(*sort_data->second)[index2]) {
return 1;
}
if (VECTOR(*sort_data->second)[index1] > VECTOR(*sort_data->second)[index2]) {
return -1;
}
return 0;
}
int igraph_i_is_graphical_degree_sequence_directed(
const igraph_vector_t *out_degrees, const igraph_vector_t *in_degrees,
igraph_bool_t *res) {
igraph_vector_long_t index_array;
long int i, j, vcount, lhs, rhs;
igraph_i_qsort_dual_vector_cmp_data_t sort_data;
/* Create an index vector that sorts the vertices by decreasing in-degree */
vcount = igraph_vector_size(out_degrees);
IGRAPH_CHECK(igraph_vector_long_init_seq(&index_array, 0, vcount - 1));
IGRAPH_FINALLY(igraph_vector_long_destroy, &index_array);
/* Set up the auxiliary struct for sorting */
sort_data.first = in_degrees;
sort_data.second = out_degrees;
/* Sort the index vector */
igraph_qsort_r(VECTOR(index_array), vcount, sizeof(long int), &sort_data,
igraph_i_qsort_dual_vector_cmp_desc);
/* Be optimistic, then check whether the Fulkerson–Chen–Anstee condition
* holds for every k. In particular, for every k in [0; n), it must be true
* that:
*
* \sum_{i=0}^k indegree[i] <=
* \sum_{i=0}^k min(outdegree[i], k) +
* \sum_{i=k+1}^{n-1} min(outdegree[i], k + 1)
*/
#define INDEGREE(x) (VECTOR(*in_degrees)[VECTOR(index_array)[x]])
#define OUTDEGREE(x) (VECTOR(*out_degrees)[VECTOR(index_array)[x]])
*res = 1;
lhs = 0;
for (i = 0; i < vcount; i++) {
lhs += INDEGREE(i);
/* It is enough to check for indexes where the in-degree is about to
* decrease in the next step; see "Stronger condition" in the Wikipedia
* entry for the Fulkerson-Chen-Anstee condition */
if (i != vcount - 1 && INDEGREE(i) == INDEGREE(i + 1)) {
continue;
}
rhs = 0;
for (j = 0; j <= i; j++) {
rhs += OUTDEGREE(j) < i ? OUTDEGREE(j) : i;
}
for (j = i + 1; j < vcount; j++) {
rhs += OUTDEGREE(j) < (i + 1) ? OUTDEGREE(j) : (i + 1);
}
if (lhs > rhs) {
*res = 0;
break;
}
}
#undef INDEGREE
#undef OUTDEGREE
igraph_vector_long_destroy(&index_array);
IGRAPH_FINALLY_CLEAN(1);
return IGRAPH_SUCCESS;
}
#undef SUCCEED
#undef FAIL
/* igraph_is_tree -- check if a graph is a tree */
/* count the number of vertices reachable from the root */
static int igraph_i_is_tree_visitor(igraph_integer_t root, const igraph_adjlist_t *al, igraph_integer_t *visited_count) {
igraph_stack_int_t stack;
igraph_vector_bool_t visited;
long i;
IGRAPH_CHECK(igraph_vector_bool_init(&visited, igraph_adjlist_size(al)));
IGRAPH_FINALLY(igraph_vector_bool_destroy, &visited);
IGRAPH_CHECK(igraph_stack_int_init(&stack, 0));
IGRAPH_FINALLY(igraph_stack_int_destroy, &stack);
*visited_count = 0;
/* push the root into the stack */
IGRAPH_CHECK(igraph_stack_int_push(&stack, root));
while (! igraph_stack_int_empty(&stack)) {
igraph_integer_t u;
igraph_vector_int_t *neighbors;
long ncount;
/* take a vertex from the stack, mark it as visited */
u = igraph_stack_int_pop(&stack);
if (IGRAPH_LIKELY(! VECTOR(visited)[u])) {
VECTOR(visited)[u] = 1;
*visited_count += 1;
}
/* register all its yet-unvisited neighbours for future processing */
neighbors = igraph_adjlist_get(al, u);
ncount = igraph_vector_int_size(neighbors);
for (i = 0; i < ncount; ++i) {
igraph_integer_t v = VECTOR(*neighbors)[i];
if (! VECTOR(visited)[v]) {
IGRAPH_CHECK(igraph_stack_int_push(&stack, v));
}
}
}
igraph_stack_int_destroy(&stack);
igraph_vector_bool_destroy(&visited);
IGRAPH_FINALLY_CLEAN(2);
return IGRAPH_SUCCESS;
}
/**
* \ingroup structural
* \function igraph_is_tree
* \brief Decides whether the graph is a tree.
*
* An undirected graph is a tree if it is connected and has no cycles.
* </para><para>
*
* In the directed case, a possible additional requirement is that all
* edges are oriented away from a root (out-tree or arborescence) or all edges
* are oriented towards a root (in-tree or anti-arborescence).
* This test can be controlled using the \p mode parameter.
* </para><para>
*
* By convention, the null graph (i.e. the graph with no vertices) is considered not to be a tree.
*
* \param graph The graph object to analyze.
* \param res Pointer to a logical variable, the result will be stored
* here.
* \param root If not \c NULL, the root node will be stored here. When \p mode
* is \c IGRAPH_ALL or the graph is undirected, any vertex can be the root
* and \p root is set to 0 (the first vertex). When \p mode is \c IGRAPH_OUT
* or \c IGRAPH_IN, the root is set to the vertex with zero in- or out-degree,
* respectively.
* \param mode For a directed graph this specifies whether to test for an
* out-tree, an in-tree or ignore edge directions. The respective
* possible values are:
* \c IGRAPH_OUT, \c IGRAPH_IN, \c IGRAPH_ALL. This argument is
* ignored for undirected graphs.
* \return Error code:
* \c IGRAPH_EINVAL: invalid mode argument.
*
* Time complexity: At most O(|V|+|E|), the
* number of vertices plus the number of edges in the graph.
*
* \sa igraph_is_weakly_connected()
*
* \example examples/simple/igraph_tree.c
*/
int igraph_is_tree(const igraph_t *graph, igraph_bool_t *res, igraph_integer_t *root, igraph_neimode_t mode) {
igraph_adjlist_t al;
igraph_integer_t iroot = 0;
igraph_integer_t visited_count;
igraph_integer_t vcount, ecount;
vcount = igraph_vcount(graph);
ecount = igraph_ecount(graph);
/* A tree must have precisely vcount-1 edges. */
/* By convention, the zero-vertex graph will not be considered a tree. */
if (ecount != vcount - 1) {
*res = 0;
return IGRAPH_SUCCESS;
}
/* The single-vertex graph is a tree, provided it has no edges (checked in the previous if (..)) */
if (vcount == 1) {
*res = 1;
if (root) {
*root = 0;
}
return IGRAPH_SUCCESS;
}
/* For higher vertex counts we cannot short-circuit due to the possibility
* of loops or multi-edges even when the edge count is correct. */
/* Ignore mode for undirected graphs. */
if (! igraph_is_directed(graph)) {
mode = IGRAPH_ALL;
}
IGRAPH_CHECK(igraph_adjlist_init(graph, &al, mode));
IGRAPH_FINALLY(igraph_adjlist_destroy, &al);
/* The main algorithm:
* We find a root and check that all other vertices are reachable from it.
* We have already checked the number of edges, so with the additional
* reachability condition we can verify if the graph is a tree.
*
* For directed graphs, the root is the node with no incoming/outgoing
* connections, depending on 'mode'. For undirected, it is arbitrary, so
* we choose 0.
*/
*res = 1; /* assume success */
switch (mode) {
case IGRAPH_ALL:
iroot = 0;
break;
case IGRAPH_IN:
case IGRAPH_OUT: {
igraph_vector_t degree;
igraph_integer_t i;
IGRAPH_CHECK(igraph_vector_init(°ree, 0));
IGRAPH_FINALLY(igraph_vector_destroy, °ree);
IGRAPH_CHECK(igraph_degree(graph, °ree, igraph_vss_all(), mode == IGRAPH_IN ? IGRAPH_OUT : IGRAPH_IN, /* loops = */ 1));
for (i = 0; i < vcount; ++i)
if (VECTOR(degree)[i] == 0) {
break;
}
/* if no suitable root is found, the graph is not a tree */
if (i == vcount) {
*res = 0;
} else {
iroot = i;
}
igraph_vector_destroy(°ree);
IGRAPH_FINALLY_CLEAN(1);
}
break;
default:
IGRAPH_ERROR("Invalid mode", IGRAPH_EINVMODE);
}
/* if no suitable root was found, skip visting vertices */
if (*res) {
IGRAPH_CHECK(igraph_i_is_tree_visitor(iroot, &al, &visited_count));
*res = visited_count == vcount;
}
if (root) {
*root = iroot;
}
igraph_adjlist_destroy(&al);
IGRAPH_FINALLY_CLEAN(1);
return IGRAPH_SUCCESS;
}