/* name of prog: polyfit
   project:      utility
   author:       Nils Bluemer
   version:      6a, 28.5.98
   interfaces:   command line: (see docu)
                 input: (depending on command line, see docu)
		      default: 3 columns (x-value, y-error, y-value)
		 parameters: (determined by command line)
		      minerror, maxerror: min/max error used for 
		                          weighting
		      maxpoints : max number of points for each fit
		      maxfits: max number of fits in a run
                 output: (depending on command line, see docu)
		      default: ?? 
		 error output: ??
   description:  *universal utility for polynomial fits up to order 9
                 *order of fit and free coefficients can be chosen on 
		  command line 
   new features: parameters can be read from data file      
                 lines beginning with # are treated as commentary 
                 lines
		 confidence level added
		 individual points can be fitted at reduced order
		 ( add #ro to input line )
		 output of desciptors with 2.3 digits
		 now output with exponent
   bug:          if all points of a fit consist have y=0, the error is 
                 wrong. workaround in line 199
   future:       !!!!! minerror must be determined dynamically !!!!	 
   */

/* 

The program reserves for each combination of descriptors ("dis[1]"..
"dis[numdis]" one value of the integer "index" and ... */
/*******************************************************************/

#include <ctype.h>
#include <math.h>       /* needed for sqrt */
#include "nrutil.c"     /* from NR: vector-definitions */
#include "nrproc.h"     /* from NR, partly modified: gaussj,covsrt,lfit
			 selfmade: ffunc (calculates monomes), test */

#define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;}
static double absarg;
#define ABS(a) ((absarg=(a)) >= 0.0 ? absarg : -absarg)
static double maxa, maxb, mina, minb;
#define MAX(a,b) ((maxa=(a)) >= (maxb=(b)) ? maxa : maxb)
#define MIN(a,b) ((mina=(a)) <= (minb=(b)) ? mina : minb)

#define maxpoints 20  /* max number of points */
#define maxfits 5000   /* max number of fits */
#define maxdis 4      /* max number of discriminators */
#define maxy 4         /* max number of functions to be fitted */
#define maxorder 9    /* max order of fit-polynomial */
#define fitall 1       
#define fiteven 2
#define fitodd 3
#define fituserdef 4
#define true 1
#define false 0

void error(char error_text[])
/* standard error handler */
{
    fprintf(stderr,"\npolyfit run-time error\n");
    fprintf(stderr,"%s\n",error_text);
    fprintf(stderr,"for general help use option -h\n");
    fprintf(stderr,"...now exiting to system...\n");
    exit(1);
}

int indexsearch (double dis[], double **disvec, int numfitpoints[], int *numfits, 
		 int numdis)
/* if combination of discriminators dis[1]..dis[numdis] has already 
been read, the corresponding index is returned, otherwise the next 
free index */
{
    int indx, i, found;
    found = false; 
    indx = 0;
    do {
	indx++;
	found = true;
	for (i=1;i<=numdis;i++)
	    if (disvec[indx][i] != dis[i])
		found = false;
    } while ((found == false) && (indx <= *numfits));
    if (indx > maxfits)
	error("indexsearch: too many comb. of determinators!");
    else {
	if (indx > *numfits) {
	    for (i=1;i<=numdis;i++)
		disvec[indx][i] = dis[i];
	    numfitpoints[indx] = 0;
	    (*numfits)++;
	}
	return indx;
    }
}

int larger (int x, int y, double **disvec, int numdis, int n)
/* compares disvec[][x] to disvec[][y] 
uses recursion, n is 1 at the original call and increases
at each recursive call */ 
{
    if (disvec[x][n] > disvec[y][n])
	return true;
    else {
	if (disvec[x][n] < disvec[y][n])
	    return false;
	else 
	    if (n < numdis)
		return larger (x,y,disvec,numdis,n+1);
	    else 
		return true;
		}
}

void sort (int a[], double **disvec, int n, int numdis)
/* shell sort algorithm (modified from NR S. 332),  sorts
according to order as defined by the function "larger" */
{
    int i,j,v,inc;
    for (i=1;i<= maxfits;i++){ /* initialisation of sort vector */
	a[i] = i;
    } 
    inc = 1;
    do {
	inc *=3;
	inc++;
    } while (inc <= n);
    do {
	inc /= 3;
	for (i=inc+1;i<=n;i++){
	    v=a[i];
	    j=i;
	    while (larger (a[j-inc],v,disvec,numdis,1)==true) {
		a[j]=a[j-inc];
		j -= inc;
		if (j <= inc) break;
	    }
	    a[j]=v;
	}
    } while (inc > 1);
}

void read (int *numfits, int *numfitpoints, double **disvec, 
	    double **xin, double **errorin, double ***yin,
	    double minerror, double maxerror, int numdis, int numy, int *flagor)
/* reads input line by line. lines beginning with # are ignored.
   each line is split in up to numdis+numy+2 strings containing
   only non-space characters
   these strings are interpreted as double numbers
 */
{
    int erg;
    double *dis, x,fehler,*y;
    int i,indx, sind;
    char c,line[100], *linep, *s, svec[maxdis+maxy+3][20];

    dis = dvector(1,numdis);
    y = dvector(1,numy);
    *numfits = 0;
    
    do {
	sind = 0;
	do {
	    s = gets (line);
	} while ((s) && (line[0]=='#'));
	if (s) {
	    linep = line;
	    strcpy(svec[numdis+numy+2],"");
	    do {
		while ((isspace(c = *linep))&& (c))
		    ++linep;
		if ((c) && (sind < numdis+numy+3))
		    sscanf(linep,"%s",svec[sind++]);
		while ((!isspace(c = *linep))&& (c))
		    ++linep;
	    } while ((c) && (sind <= numdis+numy+3));
	    if (sind >= numdis+2+numy){
		sind = 0;    
		for (i=1;i<=numdis;i++)
		    sscanf(svec[sind++], "%lf ",&dis[i]);
		sscanf(svec[sind++], "%lf",&x);
		sscanf(svec[sind++], "%lf",&fehler);
		for (i=1;i<=numy;i++)
		    erg = sscanf(svec[sind++], "%lf ",&y[i]);

		indx = indexsearch (dis, disvec, numfitpoints, numfits, numdis);
		
		if (strncmp(svec[sind++],"#ro",3)==0)
		    flagor[indx] = true;

		numfitpoints[indx]++;
		xin[indx][numfitpoints[indx]] = x;
		errorin[indx][numfitpoints[indx]] = MAX (fehler, minerror);
		for (i=1;i<=numy;i++)
		    yin[i][indx][numfitpoints[indx]] = y[i]+0.0000001;
	    }
	    else 
		if (sind != 0)  /* no empty line */
		    fprintf(stderr,"polyfit/read: input format error: not enough columns\n");
	}
    } while (s);
    free_dvector(dis,1,numdis);
    free_dvector(y,1,numy);
    }

singlefit (double x[],double y[],double sig[],int ndat, double *c, double *sigc, 
      int order, int fitvec[], int degfree, double *q)
{
    int ma, i, *ia;
    double *a, **cov, chisq;

    ma = order+1;      /* number of coefficients is order of fit +1 */

    cov = dmatrix (1,ma,1,ma);
    a = dvector (1,ma);
    ia = ivector (1,ma);

    for (i=0;i<=order;i++){
	ia[i+1] = fitvec[i];
	a[i+1] = 0;
    }

    for (i=1;i<=ndat;i++) /* data points with value 0 are errors */
	if (y[i] == 0.0)
	    sig[i] = 1;

    lfit (x,y,sig,ndat,a,ia,ma,cov,&chisq);
    for (i=0;i<=order;i++){
	c[i] = a[i+1];
	sigc[i] = sqrt(cov[i+1][i+1]);
    }
    *q = gammq(0.5*degfree,0.5*chisq);

    free_dmatrix (cov,1,ma,1,ma);
    free_dvector(a,1,ma);
    free_ivector(ia,1,ma);
}

void fit (int numfits, int *numfitpoints, double **xin, double **errorin, 
	     double ***yin,double ***coeff, double **sigcoeff, int order, int fitvec[],
	     int numy, int degfree, int reduce, double **q, int *flagor)
{
    int i,indx,temporder,tempdegfree, j;
    double y,sigy;
    for (indx=1;indx <= numfits;indx++) {
	if (numfitpoints[indx]>degfree&&(flagor[indx]==false||reduce==false))
	    tempdegfree = degfree;
	else 
	    if ((numfitpoints[indx] >2)&&(numfitpoints[indx] >= degfree-1)&&
		reduce==true)
		tempdegfree = MIN(degfree-1, numfitpoints[indx]-1);
	    else
		tempdegfree = 0; /* kein fit */	
	if (tempdegfree < degfree){ /* determine order */
	    flagor[indx] = true;
	    j = fitvec[0];
	    temporder = 0;
	    do{
		j += fitvec[++temporder];
	    } while (j <tempdegfree);
	}
	else
	    temporder = order;

	if (tempdegfree > 0){
	    for (i=1;i<=numy;i++){
		for (j=0;j<= temporder;j++)
		    coeff[i][indx][j] = sigcoeff[indx][j] = 0;
		singlefit (xin[indx],yin[i][indx],errorin[indx],
			   numfitpoints[indx], coeff[i][indx], sigcoeff[indx], 
			   temporder, fitvec, tempdegfree, &q[i][indx]);
	    }
	}
	else
	    fprintf(stderr,"not enough fitpoints\n");
    }
}
    

void printhelp ()
{
 printf("**********************************************************\n");
 printf("\npolyfit6: utility for fitting data on polynomials\n");
 printf("input:  standardin\n");
 printf("options: -l  (low: minimum error level /= 2)\n");
 printf("         -u  (up:  minimum error level *= 2)\n");
 printf("         -o# order of fit (#=0..9)\n");
 printf("         -d# number of discriminators (#=0..%d)\n",maxdis);
 printf("         -y# number of functions to be fitted (#=0..%d)\n",
	maxy);
 printf("         -fa fit all coefficients\n");
 printf("         -fe fit only even coefficients (0,2,..)\n");
 printf("         -fo fit only odd coefficients\n");
 printf("         -fu1011..0..1 userdefine fit coefficients\n");
 printf("         -h  this help\n");
 printf("         -cd interprete all expressions starting with- in first\n");
 printf("             line of data file as belonging to command line\n");
 printf("         -p  print all parameters in output\n");
 printf("         -w  unweighted fit (not impl.)\n");
 printf("         -r  reduce order if not enough fitpoints\n");
 printf("             or indicated in inputline by appending #ro\n");
 printf("Options can be concatenated (with only one -) and evoced\n");
 printf("in any order. \n\n");
 printf("format conventions:\n");
 printf("input:  dis[1]..dis[d#], x, yerror, y[y#]..y[y#]\n\n");
 printf("output: dis[1]..dis[d#], coeff[1][0]..coeff[1][o#]..coeff[y#][o#],\n");
 printf("        confidence level, #fitpoints, sigcoeff[0]..sigcoeff[o#]\n\n");

 printf("**********************************************************\n");
}

void printpar (int fittype, int order, int numdis, int numy, 
		 double minerror, int reduce)
{
    printf ("# Weighted fit of order %d in ",order);
    switch (fittype){
    case fitall: printf("all"); break;
    case fiteven: printf("even"); break;
    case fitodd: printf("odd"); break;
    case fituserdef: printf("userdefined"); break;
    }
    printf (" coefficients. Low error cutoff: %7.5lf\n", minerror);
    printf ("# Number of discriminators: %d, number of functions: %d, ",
	    numdis, numy);
    printf ("order reduction: ");
    if (reduce == true)
	printf ("on.\n");
    else
	printf ("off.\n");
}

void process (double maxerror, double minerror, int order, int fitvec[],
		 int numdis, int numy, int degfree, int fittype, int print, int reduce)
{
    int indx, numfits;
    double **disvec, F;
    int *numfitpoints, *sortvec, *flagor; 
        /*flagor: order reduction for some points */
    double **xin, **errorin;
    double ***yin, ***coeff, **sigcoeff, **q;
    int i, j, k;

    numfitpoints = ivector (1,maxfits);
    flagor = ivector (1,maxfits);
    for (i=1;i<=maxfits;i++)
	flagor[i] = false;
    disvec = dmatrix(1,maxfits,1,maxdis);
    sortvec = ivector(1,maxfits);

    xin = dmatrix(1,maxfits,1,maxpoints);
    errorin = dmatrix(1,maxfits,1,maxpoints);
    q = dmatrix(1,numy,1,maxfits);
    yin = d3tensor(1,numy,1,maxfits,1,maxpoints);
    coeff = d3tensor(1,numy,1,maxfits,0,order);
    sigcoeff = dmatrix(1,maxfits,0,order);

    read (&numfits,numfitpoints,disvec,xin,errorin,yin, minerror, maxerror, 
	   numdis, numy, flagor);
    
    fit (numfits,numfitpoints,xin,errorin,yin,coeff,sigcoeff, order, fitvec, numy, 
	 degfree, reduce, q, flagor);

    sort (sortvec, disvec, numfits, numdis);

    if (print == true)
	printpar(fittype,order, numdis, numy, minerror, reduce);
    for (j = 1; j <= numfits; j++) {
	indx = sortvec[j];
	if (sigcoeff[indx][0] == 0.0)  /* not enough points: not fitted */
	    printf("# ");
	for (i=1;i<=numdis;i++)
	    printf("%- 14.8g ",disvec[indx][i]);
	for (i=1;i<=numy;i++){
	    for (k=0;k<=order;k++)
		printf("%- 14.8g ",coeff[i][indx][k]);
	    printf (" %- 8.2g ",q[i][indx]);
	}
	printf("%2d ", numfitpoints[indx]); 
	for (k=0;k<=order;k++)
	    printf("%- 10.4g ", sigcoeff[indx][k]); 
	if (reduce==true&&flagor[indx]==true)
	    printf(" order_reduced");
	printf("\n");
    } 
    free_dmatrix(disvec,1,maxfits,1,maxdis);
    free_ivector(numfitpoints,1,maxfits);
    free_ivector(flagor,1,maxfits);
    free_ivector(sortvec,1,maxfits);
    free_dmatrix(xin,1,maxfits,1,maxpoints);
    free_dmatrix(errorin,1,maxfits,1,maxpoints);
    free_dmatrix(q,1,numy,1,maxfits);
    free_d3tensor(yin,1,numy,1,maxfits,1,maxpoints);
    free_d3tensor(coeff,1,numy,1,maxfits,0,order);
    free_dmatrix(sigcoeff,1,maxfits,0,order);
}


void parse_cl (char *s, int *fittype, int *order, char fitstring[],
	       int *numdis, int *numy, 
	       double *maxerror, int *print, int *reduce, int recurse)
/* parses command line from program call or (if option -cd has
been on command line) from the first line of the data file 
(which must in this case begin with #) */
{
    char c, line[100], s2[100], *li, *ss;
    int i;
    while (c= *++s)
	switch (c) {
	case 'c':
	    c= *++s;
	    if ((c=='d')&&(recurse == false)){
		gets (line);
		if (line[0]!='#')
		    error ("parse_cl: no valid command line in data file (first char is not a #)");
		else {
		    li = line;
		    while (c = *++li)
			if (c == '-'){
			    sscanf(li,"%s",s2);
			    parse_cl (s2, fittype, order, fitstring, numdis, numy, maxerror, 
					print, reduce, true); 
			}
		}
	    }
	    else
		error ("parse_cl: option 'c' not followed by 'd'");
	    break;
	case 'd':
	    c= *++s;
	    if ((c>='0')&&(c<='0'+maxdis))
		*numdis = (int) c - (int)'0';
	    else error("parse_cl: option 'd' not followed by number or number too large");
	    break;
	case 'y':
	    c= *++s;
	    if ((c>='1')&&(c<='0'+maxy))
		*numy = (int) c - (int)'0';
	    else error("parse_cl: option 'y' not followed by number or number too large");
	    break;
	case 'o':
	    c= *++s;
	    if ((c>='0')&&(c<='9'))
		*order = (int) c - (int)'0';
	    else error("parse_cl: option 'o' not followed by number or number too large");
	    break;
	case 'f':
	    c= *++s;
	    switch (c) {
	    case 'a': 
		*fittype = fitall;
		break;
	    case 'e': 
		*fittype = fiteven;
		break;
	    case 'o': 
		*fittype = fitodd;
		break;
	    case 'u': 
		*fittype = fituserdef;
		ss=s;
		i=0;
		while ((c=*++ss)&&((c=='1')||(c=='0'))){
		    s = ss;
		    if (i < maxorder)
			fitstring[i++]=c;
		}
		fitstring[i]='\0';
		break;
	    default: 
		error("parse_cl: option 'f' not followed by 'a','e','o' or 'u'");
		break;
	    }
	        break;
	case 'l':
	    *maxerror = *maxerror / 2;
	    break;
	case 'u':
	    *maxerror = *maxerror * 2;
	    break;
	case 'p':
	    *print = true;
	    break;
	case 'r':
	    *reduce = true;
	    break;
	case 'w':
		error("parse_cl: parameter option -w not yet implemented");
	    break;
	default:
	    printhelp();
	    exit(0);
	}
}

void determine_fitvec (int *order, int fittype, char fitstring[], int fitvec[], 
		       int *degfree)
{
    int i;
    char c;
    *degfree = 0;
    for (i=0;i<=*order;i++)
      fitvec[i]=0;
    if (fittype == fitall)
	for (i = 0; i <= *order;i++){
	    fitvec[i] = 1;
	    ++(*degfree);
	}
    if (fittype == fiteven)
	for (i = 0; i <= *order;i+=2){
	    fitvec[i] = 1;
	    ++(*degfree);
	}
    if (fittype == fitodd)
	for (i = 1; i <= *order;i+=2){
	    fitvec[i] = 1;
	    ++(*degfree);
	}
    if (fittype == fituserdef){
		*order=0;
		while ((c=*(fitstring++))&&(*order<=maxorder)){
		    fitvec[(*order)++] = (int) c - (int)'0';
		    if (c=='1')
			++(*degfree);
		}
		(*order)--;
    }
}

void main (int argc, char *argv[])
{
    double maxerror, minerror;
    int c, order, *fitvec, fittype;
    int  numdis, numy, degfree, print, reduce;
      /* number of discriminators, number of functions to be fitted
	number of degrees of freedom in fit */
    char fitstring[maxorder+1];

    maxerror = 0.002;
    order = 2;
    fittype = fitall;
    numdis = 0;
    numy = 1;
    print = false;
    reduce = false;

    while (--argc > 0 && (*++argv)[0] == '-')
	parse_cl (argv[0], &fittype,&order, fitstring,&numdis, &numy, &maxerror, &print, &reduce,
		  false);
    minerror = maxerror / 5;
    fitvec = ivector(0,maxorder); 
    determine_fitvec (&order,fittype,fitstring,fitvec,&degfree);
    process (maxerror, minerror, order, fitvec, numdis, numy, degfree, fittype,
	     print, reduce);  
    free_ivector(fitvec,0,maxorder);
    exit (0);
}