diff --git a/src/SOCadmissible.c b/src/SOCadmissible.c
index 20e9e50..a1bf59c 100644
--- a/src/SOCadmissible.c
+++ b/src/SOCadmissible.c
@@ -1,4 +1,5 @@
-// SPDX-FileCopyrightText: 2023 Ämin Baumeler <amin@indyfac.ch> and Eleftherios-Ermis Tselentis <eleftheriosermis.tselentis@oeaw.ac.at>
+// SPDX-FileCopyrightText: 2023 Ämin Baumeler <amin@indyfac.ch> and
+// Eleftherios-Ermis Tselentis <eleftheriosermis.tselentis@oeaw.ac.at>
 //
 // SPDX-License-Identifier: GPL-3.0-or-later
 
@@ -8,10 +9,10 @@
 
 /***
  * Read graphs from  file.
- * Input `filename': 	string
- * Input `cnt':		pointer to int
- * Input `maxdim':	pointer to int
- * Output:		pointer to graph data
+ * Input `filename':   string
+ * Input `cnt':    pointer to int
+ * Input `maxdim':  pointer to int
+ * Output:    pointer to graph data
  *
  * This function sets `cnt' to the total number of graphs read from `filename',
  * sets `maxdim' to the dimension (number of nodes) of the largest graph,
@@ -20,93 +21,96 @@
  * The ith graph is at position i*(maxdim*maxdim + 1).
  ***/
 int* read_graphs(char* filename, int *cnt, int *maxdim) {
-	// Open file
-	FILE *fp;
-	fp = fopen(filename, "r");
-	if(fp == NULL)
-		return NULL;
-	// Count number of lines, allocate space and initiate counter
-	char ch;
-	int n = 0;
-	*cnt = 0;               // Total number of graphs
-	*maxdim = 0;            // Maximum number of vertices
-	while(!feof(fp)) {
-		ch = fgetc(fp);
-		switch(ch) {
-		case '\n':
-			(*cnt)++;
-			n = 0;
-			break;
-		case '{':
-		case '}':
-			*maxdim = (n > *maxdim) ? n : *maxdim;
-			n = 0;
-			break;
-		case '0':
-		case '1':
-			n++;
-			break;
-		case ',':
-		case ' ':
-		case 0xffffffff:
-			break;
-		default:
-			return NULL;	// Format error
-		}
-	}
-	rewind(fp);     // Return to start of file
+  // Open file
+  FILE *fp;
+  fp = fopen(filename, "r");
+  if (fp == NULL)
+    return NULL;
+  // Count number of lines, allocate space and initiate counter
+  char ch;
+  int n = 0;
+  *cnt = 0;               // Total number of graphs
+  *maxdim = 0;            // Maximum number of vertices
+  while (!feof(fp)) {
+    ch = fgetc(fp);
+    switch (ch) {
+      case '\n':
+        (*cnt)++;
+        n = 0;
+        break;
+      case '{':
+      case '}':
+        *maxdim = (n > *maxdim) ? n : *maxdim;
+        n = 0;
+        break;
+      case '0':
+      case '1':
+        n++;
+        break;
+      case ',':
+      case ' ':
+      case 0xffffffff:
+        break;
+      default:
+        return NULL;  // Format error
+    }
+  }
+  rewind(fp);     // Return to start of file
 
-	// Allocate space for graph data
-	const int gdatasize = *maxdim * *maxdim + 1;          // Size allocation per graph
-	const int arraysize = (1+ *cnt) * gdatasize;          // Total space to allocate
-	int *graphs = (int*)malloc(sizeof(int)*arraysize);    // Allocation
-	if(!graphs) {
-		fprintf(stderr, "ERROR: Could not allocate enough memory to store the graphs.\n");
-		return NULL;
-	}
-	// Parse file
-	int col = 0;    // Column number
-	int row = 0;    // Row number
-	int i = 0;      // Current graph index
-	n = -1;         // Graph dimension
-	int j = 0;      // Running index
-	while ((ch = fgetc(fp)) != EOF) {
-		// New line encountered
-		if(ch == '\n') {
-			if(row != n) {
-				fprintf(stderr, "ERROR: File not properly formatted.\n");
-				return NULL;
-			}
-			graphs[i*gdatasize] = n;
-			i++;
-			j = 0;
-			col = 0;
-			row = 0;
-			n = -1;
-		// Matrix entry encountered
-		} else if(ch == '0' || ch == '1') {
-			int x = ch - '0';
-			graphs[i*gdatasize + 1 + j] = x;
-			col++;
-			j++;
-		// End-of-row or end-of-matrix encountered
-		} else if(ch == '}') {
-			if(n == -1) {
-				n = col;
-			} else if(col > 0 && n != col) {
-				fprintf(stderr, "ERROR: File not properly formatted.\n");
-				return NULL;
-			}
-			// End of row
-			if(col != 0) {
-				row++;
-				col = 0;
-			}
-		}
-	}
-	// Close file
-	fclose(fp);
-	return graphs;
+  // Allocate space for graph data
+  // Size allocation per graph
+  const int gdatasize = *maxdim * *maxdim + 1;
+  // Total space to allocate
+  const int arraysize = (1+ *cnt) * gdatasize;
+  // Allocation
+  int *graphs = (int*)malloc(sizeof(int)*arraysize);
+  if (!graphs) {
+    fprintf(stderr, "ERROR: Could not allocate enough memory to store the graphs.\n");
+    return NULL;
+  }
+  // Parse file
+  int col = 0;    // Column number
+  int row = 0;    // Row number
+  int i = 0;      // Current graph index
+  n = -1;         // Graph dimension
+  int j = 0;      // Running index
+  while ((ch = fgetc(fp)) != EOF) {
+    // New line encountered
+    if (ch == '\n') {
+      if (row != n) {
+        fprintf(stderr, "ERROR: File not properly formatted.\n");
+        return NULL;
+      }
+      graphs[i*gdatasize] = n;
+      i++;
+      j = 0;
+      col = 0;
+      row = 0;
+      n = -1;
+    // Matrix entry encountered
+    } else if (ch == '0' || ch == '1') {
+      int x = ch - '0';
+      graphs[i*gdatasize + 1 + j] = x;
+      col++;
+      j++;
+    // End-of-row or end-of-matrix encountered
+    } else if (ch == '}') {
+      if (n == -1) {
+        n = col;
+      } else if (col > 0 && n != col) {
+        fprintf(stderr, "ERROR: File not properly formatted.\n");
+        return NULL;
+      }
+      // End of row
+      if (col != 0) {
+        row++;
+        col = 0;
+      }
+    }
+  }
+  // Close file
+  fclose(fp);
+  return graphs;
 }
 
 /***
@@ -114,15 +118,15 @@ int* read_graphs(char* filename, int *cnt, int *maxdim) {
  * Returns 0 if we run out of interventions, and 1 otherwise.
  ***/
 int nextintervention(int n, int *interventionslen, int *interventionidx) {
-	for(int i=0; i<n; i++) {
-		if(interventionidx[i] < interventionslen[i] - 1) {
-			interventionidx[i]++;
-			for(int j=0; j<i; j++)
-				interventionidx[j] = 0;
-			return 1;
-		}
-	}
-	return 0;
+  for (int i=0; i < n; i++) {
+    if (interventionidx[i] < interventionslen[i] - 1) {
+      interventionidx[i]++;
+      for (int j=0; j < i; j++)
+        interventionidx[j] = 0;
+      return 1;
+    }
+  }
+  return 0;
 }
 
 /***
@@ -132,18 +136,21 @@ int nextintervention(int n, int *interventionslen, int *interventionidx) {
  * If `party' has no parents, then `party' will definitely receive a 1.
  * If the fixed-point value cannot be inferred, return -1.
  ***/
-int alphapre(int party, const int *parents, const int *parentslen, const int *interventions, const int *interventionidx, int maxn) {
-	if(parentslen[party] == 0)
-		return 1;
-	const int idatasize = (maxn+1)*(maxn+1)*2;		// Max nr of interventions/party
-	for(int pidx=0; pidx<parentslen[party]; pidx++) {
-		const int parent = parents[maxn*party + pidx];
-		const int f0 = interventions[idatasize*parent + 2*interventionidx[parent] + 0];
-		const int f1 = interventions[idatasize*parent + 2*interventionidx[parent] + 1];
-		if(f0 != party && f1 != party)
-			return 0;
-	}
-	return -1;
+int alphapre(int party, const int *parents, const int *parentslen, \
+      const int *interventions, const int *interventionidx, int maxn) {
+  if (parentslen[party] == 0)
+    return 1;
+  const int idatasize = (maxn+1)*(maxn+1)*2;    // Max nr of interventions/party
+  for (int pidx = 0; pidx < parentslen[party]; pidx++) {
+    const int parent = parents[maxn*party + pidx];
+    const int f0 = interventions[idatasize*parent + \
+        2*interventionidx[parent] + 0];
+    const int f1 = interventions[idatasize*parent + \
+        2*interventionidx[parent] + 1];
+    if (f0 != party && f1 != party)
+      return 0;
+  }
+  return -1;
 }
 
 /***
@@ -154,37 +161,43 @@ int alphapre(int party, const int *parents, const int *parentslen, const int *in
  *
  * Make sure to run alphapre before running this.
  ***/
-int alpha(int *path, int pathlen, const int *parents, const int *parentslen, const int *interventions, const int *interventionidx, const int *fp, int n, int maxn) {
-	const int party = path[pathlen-1];
-	for(int i=0;i<pathlen-1;i++) {
-			if(party == path[i]) {
-				return 0;
-		}
-	}
-	const int idatasize = (maxn+1)*(maxn+1)*2;		// Max nr of interventions/party
-	// Iterate over all parents
-	// If one parent does not vote for `party', immediattely return a zero.
-	// If all parents vote for `party' --- which can only be inferred after we queried all parents --- return a one.
-	for(int pidx=0; pidx<parentslen[party]; pidx++) {
-		const int parent = parents[maxn*party + pidx];
-		const int f0 = interventions[idatasize*parent + 2*interventionidx[parent] + 0];
-		const int f1 = interventions[idatasize*parent + 2*interventionidx[parent] + 1];
-		if(f0 == party && f1 == party)
-			continue;
-		// Re-use already computed fixed points, else enter recursion
-		int val;
-		if(fp[parent] != -1) {
-			val = fp[parent];
-		} else {
-			path[pathlen] = parent;
-			val = alpha(path, pathlen+1, parents, parentslen, interventions, interventionidx, fp, n, maxn);
-		}
-		// Evaluate the intervention on the input
-		const int muOfVal = (val == 0) ? f0 : f1;
-		if(muOfVal != party)
-			return 0;
-	}
-	return 1;
+int alpha(int *path, int pathlen, const int *parents, const int *parentslen, \
+      const int *interventions, const int *interventionidx, const int *fp, \
+      int n, int maxn) {
+  const int party = path[pathlen-1];
+  for (int i = 0; i < pathlen-1; i++) {
+      if (party == path[i]) {
+        return 0;
+    }
+  }
+  const int idatasize = (maxn+1)*(maxn+1)*2;    // Max nr of interventions/party
+  // Iterate over all parents
+  // If one parent does not vote for `party', immediattely return a zero.
+  // If all parents vote for `party' --- which can only be inferred after we
+  //   queried all parents --- return a one.
+  for (int pidx = 0; pidx < parentslen[party]; pidx++) {
+    const int parent = parents[maxn*party + pidx];
+    const int f0 = interventions[idatasize*parent + 2*interventionidx[parent] \
+        + 0];
+    const int f1 = interventions[idatasize*parent + 2*interventionidx[parent] \
+        + 1];
+    if (f0 == party && f1 == party)
+      continue;
+    // Re-use already computed fixed points, else enter recursion
+    int val;
+    if (fp[parent] != -1) {
+      val = fp[parent];
+    } else {
+      path[pathlen] = parent;
+      val = alpha(path, pathlen+1, parents, parentslen, interventions, \
+          interventionidx, fp, n, maxn);
+    }
+    // Evaluate the intervention on the input
+    const int muOfVal = (val == 0) ? f0 : f1;
+    if (muOfVal != party)
+      return 0;
+  }
+  return 1;
 }
 
 /***
@@ -192,12 +205,12 @@ int alpha(int *path, int pathlen, const int *parents, const int *parentslen, con
  * Helper function for pretty printing.
  ***/
 int ceillog10(int x) {
-	int res = 0;
-	while(x > 0) {
-		x /= 10;
-		res++;
-	}
-	return res;
+  int res = 0;
+  while (x > 0) {
+    x /= 10;
+    res++;
+  }
+  return res;
 }
 
 /***
@@ -212,193 +225,214 @@ int ceillog10(int x) {
  * Return -1 on error.
  ***/
 int test_graphs(int startidx, int len, int maxn, const int *graphs) {
-	const int gdatasize = maxn * maxn + 1;        // Size allocation per graph
-	const int idatasize = (maxn+1)*(maxn+1)*2;    // Max number of interventions per party
-	int range[maxn+1];                            // Range of intervention
-	int interventions[maxn*idatasize];            // Hold all interventions per party {f_0(0), f_0(1), g_0(0), g_0(1), ... , f_1(0), f_1(1), ...}
-	int interventionslen[maxn];                   // Number of interventions per party
-	int interventionidx[maxn];                    // Current intervention executed
-	int parents[maxn*maxn];                       // List all parents per party
-	int parentslen[maxn];                         // Number of parents per party
-	int children[maxn*maxn];                      // List all children per party
-	int childrenlen[maxn];                        // Numer of children per party
-	const int padlen = ceillog10(len);
-	time_t t0 = time(NULL);
-	const time_t tstart = t0;
-	for(int graphidx=startidx; graphidx<startidx+len; graphidx++) {
-		time_t t1 = time(NULL);
-		if(t1 > t0) {
-			const int count = graphidx - startidx + 1;
-			const float percentage = 100*(float)count/(float)len;
-			const int deltat = t1-tstart;
-			const float rate  = (float)count/(float)deltat;
-			fprintf(stderr, "\r%6.2f%%    %*d/%d (%*.2f graphs/s in %d seconds; current line %d)",percentage,padlen,count,len,padlen+3,rate,deltat,graphidx+1);
-			t0 = t1;
-		}
-		// Test graph with index `graphidx'
-		// `n' is the order of this graph
-		const int n = graphs[graphidx*gdatasize];
-		// Get parents and children
-		// Adj matrix is such that (row,col) = 1 <=> edge from row to col
-		for(int i=0; i<n; i++) {
-			parentslen[i] = 0;
-			childrenlen[i] = 0;
-		}
-		int idx = 0;
-		for(int row=0; row<n; row++) {
-			for(int col=0; col<n; col++) {
-				if(graphs[graphidx*gdatasize+1+idx] == 1) {
-					// Edge from row to col
-					children[maxn*row + childrenlen[row]] = col;
-					childrenlen[row]++;
-					parents[maxn*col + parentslen[col]] = row;
-					parentslen[col]++;
-				}
-				idx++;
-			}
-		}
-		// Fill in interventions
-		for(int party=0; party<n; party++) {
-			// This is a small trick to speed-up processing:
-			// If the party `party' has two or more children, ignore the `discard' intervention (-1) (tmp=0),
-			// else make use of the `discard' internvetion (tmp=1).
-			// Populate range of `party's interventions
-			// The range has childrenlen[party]+tmp entries
-			range[0] = -1;
-			const int tmp = (childrenlen[party]<=1) ? 1 : 0;
-			for(int j=0; j<childrenlen[party]; j++) {
-				range[j+tmp] = children[maxn*party+j];
-			}
-			// Enumerate, store, and count all possible interventions
-			// Each intervention is a tuple (mu(0), mu(1)), x corresponds to the first entry, y to the second
-			int cnt = 0;
-			for(int x=0; x<childrenlen[party]+tmp; x++) {
-				for(int y=0; y<childrenlen[party]+tmp; y++) {
-					interventions[party*idatasize+2*cnt] = range[x];
-					interventions[party*idatasize+2*cnt+1] = range[y];
-					cnt++;
-				}
-				interventionslen[party] = cnt;
-			}
-		}
-		// Loop over intervetions and verify fixed-point expression
-		int fp[n];                                       // Array to store alpha(party) values
-		int afterfp[n];                                  // \omega(\mu(fp))
-		int path[n+1];                                   // Current path
-		for(int i=0; i<n; i++) interventionidx[i] = 0;	 // Initialize intervention index to the first intervention
-		int running = 1;  // Flag to see whether we are still verifying, or we verified all interventions
-		while(running) {
-			// Reset fixed-point entries
-			for(int party=0; party<n; party++)
-				fp[party] = -1;
-			// Fill out fixed-points that are trivially computed
-			for(int party=0; party<n; party++)
-				alphapre(party, parents, parentslen, interventions, interventionidx, maxn);
-			// Invoke recursive function
-			for(int party=0; party<n; party++) {
-				path[0] = party;
-				const int val = alpha(path, 1, parents, parentslen, interventions, interventionidx, fp, n, maxn);
-				fp[party] = val;
-			}
+  const int gdatasize = maxn * maxn + 1;        // Size allocation per graph
+  const int idatasize = (maxn+1)*(maxn+1)*2;    // Max number of interventions
+                                                // per party
+  int range[maxn+1];                            // Range of intervention
+  int interventions[maxn*idatasize];            // Hold all interventions per
+                                                // party {f_0(0), f_0(1),
+                                                // g_0(0), g_0(1), ... ,
+                                                // f_1(0), f_1(1), ...}
+  int interventionslen[maxn];                   // Number of intervention per
+                                                // party
+  int interventionidx[maxn];                    // Current intervention executed
+  int parents[maxn*maxn];                       // List all parents per party
+  int parentslen[maxn];                         // Number of parents per party
+  int children[maxn*maxn];                      // List all children per party
+  int childrenlen[maxn];                        // Numer of children per party
+  const int padlen = ceillog10(len);
+  time_t t0 = time(NULL);
+  const time_t tstart = t0;
+  for (int graphidx = startidx; graphidx < startidx+len; graphidx++) {
+    time_t t1 = time(NULL);
+    if (t1 > t0) {
+      const int count = graphidx - startidx + 1;
+      const float percentage = 100*(float)count/(float)len;
+      const int deltat = t1-tstart;
+      const float rate  = (float)count/(float)deltat;
+      fprintf(stderr, "\r%6.2f%%    %*d/%d (%*.2f graphs/s in %d seconds; current line %d)", \
+          percentage, padlen, count, len, padlen+3, rate, deltat, graphidx+1);
+      t0 = t1;
+    }
+    // Test graph with index `graphidx'
+    // `n' is the order of this graph
+    const int n = graphs[graphidx*gdatasize];
+    // Get parents and children
+    // Adj matrix is such that (row,col) = 1 <=> edge from row to col
+    for (int i = 0; i < n; i++) {
+      parentslen[i] = 0;
+      childrenlen[i] = 0;
+    }
+    int idx = 0;
+    for (int row = 0; row < n; row++) {
+      for (int col = 0; col < n; col++) {
+        if (graphs[graphidx*gdatasize+1+idx] == 1) {
+          // Edge from row to col
+          children[maxn*row + childrenlen[row]] = col;
+          childrenlen[row]++;
+          parents[maxn*col + parentslen[col]] = row;
+          parentslen[col]++;
+        }
+        idx++;
+      }
+    }
+    // Fill in interventions
+    for (int party = 0; party < n; party++) {
+      // This is a small trick to speed-up processing:
+      // If the party `party' has two or more children,
+      // ignore the `discard' intervention (-1) (tmp=0),
+      // else make use of the `discard' internvetion (tmp=1).
+      // Populate range of `party's interventions
+      // The range has childrenlen[party]+tmp entries
+      range[0] = -1;
+      const int tmp = (childrenlen[party] <= 1) ? 1 : 0;
+      for (int j = 0; j < childrenlen[party]; j++) {
+        range[j+tmp] = children[maxn*party+j];
+      }
+      // Enumerate, store, and count all possible interventions
+      // Each intervention is a tuple (mu(0), mu(1)), x corresponds to the first
+      // entry, y to the second
+      int cnt = 0;
+      for (int x = 0; x < childrenlen[party]+tmp; x++) {
+        for (int y = 0; y < childrenlen[party]+tmp; y++) {
+          interventions[party*idatasize+2*cnt] = range[x];
+          interventions[party*idatasize+2*cnt+1] = range[y];
+          cnt++;
+        }
+        interventionslen[party] = cnt;
+      }
+    }
+    // Loop over intervetions and verify fixed-point expression
+    int fp[n];       // Array to store alpha(party) values
+    int afterfp[n];  // \omega(\mu(fp))
+    int path[n+1];   // Current path
+    // Initialize intervention index to the first intervention
+    for (int i = 0; i < n; i++) interventionidx[i] = 0;
+    // Flag to see whether we are still verifying, or we verified all
+    // interventions
+    int running = 1;
+    while (running) {
+      // Reset fixed-point entries
+      for (int party = 0; party < n; party++)
+        fp[party] = -1;
+      // Fill out fixed-points that are trivially computed
+      for (int party = 0; party < n; party++)
+        alphapre(party, parents, parentslen, interventions, interventionidx, \
+            maxn);
+      // Invoke recursive function
+      for (int party = 0; party < n; party++) {
+        path[0] = party;
+        const int val = alpha(path, 1, parents, parentslen, interventions, \
+            interventionidx, fp, n, maxn);
+        fp[party] = val;
+      }
 
-			// Compute \omega(\mu(fp))
-			for(int i=0; i<n; i++)
-				afterfp[i] = -1;
-			// Loop over parties
-			for(int party=0; party<n; party++) {
-				// Loop over parents of this party
-				for(int pidx=0; pidx<parentslen[party]; pidx++) {
-					const int parent = parents[maxn*party + pidx];
-					const int f0 = interventions[idatasize*parent + 2*interventionidx[parent] + 0];
-					const int f1 = interventions[idatasize*parent + 2*interventionidx[parent] + 1];
-					// If `parent' never votes for `party', this value is 0.
-					if(f0 != party && f1 != party) {
-						afterfp[party] = 0;
-						break;
-					}
-					// If `parent' always votes for `party', check the other parents.
-					if(f0 == party && f1 == party)
-						continue;
-					// Evaluate intervention on the input (fp[parent]) 
-					const int muOfVal = (fp[parent] == 0) ? f0 : f1;
-					// If `parent' does not vote for `party', this value is 0, else check the other parents.
-					if(muOfVal != party) {
-						afterfp[party] = 0;
-						break;
-					}
+      // Compute \omega(\mu(fp))
+      for (int i = 0; i < n; i++)
+        afterfp[i] = -1;
+      // Loop over parties
+      for (int party = 0; party < n; party++) {
+        // Loop over parents of this party
+        for (int pidx = 0; pidx < parentslen[party]; pidx++) {
+          const int parent = parents[maxn*party + pidx];
+          const int f0 = interventions[idatasize*parent + \
+              2*interventionidx[parent] + 0];
+          const int f1 = interventions[idatasize*parent + \
+              2*interventionidx[parent] + 1];
+          // If `parent' never votes for `party', this value is 0.
+          if (f0 != party && f1 != party) {
+            afterfp[party] = 0;
+            break;
+          }
+          // If `parent' always votes for `party', check the other parents.
+          if (f0 == party && f1 == party)
+            continue;
+          // Evaluate intervention on the input (fp[parent])
+          const int muOfVal = (fp[parent] == 0) ? f0 : f1;
+          // If `parent' does not vote for `party', this value is 0, else check
+          // the other parents.
+          if (muOfVal != party) {
+            afterfp[party] = 0;
+            break;
+          }
+        }
+        // All parents voted `party'
+        if (afterfp[party] == -1)
+          afterfp[party] = 1;
+      }
+      // Compare the output of the recursive alpha function (fp) with
+      //   \omega(\mu(fp)) (afterfp)
+      // If these two arrays differ, return the line number of the graph we
+      //   falsified
+      for (int party = 0; party < n; party++)
+         if (fp[party] != afterfp[party])
+          return graphidx+1;
 
-				}
-				// All parents voted `party'
-				if(afterfp[party] == -1)
-					afterfp[party] = 1;
-			}
-			// Compare the output of the recursive alpha function (fp) with \omega(\mu(fp)) (afterfp)
-			// If these two arrays differ, return the line number of the graph we falsified
-			for(int party=0;party<n;party++)
-	 			if(fp[party] != afterfp[party])
-					return graphidx+1;
-
-			// Next intervention
-			running = nextintervention(n, interventionslen, interventionidx);
-		}
-	}
-	const int deltat = time(NULL) - tstart;
-	const float rate  = (deltat == 0) ? len : (float)len/(float)deltat;
-	fprintf(stderr, "\r%6.2f%%    %*d/%d (%*.2f graphs/s in %d seconds; current line %d)\n",100.0,padlen,len,len,padlen+3,rate,deltat,startidx+len);
-	return 0;
+      // Next intervention
+      running = nextintervention(n, interventionslen, interventionidx);
+    }
+  }
+  const int deltat = time(NULL) - tstart;
+  const float rate  = (deltat == 0) ? len : (float)len/(float)deltat;
+  fprintf(stderr, "\r%6.2f%%    %*d/%d (%*.2f graphs/s in %d seconds; current line %d)\n", \
+      100.0, padlen, len, len, padlen+3, rate, deltat, startidx+len);
+  return 0;
 }
 
 int main(int argc, char *argv[]) {
-	// Parse command-line arguments
-	int START = 0;
-	int NUM = 0;
-	if(argc >= 3)
-		START = atoi(argv[2]);
-	if(argc >= 4) {
-		if(argv[3][0] == '+')
-			NUM = atoi(argv[3]+1);
-		else
-			NUM = atoi(argv[3])-START+1;
-	}
-	if(!(argc >= 2 && (argc < 3 || START > 0) && (argc < 4 || NUM > 0) && argc <= 4)) {
-		fprintf(stderr, "Usage: %s <filename> [<startline> [<endline> | +<count>]]\n", argv[0]);
-		fprintf(stderr, "       <filename>     File name with adjacency matrices of simple directed graphs\n");
-		fprintf(stderr, "       <startline>    Verify graphs starting from line `startline'\n");
-		fprintf(stderr, "       <endline>      Verify graphs up to and including line `endline'\n");
-		fprintf(stderr, "       +<count>       Verify `count' number of graphs\n");
-		fprintf(stderr, "\n");
-		fprintf(stderr, "[FILE FORMAT]\n");
-		fprintf(stderr, "       Each line in `filename' must contain the adjacency matrix of a simple directed graph in the format\n");
-		fprintf(stderr, "         {{a00,a01,...},{a10,a11,...},...}  where aij=1 if and only if the graph has the edge i -> j\n");
-		fprintf(stderr, "       The file `filename' may contain graphs with different order (number of vertices)\n");
-		fprintf(stderr, "\n");
-		fprintf(stderr, "This program verifies the admissibility of simple directed graphs.\n");
-		return -1;
-	}
+  // Parse command-line arguments
+  int START = 0;
+  int NUM = 0;
+  if (argc >= 3)
+    START = atoi(argv[2]);
+  if (argc >= 4) {
+    if (argv[3][0] == '+')
+      NUM = atoi(argv[3]+1);
+    else
+      NUM = atoi(argv[3])-START+1;
+  }
+  if (!(argc >= 2 && (argc < 3 || START > 0) && (argc < 4 || NUM > 0) && \
+      argc <= 4)) {
+    fprintf(stderr, "Usage: %s <filename> [<startline> [<endline> | +<count>]]\n", argv[0]);
+    fprintf(stderr, "       <filename>     File name with adjacency matrices of simple directed graphs\n");
+    fprintf(stderr, "       <startline>    Verify graphs starting from line `startline'\n");
+    fprintf(stderr, "       <endline>      Verify graphs up to and including line `endline'\n");
+    fprintf(stderr, "       +<count>       Verify `count' number of graphs\n");
+    fprintf(stderr, "\n");
+    fprintf(stderr, "[FILE FORMAT]\n");
+    fprintf(stderr, "       Each line in `filename' must contain the adjacency matrix of a simple directed graph in the format\n");
+    fprintf(stderr, "         {{a00,a01,...},{a10,a11,...},...}  where aij=1 if and only if the graph has the edge i -> j\n");
+    fprintf(stderr, "       The file `filename' may contain graphs with different order (number of vertices)\n");
+    fprintf(stderr, "\n");
+    fprintf(stderr, "This program verifies the admissibility of simple directed graphs.\n");
+    return -1;
+  }
 
-	// Read Graphs from File
-	int graphcount = 0;
-	int maxdimension = -1;
-	const int *graphs = read_graphs(argv[1], &graphcount, &maxdimension);
-	if(graphs == NULL) {
-		fprintf(stderr, "ERROR: Could not read graphs from file `%s'.\n", argv[1]);
-		return -1;
-	}
-	const int startidx = START ? START-1 : 0;		// Index of graph to start verifying
-	const int tot = NUM ? NUM : graphcount-startidx;	// Number of graphs to be verified
-	if(startidx >= graphcount || tot <= 0 || tot+startidx-1 >= graphcount) {
-		fprintf(stderr, "ERROR: Lines to be verified are out of range\n");
-		return -1;
-	}
+  // Read Graphs from File
+  int graphcount = 0;
+  int maxdimension = -1;
+  const int *graphs = read_graphs(argv[1], &graphcount, &maxdimension);
+  if (graphs == NULL) {
+    fprintf(stderr, "ERROR: Could not read graphs from file `%s'.\n", argv[1]);
+    return -1;
+  }
+  const int startidx = START ? START-1 : 0;         // Index of graph to start verifying
+  const int tot = NUM ? NUM : graphcount-startidx;  // Number of graphs to be verified
+  if (startidx >= graphcount || tot <= 0 || tot+startidx-1 >= graphcount) {
+    fprintf(stderr, "ERROR: Lines to be verified are out of range\n");
+    return -1;
+  }
 
-	// Test `tot' graphs, starting from index `startidx'
-	printf("Verifying the admissibility of %d graphs in the file `%s' (line %d to line %d)\n", tot, argv[1], startidx+1, startidx+tot);
-	const int result = test_graphs(startidx, tot, maxdimension, graphs);
-	if (result == 0)
-		printf("These graphs are admissible\n");
-	else if(result == -1)
-		fprintf(stderr, "ERROR: Something went wrong\n");
-	else 
-		printf("The function alpha does not represent the fixed point, or the graph on line %d is inadmissible\n", result);
-	return result;
+  // Test `tot' graphs, starting from index `startidx'
+  printf("Verifying the admissibility of %d graphs in the file `%s' (line %d to line %d)\n", \
+      tot, argv[1], startidx+1, startidx+tot);
+  const int result = test_graphs(startidx, tot, maxdimension, graphs);
+  if (result == 0)
+    printf("These graphs are admissible\n");
+  else if (result == -1)
+    fprintf(stderr, "ERROR: Something went wrong\n");
+  else
+    printf("The function alpha does not represent the fixed point, or the graph on line %d is inadmissible\n", \
+        result);
+  return result;
 }
diff --git a/src/SOCgen.c b/src/SOCgen.c
index eb9a1a9..9265369 100644
--- a/src/SOCgen.c
+++ b/src/SOCgen.c
@@ -1,4 +1,5 @@
-// SPDX-FileCopyrightText: 2023 Ämin Baumeler <amin@indyfac.ch> and Eleftherios-Ermis Tselentis <eleftheriosermis.tselentis@oeaw.ac.at>
+// SPDX-FileCopyrightText: 2023 Ämin Baumeler <amin@indyfac.ch> and
+// Eleftherios-Ermis Tselentis <eleftheriosermis.tselentis@oeaw.ac.at>
 //
 // SPDX-License-Identifier: GPL-3.0-or-later
 
@@ -11,186 +12,201 @@
  * Print the graph's adjacency matrix
  ***/
 void dumpgraph(int n, const int *children, const int *childrenlen) {
-	int adj[n*n];
-	for(int i=0; i<n*n; i++) adj[i] = 0;
-	for(int a=0; a<n; a++) {
-		for(int bidx=0; bidx<childrenlen[a]; bidx++) {
-			const int b = children[a*n+bidx];
-			adj[a*n+b] = 1;
-		}
-	}
-	// Print
-	printf("{");
-	for(int a=0; a<n; a++) {
-		if(a) printf(",");
-		printf("{");
-		printf("%d", adj[a*n+0]);
-		for(int b=1; b<n; b++) printf(",%d", adj[a*n+b]);
-		printf("}");
-	}
-	printf("}\n");
+  int adj[n*n];
+  for (int i = 0; i < n*n; i++) adj[i] = 0;
+  for (int a = 0; a < n; a++) {
+    for (int bidx = 0; bidx < childrenlen[a]; bidx++) {
+      const int b = children[a*n+bidx];
+      adj[a*n+b] = 1;
+    }
+  }
+  // Print
+  printf("{");
+  for (int a = 0; a < n; a++) {
+    if (a) printf(",");
+    printf("{");
+    printf("%d", adj[a*n+0]);
+    for (int b = 1; b < n; b++) printf(",%d", adj[a*n+b]);
+    printf("}");
+  }
+  printf("}\n");
 }
 
 /***
  * Print the graph as graphviz command
  ***/
 void dumpgv(int n, int graphnr, const int *children, const int *childrenlen) {
-	printf("strict digraph G%d {node[shape=point];", graphnr);
-	for(int v=0; v<n; v++) {
-		for(int i=0; i<childrenlen[v]; i++) {
-			const int w = children[v*n+i];
-			printf("G%dN%d -> G%dN%d", graphnr, v, graphnr, w);
-			if(childrenlen[v] >= 2)
-				printf(" [color=\"red\"]");
-			printf(";");
-		}
-	}
-	printf("}\n");
-	return;
+  printf("strict digraph G%d {node[shape=point];", graphnr);
+  for (int v = 0; v < n; v++) {
+    for (int i = 0; i < childrenlen[v]; i++) {
+      const int w = children[v*n+i];
+      printf("G%dN%d -> G%dN%d", graphnr, v, graphnr, w);
+      if (childrenlen[v] >= 2)
+        printf(" [color=\"red\"]");
+      printf(";");
+    }
+  }
+  printf("}\n");
+  return;
 }
 
 /***
  * Block size of our datastructure.
  ***/
 int blocksize(int n) {
-	int res=(n+1)*n;
-	for(int i=0; i<n; i++)
-		res *= n-i;
-	return res;
+  int res = (n+1)*n;
+  for (int i = 0; i < n; i++)
+    res *= n-i;
+  return res;
 }
 
 /***
  * Recrusive bi-directional graph traversal, used to check connectivity.
  ***/
-void traverse(int vertex, int *visited, int n, const int *parents, const int *parentslen, const int *children, const int *childrenlen) {
-	visited[vertex] = 1;
-	for(int i=0; i<parentslen[vertex]; i++) {
-		int u = parents[n*vertex + i];
-		if(visited[u] == 0)
-			traverse(u, visited, n, parents, parentslen, children, childrenlen);
-	}
-	for(int i=0; i<childrenlen[vertex]; i++) {
-		int u = children[n*vertex + i];
-		if(visited[u] == 0)
-			traverse(u, visited, n, parents, parentslen, children, childrenlen);
-	}
+void traverse(int vertex, int *visited, int n, const int *parents, \
+    const int *parentslen, const int *children, const int *childrenlen) {
+  visited[vertex] = 1;
+  for (int i = 0; i < parentslen[vertex]; i++) {
+    int u = parents[n*vertex + i];
+    if (visited[u] == 0)
+      traverse(u, visited, n, parents, parentslen, children, childrenlen);
+  }
+  for (int i = 0; i < childrenlen[vertex]; i++) {
+    int u = children[n*vertex + i];
+    if (visited[u] == 0)
+      traverse(u, visited, n, parents, parentslen, children, childrenlen);
+  }
 }
 
 /***
  * Return 1 if the graph is connected, and 0 otherwise.
  ***/
-int isconnected(int n, const int *parents, const int *parentslen, const int *children, const int *childrenlen) {
-	// Early check: Is there a node with 0 total degree?
-	for(int i=0; i<n; i++) if(parentslen[i]+childrenlen[i]==0) return 0;
-	// Check connectivity
-	int visited[n];
-	for(int i=0; i<n; i++) visited[i] = 0;
-	// Start traversal
-	traverse(0, visited, n, parents, parentslen, children, childrenlen);
-	for(int i=0; i<n; i++) if(!visited[i]) return 0;
-	return 1;
+int isconnected(int n, const int *parents, const int *parentslen, \
+    const int *children, const int *childrenlen) {
+  // Early check: Is there a node with 0 total degree?
+  for (int i = 0; i < n; i++) if (parentslen[i]+childrenlen[i] == 0) return 0;
+  // Check connectivity
+  int visited[n];
+  for (int i = 0; i < n; i++) visited[i] = 0;
+  // Start traversal
+  traverse(0, visited, n, parents, parentslen, children, childrenlen);
+  for (int i = 0; i < n; i++) if (!visited[i]) return 0;
+  return 1;
 }
 
 /***
  * Interpret graph number `graph' as a simple digraph, and populate the arrays `parents', `parentslen', `children', `childrenlen'
  * The number `graph' is understood as the adjacency matrix with the diagonal removed (no self loops).
  ***/
-void graphnrtolists(int n, unsigned long long int graph, int *parents, int *parentslen, int *children, int *childrenlen) {
-	for(int i=0; i<n; i++) {
-		parentslen[i] = 0;
-		childrenlen[i] = 0;
-	}
-	int row = 0;
-	int col = 1;    // We cannot start at (0,0); it does not exist
-	unsigned long long g = graph;
-	for(int i=0; i<n*(n-1); i++) {
-		int val = g%2;
-		g = g/2;
-		if(val) {
-			children[row*n+childrenlen[row]] = col;
-			childrenlen[row]++;
-			parents[col*n+parentslen[col]] = row;
-			parentslen[col]++;
-		}
-		col++;
-		if(col==row)
-			col++;
-		if(col==n) {
-			col = 0;
-			row++;
-		}
-
-	}
+void graphnrtolists(int n, unsigned long long int graph, int *parents, \
+    int *parentslen, int *children, int *childrenlen) {
+  for (int i = 0; i < n; i++) {
+    parentslen[i] = 0;
+    childrenlen[i] = 0;
+  }
+  int row = 0;
+  int col = 1;    // We cannot start at (0,0); it does not exist
+  unsigned long long g = graph;
+  for (int i = 0; i < n*(n-1); i++) {
+    int val = g%2;
+    g = g/2;
+    if (val) {
+      children[row*n+childrenlen[row]] = col;
+      childrenlen[row]++;
+      parents[col*n+parentslen[col]] = row;
+      parentslen[col]++;
+    }
+    col++;
+    if (col == row)
+      col++;
+    if (col == n) {
+      col = 0;
+      row++;
+    }
+  }
 }
 
 /***
  * Find the tail in `path' of length `pathlen' that includes a cycle, and save the canonical representation.
  * The canonical representation is the listing of the nodes of the cycle starting from the smallest one.
  ***/
-void canonical_cycle(int n, int *cycle, int *cyclelen, const int *path, int pathlen) {
-	int startidx;
-	int pathhascycle = 0;
-	const int last = path[pathlen-1];       // This is the last node
-	for(int k=0; k<pathlen-1; k++) {
-		if(path[k]==last) {             // We have found the first occurance of node `last' on the path
-		                           	// The cycle is the following part of path: path[startidx], path[startidx+1], ..., path[pathlen-2], and therefore has length pathlen-startidx-1
-			startidx = k;
-			(*cyclelen) = pathlen-startidx-1;
-			pathhascycle = 1;
-			break;
-		}
-	}
-	if(!pathhascycle) {                     // Sanity check. This should NEVER happen!
-		fprintf(stderr, "ERROR: You provided a path to `canonical_cycle' that does NOT contain a cycle as its tail! I'm returning now.\n");
-		return;
-	}
-	// Generate canonical representation
-	// First, find smallest node on cycle
-	int smallestidx = -1;
-	int smallest = n;
-	for(int k=startidx; k<pathlen-1; k++) {
-		if(smallest > path[k]) {
-			smallest = path[k];
-			smallestidx = k;
-		}
-	}
-	// Now, copy cycle to `cycle' starting from the smallest
-	int j=0;                                        // This is the index on where to write on the array `cycle'
-	for(int i=smallestidx;; i++) {
-		if(j==*cyclelen) break;                 // Stop filling in, if we reached the end
-		if(i==pathlen-1) i = startidx;          // If we reach the end of the path, jump to the start
-		cycle[j] = path[i];                     // Write at position `j' the element from the cycle on the path at position `i`'
-		j++;                                    // Increase index to write on `cycle'
-	}
+void canonical_cycle(int n, int *cycle, int *cyclelen, const int *path, \
+    int pathlen) {
+  int startidx;
+  int pathhascycle = 0;
+  const int last = path[pathlen-1];       // This is the last node
+  for (int k = 0; k < pathlen-1; k++) {
+    if (path[k] == last) {
+    // We have found the first occurance of node `last' on the path
+    // The cycle is the following part of path: path[startidx],
+    // path[startidx+1], ..., path[pathlen-2], and therefore has
+    // length pathlen-startidx-1
+      startidx = k;
+      (*cyclelen) = pathlen-startidx-1;
+      pathhascycle = 1;
+      break;
+    }
+  }
+  if (!pathhascycle) {  // Sanity check. This should NEVER happen!
+    fprintf(stderr, "ERROR: You provided a path to `canonical_cycle' that does NOT contain a cycle as its tail! I'm returning now.\n");
+    return;
+  }
+  // Generate canonical representation
+  // First, find smallest node on cycle
+  int smallestidx = -1;
+  int smallest = n;
+  for (int k = startidx; k < pathlen-1; k++) {
+    if (smallest > path[k]) {
+      smallest = path[k];
+      smallestidx = k;
+    }
+  }
+  // Now, copy cycle to `cycle' starting from the smallest
+  // This is the index on where to write on the array `cycle'
+  int j = 0;
+  for (int i = smallestidx;; i++) {
+    // Stop filling in, if we reached the end
+    if (j == *cyclelen) break;
+    // If we reach the end of the path, jump to the start
+    if (i == pathlen-1) i = startidx;
+    // Write at position `j' the element from the cycle on the path at
+    // position `i`'
+    cycle[j] = path[i];
+    // Increase index to write on `cycle'
+    j++;
+  }
 }
 
 /***
  * Return 1 if the cycle is not in `cycles', otherwise 0.
  ***/
-int cycleisnew(int n, const int *c, int clen, const int *cycles, const int *cyclescnt) {
-	const int cnt = cyclescnt[clen];
-	const int bs = blocksize(n);
-	for(int i=0; i<cnt; i++) {                      // Iterate over all cycles of the same length
-		int j;
-		const int startidx = clen*bs + i*clen;  // We start here
-		for(j=0; j<clen; j++)                   // Go position by position
-			if(cycles[startidx+j] != c[j])
-				break;                  // They differ at position j
-		if(j==clen)                             // We didn't break
-			return 0;
-	}
-	return 1;
+int cycleisnew(int n, const int *c, int clen, const int *cycles, \
+    const int *cyclescnt) {
+  const int cnt = cyclescnt[clen];
+  const int bs = blocksize(n);
+  // Iterate over all cycles of the same length
+  for (int i = 0; i < cnt; i++) {
+    int j;
+    const int startidx = clen*bs + i*clen;  // We start here
+    for (j = 0; j < clen; j++)              // Go position by position
+      if (cycles[startidx+j] != c[j])
+        break;                              // They differ at position j
+    if (j == clen)                          // We didn't break
+      return 0;
+  }
+  return 1;
 }
 
 /***
  * Save the given cycle `c' of length `clen' in `cycles'.
  ***/
 void savecycle(int n, const int *c, int clen, int *cycles, int *cyclescnt) {
-	const int m = cyclescnt[clen];                          // There are already `m' cycles of length `clen'
-	const int bs = blocksize(n);
-	const int startidx = clen*bs + m*clen;                  // We start writing here
-	for(int i=0; i<clen; i++) cycles[startidx+i] = c[i];
-	cyclescnt[clen]++;                                      // Increase the counter that holds the number of cycles of length `clen'
+  const int m = cyclescnt[clen];           // There are `m' `clen'-cycles
+  const int bs = blocksize(n);
+  const int startidx = clen*bs + m*clen;  // We start writing here
+  for (int i = 0; i < clen; i++) cycles[startidx+i] = c[i];
+  cyclescnt[clen]++;                      // Increase the counter that holds
+                                          // the number of `clen'-cycles
 }
 
 
@@ -202,37 +218,48 @@ void savecycle(int n, const int *c, int clen, int *cycles, int *cyclescnt) {
  * Parameter `pathlen':     Lenght of `path'
  * Parameter `visited':     n-array to mark visited nodes
  * Parameter `n':           Number of nodes
- * Parameter `found':	Pointer to integer to count the number of cycles found
+ * Parameter `found':  Pointer to integer to count the number of cycles found
  ***/
-void dfs(const int *useless, int *path, int pathlen, int *visited, int n, const int *children, const int *childrenlen, int *cycles, int *cyclescnt, unsigned int *found) {
-	int cycle[n];                                           // Temporary store for cycle
-	int cyclelen;                                           // Length of temporarilly stored cycle
-	const int vertex = path[pathlen-1];                     // We go on duing a dfs from node `vertex'
-	for(int i=0; i<childrenlen[vertex]; i++) {
-		const int u = children[n*vertex + i];           // Enter node `u'
-		if(useless[u]) continue;                        // This one is useless
-		if(visited[u] == 0) {                           // Not visited yet
-			visited[u] = 1;                         // Mark as visited
-			path[pathlen] = u;                      // Append to path, and re-enter recursion...
-			dfs(useless, path, pathlen+1, visited, n, children, childrenlen, cycles, cyclescnt, found);
-			visited[u] = 0;                         // Undo marking, we are going to take another route
-		} else {
-			// We have reached a node on the path, i.e., we have found a cycle!
-			// The cycle is somewhere along the path
-			// path[0], path[1], ... path[pathlen-1]
-			// More specifically, node `u' is somewhere there.
-			// To detect this, we first place `u' on the path.
-			// The method `canonical_cycle' will then use this information to find and return the canonical form of the cycle.
-			path[pathlen] = u;
-			canonical_cycle(n, cycle, &cyclelen, path, pathlen+1);          // Find cycle and generate the canonical representation.
-			if(cycleisnew(n, cycle, cyclelen, cycles, cyclescnt)) {         // Is it new? If yes, save
-				(*found)++;                                             // Increase the number of found cycles by 1
-				savecycle(n, cycle, cyclelen, cycles, cyclescnt);       // Save cycle
-			}
-			// else dont do anything..
-		}
-	}
-	return;
+void dfs(const int *useless, int *path, int pathlen, int *visited, int n, \
+    const int *children, const int *childrenlen, int *cycles, int *cyclescnt, \
+    unsigned int *found) {
+  int cycle[n];  // Temporary store for cycle
+  int cyclelen;  // Length of temporarilly stored cycle
+  // We go on duing a dfs from node `vertex'
+  const int vertex = path[pathlen-1];
+  for (int i = 0; i < childrenlen[vertex]; i++) {
+    const int u = children[n*vertex + i];  // Enter node `u'
+    if (useless[u]) continue;              // This one is useless
+    if (visited[u] == 0) {                 // Not visited yet
+      visited[u] = 1;                      // Mark as visited
+      path[pathlen] = u;                   // Append to path, and re-enter
+                                           // recursion...
+      dfs(useless, path, pathlen+1, visited, n, children, childrenlen, cycles, \
+          cyclescnt, found);
+      visited[u] = 0;                      // Undo marking, we are going to
+                                           // take another route
+    } else {
+      // We have reached a node on the path, i.e., we have found a cycle!
+      // The cycle is somewhere along the path
+      // path[0], path[1], ... path[pathlen-1]
+      // More specifically, node `u' is somewhere there.
+      // To detect this, we first place `u' on the path.
+      // The method `canonical_cycle' will then use this information to find
+      // and return the canonical form of the cycle.
+      path[pathlen] = u;
+      // Find cycle and generate the canonical representation.
+      canonical_cycle(n, cycle, &cyclelen, path, pathlen+1);
+      // Is it new? If yes, save
+      if (cycleisnew(n, cycle, cyclelen, cycles, cyclescnt)) {
+        // Increase the number of found cycles by 1
+        (*found)++;
+        // Save cycle
+        savecycle(n, cycle, cyclelen, cycles, cyclescnt);
+      }
+      // else dont do anything...
+    }
+  }
+  return;
 }
 
 /***
@@ -240,39 +267,46 @@ void dfs(const int *useless, int *path, int pathlen, int *visited, int n, const
  * The cycles are stored in `cycles' and `cyclescnt'.
  * This uses a recursive depth-first-search.
  ***/
-int find_cycles(int *cycles, int *cyclescnt, int n, const int *children, const int *childrenlen) {
-	int visited[n];
-	int useless[n];                 // Some nodes are useless, mark them as such
-	int path[n+1];                  // We use n+1 here because we will "close" the path to form a cycle, i.e., one node will be repeated.
-	unsigned int found = 0;         // Store total number of cycles found
-	// Initialize cyclescnt to zero
-	for(int i=0; i<n+1; i++)
-		cyclescnt[i] = 0;
-	// Mark nodes without children as useless
-	for(int i=0; i<n; i++) {
-		useless[i] = (childrenlen[i]==0) ? 1 : 0;
-	}
-	// Iterate over nodes
-	for(int vertex=0; vertex<n; vertex++) {
-		if(useless[vertex]) continue;
-		for(int i=0; i<n; i++) visited[i] = 0;  // Mark all as univisited
-		path[0] = vertex;                       // We start at `vertex'
-		visited[vertex] = 1;                    // Mark `vertex' as visited. If we revisit a node, we found a cycle.
-		// Enter depth-first-search recursion
-		dfs(useless, path, 1, visited, n, children, childrenlen, cycles, cyclescnt, &found);
-		// We processed this. Every cycle that goes through `vertex' has been found.
-		useless[vertex] = 1;
-	}
-	return found;
+int find_cycles(int *cycles, int *cyclescnt, int n, const int *children, \
+    const int *childrenlen) {
+  int visited[n];
+  // Some nodes are useless, mark them as such
+  int useless[n];
+  // We use n+1 here because we will "close" the path to form a cycle,
+  // i.e., one node will be repeated.
+  int path[n+1];
+  unsigned int found = 0;         // Store total number of cycles found
+  // Initialize cyclescnt to zero
+  for (int i = 0; i < n+1; i++)
+    cyclescnt[i] = 0;
+  // Mark nodes without children as useless
+  for (int i = 0; i < n; i++) {
+    useless[i] = (childrenlen[i] == 0) ? 1 : 0;
+  }
+  // Iterate over nodes
+  for (int vertex = 0; vertex < n; vertex++) {
+    if (useless[vertex]) continue;
+    for (int i = 0; i < n; i++) visited[i] = 0;  // Mark all as univisited
+    // We start at `vertex'
+    path[0] = vertex;
+    // Mark `vertex' as visited. If we revisit a node, we found a cycle.
+    visited[vertex] = 1;
+    // Enter depth-first-search recursion
+    dfs(useless, path, 1, visited, n, children, childrenlen, cycles, \
+        cyclescnt, &found);
+    // We processed this. Every cycle that goes through `vertex' has been found.
+    useless[vertex] = 1;
+  }
+  return found;
 }
 
 /***
  * gissoc tests wheter the graph provided is a SOC or not.
- * Parameter `n':		Number of nodes
- * Parameter `parents':	Pointer to list of parents per node
- * Parameter `parentslen':	Pointer to list of number of parents per node
- * Parameter `children':	Pointer to list of children per node
- * Parameter `childrenlen':	Pointer to list of number of children per node
+ * Parameter `n':    Number of nodes
+ * Parameter `parents':  Pointer to list of parents per node
+ * Parameter `parentslen':  Pointer to list of number of parents per node
+ * Parameter `children':  Pointer to list of children per node
+ * Parameter `childrenlen':  Pointer to list of number of children per node
  * Parameter `cycles':          Pointer to list of cycles in the graph
  * Parameter `cyclescnt':       Pointer to list of cycles count
  *
@@ -280,100 +314,120 @@ int find_cycles(int *cycles, int *cyclescnt, int n, const int *children, const i
  *
  * A SOC (Siblings-on-Cycles) is a simple digraph where every cycle has siblings, i.e., nodes with common parents.
  ***/
-int gissoc(int n, const int *parents, const int *parentslen, const int *children, const int *childrenlen, int *cycles, int *cyclescnt) {
-	// We store the cycles in groups of the cycle lenghts (length of cycle is the number of edges, or equivalentelly, the number of distinct nodes)
-	// Since the graph is a simple digraph, the smallest cycles have length 2.
-	// cycles[0], cycles[1] is the first cycle of length 2
-	// cycles[2], cycles[3] is the second cycle of length 2
-	// etc.
-	// There are at most (n choose 2) cycles of length 2.
-	// There are at most (n choose n/2) cycles of length n/2. This is upper bounded by 2^n/sqrt(pi*n/2) < 2^n; this is the largest binomial coefficient.
-	// Since generally n will be small, we just use this bound for all groups.
-	// So, the cycles with length k are saved at and after cycles[k*(2^n)].
-	// Each block has size 2^n.
-	/* Early and saftey check for self loops. SOCs dont have self-loops, and we should never generate them. Anyhow, test */
-	for(int i=0; i<n; i++) {
-		int num = childrenlen[i];
-		for(int cidx=0; cidx<num; cidx++)
-			if(i == children[n*i+cidx]) {
-				fprintf(stderr, "ERROR: THIS GRAPH HAS A SELF LOOP (node %d)\n", i);
-				return 0;
-			}
-	}
-	// Check SOC property
-	const int bs = blocksize(n);
-	// Iterate over cycle length; min length is 2, max length is n
-	// then iterate over cycle of length `len'
-	// then iterate over nodes of that cycle.
-	// In that iteration process, mark all parents.
-	// If a marked parent is marked again, we have found a common parent.
-	for(int len=2; len<n+1; len++) {
-		const int cnt = cyclescnt[len];             // Number of cycles of length `len'
-		for(int i=0; i<cnt; i++) {                  // Iterate over cycles
-			const int startidx = len*bs+len*i;  // The cycle we look at (length=`len', index=`i') is saved starting from index `startidx'
-			// Initialize marks
-			int marked[n];
-			for(int k=0; k<n; k++) marked[k] = 0;
-			// Flag: Have we found a common parent?
-			int commonparentfound = 0;
-			// Iterate over cycle
-			for(int k=0; k<len; k++) {
-				const int vertex = cycles[startidx+k];
-				for(int pidx=0; pidx<parentslen[vertex]; pidx++) {
-					const int parent = parents[n*vertex+pidx];
-					if(marked[parent]) commonparentfound = 1;    // The node `parent' is a common parent
-					else marked[parent] = 1;                     // Mark the parent node
-					if(commonparentfound) break;                 // We have found a common parent, break out of the for loop over `pidx'
-				}
-				if(commonparentfound) break;                         // We have found a common parent, break out of the for loop over `k'
-			}
-			if(!commonparentfound)                          // On cycle `i' we did NOT find a common parent, this is not a SOC
-				return 0;
-		}
-	}
-	return 1;                                                       // If all test pass, we have found a SOC
+int gissoc(int n, const int *parents, const int *parentslen, \
+    const int *children, const int *childrenlen, int *cycles, int *cyclescnt) {
+  // We store the cycles in groups of the cycle lenghts (length of cycle is the
+  // number of edges, or equivalentelly, the number of distinct nodes)
+  // Since the graph is a simple digraph, the smallest cycles have length 2.
+  // cycles[0], cycles[1] is the first cycle of length 2
+  // cycles[2], cycles[3] is the second cycle of length 2
+  // etc.
+  // There are at most (n choose 2) cycles of length 2.
+  // There are at most (n choose n/2) cycles of length n/2. This is upper
+  // bounded by 2^n/sqrt(pi*n/2) < 2^n; this is the largest binomial
+  // coefficient.
+  // Since generally n will be small, we just use this bound for all groups.
+  // So, the cycles with length k are saved at and after cycles[k*(2^n)].
+  // Each block has size 2^n.
+  /* Early and saftey check for self loops. SOCs dont have self-loops, and we
+   * should never generate them. Anyhow, test */
+  for (int i = 0; i < n; i++) {
+    int num = childrenlen[i];
+    for (int cidx = 0; cidx < num; cidx++)
+      if (i == children[n*i+cidx]) {
+        fprintf(stderr, "ERROR: THIS GRAPH HAS A SELF LOOP (node %d)\n", i);
+        return 0;
+      }
+  }
+  // Check SOC property
+  const int bs = blocksize(n);
+  // Iterate over cycle length; min length is 2, max length is n
+  // then iterate over cycle of length `len'
+  // then iterate over nodes of that cycle.
+  // In that iteration process, mark all parents.
+  // If a marked parent is marked again, we have found a common parent.
+  for (int len = 2; len < n+1; len++) {
+    // Number of cycles of length `len'
+    const int cnt = cyclescnt[len];
+    // Iterate over cycles
+    for (int i = 0; i < cnt; i++) {
+      // The cycle we look at (length=`len', index=`i') is saved starting from
+      // index `startidx'
+      const int startidx = len*bs+len*i;
+      // Initialize marks
+      int marked[n];
+      for (int k = 0; k < n; k++) marked[k] = 0;
+      // Flag: Have we found a common parent?
+      int commonparentfound = 0;
+      // Iterate over cycle
+      for (int k = 0; k < len; k++) {
+        const int vertex = cycles[startidx+k];
+        for (int pidx = 0; pidx < parentslen[vertex]; pidx++) {
+          const int parent = parents[n*vertex+pidx];
+          // Check if `parent' is a common parent, else mark
+          if (marked[parent])
+            commonparentfound = 1;
+          else
+            marked[parent] = 1;
+          // We have found a common parent, break out of the `pidx' loop
+          if (commonparentfound) break;
+        }
+        // We have found a common parent, break out of the for loop over `k'
+        if (commonparentfound) break;
+      }
+      // On cycle `i' we did NOT find a common parent, this is not a SOC
+      if (!commonparentfound)
+        return 0;
+    }
+  }
+  // If all test pass, we have found a SOC
+  return 1;
 }
 
 /***
  * Number of grandparents, and number of grandchildren
  ***/
 int grandparentslen(int node, const int *parents, const int *parentslen) {
-	int tot = 0;
-	for(int i=0; i<parentslen[node]; i++)
-		tot += parentslen[i];
-	return tot;
+  int tot = 0;
+  for (int i = 0; i < parentslen[node]; i++)
+    tot += parentslen[i];
+  return tot;
 }
 int grandchildrenlen(int node, const int *children, const int *childrenlen) {
-	int tot = 0;
-	for(int i=0; i<childrenlen[node]; i++)
-		tot += childrenlen[i];
-	return tot;
+  int tot = 0;
+  for (int i = 0; i < childrenlen[node]; i++)
+    tot += childrenlen[i];
+  return tot;
 }
 /***
  * Speed-up by ignoring all graphs that do not satisfy some ordering of out- and in-degrees.
  * If the graph does not satisfy this `canonical' form, then another ismorphic to it (it's simply a permutation of the nodes)
  * will be degree-ordered.
  ***/
-int isdegreeordered(int n, const int *parents, const int *parentslen, const int *children, const int *childrenlen) {
-	// Decreasingling order by in-degree, then out-degree, then ...
-	for(int i=0; i<n-1; i++) {
-		if(parentslen[i] > parentslen[i+1]) {
-			return 0;
-		} else if(parentslen[i] == parentslen[i+1]) {
-			if(childrenlen[i] > childrenlen[i+1]) {
-				return 0;
-			} else if(childrenlen[i] == childrenlen[i+1]) {
-				if(grandparentslen(i,parents,parentslen) > grandparentslen(i+1,parents,parentslen)) {
-					return 0;
-				} else if(grandparentslen(i,parents,parentslen) == grandparentslen(i+1,parents,parentslen)) {
-					if(grandchildrenlen(i,children,childrenlen) > grandchildrenlen(i+1,children,childrenlen)) {
-						return 0;
-					}
-				}
-			}
-		}
-	}
-	return 1;
+int isdegreeordered(int n, const int *parents, const int *parentslen, \
+    const int *children, const int *childrenlen) {
+  // Decreasingling order by in-degree, then out-degree, then ...
+  for (int i = 0; i < n-1; i++) {
+    if (parentslen[i] > parentslen[i+1]) {
+      return 0;
+    } else if (parentslen[i] == parentslen[i+1]) {
+      if (childrenlen[i] > childrenlen[i+1]) {
+        return 0;
+      } else if (childrenlen[i] == childrenlen[i+1]) {
+        if (grandparentslen(i, parents, parentslen) > \
+            grandparentslen(i+1, parents, parentslen)) {
+          return 0;
+        } else if (grandparentslen(i, parents, parentslen) == \
+            grandparentslen(i+1, parents, parentslen)) {
+          if (grandchildrenlen(i, children, childrenlen) > \
+              grandchildrenlen(i+1, children, childrenlen)) {
+            return 0;
+          }
+        }
+      }
+    }
+  }
+  return 1;
 }
 
 /* Random nubers */
@@ -390,155 +444,175 @@ int isdegreeordered(int n, const int *parents, const int *parentslen, const int
 #endif
 
 u_int64_t rand_uint64(unsigned long long max) {
-	u_int64_t r = 0;
-	for (int i=LOOP_COUNT; i > 0; i--) {
-		r = r*(RAND_MAX + (u_int64_t)1) + rand();
-	}
-	return r%max;
+  u_int64_t r = 0;
+  for (int i=LOOP_COUNT; i > 0; i--) {
+    r = r*(RAND_MAX + (u_int64_t)1) + rand();
+  }
+  return r%max;
 }
 /* End of random numbers */
 
 int main(int argc, char *argv[]) {
-	// Parse command-line arguments
-	int n = -1;
-	int NONDAGONLY = 0;
-	int NOSINK = 0;
-	int NOSOURCE = 0;
-	int RANDOM = 0;
-	int GRAPHVIZ = 0;
-	int UNKOPTION = 0;
-	for(int i=1; i<argc; i++) {
-		if(strcmp(argv[i], "-n") == 0 && i+1 < argc) {
-			n = atoi(argv[i+1]);
-			i++;
-		} else if(strcmp(argv[i], "-r") == 0 && i+1 < argc) {
-			RANDOM = atoi(argv[i+1]);
-			i++;
-		} else if(strcmp(argv[i], "-c") == 0)
-			NONDAGONLY = 1;
-		else if(strcmp(argv[i], "--no-sink") == 0)
-			NOSINK = 1;
-		else if(strcmp(argv[i], "--no-source") == 0)
-			NOSOURCE = 1;
-		else if(strcmp(argv[i], "--graphviz") == 0)
-			GRAPHVIZ = 1;
-		else {
-			UNKOPTION = 1;
-			break;
-		}
-	}
-	if(UNKOPTION || n<=1) {
-		fprintf(stderr, "Usage: %s -n <order> [-r <num>] [--grpahviz] [FILTER ...]\n", argv[0]);
-		fprintf(stderr, "       -n <order>     Generate SOCs with `order' connected nodes\n");
-		fprintf(stderr, "       -r <num>       Pick directed graphs at random, and exit after having found `num' SOCs\n");
-		fprintf(stderr, "       --graphviz     Output SOCs in Graphviz format, arcs of common parents are highlighted\n");
-		fprintf(stderr, "\n");
-		fprintf(stderr, "[FILTER]              Consider only simple directed graphs ...\n");
-		fprintf(stderr, "       -c             ... that are cyclic (i.e., not DAGs)\n");
-		fprintf(stderr, "       --no-sink      ... without sink nodes (this logically implies -c)\n");
-		fprintf(stderr, "       --no-source    ... without source nodes (also this logically implies -c)\n");
-		fprintf(stderr, "\n");
-		fprintf(stderr, "This program prints the found SOCs as adjacency matrices to stdout, unless --grpahviz has been specified.\n");
-		fprintf(stderr, "To exclude (some) of the isomorphic SOCs, it uses a degree-order filter.\n");
-		return -1;
-	}
-	// EO Parse command-line arguments
-	// Setup datastructures
-	int parents[n*n];       // Each node i can have at most n-1 parents. These are listed at parents[n*i+k] for 0<=k<n
-	int parentslen[n];      // Number of parents of node i is stored in parentslen[i]
-	int children[n*n];      // Each node i can have at most n-1 children. These are listed at parents[n*i+k] for 0<=k<n
-	int childrenlen[n];     // Number of children of node i is stored in childrenlen[i]
-	const int bs = blocksize(n);
-	int *cycles = (int*)malloc(sizeof(int)*(n+2)*bs);
-	if(!cycles) {
-		fprintf(stderr, "Failed to allocate memory to store cycles (tried to allocate %lu bytes)\n", sizeof(int)*(n+2)*bs);
-		return -1;
-	}
-	int cyclescnt[n+1];      // cyclescnt[k] stores the number of cycles with length k
-	// EO Setup datastructures
-	// Initiate graph enumeration
-	const int m = n*(n-1);
-	if(m>64) {
-		fprintf(stderr, "Too many graphs to enumarate with an unsiged long long integer (number of node pairs = %d)\n", m);
-		return -1;
-	}
-	unsigned long long max = 1L << m;		// Largest `graphnumber'
-	const int padlen = (int)((float)m/3.322)+1;	// Convert log 2 to log 10
-	int len = 0;
-	time_t t0 = time(NULL);
-	time_t t = time(NULL);
-	srand((unsigned) time(&t));
-	// EO Initiate graph enumeration
-	fprintf(stderr, "Generating SOCs with %d nodes\n", n);
-	if(RANDOM) fprintf(stderr, "  Picking graphs at random\n");
-	if(NONDAGONLY) fprintf(stderr, "  Filter: Omitting DAGs\n");
-	if(NOSINK) fprintf(stderr, "  Filter: Omitting graphs with sink nodes\n");
-	if(NOSOURCE) fprintf(stderr, "  Filter: Omitting graphs with source nodes\n");
-	unsigned long long graphnumber = 0;	// We will always omit the graph with `graphnumber' 0; it is unintersting anyway
-	unsigned long long graphschecked = 0;
-	for(;graphnumber < max && (!RANDOM || len < RANDOM);) {
-		const time_t now = time(NULL);
-		// Break if we exhaustively check all graphs or if we found enough SOCs at random
-		//if(graphnumber >= max || (RANDOM && len >= RANDOM))
-		//	break;
-		// Print status every second
-		if(now > t) {
-			t = now;
-			const int deltat = now - t0;
-			const float rate = (float)graphschecked/(float)deltat;
-			const float SOCrate = (float)len/(float)deltat;
-			const float percentage = (RANDOM==0) ? 100*(float)graphnumber/(float)max : 100*(float)len/(float)RANDOM;
-			const float ETC = (RANDOM==0) ? ((float)(max - graphschecked)/(rate))/3600 : ((float)(RANDOM - len)/(SOCrate))/3600;
-			fprintf(stderr, "\r%6.2f%%    %*llu/%llu (%i SOCs found, %d seconds, %*.2f graphs/s, %*.2f SOCs/s, %*.2fh estimated time left)",percentage,padlen,graphnumber,max,len,deltat,padlen+3,rate,padlen+3,SOCrate,padlen,ETC);
-			fflush(stderr);
-		}
-		// Convert graph index `grpahnumber' to the lists parents, children, parentslen, childrenlen
-		graphnrtolists(n, graphnumber, parents, parentslen, children, childrenlen);
-		// Increase checked counter and prepare for next iteration
-		graphschecked++;
-		if(RANDOM)
-			graphnumber = rand_uint64(max);
-		else
-			graphnumber++;
-		// EO Increase checked counter and prepare for next iteration
-		// If enabled, compare against degree options
-		if(NOSINK || NOSOURCE) {
-			// Get some degree properties
-			int minindegree = n;
-			int minoutdegree = n;
-			for(int node=0; node<n; node++) {
-				if(minindegree > parentslen[node])
-					minindegree = parentslen[node];
-				if(minoutdegree > childrenlen[node])
-					minoutdegree = childrenlen[node];
-			}
-			if(NOSINK && minoutdegree == 0) continue;
-			if(NOSOURCE && minindegree == 0) continue;
-		}
-		// Check whether this graph is canonical or not; continue if graph has not proper ordering of degress (there is an isomorphic graphs to this one that has been or will be checked)
-		if(!isdegreeordered(n, parents, parentslen, children, childrenlen)) continue;
-		// Ignore graphs that are not connected
-		if(!isconnected(n, parents, parentslen, children, childrenlen)) continue;
-		// Find cycles
-		const int num_cycles = find_cycles(cycles, cyclescnt, n, children, childrenlen);
-		// If enabled, ignore DAGs
-		if(NONDAGONLY && num_cycles > 0) continue;
-		// Test the Siblings-On-Cycles property
-		// A DAG is trivially a SOC
-		if(num_cycles > 0 && !gissoc(n, parents, parentslen, children, childrenlen, cycles, cyclescnt)) continue;
-		// We have found a SOC
-		len++;
-		if(GRAPHVIZ)
-			dumpgv(n, graphnumber, children, childrenlen);
-		else
-			dumpgraph(n, children, childrenlen);
-	}
-	const int deltat = time(NULL) - t0;
-	const float rate = (deltat == 0) ? graphschecked : (float)graphschecked/(float)deltat;
-	const float SOCrate = (deltat == 0) ? len : (float)len/(float)deltat;
-	fprintf(stderr, "\r%6.2f%%    %*llu/%llu (%i SOCs found, %d seconds, %*.2f graphs/s, %*.2f SOCs/s, %*.2fh estimated time left)",100.00,padlen,graphnumber,max,len,deltat,padlen+3,rate,padlen+3,SOCrate,padlen,0.00);
-	fprintf(stderr, "\nFound %d SOCs in %d seconds\n", len,deltat);
-	// Free manually allocated memory
-	free(cycles);
-	return 0;
+  // Parse command-line arguments
+  int n = -1;
+  int NONDAGONLY = 0;
+  int NOSINK = 0;
+  int NOSOURCE = 0;
+  int RANDOM = 0;
+  int GRAPHVIZ = 0;
+  int UNKOPTION = 0;
+  for (int i = 1; i < argc; i++) {
+    if (strcmp(argv[i], "-n") == 0 && i+1 < argc) {
+      n = atoi(argv[i+1]);
+      i++;
+    } else if (strcmp(argv[i], "-r") == 0 && i+1 < argc) {
+      RANDOM = atoi(argv[i+1]);
+      i++;
+    } else if (strcmp(argv[i], "-c") == 0) {
+      NONDAGONLY = 1;
+    } else if (strcmp(argv[i], "--no-sink") == 0) {
+      NOSINK = 1;
+    } else if (strcmp(argv[i], "--no-source") == 0) {
+      NOSOURCE = 1;
+    } else if (strcmp(argv[i], "--graphviz") == 0) {
+      GRAPHVIZ = 1;
+    } else {
+      UNKOPTION = 1;
+      break;
+    }
+  }
+  if (UNKOPTION || n <= 1) {
+    fprintf(stderr, "Usage: %s -n <order> [-r <num>] [--grpahviz] [FILTER ...]\n", argv[0]);
+    fprintf(stderr, "       -n <order>     Generate SOCs with `order' connected nodes\n");
+    fprintf(stderr, "       -r <num>       Pick directed graphs at random, and exit after having found `num' SOCs\n");
+    fprintf(stderr, "       --graphviz     Output SOCs in Graphviz format, arcs of common parents are highlighted\n");
+    fprintf(stderr, "\n");
+    fprintf(stderr, "[FILTER]              Consider only simple directed graphs ...\n");
+    fprintf(stderr, "       -c             ... that are cyclic (i.e., not DAGs)\n");
+    fprintf(stderr, "       --no-sink      ... without sink nodes (this logically implies -c)\n");
+    fprintf(stderr, "       --no-source    ... without source nodes (also this logically implies -c)\n");
+    fprintf(stderr, "\n");
+    fprintf(stderr, "This program prints the found SOCs as adjacency matrices to stdout, unless --grpahviz has been specified.\n");
+    fprintf(stderr, "To exclude (some) of the isomorphic SOCs, it uses a degree-order filter.\n");
+    return -1;
+  }
+  // EO Parse command-line arguments
+  // Setup datastructures
+  int parents[n*n];    // Each node i can have at most n-1 parents. These are listed at parents[n*i+k] for 0<=k<n
+  int parentslen[n];   // Number of parents of node i is stored in parentslen[i]
+  int children[n*n];   // Each node i can have at most n-1 children. These are listed at parents[n*i+k] for 0<=k<n
+  int childrenlen[n];  // Number of children of node i is stored in childrenlen[i]
+  const int bs = blocksize(n);
+  int *cycles = (int*)malloc(sizeof(int)*(n+2)*bs);
+  if (!cycles) {
+    fprintf(stderr, "Failed to allocate memory to store cycles (tried to allocate %lu bytes)\n", \
+        sizeof(int)*(n+2)*bs);
+    return -1;
+  }
+  int cyclescnt[n+1];  // cyclescnt[k] stores the number of cycles with length k
+  // EO Setup datastructures
+  // Initiate graph enumeration
+  const int m = n*(n-1);
+  if (m > 64) {
+    fprintf(stderr, "Too many graphs to enumarate with an unsiged long long integer (number of node pairs = %d)\n", \
+        m);
+    return -1;
+  }
+  unsigned long long max = 1L << m;            // Largest `graphnumber'
+  const int padlen = (int)((float)m/3.322)+1;  // Convert log 2 to log 10
+  int len = 0;
+  time_t t0 = time(NULL);
+  time_t t = time(NULL);
+  srand((unsigned) time(&t));
+  // EO Initiate graph enumeration
+  fprintf(stderr, "Generating SOCs with %d nodes\n", n);
+  if (RANDOM)
+    fprintf(stderr, "  Picking graphs at random\n");
+  if (NONDAGONLY)
+    fprintf(stderr, "  Filter: Omitting DAGs\n");
+  if (NOSINK)
+    fprintf(stderr, "  Filter: Omitting graphs with sink nodes\n");
+  if (NOSOURCE)
+    fprintf(stderr, "  Filter: Omitting graphs with source nodes\n");
+  // We will always omit the graph with `graphnumber' 0;
+  // it is unintersting anyway
+  unsigned long long graphnumber = 0;
+  unsigned long long graphschecked = 0;
+  for (; graphnumber < max && (!RANDOM || len < RANDOM);) {
+    const time_t now = time(NULL);
+    // Print status every second
+    if (now > t) {
+      t = now;
+      const int deltat = now - t0;
+      const float rate = (float)graphschecked/(float)deltat;
+      const float SOCrate = (float)len/(float)deltat;
+      const float percentage = (RANDOM == 0) ? \
+          100*(float)graphnumber/(float)max : 100*(float)len/(float)RANDOM;
+      const float ETC = (RANDOM == 0) ? \
+          ((float)(max - graphschecked)/(rate))/3600 : \
+          ((float)(RANDOM - len)/(SOCrate))/3600;
+      fprintf(stderr, "\r%6.2f%%    %*llu/%llu (%i SOCs found, %d seconds, %*.2f graphs/s, %*.2f SOCs/s, %*.2fh estimated time left)", \
+          percentage, padlen, graphnumber, max, len, deltat, padlen+3, rate, \
+          padlen+3, SOCrate, padlen, ETC);
+      fflush(stderr);
+    }
+    // Convert graph index `grpahnumber' to the lists parents, children,
+    // parentslen, childrenlen
+    graphnrtolists(n, graphnumber, parents, parentslen, children, childrenlen);
+    // Increase checked counter and prepare for next iteration
+    graphschecked++;
+    if (RANDOM)
+      graphnumber = rand_uint64(max);
+    else
+      graphnumber++;
+    // EO Increase checked counter and prepare for next iteration
+    // If enabled, compare against degree options
+    if (NOSINK || NOSOURCE) {
+      // Get some degree properties
+      int minindegree = n;
+      int minoutdegree = n;
+      for (int node = 0; node < n; node++) {
+        if (minindegree > parentslen[node])
+          minindegree = parentslen[node];
+        if (minoutdegree > childrenlen[node])
+          minoutdegree = childrenlen[node];
+      }
+      if (NOSINK && minoutdegree == 0) continue;
+      if (NOSOURCE && minindegree == 0) continue;
+    }
+    // Check whether this graph is canonical or not; continue if graph has not
+    // proper ordering of degress (there is an isomorphic graphs to this one
+    // that has been or will be checked)
+    if (!isdegreeordered(n, parents, parentslen, children, childrenlen))
+      continue;
+    // Ignore graphs that are not connected
+    if (!isconnected(n, parents, parentslen, children, childrenlen))
+      continue;
+    // Find cycles
+    const int num_cycles = find_cycles(cycles, cyclescnt, n, children, \
+        childrenlen);
+    // If enabled, ignore DAGs
+    if (NONDAGONLY && num_cycles > 0) continue;
+    // Test the Siblings-On-Cycles property
+    // A DAG is trivially a SOC
+    if (num_cycles > 0 && !gissoc(n, parents, parentslen, children, \
+        childrenlen, cycles, cyclescnt)) continue;
+    // We have found a SOC
+    len++;
+    if (GRAPHVIZ)
+      dumpgv(n, graphnumber, children, childrenlen);
+    else
+      dumpgraph(n, children, childrenlen);
+  }
+  const int deltat = time(NULL) - t0;
+  const float rate = (deltat == 0) ? graphschecked : \
+      (float)graphschecked/(float)deltat;
+  const float SOCrate = (deltat == 0) ? len : (float)len/(float)deltat;
+  fprintf(stderr, "\r%6.2f%%    %*llu/%llu (%i SOCs found, %d seconds, %*.2f graphs/s, %*.2f SOCs/s, %*.2fh estimated time left)", \
+      100.00, padlen, graphnumber, max, len, deltat, padlen+3, rate, padlen+3, \
+      SOCrate, padlen, 0.00);
+  fprintf(stderr, "\nFound %d SOCs in %d seconds\n", len, deltat);
+  // Free manually allocated memory
+  free(cycles);
+  return 0;
 }
diff --git a/src/SOCgraphviz.c b/src/SOCgraphviz.c
index a34d808..94db46b 100644
--- a/src/SOCgraphviz.c
+++ b/src/SOCgraphviz.c
@@ -1,80 +1,80 @@
-// SPDX-FileCopyrightText: 2023 Ämin Baumeler <amin@indyfac.ch> and Eleftherios-Ermis Tselentis <eleftheriosermis.tselentis@oeaw.ac.at>
+// SPDX-FileCopyrightText: 2023 Ämin Baumeler <amin@indyfac.ch> and
+// Eleftherios-Ermis Tselentis <eleftheriosermis.tselentis@oeaw.ac.at>
 //
 // SPDX-License-Identifier: GPL-3.0-or-later
 
 #include <stdio.h>
 
 int main(int argc, char *argv[]) {
-	// Parse command-line arguments
-	if(argc != 2) {
-		fprintf(stderr, "Usage: %s <filename>\n", argv[0]);
-		fprintf(stderr, "       <filename>     File name with adjacency matrices of simple directed graphs\n");
-		fprintf(stderr, "\n");
-		fprintf(stderr, "[FILE FORMAT]\n");
-		fprintf(stderr, "       Each line in `filename' must contain the adjacency matrix of a simple directed graph in the format\n");
-		fprintf(stderr, "         {{a00,a01,...},{a10,a11,...},...}  where aij=1 if and only if the graph has the edge i -> j\n");
-		fprintf(stderr, "       The file `filename' may contain graphs with different order (number of vertices)\n");
-		fprintf(stderr, "\n");
-		fprintf(stderr, "This program translates to adjacency matrices into the Graphviz format, and prints them to stdout.\n");
-		return -1;
-	}
+  // Parse command-line arguments
+  if (argc != 2) {
+    fprintf(stderr, "Usage: %s <filename>\n", argv[0]);
+    fprintf(stderr, "       <filename>     File name with adjacency matrices of simple directed graphs\n");
+    fprintf(stderr, "\n");
+    fprintf(stderr, "[FILE FORMAT]\n");
+    fprintf(stderr, "       Each line in `filename' must contain the adjacency matrix of a simple directed graph in the format\n");
+    fprintf(stderr, "         {{a00,a01,...},{a10,a11,...},...}  where aij=1 if and only if the graph has the edge i -> j\n");
+    fprintf(stderr, "       The file `filename' may contain graphs with different order (number of vertices)\n");
+    fprintf(stderr, "\n");
+    fprintf(stderr, "This program translates to adjacency matrices into the Graphviz format, and prints them to stdout.\n");
+    return -1;
+  }
 
-	// Open file
-	FILE *fp;
-	fp = fopen(argv[1], "r");
-	if(fp == NULL){
-		fprintf(stderr, "ERROR: Could not read graphs from file `%s'.\n", argv[1]);
-		return -1;
-	}
-	// Parse
-	int headprinted = 1;
-	int line = 1;
-	int n = -1;
-	int row = 0;
-	int col = 0;
-	char ch;
-	printf("strict digraph \"File_%s_line_%d\" {", argv[1], line);
-	while(!feof(fp)) {
-		ch = fgetc(fp);
-		switch(ch) {
-			case '{':
-				if(row==0 && col==0 && headprinted==0) {
-					headprinted = 1;
-					printf("strict digraph \"File_%s_line_%d\" {", argv[1], line);
-				}
-				break;
-			case '0':
-			case '1':
-				if(ch=='1')
-					printf("G%dN%d->G%dN%d;", line,row,line,col);
-				col++;
-				break;
-			case '}':
-				if(n==-1)
-					n = col;
-				if(col==n) {
-					row++;
-					col=0;
-				}
-				break;
-			case '\n':
-				printf("}\n");
-				headprinted = 0;
-				line++;
-				row = 0;
-				col = 0;
-				n = -1;
-				break;
-			case ',':
-			case ' ':
-			case 0xffffffff:
-				break;
-			default:
-				fprintf(stderr, "File contains wrongly formatted graph\n");
-				fclose(fp);
-				return -1;	// Format error
-		}
-	}
-	fclose(fp);
-	return 0;
+  // Open file
+  FILE *fp;
+  fp = fopen(argv[1], "r");
+  if (fp == NULL) {
+    fprintf(stderr, "ERROR: Could not read graphs from file `%s'.\n", argv[1]);
+    return -1;
+  }
+  // Parse
+  int headprinted = 1;
+  int line = 1;
+  int n = -1;
+  int row = 0;
+  int col = 0;
+  printf("strict digraph \"File_%s_line_%d\" {", argv[1], line);
+  while (!feof(fp)) {
+    char ch = fgetc(fp);
+    switch (ch) {
+      case '{':
+        if (row == 0 && col == 0 && headprinted == 0) {
+          headprinted = 1;
+          printf("strict digraph \"File_%s_line_%d\" {", argv[1], line);
+        }
+        break;
+      case '0':
+      case '1':
+        if (ch == '1')
+          printf("G%dN%d->G%dN%d;", line, row, line, col);
+        col++;
+        break;
+      case '}':
+        if (n == -1)
+          n = col;
+        if (col == n) {
+          row++;
+          col = 0;
+        }
+        break;
+      case '\n':
+        printf("}\n");
+        headprinted = 0;
+        line++;
+        row = 0;
+        col = 0;
+        n = -1;
+        break;
+      case ',':
+      case ' ':
+      case 0xffffffff:
+        break;
+      default:
+        fprintf(stderr, "File contains wrongly formatted graph\n");
+        fclose(fp);
+        return -1;  // Format error
+    }
+  }
+  fclose(fp);
+  return 0;
 }
diff --git a/tools/removeiso.py b/tools/removeiso.py
index f9337dd..d461df4 100755
--- a/tools/removeiso.py
+++ b/tools/removeiso.py
@@ -1,6 +1,7 @@
 #!/usr/bin/python3
 """
-SPDX-FileCopyrightText: 2023 Ämin Baumeler <amin@indyfac.ch> and Eleftherios-Ermis Tselentis <eleftheriosermis.tselentis@oeaw.ac.at>
+SPDX-FileCopyrightText: 2023 Ämin Baumeler <amin@indyfac.ch> and
+Eleftherios-Ermis Tselentis <eleftheriosermis.tselentis@oeaw.ac.at>
 
 SPDX-License-Identifier: GPL-3.0-or-later
 """