/* -*- mode: C -*- */
/*
IGraph R package.
Copyright (C) 2003, 2004, 2005, 2006 Gabor Csardi <csardi@rmki.kfki.hu>
MTA RMKI, Konkoly-Thege Miklos st. 29-33, Budapest 1121, Hungary
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
02110-1301 USA
*/
#include "igraph.h"
#include "random.h"
#include "memory.h"
#include <math.h>
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);
/**
* \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 graph object. This will
* contain the result and will be resized in 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 graph object. This will
* contain the result and will be resized in needed.
* \return Error code.
*
* Time complexity: O(|V|), the
* number of vertices.
*/
int igraph_layout_circle(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));
for (i=0; i<no_of_nodes; i++) {
igraph_real_t phi=2*M_PI/no_of_nodes*i;
MATRIX(*res, i, 0)=cos(phi);
MATRIX(*res, i, 1)=sin(phi);
}
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, the will be
* stored here. It will be resized.
* \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_fruchterman_reingold
* \brief Places the vertices on a plane according to the
* Fruchterman-Reingold algorithm.
*
* </para><para>
* This is a force-directed layout, see Fruchterman, T.M.J. and
* Reingold, E.M.: Graph Drawing by Force-directed Placement.
* Software -- Practice and Experience, 21/11, 1129--1164,
* 1991.
* This function was ported from the SNA R package.
* \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 in needed.
* \param niter The number of iterations to do.
* \param maxdelta The maximum distance to move a vertex in an
* iteration.
* \param area The area parameter of the algorithm.
* \param coolexp The cooling exponent of the simulated annealing.
* \param repulserad Determines the radius at which
* vertex-vertex repulsion cancels out attraction of
* adjacent vertices.
* \param use_seed Logical, if true the supplied values in the
* \p res argument are used as an initial layout, if
* false a random initial layout is used.
* \return Error code.
*
* Time complexity: O(|V|^2) in each
* iteration, |V| is the number of
* vertices in the graph.
*/
int igraph_layout_fruchterman_reingold(const igraph_t *graph, igraph_matrix_t *res,
igraph_integer_t niter, igraph_real_t maxdelta,
igraph_real_t area, igraph_real_t coolexp,
igraph_real_t repulserad, igraph_bool_t use_seed) {
igraph_real_t frk,t,ded,xd,yd;
igraph_real_t rf,af;
long int i,j,k;
long int no_of_nodes=igraph_vcount(graph);
igraph_matrix_t dxdy=IGRAPH_MATRIX_NULL;
igraph_eit_t edgeit;
igraph_integer_t from, to;
IGRAPH_CHECK(igraph_matrix_resize(res, no_of_nodes, 2));
if (!use_seed) {
IGRAPH_CHECK(igraph_layout_random(graph, res));
}
IGRAPH_MATRIX_INIT_FINALLY(&dxdy, no_of_nodes, 2);
IGRAPH_CHECK(igraph_eit_create(graph, igraph_ess_all(0), &edgeit));
IGRAPH_FINALLY(igraph_eit_destroy, &edgeit);
frk=sqrt(area/no_of_nodes);
for(i=niter;i>0;i--) {
/* Report progress in approx. every 100th step */
if (i%10 == 0)
igraph_progress("Fruchterman-Reingold layout: ",
100.0-100.0*i/niter, NULL);
/* Set the temperature (maximum move/iteration) */
t=maxdelta*pow(i/(double)niter,coolexp);
/* Clear the deltas */
igraph_matrix_null(&dxdy);
/* Increment deltas for each undirected pair */
for(j=0;j<no_of_nodes;j++) {
IGRAPH_ALLOW_INTERRUPTION();
for(k=j+1;k<no_of_nodes;k++){
/* Obtain difference vector */
xd=MATRIX(*res, j, 0)-MATRIX(*res, k, 0);
yd=MATRIX(*res, j, 1)-MATRIX(*res, k, 1);
ded=sqrt(xd*xd+yd*yd); /* Get dyadic euclidean distance */
xd/=ded; /* Rescale differences to length 1 */
yd/=ded;
/* Calculate repulsive "force" */
rf=frk*frk*(1.0/ded-ded*ded/repulserad);
MATRIX(dxdy, j, 0)+=xd*rf; /* Add to the position change vector */
MATRIX(dxdy, k, 0)-=xd*rf;
MATRIX(dxdy, j, 1)+=yd*rf;
MATRIX(dxdy, k, 1)-=yd*rf;
}
}
/* Calculate the attractive "force" */
IGRAPH_EIT_RESET(edgeit);
while (!IGRAPH_EIT_END(edgeit)) {
igraph_edge(graph, IGRAPH_EIT_GET(edgeit), &from, &to);
j=from;
k=to;
xd=MATRIX(*res, j, 0)-MATRIX(*res, k, 0);
yd=MATRIX(*res, j, 1)-MATRIX(*res, k, 1);
ded=sqrt(xd*xd+yd*yd); /* Get dyadic euclidean distance */
if (ded != 0) {
xd/=ded; /* Rescale differences to length 1 */
yd/=ded;
}
af=ded*ded/frk;
MATRIX(dxdy, j, 0)-=xd*af; /* Add to the position change vector */
MATRIX(dxdy, k, 0)+=xd*af;
MATRIX(dxdy, j, 1)-=yd*af;
MATRIX(dxdy, k, 1)+=yd*af;
IGRAPH_EIT_NEXT(edgeit);
}
/* Dampen motion, if needed, and move the points */
for(j=0;j<no_of_nodes;j++){
ded=sqrt(MATRIX(dxdy, j, 0)*MATRIX(dxdy, j, 0)+
MATRIX(dxdy, j, 1)*MATRIX(dxdy, j, 1));
if(ded>t){ /* Dampen to t */
ded=t/ded;
MATRIX(dxdy, j, 0)*=ded;
MATRIX(dxdy, j, 1)*=ded;
}
MATRIX(*res, j, 0)+=MATRIX(dxdy, j, 0); /* Update positions */
MATRIX(*res, j, 1)+=MATRIX(dxdy, j, 1);
}
}
igraph_progress("Fruchterman-Reingold layout: ", 100.0, NULL);
igraph_eit_destroy(&edgeit);
igraph_matrix_destroy(&dxdy);
IGRAPH_FINALLY_CLEAN(2);
return 0;
}
/**
* \function igraph_layout_fruchterman_reingold_3d
* \brief This is the 3D version of the force based
* Fruchterman-Reingold layout (see \ref
* igraph_layout_fruchterman_reingold for the 2D version
*
* </para><para>
* This function was ported from the SNA R package.
* \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 in needed.
* \param niter The number of iterations to do.
* \param maxdelta The maximum distance to move a vertex in an
* iteration.
* \param volume The volume parameter of the algorithm.
* \param coolexp The cooling exponent of the simulated annealing.
* \param repulserad Determines the radius at which
* vertex-vertex repulsion cancels out attraction of
* adjacent vertices.
* \param use_seed Logical, if true the supplied values in the
* \p res argument are used as an initial layout, if
* false a random initial layout is used.
* \return Error code.
*
* Added in version 0.2.</para><para>
*
* Time complexity: O(|V|^2) in each
* iteration, |V| is the number of
* vertices in the graph.
*
*/
int igraph_layout_fruchterman_reingold_3d(const igraph_t *graph,
igraph_matrix_t *res,
igraph_integer_t niter, igraph_real_t maxdelta,
igraph_real_t volume, igraph_real_t coolexp,
igraph_real_t repulserad,
igraph_bool_t use_seed) {
igraph_real_t frk, t, ded, xd, yd, zd;
igraph_matrix_t dxdydz;
igraph_real_t rf, af;
long int i, j, k;
long int no_of_nodes=igraph_vcount(graph);
igraph_eit_t edgeit;
igraph_integer_t from, to;
IGRAPH_CHECK(igraph_matrix_init(&dxdydz, no_of_nodes, 3));
IGRAPH_FINALLY(igraph_matrix_destroy, &dxdydz);
IGRAPH_CHECK(igraph_matrix_resize(res, no_of_nodes, 3));
if (!use_seed) {
IGRAPH_CHECK(igraph_layout_random_3d(graph, res));
}
IGRAPH_CHECK(igraph_eit_create(graph, igraph_ess_all(0), &edgeit));
IGRAPH_FINALLY(igraph_eit_destroy, &edgeit);
frk=pow(volume/(double)no_of_nodes,1.0/3.0); /*Define the F-R constant*/
/*Run the annealing loop*/
for(i=niter;i>=0;i--){
if (i%10 == 0)
igraph_progress("3D Fruchterman-Reingold layout: ",
100.0-100.0*i/niter, NULL);
/*Set the temperature (maximum move/iteration)*/
t=maxdelta*pow(i/(double)niter,coolexp);
/*Clear the deltas*/
igraph_matrix_null(&dxdydz);
/*Increment deltas for each undirected pair*/
for(j=0;j<no_of_nodes;j++) {
IGRAPH_ALLOW_INTERRUPTION();
for(k=j+1;k<no_of_nodes;k++){
/*Obtain difference vector*/
xd=MATRIX(*res, j, 0)-MATRIX(*res, k, 0);
yd=MATRIX(*res, j, 1)-MATRIX(*res, k, 1);
zd=MATRIX(*res, j, 2)-MATRIX(*res, k, 2);
ded=sqrt(xd*xd+yd*yd+zd*zd); /*Get dyadic euclidean distance*/
if (ded != 0) {
xd/=ded; /*Rescale differences to length 1*/
yd/=ded;
zd/=ded;
}
/*Calculate repulsive "force"*/
rf=frk*frk*(1.0/ded-ded*ded/repulserad);
MATRIX(dxdydz, j, 0)+=xd*rf; /*Add to the position change vector*/
MATRIX(dxdydz, k, 0)-=xd*rf;
MATRIX(dxdydz, j, 1)+=yd*rf;
MATRIX(dxdydz, k, 1)-=yd*rf;
MATRIX(dxdydz, j, 2)+=zd*rf;
MATRIX(dxdydz, k, 2)-=zd*rf;
}
}
/*Calculate the attractive "force"*/
IGRAPH_EIT_RESET(edgeit);
while (!IGRAPH_EIT_END(edgeit)) {
igraph_edge(graph, IGRAPH_EIT_GET(edgeit), &from, &to);
j=from;
k=to;
xd=MATRIX(*res, j, 0)-MATRIX(*res, k, 0);
yd=MATRIX(*res, j, 1)-MATRIX(*res, k, 1);
zd=MATRIX(*res, j, 2)-MATRIX(*res, k, 2);
ded=sqrt(xd*xd+yd*yd+zd*zd); /*Get dyadic euclidean distance*/
if (ded != 0) {
xd/=ded; /*Rescale differences to length 1*/
yd/=ded;
zd/=ded;
}
af=ded*ded/frk;
MATRIX(dxdydz, j, 0)-=xd*af; /*Add to the position change vector*/
MATRIX(dxdydz, k, 0)+=xd*af;
MATRIX(dxdydz, j, 1)-=yd*af;
MATRIX(dxdydz, k, 1)+=yd*af;
MATRIX(dxdydz, j, 2)-=zd*af;
MATRIX(dxdydz, k, 2)+=zd*af;
IGRAPH_EIT_NEXT(edgeit);
}
/*Dampen motion, if needed, and move the points*/
for(j=0;j<no_of_nodes;j++){
ded=sqrt(MATRIX(dxdydz, j, 0)*MATRIX(dxdydz, j, 0)+
MATRIX(dxdydz, j, 1)*MATRIX(dxdydz, j, 1)+
MATRIX(dxdydz, j, 2)*MATRIX(dxdydz, j, 2));
if(ded>t){ /*Dampen to t*/
ded=t/ded;
MATRIX(dxdydz, j, 0)*=ded;
MATRIX(dxdydz, j, 1)*=ded;
MATRIX(dxdydz, j, 2)*=ded;
}
MATRIX(*res, j, 0)+=MATRIX(dxdydz, j, 0); /*Update positions*/
MATRIX(*res, j, 1)+=MATRIX(dxdydz, j, 1);
MATRIX(*res, j, 2)+=MATRIX(dxdydz, j, 2);
}
}
igraph_progress("3D Fruchterman-Reingold layout: ", 100.0, NULL);
igraph_matrix_destroy(&dxdydz);
igraph_eit_destroy(&edgeit);
IGRAPH_FINALLY_CLEAN(2);
return 0;
}
/**
* \ingroup layout
* \function igraph_layout_kamada_kawai
* \brief Places the vertices on a plane according the Kamada-Kawai
* algorithm.
*
* </para><para>
* This is a force directed layout, see Kamada, T. and Kawai, S.: An
* Algorithm for Drawing General Undirected Graphs. Information
* Processing Letters, 31/1, 7--15, 1989.
* This function was ported from the SNA R package.
* \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 niter The number of iterations to perform.
* \param sigma Sets the base standard deviation of position
* change proposals.
* \param initemp Sets the initial temperature for the annealing.
* \param coolexp The cooling exponent of the annealing.
* \param kkconst The Kamada-Kawai vertex attraction constant.
* \return Error code.
*
* Time complexity: O(|V|^2) for each
* iteration, |V| is the number of
* vertices in the graph.
*/
int igraph_layout_kamada_kawai(const igraph_t *graph, igraph_matrix_t *res,
igraph_integer_t niter, igraph_real_t sigma,
igraph_real_t initemp, igraph_real_t coolexp,
igraph_real_t kkconst) {
igraph_real_t temp, candx, candy;
igraph_real_t dpot, odis, ndis, osqd, nsqd;
long int n,i,j,k;
igraph_matrix_t elen;
/* Define various things */
n=igraph_vcount(graph);
/* Calculate elen, initial x & y */
RNG_BEGIN();
IGRAPH_CHECK(igraph_matrix_resize(res, n, 2));
IGRAPH_MATRIX_INIT_FINALLY(&elen, n, n);
IGRAPH_CHECK(igraph_shortest_paths(graph, &elen, igraph_vss_all(),
IGRAPH_ALL));
for (i=0; i<n; i++) {
MATRIX(elen, i, i) = sqrt(n);
MATRIX(*res, i, 0) = RNG_NORMAL(0, n/4.0);
MATRIX(*res, i, 1) = RNG_NORMAL(0, n/4.0);
}
/*Perform the annealing loop*/
temp=initemp;
for(i=0;i<niter;i++){
/* Report progress in approx. every 100th step */
if (i%10 == 0)
igraph_progress("Kamada-Kawai layout: ",
100.0*i/niter, NULL);
/*Update each vertex*/
for(j=0;j<n;j++){
IGRAPH_ALLOW_INTERRUPTION();
/*Draw the candidate via a gaussian perturbation*/
candx=RNG_NORMAL(MATRIX(*res, j, 0),sigma*temp/initemp);
candy=RNG_NORMAL(MATRIX(*res, j, 1),sigma*temp/initemp);
/*Calculate the potential difference for the new position*/
dpot=0.0;
for(k=0;k<n;k++) /*Potential differences for pairwise effects*/
if(j!=k){
odis=sqrt((MATRIX(*res, j, 0)-MATRIX(*res, k, 0))*
(MATRIX(*res, j, 0)-MATRIX(*res, k, 0))+
(MATRIX(*res, j, 1)-MATRIX(*res, k, 1))*
(MATRIX(*res, j, 1)-MATRIX(*res, k, 1)));
ndis=sqrt((candx-MATRIX(*res, k, 0))*(candx-MATRIX(*res, k, 0))+
(candy-MATRIX(*res, k, 1))*(candy-MATRIX(*res, k, 1)));
osqd=(odis-MATRIX(elen, j, k))*(odis-MATRIX(elen, j, k));
nsqd=(ndis-MATRIX(elen, j, k))*(ndis-MATRIX(elen, j, k));
dpot+=kkconst*(osqd-nsqd)/(MATRIX(elen, j, k)*MATRIX(elen, j, k));
}
/*Make a keep/reject decision*/
if(log(RNG_UNIF(0.0,1.0))<dpot/temp){
MATRIX(*res, j, 0)=candx;
MATRIX(*res, j, 1)=candy;
}
}
/*Cool the system*/
temp*=coolexp;
}
igraph_progress("Kamada-Kawai layout: ", 100.0, NULL);
RNG_END();
igraph_matrix_destroy(&elen);
IGRAPH_FINALLY_CLEAN(1);
return 0;
}
/**
* \function igraph_layout_kamada_kawai_3d
* \brief 3D version of the force based Kamada-Kawai layout, the pair
* of the \ref igraph_layout_kamada_kawai 2D layout generator
*
* </para><para>
* This function was ported from the SNA R package.
* \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 niter The number of iterations to perform.
* \param sigma Sets the base standard deviation of position
* change proposals.
* \param initemp Sets the initial temperature for the annealing.
* \param coolexp The cooling exponent of the annealing.
* \param kkconst The Kamada-Kawai vertex attraction constant.
* \return Error code.
*
* Added in version 0.2.</para><para>
*
* Time complexity: O(|V|^2) for each
* iteration, |V| is the number of
* vertices in the graph.
*/
int igraph_layout_kamada_kawai_3d(const igraph_t *graph, igraph_matrix_t *res,
igraph_integer_t niter, igraph_real_t sigma,
igraph_real_t initemp, igraph_real_t coolexp,
igraph_real_t kkconst) {
igraph_real_t temp, candx, candy, candz;
igraph_real_t dpot, odis, ndis, osqd, nsqd;
long int i,j,k;
long int no_of_nodes=igraph_vcount(graph);
igraph_matrix_t elen;
RNG_BEGIN();
IGRAPH_CHECK(igraph_matrix_resize(res, no_of_nodes, 3));
IGRAPH_MATRIX_INIT_FINALLY(&elen, no_of_nodes, no_of_nodes);
IGRAPH_CHECK(igraph_shortest_paths(graph, &elen, igraph_vss_all(),
IGRAPH_ALL));
temp=initemp;
for(i=0;i<niter;i++){
if (i%10 == 0)
igraph_progress("3D Kamada-Kawai layout: ",
100.0*i/niter, NULL);
/*Update each vertex*/
for(j=0;j<no_of_nodes;j++){
IGRAPH_ALLOW_INTERRUPTION();
/*Draw the candidate via a gaussian perturbation*/
candx=RNG_NORMAL(MATRIX(*res, j, 0) ,sigma*temp/initemp);
candy=RNG_NORMAL(MATRIX(*res, j, 1) ,sigma*temp/initemp);
candz=RNG_NORMAL(MATRIX(*res, j, 2) ,sigma*temp/initemp);
/*Calculate the potential difference for the new position*/
dpot=0.0;
for(k=0;k<no_of_nodes;k++) /*Potential differences for pairwise effects*/
if(j!=k){
odis=sqrt((MATRIX(*res, j, 0)-MATRIX(*res, k, 0))*
(MATRIX(*res, j, 0)-MATRIX(*res, k, 0))+
(MATRIX(*res, j, 1)-MATRIX(*res, k, 1))*
(MATRIX(*res, j, 1)-MATRIX(*res, k, 1))+
(MATRIX(*res, j, 2)-MATRIX(*res, k, 2))*
(MATRIX(*res, j, 2)-MATRIX(*res, k, 2)));
ndis=sqrt((candx-MATRIX(*res, k, 0))*(candx-MATRIX(*res, k, 0))+
(candy-MATRIX(*res, k, 1))*(candy-MATRIX(*res, k, 1))+
(candz-MATRIX(*res, k, 2))*(candz-MATRIX(*res, k, 2)));
osqd=(odis-MATRIX(elen, j, k))*(odis-MATRIX(elen, j, k));
nsqd=(ndis-MATRIX(elen, j, k))*(ndis-MATRIX(elen, j, k));
dpot+=kkconst*(osqd-nsqd)/(MATRIX(elen, j, k)*MATRIX(elen, j, k));
}
/*Make a keep/reject decision*/
if(log(RNG_UNIF(0.0,1.0))<dpot/temp){
MATRIX(*res, j, 0)=candx;
MATRIX(*res, j, 1)=candy;
MATRIX(*res, j, 2)=candz;
}
}
/*Cool the system*/
temp*=coolexp;
}
igraph_progress("3D Kamada-Kawai layout: ", 100.0, NULL);
RNG_END();
igraph_matrix_destroy(&elen);
IGRAPH_FINALLY_CLEAN(1);
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) {
/* TODO */
return 0;
}
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://bioinformatics.icmb.utexas.edu/lgl/). But unlike LGL, this
* version uses a Fruchterman-Reingold style simulated annealing
* algorithm for placing the vertices. The speedup is achived 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.
* \param maxdelta The maximum length of the move allowed for a vertex
* in a single iteration.
* \param area This parameter gives the area of the square on which
* the vertices will be placed.
* \param coolexp The cooling exponent.
* \param repulserad Determines the radius at which vertex-vertex
* repulsion cancels out attraction of adjacenct vertices.
* \param cellsize The size of the grid cells, one side of the
* square.
* \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_bfs(&mst, 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_multiply(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=VECTOR(layers)[actlayer];
for (i=VECTOR(layers)[actlayer-1]; i<VECTOR(layers)[actlayer]; i++) {
long int vid=VECTOR(vids)[i];
long int par=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, VECTOR(vids)[j], sx+rx, sy+ry);
j++;
}
}
RNG_END();
/*-----------------------------------------*/
/* Step 2: add the edges of the next layer */
/*-----------------------------------------*/
for (j=VECTOR(layers)[actlayer]; j<VECTOR(layers)[actlayer+1]; j++) {
long int vid=VECTOR(vids)[j];
long int k;
IGRAPH_ALLOW_INTERRUPTION();
IGRAPH_CHECK(igraph_adjacent(graph, &eids, vid, IGRAPH_ALL));
for (k=0;k<igraph_vector_size(&eids);k++) {
long int eid=VECTOR(eids)[k];
igraph_integer_t from, to;
igraph_edge(graph, 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 j;
igraph_real_t t=maxdelta*pow((maxit-it)/(double)maxit, coolexp);
long int vid, nei;
/* init */
igraph_vector_null(&forcex);
igraph_vector_null(&forcey);
maxchange=0;
/* attractive "forces" along the edges */
for (j=0; j<igraph_vector_size(&edges); j++) {
igraph_integer_t from, to;
igraph_real_t xd, yd, dist, force;
IGRAPH_ALLOW_INTERRUPTION();
igraph_edge(graph, VECTOR(edges)[j], &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 (j=0; j<VECTOR(layers)[actlayer+1]; j++) {
long int vid=VECTOR(vids)[j];
igraph_real_t fx=VECTOR(forcex)[vid];
igraph_real_t fy=VECTOR(forcey)[vid];
igraph_real_t ded=sqrt(fx*fx+fy*fy);
if (ded > t) {
ded=t/ded;
fx*=ded; fy *=ded;
}
igraph_2dgrid_move(&grid, vid, fx, fy);
if (fx > maxchange) { maxchange=fx; }
if (fy > maxchange) { maxchange=fy; }
}
it++;
/* printf("%li iterations, maxchange: %f\n", it, (double)maxchange); */
}
}
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;
}
/**
* \function igraph_layout_grid_fruchterman_reingold
* \brief Force based layout generator for large graphs.
*
* </para><para>
* This algorithm is the same as the Fruchterman-Reingold layout
* generator, but it partitions the 2d space to a grid and and vertex
* repulsion is calculated only for vertices nearby.
*
* \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 niter Number of iterations to perform.
* \param maxdelta Maximum distance to move a vertex in an iteration.
* \param area The area of the square on which the vertices will be
* placed.
* \param coolexp The cooling exponent.
* \param repulserad Determines the radius at which vertex-vertex
* repulsion cancels out attraction of adjacenct vertices.
* \param cellsize The size of the grid cells.
* \param use_seed Logical, if true, the coordinates passed in \p res
* (should have the appropriate size) will be used for the first
* iteration.
* \return Error code.
*
* Added in version 0.2.</para><para>
*
* Time complexity: ideally (constant number of vertices in each cell)
* O(niter*(|V|+|E|)), in the worst case O(niter*(|V|^2+|E|)).
*/
int igraph_layout_grid_fruchterman_reingold(const igraph_t *graph,
igraph_matrix_t *res,
igraph_integer_t niter, igraph_real_t maxdelta,
igraph_real_t area, igraph_real_t coolexp,
igraph_real_t repulserad,
igraph_real_t cellsize,
igraph_bool_t use_seed) {
long int no_of_nodes=igraph_vcount(graph);
long int no_of_edges=igraph_ecount(graph);
igraph_2dgrid_t grid;
igraph_vector_t forcex;
igraph_vector_t forcey;
long int i, it=0;
igraph_2dgrid_iterator_t vidit;
igraph_real_t frk=sqrt(area/no_of_nodes);
IGRAPH_CHECK(igraph_matrix_resize(res, no_of_nodes, 2));
IGRAPH_VECTOR_INIT_FINALLY(&forcex, no_of_nodes);
IGRAPH_VECTOR_INIT_FINALLY(&forcey, no_of_nodes);
/* initial layout */
if (!use_seed) {
IGRAPH_CHECK(igraph_layout_random(graph, res));
igraph_matrix_multiply(res, sqrt(area/M_PI));
}
/* make grid */
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 vertices on grid */
for (i=0; i<no_of_nodes; i++) {
igraph_2dgrid_add2(&grid, i);
}
while (it<niter) {
long int j;
igraph_real_t t=maxdelta*pow((niter-it)/(double)niter, coolexp);
long int vid, nei;
/* Report progress */
if (it%10 == 0) {
igraph_progress("Grid based Fruchterman-Reingold layout: ",
(100.0*it)/niter, NULL);
}
igraph_vector_null(&forcex);
igraph_vector_null(&forcey);
/* attraction */
for (j=0; j<no_of_edges; j++) {
igraph_integer_t from, to;
igraph_real_t xd, yd, dist, force;
igraph_edge(graph, j, &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;
}
/* repulsion */
igraph_2dgrid_reset(&grid, &vidit);
while ( (vid=igraph_2dgrid_next(&grid, &vidit)-1) != -1) {
IGRAPH_ALLOW_INTERRUPTION();
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) {
if (dist==0) { dist=1e-6; };
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;
}
}
}
/* update */
for (j=0; j<no_of_nodes; j++) {
long int vid=j;
igraph_real_t fx=VECTOR(forcex)[vid];
igraph_real_t fy=VECTOR(forcey)[vid];
igraph_real_t ded=sqrt(fx*fx+fy*fy);
if (ded > t) {
ded=t/ded;
fx*=ded; fy *=ded;
}
igraph_2dgrid_move(&grid, vid, fx, fy);
}
it++;
} /* it<niter */
igraph_progress("Grid based Fruchterman-Reingold layout: ",
100.0, NULL);
igraph_vector_destroy(&forcex);
igraph_vector_destroy(&forcey);
igraph_2dgrid_destroy(&grid);
IGRAPH_FINALLY_CLEAN(3);
return 0;
}
/* Internal structure for Reingold-Tilford layout */
struct igraph_i_reingold_tilford_vertex {
int parent; /* Parent node index */
int level; /* Level of the node */
igraph_real_t offset; /* X offset from parent node */
int left_contour; /* Next left node of the contour
of the subtree rooted at this node */
int right_contour; /* Next right node of the contour
of the subtree rooted at this node */
};
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);
/**
* \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 root The index of the root vertex.
* \return Error code.
*
* Added in version 0.2.
*
* </para><para>
* TODO: decompose and merge for not fully connected graphs
* 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)
*/
int igraph_layout_reingold_tilford(const igraph_t *graph,
igraph_matrix_t *res, long int root) {
long int no_of_nodes=igraph_vcount(graph);
long int i, n, j;
igraph_dqueue_t q=IGRAPH_DQUEUE_NULL;
igraph_i_adjlist_t allneis;
igraph_vector_t *neis;
struct igraph_i_reingold_tilford_vertex *vdata;
if (root<0 || root>=no_of_nodes) {
IGRAPH_ERROR("invalid vertex id", IGRAPH_EINVVID);
}
IGRAPH_CHECK(igraph_matrix_resize(res, no_of_nodes, 2));
IGRAPH_DQUEUE_INIT_FINALLY(&q, 100);
IGRAPH_CHECK(igraph_i_adjlist_init(graph, &allneis, IGRAPH_ALL));
IGRAPH_FINALLY(igraph_i_adjlist_destroy, &allneis);
vdata=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[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=igraph_dqueue_pop(&q);
long int actdist=igraph_dqueue_pop(&q);
neis=igraph_i_adjlist_get(&allneis, actnode);
n=igraph_vector_size(neis);
for (j=0; j<n; j++) {
long int neighbor=VECTOR(*neis)[j];
if (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_i_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, n;
MATRIX(*res, node, 0) = xpos;
for (i=0, n=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 left contour */
if (vdata[lnode].right_contour >= 0) {
lnode = vdata[lnode].right_contour;
loffset += vdata[lnode].offset;
} else {
lnode = vdata[lnode].left_contour;
if (lnode >= 0) loffset += vdata[lnode].offset;
}
/* Step to the next level on the right contour */
if (vdata[rnode].left_contour >= 0) {
rnode = vdata[rnode].left_contour;
roffset += vdata[rnode].offset;
} else if (vdata[rnode].right_contour >= 0) {
rnode = vdata[rnode].right_contour;
roffset += vdata[rnode].offset;
} else {
/* Right subtree ended. If the left subtree is deeper, the
* left contour will continue on the left subtree. We can
* use only lnode here, since it has already been advanced
* to the next level in the previous if statement */
vdata[rnode].left_contour = lnode;
vdata[rnode].right_contour = lnode;
rnode = -1;
/* printf(" Right subtree ended, continuing its contours to %d\n", vdata[rnode].left_contour); */
}
/* 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;
avg = (avg*j)/(j+1) + rootsep/(j+1);
leftrootidx=j;
leftroot=i;
} else {
leftrootidx=j;
leftroot=i;
vdata[node].left_contour=i;
avg = vdata[i].offset;
}
j++;
}
}
/* printf("Shifting node to be centered above children. Shift amount: %lf\n", avg); */
vdata[node].offset += avg;
for (i=0, j=0; i<vcount; i++) {
if (i == node) continue;
if (vdata[i].parent == node) vdata[i].offset -= avg;
}
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);
/**
* \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;
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);
for (i=0; i<igraph_vector_ptr_size(coords); i++) {
igraph_matrix_t *mat=VECTOR(*coords)[i];
long int size=igraph_matrix_nrow(mat);
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=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=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-maxr, maxx-maxr+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;
}
}
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;
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;
}
syntax highlighted by Code2HTML, v. 0.9.1