// 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;
}