#ifndef lint
static char *RCSid = "$Id: fit.c,v 1.33 1995/12/11 23:14:04 drd Exp $";
#endif

/*
 *	Nonlinear least squares fit according to the
 *	Marquardt-Levenberg-algorithm
 *
 *	added as Patch to Gnuplot (v3.2 and higher)
 *	by Carsten Grammes
 *	Experimental Physics, University of Saarbruecken, Germany
 *
 *	Internet address: cagr@rz.uni-sb.de
 *
 *	Copyright of this module:   Carsten Grammes, 1993
 *
 *	Permission to use, copy, and distribute this software and its
 *	documentation for any purpose with or without fee is hereby granted,
 *	provided that the above copyright notice appear in all copies and
 *	that both that copyright notice and this permission notice appear
 *	in supporting documentation.
 *
 *	This software is provided "as is" without express or implied warranty.
 *
 *	930726:     Recoding of the Unix-like raw console I/O routines by:
 *		    Michele Marziani (marziani@ferrara.infn.it)
 * drd: start unitialised variables at 1 rather than NEARLY_ZERO
 *  (fit is more likely to converge if started from 1 than 1e-30 ?)
 */


#define FIT_MAIN

#include <math.h>
#include <ctype.h>
#include <signal.h>

#if defined(MSDOS) || defined(DOS386)	     /* non-blocking IO stuff */
#include <io.h>
#include <conio.h>
#include <dos.h>
#else
#if defined(AMIGA_SC_6_1)
#include <fcntl.h>
#include <ios1.h>
#include <exec/types.h>
#include <exec/tasks.h>
#include <dos/dos.h>
#include <proto/exec.h>
#include <proto/dos.h>
#else
#ifdef OSK
#include <stdio.h>
#else
#ifndef VMS
#include <fcntl.h>
#endif	/* VMS */
#endif /* OSK */
#endif	/* AMIGA */
#endif	/* DOS */

#if defined(ATARI) || defined(MTOS)
#include <osbind.h>
#include <stdio.h>
#include <time.h>
#define getchx() Crawcin()
int kbhit(void);
#endif

#define STANDARD    stderr	/* compatible with gnuplot philosophy */

#include "plot.h"
#include "type.h"               /* own types */
#include "matrix.h"
#include "fit.h"
#include "setshow.h"   /* for load_range */

enum marq_res { OK, ERROR, BETTER, WORSE };
typedef enum marq_res marq_res_t;


#define INFINITY    1e30
#define NEARLY_ZERO 1e-30
#define INITIAL_VALUE 1.0  /* create new variables with this value (was NEARLY_ZERO) */
#define DELTA	    0.001
#define MAX_DATA    2048
#define MAX_PARAMS  32
#define MAX_LAMBDA  1e20
#define MAX_VALUES_PER_LINE	32
#define MAX_VARLEN  32
#define START_LAMBDA 0.01
#if defined(MSDOS) || defined(OS2) || defined(DOS386)
#define PLUSMINUS   "\xF1"                      /* plusminus sign */
#else
#define PLUSMINUS   "+/-"
#endif
#define TEMPNAME "gfttemp.$$$"			/* used for msdos and os/2 */

static double epsilon	   = 1e-5;		/* convergence criteria */

static char fit_script[127];
static char logfile[128]    = "fit.log";
static char *fit_index_def  = "FIT_INDEX";
static char *FIXED	    = "# FIXED";
static char *GNUFITLOG	    = "FIT_LOG";
static char *FITLIMIT	    = "FIT_LIMIT";
static char *FITSCRIPT	    = "FIT_SCRIPT";
static char *DEFAULT_CMD    = "replot";         /* if no fitscript spec. */


static FILE	    *log_f = NULL;

static word	    num_data,
		    num_params;
static double	    *fit_x,
		    *fit_y,
		    *err_data,
		    *a;
static boolean	    ctrlc_flag = FALSE;

static struct udft_entry func;

typedef char fixstr[MAX_VARLEN+1];
static fixstr	    *par_name;


/*****************************************************************
			 internal Prototypes
*****************************************************************/

static RETSIGTYPE ctrlc_handle __P((int an_int));
static void ctrlc_setup __P((void));
static marq_res_t marquardt __P((double a[], double **alpha, double *chisq,
			        double *lambda, double *varianz));
static boolean analyze __P((double a[], double **alpha, double beta[],
			    double *chisq, double *varianz));
static void calculate __P((double *yfunc, double **dyda, double a[]));
static void call_gnuplot __P((double *par, double *data));
static void show_fit __P((int i, double chisq, double last_chisq, double *a,
		          double lambda, FILE *device));
static boolean regress __P((double a[]));
#if 0 /* not used */
static int scan_num_entries __P((char *s, double vlist[]));
#endif
static boolean is_variable __P((char *s));
static void copy_max __P((char *d, char *s, int n));
static double getdvar __P((char *varname));
static double createdvar __P((char *varname, double value));
static void splitpath __P((char *s, char *p, char *f));
#if defined(MSDOS) || defined(DOS386)
static void subst __P((char *s, char from, char to));
#endif
static boolean fit_interrupt __P((void));
static boolean is_empty __P((char *s));

/*****************************************************************
    Useful macros
    We avoid any use of varargs/stdargs (not good style but portable)
*****************************************************************/
#define Dblf(a) 	{fprintf (STANDARD,a); fprintf (log_f,a);}
#define Dblf2(a,b)	{fprintf (STANDARD,a,b); fprintf (log_f,a,b);}
#define Dblf5(a,b,c,d,e) \
		{fprintf (STANDARD,a,b,c,d,e); fprintf (log_f,a,b,c,d,e);}


/*****************************************************************
    This is called when a SIGINT occurs during fit
*****************************************************************/

static RETSIGTYPE ctrlc_handle (an_int)
int an_int;
{
#ifdef OS2
    (void) signal(an_int, SIG_ACK);
#else
    /* reinstall signal handler (necessary on SysV) */
    (void) signal(SIGINT, (sigfunc)ctrlc_handle);
#endif
    ctrlc_flag = TRUE;
}


/*****************************************************************
    setup the ctrl_c signal handler
*****************************************************************/
static void ctrlc_setup()
{
/*
 *  MSDOS defines signal(SIGINT) but doesn't handle it through
 *  real interrupts. So there remain cases in which a ctrl-c may
 *  be uncaught by signal. We must use kbhit() instead that really
 *  serves the keyboard interrupt (or write an own interrupt func
 *  which also generates #ifdefs)
 *
 *  I hope that other OSes do it better, if not... add #ifdefs :-(
 */
#if (defined(__EMX__) || !defined(MSDOS) && !defined(DOS386)) && !defined(ATARI) && !defined(MTOS)
	(void) signal(SIGINT, (sigfunc)ctrlc_handle);
#endif
}


/*****************************************************************
    getch that handles also function keys etc.
*****************************************************************/
#if defined(MSDOS) || defined(DOS386)
int getchx ()
{
    register int c = getch();
    if ( !c  ||  c == 0xE0 ) {
	c <<= 8;
	c |= getch ();
    }
    return c;
}
#endif

/*****************************************************************
    in case of fatal errors
*****************************************************************/
void error_ex ()
{
    char *sp;

    strncpy (fitbuf, "         ", 9);         /* start after GNUPLOT> */
    sp = strchr (fitbuf, '\0');
    while ( *--sp == '\n' )
	;
    strcpy (sp+1, "\n\n");              /* terminate with exactly 2 newlines */
    fprintf (STANDARD, fitbuf);
      if ( log_f ) {
	fprintf (log_f, "BREAK: %s", fitbuf);
  	fclose (log_f);
  	log_f = NULL;
      }
    if ( func.at ) {
	free (func.at); 		    /* release perm. action table */
	func.at = (struct at_type *) NULL;
    }

    /* restore original SIGINT function */
#if !defined(ATARI) && !defined(MTOS)
    interrupt_setup ();
#endif

    longjmp(env, TRUE); 	/* bail out to command line */
}


/*****************************************************************
    Marquardt's nonlinear least squares fit
*****************************************************************/
static marq_res_t marquardt (a, alpha, chisq, lambda, varianz)
double a[];
double **alpha;
double *chisq;
double *lambda;
double *varianz;
{
    int 	    j;
    static double   *da,	    /* delta-step of the parameter */
		    *temp_a,	    /* temptative new params set   */
		    **one_da,
		    *beta,
		    *tmp_beta,
		    **tmp_alpha,
		    **covar,
		    old_chisq;

    /* Initialization when lambda < 0 */

    if ( *lambda < 0 ) {	/* Get first chi-square check */
	one_da	    = matr (num_params, 1);
	temp_a	    = vec (num_params);
	beta	    = vec (num_params);
	tmp_beta    = vec (num_params);
	da	    = vec (num_params);
	covar	    = matr (num_params, num_params);
	tmp_alpha   = matr (num_params, num_params);

	*lambda  = -*lambda;	  /* make lambda positive */
	return analyze (a, alpha, beta, chisq, varianz) ? OK : ERROR;
    }
    old_chisq = *chisq;

    for ( j=0 ; j<num_params ; j++ ) {
	memcpy (covar[j], alpha[j], num_params*sizeof(double));
	covar[j][j] *= 1 + *lambda;
	one_da[j][0] = beta [j];
    }

    solve (covar, num_params, one_da, 1);	/* Equation solution */

    for ( j=0 ; j<num_params ; j++ )
	da[j] = one_da[j][0];			/* changes in paramss */

    /* once converged, free dynamic allocated vars */

    if ( *lambda == 0.0 ) {
	free (beta);
	free (tmp_beta);
	free (da);
	free (temp_a);
	free_matr (one_da, num_params);
	free_matr (tmp_alpha, num_params);
	free_matr (covar, num_params);
        return OK;
    }

    /* check if trial did ameliorate sum of squares */

    for ( j=0 ; j<num_params ; j++ ) {
	temp_a[j] = a[j] + da[j];
	tmp_beta[j] = beta [j];
	memcpy (tmp_alpha[j], alpha[j], num_params*sizeof(double));
    }

    if ( !analyze (temp_a, tmp_alpha, tmp_beta, chisq, varianz) )
	return ERROR;

    if ( *chisq < old_chisq ) {     /* Success, accept new solution */
	if ( *lambda > 1e-20 )
	    *lambda /= 10;
	old_chisq = *chisq;
	for ( j=0 ; j<num_params ; j++ ) {
	    memcpy (alpha[j], tmp_alpha[j], num_params*sizeof(double));
	    beta[j] = tmp_beta[j];
	    a[j] = temp_a[j];
	}
	return BETTER;
    }
    else {			    /* failure, increase lambda and return */
	*lambda *= 10;
	*chisq = old_chisq;
	return WORSE;
    }
}


/*****************************************************************
    compute chi-square and numeric derivations
*****************************************************************/
static boolean analyze (a, alpha, beta, chisq, varianz)
double a[];
double **alpha;
double beta[];
double *chisq;
double *varianz;
{

/*
 *  used by marquardt to evaluate the linearized fitting matrix alpha
 *  and vector beta
 */

    word    k, j, i;
    double  wt, sig2i, dy, **dyda, *yfunc;
#if !defined(ATARI) && !defined(MTOS)
    double *edp, *yp, *fp, **dr, **ar, *dc, *bp, *ac;
#endif

    yfunc = vec (num_data);
    dyda = matr (num_data, num_params);

    /* initialize alpha beta */
    for ( j=0 ; j<num_params ; j++ ) {
	for ( k=0 ; k<=j ; k++ )
	    alpha[j][k] = 0;
	beta[j] = 0;
    }

    *chisq = *varianz = 0;
    calculate (yfunc, dyda, a);

#if !defined(ATARI) && !defined(MTOS)
    edp = err_data;
    yp = fit_y;
    fp = yfunc;
    dr = dyda;

    /* Summation loop over all data */
    for ( i=0 ; i<num_data ; i++, dr++ ) {
	/* The original code read: sig2i = 1/(*edp * *edp++); */
	/* There are some compilers that evaluate the operation */
	/* from right to left, although it is an error to do so. */
	/* Hence the following modification: */
	sig2i = 1/(*edp * *edp);
	edp++;
	dy = *yp++ - *fp++;
	*varianz += dy*dy;
	ar = alpha;
        dc = *dr;
	bp = beta;
	for ( j=0 ; j<num_params ; j++ ) {
	    wt = *dc++ * sig2i;
	    ac = *ar++;
            for ( k=0 ; k<=j ; k++ )
		*ac++ += wt * (*dr)[k];
	    *bp++ += dy * wt;
	}
	*chisq += dy*dy*sig2i;		    /* and find chi-square */
    }
#else
     /* Summation loop over all data */
    for ( i=0 ; i<num_data ; i++) {
        sig2i = 1/(err_data[i] * err_data[i]);
        dy = fit_y[i] - yfunc[i];
	*varianz += dy*dy;
        for ( j=0 ; j<num_params ; j++ ) {
	    wt = dyda[i][j] * sig2i;
            for ( k=0 ; k<=j ; k++ )
		alpha[j][k] += wt * dyda[i][k];
            beta[j] += dy * wt;
	}
	*chisq += dy*dy*sig2i;		    /* and find chi-square */
    }
#endif
    *varianz /= num_data;
    for ( j=1 ; j<num_params ; j++ )	    /* fill in the symmetric side */
	for ( k=0 ; k<=j-1 ; k++ )
	    alpha[k][j] = alpha [j][k];
    free (yfunc);
    free_matr (dyda, num_data);
    return TRUE;
}


/*****************************************************************
    compute function values and partial derivatives of chi-square
*****************************************************************/
static void calculate (yfunc, dyda, a)
double *yfunc;
double **dyda;
double a[];
{
    word    k, p;
    double  tmp_a;
    double  *tmp_low,
	    *tmp_high,
	    *tmp_pars;

    tmp_low = vec (num_data);	      /* numeric derivations */
    tmp_high = vec (num_data);
    tmp_pars = vec (num_params);

    /* first function values */

    call_gnuplot (a, yfunc);

    /* then derivatives */

    for ( p=0 ; p<num_params ; p++ )
	tmp_pars[p] = a[p];
    for ( p=0 ; p<num_params ; p++ ) {
	tmp_pars[p] = tmp_a = fabs(a[p]) < NEARLY_ZERO ? NEARLY_ZERO : a[p];
/*
 *  the more exact method costs double execution time and is therefore
 *  commented out here. Change if you like!
 */
/*	  tmp_pars[p] *= 1-DELTA;
	call_gnuplot (tmp_pars, tmp_low);  */

	tmp_pars[p] = tmp_a * (1+DELTA);
	call_gnuplot (tmp_pars, tmp_high);
	for ( k=0 ; k<num_data ; k++ )

/*            dyda[k][p] = (tmp_high[k] - tmp_low[k])/(2*tmp_a*DELTA); */

	    dyda[k][p] = (tmp_high[k] - yfunc[k])/(tmp_a*DELTA);
        tmp_pars[p] = a[p];
    }

    free (tmp_low);
    free (tmp_high);
    free (tmp_pars);
}


/*****************************************************************
    call internal gnuplot functions
*****************************************************************/
static void call_gnuplot (par, data)
double *par;
double *data;
{
    word    i;
    struct value v;

    /* set parameters first */
    for ( i=0 ; i<num_params ; i++ ) {
	Gcomplex (&v, par[i], 0.0);
	setvar (par_name[i], v);
    }

    for (i = 0; i < num_data; i++) {
	/* fit_index for multi-branch functions */
	(void) Ginteger (&func.dummy_values[0], (i+1));
	setvar (fit_index, func.dummy_values[0]);
	/* calculate fit-function value */
	(void) Gcomplex (&func.dummy_values[0], fit_x[i], 0.0);
	evaluate_at (func.at,&v);
	if ( undefined )
	    Eex ("Undefined value during function evaluation")
	data[i] = real(&v);
    }
    (void) Ginteger (&func.dummy_values[0], 0);
    setvar (fit_index, func.dummy_values[0]);
}


/*****************************************************************
    handle user interrupts during fit
*****************************************************************/
static boolean fit_interrupt ()
{
    char c;

    while ( TRUE ) {
	fprintf (STANDARD, "\n\n(S)top fit, (C)ontinue, (E)xecute:  ");
        /* toupper is a macro that might call fgetc twice on the Amiga */
        c = fgetc (stdin);
	switch (toupper(c)) {
	case 'S' :  fprintf (STANDARD, "Stop.");
		    return FALSE;
	case 'C' :  fprintf (STANDARD, "Continue.");
		    return TRUE;
	case 'E' :  {
			int i;
			struct value v;
			char *tmp;

			tmp = *fit_script ? fit_script : DEFAULT_CMD;
			fprintf (STANDARD, "executing: %s", tmp);
			/* set parameters visible to gnuplot */
			for ( i=0 ; i<num_params ; i++ ) {
			    Gcomplex (&v, a[i], 0.0);
			    setvar (par_name[i], v);
			}
			sprintf (input_line, tmp);
			(void) do_line ();
		    }
	}
    }
}


/*****************************************************************
    frame routine for the marquardt-fit
*****************************************************************/
static boolean regress (a)
double a[];
{
    double	**covar,
		**correl,
		*dpar,
		**alpha,
		chisq,
		last_chisq,
		varianz,
		lambda;
    word	i, j;
    marq_res_t	res;
    struct value val;

    chisq   = last_chisq = INFINITY;
    alpha   = matr (num_params, num_params);
    lambda  = -START_LAMBDA;			/* use sign as flag */
    i = 0;					/* iteration counter  */

    /* ctrlc now serves as Hotkey */
    ctrlc_setup ();

    /* Initialize internal variables and 1st chi-square check */
    if ( (res = marquardt (a, alpha, &chisq, &lambda, &varianz)) == ERROR )
	Eex ("FIT: error occured during fit")
    res = BETTER;

    Dblf ("\nInitial set of free parameters:\n")
    show_fit (i, chisq, chisq, a, lambda, STANDARD);
    show_fit (i, chisq, chisq, a, lambda, log_f);

    /* MAIN FIT LOOP: do the regression iteration */

    do {
/*
 *  MSDOS defines signal(SIGINT) but doesn't handle it through
 *  real interrupts. So there remain cases in which a ctrl-c may
 *  be uncaught by signal. We must use kbhit() instead that really
 *  serves the keyboard interrupt (or write an own interrupt func
 *  which also generates #ifdefs)
 *
 *  I hope that other OSes do it better, if not... add #ifdefs :-(
 *  EMX does not have kbhit.
 */
#if ((defined(MSDOS) || defined(DOS386)) && !defined(__EMX__)) || defined(ATARI) || defined(MTOS)
 	if ( kbhit () ) {
	    do {getchx ();} while ( kbhit() );
	    ctrlc_flag = TRUE;
	}
#endif

        if ( ctrlc_flag ) {
	    show_fit (i, chisq, last_chisq, a, lambda, STANDARD);
	    ctrlc_flag = FALSE;
	    if ( !fit_interrupt () )	       /* handle keys */
		break;
	}
	if ( res == BETTER ) {
            i++;
            last_chisq = chisq;
	}
        if ( (res = marquardt (a, alpha, &chisq, &lambda, &varianz)) == BETTER )
	    show_fit (i, chisq, last_chisq, a, lambda, STANDARD);
    }
    while ( res != ERROR  &&  lambda < MAX_LAMBDA  &&
	    (res == WORSE  ||  ( chisq > NEARLY_ZERO ? (last_chisq-chisq)/chisq : (last_chisq-chisq)) > epsilon) );

    /* restore original SIGINT function */
#if !defined(ATARI) && !defined(MTOS)
    interrupt_setup ();
#endif
    Dblf2 ("\nAfter %d iterations the fit converged.\n", i)
    Dblf2 ("final sum of squares residuals    : %g\n", chisq)
    if (chisq > NEARLY_ZERO) {
	    Dblf2 ("rel. change during last iteration : %g\n\n", (chisq-last_chisq)/chisq)
	 } else {
	    Dblf2 ("abs. change during last iteration : %g\n\n", (chisq-last_chisq))
	 }

    if ( res == ERROR )
	Eex ("FIT: error occured during fit")

    /* compute errors in the parameters */

    if (varianz > NEARLY_ZERO) { /* what should we do if varianz=0 (exact fit) ? */
    for ( i=0 ; i<num_params; i++ )
	for ( j=0 ; j<num_params; j++ )
	    alpha[i][j] /= varianz;
	 }
	 
    /* get covariance-, Korrelations- and Kurvature-Matrix */
    /* and errors in the parameters			  */

    covar   = matr (num_params, num_params);
    correl  = matr (num_params, num_params);
    dpar = vec (num_params);

    inverse (alpha, covar, num_params);

    Dblf ("Final set of parameters            68.3%% confidence interval (one at a time)\n")
    Dblf ("=======================            =========================================\n\n")
    for ( i=0 ; i<num_params; i++ ) {
	dpar[i] = sqrt (covar[i][i]);
	Dblf5 ("%-15.15s = %-15g  %-3.3s %g\n", par_name[i], a[i], PLUSMINUS, dpar[i])
    }

    Dblf ("\n\ncorrelation matrix of the fit parameters:\n\n")
    Dblf ("               ")
    for ( j=0 ; j<num_params ; j++ )
	Dblf2 ("%-6.6s ", par_name[j])
    Dblf ("\n")

    for ( i=0 ; i<num_params; i++ ) {
	Dblf2 ("%-15.15s", par_name[i])
	for ( j=0 ; j<num_params; j++ ) {
	    correl[i][j] = covar[i][j] / (dpar[i]*dpar[j]);
	    Dblf2 ("%6.3f ", correl[i][j])
        }
	Dblf ("\n")
    }

    /* restore last parameter's value (not done by calculate) */
    Gcomplex (&val, a[num_params-1], 0.0);
    setvar (par_name[num_params-1], val);

    /* call destruktor for allocated vars */
    lambda = 0;
    marquardt (a, alpha, &chisq, &lambda, &varianz);

    free_matr (covar, num_params);
    free_matr (alpha, num_params);
    free_matr (correl, num_params);
    free (dpar);

    return TRUE;
}


/*****************************************************************
    display actual state of the fit
*****************************************************************/
static void show_fit (i, chisq, last_chisq, a, lambda, device)
int i;
double chisq;
double last_chisq;
double *a;
double lambda;
FILE *device;
{
    int k;

    fprintf (device, "\n\n\
Iteration %d\n\
chisquare : %-15g   relative deltachi2 : %g\n\
deltachi2 : %-15g   limit for stopping : %g\n\
lambda	  : %g\n\nactual set of parameters\n\n",
i, chisq, chisq > NEARLY_ZERO ? (chisq - last_chisq)/chisq : 0.0, chisq - last_chisq, epsilon, lambda);
    for ( k=0 ; k<num_params ; k++ )
	fprintf (device, "%-15.15s = %g\n", par_name[k], a[k]);
}



/*****************************************************************
    is_empty: check for valid string entries
*****************************************************************/
static boolean is_empty (s)
char *s;
{
    while ( *s == ' '  ||  *s == '\t'  ||  *s == '\n' )
	s++;
    return (boolean)( *s == '#'  ||  *s == '\0' );
}


/*****************************************************************
    get next word of a multi-word string, advance pointer
*****************************************************************/
char *get_next_word (s, subst)
char **s;
char *subst;
{
    char *tmp = *s;

    while ( *tmp==' ' || *tmp=='\t' || *tmp=='=' )
	tmp++;
    if ( *tmp=='\n' || *tmp=='\0' )                    /* not found */
	return NULL;
    if ( (*s = strpbrk (tmp, " =\t\n")) == NULL )
	*s = tmp + strlen(tmp);
    *subst = **s;
    *(*s)++ = '\0';
    return tmp;
}



/*****************************************************************
    get valid numerical entries
*****************************************************************/
#if 0 /* not used */
static int scan_num_entries (s, vlist)
char *s;
double vlist[];
{
    int num = 0;
    char c, *tmp;

    while ( (tmp = get_next_word (&s, &c)) != NULL
					&&  num <= MAX_VALUES_PER_LINE )
	if ( !sscanf (tmp, "%lf", &vlist[++num]) )
	    Eex ("syntax error in data file")
    return num;
}
#endif

/*****************************************************************
    check for variable identifiers
*****************************************************************/
static boolean is_variable (s)
char *s;
{
    while ( *s ) {
	if ( !isalnum(*s) && *s!='_' )
	    return FALSE;
	s++;
    }
    return TRUE;
}


/*****************************************************************
    strcpy but max n chars
*****************************************************************/
static void copy_max (d, s, n)
char *d;
char *s;
int n;
{
    strncpy (d, s, n);
    if ( strlen(s) >= n )
	d[n-1] = '\0';
}


/*****************************************************************
    portable implementation of strnicmp (hopefully)
*****************************************************************/

#ifdef NEED_STRNICMP
int strnicmp (s1, s2, n)
char *s1;
char *s2;
int n;
{
    char c1,c2;

    if(n==0) return 0;

    do {
	c1 = *s1++; if(islower(c1)) c1=toupper(c1);
	c2 = *s2++; if(islower(c2)) c2=toupper(c2);
    } while(c1==c2 && c1 && c2 && --n>0);

    if(n==0 || c1==c2) return 0;
    if(c1=='\0' || c1<c2) return 1;
    return -1;
}
#endif


/*****************************************************************
    first time settings
*****************************************************************/
void init_fit ()
{
    struct value val;

    /* index used for multi-branch functions */
    fit_index = fit_index_def;
    val.type = INTGR;
    val.v.int_val = 0;
    setvar (fit_index, val);
    func.at = (struct at_type *) NULL;      /* need to parse 1 time */
}


/*****************************************************************
	    Set a GNUPLOT user-defined variable
******************************************************************/

void setvar (varname, data)
char *varname;
struct value data;
{
    register struct udvt_entry *udv_ptr = first_udv,
			       *last = first_udv;

    /* check if it's already in the table... */

    while (udv_ptr) {
	    last = udv_ptr;
	    if ( !strcmp (varname, udv_ptr->udv_name))
		break;
	    udv_ptr = udv_ptr->next_udv;
    }

    if ( !udv_ptr ) {			    /* generate new entry */
	last->next_udv = udv_ptr = (struct udvt_entry *)
	  alloc ((unsigned int)sizeof(struct udvt_entry), "setvar");
	udv_ptr->next_udv = NULL;
    }
    copy_max (udv_ptr->udv_name, varname, sizeof(udv_ptr->udv_name));
    udv_ptr->udv_value = data;
    udv_ptr->udv_undef = FALSE;
}



/*****************************************************************
    Read INTGR Variable value, return 0 if undefined or wrong type
*****************************************************************/
int getivar (varname)
char *varname;
{
    register struct udvt_entry *udv_ptr = first_udv;

    while (udv_ptr) {
	    if ( !strcmp (varname, udv_ptr->udv_name))
		return udv_ptr->udv_value.type == INTGR
		       ? udv_ptr->udv_value.v.int_val	/* valid */
		       : 0;				/* wrong type */
            udv_ptr = udv_ptr->next_udv;
    }
    return 0;						/* not in table */
}


/*****************************************************************
    Read DOUBLE Variable value, return 0 if undefined or wrong type
   I dont think it's a problem that it's an integer - div
*****************************************************************/
static double getdvar (varname)
char *varname;
{
    register struct udvt_entry *udv_ptr = first_udv;

    for ( ; udv_ptr ; udv_ptr=udv_ptr->next_udv)
	    if ( strcmp (varname, udv_ptr->udv_name)==0)
	      return real(&(udv_ptr->udv_value));

    /* get here => not found */
    return 0;
}

/* like getdvar, but
 * - convert it from integer to real if necessary
 * - create it with value INITIAL_VALUE if not found or undefined
 */
static double createdvar(varname,value)
char *varname;
double value;
{
    register struct udvt_entry *udv_ptr = first_udv;

    for ( ; udv_ptr ; udv_ptr=udv_ptr->next_udv)
	    if ( strcmp (varname, udv_ptr->udv_name)==0) {
	    	if (udv_ptr->udv_undef) {
	    		udv_ptr->udv_undef=0;
	    		Gcomplex(&udv_ptr->udv_value, value, 0.0);
	    	} else if (udv_ptr->udv_value.type==INTGR) {
	    		Gcomplex(&udv_ptr->udv_value, (double)udv_ptr->udv_value.v.int_val, 0.0);
	    	}
		   return real(&(udv_ptr->udv_value));
		 }

    /* get here => not found */

    {
        struct value a;
        Gcomplex(&a, value, 0.0);
        setvar(varname, a);
    }
    
    return value;
}
	

/*****************************************************************
    Split Identifier into path and filename
*****************************************************************/
static void splitpath (s, p, f)
char *s;
char *p;
char *f;
{
    register char *tmp;
    tmp = s + strlen(s) - 1;
    while ( *tmp != '\\'  &&  *tmp != '/'  &&  *tmp != ':'  &&  tmp-s >= 0 )
	tmp--;
    strcpy (f, tmp+1);
    strncpy (p, s, tmp-s+1);
    p[tmp-s+1] = '\0';
}


#if defined(MSDOS) || defined(DOS386)
/*****************************************************************
    Character substitution
*****************************************************************/
#ifdef ANSI_C
static void subst (char *s, char from, char to)
#else
static void subst (s, from, to)
char *s;
char from;
char to;
#endif
{
    while ( *s ) {
	if ( *s == from )
	    *s = to;
	s++;
    }
}
#endif


/*****************************************************************
    write the actual parameters to start parameter file
*****************************************************************/
void update (pfile)
char *pfile;
{
    char	fnam[256],
		path[256],
		tmpname[256],
		sstr[256],
		pname[64],
		tail[127],
		*s = sstr,
#if !defined(AMIGA_SC_6_1)
		*tmpn,
#endif
		*tmp,
		c;
    FILE        *of,
		*nf;
    double	pval;

    /* split into path and filename */

#if defined(MSDOS) || defined(DOS386)
    /* don't get problems during system() call */
    subst (pfile, '/', '\\');
#endif
    splitpath (pfile, path, fnam);
    if ( !(of = fopen (pfile, "r")) )
	Eex2 ("parameter file %s could not be read", pfile)
#if defined(AMIGA_SC_6_1)
    sprintf (tmpname, "%sgf%08lx", path, FindTask (NULL));
#else /* !AMIGA_SC_6_1 */
    /* Under MSDOS rename will not work properly if TMP points to another
     * drive/dir. In this case use 'path' directory for tempfile
     */
#if defined (MSDOS)  ||  defined (OS2)
    sprintf (tmpname, "%s%s", path, TEMPNAME);
#else
#ifdef HAVE_TEMPNAM
    tmpn = tempnam (path, "gnuft");
#else
    tmpn = tmpnam(NULL);
#endif

    if ( tmpn == NULL )
	Eex ("no more unique names possible during updating")
    strcpy (tmpname, tmpn);
    free (tmpn);
#endif
#endif /* !AMIGA_SC_6_1 */
    if ( !(nf = fopen (tmpname, "w")) )
	Eex ("could not open temporary file")
    while ( TRUE ) {
	if ( !fgets (s = sstr, sizeof(sstr), of) )	/* EOF found */
	    break;
	if ( is_empty(s) ) {
	    fprintf (nf, s);				/* preserve comments */
            continue;
        }
        if ( (tmp = strchr (s, '#')) != NULL ) {
	    copy_max (tail, tmp, sizeof(tail));
	    *tmp = '\0';
	}
	else
	    strcpy (tail, "\n");
	tmp = get_next_word (&s, &c);
	if ( !is_variable (tmp)  ||  strlen(tmp) > MAX_VARLEN ) {
	    fclose (of);
	    Eex2 ("syntax error in parameter file %s", fnam)
	}
	copy_max (pname, tmp, sizeof(pname));
	/* next must be '=' */
	if ( c != '=' ) {
	    tmp = strchr (s, '=');
	    if ( tmp == NULL ) {
		fclose (of);
		Eex2 ("syntax error in parameter file %s", fnam)
	    }
	    s = tmp+1;
	}
	tmp = get_next_word (&s, &c);
	if ( !sscanf (tmp, "%lf", &pval) ) {
	    fclose (of);
	    Eex2 ("syntax error in parameter file %s", fnam)
	}
	if ( (tmp = get_next_word (&s, &c)) != NULL ) {
	    fclose (of);
	    Eex2 ("syntax error in parameter file %s", fnam)
	}

        /* now modify */

	if ( (pval = getdvar (pname)) == 0 )
	    pval = (double) getivar (pname);

	sprintf (sstr, "%g", pval);
	if ( !strchr (sstr, '.')  &&  !strchr (sstr, 'e') )
	    strcat (sstr, ".0");                /* assure CMPLX-type */
	fprintf (nf, "%-15.15s = %-15.15s   %s", pname, sstr, tail);
    }
    if ( fclose (of)  ||  fclose (nf) )
	Eex ("I/O error during update")

#ifdef VMS
    if ( delete (pfile) )     /* implies write permission !!!! */
#else	/* VMS */
    if ( unlink (pfile) )     /* implies write permission !!!! */
#endif	/* VMS */
      Eex2 ("parameter file %s could not be deleted", pfile)

#if defined(MSDOS) || defined(OS2)
    sprintf (sstr, "ren %s %s", tmpname, fnam);
#ifdef _Windows
    if ( winsystem (sstr) )
#else
    if ( system (sstr) )
#endif
#else                   /* implies write permission !!!! */
    if ( rename (tmpname, pfile) )
#endif
      Eex3 ("unable to rename %s to %s", tmpname, pfile)
}



/*****************************************************************
    Interface to the classic gnuplot-software
*****************************************************************/

void do_fit ()
{
    int autorange=3;  /* yes */
    double min, max;  /* range to fit */
    int dummy_token=-1; /* local variable name */

	int 	i;
	int columns;
	double	v[3];
	double tmpd;
	time_t	timer;
	int token1, token2, token3;
	char *tmp;
	extern struct udft_entry *dummy_func;
	extern char dummy_var[MAX_NUM_VAR][MAX_ID_LEN+1];
	extern char c_dummy_var[MAX_NUM_VAR][MAX_ID_LEN+1];
	extern int c_token;
   extern int df_line_number;
   
/* first look for a restricted fit range */

    if (equals(c_token, "[")) {
	c_token++;
	if (isletter(c_token)) {
	    if (equals(c_token + 1, "=")) {
		dummy_token = c_token;
		c_token += 2;
	    }
			/* else parse it as a xmin expression */
	}
	autorange = load_range(FIRST_X_AXIS,&min, &max, autorange);
	if (!equals(c_token, "]"))
		int_error("']' expected", c_token);
	c_token++;
    }


/* now compile the function */

    token1 = c_token;

    if (func.at) {
	free(func.at);
	func.at=NULL;  /* in case perm_at() does int_error */
    }

    dummy_func = &func;

    if (dummy_token >= 0)
	copy_str(c_dummy_var[0], dummy_token, MAX_ID_LEN);
    else
	strcpy(c_dummy_var[0], dummy_var[0]);
	
    func.at = perm_at();
    dummy_func = NULL;

    token2 = c_token;

/* use datafile module to parse the datafile and qualifiers */

    columns = df_open(3);  /* up to 3 using specs allowed */
    if (columns==1)
	int_error("Need 2 or 3 using specs", c_token);
    /* defer actually reading the data until we have parsed the rest of the line */
    
	token3 = c_token;

    tmpd = getdvar (FITLIMIT);		  /* get epsilon if given explicitly */
    if ( tmpd < 1.0  &&  tmpd > 0.0 )
        epsilon = tmpd;

    *fit_script = '\0';                   /* FIT_SKIP 1 means each 2nd value */
    if ( (tmp = getenv (FITSCRIPT)) != NULL )
	strcpy (fit_script, tmp);

    tmp = getenv (GNUFITLOG);		  /* open logfile */
    if ( tmp != NULL ) {
        char *tmp2 = &tmp[strlen(tmp)-1];
        if ( *tmp2 == '/'  ||  *tmp2 == '\\' )
            sprintf (logfile, "%s%s", tmp, logfile);
        else
            strcpy (logfile, tmp);
    }
    if ( !log_f && /* div */ !(log_f = fopen (logfile, "a")) )
	Eex2 ("could not open log-file %s", logfile)
    fprintf (log_f, "\n\n*******************************************************************************\n");
    time (&timer);
    fprintf (log_f, "%s\n\n", ctime (&timer));
	{	char line[MAX_LINE_LEN];
		capture(line, token2,token3-1,MAX_LINE_LEN);
			fprintf (log_f, "FIT:    data read from %s\n", line);
	}

    fit_x = vec (MAX_DATA);		 /* start with max. value */
    fit_y = vec (MAX_DATA);

/* first read in experimental data */

    err_data = vec (MAX_DATA);
    num_data = 0;

	while ( (i=df_readline(v, 3)) != EOF ) {

		if (i==0)
			continue;

		if ( !(autorange&1) && (v[0] < min))
			continue;  /* skip this one - div */
		if ( !(autorange&2) && (v[0] > max))
			continue;  /* and again */

		if ( num_data==MAX_DATA ) {
			df_close();
			Eex2 ("max. # of datapoints %d exceeded", MAX_DATA);
		}

		switch (i)
		{
			case DF_UNDEFINED:
			case DF_FIRST_BLANK:
			case DF_SECOND_BLANK:
				continue;
			case 0:
				Eex2("bad data on line %d", df_line_number);
				break;
			case 1:
				v[1]=v[0];
				v[0]=df_line_number;
				break;
		}
				
		fit_x[num_data] = v[0];
		fit_y[num_data] = v[1];
		/* datafile module is such that if columns==3, it
		 * will only return lines with 3 data read
		 */
		err_data[num_data++] = (columns==3 ? v[2] : 1);
	}
	df_close();

    /* now resize fields */
    redim_vec (&fit_x, num_data);
    redim_vec (&fit_y, num_data);
    redim_vec (&err_data, num_data);
    
    fprintf (log_f, "        #datapoints = %d\n", num_data);
    if ( columns==2 )
	fprintf (log_f, "        y-errors assumed equally distributed\n\n");

	{	char line[MAX_LINE_LEN];
		capture(line, token1, token2-1,MAX_LINE_LEN);		
		fprintf (log_f, "function used for fitting: %s\n", line);
	}

    /* read in parameters */

    a = vec (MAX_PARAMS);
    if ((par_name = (fixstr *) malloc ((MAX_PARAMS+1)*sizeof(fixstr))) == NULL)
	Eex ("Memory running low during fit")
    num_params = 0;

	if (!equals(c_token, "via"))
		int_error("Need via and either parameter list or file", c_token);
	if (isstring(++c_token)) {
		boolean fixed;
		double	tmp_par;
		char c, *s;
		char sstr[MAX_LINE_LEN];
		FILE *f;

		quote_str(sstr, c_token++, MAX_LINE_LEN);

		/* get parameters and values out of file and ignore fixed ones */

		fprintf (log_f, "take parameters from file: %s\n\n", sstr);
		if ( !(f = fopen (sstr, "r")) )
			Eex2 ("could not read parameter-file %s", sstr);
		while ( TRUE ) {
			if ( !fgets (s = sstr, sizeof(sstr), f) )		/* EOF found */
				break;
			if ( (tmp = strstr (s, FIXED)) != NULL ) {	      /* ignore fixed params */
				*tmp = '\0';
				fprintf (log_f, "FIXED:  %s\n", s);
				fixed = TRUE;
			}
			else
				fixed = FALSE;
			if ( (tmp = strchr (s, '#')) != NULL )
				*tmp = '\0';
			if ( is_empty(s) )
				continue;
			tmp = get_next_word (&s, &c);
			if ( !is_variable (tmp)  ||  strlen(tmp) > MAX_VARLEN ) {
				fclose (f);
				Eex ("syntax error in parameter file");
			}
			
			copy_max (par_name[num_params], tmp, sizeof(fixstr));
			/* next must be '=' */
			if ( c != '=' ) {
				tmp = strchr (s, '=');
				if ( tmp == NULL ) {
					fclose (f);
					Eex ("syntax error in parameter file");
				}
				s = tmp+1;
			}
			tmp = get_next_word (&s, &c);
			if ( sscanf (tmp, "%lf", &tmp_par) != 1 ) {
				fclose (f);
				Eex ("syntax error in parameter file");
			}
			
			/* make fixed params visible to GNUPLOT */
			if ( fixed ) {
				struct value v;
				Gcomplex (&v, tmp_par, 0.0);
				setvar (par_name[num_params], v);   /* use parname as temp */
			}
			else
				a[num_params++] = tmp_par;
			
			if ( (tmp = get_next_word (&s, &c)) != NULL ) {
				fclose (f);
				Eex ("syntax error in parameter file");
			}
		}
		fclose (f);
	} else { /* not a string after via */
		fprintf (log_f, "use actual parameter values\n\n");
		do {
			if (!isletter(c_token))
				Eex ("no parameter specified");
			capture(par_name[num_params], c_token, c_token, sizeof(par_name[0]));
			/* create variable if it doesn't exist */
			a[num_params] = createdvar (par_name[num_params], INITIAL_VALUE);
			num_params++;
		}	while (equals(++c_token, ",") && ++c_token);
	}
    redim_vec (&a, num_params);
    par_name = (fixstr *) realloc (par_name, (num_params+1)*sizeof(fixstr));

    /* avoid parameters being equal to zero */
    for ( i=0 ; i<num_params ; i++ )
	if ( a[i] == 0 )
	    a[i] = NEARLY_ZERO;

    regress (a);

    fclose (log_f);
    log_f = NULL;
    free (fit_x);
    free (fit_y);
    free (err_data);
    free (a);
    free (par_name);
    if ( func.at ) {
        free (func.at);                     /* release perm. action table */
        func.at = (struct at_type *) NULL;
    }

}

#if defined(ATARI) || defined(MTOS)
int kbhit()
{
	fd_set rfds;
	struct timeval timeout;
	
	timeout.tv_sec = 0;
	timeout.tv_usec = 0;
	rfds = (fd_set)(1L << fileno(stdin));
	return((select(0, &rfds, NULL, NULL, &timeout) > 0) ? 1 : 0);
}
#endif