/* -*- mode: C -*- */
/* vim:set ts=4 sw=4 sts=4 et: */
/*
IGraph R package.
Copyright (C) 2003-2014 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_layout.h"
#include "igraph_random.h"
#include "igraph_memory.h"
#include "igraph_iterators.h"
#include "igraph_interface.h"
#include "igraph_adjlist.h"
#include "igraph_progress.h"
#include "igraph_interrupt_internal.h"
#include "igraph_paths.h"
#include "igraph_structural.h"
#include "igraph_visitor.h"
#include "igraph_topology.h"
#include "igraph_components.h"
#include "igraph_types_internal.h"
#include "igraph_dqueue.h"
#include "igraph_arpack.h"
#include "igraph_blas.h"
#include "igraph_centrality.h"
#include "igraph_eigen.h"
#include "config.h"
#include <math.h>
#include "igraph_math.h"
#include <stdio.h> /* FIXME */
/**
* \section about_layouts
*
* <para>Layout generator functions (or at least most of them) try to place the
* vertices and edges of a graph on a 2D plane or in 3D space in a way
* which visually pleases the human eye.</para>
*
* <para>They take a graph object and a number of parameters as arguments
* and return an \type igraph_matrix_t, in which each row gives the
* coordinates of a vertex.</para>
*/
/**
* \ingroup layout
* \function igraph_layout_random
* \brief Places the vertices uniform randomly on a plane.
*
* \param graph Pointer to an initialized graph object.
* \param res Pointer to an initialized matrix object. This will
* contain the result and will be resized as needed.
* \return Error code. The current implementation always returns with
* success.
*
* Time complexity: O(|V|), the
* number of vertices.
*/
int igraph_layout_random(const igraph_t *graph, igraph_matrix_t *res) {
long int no_of_nodes = igraph_vcount(graph);
long int i;
IGRAPH_CHECK(igraph_matrix_resize(res, no_of_nodes, 2));
RNG_BEGIN();
for (i = 0; i < no_of_nodes; i++) {
MATRIX(*res, i, 0) = RNG_UNIF(-1, 1);
MATRIX(*res, i, 1) = RNG_UNIF(-1, 1);
}
RNG_END();
return 0;
}
/**
* \function igraph_layout_random_3d
* \brief Random layout in 3D
*
* \param graph The graph to place.
* \param res Pointer to an initialized matrix object. It will be
* resized to hold the result.
* \return Error code. The current implementation always returns with
* success.
*
* Added in version 0.2.</para><para>
*
* Time complexity: O(|V|), the number of vertices.
*/
int igraph_layout_random_3d(const igraph_t *graph, igraph_matrix_t *res) {
long int no_of_nodes = igraph_vcount(graph);
long int i;
IGRAPH_CHECK(igraph_matrix_resize(res, no_of_nodes, 3));
RNG_BEGIN();
for (i = 0; i < no_of_nodes; i++) {
MATRIX(*res, i, 0) = RNG_UNIF(-1, 1);
MATRIX(*res, i, 1) = RNG_UNIF(-1, 1);
MATRIX(*res, i, 2) = RNG_UNIF(-1, 1);
}
RNG_END();
return 0;
}
/**
* \ingroup layout
* \function igraph_layout_circle
* \brief Places the vertices uniformly on a circle, in the order of vertex ids.
*
* \param graph Pointer to an initialized graph object.
* \param res Pointer to an initialized matrix object. This will
* contain the result and will be resized as needed.
* \param order The order of the vertices on the circle. The vertices
* not included here, will be placed at (0,0). Supply
* \ref igraph_vss_all() here for all vertices, in the order of
* their vertex ids.
* \return Error code.
*
* Time complexity: O(|V|), the
* number of vertices.
*/
int igraph_layout_circle(const igraph_t *graph, igraph_matrix_t *res,
igraph_vs_t order) {
long int no_of_nodes = igraph_vcount(graph);
igraph_integer_t vs_size;
long int i;
igraph_vit_t vit;
IGRAPH_CHECK(igraph_vs_size(graph, &order, &vs_size));
IGRAPH_CHECK(igraph_matrix_resize(res, no_of_nodes, 2));
igraph_matrix_null(res);
igraph_vit_create(graph, order, &vit);
for (i = 0; !IGRAPH_VIT_END(vit); IGRAPH_VIT_NEXT(vit), i++) {
igraph_real_t phi = 2 * M_PI / vs_size * i;
int idx = IGRAPH_VIT_GET(vit);
MATRIX(*res, idx, 0) = cos(phi);
MATRIX(*res, idx, 1) = sin(phi);
}
igraph_vit_destroy(&vit);
return 0;
}
/**
* \function igraph_layout_star
* Generate a star-like layout
*
* \param graph The input graph.
* \param res Pointer to an initialized matrix object. This will
* contain the result and will be resized as needed.
* \param center The id of the vertex to put in the center.
* \param order A numeric vector giving the order of the vertices
* (including the center vertex!). If a null pointer, then the
* vertices are placed in increasing vertex id order.
* \return Error code.
*
* Time complexity: O(|V|), linear in the number of vertices.
*
* \sa \ref igraph_layout_circle() and other layout generators.
*/
int igraph_layout_star(const igraph_t *graph, igraph_matrix_t *res,
igraph_integer_t center, const igraph_vector_t *order) {
long int no_of_nodes = igraph_vcount(graph);
long int c = center;
long int i;
igraph_real_t step;
igraph_real_t phi;
if (order && igraph_vector_size(order) != no_of_nodes) {
IGRAPH_ERROR("Invalid order vector length", IGRAPH_EINVAL);
}
IGRAPH_CHECK(igraph_matrix_resize(res, no_of_nodes, 2));
if (no_of_nodes == 1) {
MATRIX(*res, 0, 0) = MATRIX(*res, 0, 1) = 0.0;
} else {
for (i = 0, step = 2 * M_PI / (no_of_nodes - 1), phi = 0;
i < no_of_nodes; i++) {
long int node = order ? (long int) VECTOR(*order)[i] : i;
if (node != c) {
MATRIX(*res, node, 0) = cos(phi);
MATRIX(*res, node, 1) = sin(phi);
phi += step;
} else {
MATRIX(*res, node, 0) = MATRIX(*res, node, 1) = 0.0;
}
}
}
return 0;
}
/**
* \function igraph_layout_sphere
* \brief Places vertices (more or less) uniformly on a sphere.
*
* </para><para>
* The algorithm was described in the following paper:
* Distributing many points on a sphere by E.B. Saff and
* A.B.J. Kuijlaars, \emb Mathematical Intelligencer \eme 19.1 (1997)
* 5--11.
*
* \param graph Pointer to an initialized graph object.
* \param res Pointer to an initialized matrix object. This will
* contain the result and will be resized as needed.
* \return Error code. The current implementation always returns with
* success.
*
* Added in version 0.2.</para><para>
*
* Time complexity: O(|V|), the number of vertices in the graph.
*/
int igraph_layout_sphere(const igraph_t *graph, igraph_matrix_t *res) {
long int no_of_nodes = igraph_vcount(graph);
long int i;
igraph_real_t h;
IGRAPH_CHECK(igraph_matrix_resize(res, no_of_nodes, 3));
if (no_of_nodes != 0) {
MATRIX(*res, 0, 0) = M_PI;
MATRIX(*res, 0, 1) = 0;
}
for (i = 1; i < no_of_nodes - 1; i++) {
h = -1 + 2 * i / (double)(no_of_nodes - 1);
MATRIX(*res, i, 0) = acos(h);
MATRIX(*res, i, 1) = fmod((MATRIX(*res, i - 1, 1) +
3.6 / sqrt(no_of_nodes * (1 - h * h))), 2 * M_PI);
IGRAPH_ALLOW_INTERRUPTION();
}
if (no_of_nodes >= 2) {
MATRIX(*res, no_of_nodes - 1, 0) = 0;
MATRIX(*res, no_of_nodes - 1, 1) = 0;
}
for (i = 0; i < no_of_nodes; i++) {
igraph_real_t x = cos(MATRIX(*res, i, 1)) * sin(MATRIX(*res, i, 0));
igraph_real_t y = sin(MATRIX(*res, i, 1)) * sin(MATRIX(*res, i, 0));
igraph_real_t z = cos(MATRIX(*res, i, 0));
MATRIX(*res, i, 0) = x;
MATRIX(*res, i, 1) = y;
MATRIX(*res, i, 2) = z;
IGRAPH_ALLOW_INTERRUPTION();
}
return 0;
}
/**
* \ingroup layout
* \function igraph_layout_grid
* \brief Places the vertices on a regular grid on the plane.
*
* \param graph Pointer to an initialized graph object.
* \param res Pointer to an initialized matrix object. This will
* contain the result and will be resized as needed.
* \param width The number of vertices in a single row of the grid.
* When zero or negative, the width of the grid will be the
* square root of the number of vertices, rounded up if needed.
* \return Error code. The current implementation always returns with
* success.
*
* Time complexity: O(|V|), the number of vertices.
*/
int igraph_layout_grid(const igraph_t *graph, igraph_matrix_t *res, long int width) {
long int i, no_of_nodes = igraph_vcount(graph);
igraph_real_t x, y;
IGRAPH_CHECK(igraph_matrix_resize(res, no_of_nodes, 2));
if (width <= 0) {
width = (long int) ceil(sqrt(no_of_nodes));
}
x = y = 0;
for (i = 0; i < no_of_nodes; i++) {
MATRIX(*res, i, 0) = x++;
MATRIX(*res, i, 1) = y;
if (x == width) {
x = 0; y++;
}
}
return 0;
}
/**
* \ingroup layout
* \function igraph_layout_grid_3d
* \brief Places the vertices on a regular grid in the 3D space.
*
* \param graph Pointer to an initialized graph object.
* \param res Pointer to an initialized matrix object. This will
* contain the result and will be resized as needed.
* \param width The number of vertices in a single row of the grid. When
* zero or negative, the width is determined automatically.
* \param height The number of vertices in a single column of the grid. When
* zero or negative, the height is determined automatically.
*
* \return Error code. The current implementation always returns with
* success.
*
* Time complexity: O(|V|), the number of vertices.
*/
int igraph_layout_grid_3d(const igraph_t *graph, igraph_matrix_t *res,
long int width, long int height) {
long int i, no_of_nodes = igraph_vcount(graph);
igraph_real_t x, y, z;
IGRAPH_CHECK(igraph_matrix_resize(res, no_of_nodes, 3));
if (width <= 0 && height <= 0) {
width = height = (long int) ceil(pow(no_of_nodes, 1.0 / 3));
} else if (width <= 0) {
width = (long int) ceil(sqrt(no_of_nodes / (double)height));
} else if (height <= 0) {
height = (long int) ceil(sqrt(no_of_nodes / (double)width));
}
x = y = z = 0;
for (i = 0; i < no_of_nodes; i++) {
MATRIX(*res, i, 0) = x++;
MATRIX(*res, i, 1) = y;
MATRIX(*res, i, 2) = z;
if (x == width) {
x = 0; y++;
if (y == height) {
y = 0; z++;
}
}
}
return 0;
}
int igraph_layout_springs(const igraph_t *graph, igraph_matrix_t *res,
igraph_real_t mass, igraph_real_t equil, igraph_real_t k,
igraph_real_t repeqdis, igraph_real_t kfr, igraph_bool_t repulse) {
IGRAPH_UNUSED(graph); IGRAPH_UNUSED(res); IGRAPH_UNUSED(mass);
IGRAPH_UNUSED(equil); IGRAPH_UNUSED(k); IGRAPH_UNUSED(repeqdis);
IGRAPH_UNUSED(kfr); IGRAPH_UNUSED(repulse);
IGRAPH_ERROR("Springs layout not implemented", IGRAPH_UNIMPLEMENTED);
/* TODO */
return 0;
}
void igraph_i_norm2d(igraph_real_t *x, igraph_real_t *y);
void igraph_i_norm2d(igraph_real_t *x, igraph_real_t *y) {
igraph_real_t len = sqrt((*x) * (*x) + (*y) * (*y));
if (len != 0) {
*x /= len;
*y /= len;
}
}
/**
* \function igraph_layout_lgl
* \brief Force based layout algorithm for large graphs.
*
* </para><para>
* This is a layout generator similar to the Large Graph Layout
* algorithm and program
* (http://lgl.sourceforge.net/). But unlike LGL, this
* version uses a Fruchterman-Reingold style simulated annealing
* algorithm for placing the vertices. The speedup is achieved by
* placing the vertices on a grid and calculating the repulsion only
* for vertices which are closer to each other than a limit.
*
* \param graph The (initialized) graph object to place.
* \param res Pointer to an initialized matrix object to hold the
* result. It will be resized if needed.
* \param maxit The maximum number of cooling iterations to perform
* for each layout step. A reasonable default is 150.
* \param maxdelta The maximum length of the move allowed for a vertex
* in a single iteration. A reasonable default is the number of
* vertices.
* \param area This parameter gives the area of the square on which
* the vertices will be placed. A reasonable default value is the
* number of vertices squared.
* \param coolexp The cooling exponent. A reasonable default value is
* 1.5.
* \param repulserad Determines the radius at which vertex-vertex
* repulsion cancels out attraction of adjacent vertices. A
* reasonable default value is \p area times the number of vertices.
* \param cellsize The size of the grid cells, one side of the
* square. A reasonable default value is the fourth root of
* \p area (or the square root of the number of vertices if \p area
* is also left at its default value).
* \param proot The root vertex, this is placed first, its neighbors
* in the first iteration, second neighbors in the second, etc. If
* negative then a random vertex is chosen.
* \return Error code.
*
* Added in version 0.2.</para><para>
*
* Time complexity: ideally O(dia*maxit*(|V|+|E|)), |V| is the number
* of vertices,
* dia is the diameter of the graph, worst case complexity is still
* O(dia*maxit*(|V|^2+|E|)), this is the case when all vertices happen to be
* in the same grid cell.
*/
int igraph_layout_lgl(const igraph_t *graph, igraph_matrix_t *res,
igraph_integer_t maxit, igraph_real_t maxdelta,
igraph_real_t area, igraph_real_t coolexp,
igraph_real_t repulserad, igraph_real_t cellsize,
igraph_integer_t proot) {
long int no_of_nodes = igraph_vcount(graph);
long int no_of_edges = igraph_ecount(graph);
igraph_t mst;
long int root;
long int no_of_layers, actlayer = 0;
igraph_vector_t vids;
igraph_vector_t layers;
igraph_vector_t parents;
igraph_vector_t edges;
igraph_2dgrid_t grid;
igraph_vector_t eids;
igraph_vector_t forcex;
igraph_vector_t forcey;
igraph_real_t frk = sqrt(area / no_of_nodes);
igraph_real_t H_n = 0;
IGRAPH_CHECK(igraph_minimum_spanning_tree_unweighted(graph, &mst));
IGRAPH_FINALLY(igraph_destroy, &mst);
IGRAPH_CHECK(igraph_matrix_resize(res, no_of_nodes, 2));
/* Determine the root vertex, random pick right now */
if (proot < 0) {
root = RNG_INTEGER(0, no_of_nodes - 1);
} else {
root = proot;
}
/* Assign the layers */
IGRAPH_VECTOR_INIT_FINALLY(&vids, 0);
IGRAPH_VECTOR_INIT_FINALLY(&layers, 0);
IGRAPH_VECTOR_INIT_FINALLY(&parents, 0);
IGRAPH_CHECK(igraph_i_bfs(&mst, (igraph_integer_t) root, IGRAPH_ALL, &vids,
&layers, &parents));
no_of_layers = igraph_vector_size(&layers) - 1;
/* We don't need the mst any more */
igraph_destroy(&mst);
igraph_empty(&mst, 0, IGRAPH_UNDIRECTED); /* to make finalization work */
IGRAPH_VECTOR_INIT_FINALLY(&edges, 0);
IGRAPH_CHECK(igraph_vector_reserve(&edges, no_of_edges));
IGRAPH_VECTOR_INIT_FINALLY(&eids, 0);
IGRAPH_VECTOR_INIT_FINALLY(&forcex, no_of_nodes);
IGRAPH_VECTOR_INIT_FINALLY(&forcey, no_of_nodes);
/* Place the vertices randomly */
IGRAPH_CHECK(igraph_layout_random(graph, res));
igraph_matrix_scale(res, 1e6);
/* This is the grid for calculating the vertices near to a given vertex */
IGRAPH_CHECK(igraph_2dgrid_init(&grid, res,
-sqrt(area / M_PI), sqrt(area / M_PI), cellsize,
-sqrt(area / M_PI), sqrt(area / M_PI), cellsize));
IGRAPH_FINALLY(igraph_2dgrid_destroy, &grid);
/* Place the root vertex */
igraph_2dgrid_add(&grid, root, 0, 0);
for (actlayer = 1; actlayer < no_of_layers; actlayer++) {
H_n += 1.0 / actlayer;
}
for (actlayer = 1; actlayer < no_of_layers; actlayer++) {
igraph_real_t c = 1;
long int i, j;
igraph_real_t massx, massy;
igraph_real_t px, py;
igraph_real_t sx, sy;
long int it = 0;
igraph_real_t epsilon = 10e-6;
igraph_real_t maxchange = epsilon + 1;
long int pairs;
igraph_real_t sconst = sqrt(area / M_PI) / H_n;
igraph_2dgrid_iterator_t vidit;
/* printf("Layer %li:\n", actlayer); */
/*-----------------------------------------*/
/* Step 1: place the next layer on spheres */
/*-----------------------------------------*/
RNG_BEGIN();
j = (long int) VECTOR(layers)[actlayer];
for (i = (long int) VECTOR(layers)[actlayer - 1];
i < VECTOR(layers)[actlayer]; i++) {
long int vid = (long int) VECTOR(vids)[i];
long int par = (long int) VECTOR(parents)[vid];
IGRAPH_ALLOW_INTERRUPTION();
igraph_2dgrid_getcenter(&grid, &massx, &massy);
igraph_i_norm2d(&massx, &massy);
px = MATRIX(*res, vid, 0) - MATRIX(*res, par, 0);
py = MATRIX(*res, vid, 1) - MATRIX(*res, par, 1);
igraph_i_norm2d(&px, &py);
sx = c * (massx + px) + MATRIX(*res, vid, 0);
sy = c * (massy + py) + MATRIX(*res, vid, 1);
/* The neighbors of 'vid' */
while (j < VECTOR(layers)[actlayer + 1] &&
VECTOR(parents)[(long int)VECTOR(vids)[j]] == vid) {
igraph_real_t rx, ry;
if (actlayer == 1) {
igraph_real_t phi = 2 * M_PI / (VECTOR(layers)[2] - 1) * (j - 1);
rx = cos(phi);
ry = sin(phi);
} else {
rx = RNG_UNIF(-1, 1);
ry = RNG_UNIF(-1, 1);
}
igraph_i_norm2d(&rx, &ry);
rx = rx / actlayer * sconst;
ry = ry / actlayer * sconst;
igraph_2dgrid_add(&grid, (long int) VECTOR(vids)[j], sx + rx, sy + ry);
j++;
}
}
RNG_END();
/*-----------------------------------------*/
/* Step 2: add the edges of the next layer */
/*-----------------------------------------*/
for (j = (long int) VECTOR(layers)[actlayer];
j < VECTOR(layers)[actlayer + 1]; j++) {
long int vid = (long int) VECTOR(vids)[j];
long int k;
IGRAPH_ALLOW_INTERRUPTION();
IGRAPH_CHECK(igraph_incident(graph, &eids, (igraph_integer_t) vid,
IGRAPH_ALL));
for (k = 0; k < igraph_vector_size(&eids); k++) {
long int eid = (long int) VECTOR(eids)[k];
igraph_integer_t from, to;
igraph_edge(graph, (igraph_integer_t) eid, &from, &to);
if ((from != vid && igraph_2dgrid_in(&grid, from)) ||
(to != vid && igraph_2dgrid_in(&grid, to))) {
igraph_vector_push_back(&edges, eid);
}
}
}
/*-----------------------------------------*/
/* Step 3: let the springs spring */
/*-----------------------------------------*/
maxchange = epsilon + 1;
while (it < maxit && maxchange > epsilon) {
long int jj;
igraph_real_t t = maxdelta * pow((maxit - it) / (double)maxit, coolexp);
long int vid, nei;
IGRAPH_PROGRESS("Large graph layout",
100.0 * ((actlayer - 1.0) / (no_of_layers - 1.0) + ((float)it) / (maxit * (no_of_layers - 1.0))),
0);
/* init */
igraph_vector_null(&forcex);
igraph_vector_null(&forcey);
maxchange = 0;
/* attractive "forces" along the edges */
for (jj = 0; jj < igraph_vector_size(&edges); jj++) {
igraph_integer_t from, to;
igraph_real_t xd, yd, dist, force;
IGRAPH_ALLOW_INTERRUPTION();
igraph_edge(graph, (igraph_integer_t) VECTOR(edges)[jj], &from, &to);
xd = MATRIX(*res, (long int)from, 0) - MATRIX(*res, (long int)to, 0);
yd = MATRIX(*res, (long int)from, 1) - MATRIX(*res, (long int)to, 1);
dist = sqrt(xd * xd + yd * yd);
if (dist != 0) {
xd /= dist;
yd /= dist;
}
force = dist * dist / frk;
VECTOR(forcex)[(long int)from] -= xd * force;
VECTOR(forcex)[(long int)to] += xd * force;
VECTOR(forcey)[(long int)from] -= yd * force;
VECTOR(forcey)[(long int)to] += yd * force;
}
/* repulsive "forces" of the vertices nearby */
pairs = 0;
igraph_2dgrid_reset(&grid, &vidit);
while ( (vid = igraph_2dgrid_next(&grid, &vidit) - 1) != -1) {
while ( (nei = igraph_2dgrid_next_nei(&grid, &vidit) - 1) != -1) {
igraph_real_t xd = MATRIX(*res, (long int)vid, 0) -
MATRIX(*res, (long int)nei, 0);
igraph_real_t yd = MATRIX(*res, (long int)vid, 1) -
MATRIX(*res, (long int)nei, 1);
igraph_real_t dist = sqrt(xd * xd + yd * yd);
igraph_real_t force;
if (dist < cellsize) {
pairs++;
if (dist == 0) {
dist = epsilon;
};
xd /= dist; yd /= dist;
force = frk * frk * (1.0 / dist - dist * dist / repulserad);
VECTOR(forcex)[(long int)vid] += xd * force;
VECTOR(forcex)[(long int)nei] -= xd * force;
VECTOR(forcey)[(long int)vid] += yd * force;
VECTOR(forcey)[(long int)nei] -= yd * force;
}
}
}
/* printf("verties: %li iterations: %li\n", */
/* (long int) VECTOR(layers)[actlayer+1], pairs); */
/* apply the changes */
for (jj = 0; jj < VECTOR(layers)[actlayer + 1]; jj++) {
long int vvid = (long int) VECTOR(vids)[jj];
igraph_real_t fx = VECTOR(forcex)[vvid];
igraph_real_t fy = VECTOR(forcey)[vvid];
igraph_real_t ded = sqrt(fx * fx + fy * fy);
if (ded > t) {
ded = t / ded;
fx *= ded; fy *= ded;
}
igraph_2dgrid_move(&grid, vvid, fx, fy);
if (fx > maxchange) {
maxchange = fx;
}
if (fy > maxchange) {
maxchange = fy;
}
}
it++;
/* printf("%li iterations, maxchange: %f\n", it, (double)maxchange); */
}
}
IGRAPH_PROGRESS("Large graph layout", 100.0, 0);
igraph_destroy(&mst);
igraph_vector_destroy(&vids);
igraph_vector_destroy(&layers);
igraph_vector_destroy(&parents);
igraph_vector_destroy(&edges);
igraph_2dgrid_destroy(&grid);
igraph_vector_destroy(&eids);
igraph_vector_destroy(&forcex);
igraph_vector_destroy(&forcey);
IGRAPH_FINALLY_CLEAN(9);
return 0;
}
int igraph_i_layout_reingold_tilford_unreachable(
const igraph_t *graph,
igraph_neimode_t mode,
long int real_root,
long int no_of_nodes,
igraph_vector_t *pnewedges);
int igraph_i_layout_reingold_tilford_unreachable(
const igraph_t *graph,
igraph_neimode_t mode,
long int real_root,
long int no_of_nodes,
igraph_vector_t *pnewedges) {
long int no_of_newedges;
igraph_vector_t visited;
long int i, j, n;
igraph_dqueue_t q = IGRAPH_DQUEUE_NULL;
igraph_adjlist_t allneis;
igraph_vector_int_t *neis;
igraph_vector_resize(pnewedges, 0);
/* traverse from real_root and see what nodes you cannot reach */
no_of_newedges = 0;
IGRAPH_VECTOR_INIT_FINALLY(&visited, no_of_nodes);
IGRAPH_DQUEUE_INIT_FINALLY(&q, 100);
IGRAPH_CHECK(igraph_adjlist_init(graph, &allneis, mode));
IGRAPH_FINALLY(igraph_adjlist_destroy, &allneis);
/* start from real_root and go BFS */
IGRAPH_CHECK(igraph_dqueue_push(&q, real_root));
while (!igraph_dqueue_empty(&q)) {
long int actnode = (long int) igraph_dqueue_pop(&q);
neis = igraph_adjlist_get(&allneis, actnode);
n = igraph_vector_int_size(neis);
VECTOR(visited)[actnode] = 1;
for (j = 0; j < n; j++) {
long int neighbor = (long int) VECTOR(*neis)[j];
if (!(long int)VECTOR(visited)[neighbor]) {
IGRAPH_CHECK(igraph_dqueue_push(&q, neighbor));
}
}
}
for (j = 0; j < no_of_nodes; j++) {
no_of_newedges += 1 - VECTOR(visited)[j];
}
/* if any nodes are unreachable, add edges between them and real_root */
if (no_of_newedges != 0) {
igraph_vector_resize(pnewedges, no_of_newedges * 2);
j = 0;
for (i = 0; i < no_of_nodes; i++) {
if (!VECTOR(visited)[i]) {
if (mode != IGRAPH_IN) {
VECTOR(*pnewedges)[2 * j] = real_root;
VECTOR(*pnewedges)[2 * j + 1] = i;
} else {
VECTOR(*pnewedges)[2 * j] = i;
VECTOR(*pnewedges)[2 * j + 1] = real_root;
}
j++;
}
}
}
igraph_dqueue_destroy(&q);
igraph_adjlist_destroy(&allneis);
igraph_vector_destroy(&visited);
IGRAPH_FINALLY_CLEAN(3);
return IGRAPH_SUCCESS;
}
/* Internal structure for Reingold-Tilford layout */
struct igraph_i_reingold_tilford_vertex {
long int parent; /* Parent node index */
long int level; /* Level of the node */
igraph_real_t offset; /* X offset from parent node */
long int left_contour; /* Next left node of the contour
of the subtree rooted at this node */
long int right_contour; /* Next right node of the contour
of the subtree rooted at this node */
igraph_real_t offset_follow_lc; /* X offset when following the left contour */
igraph_real_t offset_follow_rc; /* X offset when following the right contour */
};
int igraph_i_layout_reingold_tilford_postorder(struct igraph_i_reingold_tilford_vertex *vdata,
long int node, long int vcount);
int igraph_i_layout_reingold_tilford_calc_coords(struct igraph_i_reingold_tilford_vertex *vdata,
igraph_matrix_t *res, long int node,
long int vcount, igraph_real_t xpos);
int igraph_i_layout_reingold_tilford(const igraph_t *graph,
igraph_matrix_t *res,
igraph_neimode_t mode,
long int root);
int igraph_i_layout_reingold_tilford(const igraph_t *graph,
igraph_matrix_t *res,
igraph_neimode_t mode,
long int root) {
long int no_of_nodes = igraph_vcount(graph);
long int i, n, j;
igraph_dqueue_t q = IGRAPH_DQUEUE_NULL;
igraph_adjlist_t allneis;
igraph_vector_int_t *neis;
struct igraph_i_reingold_tilford_vertex *vdata;
IGRAPH_CHECK(igraph_matrix_resize(res, no_of_nodes, 2));
IGRAPH_DQUEUE_INIT_FINALLY(&q, 100);
IGRAPH_CHECK(igraph_adjlist_init(graph, &allneis, mode));
IGRAPH_FINALLY(igraph_adjlist_destroy, &allneis);
vdata = igraph_Calloc(no_of_nodes, struct igraph_i_reingold_tilford_vertex);
if (vdata == 0) {
IGRAPH_ERROR("igraph_layout_reingold_tilford failed", IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(igraph_free, vdata);
for (i = 0; i < no_of_nodes; i++) {
vdata[i].parent = -1;
vdata[i].level = -1;
vdata[i].offset = 0.0;
vdata[i].left_contour = -1;
vdata[i].right_contour = -1;
vdata[i].offset_follow_lc = 0.0;
vdata[i].offset_follow_rc = 0.0;
}
vdata[root].parent = root;
vdata[root].level = 0;
MATRIX(*res, root, 1) = 0;
/* Step 1: assign Y coordinates based on BFS and setup parents vector */
IGRAPH_CHECK(igraph_dqueue_push(&q, root));
IGRAPH_CHECK(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);
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 (vdata[neighbor].parent >= 0) {
continue;
}
MATRIX(*res, neighbor, 1) = actdist + 1;
IGRAPH_CHECK(igraph_dqueue_push(&q, neighbor));
IGRAPH_CHECK(igraph_dqueue_push(&q, actdist + 1));
vdata[neighbor].parent = actnode;
vdata[neighbor].level = actdist + 1;
}
}
/* Step 2: postorder tree traversal, determines the appropriate X
* offsets for every node */
igraph_i_layout_reingold_tilford_postorder(vdata, root, no_of_nodes);
/* Step 3: calculate real coordinates based on X offsets */
igraph_i_layout_reingold_tilford_calc_coords(vdata, res, root, no_of_nodes, vdata[root].offset);
igraph_dqueue_destroy(&q);
igraph_adjlist_destroy(&allneis);
igraph_free(vdata);
IGRAPH_FINALLY_CLEAN(3);
IGRAPH_PROGRESS("Reingold-Tilford tree layout", 100.0, NULL);
return 0;
}
int igraph_i_layout_reingold_tilford_calc_coords(struct igraph_i_reingold_tilford_vertex *vdata,
igraph_matrix_t *res, long int node,
long int vcount, igraph_real_t xpos) {
long int i;
MATRIX(*res, node, 0) = xpos;
for (i = 0; i < vcount; i++) {
if (i == node) {
continue;
}
if (vdata[i].parent == node) {
igraph_i_layout_reingold_tilford_calc_coords(vdata, res, i, vcount,
xpos + vdata[i].offset);
}
}
return 0;
}
int igraph_i_layout_reingold_tilford_postorder(struct igraph_i_reingold_tilford_vertex *vdata,
long int node, long int vcount) {
long int i, j, childcount, leftroot, leftrootidx;
igraph_real_t avg;
/* printf("Starting visiting node %d\n", node); */
/* Check whether this node is a leaf node */
childcount = 0;
for (i = 0; i < vcount; i++) {
if (i == node) {
continue;
}
if (vdata[i].parent == node) {
/* Node i is a child, so visit it recursively */
childcount++;
igraph_i_layout_reingold_tilford_postorder(vdata, i, vcount);
}
}
if (childcount == 0) {
return 0;
}
/* Here we can assume that all of the subtrees have been placed and their
* left and right contours are calculated. Let's place them next to each
* other as close as we can.
* We will take each subtree in an arbitrary order. The root of the
* first one will be placed at offset 0, the next ones will be placed
* as close to each other as possible. leftroot stores the root of the
* rightmost subtree of the already placed subtrees - its right contour
* will be checked against the left contour of the next subtree */
leftroot = leftrootidx = -1;
avg = 0.0;
/*printf("Visited node %d and arranged its subtrees\n", node);*/
for (i = 0, j = 0; i < vcount; i++) {
if (i == node) {
continue;
}
if (vdata[i].parent == node) {
/*printf(" Placing child %d on level %d\n", i, vdata[i].level);*/
if (leftroot >= 0) {
/* Now we will follow the right contour of leftroot and the
* left contour of the subtree rooted at i */
long lnode, rnode;
igraph_real_t loffset, roffset, minsep, rootsep;
lnode = leftroot; rnode = i;
minsep = 1;
rootsep = vdata[leftroot].offset + minsep;
loffset = 0; roffset = minsep;
/*printf(" Contour: [%d, %d], offsets: [%lf, %lf], rootsep: %lf\n",
lnode, rnode, loffset, roffset, rootsep);*/
while ((lnode >= 0) && (rnode >= 0)) {
/* Step to the next level on the right contour of the left subtree */
if (vdata[lnode].right_contour >= 0) {
loffset += vdata[lnode].offset_follow_rc;
lnode = vdata[lnode].right_contour;
} else {
/* Left subtree ended there. The right contour of the left subtree
* will continue to the next step on the right subtree. */
if (vdata[rnode].left_contour >= 0) {
/*printf(" Left subtree ended, continuing left subtree's left and right contour on right subtree (node %ld)\n", vdata[rnode].left_contour);*/
vdata[lnode].left_contour = vdata[rnode].left_contour;
vdata[lnode].right_contour = vdata[rnode].left_contour;
vdata[lnode].offset_follow_lc = vdata[lnode].offset_follow_rc =
(roffset - loffset) + vdata[rnode].offset_follow_lc;
/*printf(" vdata[lnode].offset_follow_* = %.4f\n", vdata[lnode].offset_follow_lc);*/
}
lnode = -1;
}
/* Step to the next level on the left contour of the right subtree */
if (vdata[rnode].left_contour >= 0) {
roffset += vdata[rnode].offset_follow_lc;
rnode = vdata[rnode].left_contour;
} else {
/* Right subtree ended here. The left contour of the right
* subtree will continue to the next step on the left subtree.
* Note that lnode has already been advanced here */
if (lnode >= 0) {
/*printf(" Right subtree ended, continuing right subtree's left and right contour on left subtree (node %ld)\n", lnode);*/
vdata[rnode].left_contour = lnode;
vdata[rnode].right_contour = lnode;
vdata[rnode].offset_follow_lc = vdata[rnode].offset_follow_rc =
(loffset - roffset); /* loffset has also been increased earlier */
/*printf(" vdata[rnode].offset_follow_* = %.4f\n", vdata[rnode].offset_follow_lc);*/
}
rnode = -1;
}
/*printf(" Contour: [%d, %d], offsets: [%lf, %lf], rootsep: %lf\n",
lnode, rnode, loffset, roffset, rootsep);*/
/* Push subtrees away if necessary */
if ((lnode >= 0) && (rnode >= 0) && (roffset - loffset < minsep)) {
/*printf(" Pushing right subtree away by %lf\n", minsep-roffset+loffset);*/
rootsep += minsep - roffset + loffset;
roffset = loffset + minsep;
}
}
/*printf(" Offset of subtree with root node %d will be %lf\n", i, rootsep);*/
vdata[i].offset = rootsep;
vdata[node].right_contour = i;
vdata[node].offset_follow_rc = rootsep;
avg = (avg * j) / (j + 1) + rootsep / (j + 1);
leftrootidx = j;
leftroot = i;
} else {
leftrootidx = j;
leftroot = i;
vdata[node].left_contour = i;
vdata[node].right_contour = i;
vdata[node].offset_follow_lc = 0.0;
vdata[node].offset_follow_rc = 0.0;
avg = vdata[i].offset;
}
j++;
}
}
/*printf("Shifting node to be centered above children. Shift amount: %lf\n", avg);*/
vdata[node].offset_follow_lc -= avg;
vdata[node].offset_follow_rc -= avg;
for (i = 0, j = 0; i < vcount; i++) {
if (i == node) {
continue;
}
if (vdata[i].parent == node) {
vdata[i].offset -= avg;
}
}
return 0;
}
/**
* \function igraph_layout_reingold_tilford
* \brief Reingold-Tilford layout for tree graphs
*
* </para><para>
* Arranges the nodes in a tree where the given node is used as the root.
* The tree is directed downwards and the parents are centered above its
* children. For the exact algorithm, see:
*
* </para><para>
* Reingold, E and Tilford, J: Tidier drawing of trees.
* IEEE Trans. Softw. Eng., SE-7(2):223--228, 1981
*
* </para><para>
* If the given graph is not a tree, a breadth-first search is executed
* first to obtain a possible spanning tree.
*
* \param graph The graph object.
* \param res The result, the coordinates in a matrix. The parameter
* should point to an initialized matrix object and will be resized.
* \param mode Specifies which edges to consider when building the tree.
* If it is \c IGRAPH_OUT then only the outgoing, if it is \c IGRAPH_IN
* then only the incoming edges of a parent are considered. If it is
* \c IGRAPH_ALL then all edges are used (this was the behavior in
* igraph 0.5 and before). This parameter also influences how the root
* vertices are calculated, if they are not given. See the \p roots parameter.
* \param roots The index of the root vertex or root vertices.
* If this is a non-empty vector then the supplied vertex ids are used
* as the roots of the trees (or a single tree if the graph is connected).
* If it is a null pointer of a pointer to an empty vector, then the root
* vertices are automatically calculated based on topological sorting,
* performed with the opposite mode than the \p mode argument.
* After the vertices have been sorted, one is selected from each component.
* \param rootlevel This argument can be useful when drawing forests which are
* not trees (i.e. they are unconnected and have tree components). It specifies
* the level of the root vertices for every tree in the forest. It is only
* considered if not a null pointer and the \p roots argument is also given
* (and it is not a null pointer of an empty vector).
* \return Error code.
*
* Added in version 0.2.
*
* \sa \ref igraph_layout_reingold_tilford_circular().
*
* \example examples/simple/igraph_layout_reingold_tilford.c
*/
int igraph_layout_reingold_tilford(const igraph_t *graph,
igraph_matrix_t *res,
igraph_neimode_t mode,
const igraph_vector_t *roots,
const igraph_vector_t *rootlevel) {
long int no_of_nodes_orig = igraph_vcount(graph);
long int no_of_nodes = no_of_nodes_orig;
long int real_root;
igraph_t extended;
const igraph_t *pextended = graph;
igraph_vector_t myroots;
const igraph_vector_t *proots = roots;
igraph_neimode_t mode2;
long int i;
igraph_vector_t newedges;
/* TODO: possible speedup could be achieved if we use a table for storing
* the children of each node in the tree. (Now the implementation uses a
* single array containing the parent of each node and a node's children
* are determined by looking for other nodes that have this node as parent)
*/
/* at various steps it might be necessary to add edges to the graph */
IGRAPH_VECTOR_INIT_FINALLY(&newedges, 0);
if (!igraph_is_directed(graph)) {
mode = IGRAPH_ALL;
}
if ( (!roots || igraph_vector_size(roots) == 0) &&
rootlevel && igraph_vector_size(rootlevel) != 0 ) {
IGRAPH_WARNING("Reingold-Tilford layout: 'rootlevel' ignored");
}
/* ----------------------------------------------------------------------- */
/* If root vertices are not given, then do a topological sort and take
the last element from every component for directed graphs and mode == out,
or the first element from every component for directed graphs and mode ==
in,or select the vertex with the maximum degree from each component for
undirected graphs */
if (!roots || igraph_vector_size(roots) == 0) {
igraph_vector_t order, membership;
igraph_integer_t no_comps;
long int i, noseen = 0;
IGRAPH_VECTOR_INIT_FINALLY(&myroots, 0);
IGRAPH_VECTOR_INIT_FINALLY(&order, no_of_nodes);
IGRAPH_VECTOR_INIT_FINALLY(&membership, no_of_nodes);
if (mode != IGRAPH_ALL) {
/* look for roots by swimming against the stream */
mode2 = (mode == IGRAPH_IN) ? IGRAPH_OUT : IGRAPH_IN;
IGRAPH_CHECK(igraph_topological_sorting(graph, &order, mode2));
IGRAPH_CHECK(igraph_clusters(graph, &membership, /*csize=*/ 0,
&no_comps, IGRAPH_WEAK));
} else {
IGRAPH_CHECK(igraph_sort_vertex_ids_by_degree(graph, &order,
igraph_vss_all(), IGRAPH_ALL, 0, IGRAPH_ASCENDING, 0));
IGRAPH_CHECK(igraph_clusters(graph, &membership, /*csize=*/ 0,
&no_comps, IGRAPH_WEAK));
}
IGRAPH_CHECK(igraph_vector_resize(&myroots, no_comps));
/* go backwards and fill the roots vector with indices [1, no_of_nodes]
The index 0 is used to signal this root has not been found yet:
all indices are then decreased by one to [0, no_of_nodes - 1] */
igraph_vector_null(&myroots);
proots = &myroots;
for (i = no_of_nodes - 1; noseen < no_comps && i >= 0; i--) {
long int v = (long int) VECTOR(order)[i];
long int mem = (long int) VECTOR(membership)[v];
if (VECTOR(myroots)[mem] == 0) {
noseen += 1;
VECTOR(myroots)[mem] = v + 1;
}
}
for (i = 0; i < no_comps; i++) {
VECTOR(myroots)[i] -= 1;
}
igraph_vector_destroy(&membership);
igraph_vector_destroy(&order);
IGRAPH_FINALLY_CLEAN(2);
} else if (rootlevel && igraph_vector_size(rootlevel) > 0 &&
igraph_vector_size(roots) > 1) {
/* ----------------------------------------------------------------------- */
/* Many roots were given to us, check 'rootlevel' */
long int plus_levels = 0;
long int i;
if (igraph_vector_size(roots) != igraph_vector_size(rootlevel)) {
IGRAPH_ERROR("Reingold-Tilford: 'roots' and 'rootlevel' lengths differ",
IGRAPH_EINVAL);
}
/* count the rootlevels that are not zero */
for (i = 0; i < igraph_vector_size(roots); i++) {
plus_levels += VECTOR(*rootlevel)[i];
}
/* make copy of graph, add vertices/edges */
if (plus_levels != 0) {
long int edgeptr = 0;
pextended = &extended;
IGRAPH_CHECK(igraph_copy(&extended, graph));
IGRAPH_FINALLY(igraph_destroy, &extended);
IGRAPH_CHECK(igraph_add_vertices(&extended,
(igraph_integer_t) plus_levels, 0));
igraph_vector_resize(&newedges, plus_levels * 2);
for (i = 0; i < igraph_vector_size(roots); i++) {
long int rl = (long int) VECTOR(*rootlevel)[i];
long int rn = (long int) VECTOR(*roots)[i];
long int j;
/* zero-level roots don't get anything special */
if (rl == 0) {
continue;
}
/* for each nonzero-level root, add vertices
and edges at all levels [1, 2, .., rl]
piercing through the graph. If mode=="in"
they pierce the other way */
if (mode != IGRAPH_IN) {
VECTOR(newedges)[edgeptr++] = no_of_nodes;
VECTOR(newedges)[edgeptr++] = rn;
for (j = 0; j < rl - 1; j++) {
VECTOR(newedges)[edgeptr++] = no_of_nodes + 1;
VECTOR(newedges)[edgeptr++] = no_of_nodes;
no_of_nodes++;
}
} else {
VECTOR(newedges)[edgeptr++] = rn;
VECTOR(newedges)[edgeptr++] = no_of_nodes;
for (j = 0; j < rl - 1; j++) {
VECTOR(newedges)[edgeptr++] = no_of_nodes;
VECTOR(newedges)[edgeptr++] = no_of_nodes + 1;
no_of_nodes++;
}
}
/* move on to the next root */
VECTOR(*roots)[i] = no_of_nodes++;
}
/* actually add the edges to the graph */
IGRAPH_CHECK(igraph_add_edges(&extended, &newedges, 0));
}
}
/* We have root vertices now. If one or more nonzero-level roots were
chosen by the user, we have copied the graph and added a few vertices
and (directed) edges to connect those floating roots to nonfloating,
zero-level equivalent roots.
Below, the function
igraph_i_layout_reingold_tilford(pextended, res, mode, real_root)
calculates the actual rt coordinates of the graph. However, for
simplicity that function requires a connected graph and a single root.
For directed graphs, it needs not be strongly connected, however all
nodes must be reachable from the root following the stream (i.e. the
root must be a "mother vertex").
So before we call that function we have to make sure the (copied) graph
satisfies that condition. That requires:
1. if there is more than one root, defining a single real_root
2. if a real_root is defined, adding edges to connect all roots to it
3. ensure real_root is mother of the whole graph. If it is not,
add shortcut edges from real_root to any disconnected node for now.
NOTE: 3. could be done better, e.g. by topological sorting of some kind.
But for now it's ok like this.
*/
/* if there is only one root, no need for real_root */
if (igraph_vector_size(proots) == 1) {
real_root = (long int) VECTOR(*proots)[0];
if (real_root < 0 || real_root >= no_of_nodes) {
IGRAPH_ERROR("invalid vertex id", IGRAPH_EINVVID);
}
/* else, we need to make real_root */
} else {
long int no_of_newedges;
/* Make copy of the graph unless it exists already */
if (pextended == graph) {
pextended = &extended;
IGRAPH_CHECK(igraph_copy(&extended, graph));
IGRAPH_FINALLY(igraph_destroy, &extended);
}
/* add real_root to the vertices */
real_root = no_of_nodes;
IGRAPH_CHECK(igraph_add_vertices(&extended, 1, 0));
no_of_nodes++;
/* add edges from the roots to real_root */
no_of_newedges = igraph_vector_size(proots);
igraph_vector_resize(&newedges, no_of_newedges * 2);
for (i = 0; i < no_of_newedges; i++) {
VECTOR(newedges)[2 * i] = no_of_nodes - 1;
VECTOR(newedges)[2 * i + 1] = VECTOR(*proots)[i];
}
IGRAPH_CHECK(igraph_add_edges(&extended, &newedges, 0));
}
/* prepare edges to unreachable parts of the graph */
IGRAPH_CHECK(igraph_i_layout_reingold_tilford_unreachable(pextended, mode, real_root, no_of_nodes, &newedges));
if (igraph_vector_size(&newedges) != 0) {
/* Make copy of the graph unless it exists already */
if (pextended == graph) {
pextended = &extended;
IGRAPH_CHECK(igraph_copy(&extended, graph));
IGRAPH_FINALLY(igraph_destroy, &extended);
}
IGRAPH_CHECK(igraph_add_edges(&extended, &newedges, 0));
}
igraph_vector_destroy(&newedges);
IGRAPH_FINALLY_CLEAN(1);
/* ----------------------------------------------------------------------- */
/* Layout */
IGRAPH_CHECK(igraph_i_layout_reingold_tilford(pextended, res, mode, real_root));
/* Remove the new vertices from the layout */
if (no_of_nodes != no_of_nodes_orig) {
if (no_of_nodes - 1 == no_of_nodes_orig) {
IGRAPH_CHECK(igraph_matrix_remove_row(res, no_of_nodes_orig));
} else {
igraph_matrix_t tmp;
long int i;
IGRAPH_MATRIX_INIT_FINALLY(&tmp, no_of_nodes_orig, 2);
for (i = 0; i < no_of_nodes_orig; i++) {
MATRIX(tmp, i, 0) = MATRIX(*res, i, 0);
MATRIX(tmp, i, 1) = MATRIX(*res, i, 1);
}
IGRAPH_CHECK(igraph_matrix_update(res, &tmp));
igraph_matrix_destroy(&tmp);
IGRAPH_FINALLY_CLEAN(1);
}
}
if (pextended != graph) {
igraph_destroy(&extended);
IGRAPH_FINALLY_CLEAN(1);
}
/* Remove the roots vector if it was created by us */
if (proots != roots) {
igraph_vector_destroy(&myroots);
IGRAPH_FINALLY_CLEAN(1);
}
return 0;
}
/**
* \function igraph_layout_reingold_tilford_circular
* \brief Circular Reingold-Tilford layout for trees
*
* </para><para>
* This layout is almost the same as \ref igraph_layout_reingold_tilford(), but
* the tree is drawn in a circular way, with the root vertex in the center.
*
* \param graph The graph object.
* \param res The result, the coordinates in a matrix. The parameter
* should point to an initialized matrix object and will be resized.
* \param mode Specifies which edges to consider when building the tree.
* If it is \c IGRAPH_OUT then only the outgoing, if it is \c IGRAPH_IN
* then only the incoming edges of a parent are considered. If it is
* \c IGRAPH_ALL then all edges are used (this was the behavior in
* igraph 0.5 and before). This parameter also influences how the root
* vertices are calculated, if they are not given. See the \p roots parameter.
* \param roots The index of the root vertex or root vertices.
* If this is a non-empty vector then the supplied vertex ids are used
* as the roots of the trees (or a single tree if the graph is connected).
* If it is a null pointer of a pointer to an empty vector, then the root
* vertices are automatically calculated based on topological sorting,
* performed with the opposite mode than the \p mode argument.
* After the vertices have been sorted, one is selected from each component.
* \param rootlevel This argument can be useful when drawing forests which are
* not trees (i.e. they are unconnected and have tree components). It specifies
* the level of the root vertices for every tree in the forest. It is only
* considered if not a null pointer and the \p roots argument is also given
* (and it is not a null pointer of an empty vector). Note that if you supply
* a null pointer here and the graph has multiple components, all of the root
* vertices will be mapped to the origin of the coordinate system, which does
* not really make sense.
* \return Error code.
*
* \sa \ref igraph_layout_reingold_tilford().
*/
int igraph_layout_reingold_tilford_circular(const igraph_t *graph,
igraph_matrix_t *res,
igraph_neimode_t mode,
const igraph_vector_t *roots,
const igraph_vector_t *rootlevel) {
long int no_of_nodes = igraph_vcount(graph);
long int i;
igraph_real_t ratio = 2 * M_PI * (no_of_nodes - 1.0) / no_of_nodes;
igraph_real_t minx, maxx;
IGRAPH_CHECK(igraph_layout_reingold_tilford(graph, res, mode, roots, rootlevel));
if (no_of_nodes == 0) {
return 0;
}
minx = maxx = MATRIX(*res, 0, 0);
for (i = 1; i < no_of_nodes; i++) {
if (MATRIX(*res, i, 0) > maxx) {
maxx = MATRIX(*res, i, 0);
}
if (MATRIX(*res, i, 0) < minx) {
minx = MATRIX(*res, i, 0);
}
}
if (maxx > minx) {
ratio /= (maxx - minx);
}
for (i = 0; i < no_of_nodes; i++) {
igraph_real_t phi = (MATRIX(*res, i, 0) - minx) * ratio;
igraph_real_t r = MATRIX(*res, i, 1);
MATRIX(*res, i, 0) = r * cos(phi);
MATRIX(*res, i, 1) = r * sin(phi);
}
return 0;
}
#define COULOMBS_CONSTANT 8987500000.0
igraph_real_t igraph_i_distance_between(const igraph_matrix_t *c, long int a,
long int b);
int igraph_i_determine_electric_axal_forces(const igraph_matrix_t *pos,
igraph_real_t *x,
igraph_real_t *y,
igraph_real_t directed_force,
igraph_real_t distance,
long int other_node,
long int this_node);
int igraph_i_apply_electrical_force(const igraph_matrix_t *pos,
igraph_vector_t *pending_forces_x,
igraph_vector_t *pending_forces_y,
long int other_node, long int this_node,
igraph_real_t node_charge,
igraph_real_t distance);
int igraph_i_determine_spring_axal_forces(const igraph_matrix_t *pos,
igraph_real_t *x, igraph_real_t *y,
igraph_real_t directed_force,
igraph_real_t distance,
int spring_length,
long int other_node,
long int this_node);
int igraph_i_apply_spring_force(const igraph_matrix_t *pos,
igraph_vector_t *pending_forces_x,
igraph_vector_t *pending_forces_y,
long int other_node,
long int this_node, int spring_length,
igraph_real_t spring_constant);
int igraph_i_move_nodes(igraph_matrix_t *pos,
const igraph_vector_t *pending_forces_x,
const igraph_vector_t *pending_forces_y,
igraph_real_t node_mass,
igraph_real_t max_sa_movement);
igraph_real_t igraph_i_distance_between(const igraph_matrix_t *c, long int a,
long int b) {
igraph_real_t diffx = MATRIX(*c, a, 0) - MATRIX(*c, b, 0);
igraph_real_t diffy = MATRIX(*c, a, 1) - MATRIX(*c, b, 1);
return sqrt( diffx * diffx + diffy * diffy );
}
int igraph_i_determine_electric_axal_forces(const igraph_matrix_t *pos,
igraph_real_t *x,
igraph_real_t *y,
igraph_real_t directed_force,
igraph_real_t distance,
long int other_node,
long int this_node) {
// We know what the directed force is. We now need to translate it
// into the appropriate x and y components.
// First, assume:
// other_node
// /|
// directed_force / |
// / | y
// /______|
// this_node x
//
// other_node.x > this_node.x
// other_node.y > this_node.y
// the force will be on this_node away from other_node
// the proportion (distance/y_distance) is equal to the proportion
// (directed_force/y_force), as the two triangles are similar.
// therefore, the magnitude of y_force = (directed_force*y_distance)/distance
// the sign of y_force is negative, away from other_node
igraph_real_t x_distance, y_distance;
y_distance = MATRIX(*pos, other_node, 1) - MATRIX(*pos, this_node, 1);
if (y_distance < 0) {
y_distance = -y_distance;
}
*y = -1 * ((directed_force * y_distance) / distance);
// the x component works in exactly the same way.
x_distance = MATRIX(*pos, other_node, 0) - MATRIX(*pos, this_node, 0);
if (x_distance < 0) {
x_distance = -x_distance;
}
*x = -1 * ((directed_force * x_distance) / distance);
// Now we need to reverse the polarity of our answers based on the falsness
// of our assumptions.
if (MATRIX(*pos, other_node, 0) < MATRIX(*pos, this_node, 0)) {
*x = *x * -1;
}
if (MATRIX(*pos, other_node, 1) < MATRIX(*pos, this_node, 1)) {
*y = *y * -1;
}
return 0;
}
int igraph_i_apply_electrical_force(const igraph_matrix_t *pos,
igraph_vector_t *pending_forces_x,
igraph_vector_t *pending_forces_y,
long int other_node, long int this_node,
igraph_real_t node_charge,
igraph_real_t distance) {
igraph_real_t directed_force = COULOMBS_CONSTANT *
((node_charge * node_charge) / (distance * distance));
igraph_real_t x_force, y_force;
igraph_i_determine_electric_axal_forces(pos, &x_force, &y_force,
directed_force, distance,
other_node, this_node);
VECTOR(*pending_forces_x)[this_node] += x_force;
VECTOR(*pending_forces_y)[this_node] += y_force;
VECTOR(*pending_forces_x)[other_node] -= x_force;
VECTOR(*pending_forces_y)[other_node] -= y_force;
return 0;
}
int igraph_i_determine_spring_axal_forces(const igraph_matrix_t *pos,
igraph_real_t *x, igraph_real_t *y,
igraph_real_t directed_force,
igraph_real_t distance,
int spring_length,
long int other_node, long int this_node) {
// if the spring is just the right size, the forces will be 0, so we can
// skip the computation.
//
// if the spring is too long, our forces will be identical to those computed
// by determine_electrical_axal_forces() (this_node will be pulled toward
// other_node).
//
// if the spring is too short, our forces will be the opposite of those
// computed by determine_electrical_axal_forces() (this_node will be pushed
// away from other_node)
//
// finally, since both nodes are movable, only one-half of the total force
// should be applied to each node, so half the forces for our answer.
if (distance == spring_length) {
*x = 0.0;
*y = 0.0;
} else {
igraph_i_determine_electric_axal_forces(pos, x, y, directed_force, distance,
other_node, this_node);
if (distance < spring_length) {
*x = -1 * *x;
*y = -1 * *y;
}
*x = 0.5 * *x;
*y = 0.5 * *y;
}
return 0;
}
int igraph_i_apply_spring_force(const igraph_matrix_t *pos,
igraph_vector_t *pending_forces_x,
igraph_vector_t *pending_forces_y,
long int other_node,
long int this_node, int spring_length,
igraph_real_t spring_constant) {
// determined using Hooke's Law:
// force = -kx
// where:
// k = spring constant
// x = displacement from ideal length in meters
igraph_real_t distance, displacement, directed_force, x_force, y_force;
distance = igraph_i_distance_between(pos, other_node, this_node);
// let's protect ourselves from division by zero by ignoring two nodes that
// happen to be in the same place. Since we separate all nodes before we
// work on any of them, this will only happen in extremely rare circumstances,
// and when it does, electrical force will probably push one or both of them
// one way or another anyway.
if (distance == 0.0) {
return 0;
}
displacement = distance - spring_length;
if (displacement < 0) {
displacement = -displacement;
}
directed_force = -1 * spring_constant * displacement;
// remember, this is force directed away from the spring;
// a negative number is back towards the spring (or, in our case, back towards
// the other node)
// get the force that should be applied to >this< node
igraph_i_determine_spring_axal_forces(pos, &x_force, &y_force,
directed_force, distance, spring_length,
other_node, this_node);
VECTOR(*pending_forces_x)[this_node] += x_force;
VECTOR(*pending_forces_y)[this_node] += y_force;
VECTOR(*pending_forces_x)[other_node] -= x_force;
VECTOR(*pending_forces_y)[other_node] -= y_force;
return 0;
}
int igraph_i_move_nodes(igraph_matrix_t *pos,
const igraph_vector_t *pending_forces_x,
const igraph_vector_t *pending_forces_y,
igraph_real_t node_mass,
igraph_real_t max_sa_movement) {
// Since each iteration is isolated, time is constant at 1.
// Therefore:
// Force effects acceleration.
// acceleration (d(velocity)/time) = velocity
// velocity (d(displacement)/time) = displacement
// displacement = acceleration
// determined using Newton's second law:
// sum(F) = ma
// therefore:
// acceleration = force / mass
// velocity = force / mass
// displacement = force / mass
long int this_node, no_of_nodes = igraph_vector_size(pending_forces_x);
for (this_node = 0; this_node < no_of_nodes; this_node++) {
igraph_real_t x_movement, y_movement;
x_movement = VECTOR(*pending_forces_x)[this_node] / node_mass;
if (x_movement > max_sa_movement) {
x_movement = max_sa_movement;
} else if (x_movement < -max_sa_movement) {
x_movement = -max_sa_movement;
}
y_movement = VECTOR(*pending_forces_y)[this_node] / node_mass;
if (y_movement > max_sa_movement) {
y_movement = max_sa_movement;
} else if (y_movement < -max_sa_movement) {
y_movement = -max_sa_movement;
}
MATRIX(*pos, this_node, 0) += x_movement;
MATRIX(*pos, this_node, 1) += y_movement;
}
return 0;
}
/**
* \function igraph_layout_graphopt
* \brief Optimizes vertex layout via the graphopt algorithm.
*
* </para><para>
* This is a port of the graphopt layout algorithm by Michael Schmuhl.
* graphopt version 0.4.1 was rewritten in C and the support for
* layers was removed (might be added later) and a code was a bit
* reorganized to avoid some unnecessary steps is the node charge (see below)
* is zero.
*
* </para><para>
* graphopt uses physical analogies for defining attracting and repelling
* forces among the vertices and then the physical system is simulated
* until it reaches an equilibrium. (There is no simulated annealing or
* anything like that, so a stable fixed point is not guaranteed.)
*
* </para><para>
* See also http://www.schmuhl.org/graphopt/ for the original graphopt.
* \param graph The input graph.
* \param res Pointer to an initialized matrix, the result will be stored here
* and its initial contents is used the starting point of the simulation
* if the \p use_seed argument is true. Note that in this case the
* matrix should have the proper size, otherwise a warning is issued and
* the supplied values are ignored. If no starting positions are given
* (or they are invalid) then a random staring position is used.
* The matrix will be resized if needed.
* \param niter Integer constant, the number of iterations to perform.
* Should be a couple of hundred in general. If you have a large graph
* then you might want to only do a few iterations and then check the
* result. If it is not good enough you can feed it in again in
* the \p res argument. The original graphopt default if 500.
* \param node_charge The charge of the vertices, used to calculate electric
* repulsion. The original graphopt default is 0.001.
* \param node_mass The mass of the vertices, used for the spring forces.
* The original graphopt defaults to 30.
* \param spring_length The length of the springs, an integer number.
* The original graphopt defaults to zero.
* \param spring_constant The spring constant, the original graphopt defaults
* to one.
* \param max_sa_movement Real constant, it gives the maximum amount of movement
* allowed in a single step along a single axis. The original graphopt
* default is 5.
* \param use_seed Logical scalar, whether to use the positions in \p res as
* a starting configuration. See also \p res above.
* \return Error code.
*
* Time complexity: O(n (|V|^2+|E|) ), n is the number of iterations,
* |V| is the number of vertices, |E| the number
* of edges. If \p node_charge is zero then it is only O(n|E|).
*/
int igraph_layout_graphopt(const igraph_t *graph, igraph_matrix_t *res,
igraph_integer_t niter,
igraph_real_t node_charge, igraph_real_t node_mass,
igraph_real_t spring_length,
igraph_real_t spring_constant,
igraph_real_t max_sa_movement,
igraph_bool_t use_seed) {
long int no_of_nodes = igraph_vcount(graph);
long int no_of_edges = igraph_ecount(graph);
int my_spring_length = (int) spring_length;
igraph_vector_t pending_forces_x, pending_forces_y;
/* Set a flag to calculate (or not) the electrical forces that the nodes */
/* apply on each other based on if both node types' charges are zero. */
igraph_bool_t apply_electric_charges = (node_charge != 0);
long int this_node, other_node, edge;
igraph_real_t distance;
long int i;
IGRAPH_VECTOR_INIT_FINALLY(&pending_forces_x, no_of_nodes);
IGRAPH_VECTOR_INIT_FINALLY(&pending_forces_y, no_of_nodes);
if (use_seed) {
if (igraph_matrix_nrow(res) != no_of_nodes ||
igraph_matrix_ncol(res) != 2) {
IGRAPH_WARNING("Invalid size for initial matrix, starting from random layout");
IGRAPH_CHECK(igraph_layout_random(graph, res));
}
} else {
IGRAPH_CHECK(igraph_layout_random(graph, res));
}
IGRAPH_PROGRESS("Graphopt layout", 0, NULL);
for (i = niter; i > 0; i--) {
/* Report progress in approx. every 100th step */
if (i % 10 == 0) {
IGRAPH_PROGRESS("Graphopt layout", 100.0 - 100.0 * i / niter, NULL);
}
/* Clear pending forces on all nodes */
igraph_vector_null(&pending_forces_x);
igraph_vector_null(&pending_forces_y);
// Apply electrical force applied by all other nodes
if (apply_electric_charges) {
// Iterate through all nodes
for (this_node = 0; this_node < no_of_nodes; this_node++) {
IGRAPH_ALLOW_INTERRUPTION();
for (other_node = this_node + 1;
other_node < no_of_nodes;
other_node++) {
distance = igraph_i_distance_between(res, this_node, other_node);
// let's protect ourselves from division by zero by ignoring
// two nodes that happen to be in the same place. Since we
// separate all nodes before we work on any of them, this
// will only happen in extremely rare circumstances, and when
// it does, springs will probably pull them apart anyway.
// also, if we are more than 50 away, the electric force
// will be negligible.
// ***** may not always be desirable ****
if ((distance != 0.0) && (distance < 500.0)) {
// if (distance != 0.0) {
// Apply electrical force from node(counter2) on
// node(counter)
igraph_i_apply_electrical_force(res, &pending_forces_x,
&pending_forces_y,
other_node, this_node,
node_charge,
distance);
}
}
}
}
// Apply force from springs
for (edge = 0; edge < no_of_edges; edge++) {
long int tthis_node = IGRAPH_FROM(graph, edge);
long int oother_node = IGRAPH_TO(graph, edge);
// Apply spring force on both nodes
igraph_i_apply_spring_force(res, &pending_forces_x, &pending_forces_y,
oother_node, tthis_node, my_spring_length,
spring_constant);
}
// Effect the movement of the nodes based on all pending forces
igraph_i_move_nodes(res, &pending_forces_x, &pending_forces_y, node_mass,
max_sa_movement);
}
IGRAPH_PROGRESS("Graphopt layout", 100, NULL);
igraph_vector_destroy(&pending_forces_y);
igraph_vector_destroy(&pending_forces_x);
IGRAPH_FINALLY_CLEAN(2);
return 0;
}
int igraph_i_layout_merge_dla(igraph_i_layout_mergegrid_t *grid,
long int actg, igraph_real_t *x, igraph_real_t *y, igraph_real_t r,
igraph_real_t cx, igraph_real_t cy, igraph_real_t startr,
igraph_real_t killr);
int igraph_i_layout_sphere_2d(igraph_matrix_t *coords, igraph_real_t *x,
igraph_real_t *y, igraph_real_t *r);
int igraph_i_layout_sphere_3d(igraph_matrix_t *coords, igraph_real_t *x,
igraph_real_t *y, igraph_real_t *z,
igraph_real_t *r);
/**
* \function igraph_layout_merge_dla
* \brief Merge multiple layouts by using a DLA algorithm
*
* </para><para>
* First each layout is covered by a circle. Then the layout of the
* largest graph is placed at the origin. Then the other layouts are
* placed by the DLA algorithm, larger ones first and smaller ones
* last.
* \param thegraphs Pointer vector containing the graph object of
* which the layouts will be merged.
* \param coords Pointer vector containing matrix objects with the 2d
* layouts of the graphs in \p thegraphs.
* \param res Pointer to an initialized matrix object, the result will
* be stored here. It will be resized if needed.
* \return Error code.
*
* Added in version 0.2. This function is experimental.
*
* </para><para>
* Time complexity: TODO.
*/
int igraph_layout_merge_dla(igraph_vector_ptr_t *thegraphs,
igraph_vector_ptr_t *coords,
igraph_matrix_t *res) {
long int graphs = igraph_vector_ptr_size(coords);
igraph_vector_t sizes;
igraph_vector_t x, y, r;
igraph_vector_t nx, ny, nr;
long int allnodes = 0;
long int i, j;
long int actg;
igraph_i_layout_mergegrid_t grid;
long int jpos = 0;
igraph_real_t minx, maxx, miny, maxy;
igraph_real_t area = 0;
igraph_real_t maxr = 0;
long int respos;
/* Graphs are currently not used, only the coordinates */
IGRAPH_UNUSED(thegraphs);
IGRAPH_VECTOR_INIT_FINALLY(&sizes, graphs);
IGRAPH_VECTOR_INIT_FINALLY(&x, graphs);
IGRAPH_VECTOR_INIT_FINALLY(&y, graphs);
IGRAPH_VECTOR_INIT_FINALLY(&r, graphs);
IGRAPH_VECTOR_INIT_FINALLY(&nx, graphs);
IGRAPH_VECTOR_INIT_FINALLY(&ny, graphs);
IGRAPH_VECTOR_INIT_FINALLY(&nr, graphs);
RNG_BEGIN();
for (i = 0; i < igraph_vector_ptr_size(coords); i++) {
igraph_matrix_t *mat = VECTOR(*coords)[i];
long int size = igraph_matrix_nrow(mat);
if (igraph_matrix_ncol(mat) != 2) {
IGRAPH_ERROR("igraph_layout_merge_dla works for 2D layouts only",
IGRAPH_EINVAL);
}
IGRAPH_ALLOW_INTERRUPTION();
allnodes += size;
VECTOR(sizes)[i] = size;
VECTOR(r)[i] = pow(size, .75);
area += VECTOR(r)[i] * VECTOR(r)[i];
if (VECTOR(r)[i] > maxr) {
maxr = VECTOR(r)[i];
}
igraph_i_layout_sphere_2d(mat,
igraph_vector_e_ptr(&nx, i),
igraph_vector_e_ptr(&ny, i),
igraph_vector_e_ptr(&nr, i));
}
igraph_vector_order2(&sizes); /* largest first */
/* 0. create grid */
minx = miny = -sqrt(5 * area);
maxx = maxy = sqrt(5 * area);
igraph_i_layout_mergegrid_init(&grid, minx, maxx, 200,
miny, maxy, 200);
IGRAPH_FINALLY(igraph_i_layout_mergegrid_destroy, &grid);
/* fprintf(stderr, "Ok, starting DLA\n"); */
/* 1. place the largest */
actg = (long int) VECTOR(sizes)[jpos++];
igraph_i_layout_merge_place_sphere(&grid, 0, 0, VECTOR(r)[actg], actg);
IGRAPH_PROGRESS("Merging layouts via DLA", 0.0, NULL);
while (jpos < graphs) {
IGRAPH_ALLOW_INTERRUPTION();
/* fprintf(stderr, "comp: %li", jpos); */
IGRAPH_PROGRESS("Merging layouts via DLA", (100.0 * jpos) / graphs, NULL);
actg = (long int) VECTOR(sizes)[jpos++];
/* 2. random walk, TODO: tune parameters */
igraph_i_layout_merge_dla(&grid, actg,
igraph_vector_e_ptr(&x, actg),
igraph_vector_e_ptr(&y, actg),
VECTOR(r)[actg], 0, 0,
maxx, maxx + 5);
/* 3. place sphere */
igraph_i_layout_merge_place_sphere(&grid, VECTOR(x)[actg], VECTOR(y)[actg],
VECTOR(r)[actg], actg);
}
IGRAPH_PROGRESS("Merging layouts via DLA", 100.0, NULL);
/* Create the result */
IGRAPH_CHECK(igraph_matrix_resize(res, allnodes, 2));
respos = 0;
for (i = 0; i < graphs; i++) {
long int size = igraph_matrix_nrow(VECTOR(*coords)[i]);
igraph_real_t xx = VECTOR(x)[i];
igraph_real_t yy = VECTOR(y)[i];
igraph_real_t rr = VECTOR(r)[i] / VECTOR(nr)[i];
igraph_matrix_t *mat = VECTOR(*coords)[i];
IGRAPH_ALLOW_INTERRUPTION();
if (VECTOR(nr)[i] == 0) {
rr = 1;
}
for (j = 0; j < size; j++) {
MATRIX(*res, respos, 0) = rr * (MATRIX(*mat, j, 0) - VECTOR(nx)[i]);
MATRIX(*res, respos, 1) = rr * (MATRIX(*mat, j, 1) - VECTOR(ny)[i]);
MATRIX(*res, respos, 0) += xx;
MATRIX(*res, respos, 1) += yy;
++respos;
}
}
RNG_END();
igraph_i_layout_mergegrid_destroy(&grid);
igraph_vector_destroy(&sizes);
igraph_vector_destroy(&x);
igraph_vector_destroy(&y);
igraph_vector_destroy(&r);
igraph_vector_destroy(&nx);
igraph_vector_destroy(&ny);
igraph_vector_destroy(&nr);
IGRAPH_FINALLY_CLEAN(8);
return 0;
}
int igraph_i_layout_sphere_2d(igraph_matrix_t *coords, igraph_real_t *x, igraph_real_t *y,
igraph_real_t *r) {
long int nodes = igraph_matrix_nrow(coords);
long int i;
igraph_real_t xmin, xmax, ymin, ymax;
xmin = xmax = MATRIX(*coords, 0, 0);
ymin = ymax = MATRIX(*coords, 0, 1);
for (i = 1; i < nodes; i++) {
if (MATRIX(*coords, i, 0) < xmin) {
xmin = MATRIX(*coords, i, 0);
} else if (MATRIX(*coords, i, 0) > xmax) {
xmax = MATRIX(*coords, i, 0);
}
if (MATRIX(*coords, i, 1) < ymin) {
ymin = MATRIX(*coords, i, 1);
} else if (MATRIX(*coords, i, 1) > ymax) {
ymax = MATRIX(*coords, i, 1);
}
}
*x = (xmin + xmax) / 2;
*y = (ymin + ymax) / 2;
*r = sqrt( (xmax - xmin) * (xmax - xmin) + (ymax - ymin) * (ymax - ymin) ) / 2;
return 0;
}
int igraph_i_layout_sphere_3d(igraph_matrix_t *coords, igraph_real_t *x, igraph_real_t *y,
igraph_real_t *z, igraph_real_t *r) {
long int nodes = igraph_matrix_nrow(coords);
long int i;
igraph_real_t xmin, xmax, ymin, ymax, zmin, zmax;
xmin = xmax = MATRIX(*coords, 0, 0);
ymin = ymax = MATRIX(*coords, 0, 1);
zmin = zmax = MATRIX(*coords, 0, 2);
for (i = 1; i < nodes; i++) {
if (MATRIX(*coords, i, 0) < xmin) {
xmin = MATRIX(*coords, i, 0);
} else if (MATRIX(*coords, i, 0) > xmax) {
xmax = MATRIX(*coords, i, 0);
}
if (MATRIX(*coords, i, 1) < ymin) {
ymin = MATRIX(*coords, i, 1);
} else if (MATRIX(*coords, i, 1) > ymax) {
ymax = MATRIX(*coords, i, 1);
}
if (MATRIX(*coords, i, 2) < zmin) {
zmin = MATRIX(*coords, i, 2);
} else if (MATRIX(*coords, i, 2) > zmax) {
zmax = MATRIX(*coords, i, 2);
}
}
*x = (xmin + xmax) / 2;
*y = (ymin + ymax) / 2;
*z = (zmin + zmax) / 2;
*r = sqrt( (xmax - xmin) * (xmax - xmin) + (ymax - ymin) * (ymax - ymin) +
(zmax - zmin) * (zmax - zmin) ) / 2;
return 0;
}
#define DIST(x,y) (sqrt(pow((x)-cx,2)+pow((y)-cy,2)))
int igraph_i_layout_merge_dla(igraph_i_layout_mergegrid_t *grid,
long int actg, igraph_real_t *x, igraph_real_t *y, igraph_real_t r,
igraph_real_t cx, igraph_real_t cy, igraph_real_t startr,
igraph_real_t killr) {
long int sp = -1;
igraph_real_t angle, len;
long int steps = 0;
/* The graph is not used, only its coordinates */
IGRAPH_UNUSED(actg);
while (sp < 0) {
/* start particle */
do {
steps++;
angle = RNG_UNIF(0, 2 * M_PI);
len = RNG_UNIF(.5 * startr, startr);
*x = cx + len * cos(angle);
*y = cy + len * sin(angle);
sp = igraph_i_layout_mergegrid_get_sphere(grid, *x, *y, r);
} while (sp >= 0);
while (sp < 0 && DIST(*x, *y) < killr) {
igraph_real_t nx, ny;
steps++;
angle = RNG_UNIF(0, 2 * M_PI);
len = RNG_UNIF(0, startr / 100);
nx = *x + len * cos(angle);
ny = *y + len * sin(angle);
sp = igraph_i_layout_mergegrid_get_sphere(grid, nx, ny, r);
if (sp < 0) {
*x = nx; *y = ny;
}
}
}
/* fprintf(stderr, "%li ", steps); */
return 0;
}
int igraph_i_layout_mds_step(igraph_real_t *to, const igraph_real_t *from,
int n, void *extra);
int igraph_i_layout_mds_single(const igraph_t* graph, igraph_matrix_t *res,
igraph_matrix_t *dist, long int dim);
int igraph_i_layout_mds_step(igraph_real_t *to, const igraph_real_t *from,
int n, void *extra) {
igraph_matrix_t* matrix = (igraph_matrix_t*)extra;
IGRAPH_UNUSED(n);
igraph_blas_dgemv_array(0, 1, matrix, from, 0, to);
return 0;
}
/* MDS layout for a connected graph, with no error checking on the
* input parameters. The distance matrix will be modified in-place. */
int igraph_i_layout_mds_single(const igraph_t* graph, igraph_matrix_t *res,
igraph_matrix_t *dist, long int dim) {
long int no_of_nodes = igraph_vcount(graph);
long int nev = dim;
igraph_matrix_t vectors;
igraph_vector_t values, row_means;
igraph_real_t grand_mean;
long int i, j, k;
igraph_eigen_which_t which;
/* Handle the trivial cases */
if (no_of_nodes == 1) {
IGRAPH_CHECK(igraph_matrix_resize(res, 1, dim));
igraph_matrix_fill(res, 0);
return IGRAPH_SUCCESS;
}
if (no_of_nodes == 2) {
IGRAPH_CHECK(igraph_matrix_resize(res, 2, dim));
igraph_matrix_fill(res, 0);
for (j = 0; j < dim; j++) {
MATRIX(*res, 1, j) = 1;
}
return IGRAPH_SUCCESS;
}
/* Initialize some stuff */
IGRAPH_VECTOR_INIT_FINALLY(&values, no_of_nodes);
IGRAPH_CHECK(igraph_matrix_init(&vectors, no_of_nodes, dim));
IGRAPH_FINALLY(igraph_matrix_destroy, &vectors);
/* Take the square of the distance matrix */
for (i = 0; i < no_of_nodes; i++) {
for (j = 0; j < no_of_nodes; j++) {
MATRIX(*dist, i, j) *= MATRIX(*dist, i, j);
}
}
/* Double centering of the distance matrix */
IGRAPH_VECTOR_INIT_FINALLY(&row_means, no_of_nodes);
igraph_vector_fill(&values, 1.0 / no_of_nodes);
igraph_blas_dgemv(0, 1, dist, &values, 0, &row_means);
grand_mean = igraph_vector_sum(&row_means) / no_of_nodes;
igraph_matrix_add_constant(dist, grand_mean);
for (i = 0; i < no_of_nodes; i++) {
for (j = 0; j < no_of_nodes; j++) {
MATRIX(*dist, i, j) -= VECTOR(row_means)[i] + VECTOR(row_means)[j];
MATRIX(*dist, i, j) *= -0.5;
}
}
igraph_vector_destroy(&row_means);
IGRAPH_FINALLY_CLEAN(1);
/* Calculate the top `dim` eigenvectors. */
which.pos = IGRAPH_EIGEN_LA;
which.howmany = (int) nev;
IGRAPH_CHECK(igraph_eigen_matrix_symmetric(/*A=*/ 0, /*sA=*/ 0,
/*fun=*/ igraph_i_layout_mds_step,
/*n=*/ (int) no_of_nodes, /*extra=*/ dist,
/*algorithm=*/ IGRAPH_EIGEN_LAPACK,
&which, /*options=*/ 0, /*storage=*/ 0,
&values, &vectors));
/* Calculate and normalize the final coordinates */
for (j = 0; j < nev; j++) {
VECTOR(values)[j] = sqrt(fabs(VECTOR(values)[j]));
}
IGRAPH_CHECK(igraph_matrix_resize(res, no_of_nodes, dim));
for (i = 0; i < no_of_nodes; i++) {
for (j = 0, k = nev - 1; j < nev; j++, k--) {
MATRIX(*res, i, k) = VECTOR(values)[j] * MATRIX(vectors, i, j);
}
}
igraph_matrix_destroy(&vectors);
igraph_vector_destroy(&values);
IGRAPH_FINALLY_CLEAN(2);
return IGRAPH_SUCCESS;
}
/**
* \function igraph_layout_mds
* \brief Place the vertices on a plane using multidimensional scaling.
*
* </para><para>
* This layout requires a distance matrix, where the intersection of
* row i and column j specifies the desired distance between vertex i
* and vertex j. The algorithm will try to place the vertices in a
* space having a given number of dimensions in a way that approximates
* the distance relations prescribed in the distance matrix. igraph
* uses the classical multidimensional scaling by Torgerson; for more
* details, see Cox & Cox: Multidimensional Scaling (1994), Chapman
* and Hall, London.
*
* </para><para>
* If the input graph is disconnected, igraph will decompose it
* first into its subgraphs, lay out the subgraphs one by one
* using the appropriate submatrices of the distance matrix, and
* then merge the layouts using \ref igraph_layout_merge_dla.
* Since \ref igraph_layout_merge_dla works for 2D layouts only,
* you cannot run the MDS layout on disconnected graphs for
* more than two dimensions.
*
* </para><para>
* Warning: if the graph is symmetric to the exchange of two vertices
* (as is the case with leaves of a tree connecting to the same parent),
* classical multidimensional scaling may assign the same coordinates to
* these vertices.
*
* \param graph A graph object.
* \param res Pointer to an initialized matrix object. This will
* contain the result and will be resized if needed.
* \param dist The distance matrix. It must be symmetric and this
* function does not check whether the matrix is indeed
* symmetric. Results are unspecified if you pass a non-symmetric
* matrix here. You can set this parameter to null; in this
* case, the shortest path lengths between vertices will be
* used as distances.
* \param dim The number of dimensions in the embedding space. For
* 2D layouts, supply 2 here.
* \param options This argument is currently ignored, it was used for
* ARPACK, but LAPACK is used now for calculating the eigenvectors.
* \return Error code.
*
* Added in version 0.6.
*
* </para><para>
* Time complexity: usually around O(|V|^2 dim).
*/
int igraph_layout_mds(const igraph_t* graph, igraph_matrix_t *res,
const igraph_matrix_t *dist, long int dim,
igraph_arpack_options_t *options) {
long int i, no_of_nodes = igraph_vcount(graph);
igraph_matrix_t m;
igraph_bool_t conn;
RNG_BEGIN();
/* Check the distance matrix */
if (dist && (igraph_matrix_nrow(dist) != no_of_nodes ||
igraph_matrix_ncol(dist) != no_of_nodes)) {
IGRAPH_ERROR("invalid distance matrix size", IGRAPH_EINVAL);
}
/* Check the number of dimensions */
if (dim <= 1) {
IGRAPH_ERROR("dim must be positive", IGRAPH_EINVAL);
}
if (dim > no_of_nodes) {
IGRAPH_ERROR("dim must be less than the number of nodes", IGRAPH_EINVAL);
}
/* Copy or obtain the distance matrix */
if (dist == 0) {
IGRAPH_CHECK(igraph_matrix_init(&m, no_of_nodes, no_of_nodes));
IGRAPH_FINALLY(igraph_matrix_destroy, &m);
IGRAPH_CHECK(igraph_shortest_paths(graph, &m,
igraph_vss_all(), igraph_vss_all(), IGRAPH_ALL));
} else {
IGRAPH_CHECK(igraph_matrix_copy(&m, dist));
IGRAPH_FINALLY(igraph_matrix_destroy, &m);
/* Make sure that the diagonal contains zeroes only */
for (i = 0; i < no_of_nodes; i++) {
MATRIX(m, i, i) = 0.0;
}
}
/* Check whether the graph is connected */
IGRAPH_CHECK(igraph_is_connected(graph, &conn, IGRAPH_WEAK));
if (conn) {
/* Yes, it is, just do the MDS */
IGRAPH_CHECK(igraph_i_layout_mds_single(graph, res, &m, dim));
} else {
/* The graph is not connected, lay out the components one by one */
igraph_vector_ptr_t layouts;
igraph_vector_t comp, vertex_order;
igraph_t subgraph;
igraph_matrix_t *layout;
igraph_matrix_t dist_submatrix;
igraph_bool_t *seen_vertices;
long int j, n, processed_vertex_count = 0;
IGRAPH_VECTOR_INIT_FINALLY(&comp, 0);
IGRAPH_VECTOR_INIT_FINALLY(&vertex_order, no_of_nodes);
IGRAPH_CHECK(igraph_vector_ptr_init(&layouts, 0));
IGRAPH_FINALLY(igraph_vector_ptr_destroy_all, &layouts);
igraph_vector_ptr_set_item_destructor(&layouts, (igraph_finally_func_t*)igraph_matrix_destroy);
IGRAPH_CHECK(igraph_matrix_init(&dist_submatrix, 0, 0));
IGRAPH_FINALLY(igraph_matrix_destroy, &dist_submatrix);
seen_vertices = igraph_Calloc(no_of_nodes, igraph_bool_t);
if (seen_vertices == 0) {
IGRAPH_ERROR("cannot calculate MDS layout", IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(igraph_free, seen_vertices);
for (i = 0; i < no_of_nodes; i++) {
if (seen_vertices[i]) {
continue;
}
/* This is a vertex whose component we did not lay out so far */
IGRAPH_CHECK(igraph_subcomponent(graph, &comp, i, IGRAPH_ALL));
/* Take the subgraph */
IGRAPH_CHECK(igraph_induced_subgraph(graph, &subgraph, igraph_vss_vector(&comp),
IGRAPH_SUBGRAPH_AUTO));
IGRAPH_FINALLY(igraph_destroy, &subgraph);
/* Calculate the submatrix of the distances */
IGRAPH_CHECK(igraph_matrix_select_rows_cols(&m, &dist_submatrix,
&comp, &comp));
/* Allocate a new matrix for storing the layout */
layout = igraph_Calloc(1, igraph_matrix_t);
if (layout == 0) {
IGRAPH_ERROR("cannot calculate MDS layout", IGRAPH_ENOMEM);
}
IGRAPH_FINALLY(igraph_free, layout);
IGRAPH_CHECK(igraph_matrix_init(layout, 0, 0));
IGRAPH_FINALLY(igraph_matrix_destroy, layout);
/* Lay out the subgraph */
IGRAPH_CHECK(igraph_i_layout_mds_single(&subgraph, layout, &dist_submatrix, dim));
/* Store the layout */
IGRAPH_CHECK(igraph_vector_ptr_push_back(&layouts, layout));
IGRAPH_FINALLY_CLEAN(2); /* ownership of layout taken by layouts */
/* Free the newly created subgraph */
igraph_destroy(&subgraph);
IGRAPH_FINALLY_CLEAN(1);
/* Mark all the vertices in the component as visited */
n = igraph_vector_size(&comp);
for (j = 0; j < n; j++) {
seen_vertices[(long int)VECTOR(comp)[j]] = 1;
VECTOR(vertex_order)[(long int)VECTOR(comp)[j]] = processed_vertex_count++;
}
}
/* Merge the layouts - reusing dist_submatrix here */
IGRAPH_CHECK(igraph_layout_merge_dla(0, &layouts, &dist_submatrix));
/* Reordering the rows of res to match the original graph */
IGRAPH_CHECK(igraph_matrix_select_rows(&dist_submatrix, res, &vertex_order));
igraph_free(seen_vertices);
igraph_matrix_destroy(&dist_submatrix);
igraph_vector_ptr_destroy_all(&layouts);
igraph_vector_destroy(&vertex_order);
igraph_vector_destroy(&comp);
IGRAPH_FINALLY_CLEAN(5);
}
RNG_END();
igraph_matrix_destroy(&m);
IGRAPH_FINALLY_CLEAN(1);
return IGRAPH_SUCCESS;
}
/**
* \function igraph_layout_bipartite
* Simple layout for bipartite graphs
*
* The layout is created by first placing the vertices in two rows,
* according to their types. Then the positions within the rows are
* optimized to minimize edge crossings, by calling \ref
* igraph_layout_sugiyama().
*
* \param graph The input graph.
* \param types A boolean vector containing ones and zeros, the vertex
* types. Its length must match the number of vertices in the graph.
* \param res Pointer to an initialized matrix, the result, the x and
* y coordinates are stored here.
* \param hgap The preferred minimum horizontal gap between vertices
* in the same layer (i.e. vertices of the same type).
* \param vgap The distance between layers.
* \param maxiter Maximum number of iterations in the crossing
* minimization stage. 100 is a reasonable default; if you feel
* that you have too many edge crossings, increase this.
* \return Error code.
*
* \sa \ref igraph_layout_sugiyama().
*/
int igraph_layout_bipartite(const igraph_t *graph,
const igraph_vector_bool_t *types,
igraph_matrix_t *res, igraph_real_t hgap,
igraph_real_t vgap, long int maxiter) {
long int i, no_of_nodes = igraph_vcount(graph);
igraph_vector_t layers;
if (igraph_vector_bool_size(types) != no_of_nodes) {
IGRAPH_ERROR("Invalid vertex type vector size", IGRAPH_EINVAL);
}
IGRAPH_VECTOR_INIT_FINALLY(&layers, no_of_nodes);
for (i = 0; i < no_of_nodes; i++) {
VECTOR(layers)[i] = 1 - VECTOR(*types)[i];
}
IGRAPH_CHECK(igraph_layout_sugiyama(graph, res, /*extd_graph=*/ 0,
/*extd_to_orig_eids=*/ 0, &layers, hgap,
vgap, maxiter, /*weights=*/ 0));
igraph_vector_destroy(&layers);
IGRAPH_FINALLY_CLEAN(1);
return 0;
}