/*
 *	elem 0.5 (C) 2001 Ivan Tubert (itubert@hotmail.com)
 *	
 *	You can do whatever you want with this program except claiming you wrote it.
 *	
 *	This program has ABSOLUTELY NO WARRANTY
 *
 */


#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "eprintf.h"

enum{
	TRUE = 1,
	FALSE = 0,
	MAXBLOCKS = 20,
	SYMLENGTH = 8,
	MAXELEM = 256,
	LINELENGTH = 128,
	MAXBEST = 1000
};

typedef int BOOL;

typedef struct {
	char symbol[SYMLENGTH];
	double mass;
} ELEMENT;

typedef struct {
	int type;
	double percentage;
} COMPOSITION_ELEMENT;

typedef struct {
	int n;
	COMPOSITION_ELEMENT *elem;
} COMPOSITION;

typedef struct {
	int elem;
	int index;
} FORMULA_ITEM;

typedef struct {
	double min, max, step, curr;
	FORMULA_ITEM *formula;
	int n_elem;		// number of different elements in the formula
	double charge;
	char *formula_string;
	char *name;
} BLOCK;

typedef struct {
	COMPOSITION *comp,
		*rel_err;	// rel_err[i] = (exp[i] - theo[i]) / theo[i]
	double sq_err;  // sq_err = sum(sq(exp-theo))
	double charge;
	double vec[MAXBLOCKS];  // vector describing the composition
} BLOCKSUM;


/*** GLOBAL VARIABLES ***/
ELEMENT elements[MAXELEM];
int NELEM;				// Number of elements
int NBLOCK = 0;		//  The number of components to be varied (aka blocks)
int CHARGEREQ = FALSE;	// Requirement to find only electrically neutral combinations;
int NBEST = 10;
int VERBOSITY = 0;
/* Verbosity levels:
0	normal
1	verbose
2	very verbose
3	extremely verbose
-1	silent
*/

/* Returns number of elements read from file fname
 * Entries are in the format "symbol mass ..." */
int read_elements(char *fname)
{
	FILE *fp;
	char s[SYMLENGTH];
	int n = 0;
	float f;
	
	memset(elements, 0, sizeof(elements));

	fp = fopen(fname, "r");

	if(fp) {
		while(fscanf(fp, "%s %f", s, &f) == 2) { // in God we trust
			elements[n].mass = f;
			strcpy(elements[n].symbol, s);
			++n;
		}
	} else {
		printf("Could not open element file %s\n", fname);
	}
	return(n);
}

void parse(BLOCK *b)
{
	int n;
	char *p, *endptr;
	int i, len, maxlen, z, charge = 0;

	p = b->formula_string;
	while(*p){
		/*** Read symbol ***/
		
		for(i = maxlen = 0; (i < NELEM); ++i){
			len = strlen(elements[i].symbol);
			if((len>maxlen) && (strncmp(elements[i].symbol, p, len) == 0)){
				maxlen = len;
				z = i;
			}
		}
		if(!maxlen) {
			printf("Unidentified symbol in %s\n", b->formula_string);
			printf("                       ");
			for(i=0; i < p - b->formula_string; ++i)
				printf(" ");
			printf("^\n");
			exit(1);
		}
		p += maxlen;

		if(VERBOSITY > 1)
			printf("    Symbol: %s\n", elements[z].symbol);

		/*** Read index ***/
		if(*p != '+' && *p != '-') { // make sure that number is not intended as a charge
			n = strtol(p, &endptr, 10);
			if(endptr == p) {
				n = 1;
			} else
				p = endptr;
		} else
			n = 1;

		if(VERBOSITY > 1)
			printf("    Index: %i\n", n);


		if(*p == '+' || *p == '-') {
			charge = strtol(p, &endptr, 10);
			if(endptr == p) {
				charge = (*p == '+') ? 1 : -1;
				++p;
			} else
				p = endptr;
		}
		if(VERBOSITY > 1)
			printf("    Charge: %i\n", charge);
		
		// Add to formula
		b->charge = charge;
		for(i=0; i<b->n_elem; ++i){
			if(b->formula[i].elem == z){
				b->formula[i].index += n;
				break;
			}
		}
		if(i==b->n_elem){ // add to list
			b->formula = realloc(b->formula, (b->n_elem+1)*sizeof(FORMULA_ITEM));
			b->formula[i].elem = z;
			b->formula[i].index = n;
			++b->n_elem;
		}
	}
}



int inc(BLOCK *b, int n)
{
	if(n<0)
		return(FALSE);

	b[n].curr += b[n].step;
	if(b[n].curr - 1.0E-6 > b[n].max){
		b[n].curr = b[n].min;
		return(inc(b, n-1));
	}
	return(TRUE);
}

/* Calculates composition and returns newly allocated structure */
COMPOSITION *calc_comp(BLOCK *b)
{
	COMPOSITION *comp;
	int i, j;
	double MW = 0.0;
	double fcomp[MAXELEM];
	double f;

	comp = calloc(1, sizeof(COMPOSITION));
	memset(&fcomp, 0, sizeof(fcomp));

	for(i=0; i < NBLOCK; ++i){
		for(j=0; j < b[i].n_elem; ++j) {
			f = b[i].formula[j].index * elements[b[i].formula[j].elem].mass * b[i].curr;
			MW += f;
			fcomp[b[i].formula[j].elem] += f; 
		}
	}
	// compress the list
	for(i = 0; i < NELEM; ++i) {
		if(fcomp[i]){
			comp->elem = realloc(comp->elem, (comp->n + 1) * sizeof(COMPOSITION_ELEMENT));
			comp->elem[comp->n].type = i;
			comp->elem[comp->n].percentage = fcomp[i] / MW * 100;
			++comp->n;
		}
	}

	return(comp);
}


double calc_charge(BLOCK *b)
{
	double charge = 0.0;
	int i;
	
	for(i=0; i < NBLOCK; ++i)
		charge += b[i].charge * b[i].curr;
	
	return(charge);
}


COMPOSITION *error(COMPOSITION *exp, COMPOSITION *theo, double *sq_err)
{
	double c1[MAXELEM];
	double c2[MAXELEM];
	double d = 0.0;
	int i;
	COMPOSITION *comp;

	comp = calloc(1, sizeof(COMPOSITION));

	for(i = 0; i < NELEM; ++i)
		c1[i] = -1.0; // negative number indicates undefined element
	for(i = 0; i < exp->n; ++i)
		c1[exp->elem[i].type] = exp->elem[i].percentage;

	memset(&c2, 0, sizeof(c2));
	for(i = 0; i < theo->n; ++i)
		c2[theo->elem[i].type] = theo->elem[i].percentage;

	for(i = 0; i < NELEM; ++i)
		if(c1[i] >= 0.0) {
			d += (c1[i] - c2[i]) * (c1[i] - c2[i]);
			comp->elem = realloc(comp->elem, (comp->n + 1) * sizeof(COMPOSITION_ELEMENT));
			comp->elem[comp->n].type = i;
			comp->elem[comp->n].percentage = c1[i] = (c1[i] - c2[i]) / c2[i] * 100;
			++comp->n;
		}
	
	if(sq_err)
		*sq_err = d;

	return(comp);
}


BOOL addtobest(BLOCKSUM *best, BLOCK *b, COMPOSITION *compteo, COMPOSITION *compexp)
{
	double sq_err, charge;
	int i, j;
	COMPOSITION *rel_err;

	rel_err = error(compexp, compteo, &sq_err); // Theological line
	charge = calc_charge(b);

	for(i=0; i<NBEST; ++i) {
		if((sq_err < best[i].sq_err) && (!CHARGEREQ || fabs(charge) < 1E-6))
		{
			memmove(&best[i+1], &best[i], sizeof(BLOCKSUM) * (NBEST - i - 1));
			best[i].sq_err = sq_err;
			best[i].rel_err = rel_err;
			best[i].comp = compteo;
			best[i].charge = charge;
			for(j = 0; j < NBLOCK; ++j){
				best[i].vec[j] = b[j].curr;
			}
			return(TRUE);
		}
	}
	return(FALSE);
}

void *destroy_comp(COMPOSITION *c)
{
	free(c->elem);
	free(c);
	return(NULL);
}


char *readline(char *line, FILE *fp)
{
	int len;

	fgets(line, LINELENGTH, fp);
	if(feof(fp))
		eprintf("Error: Unexpected end of file\n");

	len = strlen(line);
	if((len == LINELENGTH) && (line[len-1] != '\n'))
		eprintf("Error: line too long.\n");
	return(line);
}

void printcomp(COMPOSITION *c)
{
	int i;

	for(i = 0; i < c->n; ++i){
		printf("%s=%5.2f%%; ", elements[c->elem[i].type].symbol, c->elem[i].percentage);

	}

}

int main(int argc, char **argv)
{
	FILE *fp;
	char *fname = NULL;
	COMPOSITION *exp, *curr;
	BLOCK b[MAXBLOCKS];
	int i, j;
	BLOCKSUM *best;
	float percentage, min, max, step;
	char line[LINELENGTH];
	char symbol[SYMLENGTH];
	char name[100];
	char formula_string[100];
	int n = 0;	/* Number of trials */
	char *arg;
		
	memset(b, 0, sizeof(b));


/*** READ COMMAND LINE ARGUMENTS ***/
	for(i = 1; i < argc; ++i) {
		arg = argv[i];

		if(arg[0] == '-') { //options
			char *p;
			for(p = arg+1; *p; ++p)
				switch(*p){
					case '0': // zero
						CHARGEREQ = TRUE;
						break;

					case 'h': 
						printf("elem v. 0.5 by Ivan Tubert\n");
						printf("use: elem [options...] [filename]\n");
						printf("\nOptions:\n");
						printf("-niii Display 'iii' results\n");
						printf("-0    Show only electrically neutral results\n");
						printf("-h    Show this help and exit\n");
						printf("-v    Verbose mode\n");
						printf("-vv   Very verbose mode\n");
						printf("-s    Silent mode\n");
						exit(0);
						break;

					case 'v':
						++VERBOSITY;
						break;

					case 's':
						--VERBOSITY;
						break;

					case 'n': {
							int m;
							char *endptr;

							if(p[1]){
								m = strtol(p+1, &endptr, 10);
								p = endptr;
								if(!*p)
									--p;
							} else {
								++i;
								if(i<argc){
									p = arg = argv[i];
									m = strtol(p, &endptr, 10);
									p += strlen(arg)-1;
								}

							}
							if(m > 0 && m <= MAXBEST)
								NBEST = m;
							else
								weprintf("Expected a number in [1, %i] after 'n' option. Using default of %i.", MAXBEST, NBEST);
						}break;

					default:
						weprintf("Warning: unrecognized option '%c'\n", *p);
						break;
				}
		} else
			if(!fname)
				fname = arg;
			else
				weprintf("Ignoring multiple files: %s", arg);
	}

	best = calloc(NBEST, sizeof(BLOCKSUM));
	for (i = 0; i < NBEST; ++i)
		best[i].sq_err = 1E20;

	NELEM = read_elements("elem.txt");

/*** READ INPUT FILE ***/

/* Open file */
	if(fname) {
		if((fp = fopen(fname, "r")) == 0)
			eprintf("Error opening file %s\n", fname);
	} else {
		fp = stdin;
		fname = "STDIN";
	}

	if(VERBOSITY > -1)
		printf("Input file: %s\n", fname);

	do{ /* Ignore everything up to first dot ".exp" */
		readline(line, fp);
	} while(line[0] != '.');


 /* Read experimental composition */

	if(VERBOSITY > -1)
		printf("Experimental composition:\n");
	exp = calloc(1, sizeof(COMPOSITION));
	do{
		readline(line, fp);
		if(line[0] != '.') {
			if(sscanf(line, "%s %f", symbol, &percentage) != 2)  // In God we trust
				eprintf("Error scanning in .EXP block");
			for(i = 0; i < NELEM; ++i){
				if(strncmp(symbol, elements[i].symbol, SYMLENGTH) == 0) {
					exp->elem = erealloc(exp->elem, (exp->n + 1) * sizeof(COMPOSITION_ELEMENT));
					exp->elem[exp->n].type = i;
					exp->elem[exp->n].percentage = percentage;
					++exp->n;
					if(VERBOSITY > -1)
						printf("%s: %5.2f\n", elements[i].symbol, percentage);
					break;
				}
			}
			if(i == NELEM)
				eprintf("Unidentified element in .EXP block");
		}
	} while(line[0] != '.');

	
	if(VERBOSITY > 0) {
		printf(".EXP block done\n");
	}

	while(fscanf(fp, "%s %s %f %f %f", name, formula_string, &min, &max, &step) == 5) {
		b[NBLOCK].min = min;
		b[NBLOCK].max = max;
		b[NBLOCK].step = step;
		b[NBLOCK].name = estrdup(name);
		b[NBLOCK].formula_string = estrdup(formula_string);
		if(VERBOSITY > -1)
			printf("%s (%s): from %5.2f to %5.2f step %5.2f\n", b[NBLOCK].name, b[NBLOCK].formula_string, b[NBLOCK].min, b[NBLOCK].max, b[NBLOCK].step);
		parse(&b[NBLOCK]);
		b[NBLOCK].curr = b[NBLOCK].min;

		++NBLOCK;
	}


	if(VERBOSITY > 0)
		printf("N = %i\n\n", NBLOCK);

	if(VERBOSITY > 0)
		printf("\nStarting...\n");

/*** PERFORM SEARCH ***/
	do{
		curr = calc_comp(b);
		if(!addtobest(best, b, curr, exp))
			destroy_comp(curr);
		++n;
	} while(inc(b, NBLOCK-1));


/*** DISPLAY RESULTS ***/
	
	if(VERBOSITY > 0)
		printf("Tried %i combinations.\n", n);

	if(VERBOSITY > -1)
		printf("BEST %i RESULTS\n", NBEST);
	
	for(i = 0; (i < NBEST) && (best[i].sq_err < 1E19); ++i) {
		for(j=0; j < NBLOCK; ++j){
			printf("%s: %4.2f; ", b[j].name, best[i].vec[j]);
		}
		printf("Charge=%4.2f\n", best[i].charge);
		printcomp(best[i].comp);
		printf("\nTotal squared absolute error = %.4g\n", best[i].sq_err);
		printf("Relative error by element: "); 
		printcomp(best[i].rel_err);
		printf("\n\n");
	}

	return 0;
}