/*
	maximum-bipartite-matching.c

	To compile with GCC, execute:
	gcc -Wall -std=c99 -pedantic -o match maximum-bipartite-matching.c

	To run, execute: ./match graph.txt

	"Everything" has been extracted from:
	https://github.com/rhyscitlema/algorithms

	---------------------------------------------------

	To convert tabs to spaces use the expand command.
	See https://www.rhyscitlema.com/about/codingstyle.html

	Provided by Rhyscitlema
	https://www.rhyscitlema.com/algorithms/maximum-bipartite-matching

	USE AT YOUR OWN RISK!
*/
#include <stdio.h>
#include <string.h>
#include <stddef.h>
#include <stdbool.h>
#include <assert.h>
#include <malloc.h>


typedef struct { char* data; long size; } Array;
Array file_to_array (const char* filename, Array old);

typedef long cost; // must match with load_edge_list()
typedef struct { int u, v; cost c; } Edge, *EdgeP;
typedef struct { int    v; cost c; } Adge, *Adja, **AdjaP;
typedef const Edge* const_EdgeP;
typedef const Adja* const_AdjaP;

AdjaP load_adja_list (const char* str); // parses and returns adjacency list
char* print_edge_list (char* out, const_EdgeP edge); // returns end of string
EdgeP maximum_matching_unweighted (const_AdjaP adja, bool skipPart1);


int main (int argc, char** argv)
{
	if(argc!=2)
		printf("Usage: <program> graph.txt\n");
	else
	{
		Array content = {0};
		content = file_to_array(argv[1], content);
		if(content.size<=0)
			printf("Error: failed to open file '%s'.\n", argv[1]);
		else
		{
			AdjaP adjaList = load_adja_list(content.data);
			EdgeP edgeList = maximum_matching_unweighted(adjaList, true);

			char str[10000];
			print_edge_list(str, edgeList);
			puts(str);

			free(adjaList);
			free(edgeList);
		}
		free(content.data);
	}
	return 0;
}


Array file_to_array (const char* filename, Array old)
{
	Array out={0};

	FILE* file = fopen (filename, "rb");
	if(file==NULL) out.size = -1;
	else
	{
		fseek (file, 0, SEEK_END);
		out.size = ftell (file);
		fseek (file, 0, SEEK_SET);

		if(out.size == 0x7FFFFFFF) out.size = -2;
		else if(old.size != -1)
		{
			out.data = (char*) realloc (old.data, out.size+1);

			if(fread(out.data, 1, out.size, file) != out.size)
			{
				free(out.data);
				out.data = NULL;
			}
			else out.data[out.size] = 0;
		}
		fclose(file);
	}
	return out;
}



AdjaP edge_list_to_adja_list (const_EdgeP edge)
{
	if(!edge) return NULL;
	int V = edge[0].u;
	int E = edge[0].v;
	bool directed = (bool)edge[0].c;

	Adja *next, *prev;
	Adja dnex, dpre;
	int i, u, v;

	u = 1+V;     // for adjacency list pointers.
	v = 1+V+E*2; // for the data at list pointers.
	if(directed) { u+=V; v+=V; }

	next = (AdjaP)malloc(u*sizeof(Adja) + v*sizeof(Adge));

	dnex = (Adja)(next+u);
	prev = next; if(directed) prev += V;
	dpre = dnex; if(directed) dpre += V+E;

	*(int*)next = V;
	dnex->v = E;
	dnex->c = directed;

	// get ready to count adjacent vertices
	for(i=1; i<=V; i++)
	{
		dnex[i].v=0;
		dpre[i].v=0;
	}
	for(i=1; i<=E; i++)
	{
		u = edge[i].u;
		v = edge[i].v;
		dnex[u].v++;    // count adjacent vertices
		dpre[v].v++;
	}
	u=1; v=1;
	for(i=1; i<=V; i++)
	{
		next[i] = dnex+u;   // assign pointers of
		prev[i] = dpre+v;   // adjacency list
		u += 1+dnex[i].v;
		v += 1+dpre[i].v;
	}
	for(i=1; i<=V; i++)
	{
		next[i][0].v=0;
		prev[i][0].v=0;
	}
	for(i=1; i<=E; i++)
	{
		cost c;
		int j;
		u = edge[i].u;
		v = edge[i].v;
		c = edge[i].c;

		j = ++next[u][0].v;
		next[u][j].v = v;
		next[u][j].c = c;

		j = ++prev[v][0].v;
		prev[v][j].v = u;
		prev[v][j].c = c;
	}
	return next;
}

#define AdjaU(adja,u) (adja[u]) // get list of vertexes adjacent to u
#define AdjaV(adja) (*(int*)(adja))         // get number of vertexes
#define AdjaE(adja) (AdjaU(adja,1)-1)->v    // get number of edges
#define AdjaD(adja) (AdjaU(adja,1)-1)->c    // get graph-is-directed
#define AdjaPrev(adja) (adja + (AdjaD(adja) ? AdjaV(adja) : 0))



const char* strNext1 (const char* str)
{
	while(true) // not a loop
	{
		if(!str || !*str) break;
		if(*str != '#')
		{
			str++;
			if(!*str) break;
			if(*str != '#') break;
		}
		str++;
		if(!*str) break;
		if(*str != '{')
		{
			// skip till '\n' is found, excluding it
			while(*str && *str != '\n') str++;
			break;
		}
		// else: skip till "}#" is found, including it
		char a, b;
		int level=1;
		str++;
		b = *str;
		while(true)
		{
			str++;
			if(!*str) break;
			a = b;
			b = *str;
			if(a=='#' && b=='{') level++;
			if(a=='}' && b=='#')
			{
				if(--level==0)
				{ str++; break; }
			}
		}
		break;
	}
	return str;
}

static inline bool isSpace(char c)
{ return (c==' ' || c=='\t' || c=='\r' || c=='\n'); }

const char* nextWord (const char* str)
{
	while(!isSpace(*str) && *str) str = strNext1(str);
	while( isSpace(*str)) str = strNext1(str);
	return str;
}



EdgeP load_edge_list (const char* str)
{
	if(!str) return NULL;
	int i, V, E;
	bool directed;
	bool weighted;

	if(*str=='#') str = strNext1(str);
	while(isSpace(*str)) str = strNext1(str);

	if(!*str || sscanf(str, "%d", &V)!=1) return NULL; str=nextWord(str);
	if(!*str || sscanf(str, "%d", &E)!=1) return NULL; str=nextWord(str);
	if(!*str || sscanf(str, "%d", &i)!=1) return NULL; str=nextWord(str); directed=i;
	if(!*str || sscanf(str, "%d", &i)!=1) return NULL; str=nextWord(str); weighted=i;
	if(!*str || sscanf(str, "%d", &i)!=1) return NULL; str=nextWord(str);
	if(V<0 || E<0) return NULL;

	EdgeP edge = (EdgeP)malloc((1+E)*sizeof(Edge));
	edge[0].u = V;  // V = number of vertexes
	edge[0].v = E;  // E = number of edges
	edge[0].c = directed;

	if(i==0) // if an edge_list
	{
		for(i=1; i<=E; i++)
		{
			sscanf(str, "%d", &edge[i].u); str=nextWord(str); // get source vertex u
			sscanf(str, "%d", &edge[i].v); str=nextWord(str); // get sink vertex v
			if(weighted)
			{
				sscanf(str, "%ld", &edge[i].c); // get edge cost c
				str=nextWord(str);
			}
			else edge[i].c = 1; // else mark as unweighted
		}
	}
	else if(i==1) // else an adja_list of source-to-sink vertexes
	{
		int t, u, m;
		for(E=0, t=1; t<=V; t++)
		{
			sscanf(str, "%d", &u); str=nextWord(str); // get source vertex u
			sscanf(str, "%d", &m); str=nextWord(str); // get number of adjacent
			for(i=1; i<=m; i++)
			{
				cost c=1;
				int v;
				sscanf(str, "%d", &v); // get sink vertex v
				str=nextWord(str);

				if(weighted)
				{
					sscanf(str, "%ld", &edge[i].c); // get edge cost c
					str=nextWord(str);
				}
				else edge[i].c = 1; // else mark as unweighted

				if(!directed && u>v) continue; // avoid redundant edges
				E++;
				edge[E].u = u;  // record loaded edge
				edge[E].v = v;
				edge[E].c = c;
			}
		}
	}
	else assert(false && "Graph neither edge list nor adja list");
	return edge;
}

AdjaP load_adja_list (const char* str)
{
	EdgeP edge = load_edge_list (str);
	AdjaP adja = edge_list_to_adja_list (edge);
	free(edge);
	return adja;
}



#ifndef HEADER
#define graph_header "\n%d %d %d %d %d\n\n"
#else
#define graph_header
	"#{\n"
	"\tThis is a block comment: #{...}#\n"
	"\tA single-line comment has form: #...\n"
	"\tVertex numbering starts from 1 not 0.\n"
	"\tNote the content structure as described.\n"
	"}#\n"
	"\n"
	"%-5d # number of vertexes\n"
	"%-5d # number of edges\n"
	"\n"
	"%d     # 1 if graph is directed, 0 otherwise\n"
	"%d     # 1 if graph is weighted, 0 otherwise\n"
	"\n"
	"%d     # 0 if graph is given as an edge list, or,\n"
	"      # 1 if graph is given as an adjacency list\n"
	"      #      of source-to-sink vertexes.\n"
	"\n";
#endif


char* print_edge_list (char* out, const_EdgeP edge)
{
	if(!out) return out;
	if(!edge) { *out=0; return out; }

	int V = edge[0].u;
	int E = edge[0].v;
	bool directed = (bool)edge[0].c;
	bool weighted = 0;
	int i;

	// check if any edge cost is not 1
	for(i=1; i<=E; i++) if(edge[i].c!=1) { weighted=1; break; }

	out += sprintf(out, graph_header, V, E, directed, weighted, 0);

	for(i=1; i<=E; i++)
	{
		out += sprintf(out, "%d %d", edge[i].u, edge[i].v);
		if(weighted) out += sprintf(out, " %ld", edge[i].c);
		*out++ = '\n';
	}
	*out=0;
	return out;
}



/**
	First initialise heap.size to = 0
	and initialise heap.data to =
	(void**) malloc (heap.maxSize * sizeof(void*));

	The functions will only change these two.
	Set heap.node_compare to the callback function.
	heap.arg is only passed to heap.node_compare().

	** Indexed heap:
	This is only necessary if heap_update() or heap_remove() will be used.
	If heap is indexed then every "void* node" must actually be an int* node
	pointing to FREE int* memory to only be used by these heap functions.
*/
typedef struct _Heap
{   void** data;
	int size;
	int maxSize;
	bool indexed;
	const void* arg;
	int (*node_compare) (const void* a, const void* b, const void* arg);
} Heap;



#define INDEX_UPDATE if(heap->indexed) *(int*)data[i] = i;

static bool update_downward (Heap* heap, void* node, int i)
{
	void** data = heap->data;
	while(i>1)
	{
		if(heap->node_compare(node, data[i/2], heap->arg) < 0)
		{
			data[i] = data[i/2];
			INDEX_UPDATE
			i = i/2;
		}
		else break;
	}
	if(data[i]==node) return 0; // if no update
	data[i] = node;
	INDEX_UPDATE
	return 1;
}

static bool update_upward (Heap* heap, void* node, int i)
{
	void** data = heap->data;
	int size = heap->size;
	while(1)
	{
		if(i*2 > size) break;

		if(i*2 == size
		|| heap->node_compare(data[i*2], data[i*2+1], heap->arg)<0)
		{
			if(heap->node_compare(data[i*2], node, heap->arg)<=0)
			{
				// shift left child
				data[i] = data[i*2];
				INDEX_UPDATE
				i = i*2;
				continue;
			}
		}
		else
		{
			if(heap->node_compare(data[i*2+1], node, heap->arg)<=0)
			{
				// shift right child
				data[i] = data[i*2+1];
				INDEX_UPDATE
				i = i*2+1;
				continue;
			}
		}
		break;
	}
	if(data[i]==node) return 0; // if no update
	data[i] = node;
	INDEX_UPDATE
	return 1;
}

void heap_push (Heap* heap, void* node)
{
	assert(heap!=NULL);
	if(!heap) return;
	int i = ++heap->size;
	heap->data[i] = node;
	if(heap->indexed) *(int*)node = i;
	update_downward(heap, node, i);
}

void* heap_pop (Heap* heap)
{
	assert(heap!=NULL);
	if(!heap) return NULL;
	if(heap->size <= 0) return NULL;
	void* root = heap->data[1];
	if(heap->indexed) *(int*)root = -1;
	int size = heap->size--;
	void* node = heap->data[size];
	update_upward(heap, node, 1);
	return root;
}

void heap_remove (Heap* heap, void* node)
{
	assert(heap && node && heap->indexed);
	if(!(heap && node && heap->indexed)) return;
	int i = *(int*)node;
	if(i<0) return;
	*(int*)node = -1;
	int size = heap->size--;
	node = heap->data[size];
	update_upward(heap, node, i);
}

void heap_update (Heap* heap, void* node)
{
	assert(heap && node && heap->indexed);
	if(!(heap && node && heap->indexed)) return;
	int i = *(int*)node;
	if(i<0) return;
	if(!update_downward(heap, node, i))
		update_upward(heap, node, i);
}



typedef struct {
	int i;  // the free-memory needed by heap.h
	int c;  // count of free adjacent vertexes
	// i must be declared first so to have:
	// vertex = (mmug*)heap_node - mmug_array;
} mmug;

/* compare count of free adjacent vertexes */
static int mmug_heap_node_compare (const void* a, const void* b, const void* arg)
{ return ( ((const mmug*)a)->c - ((const mmug*)b)->c ); }


EdgeP maximum_matching_unweighted (const_AdjaP adja, bool skipPart1)
{
	if(!adja) return NULL;
	int i, m, u, v;
	int V = AdjaV(adja); // get V = number of vertexes

	const_AdjaP graph, prev = AdjaPrev(adja);
	// adja = adjacency list of source-to-sink vertexes
	// prev = adjacency list of sink-from-source vertexes
	// if adja==prev then graph is an undirected graph

	// get m = total memory needed
	m = (1+V)*(sizeof(int)+sizeof(mmug)+sizeof(void*));
	int* match = (int*)malloc(m);
	mmug* count = (mmug*)(match+1+V);
	void** hdata = (void**)(count+1+V);

	if(match==NULL) return NULL;    // if failed to allocate memory
	for(u=1; u<=V; u++) match[u]=0; // else initialise main array


	//*****************************************************************
	//*** Part 1 of algorithm : is often enough to solve optimally! ***
	//*****************************************************************
	if(!skipPart1) {

	// initialise the heap data structure
	Heap _heap = { hdata, 0, V, true, NULL, mmug_heap_node_compare };
	Heap *heap = &_heap;

	// initialise the heap with all vertexes
	for(u=1; u<=V; u++)
	{
		// get m = total number of adjacent vertexes
		m = AdjaU(adja,u)[0].v;
		if(adja!=prev) // if a directed graph
			m += AdjaU(prev,u)[0].v;
		count[u].c = m;
		if(m) // if at least one adjacent vertex
			heap_push(heap, &count[u].i);
	}

	while(true)
	{
		void* hnode = heap_pop(heap);
		if(hnode==NULL) break; // if heap is empty then quit

		// get u = vertex with minimum count of free adjacent vertexes
		u = (mmug*)hnode - count;

		if(count[u].c==0) continue; // if no free adjacent vertex then skip,
		count[u].c=0;               // else mark u as removed from the heap.

		int t=0;
		graph = adja;               // start with 'adja' (source-to-sink) list
		while(true)
		{
			Adja U = AdjaU(graph,u);
			m = U[0].v;             // get number of adjacent vertexes
			for(i=1; i<=m; i++)     // for every adjacent vertex
			{
				v = U[i].v;         // get adjacent vertex
				if(count[v].c>0)    // if is inside heap
				{
					count[v].c--;   // decrement count of free adjacent vertexes
					heap_update(heap, &count[v].i); // update position in heap

					// get free vertex with minimum count of free adjacent vertexes
					if(t==0 || count[t].c > count[v].c) t=v;
				}
			}
			if(graph==prev) break;
			else graph = prev;      // switch to 'prev' (sink-from-source) list
		}
		v=t;    // get v = new free vertex

		match[u] = v;   // match free vertex
		match[v] = u;

		count[v].c = 0; // remove v from heap
		heap_remove(heap, &count[v].i);

		// now tell all adjacent to v that v is monopolised.
		u=v; // set u=v so to use exact same code as before!

		graph = adja;               // start with 'adja' (source-to-sink) list
		while(true)
		{
			Adja U = AdjaU(graph,u);
			m = U[0].v;             // get number of adjacent vertexes
			for(i=1; i<=m; i++)     // for every adjacent vertex
			{
				v = U[i].v;         // get adjacent vertex
				if(count[v].c>0)    // if is inside heap
				{
					count[v].c--;   // decrement count of free adjacent vertexes
					heap_update(heap, &count[v].i); // update position in heap
				}
			}
			if(graph==prev) break;
			else graph = prev;      // switch to 'prev' (sink-from-source) list
		}
	}
	} // end of if(!skipPart1)


	//****************************************************************
	//*** Part 2 of algorithm : may be faster when executed alone! ***
	//****************************************************************

	int E; // number of edges of the result graph
	while(true)
	{
		int beg=0, end=0;
		int* queue = match+(1+V)*1; // queue used by the breadth-first search.
		int* claim = match+(1+V)*2; // claim[v] is the claimer of the claimed v.
		int* vfree = match+(1+V)*3; // vfree[v] is the free vertex to start claiming

		// initialisations, push all unmatched vertexes to queue
		for(u=1; u<=V; u++)
		{
			if(match[u]) // if u is a matched vertex
			     { claim[u] = vfree[u] = 0; }
			else { claim[u] = vfree[u] = queue[end++] = u; }
		}

		// get E = count of matched vertex pairs
		E = (V-end)/2;
		//printf("-------------------------------------- E=%d\n", E);

		for( ; beg!=end; beg++)     // do the queue-based backtracking-search
		{
			u = queue[beg];
			v = vfree[u];           // get v = free vertex which u is derived from
			assert(v>0);
			if(match[v]) continue;  // check whether this v is still unmatched

			graph = adja;           // start with 'adja' (source-to-sink) list
			while(true)
			{
				Adja U = AdjaU(graph,u);
				m = U[0].v;  // get number of adjacent vertexes

				for(i=1; i<=m; i++)     // for every adjacent vertex
				{
					v = U[i].v;         // get adjacent vertex
					int t = vfree[v];   // get t = free vertex which v is derived from

					// check whether a new re-matching has been uncovered
					if(t && !match[t] && t != vfree[u])
					{
						claim[v] = u;       // ensure v is never claimed again
						int tv = match[v];  // prepare for 2nd re-matching
						while(true)
						{
							claim[u] = u;   // ensure u is never claimed again
							t = match[u];   // prepare for next iteration
							match[u] = v;
							match[v] = u;   // uncover the augmenting path
							if(t==0){
								if(tv==0) break; // if no further re-matching
								t=tv; tv=0;      // prevent any further re-matching
							}
							v = t;          // get v = the next claimed vertex
							u = claim[t];   // get u = the claimer to this v
						}
						E=-1; // notify that a re-matching process was performed
						break;
					}
					if(!claim[v])   // if not yet claimed
					{
						claim[v] = u;       // Mark u as claiming v
						t = match[v];       // Get who monopolises v
						vfree[t] = vfree[u];
						queue[end++] = t;   // Push it to queue
					}
				}
				if(i<=m) break; // if re-matching was performed then break

				if(graph==prev) break;
				else graph = prev;  // switch to 'prev' (sink-from-source) list
			}
		}
		if(E!=-1) break; // no re-matching performed means solution is optimal
	}

	// allocate memory for the result graph
	EdgeP edge = (EdgeP) malloc ((1+E)*sizeof(Edge));

	// build the edge-list of the result graph
	for(E=0, u=1; u<=V; u++)
	{
		v = match[u];
		if(v<u) continue; // avoid repetitions
		++E;
		edge[E].u = u;
		edge[E].v = v;
		edge[E].c = 1; // set weight to = 1
	}
	edge[0].u = V;
	edge[0].v = E;
	edge[0].c = 0; // mark as undirected

	// free allocated memory
	free(match);
	return edge;
}