/**
* 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 3 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, see
* <http://www.gnu.org/licenses/>.
*
* (c) Vincenzo Nicosia 2009-2017 -- <v.nicosia@qmul.ac.uk>
*
* This file is part of NetBunch, a package for complex network
* analysis and modelling. For more information please visit:
*
* http://www.complex-networks.net/
*
* If you use this software, please add a reference to
*
* V. Latora, V. Nicosia, G. Russo
* "Complex Networks: Principles, Methods and Applications"
* Cambridge University Press (2017)
* ISBN: 9781107103184
*
***********************************************************************
*
* Compute all the simple cycles of an undirected graph provided as
* input (optionally, only up to a certain length) using the Johnson
* algorithm.
*
* References:
*
* [1] D. B. Johnson. "Finding All the Elementary Circuits of a
* Directed Graph". SIAM J. Comput. 4 (1975), 77-84.
*
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "utils.h"
#include "stack.h"
#define FALSE 0
#define TRUE 1
void usage(char *argv[]){
printf("********************************************************************\n"
"** **\n"
"** -*- johnson_cycles -*- **\n"
"** **\n"
"** Enumerates all the elementary cycles of a graph on input. **\n"
"** **\n"
"** If 'graph_in' is equal to '-' (dash), read the edge list **\n"
"** from standard input (STDIN). **\n"
"** **\n"
"** If 'max_length' is provided, ignore all cycles whose length **\n"
"** is larger than 'max_length' **\n"
"** **\n"
"** The program prints on output a row for each cycle length, **\n"
"** in the format: **\n"
"** **\n"
"** 2 N_2 **\n"
"** 3 N_3 **\n"
"** 4 N_4 **\n"
"** 5 N_5 **\n"
"** .... **\n"
"** **\n"
"** where 2, 3, 4, 5, ... is the length of the cycle, and N_2, **\n"
"** N_3, N_4, N_5, ... is the number of cycles of that length. **\n"
"** Notice that in undirected graphs each cycle is reported **\n"
"** (and counted) twice, once for each possible direction in **\n"
"** which the cycle can be traversed. **\n"
"** **\n"
"** If the third parameter is equal to 'SHOW', the program will **\n"
"** dump all the found cycles on STDERR, using the format: **\n"
"** **\n"
"** N_(l-1) N_(l-2)... N_0 **\n"
"** **\n"
"** where N_(l-i) is the i-th node in the cycle starting at node **\n"
"** N_0. Notice that according to this notation we report only **\n"
"** 'l' nodes belonging to the cycle, in reversed traversal **\n"
"** order, avoiding to report the starting node (N_0) twice. **\n"
"** **\n"
"** The program shows also cycles of length '2', that in **\n"
"** undirected graphs correspond to the degenerate cycle obtained **\n"
"** by traversing the same edge in the two directions. **\n"
"** **\n"
"********************************************************************\n"
" This is Free Software - You can use and distribute it under \n"
" the terms of the GNU General Public License, version 3 or later\n\n"
" Please visit http://www.complex-networks.net for more information\n\n"
" (c) Vincenzo Nicosia 2009-2017 (v.nicosia@qmul.ac.uk)\n"
"********************************************************************\n\n"
);
printf("Usage: %s <graph_in> [<max_length> [SHOW]]\n", argv[0]);
}
typedef struct{
unsigned int num;
double *c;
} count_t;
/**
*
* This structure contains all the information needed by the algorithm
* to run:
*
* - N is the number of nodes
*
* - K is twice the number of edges
*
* - J_slap and r_slap are the usual SLAP representation of the
* adjacency matrix
*
* - B_len is an array with the length of the set B(u) for each node u
*
* - B is a vector of character of the same size of J_slap, which
* stores the values of the indicator vector B(u) for each node u
*
* - blocked is an indicator vector of size N
*
* - v is the stack of nodes
*
* - stats maintains the counts for all the cycle lengths
*
*/
typedef struct{
unsigned int N;
unsigned int K;
unsigned int *J_slap;
unsigned int *r_slap;
my_stack_t *stack;
int max_depth;
count_t stats;
unsigned int *B_len;
unsigned int *B;
char *blocked;
char show;
} algo_info_t;
void unblock(algo_info_t *info, unsigned int u){
unsigned int w, idx;
info->blocked[u] = FALSE;
while(info->B_len[u] >0){
info->B_len[u] -= 1;
idx = info->r_slap[u] + info->B_len[u];
w = info->B[idx];
if (info->blocked[w] == TRUE){
unblock(info, w);
}
}
}
/* add v to B(w) */
void set_B(algo_info_t *info, unsigned int v, unsigned int w){
unsigned int idx, i;
idx = info->r_slap[w];
for(i=0; i <info->B_len[w]; i ++){
if (info->B[idx + i] == v)
return;
}
info->B[idx + info->B_len[w]] = v;
info->B_len[w] += 1;
}
char circuit(algo_info_t *info, unsigned int s, unsigned int v){
unsigned int i, w;
char f = FALSE;
if (v < s)
return FALSE;
/* check if maximum depth has been reached */
if(stack_size(info->stack) == info->max_depth){
return TRUE;
}
stack_push(info->stack, v);
info->blocked[v] = TRUE;
for(i=info->r_slap[v]; i< info->r_slap[v+1]; i++){
w = info->J_slap[i];
if (w == s){
/* output a cycle here */
if (info->show)
stack_dump(info->stack, stderr);
info->stats.c[stack_size(info->stack)] += 1;
f = TRUE;
}
if (w < s)
continue; /* We only consider nodes with higher IDs */
else
if (! (info->blocked[w])){
if (circuit(info, s, w) == TRUE)
f = TRUE;
}
}
if (f == TRUE) {
unblock(info, v);
}
else{
for(i=info->r_slap[v]; i < info->r_slap[v+1]; i++){
w = info->J_slap[i];
if (w > s)
set_B(info, v, w); /* add v to B(w) only if v is not already in B(w) */
}
}
stack_pop(info->stack, &w);
return f;
}
void johnson_cycles(algo_info_t *info){
unsigned int i, N;
N = info->N;
for(i=0; i<N; i++){
/* clear the vector of blocked nodes */
memset(info->blocked, 0, N * sizeof(char));
/* clear the indicator vector B */
/* we reset the length of each set B*/
memset(info->B_len, 0, N * sizeof(unsigned int));
circuit(info, i, i);
}
}
int main(int argc, char *argv[]){
algo_info_t info;
FILE *filein;
unsigned int i;
if (argc < 2){
usage(argv);
exit(1);
}
info.J_slap = NULL;
info.r_slap = NULL;
if(!strcmp(argv[1], "-")){
filein = stdin;
}
else{
filein = openfile_or_exit(argv[1], "r", 2);
}
read_slap(filein, &(info.K), &(info.N), &(info.J_slap), &(info.r_slap));
fclose(filein);
info.B = malloc(info.r_slap[info.N] * sizeof(unsigned int)); /* This is the indicator vector B(u) */
info.B_len = malloc(info.N * sizeof(unsigned int)); /* lengths of the B(u) */
info.blocked = malloc(info.N * sizeof(char)); /* This is the indicator vector
of blocked nodes */
info.stack = malloc(sizeof(my_stack_t));
stack_create(info.stack);
info.show = FALSE;
if (argc > 2){
info.max_depth = atoi(argv[2]);
if (argc > 3){
if (!my_strcasecmp(argv[3], "show")){
info.show = TRUE;
}
}
}
else{
info.max_depth = info.N;
}
/* Here we initialise the counters */
info.stats.num = info.max_depth + 1;
info.stats.c = malloc(info.stats.num * sizeof(double));
/* Here we set the counters to zero */
memset(info.stats.c, 0, info.stats.num * sizeof(double));
johnson_cycles(&info);
for(i=2; i<info.stats.num ; i ++){
printf("%d %1.20g\n", i, info.stats.c[i]);
}
free(info.B);
free(info.B_len);
free(info.blocked);
free(info.stack->v);
free(info.stack);
free(info.stats.c);
free(info.r_slap);
free(info.J_slap);
}