diff --git a/SOCadmissible.c b/SOCadmissible.c
deleted file mode 100644
index 424f94a..0000000
--- a/SOCadmissible.c
+++ /dev/null
@@ -1,404 +0,0 @@
-// 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>
-#include <stdlib.h>
-#include <time.h>
-
-/***
- * Read graphs from  file.
- * 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,
- * and returns the graph data. Graph data is structured as follows:
- * {dimension graph 1 | (0,0) | (0,1) | ... | (1,0) | (1,1) | ... | dimension graph 2 | ...}
- * 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
-
-	// 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;
-}
-
-/***
- * Increase intervention counter by one.
- * 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;
-}
-
-/***
- * Early check
- * With this, we might find some fixed-points entries without the need of entering any recursion.
- * If `party' is not in the range of a parent's intervention, then `party' will definitely receive a 0.
- * If `party' has no parents, then `party' will definitely receive a 1.
- * If the fixed-point value cannot be inferred, return -1, otherwise return the fixed-point entry.
- ***/
-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;
-}
-
-/***
- * Recursively compute the fixed point
- * This methods reuses already computed fixed points.
- * `fp'    pointer to fixed-point array
- * `fplen' the first `fplen' parties have already computed the fixed points
- *
- * 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;
-}
-
-/***
- * ceil(log(x, base=10))
- * Helper function for pretty printing.
- ***/
-int ceillog10(int x) {
-	int res = 0;
-	while(x > 0) {
-		x /= 10;
-		res++;
-	}
-	return res;
-}
-
-/***
- * Verify the admissibility of the graphs
- * Go through all graphs and test them
- * `startidx': Index from which to start, assumed to be non-negative
- * `len':      Number of graphs to test, assumed to be non-negative and such that startidx+len does not exceed the length of the `graphs' array
- * `maxn':     Max. number of nodes
- * `graphs':   Pointer to graph structure
- *
- * Return 0 if test succeeds, else return graphidx+1 of failed graph (this is the line number).
- * 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;
-			}
-
-			// 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;
-
-			// 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;
-	}
-
-	// 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;
-}
diff --git a/SOCgen.c b/SOCgen.c
deleted file mode 100644
index 6e53bcc..0000000
--- a/SOCgen.c
+++ /dev/null
@@ -1,519 +0,0 @@
-// 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>
-#include <stdlib.h>
-#include <string.h>
-#include <time.h>
-
-/***
- * 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");
-}
-
-/***
- * 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;
-}
-
-/***
- * 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);
-	}
-}
-
-/***
- * 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;
-}
-
-/***
- * 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++;
-		}
-
-	}
-}
-
-/***
- * 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 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'
-	}
-}
-
-/***
- * 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;
-}
-
-/***
- * 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'
-}
-
-
-/***
- * Recursive depth-first-search to find cycles in the graph specified by `children' and `childrenlen'.
- * The cycles are stored in `cycles' and `cyclescnt'.
- * Parameter `useless':     n-array where useless nodes are marked
- * Parameter `path':        Array of nodes that represent a path on the graph
- * 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
- ***/
-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;
-}
-
-/***
- * Find cycles in the `n'-nodes graph specified by the arrays `children' and `childrenlen', and return the number of cycles found.
- * 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;
-}
-
-/***
- * 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 `cycles':          Pointer to list of cycles in the graph
- * Parameter `cyclescnt':       Pointer to list of cycles count
- *
- * This function retruns 0 if the graph is NOT a SOC, and 1 otherwise.
- *
- * 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
-}
-
-/***
- * 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 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;
-}
-/***
- * 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;
-}
-
-/* Random nubers */
-#if RAND_MAX/256 >= 0xFFFFFFFFFFFFFF
-#define LOOP_COUNT 1
-#elif RAND_MAX/256 >= 0xFFFFFF
-#define LOOP_COUNT 2
-#elif RAND_MAX/256 >= 0x3FFFF
-#define LOOP_COUNT 3
-#elif RAND_MAX/256 >= 0x1FF
-#define LOOP_COUNT 4
-#else
-#define LOOP_COUNT 5
-#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;
-}
-/* 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 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 {
-			UNKOPTION = 1;
-			break;
-		}
-	}
-	if(UNKOPTION || n<=1) {
-		fprintf(stderr, "Usage: %s -n <order> [-r <num> ] [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, "\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.\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;
-	}
-	const 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++;
-		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;
-}