First commit on github -- NetBunch 1.0

1 files changed, 328 insertions, 0 deletions

diff --git a/src/johnson/johnson_cycles.c b/src/johnson/johnson_cycles.c
new file mode 100644
index 0000000..5135b20
--- /dev/null
+++ b/src/johnson/johnson_cycles.c
@@ -0,0 +1,328 @@
+/**
+ *  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);
+}