extened Makefile, added ./tools/removeiso.py, added option --graphviz, readme Makefile

Makefile

@@ -6,5 +6,8 @@ CFLAGS = -Wall -Ofast
 
 all: SOCgen SOCadmissible SOCgraphviz
 
+clean:
+	rm -v SOCgen SOCadmissible SOCgraphviz
+
 %: src/%.c
 	gcc $(CFLAGS) -o $@ $^

README.md

@@ -7,11 +7,12 @@ SPDX-License-Identifier: GPL-3.0-or-later
 # SOC Observation Code
 
 ## Description
-SOC stands for **S**OC **O**bservation **C**ode, and is composed of three C programs:
+SOC stands for **S**OC **O**bservation **C**ode, and is composed of three C programs and a python script:
 
 - One, `SOCgen`, to generate SOC graphs (here, SOC stands for Siblings-on-Cycles),
 - and another, `SOCadmissible`, to verify the admissibility of these graphs as quantum causal structures,
-- and the last, `SOCgraphviz`, to translate adjancency matrices into the [Graphviz language](https://graphviz.org/).
+- and the last, `SOCgraphviz`, to translate adjacency matrices into the [Graphviz language](https://graphviz.org/).
+- The python script `/tools/removeiso.py` removes isomorphic graphs.
 
 The first two programs are used in support of Conjecture 1 in the article [Admissible Causal Structures and Correlations, arXiv:2210.12796 \[quant-ph\]](https://arxiv.org/abs/2210.12796).
 
@@ -31,16 +32,17 @@ To display help and exit, run the respective program without command-line argume
 ### SOCgen
 ```
 $ ./SOCgen 
-Usage: ./SOCgen -n <order> [-r <num> ] [FILTER ...]
+Usage: ./SOCgen -n <order> [-r <num>] [--grpahviz] [FILTER ...]
        -n <order>     Generate SOCs with `order' connected nodes
        -r <num>       Pick directed graphs at random, and exit after having found `num' SOCs
+       --graphviz     Output SOCs in Graphviz format, arcs of common parents are highlighted
 
 [FILTER]              Consider only simple directed graphs ...
        -c             ... that are cyclic (i.e., not DAGs)
        --no-sink      ... without sink nodes (this logically implies -c)
        --no-source    ... without source nodes (also this logically implies -c)
 
-This program prints the found SOCs as adjacency matrices to stdout.
+This program prints the found SOCs as adjacency matrices to stdout, unless --grpahviz has been specified.
 To exclude (some) of the isomorphic SOCs, it uses a degree-order filter.
 ```
 
@@ -75,8 +77,17 @@ Usage: ./SOCgraphviz <filename>
 This program translates to adjacency matrices into the Graphviz format, and prints them to stdout.
 ```
 
+### removeiso.py
+```
+$ ./tools/removeiso.py  
+Usage: ./tools/removeiso.py -f <filename>
+
+Filters out the isomorphic graphs in `filename'.
+```
+This script requires [NumPy](https://numpy.org/) and [NetworkX](https://networkx.org/).
+
 ### Examples
-To generate all SOCs with three nodes, and save them in the file `3.soc`, run:
+To generate all SOCs with three nodes, and save them in the file `3.soc`, you may run:
 ```
 $ ./SOCgen -n 3 > 3.soc
 Generating SOCs with 3 nodes
@@ -93,7 +104,11 @@ These graphs are admissible
 
 The SOCs generated can easily be displayed with the [Graphviz](https://graphviz.org/) tools:
 ```
-./SOCgraphviz 3.soc | dot | gvpack | circo -Nshape=point -Tx11
+$ ./SOCgraphviz 3.soc | dot | gvpack | circo -Nshape=point -Tx11
+```
+also with the `--graphviz` option:
+```
+$ ./SOCgen -n 3 --graphviz | dot | gvpack | neato -Tx11
 ```
 or in [Wolfram Mathematica](https://www.wolfram.com/mathematica/) using:
 ```
@@ -102,10 +117,11 @@ SOCs = DeleteDuplicatesBy[SOCs, CanonicalGraph];
 SOCs
 ```
 
-To generate all SOCs with four nodes, and to display them using Graphviz, you may run:
+To generate all SOCs with four nodes, save them in the file `4.soc`, remove isomorphic ones (as stated below, `SOCgen` only rudimentary checks for isomorphic graphs), and to display them using Graphviz, you may run:
 ```
-./SOCgen -n 4 | ./SOCgraphviz /dev/stdin | dot | gvpack | circo -Nshape=point -Tx11
+$ ./SOCgen -n 4 | tee 4.soc | ./tools/removeiso.py -f 4.soc | ./SOCgraphviz /dev/stdin | dot | gvpack | circo -Nshape=point -Tx11
 ```
+
 ![All SOCs with four nodes](./example.png "SOCs with four nodes")
 
 ## Limitations

src/SOCadmissible.c

@@ -130,7 +130,7 @@ int nextintervention(int n, int *interventionslen, int *interventionidx) {
  * 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.
+ * 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)

src/SOCgen.c

@@ -31,6 +31,24 @@ void dumpgraph(int n, const int *children, const int *childrenlen) {
 	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;
+}
+
 /***
  * Block size of our datastructure.
  ***/
@@ -112,10 +130,10 @@ void graphnrtolists(int n, unsigned long long int graph, int *parents, int *pare
 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
+	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
+		                           	// 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;
@@ -187,9 +205,9 @@ void savecycle(int n, const int *c, int clen, int *cycles, int *cyclescnt) {
  * 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 cycle[n];                                           // Temporary store for cycle
-	int cyclelen;                           // Length of temporarilly stored cycle
+	int cyclelen;                                           // Length of temporarilly stored cycle
-	const int vertex = path[pathlen-1];     // We go on duing a dfs from node `vertex'
+	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
@@ -387,6 +405,7 @@ int main(int argc, char *argv[]) {
 	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) {
@@ -401,22 +420,25 @@ int main(int argc, char *argv[]) {
 			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> ] [FILTER ...]\n", argv[0]);
+		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.\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;
 	}
@@ -440,7 +462,7 @@ int main(int argc, char *argv[]) {
 		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'
+	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);
@@ -506,7 +528,10 @@ int main(int argc, char *argv[]) {
 		if(num_cycles > 0 && !gissoc(n, parents, parentslen, children, childrenlen, cycles, cyclescnt)) continue;
 		// We have found a SOC
 		len++;
-		dumpgraph(n, children, childrenlen);
+		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;

tools/removeiso.py

@@ -0,0 +1,77 @@
+#!/usr/bin/python3
+"""
+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
+"""
+import sys
+import getopt
+import numpy as np
+import networkx as nx
+
+"""Parse command line arguments first."""
+helpstr = """Usage: """+sys.argv[0]+""" -f <filename>
+
+Filters out the isomorphic graphs in `filename'."""
+
+FILENAME = False
+
+try:
+    opts, args = getopt.getopt(sys.argv[1:], "hf:")
+except getopt.GetoptError:
+    print(helpstr, file=sys.stderr)
+    sys.exit(2)
+for opt, arg in opts:
+    if opt == '-h':
+        print(helpstr)
+        sys.exit()
+    elif opt == '-f':
+        FILENAME = arg
+
+if not FILENAME:
+    print(helpstr, file=sys.stderr)
+    sys.exit(2)
+
+
+with open(FILENAME, 'rb') as f:
+    num_lines = sum(1 for _ in f)
+print(f"Total number of digraphs is {num_lines}", file=sys.stderr)
+
+
+nonisodigraphs = set()
+
+cnt = 0
+fnd = 0
+with open(FILENAME, 'r') as f:
+    while 1:
+        cnt = cnt + 1
+        print("\r{:.2f}%".format(100*cnt/num_lines), end='', file=sys.stderr)
+        line = f.readline().strip('{}\n')
+        if not line:
+            break
+        G = nx.from_numpy_matrix(np.vectorize(int)(np.matrix([x.split(',') for x in line.split('},{')])), create_using=nx.DiGraph)
+
+        couldbeold = False
+        for H in nonisodigraphs:
+            # Returns False if graphs are definitely not isomorphic.
+            # True does NOT guarantee isomorphism.
+            if nx.faster_could_be_isomorphic(G, H):
+                couldbeold = True
+                break
+
+        if not couldbeold:
+            nonisodigraphs.add(G)
+            fnd = fnd + 1
+            print(f"{{{{{line}}}}}")
+            continue
+
+        new = True
+        for H in nonisodigraphs:
+            if nx.is_isomorphic(G, H):
+                new = False
+                break
+        if new:
+            nonisodigraphs.add(G)
+            fnd = fnd + 1
+            print(f"{{{{{line}}}}}")
+    print(f"\r100%  Found {fnd} non-isomoric graphs", file=sys.stderr)