/*
* OPTIMIZE - Make the cross-identification between CD and PPMX
* Made by Daniel E. Severin (CONICET) in 2015-2017
*/

#include "common.h"
#include <gurobi_c++.h>
#include <time.h>
#ifdef LINUX
#include <unistd.h>
#include <sys/times.h>
#endif


#define PROB0 1e-10
#define PROBSLG 1e-4
#define DISTTOLDPL 80.0/3600.0
#define DSEPPARAM 34.9/3600.0
#define SIGMADSEPPARAM 13.65/3600.0
#define SIGMAMAGPARAM 0.915
#define MAXSIZEPA 1000
#define MAXSIZEACOMP 100
#define MAXSIZEBCOMP 200
#define MAXVARIABLES 1000
//#define GENERATE_CROSS23 /* just for comparing with previous cross-identification */

struct CDstar_struct *CDstar;
struct Pa_struct *Pa; /* sets P_a for each CD star */
int CDstars; /* number of CD stars */

struct PPMXstar_struct *PPMXstar;
int PPMXstars; /* number of PPMX stars */

#define MAXEDGES 150000
int edge_A[MAXEDGES], edge_B[MAXEDGES]; /* "edges" of the initial matching */
int edges; /* number of "edges" */
int *degree_A, *degree_B; /* "degrees of vertices" of the initial matching */
int *A1; /* <-- A1[v] = first element of component having v as an element */
int *A2; /* <-- A2[v] = cardinal of component having v as an element */
int *A3; /* <-- A3[v] = next element in component after v, -1 = if v is the last one */
int *A_comp; /* active stars of first catalog in a component (A') */
int *B_comp; /* active stars of second catalog in a component (B') */
int hardest_variables, hardest_constraints; /* Hardest instance: variables and constraints */
double hardest_nodes, hardest_time; /* Hardest instance: nodes of B&B tree and CPU time */

/*
* Tclock - obtain a timestamp (in seconds)
*/
double Tclock() {
#ifdef LINUX
	/* measure user-time: use it on single-threaded system in Linux (is more accurate) */
	struct tms buf;
	times(&buf);
	return ((double)buf.tms_utime) / (double)sysconf(_SC_CLK_TCK);
#else
	/* measure standard wall-clock: use it on Windows
	* (Linux clock() function overflows every 72 minutes) */
	return ((double)clock()) / (double)CLOCKS_PER_SEC;
#endif
}

/*
* optimize - perform the optimization over a component
*/
void optimize(int v)
{
	try {
		GRBEnv env = GRBEnv();
		GRBModel model = GRBModel(env);

		/* obtain sets A' and B' (sets of stars restricted to the component) */
		int Astars = 0, Bstars = 0;
		int w = v;
		for (int i = 0; i < A2[v]; i++) {
			if (w < CDstars) {
				if (Astars == MAXSIZEACOMP) bye("Maximum size of Acomp reached!\n");
				A_comp[Astars++] = w;
			}
			else {
				if (Bstars == MAXSIZEBCOMP) bye("Maximum size of Bcomp reached!\n");
				B_comp[Bstars++] = w - CDstars;
			}
			w = A3[w];
		}

		/* set some parameters */
		model.getEnv().resetParams();
		model.getEnv().set(GRB_DoubleParam_NodefileStart, 0.5);
		model.getEnv().set(GRB_DoubleParam_TimeLimit, 600.0); /* limit: 10 minutes of CPU time */
		model.getEnv().set(GRB_DoubleParam_MIPGap, 0.0);
		model.getEnv().set(GRB_DoubleParam_MIPGapAbs, 0.0);
		model.getEnv().set(GRB_DoubleParam_IntFeasTol, 0.000001);
		model.getEnv().set(GRB_IntParam_Threads, 1);
		model.getEnv().set(GRB_IntParam_Seed, 1);
		model.getEnv().set(GRB_IntParam_OutputFlag, 0);

		/* create binary variables x(i,j), i in A and j in P_a */
		GRBVar *vars = new GRBVar[MAXVARIABLES];
		int variables = 0;
		char name[128];
		for (int i = 0; i < Astars; i++) {
			int a = A_comp[i];
			for (int j = 0; j < Pa[a].cardinal; j++) {
				if (variables == MAXVARIABLES) bye("Maximum number of variables reached!\n");
				sprintf(name, "x(%d,%d,%d)", a, Pa[a].b1[j], Pa[a].b2[j]);
				vars[variables++] = model.addVar(0.0, 1.0, 0.0, GRB_BINARY, name);
			}
		}
		model.update();

		/* generate objective function */
		GRBLinExpr fobj = 0.0;
		int index = 0;
		for (int i = 0; i < Astars; i++) {
			int a = A_comp[i];
			for (int j = 0; j < Pa[a].cardinal; j++) fobj += Pa[a].weight[j] * vars[index++];
		}
		model.setObjective(fobj, GRB_MINIMIZE);

		/* generate constraints (1) */
		int constraints = 0;
		index = 0;
		for (int i = 0; i < Astars; i++) {
			int a = A_comp[i];
			GRBLinExpr restr = 0.0;
			for (int j = 0; j < Pa[a].cardinal; j++) restr += vars[index++];
			model.addConstr(restr == 1);
			constraints++;
		}

		/* generate constraints (2) */
		for (int k = 0; k < Bstars; k++) {
			int b = B_comp[k];
			GRBLinExpr restr = 0.0;
			int varadded = 0;
			index = 0;
			for (int i = 0; i < Astars; i++) {
				int a = A_comp[i];
				for (int j = 0; j < Pa[a].cardinal; j++) {
					if (Pa[a].b1[j] == b || Pa[a].b2[j] == b) {
						restr += vars[index];
						varadded++;
					}
					index++;
				}
			}
			if (varadded >= 2) {
				model.addConstr(restr <= 1);
				constraints++;
			}
		}

		/* perform MIP optimization */
		model.update();
		//model.write("model.lp");
		model.optimize();

		int status = model.get(GRB_IntAttr_Status);
		if (status != GRB_OPTIMAL) {
			switch (status) {
			case GRB_INF_OR_UNBD:
			case GRB_INFEASIBLE: bye("Infeasible :(");
			case GRB_TIME_LIMIT: bye("Time limit reached :(");
			default: bye("Unexpected error :(");
			}
		}

		/* obtain results from optimization */
		double nodes = model.get(GRB_DoubleAttr_NodeCount);
		if (nodes > 0.5) printf("%.0lf nodes of B&B tree were explored\n", nodes);
		double runtime = model.get(GRB_DoubleAttr_Runtime);
		if (runtime > hardest_time) {
			hardest_variables = variables;
			hardest_constraints = constraints;
			hardest_nodes = nodes;
			hardest_time = runtime;
		}

		/* make assignment */
		index = 0;
		for (int i = 0; i < Astars; i++) {
			int a = A_comp[i];
			for (int j = 0; j < Pa[a].cardinal; j++) {
				if (vars[index].get(GRB_DoubleAttr_X) > 0.5) Pa[a].chosen = j;
				index++;
			}
		}
		return;
	}
	catch (GRBException e) {
		printf("Error code = %d\n", e.getErrorCode());
	}
	catch (...) {
		printf("An exception has occured during optimization\n");
	}
	bye("Internal error!\n");
}

/*
 * write_assignment - write one register of the new cd.txt
 */
void write_assignment(FILE *stream, int a, int b1, int b2)
{
	char line[128];
	sprintf(&line[0], "CD%s", CDstar[a].reg);
	sprintf(&line[30], "%c%c", CDstar[a].color ? 'C' : ' ', CDstar[a].dpl ? 'D' : ' ');
	if (b1 == -1) {
		/* an unassigned star */
		sprintf(&line[32], "0                                      ");
	}
	else {
		if (b2 == -1) {
			/* a single assigned star */
			sprintf(&line[32], "1PPMX%s                   ", PPMXstar[b1].ppmx_name);
		}
		else {
			/* a double assigned star */
			sprintf(&line[32], "2PPMX%sPPMX%s", PPMXstar[b1].ppmx_name, PPMXstar[b2].ppmx_name);
		}
	}
	fprintf(stream, "%s\n", line);
}
#ifdef GENERATE_CROSS23
void write_assignment(FILE *stream, int a, int b1, int b2)
{
	char line[128];
	if (CDstar[a].decl_ref != -23) return;
	sprintf(&line[0], "CD%s", CDstar[a].reg);
	if (b1 == -1) return;
	double x1, y1, z1, x2, y2, z2;
	sph2rec(CDstar[a].RA1875, CDstar[a].DE1875, &x1, &y1, &z1);
	if (b2 == -1) {
		/* a single assigned star */
		sph2rec(PPMXstar[b1].RA1875, PPMXstar[b1].DE1875, &x2, &y2, &z2);
		float dist = calcdist(calcscalar(x1, y1, z1, x2, y2, z2));
		sprintf(&line[11], "%c %6.2fPPMX%s", CDstar[a].dpl ? 'U' : 'S', dist * 3600.0, PPMXstar[b1].ppmx_name);
	}
	else {
		/* a double assigned star */
		sph2rec(PPMXstar[b1].RA1875, PPMXstar[b1].DE1875, &x2, &y2, &z2);
		float dist = calcdist(calcscalar(x1, y1, z1, x2, y2, z2));
		sprintf(&line[11], "A %6.2fPPMX%s", dist * 3600.0, PPMXstar[b1].ppmx_name);
		fprintf(stream, "%s\n", line);
		sph2rec(PPMXstar[b2].RA1875, PPMXstar[b2].DE1875, &x2, &y2, &z2);
		dist = calcdist(calcscalar(x1, y1, z1, x2, y2, z2));
		sprintf(&line[11], "B %6.2fPPMX%s", dist * 3600.0, PPMXstar[b2].ppmx_name);
	}
	fprintf(stream, "%s\n", line);
}
#endif

/*
 * main
 */
int main(int argc, char** argv)
{
	FILE *stream;
	char buffer[1024];
	char cell[256];

	printf("OPTIMIZE - Make the cross-identification between CD and PPMX.\n");
	printf("Made by Daniel E. Severin (CONICET) in 2015-2017.\n");

	/* read CD data */
	stream = fopen("cd_data.bin", "rb");
	if (stream == NULL) {
		perror("Error reading bin file");
		exit(1);
	}
	fread(&CDstars, sizeof(int), 1, stream);
	CDstar = (struct CDstar_struct *)calloc(CDstars, sizeof(struct CDstar_struct));
	fread(CDstar, sizeof(struct CDstar_struct), CDstars, stream);
	fclose(stream);

	/* read PPMX data */
	stream = fopen("ppmx_data.bin", "rb");
	if (stream == NULL) {
		perror("Error reading bin file");
		exit(1);
	}
	fread(&PPMXstars, sizeof(int), 1, stream);
	PPMXstar = (struct PPMXstar_struct *)calloc(PPMXstars, sizeof(struct PPMXstar_struct));
	fread(PPMXstar, sizeof(struct PPMXstar_struct), PPMXstars, stream);
	fclose(stream);

	/* read initial cross-identification */
	stream = fopen("tables\\initial_cross.txt", "rt");
	if (stream == NULL) {
		perror("Error reading initial cross file");
		exit(1);
	}

	degree_A = (int*)calloc(CDstars, sizeof(int));
	degree_B = (int*)calloc(PPMXstars, sizeof(int));
	for (int i = 0; i < CDstars; i++) degree_A[i] = 0;
	for (int i = 0; i < PPMXstars; i++) degree_B[i] = 0;

	edges = 0;
	while (fgets(buffer, 1023, stream) != NULL) {
		int CDid, PPMXid;
		sscanf(buffer, "%d,%d\n", &CDid, &PPMXid);
		if (CDid < 0 || CDid >= CDstars || PPMXid < 0 || PPMXid >= PPMXstars) {
			printf("Error reading edge: (%d,%d)\n", CDid, PPMXid);
			exit(1);
		}
		if (CDstar[CDid].deleted > 0) bye("Internal error!\n");
		if (PPMXstar[PPMXid].vmag > 13.5) continue; /* omit faint stars of PPMX */
		/* if the distance between both stars is more than 2.2 arcmin, warn about it
		 * (there are few cases due to the proper motion of PPMX stars) */
		double x1, y1, z1, x2, y2, z2;
		sph2rec(CDstar[CDid].RA1875, CDstar[CDid].DE1875, &x1, &y1, &z1);
		sph2rec(PPMXstar[PPMXid].RA1875, PPMXstar[PPMXid].DE1875, &x2, &y2, &z2);
		double dist = calcdist(calcscalar(x1, y1, z1, x2, y2, z2));
		if (dist > DISTTOL*1.1) printf("Note: edge (%d,%d) has angular distance = %.4f arcmin\n", CDid, PPMXid, dist*60.0);
		if (edges == MAXEDGES) bye("Maximum number of edges reached!\n");
		edge_A[edges] = CDid;
		edge_B[edges] = PPMXid;
		edges++;
		degree_A[CDid]++;
		degree_B[PPMXid]++;
	}
	fclose(stream);

	printf("Number of edges processed: %d\n", edges);

	double start_t = Tclock();
	/* compute effective cardinal of A and B */
	int card_A = 0, card_A2 = 0, max_degree_A = 0;
	for (int i = 0; i < CDstars; i++) {
		int d = degree_A[i];
		if (d > 0) {
			card_A++;
			if (CDstar[i].dpl) card_A2++;
			if (max_degree_A < d) max_degree_A = d;
		}
	}
	printf("|A| = %d,   |A_2| = %d,   max_degree = %d\n", card_A, card_A2, max_degree_A);
	int card_B = 0, max_degree_B = 0;
	for (int i = 0; i < PPMXstars; i++) {
		int d = degree_B[i];
		if (d > 0) {
			card_B++;
			if (max_degree_B < d) max_degree_B = d;
		}
	}
	printf("|B| = %d,   max_degree = %d\n", card_B, max_degree_B);

	/* ask for memory */
	Pa = (struct Pa_struct *)calloc(CDstars, sizeof(struct Pa_struct));
	struct Pa_struct *Pa_current = (struct Pa_struct *)calloc(1, sizeof(struct Pa_struct));
	Pa_current->b1 = (int*)calloc(MAXSIZEPA, sizeof(int));
	Pa_current->b2 = (int*)calloc(MAXSIZEPA, sizeof(int));
	Pa_current->weight = (long double*)calloc(MAXSIZEPA, sizeof(long double));

	/* generate sets P_a for each CD star */
	printf("Generating P_a sets...\n");
	int max_Pa_card = 0;
	int *neighbor = (int*)calloc(max_degree_A, sizeof(int));
	for (int i = 0; i < CDstars; i++) {
		int d = degree_A[i];
		Pa[i].chosen = -1;
		if (d == 0) {
			/* unassigned star */
			Pa[i].cardinal = 0;
		}
		else {
			/* obtain alpha*, delta and m */
			double delta1 = CDstar[i].DE1875;
			double alpha1 = CDstar[i].RA1875;
			double mag1 = CDstar[i].mag;
			
			/* obtain the set of "near" stars in B */
			int neighbors = 0;
			for (int e = 0; e < edges; e++) if (edge_A[e] == i) neighbor[neighbors++] = edge_B[e];
			if (neighbors != d) bye("Internal error!\n");

			if (d == 1) {
				/* If only one star is near (independently whether CD star is double or not), the weight is obvious */
				Pa_current->b1[0] = neighbor[0];
				Pa_current->b2[0] = -1;
				Pa_current->weight[0] = -log(1 - PROB0);
				Pa_current->b1[1] = -1;
				Pa_current->b2[1] = -1;
				Pa_current->weight[1] = -log(PROB0);
				Pa_current->cardinal = 2;
			}
			else {
				if (CDstar[i].dpl == false) {
					/* TREATMENT OF A SINGLE STAR */
					long double divisor = 0.0;
					for (int s = 0; s < d; s++) {
						int j = neighbor[s];
						double delta2 = PPMXstar[j].DE1875;
						double alpha2 = PPMXstar[j].RA1875;
						double mag2 = PPMXstar[j].mag;

						/* here, we just compute PDF1(0,0,0; a, b) */
						long double PDF1 = 1.0;
						PDF1 *= (long double)pdfn(alpha_difference(alpha1, alpha2, (delta1 + delta2) / 2.0), SIGMAALPHA);
						PDF1 *= (long double)pdfn(delta1 - delta2, SIGMADELTA);
						PDF1 *= (long double)pdfn(mag1 - mag2, SIGMAMAG);

						Pa_current->b1[s] = j;
						Pa_current->b2[s] = -1;
						Pa_current->weight[s] = PDF1;
						divisor += PDF1;
					}
					/* now, we normalize probabilities and compute weights */
					for (int s = 0; s < d; s++) {
						long double PDF1 = Pa_current->weight[s];
						long double probability = (PDF1 / divisor) * (1 - PROB0);
						Pa_current->weight[s] = -log(probability);
					}

					/* add the "emptyset" element */
					Pa_current->b1[d] = -1;
					Pa_current->b2[d] = -1;
					Pa_current->weight[d] = -log(PROB0);
					Pa_current->cardinal = d + 1;
				}
				else {
					/* TREATMENT OF A DOUBLE STAR */
					long double divisor = 0.0;
					int Pa_card = 0;
					for (int s = 0; s < d; s++) {
						int j = neighbor[s];
						double delta2 = PPMXstar[j].DE1875;
						double alpha2 = PPMXstar[j].RA1875;
						double mag2 = PPMXstar[j].mag;

						/* here, we just compute PDF1(0,0,0; a, b) */
						long double PDF1 = 1.0;
						PDF1 *= (long double)pdfn(alpha_difference(alpha1, alpha2, (delta1 + delta2) / 2.0), SIGMAALPHA);
						PDF1 *= (long double)pdfn(delta1 - delta2, SIGMADELTA);
						PDF1 *= (long double)pdfn(mag1 - mag2, SIGMAMAG);

						Pa_current->b1[s] = j;
						Pa_current->b2[s] = -1;
						Pa_current->weight[s] = PDF1;
						divisor += PDF1;
						Pa_card++;
					}
					/* now, we normalize probabilities and compute weights */
					for (int s = 0; s < d; s++) {
						long double PDF1 = Pa_current->weight[s];
						long double probability = (PDF1 / divisor) * (1 - PROB0) * PROBSLG;
						Pa_current->weight[s] = -log(probability);
					}

					/* compute candidate pairs, the criterion is to take components (b1,b2) such that:
					 * 1) the distance between b1 and b2 is less than 80 arcsec
					 * 2) component b1 can not be dimmer than b2 in 1 magnitude
					 * the item (1) is deactivated only if there are very few PPMX candidate stars (less than 4) */
					divisor = 0.0;
					for (int s2 = 0; s2 < d; s2++) {
						int j2 = neighbor[s2]; /* <--- component b1 */
						double delta2 = PPMXstar[j2].DE1875;
						double alpha2 = PPMXstar[j2].RA1875;
						double mag2 = PPMXstar[j2].mag;
						double x2, y2, z2;
						sph2rec(alpha2, delta2, &x2, &y2, &z2);

						for (int s3 = 0; s3 < d; s3++) {
							if (s3 == s2) continue;

							int j3 = neighbor[s3]; /* <--- component b2 */
							double delta3 = PPMXstar[j3].DE1875;
							double alpha3 = PPMXstar[j3].RA1875;
							double mag3 = PPMXstar[j3].mag;
							double x3, y3, z3;
							sph2rec(alpha3, delta3, &x3, &y3, &z3);
							double dist = calcdist(calcscalar(x2, y2, z2, x3, y3, z3));

							/* apply criterion */
							if (d > 3) if (dist > DISTTOLDPL) continue;
							if (mag2 - mag3 > 1.0) continue;
							if (Pa_card == MAXSIZEPA) bye("Memory exhausted for Pa structure!\n");

							/* here, we compute PDF2(0,0,0,0; a, b1, b2) */
							long double PDF2 = 1.0;
							PDF2 *= (long double)pdfn(alpha_difference(alpha1, alpha2, (delta1 + delta2) / 2.0), SIGMAALPHA);
							PDF2 *= (long double)pdfn(delta1 - delta2, SIGMADELTA);
							PDF2 *= (long double)pdfn(mag1 - mag2, SIGMAMAG);
							PDF2 *= (long double)pdfn(dist - DSEPPARAM, SIGMADSEPPARAM);
							PDF2 *= (long double)pdfn(mag2 - mag3, SIGMAMAGPARAM);

							Pa_current->b1[Pa_card] = j2;
							Pa_current->b2[Pa_card] = j3;
							Pa_current->weight[Pa_card] = PDF2;
							divisor += PDF2;
							Pa_card++;
						}
					}
					if (d == Pa_card) bye("No pair were found!\n");
					/* now, we normalize probabilities and compute weights */
					for (int s = d; s < Pa_card; s++) {
						long double PDF2 = Pa_current->weight[s];
						long double probability = (PDF2 / divisor) * (1 - PROB0) * (1 - PROBSLG);
						Pa_current->weight[s] = -log(probability);
					}

					/* add the "emptyset" element */
					Pa_current->b1[Pa_card] = -1;
					Pa_current->b2[Pa_card] = -1;
					Pa_current->weight[Pa_card] = -log(PROB0);
					Pa_card++;
					Pa_current->cardinal = Pa_card;
				}
			}

			/* check if weights define a probability distribution */
			long double sum = 0.0;
			for (int s = 0; s < Pa_current->cardinal; s++) sum += exp(-Pa_current->weight[s]);
			if (sum < 0.999999 || sum > 1.000001) {
				printf("Sum of exp(-weights) gives: %.6f\n", sum);
				bye("Internal error!\n");
			}

			/* copy Pa_current to Pa[i] */
			Pa[i].cardinal = Pa_current->cardinal;
			Pa[i].b1 = (int*)calloc(Pa_current->cardinal, sizeof(int));
			Pa[i].b2 = (int*)calloc(Pa_current->cardinal, sizeof(int));
			Pa[i].weight = (long double*)calloc(Pa_current->cardinal, sizeof(long double));
			for (int s = 0; s < Pa_current->cardinal; s++) {
				Pa[i].b1[s] = Pa_current->b1[s];
				Pa[i].b2[s] = Pa_current->b2[s];
				Pa[i].weight[s] = Pa_current->weight[s];
			}
			if (max_Pa_card < Pa_current->cardinal) max_Pa_card = Pa_current->cardinal;
		}
	}
	printf("Maximum cardinal of P_a: %d\n", max_Pa_card);

	/* generate components by using Kruskal-like algorithm based on disjoint set union structure */
	int cardinal = CDstars + PPMXstars;
	A1 = (int*)calloc(cardinal, sizeof(int));
	A2 = (int*)calloc(cardinal, sizeof(int));
	A3 = (int*)calloc(cardinal, sizeof(int));

	printf("Computing connected components...\n");
	for (int v = 0; v < cardinal; v++) {
		/* at first, every singleton is a component */
		A1[v] = v;
		A2[v] = 1;
		A3[v] = -1;
	}

	/* process edges */
	for (int e = 0; e < edges; e++) {
		int u = edge_A[e];            /* stars of CD are 0...CDstars-1 */
		int v = edge_B[e] + CDstars;  /* stars of PPMX are CDstars...CDstars+PPMXstars-1 */

		if (A1[u] != A1[v]) {
			/* if they belong to different components, make an union of them */
			int v1, v2;
			if (A2[A1[u]] < A2[A1[v]]) { /* choose the one with least elements */
				v1 = A1[u];
				v2 = A1[v];
			}
			else {
				v1 = A1[v];
				v2 = A1[u];
			}
			A2[v2] = A2[v1] + A2[v2]; /* add the size of both components */
			int r = A3[v2];
			A3[v2] = v1;
			int w = v2;
			/* modify the first vertices of the component with least elements */
			for (int i = 1; i <= A2[v1]; i++) {
				w = A3[w];
				A1[w] = v2;
			}
			A3[w] = r;
		}
	}
	
	/* count components */
	int unique_number = 0;
	int single_number = 0;
	int multiple_number = 0;
	int largest_single_component = 0;
	int largest_multiple_component = 0;
	for (int v = 0; v < cardinal; v++) {
		if (A2[v] > 1 && A1[v] == v) {
			/* we scan those components starting with "v" and having more than one element */
			int CDstars_in_component = 0;
			bool some_multiple = false;
			if (v < CDstars) {
				CDstars_in_component++;
				if (CDstar[v].dpl) some_multiple = true;
			}
			/* enter recursively into the component */
			int w = v;
			for (int i = 1; i < A2[v]; i++) {
				w = A3[w];
				if (w < CDstars) {
					CDstars_in_component++;
					if (CDstar[w].dpl) some_multiple = true;
				}
			}
			if (CDstars_in_component == 0) bye("Internal error!\n");
			if (CDstars_in_component == 1) unique_number++;
			else {
				if (some_multiple) {
					multiple_number++;
					if (largest_multiple_component < CDstars_in_component) largest_multiple_component = CDstars_in_component;
				}
				else {
					single_number++;
					if (largest_single_component < CDstars_in_component) largest_single_component = CDstars_in_component;
				}
			}
		}
	}
	printf("Number of \"unique star\" components: %d\n", unique_number);
	printf("Number of \"single stars\" components: %d,    largest one has %d CD stars\n", single_number, largest_single_component);
	printf("Number of \"multiple stars\" components: %d,    largest one has %d CD stars\n", multiple_number, largest_multiple_component);

	/* perform optimization */
	A_comp = (int*)calloc(MAXSIZEACOMP, sizeof(int));
	B_comp = (int*)calloc(MAXSIZEBCOMP, sizeof(int));
	hardest_variables = 0;
	hardest_constraints = 0;
	hardest_nodes = 0.0;
	hardest_time = 0.0;
	printf("Optimizing...\n");
	for (int v = 0; v < cardinal; v++) {
		if (A2[v] > 1 && A1[v] == v) {
			/* we optimize over those components starting with "v" and having more than one element */
			optimize(v);
		}
	}
	printf("Total time elapsed in optimization: %.4lf seconds\n", Tclock() - start_t);
	printf("Hardest model:\n  variables = %d\n  constraints = %d\n  nodes = %.0lf\n  time = %.4lf sec.\n",
		hardest_variables, hardest_constraints, hardest_nodes, hardest_time);

	/* check solution */
	bool *ppmx_used = (bool*)calloc(PPMXstars, sizeof(bool));
	int unassigned = 0; /* stars of CD not assigned (even having a neighbor in PPMX) */
	int singleassigned = 0; /* double stars of CD assigned to 1 PPMX star */
	int doubleassigned = 0; /* double stars of CD assigned to 2 PPMX star */
	int normalassigned = 0; /* normal stars of CD assigned to a PPMX star */
	for (int j = 0; j < PPMXstars; j++) ppmx_used[j] = false;
	for (int i = 0; i < CDstars; i++) {
		if (Pa[i].cardinal > 0) {
			int chosen = Pa[i].chosen;
			if (chosen < 0 || chosen >= Pa[i].cardinal) bye("Internal error!\n");
			int b1 = Pa[i].b1[chosen];
			int b2 = Pa[i].b2[chosen];
			if (b1 == -1) {
				/* the star of CD is not assigned because this is the most probable fact */
				unassigned++;
			}
			else {
				/* at least 1 star of PPMX was assigned to the CD star */
				if (ppmx_used[b1]) bye("PPMX star previously matched! Error in solution!\n");
				else ppmx_used[b1] = true;

				if (b2 == -1) {
					if (CDstar[i].dpl == true) singleassigned++;
					else normalassigned++;
				}
				else {
					/* 2 stars of PPMX were assigned */
					if (ppmx_used[b2]) bye("PPMX star previously matched! Error in solution!\n");
					else ppmx_used[b2] = true;

					doubleassigned++;
				}
			}
		}
	}
	printf("Statistics:\n  %d unassigned stars\n  %d single stars assigned\n  %d double stars assigned\n  %d double stars assigned to a single one in PPMX\n",
		unassigned, normalassigned, doubleassigned, singleassigned);

	/* save solution */
	stream = fopen("tables\\final_cross.txt", "wt");
	if (stream == NULL) {
		perror("Error writing final cross file");
		exit(1);
	}
	for (int i = 0; i < CDstars; i++) {
		if (Pa[i].cardinal == 0) {
			write_assignment(stream, i, -1, -1);
		}
		else {
			int chosen = Pa[i].chosen;
			write_assignment(stream, i, Pa[i].b1[chosen], Pa[i].b2[chosen]);
		}
	}
	fclose(stream);

	/* free memory */
	free(ppmx_used);
	free(B_comp);
	free(A_comp);
	free(A3);
	free(A2);
	free(A1);
	for (int i = 0; i < CDstars; i++) {
		if (Pa[i].cardinal > 0) {
			free(Pa[i].weight);
			free(Pa[i].b2);
			free(Pa[i].b1);
		}
	}
	free(neighbor);
	free(Pa_current->weight);
	free(Pa_current->b2);
	free(Pa_current->b1);
	free(Pa_current);
	free(Pa);
	free(degree_B);
	free(degree_A);
	free(PPMXstar);
	free(CDstar);
	return 0;
}