amin/SOC-Observation-Code
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

replaced tabs with whitespace, linted c code, reduced line lengths if possible for better readability

Browse Source
This commit is contained in:
Ämin Baumeler 2023-09-07 13:36:33 +02:00
parent ab5a90dc38
commit 642b78c900
4 changed files with 935 additions and 826 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

206
src/SOCadmissible.c
View File

@ -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 // SPDX-License-Identifier: GPL-3.0-or-later
@ -23,16 +24,16 @@ int* read_graphs(char* filename, int *cnt, int *maxdim) {
// Open file // Open file
FILE *fp; FILE *fp;
fp = fopen(filename, "r"); fp = fopen(filename, "r");
if(fp == NULL) if (fp == NULL)
return NULL; return NULL;
// Count number of lines, allocate space and initiate counter // Count number of lines, allocate space and initiate counter
char ch; char ch;
int n = 0; int n = 0;
*cnt = 0; // Total number of graphs *cnt = 0; // Total number of graphs
*maxdim = 0; // Maximum number of vertices *maxdim = 0; // Maximum number of vertices
while(!feof(fp)) { while (!feof(fp)) {
ch = fgetc(fp); ch = fgetc(fp);
switch(ch) { switch (ch) {
case '\n': case '\n':
(*cnt)++; (*cnt)++;
n = 0; n = 0;
@ -57,10 +58,13 @@ int* read_graphs(char* filename, int *cnt, int *maxdim) {
rewind(fp); // Return to start of file rewind(fp); // Return to start of file
// Allocate space for graph data // Allocate space for graph data
const int gdatasize = *maxdim * *maxdim + 1; // Size allocation per graph // Size allocation per graph
const int arraysize = (1+ *cnt) * gdatasize; // Total space to allocate const int gdatasize = *maxdim * *maxdim + 1;
int *graphs = (int*)malloc(sizeof(int)*arraysize); // Allocation // Total space to allocate
if(!graphs) { 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"); fprintf(stderr, "ERROR: Could not allocate enough memory to store the graphs.\n");
return NULL; return NULL;
} }
@ -72,8 +76,8 @@ int* read_graphs(char* filename, int *cnt, int *maxdim) {
int j = 0; // Running index int j = 0; // Running index
while ((ch = fgetc(fp)) != EOF) { while ((ch = fgetc(fp)) != EOF) {
// New line encountered // New line encountered
if(ch == '\n') { if (ch == '\n') {
if(row != n) { if (row != n) {
fprintf(stderr, "ERROR: File not properly formatted.\n"); fprintf(stderr, "ERROR: File not properly formatted.\n");
return NULL; return NULL;
} }
@ -84,21 +88,21 @@ int* read_graphs(char* filename, int *cnt, int *maxdim) {
row = 0; row = 0;
n = -1; n = -1;
// Matrix entry encountered // Matrix entry encountered
} else if(ch == '0' || ch == '1') { } else if (ch == '0' || ch == '1') {
int x = ch - '0'; int x = ch - '0';
graphs[i*gdatasize + 1 + j] = x; graphs[i*gdatasize + 1 + j] = x;
col++; col++;
j++; j++;
// End-of-row or end-of-matrix encountered // End-of-row or end-of-matrix encountered
} else if(ch == '}') { } else if (ch == '}') {
if(n == -1) { if (n == -1) {
n = col; n = col;
} else if(col > 0 && n != col) { } else if (col > 0 && n != col) {
fprintf(stderr, "ERROR: File not properly formatted.\n"); fprintf(stderr, "ERROR: File not properly formatted.\n");
return NULL; return NULL;
} }
// End of row // End of row
if(col != 0) { if (col != 0) {
row++; row++;
col = 0; col = 0;
} }
@ -114,10 +118,10 @@ int* read_graphs(char* filename, int *cnt, int *maxdim) {
* Returns 0 if we run out of interventions, and 1 otherwise. * Returns 0 if we run out of interventions, and 1 otherwise.
***/ ***/
int nextintervention(int n, int *interventionslen, int *interventionidx) { int nextintervention(int n, int *interventionslen, int *interventionidx) {
for(int i=0; i<n; i++) { for (int i=0; i < n; i++) {
if(interventionidx[i] < interventionslen[i] - 1) { if (interventionidx[i] < interventionslen[i] - 1) {
interventionidx[i]++; interventionidx[i]++;
for(int j=0; j<i; j++) for (int j=0; j < i; j++)
interventionidx[j] = 0; interventionidx[j] = 0;
return 1; return 1;
} }
@ -132,15 +136,18 @@ int nextintervention(int n, int *interventionslen, int *interventionidx) {
* If `party' has no parents, then `party' will definitely receive a 1. * If `party' has no parents, then `party' will definitely receive a 1.
* If the fixed-point value cannot be inferred, return -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) { int alphapre(int party, const int *parents, const int *parentslen, \
if(parentslen[party] == 0) const int *interventions, const int *interventionidx, int maxn) {
if (parentslen[party] == 0)
return 1; return 1;
const int idatasize = (maxn+1)*(maxn+1)*2; // Max nr of interventions/party const int idatasize = (maxn+1)*(maxn+1)*2; // Max nr of interventions/party
for(int pidx=0; pidx<parentslen[party]; pidx++) { for (int pidx = 0; pidx < parentslen[party]; pidx++) {
const int parent = parents[maxn*party + pidx]; const int parent = parents[maxn*party + pidx];
const int f0 = interventions[idatasize*parent + 2*interventionidx[parent] + 0]; const int f0 = interventions[idatasize*parent + \
const int f1 = interventions[idatasize*parent + 2*interventionidx[parent] + 1]; 2*interventionidx[parent] + 0];
if(f0 != party && f1 != party) const int f1 = interventions[idatasize*parent + \
2*interventionidx[parent] + 1];
if (f0 != party && f1 != party)
return 0; return 0;
} }
return -1; return -1;
@ -154,34 +161,40 @@ int alphapre(int party, const int *parents, const int *parentslen, const int *in
* *
* Make sure to run alphapre before running this. * 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) { 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]; const int party = path[pathlen-1];
for(int i=0;i<pathlen-1;i++) { for (int i = 0; i < pathlen-1; i++) {
if(party == path[i]) { if (party == path[i]) {
return 0; return 0;
} }
} }
const int idatasize = (maxn+1)*(maxn+1)*2; // Max nr of interventions/party const int idatasize = (maxn+1)*(maxn+1)*2; // Max nr of interventions/party
// Iterate over all parents // Iterate over all parents
// If one parent does not vote for `party', immediattely return a zero. // 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. // If all parents vote for `party' --- which can only be inferred after we
for(int pidx=0; pidx<parentslen[party]; pidx++) { // queried all parents --- return a one.
for (int pidx = 0; pidx < parentslen[party]; pidx++) {
const int parent = parents[maxn*party + pidx]; const int parent = parents[maxn*party + pidx];
const int f0 = interventions[idatasize*parent + 2*interventionidx[parent] + 0]; const int f0 = interventions[idatasize*parent + 2*interventionidx[parent] \
const int f1 = interventions[idatasize*parent + 2*interventionidx[parent] + 1]; + 0];
if(f0 == party && f1 == party) const int f1 = interventions[idatasize*parent + 2*interventionidx[parent] \
+ 1];
if (f0 == party && f1 == party)
continue; continue;
// Re-use already computed fixed points, else enter recursion // Re-use already computed fixed points, else enter recursion
int val; int val;
if(fp[parent] != -1) { if (fp[parent] != -1) {
val = fp[parent]; val = fp[parent];
} else { } else {
path[pathlen] = parent; path[pathlen] = parent;
val = alpha(path, pathlen+1, parents, parentslen, interventions, interventionidx, fp, n, maxn); val = alpha(path, pathlen+1, parents, parentslen, interventions, \
interventionidx, fp, n, maxn);
} }
// Evaluate the intervention on the input // Evaluate the intervention on the input
const int muOfVal = (val == 0) ? f0 : f1; const int muOfVal = (val == 0) ? f0 : f1;
if(muOfVal != party) if (muOfVal != party)
return 0; return 0;
} }
return 1; return 1;
@ -193,7 +206,7 @@ int alpha(int *path, int pathlen, const int *parents, const int *parentslen, con
***/ ***/
int ceillog10(int x) { int ceillog10(int x) {
int res = 0; int res = 0;
while(x > 0) { while (x > 0) {
x /= 10; x /= 10;
res++; res++;
} }
@ -213,10 +226,15 @@ int ceillog10(int x) {
***/ ***/
int test_graphs(int startidx, int len, int maxn, const int *graphs) { int test_graphs(int startidx, int len, int maxn, const int *graphs) {
const int gdatasize = maxn * maxn + 1; // Size allocation per graph const int gdatasize = maxn * maxn + 1; // Size allocation per graph
const int idatasize = (maxn+1)*(maxn+1)*2; // Max number of interventions per party const int idatasize = (maxn+1)*(maxn+1)*2; // Max number of interventions
// per party
int range[maxn+1]; // Range of intervention 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 interventions[maxn*idatasize]; // Hold all interventions per
int interventionslen[maxn]; // Number of interventions per party // 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 interventionidx[maxn]; // Current intervention executed
int parents[maxn*maxn]; // List all parents per party int parents[maxn*maxn]; // List all parents per party
int parentslen[maxn]; // Number of parents per party int parentslen[maxn]; // Number of parents per party
@ -225,14 +243,15 @@ int test_graphs(int startidx, int len, int maxn, const int *graphs) {
const int padlen = ceillog10(len); const int padlen = ceillog10(len);
time_t t0 = time(NULL); time_t t0 = time(NULL);
const time_t tstart = t0; const time_t tstart = t0;
for(int graphidx=startidx; graphidx<startidx+len; graphidx++) { for (int graphidx = startidx; graphidx < startidx+len; graphidx++) {
time_t t1 = time(NULL); time_t t1 = time(NULL);
if(t1 > t0) { if (t1 > t0) {
const int count = graphidx - startidx + 1; const int count = graphidx - startidx + 1;
const float percentage = 100*(float)count/(float)len; const float percentage = 100*(float)count/(float)len;
const int deltat = t1-tstart; const int deltat = t1-tstart;
const float rate = (float)count/(float)deltat; 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); 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; t0 = t1;
} }
// Test graph with index `graphidx' // Test graph with index `graphidx'
@ -240,14 +259,14 @@ int test_graphs(int startidx, int len, int maxn, const int *graphs) {
const int n = graphs[graphidx*gdatasize]; const int n = graphs[graphidx*gdatasize];
// Get parents and children // Get parents and children
// Adj matrix is such that (row,col) = 1 <=> edge from row to col // Adj matrix is such that (row,col) = 1 <=> edge from row to col
for(int i=0; i<n; i++) { for (int i = 0; i < n; i++) {
parentslen[i] = 0; parentslen[i] = 0;
childrenlen[i] = 0; childrenlen[i] = 0;
} }
int idx = 0; int idx = 0;
for(int row=0; row<n; row++) { for (int row = 0; row < n; row++) {
for(int col=0; col<n; col++) { for (int col = 0; col < n; col++) {
if(graphs[graphidx*gdatasize+1+idx] == 1) { if (graphs[graphidx*gdatasize+1+idx] == 1) {
// Edge from row to col // Edge from row to col
children[maxn*row + childrenlen[row]] = col; children[maxn*row + childrenlen[row]] = col;
childrenlen[row]++; childrenlen[row]++;
@ -258,22 +277,24 @@ int test_graphs(int startidx, int len, int maxn, const int *graphs) {
} }
} }
// Fill in interventions // Fill in interventions
for(int party=0; party<n; party++) { for (int party = 0; party < n; party++) {
// This is a small trick to speed-up processing: // 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), // 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). // else make use of the `discard' internvetion (tmp=1).
// Populate range of `party's interventions // Populate range of `party's interventions
// The range has childrenlen[party]+tmp entries // The range has childrenlen[party]+tmp entries
range[0] = -1; range[0] = -1;
const int tmp = (childrenlen[party]<=1) ? 1 : 0; const int tmp = (childrenlen[party] <= 1) ? 1 : 0;
for(int j=0; j<childrenlen[party]; j++) { for (int j = 0; j < childrenlen[party]; j++) {
range[j+tmp] = children[maxn*party+j]; range[j+tmp] = children[maxn*party+j];
} }
// Enumerate, store, and count all possible interventions // 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 // Each intervention is a tuple (mu(0), mu(1)), x corresponds to the first
// entry, y to the second
int cnt = 0; int cnt = 0;
for(int x=0; x<childrenlen[party]+tmp; x++) { for (int x = 0; x < childrenlen[party]+tmp; x++) {
for(int y=0; y<childrenlen[party]+tmp; y++) { for (int y = 0; y < childrenlen[party]+tmp; y++) {
interventions[party*idatasize+2*cnt] = range[x]; interventions[party*idatasize+2*cnt] = range[x];
interventions[party*idatasize+2*cnt+1] = range[y]; interventions[party*idatasize+2*cnt+1] = range[y];
cnt++; cnt++;
@ -285,57 +306,66 @@ int test_graphs(int startidx, int len, int maxn, const int *graphs) {
int fp[n]; // Array to store alpha(party) values int fp[n]; // Array to store alpha(party) values
int afterfp[n]; // \omega(\mu(fp)) int afterfp[n]; // \omega(\mu(fp))
int path[n+1]; // Current path int path[n+1]; // Current path
for(int i=0; i<n; i++) interventionidx[i] = 0; // Initialize intervention index to the first intervention // Initialize intervention index to the first intervention
int running = 1; // Flag to see whether we are still verifying, or we verified all interventions for (int i = 0; i < n; i++) interventionidx[i] = 0;
while(running) { // Flag to see whether we are still verifying, or we verified all
// interventions
int running = 1;
while (running) {
// Reset fixed-point entries // Reset fixed-point entries
for(int party=0; party<n; party++) for (int party = 0; party < n; party++)
fp[party] = -1; fp[party] = -1;
// Fill out fixed-points that are trivially computed // Fill out fixed-points that are trivially computed
for(int party=0; party<n; party++) for (int party = 0; party < n; party++)
alphapre(party, parents, parentslen, interventions, interventionidx, maxn); alphapre(party, parents, parentslen, interventions, interventionidx, \
maxn);
// Invoke recursive function // Invoke recursive function
for(int party=0; party<n; party++) { for (int party = 0; party < n; party++) {
path[0] = party; path[0] = party;
const int val = alpha(path, 1, parents, parentslen, interventions, interventionidx, fp, n, maxn); const int val = alpha(path, 1, parents, parentslen, interventions, \
interventionidx, fp, n, maxn);
fp[party] = val; fp[party] = val;
} }
// Compute \omega(\mu(fp)) // Compute \omega(\mu(fp))
for(int i=0; i<n; i++) for (int i = 0; i < n; i++)
afterfp[i] = -1; afterfp[i] = -1;
// Loop over parties // Loop over parties
for(int party=0; party<n; party++) { for (int party = 0; party < n; party++) {
// Loop over parents of this party // Loop over parents of this party
for(int pidx=0; pidx<parentslen[party]; pidx++) { for (int pidx = 0; pidx < parentslen[party]; pidx++) {
const int parent = parents[maxn*party + pidx]; const int parent = parents[maxn*party + pidx];
const int f0 = interventions[idatasize*parent + 2*interventionidx[parent] + 0]; const int f0 = interventions[idatasize*parent + \
const int f1 = interventions[idatasize*parent + 2*interventionidx[parent] + 1]; 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 `parent' never votes for `party', this value is 0.
if(f0 != party && f1 != party) { if (f0 != party && f1 != party) {
afterfp[party] = 0; afterfp[party] = 0;
break; break;
} }
// If `parent' always votes for `party', check the other parents. // If `parent' always votes for `party', check the other parents.
if(f0 == party && f1 == party) if (f0 == party && f1 == party)
continue; continue;
// Evaluate intervention on the input (fp[parent]) // Evaluate intervention on the input (fp[parent])
const int muOfVal = (fp[parent] == 0) ? f0 : f1; 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 `parent' does not vote for `party', this value is 0, else check
if(muOfVal != party) { // the other parents.
if (muOfVal != party) {
afterfp[party] = 0; afterfp[party] = 0;
break; break;
} }
} }
// All parents voted `party' // All parents voted `party'
if(afterfp[party] == -1) if (afterfp[party] == -1)
afterfp[party] = 1; afterfp[party] = 1;
} }
// Compare the output of the recursive alpha function (fp) with \omega(\mu(fp)) (afterfp) // Compare the output of the recursive alpha function (fp) with
// If these two arrays differ, return the line number of the graph we falsified // \omega(\mu(fp)) (afterfp)
for(int party=0;party<n;party++) // If these two arrays differ, return the line number of the graph we
if(fp[party] != afterfp[party]) // falsified
for (int party = 0; party < n; party++)
if (fp[party] != afterfp[party])
return graphidx+1; return graphidx+1;
// Next intervention // Next intervention
@ -344,7 +374,8 @@ int test_graphs(int startidx, int len, int maxn, const int *graphs) {
} }
const int deltat = time(NULL) - tstart; const int deltat = time(NULL) - tstart;
const float rate = (deltat == 0) ? len : (float)len/(float)deltat; 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); 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; return 0;
} }
@ -352,15 +383,16 @@ int main(int argc, char *argv[]) {
// Parse command-line arguments // Parse command-line arguments
int START = 0; int START = 0;
int NUM = 0; int NUM = 0;
if(argc >= 3) if (argc >= 3)
START = atoi(argv[2]); START = atoi(argv[2]);
if(argc >= 4) { if (argc >= 4) {
if(argv[3][0] == '+') if (argv[3][0] == '+')
NUM = atoi(argv[3]+1); NUM = atoi(argv[3]+1);
else else
NUM = atoi(argv[3])-START+1; NUM = atoi(argv[3])-START+1;
} }
if(!(argc >= 2 && (argc < 3 || START > 0) && (argc < 4 || NUM > 0) && argc <= 4)) { 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, "Usage: %s <filename> [<startline> [<endline> | +<count>]]\n", argv[0]);
fprintf(stderr, " <filename> File name with adjacency matrices of simple directed graphs\n"); 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, " <startline> Verify graphs starting from line `startline'\n");
@ -380,25 +412,27 @@ int main(int argc, char *argv[]) {
int graphcount = 0; int graphcount = 0;
int maxdimension = -1; int maxdimension = -1;
const int *graphs = read_graphs(argv[1], &graphcount, &maxdimension); const int *graphs = read_graphs(argv[1], &graphcount, &maxdimension);
if(graphs == NULL) { if (graphs == NULL) {
fprintf(stderr, "ERROR: Could not read graphs from file `%s'.\n", argv[1]); fprintf(stderr, "ERROR: Could not read graphs from file `%s'.\n", argv[1]);
return -1; return -1;
} }
const int startidx = START ? START-1 : 0; // Index of graph to start verifying 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 const int tot = NUM ? NUM : graphcount-startidx; // Number of graphs to be verified
if(startidx >= graphcount || tot <= 0 || tot+startidx-1 >= graphcount) { if (startidx >= graphcount || tot <= 0 || tot+startidx-1 >= graphcount) {
fprintf(stderr, "ERROR: Lines to be verified are out of range\n"); fprintf(stderr, "ERROR: Lines to be verified are out of range\n");
return -1; return -1;
} }
// Test `tot' graphs, starting from index `startidx' // 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); 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); const int result = test_graphs(startidx, tot, maxdimension, graphs);
if (result == 0) if (result == 0)
printf("These graphs are admissible\n"); printf("These graphs are admissible\n");
else if(result == -1) else if (result == -1)
fprintf(stderr, "ERROR: Something went wrong\n"); fprintf(stderr, "ERROR: Something went wrong\n");
else else
printf("The function alpha does not represent the fixed point, or the graph on line %d is inadmissible\n", result); printf("The function alpha does not represent the fixed point, or the graph on line %d is inadmissible\n", \
result);
return result; return result;
} }

374
src/SOCgen.c
View File

@ -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 // SPDX-License-Identifier: GPL-3.0-or-later
@ -12,20 +13,20 @@
***/ ***/
void dumpgraph(int n, const int *children, const int *childrenlen) { void dumpgraph(int n, const int *children, const int *childrenlen) {
int adj[n*n]; int adj[n*n];
for(int i=0; i<n*n; i++) adj[i] = 0; for (int i = 0; i < n*n; i++) adj[i] = 0;
for(int a=0; a<n; a++) { for (int a = 0; a < n; a++) {
for(int bidx=0; bidx<childrenlen[a]; bidx++) { for (int bidx = 0; bidx < childrenlen[a]; bidx++) {
const int b = children[a*n+bidx]; const int b = children[a*n+bidx];
adj[a*n+b] = 1; adj[a*n+b] = 1;
} }
} }
// Print // Print
printf("{"); printf("{");
for(int a=0; a<n; a++) { for (int a = 0; a < n; a++) {
if(a) printf(","); if (a) printf(",");
printf("{"); printf("{");
printf("%d", adj[a*n+0]); printf("%d", adj[a*n+0]);
for(int b=1; b<n; b++) printf(",%d", adj[a*n+b]); for (int b = 1; b < n; b++) printf(",%d", adj[a*n+b]);
printf("}"); printf("}");
} }
printf("}\n"); printf("}\n");
@ -36,11 +37,11 @@ void dumpgraph(int n, const int *children, const int *childrenlen) {
***/ ***/
void dumpgv(int n, int graphnr, const int *children, const int *childrenlen) { void dumpgv(int n, int graphnr, const int *children, const int *childrenlen) {
printf("strict digraph G%d {node[shape=point];", graphnr); printf("strict digraph G%d {node[shape=point];", graphnr);
for(int v=0; v<n; v++) { for (int v = 0; v < n; v++) {
for(int i=0; i<childrenlen[v]; i++) { for (int i = 0; i < childrenlen[v]; i++) {
const int w = children[v*n+i]; const int w = children[v*n+i];
printf("G%dN%d -> G%dN%d", graphnr, v, graphnr, w); printf("G%dN%d -> G%dN%d", graphnr, v, graphnr, w);
if(childrenlen[v] >= 2) if (childrenlen[v] >= 2)
printf(" [color=\"red\"]"); printf(" [color=\"red\"]");
printf(";"); printf(";");
} }
@ -53,8 +54,8 @@ void dumpgv(int n, int graphnr, const int *children, const int *childrenlen) {
* Block size of our datastructure. * Block size of our datastructure.
***/ ***/
int blocksize(int n) { int blocksize(int n) {
int res=(n+1)*n; int res = (n+1)*n;
for(int i=0; i<n; i++) for (int i = 0; i < n; i++)
res *= n-i; res *= n-i;
return res; return res;
} }
@ -62,16 +63,17 @@ int blocksize(int n) {
/*** /***
* Recrusive bi-directional graph traversal, used to check connectivity. * 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) { void traverse(int vertex, int *visited, int n, const int *parents, \
const int *parentslen, const int *children, const int *childrenlen) {
visited[vertex] = 1; visited[vertex] = 1;
for(int i=0; i<parentslen[vertex]; i++) { for (int i = 0; i < parentslen[vertex]; i++) {
int u = parents[n*vertex + i]; int u = parents[n*vertex + i];
if(visited[u] == 0) if (visited[u] == 0)
traverse(u, visited, n, parents, parentslen, children, childrenlen); traverse(u, visited, n, parents, parentslen, children, childrenlen);
} }
for(int i=0; i<childrenlen[vertex]; i++) { for (int i = 0; i < childrenlen[vertex]; i++) {
int u = children[n*vertex + i]; int u = children[n*vertex + i];
if(visited[u] == 0) if (visited[u] == 0)
traverse(u, visited, n, parents, parentslen, children, childrenlen); traverse(u, visited, n, parents, parentslen, children, childrenlen);
} }
} }
@ -79,15 +81,16 @@ void traverse(int vertex, int *visited, int n, const int *parents, const int *pa
/*** /***
* Return 1 if the graph is connected, and 0 otherwise. * 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) { 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? // 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; for (int i = 0; i < n; i++) if (parentslen[i]+childrenlen[i] == 0) return 0;
// Check connectivity // Check connectivity
int visited[n]; int visited[n];
for(int i=0; i<n; i++) visited[i] = 0; for (int i = 0; i < n; i++) visited[i] = 0;
// Start traversal // Start traversal
traverse(0, visited, n, parents, parentslen, children, childrenlen); traverse(0, visited, n, parents, parentslen, children, childrenlen);
for(int i=0; i<n; i++) if(!visited[i]) return 0; for (int i = 0; i < n; i++) if (!visited[i]) return 0;
return 1; return 1;
} }
@ -95,31 +98,31 @@ int isconnected(int n, const int *parents, const int *parentslen, const int *chi
* Interpret graph number `graph' as a simple digraph, and populate the arrays `parents', `parentslen', `children', `childrenlen' * 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). * 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) { void graphnrtolists(int n, unsigned long long int graph, int *parents, \
for(int i=0; i<n; i++) { int *parentslen, int *children, int *childrenlen) {
for (int i = 0; i < n; i++) {
parentslen[i] = 0; parentslen[i] = 0;
childrenlen[i] = 0; childrenlen[i] = 0;
} }
int row = 0; int row = 0;
int col = 1; // We cannot start at (0,0); it does not exist int col = 1; // We cannot start at (0,0); it does not exist
unsigned long long g = graph; unsigned long long g = graph;
for(int i=0; i<n*(n-1); i++) { for (int i = 0; i < n*(n-1); i++) {
int val = g%2; int val = g%2;
g = g/2; g = g/2;
if(val) { if (val) {
children[row*n+childrenlen[row]] = col; children[row*n+childrenlen[row]] = col;
childrenlen[row]++; childrenlen[row]++;
parents[col*n+parentslen[col]] = row; parents[col*n+parentslen[col]] = row;
parentslen[col]++; parentslen[col]++;
} }
col++; col++;
if(col==row) if (col == row)
col++; col++;
if(col==n) { if (col == n) {
col = 0; col = 0;
row++; row++;
} }
} }
} }
@ -127,20 +130,24 @@ void graphnrtolists(int n, unsigned long long int graph, int *parents, int *pare
* Find the tail in `path' of length `pathlen' that includes a cycle, and save the canonical representation. * 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. * 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) { void canonical_cycle(int n, int *cycle, int *cyclelen, const int *path, \
int pathlen) {
int startidx; int startidx;
int pathhascycle = 0; 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++) { for (int k = 0; k < pathlen-1; k++) {
if(path[k]==last) { // We have found the first occurance of node `last' on the path if (path[k] == last) {
// The cycle is the following part of path: path[startidx], path[startidx+1], ..., path[pathlen-2], and therefore has length pathlen-startidx-1 // 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; startidx = k;
(*cyclelen) = pathlen-startidx-1; (*cyclelen) = pathlen-startidx-1;
pathhascycle = 1; pathhascycle = 1;
break; break;
} }
} }
if(!pathhascycle) { // Sanity check. This should NEVER happen! 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"); 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; return;
} }
@ -148,35 +155,43 @@ void canonical_cycle(int n, int *cycle, int *cyclelen, const int *path, int path
// First, find smallest node on cycle // First, find smallest node on cycle
int smallestidx = -1; int smallestidx = -1;
int smallest = n; int smallest = n;
for(int k=startidx; k<pathlen-1; k++) { for (int k = startidx; k < pathlen-1; k++) {
if(smallest > path[k]) { if (smallest > path[k]) {
smallest = path[k]; smallest = path[k];
smallestidx = k; smallestidx = k;
} }
} }
// Now, copy cycle to `cycle' starting from the smallest // Now, copy cycle to `cycle' starting from the smallest
int j=0; // This is the index on where to write on the array `cycle' // This is the index on where to write on the array `cycle'
for(int i=smallestidx;; i++) { int j = 0;
if(j==*cyclelen) break; // Stop filling in, if we reached the end for (int i = smallestidx;; i++) {
if(i==pathlen-1) i = startidx; // If we reach the end of the path, jump to the start // Stop filling in, if we reached the end
cycle[j] = path[i]; // Write at position `j' the element from the cycle on the path at position `i`' if (j == *cyclelen) break;
j++; // Increase index to write on `cycle' // 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. * 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) { int cycleisnew(int n, const int *c, int clen, const int *cycles, \
const int *cyclescnt) {
const int cnt = cyclescnt[clen]; const int cnt = cyclescnt[clen];
const int bs = blocksize(n); const int bs = blocksize(n);
for(int i=0; i<cnt; i++) { // Iterate over all cycles of the same length // Iterate over all cycles of the same length
for (int i = 0; i < cnt; i++) {
int j; int j;
const int startidx = clen*bs + i*clen; // We start here const int startidx = clen*bs + i*clen; // We start here
for(j=0; j<clen; j++) // Go position by position for (j = 0; j < clen; j++) // Go position by position
if(cycles[startidx+j] != c[j]) if (cycles[startidx+j] != c[j])
break; // They differ at position j break; // They differ at position j
if(j==clen) // We didn't break if (j == clen) // We didn't break
return 0; return 0;
} }
return 1; return 1;
@ -186,11 +201,12 @@ int cycleisnew(int n, const int *c, int clen, const int *cycles, const int *cycl
* Save the given cycle `c' of length `clen' in `cycles'. * Save the given cycle `c' of length `clen' in `cycles'.
***/ ***/
void savecycle(int n, const int *c, int clen, int *cycles, int *cyclescnt) { 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 m = cyclescnt[clen]; // There are `m' `clen'-cycles
const int bs = blocksize(n); const int bs = blocksize(n);
const int startidx = clen*bs + m*clen; // We start writing here const int startidx = clen*bs + m*clen; // We start writing here
for(int i=0; i<clen; i++) cycles[startidx+i] = c[i]; 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' cyclescnt[clen]++; // Increase the counter that holds
// the number of `clen'-cycles
} }
@ -204,32 +220,43 @@ void savecycle(int n, const int *c, int clen, int *cycles, int *cyclescnt) {
* Parameter `n': Number of 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) { 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' // We go on duing a dfs from node `vertex'
for(int i=0; i<childrenlen[vertex]; i++) { const int vertex = path[pathlen-1];
for (int i = 0; i < childrenlen[vertex]; i++) {
const int u = children[n*vertex + i]; // Enter node `u' const int u = children[n*vertex + i]; // Enter node `u'
if(useless[u]) continue; // This one is useless if (useless[u]) continue; // This one is useless
if(visited[u] == 0) { // Not visited yet if (visited[u] == 0) { // Not visited yet
visited[u] = 1; // Mark as visited visited[u] = 1; // Mark as visited
path[pathlen] = u; // Append to path, and re-enter recursion... path[pathlen] = u; // Append to path, and re-enter
dfs(useless, path, pathlen+1, visited, n, children, childrenlen, cycles, cyclescnt, found); // recursion...
visited[u] = 0; // Undo marking, we are going to take another route 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 { } else {
// We have reached a node on the path, i.e., we have found a cycle! // We have reached a node on the path, i.e., we have found a cycle!
// The cycle is somewhere along the path // The cycle is somewhere along the path
// path[0], path[1], ... path[pathlen-1] // path[0], path[1], ... path[pathlen-1]
// More specifically, node `u' is somewhere there. // More specifically, node `u' is somewhere there.
// To detect this, we first place `u' on the path. // 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. // The method `canonical_cycle' will then use this information to find
// and return the canonical form of the cycle.
path[pathlen] = u; path[pathlen] = u;
canonical_cycle(n, cycle, &cyclelen, path, pathlen+1); // Find cycle and generate the canonical representation. // Find cycle and generate the canonical representation.
if(cycleisnew(n, cycle, cyclelen, cycles, cyclescnt)) { // Is it new? If yes, save canonical_cycle(n, cycle, &cyclelen, path, pathlen+1);
(*found)++; // Increase the number of found cycles by 1 // Is it new? If yes, save
savecycle(n, cycle, cyclelen, cycles, cyclescnt); // Save cycle 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.. // else dont do anything...
} }
} }
return; return;
@ -240,26 +267,33 @@ void dfs(const int *useless, int *path, int pathlen, int *visited, int n, const
* The cycles are stored in `cycles' and `cyclescnt'. * The cycles are stored in `cycles' and `cyclescnt'.
* This uses a recursive depth-first-search. * This uses a recursive depth-first-search.
***/ ***/
int find_cycles(int *cycles, int *cyclescnt, int n, const int *children, const int *childrenlen) { int find_cycles(int *cycles, int *cyclescnt, int n, const int *children, \
const int *childrenlen) {
int visited[n]; int visited[n];
int useless[n]; // Some nodes are useless, mark them as such // 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. 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 unsigned int found = 0; // Store total number of cycles found
// Initialize cyclescnt to zero // Initialize cyclescnt to zero
for(int i=0; i<n+1; i++) for (int i = 0; i < n+1; i++)
cyclescnt[i] = 0; cyclescnt[i] = 0;
// Mark nodes without children as useless // Mark nodes without children as useless
for(int i=0; i<n; i++) { for (int i = 0; i < n; i++) {
useless[i] = (childrenlen[i]==0) ? 1 : 0; useless[i] = (childrenlen[i] == 0) ? 1 : 0;
} }
// Iterate over nodes // Iterate over nodes
for(int vertex=0; vertex<n; vertex++) { for (int vertex = 0; vertex < n; vertex++) {
if(useless[vertex]) continue; if (useless[vertex]) continue;
for(int i=0; i<n; i++) visited[i] = 0; // Mark all as univisited for (int i = 0; i < n; i++) visited[i] = 0; // Mark all as univisited
path[0] = vertex; // We start at `vertex' // We start at `vertex'
visited[vertex] = 1; // Mark `vertex' as visited. If we revisit a node, we found a cycle. path[0] = vertex;
// Mark `vertex' as visited. If we revisit a node, we found a cycle.
visited[vertex] = 1;
// Enter depth-first-search recursion // Enter depth-first-search recursion
dfs(useless, path, 1, visited, n, children, childrenlen, cycles, cyclescnt, &found); dfs(useless, path, 1, visited, n, children, childrenlen, cycles, \
cyclescnt, &found);
// We processed this. Every cycle that goes through `vertex' has been found. // We processed this. Every cycle that goes through `vertex' has been found.
useless[vertex] = 1; useless[vertex] = 1;
} }
@ -280,22 +314,27 @@ 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. * 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) { int gissoc(int n, const int *parents, const int *parentslen, \
// 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) 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. // 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[0], cycles[1] is the first cycle of length 2
// cycles[2], cycles[3] is the second cycle of length 2 // cycles[2], cycles[3] is the second cycle of length 2
// etc. // etc.
// There are at most (n choose 2) cycles of length 2. // 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. // 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. // 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)]. // So, the cycles with length k are saved at and after cycles[k*(2^n)].
// Each block has size 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 */ /* Early and saftey check for self loops. SOCs dont have self-loops, and we
for(int i=0; i<n; i++) { * should never generate them. Anyhow, test */
for (int i = 0; i < n; i++) {
int num = childrenlen[i]; int num = childrenlen[i];
for(int cidx=0; cidx<num; cidx++) for (int cidx = 0; cidx < num; cidx++)
if(i == children[n*i+cidx]) { if (i == children[n*i+cidx]) {
fprintf(stderr, "ERROR: THIS GRAPH HAS A SELF LOOP (node %d)\n", i); fprintf(stderr, "ERROR: THIS GRAPH HAS A SELF LOOP (node %d)\n", i);
return 0; return 0;
} }
@ -307,31 +346,42 @@ int gissoc(int n, const int *parents, const int *parentslen, const int *children
// then iterate over nodes of that cycle. // then iterate over nodes of that cycle.
// In that iteration process, mark all parents. // In that iteration process, mark all parents.
// If a marked parent is marked again, we have found a common parent. // If a marked parent is marked again, we have found a common parent.
for(int len=2; len<n+1; len++) { for (int len = 2; len < n+1; len++) {
const int cnt = cyclescnt[len]; // Number of cycles of length `len' // Number of cycles of length `len'
for(int i=0; i<cnt; i++) { // Iterate over cycles const int cnt = cyclescnt[len];
const int startidx = len*bs+len*i; // The cycle we look at (length=`len', index=`i') is saved starting from index `startidx' // 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 // Initialize marks
int marked[n]; int marked[n];
for(int k=0; k<n; k++) marked[k] = 0; for (int k = 0; k < n; k++) marked[k] = 0;
// Flag: Have we found a common parent? // Flag: Have we found a common parent?
int commonparentfound = 0; int commonparentfound = 0;
// Iterate over cycle // Iterate over cycle
for(int k=0; k<len; k++) { for (int k = 0; k < len; k++) {
const int vertex = cycles[startidx+k]; const int vertex = cycles[startidx+k];
for(int pidx=0; pidx<parentslen[vertex]; pidx++) { for (int pidx = 0; pidx < parentslen[vertex]; pidx++) {
const int parent = parents[n*vertex+pidx]; const int parent = parents[n*vertex+pidx];
if(marked[parent]) commonparentfound = 1; // The node `parent' is a common parent // Check if `parent' is a common parent, else mark
else marked[parent] = 1; // Mark the parent node if (marked[parent])
if(commonparentfound) break; // We have found a common parent, break out of the for loop over `pidx' commonparentfound = 1;
else
marked[parent] = 1;
// We have found a common parent, break out of the `pidx' loop
if (commonparentfound) break;
} }
if(commonparentfound) break; // We have found a common parent, break out of the for loop over `k' // We have found a common parent, break out of the for loop over `k'
if (commonparentfound) break;
} }
if(!commonparentfound) // On cycle `i' we did NOT find a common parent, this is not a SOC // On cycle `i' we did NOT find a common parent, this is not a SOC
if (!commonparentfound)
return 0; return 0;
} }
} }
return 1; // If all test pass, we have found a SOC // If all test pass, we have found a SOC
return 1;
} }
/*** /***
@ -339,13 +389,13 @@ int gissoc(int n, const int *parents, const int *parentslen, const int *children
***/ ***/
int grandparentslen(int node, const int *parents, const int *parentslen) { int grandparentslen(int node, const int *parents, const int *parentslen) {
int tot = 0; int tot = 0;
for(int i=0; i<parentslen[node]; i++) for (int i = 0; i < parentslen[node]; i++)
tot += parentslen[i]; tot += parentslen[i];
return tot; return tot;
} }
int grandchildrenlen(int node, const int *children, const int *childrenlen) { int grandchildrenlen(int node, const int *children, const int *childrenlen) {
int tot = 0; int tot = 0;
for(int i=0; i<childrenlen[node]; i++) for (int i = 0; i < childrenlen[node]; i++)
tot += childrenlen[i]; tot += childrenlen[i];
return tot; return tot;
} }
@ -354,19 +404,23 @@ int grandchildrenlen(int node, const int *children, const int *childrenlen) {
* If the graph does not satisfy this `canonical' form, then another ismorphic to it (it's simply a permutation of the nodes) * 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. * will be degree-ordered.
***/ ***/
int isdegreeordered(int n, const int *parents, const int *parentslen, const int *children, const int *childrenlen) { 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 ... // Decreasingling order by in-degree, then out-degree, then ...
for(int i=0; i<n-1; i++) { for (int i = 0; i < n-1; i++) {
if(parentslen[i] > parentslen[i+1]) { if (parentslen[i] > parentslen[i+1]) {
return 0; return 0;
} else if(parentslen[i] == parentslen[i+1]) { } else if (parentslen[i] == parentslen[i+1]) {
if(childrenlen[i] > childrenlen[i+1]) { if (childrenlen[i] > childrenlen[i+1]) {
return 0; return 0;
} else if(childrenlen[i] == childrenlen[i+1]) { } else if (childrenlen[i] == childrenlen[i+1]) {
if(grandparentslen(i,parents,parentslen) > grandparentslen(i+1,parents,parentslen)) { if (grandparentslen(i, parents, parentslen) > \
grandparentslen(i+1, parents, parentslen)) {
return 0; return 0;
} else if(grandparentslen(i,parents,parentslen) == grandparentslen(i+1,parents,parentslen)) { } else if (grandparentslen(i, parents, parentslen) == \
if(grandchildrenlen(i,children,childrenlen) > grandchildrenlen(i+1,children,childrenlen)) { grandparentslen(i+1, parents, parentslen)) {
if (grandchildrenlen(i, children, childrenlen) > \
grandchildrenlen(i+1, children, childrenlen)) {
return 0; return 0;
} }
} }
@ -407,27 +461,27 @@ int main(int argc, char *argv[]) {
int RANDOM = 0; int RANDOM = 0;
int GRAPHVIZ = 0; int GRAPHVIZ = 0;
int UNKOPTION = 0; int UNKOPTION = 0;
for(int i=1; i<argc; i++) { for (int i = 1; i < argc; i++) {
if(strcmp(argv[i], "-n") == 0 && i+1 < argc) { if (strcmp(argv[i], "-n") == 0 && i+1 < argc) {
n = atoi(argv[i+1]); n = atoi(argv[i+1]);
i++; i++;
} else if(strcmp(argv[i], "-r") == 0 && i+1 < argc) { } else if (strcmp(argv[i], "-r") == 0 && i+1 < argc) {
RANDOM = atoi(argv[i+1]); RANDOM = atoi(argv[i+1]);
i++; i++;
} else if(strcmp(argv[i], "-c") == 0) } else if (strcmp(argv[i], "-c") == 0) {
NONDAGONLY = 1; NONDAGONLY = 1;
else if(strcmp(argv[i], "--no-sink") == 0) } else if (strcmp(argv[i], "--no-sink") == 0) {
NOSINK = 1; NOSINK = 1;
else if(strcmp(argv[i], "--no-source") == 0) } else if (strcmp(argv[i], "--no-source") == 0) {
NOSOURCE = 1; NOSOURCE = 1;
else if(strcmp(argv[i], "--graphviz") == 0) } else if (strcmp(argv[i], "--graphviz") == 0) {
GRAPHVIZ = 1; GRAPHVIZ = 1;
else { } else {
UNKOPTION = 1; UNKOPTION = 1;
break; break;
} }
} }
if(UNKOPTION || n<=1) { if (UNKOPTION || n <= 1) {
fprintf(stderr, "Usage: %s -n <order> [-r <num>] [--grpahviz] [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, " -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, " -r <num> Pick directed graphs at random, and exit after having found `num' SOCs\n");
@ -450,16 +504,18 @@ int main(int argc, char *argv[]) {
int childrenlen[n]; // Number of children of node i is stored in childrenlen[i] int childrenlen[n]; // Number of children of node i is stored in childrenlen[i]
const int bs = blocksize(n); const int bs = blocksize(n);
int *cycles = (int*)malloc(sizeof(int)*(n+2)*bs); int *cycles = (int*)malloc(sizeof(int)*(n+2)*bs);
if(!cycles) { if (!cycles) {
fprintf(stderr, "Failed to allocate memory to store cycles (tried to allocate %lu bytes)\n", sizeof(int)*(n+2)*bs); fprintf(stderr, "Failed to allocate memory to store cycles (tried to allocate %lu bytes)\n", \
sizeof(int)*(n+2)*bs);
return -1; return -1;
} }
int cyclescnt[n+1]; // cyclescnt[k] stores the number of cycles with length k int cyclescnt[n+1]; // cyclescnt[k] stores the number of cycles with length k
// EO Setup datastructures // EO Setup datastructures
// Initiate graph enumeration // Initiate graph enumeration
const int m = n*(n-1); const int m = n*(n-1);
if(m>64) { if (m > 64) {
fprintf(stderr, "Too many graphs to enumarate with an unsiged long long integer (number of node pairs = %d)\n", m); fprintf(stderr, "Too many graphs to enumarate with an unsiged long long integer (number of node pairs = %d)\n", \
m);
return -1; return -1;
} }
unsigned long long max = 1L << m; // Largest `graphnumber' unsigned long long max = 1L << m; // Largest `graphnumber'
@ -470,74 +526,92 @@ int main(int argc, char *argv[]) {
srand((unsigned) time(&t)); srand((unsigned) time(&t));
// EO Initiate graph enumeration // EO Initiate graph enumeration
fprintf(stderr, "Generating SOCs with %d nodes\n", n); fprintf(stderr, "Generating SOCs with %d nodes\n", n);
if(RANDOM) fprintf(stderr, " Picking graphs at random\n"); if (RANDOM)
if(NONDAGONLY) fprintf(stderr, " Filter: Omitting DAGs\n"); fprintf(stderr, " Picking graphs at random\n");
if(NOSINK) fprintf(stderr, " Filter: Omitting graphs with sink nodes\n"); if (NONDAGONLY)
if(NOSOURCE) fprintf(stderr, " Filter: Omitting graphs with source nodes\n"); fprintf(stderr, " Filter: Omitting DAGs\n");
unsigned long long graphnumber = 0; // We will always omit the graph with `graphnumber' 0; it is unintersting anyway 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; unsigned long long graphschecked = 0;
for(;graphnumber < max && (!RANDOM || len < RANDOM);) { for (; graphnumber < max && (!RANDOM || len < RANDOM);) {
const time_t now = time(NULL); 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 // Print status every second
if(now > t) { if (now > t) {
t = now; t = now;
const int deltat = now - t0; const int deltat = now - t0;
const float rate = (float)graphschecked/(float)deltat; const float rate = (float)graphschecked/(float)deltat;
const float SOCrate = (float)len/(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 percentage = (RANDOM == 0) ? \
const float ETC = (RANDOM==0) ? ((float)(max - graphschecked)/(rate))/3600 : ((float)(RANDOM - len)/(SOCrate))/3600; 100*(float)graphnumber/(float)max : 100*(float)len/(float)RANDOM;
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); 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); fflush(stderr);
} }
// Convert graph index `grpahnumber' to the lists parents, children, parentslen, childrenlen // Convert graph index `grpahnumber' to the lists parents, children,
// parentslen, childrenlen
graphnrtolists(n, graphnumber, parents, parentslen, children, childrenlen); graphnrtolists(n, graphnumber, parents, parentslen, children, childrenlen);
// Increase checked counter and prepare for next iteration // Increase checked counter and prepare for next iteration
graphschecked++; graphschecked++;
if(RANDOM) if (RANDOM)
graphnumber = rand_uint64(max); graphnumber = rand_uint64(max);
else else
graphnumber++; graphnumber++;
// EO Increase checked counter and prepare for next iteration // EO Increase checked counter and prepare for next iteration
// If enabled, compare against degree options // If enabled, compare against degree options
if(NOSINK || NOSOURCE) { if (NOSINK || NOSOURCE) {
// Get some degree properties // Get some degree properties
int minindegree = n; int minindegree = n;
int minoutdegree = n; int minoutdegree = n;
for(int node=0; node<n; node++) { for (int node = 0; node < n; node++) {
if(minindegree > parentslen[node]) if (minindegree > parentslen[node])
minindegree = parentslen[node]; minindegree = parentslen[node];
if(minoutdegree > childrenlen[node]) if (minoutdegree > childrenlen[node])
minoutdegree = childrenlen[node]; minoutdegree = childrenlen[node];
} }
if(NOSINK && minoutdegree == 0) continue; if (NOSINK && minoutdegree == 0) continue;
if(NOSOURCE && minindegree == 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) // Check whether this graph is canonical or not; continue if graph has not
if(!isdegreeordered(n, parents, parentslen, children, childrenlen)) continue; // 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 // Ignore graphs that are not connected
if(!isconnected(n, parents, parentslen, children, childrenlen)) continue; if (!isconnected(n, parents, parentslen, children, childrenlen))
continue;
// Find cycles // Find cycles
const int num_cycles = find_cycles(cycles, cyclescnt, n, children, childrenlen); const int num_cycles = find_cycles(cycles, cyclescnt, n, children, \
childrenlen);
// If enabled, ignore DAGs // If enabled, ignore DAGs
if(NONDAGONLY && num_cycles > 0) continue; if (NONDAGONLY && num_cycles > 0) continue;
// Test the Siblings-On-Cycles property // Test the Siblings-On-Cycles property
// A DAG is trivially a SOC // A DAG is trivially a SOC
if(num_cycles > 0 && !gissoc(n, parents, parentslen, children, childrenlen, cycles, cyclescnt)) continue; if (num_cycles > 0 && !gissoc(n, parents, parentslen, children, \
childrenlen, cycles, cyclescnt)) continue;
// We have found a SOC // We have found a SOC
len++; len++;
if(GRAPHVIZ) if (GRAPHVIZ)
dumpgv(n, graphnumber, children, childrenlen); dumpgv(n, graphnumber, children, childrenlen);
else else
dumpgraph(n, children, childrenlen); dumpgraph(n, children, childrenlen);
} }
const int deltat = time(NULL) - t0; const int deltat = time(NULL) - t0;
const float rate = (deltat == 0) ? graphschecked : (float)graphschecked/(float)deltat; const float rate = (deltat == 0) ? graphschecked : \
(float)graphschecked/(float)deltat;
const float SOCrate = (deltat == 0) ? len : (float)len/(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, "\r%6.2f%% %*llu/%llu (%i SOCs found, %d seconds, %*.2f graphs/s, %*.2f SOCs/s, %*.2fh estimated time left)", \
fprintf(stderr, "\nFound %d SOCs in %d seconds\n", len,deltat); 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 manually allocated memory
free(cycles); free(cycles);
return 0; return 0;

26
src/SOCgraphviz.c
View File

@ -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 // SPDX-License-Identifier: GPL-3.0-or-later
@ -6,7 +7,7 @@
int main(int argc, char *argv[]) { int main(int argc, char *argv[]) {
// Parse command-line arguments // Parse command-line arguments
if(argc != 2) { if (argc != 2) {
fprintf(stderr, "Usage: %s <filename>\n", argv[0]); fprintf(stderr, "Usage: %s <filename>\n", argv[0]);
fprintf(stderr, " <filename> File name with adjacency matrices of simple directed graphs\n"); fprintf(stderr, " <filename> File name with adjacency matrices of simple directed graphs\n");
fprintf(stderr, "\n"); fprintf(stderr, "\n");
@ -22,7 +23,7 @@ int main(int argc, char *argv[]) {
// Open file // Open file
FILE *fp; FILE *fp;
fp = fopen(argv[1], "r"); fp = fopen(argv[1], "r");
if(fp == NULL){ if (fp == NULL) {
fprintf(stderr, "ERROR: Could not read graphs from file `%s'.\n", argv[1]); fprintf(stderr, "ERROR: Could not read graphs from file `%s'.\n", argv[1]);
return -1; return -1;
} }
@ -32,29 +33,28 @@ int main(int argc, char *argv[]) {
int n = -1; int n = -1;
int row = 0; int row = 0;
int col = 0; int col = 0;
char ch;
printf("strict digraph \"File_%s_line_%d\" {", argv[1], line); printf("strict digraph \"File_%s_line_%d\" {", argv[1], line);
while(!feof(fp)) { while (!feof(fp)) {
ch = fgetc(fp); char ch = fgetc(fp);
switch(ch) { switch (ch) {
case '{': case '{':
if(row==0 && col==0 && headprinted==0) { if (row == 0 && col == 0 && headprinted == 0) {
headprinted = 1; headprinted = 1;
printf("strict digraph \"File_%s_line_%d\" {", argv[1], line); printf("strict digraph \"File_%s_line_%d\" {", argv[1], line);
} }
break; break;
case '0': case '0':
case '1': case '1':
if(ch=='1') if (ch == '1')
printf("G%dN%d->G%dN%d;", line,row,line,col); printf("G%dN%d->G%dN%d;", line, row, line, col);
col++; col++;
break; break;
case '}': case '}':
if(n==-1) if (n == -1)
n = col; n = col;
if(col==n) { if (col == n) {
row++; row++;
col=0; col = 0;
} }
break; break;
case '\n': case '\n':

3
tools/removeiso.py
View File

@ -1,6 +1,7 @@
#!/usr/bin/python3 #!/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 SPDX-License-Identifier: GPL-3.0-or-later
""" """