SOC-Observation-Code/SOCadmissible.c

406 lines
14 KiB
C
Raw Normal View History

2023-05-04 18:10:31 +02:00
// 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 (verification rate %*.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 (verification rate %*.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-START+1; // Number of graphs to be verified
printf("startidx=%d, tot=%d\n", startidx, tot);
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;
}