#ifndef lint
static char *RCSid = "$Id: fit.c 1.3 1993/08/03 07:14:17 cg Exp cg $";
#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)
 */


#define FIT_MAIN

#include <stdlib.h>
#include <stdio.h>

#include <string.h>
#include <math.h>
#include <stdarg.h>
#include <time.h>
#include <ctype.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
#include <fcntl.h>
#include <termio.h>
#endif
#endif

#include <setjmp.h>
#define STANDARD    stderr	/* compatible with gnuplot philosophy */

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

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


#define INFINITY    1e30
#define NEARLY_ZERO 1e-30
#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

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

static char fitfunction[80];
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 char *TEMPNAME	    = "gfttemp.$$$";


static FILE	    *log_f = NULL;

static word	    num_data,
		    num_params,
		    skip;
static double	    *fit_x,
		    *fit_y,
		    *err_data,
		    *a;

static struct udft_entry func;

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

extern jmp_buf env;			    /* from plot.c */
extern char input_line[];		    /* command.c */

#if !defined(MSDOS) && !defined(DOS386) && !defined(AMIGA_SC_6_1)
static	char buff = 0;	/* Unix, Linux, OS/2: buffer for getch(), kbhit() */
#ifndef OS2
#define getchx	getch	/* Unix or Linux */
#endif
#endif

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

static marq_res_t marquardt (double a[], double **alpha, double *chisq,
			     double *lambda, double *varianz);
static boolean analyze (double a[], double **alpha, double beta[],
			   double *chisq, double *varianz);
static void calculate (double *yfunc, double **dyda, double a[]);
static void call_gnuplot (double *par, double *data);
static void show_fit (int i, double chisq, double last_chisq, double *a,
		      double lambda, FILE *device);
static boolean regress (double a[]);
static int scan_num_entries (char *s, double vlist[]);
static boolean is_variable (char *s);
static void copy_max (char *d, char *s, int n);
static double getdvar (char *varname);
static void splitpath (char *s, char *p, char *f);
static void subst (char *s, char from, char to);
static void dbl_fprintf (FILE *f1, FILE *f2, char *s, ...);
static boolean fit_interrupt ();
static boolean is_empty (char *s);


extern struct value *Gcomplex(struct value *v, double re, double im);
extern struct value *Ginteger(struct value *a, int i);
extern double real(struct value *val);
extern void evaluate_at(struct at_type *at_ptr, struct value *val_ptr);
extern int do_line ();

#if !defined(MSDOS) && !defined(DOS386) && !defined(AMIGA_SC_6_1)
/* non-blocking IO stuff */
static int kbhit();
static int getch();
#endif

/*****************************************************************
    in case of fatal errors
*****************************************************************/
void error_ex (char *s, ...)
{
    va_list p;
    char    tmp[128],
	    *sp;

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

/* restore raw or cooked */
#if defined(AMIGA_SC_6_1)
    (void) rawcon (0);  /* Since we're returning to the command level, */
                        /* we might as well assume non-raw mode. */
#endif

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



/*****************************************************************
    fprintf to 2 streams at a time
*****************************************************************/
static void dbl_fprintf (FILE *f1, FILE *f2, char *s, ...)
{
    va_list p;

    va_start (p, s);
    vfprintf (f1, s, p);
    va_start (p, s);
    vfprintf (f2, s, p);
    va_end (p);
}


/*****************************************************************
    Marquardt's nonlinear least squares fit
*****************************************************************/
static marq_res_t marquardt (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 (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, *edp, *yp, *fp, **dr,
	    **ar, *dc, *bp, *ac;

    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);
    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 */
    }

    *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 (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 (double *par, double *data)
{
    word    i;
    struct value v;

    extern struct at_type *perm_at();
    extern int scanner (char *);

    /* all in command.c */
    extern char input_line[];
    extern int num_tokens, c_token;
    extern char dummy_var[MAX_NUM_VAR][MAX_ID_LEN+1];
    extern char c_dummy_var[MAX_NUM_VAR][MAX_ID_LEN+1];
    extern struct udft_entry *dummy_func;

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

    (void) strcpy (c_dummy_var[0], dummy_var[0]);

    if ( !func.at ) {		    /* parse fitfunction, build action table */
	strcpy (input_line, fitfunction);
	num_tokens = scanner (input_line);	/* divide into tokens */
	c_token = 0;
	dummy_func = &func;
	func.at = perm_at();			/* parse function */
    }
    for (i = 0; i < MAX_NUM_VAR; i++)
	c_dummy_var[i][0] = '\0';               /* no dummy variables */

    dummy_func = &func; 			/* by caution */
    for (i = 0; i < num_data; i++) {
	/* fit_index for multi-branch functions */
	(void) Ginteger (&func.dummy_values[0], (i+1)*skip);
	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);
	data[i] = real(&v);
    }
    (void) Ginteger (&func.dummy_values[0], 0);
    setvar (fit_index, func.dummy_values[0]);
}


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


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

    while ( TRUE ) {
	fprintf (STANDARD, "\n\n(S)top fit, (C)ontinue, (E)xecute:  ");
#if defined(AMIGA_SC_6_1)
	c = getch ();
#else /* !AMIGA_SC_6_1 */
	do {c = getchx ();} while ( kbhit() );
#endif /* !AMIGA_SC_6_1 */
        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 (double a[])
{
    double	**covar,
		**correl,
		*dpar,
		**alpha,
		chisq,
		last_chisq,
		varianz,
		lambda;
    word	i, j;
    int 	c;
    marq_res_t	res;
    struct value val;
#if defined(AMIGA_SC_6_1)
    BPTR StdinHandle;
#endif

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

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

    dbl_fprintf (STANDARD, log_f, "\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 {
#if defined(AMIGA_SC_6_1)
	StdinHandle = chkufb (fileno (stdin))->ufbfh;
	if (IsInteractive (StdinHandle) && 
	    WaitForChar (StdinHandle, 0) == DOSTRUE) {
	    fread (&c, 1, 1, stdin);
#else /* !AMIGA_SC_6_1 */
	if ( kbhit () ) {
	    do { getchx (); } while ( kbhit() );
#endif /* !AMIGA_SC_6_1 */
	    show_fit (i, chisq, last_chisq, a, lambda, STANDARD);
	    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  ||  (last_chisq-chisq)/chisq > epsilon) );

    dbl_fprintf (STANDARD, log_f, "\nAfter %d iterations the fit converged.\n", i);
    dbl_fprintf (STANDARD, log_f, "final sum of squares residuals    : %g\n", chisq);
    dbl_fprintf (STANDARD, log_f, "rel. change during last iteration : %g\n\n", (chisq-last_chisq)/chisq);

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

    /* compute errors in the parameters */

    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);

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

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

    for ( i=0 ; i<num_params; i++ ) {
	dbl_fprintf (STANDARD, log_f, "%-15.15s", par_name[i]);
	for ( j=0 ; j<num_params; j++ ) {
	    correl[i][j] = covar[i][j] / (dpar[i]*dpar[j]);
	    dbl_fprintf (STANDARD, log_f, "%6.3f ", correl[i][j]);
        }
	dbl_fprintf (STANDARD, log_f, "\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 (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 - last_chisq)/chisq, 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 (char *s)
{
    while ( *s == ' '  ||  *s == '\t'  ||  *s == '\n' )
	s++;
    return ( *s == '#'  ||  *s == '\0' );
}


/*****************************************************************
    get next word of a multi-word string, advance pointer
*****************************************************************/
char *get_next_word (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
*****************************************************************/
static int scan_num_entries (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, "%lg", &vlist[++num]) )
	    error_ex ("syntax error in data file");
    return num;
}


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


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


/*****************************************************************
    Slow but portable implementation of strnicmp
*****************************************************************/
/* 
 * please don't add further machine-specific macros, use -DSTRNICMP
 * as compiler switch if necessary
 */

#ifdef STRNICMP
int strnicmp (char *s1, char *s2, int n)
{
    static char ss1[256], ss2[256];
    register char *p;

    copy_max (ss1, s1, sizeof(ss1));
    copy_max (ss2, s2, sizeof(ss2));
    p = ss1-1;
    while ( *++p )
	toupper ( *p );
    p = ss2-1;
    while ( *++p )
        toupper ( *p );
    return strncmp (ss1, ss2, n);
}
#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
******************************************************************/
extern struct udvt_entry *first_udv;

void setvar (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 (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
*****************************************************************/
static double getdvar (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 == CMPLX
		       ? udv_ptr->udv_value.v.cmplx_val.real	/* valid */
		       : 0;				   /* wrong type */
            udv_ptr = udv_ptr->next_udv;
    }
    return 0;						 /* not in table */
}


/*****************************************************************
    Split Identifier into path and filename
*****************************************************************/
static void splitpath (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';
}


/*****************************************************************
    Character substitution
*****************************************************************/
static void subst (char *s, char from, char to)
{
    while ( *s ) {
	if ( *s == from )
	    *s = to;
	s++;
    }
}


/*****************************************************************
    write the actual parameters to start parameter file
*****************************************************************/
void update (char *pfile)
{
    char	fnam[256],
		path[256],
		tmpnam[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")) )
	error_ex ("parameter file %s could not be read", pfile);
#if defined(AMIGA_SC_6_1)
    sprintf (tmpnam, "%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 (tmpnam, "%s%s", path, TEMPNAME);
#else
    tmpn = tempnam (path, "gnuft");
    if ( tmpn == NULL )
        error_ex ("no more unique names possible during updating");
    strcpy (tmpnam, tmpn);
    free (tmpn);
#endif
#endif /* !AMIGA_SC_6_1 */
    if ( !(nf = fopen (tmpnam, "w")) )
	error_ex ("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, '#') ) {
	    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);
	    error_ex ("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);
		error_ex ("syntax error in parameter file %s", fnam);
	    }
	    s = tmp+1;
	}
	tmp = get_next_word (&s, &c);
	if ( !sscanf (tmp, "%lg", &pval) ) {
	    fclose (of);
	    error_ex ("syntax error in parameter file %s", fnam);
	}
	if ( (tmp = get_next_word (&s, &c)) != NULL ) {
	    fclose (of);
	    error_ex ("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) )
	error_ex ("I/O error during update");

    if ( unlink (pfile) )     /* implies write permission !!!! */
      error_ex ("parameter file %s could not be deleted", pfile);

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



/*****************************************************************
    Interface to the classic gnuplot-software
*****************************************************************/
void do_fit (char *fitfunc, char *datafile, int xcol, int ycol, int dycol,
	     char *paramsfile)
{
    FILE	*f;
    int 	i, num, num_rawdata;
    static char sstr[1024];
    char	c,
		*s = sstr,
		*tmp;
    double	vlist[MAX_VALUES_PER_LINE],
		tmpd;
    time_t	timer;


/* Amiga: force terminal to raw mode to handle user interrupts */
#if defined(AMIGA_SC_6_1)
    (void) rawcon (1);
#endif

    strcpy (fitfunction, fitfunc);
    if ( !(f = fopen (datafile, "r")) )
	error_ex ("datafile %s could not be read", datafile);

    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) )
	strcpy (fit_script, tmp);
    skip = 1 + (word) getivar (FIT_SKIP);

    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 = fopen (logfile, "a")) )
        error_ex ("could not open log-file %s", logfile);
    fprintf (log_f, "\n\n*******************************************************************************\n");
    time (&timer);
    fprintf (log_f, "%s\n\n", ctime (&timer));
    fprintf (log_f, "FIT:    data from %s\n        x = column #%d\n        y = column #%d\n",
                datafile, xcol, ycol);

    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_rawdata = num_data = 0;
    while ( TRUE ) {
	if ( ! fgets (s = sstr, sizeof(sstr), f) )
	    break;				    /* EOF found */
	if ( is_empty(s) )
	    continue;
	if ( num_data==MAX_DATA ) {
	    fclose (f);
	    error_ex ("max. # of datapoints %d exceeded", MAX_DATA);
	}
	/* expecting n value entries (x, y (, yerror) ) */
	num = scan_num_entries (s, vlist);
	if ( xcol>num || ycol>num || dycol>num ) {
	    fclose (f);
	    error_ex ("columns in %s do not match specification", datafile);
	}
	if ( ++num_rawdata % skip )		    /* ignore this one */
	    continue;
	fit_x[num_data] = vlist[xcol];
	fit_y[num_data]   = vlist[ycol];
	err_data[num_data++] = dycol ? vlist[dycol] : 1;
    }
    fclose (f);
    /* 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", num_rawdata);
    if ( skip > 1 )
	fprintf (log_f, "   use each %d. value\n", skip);
    else
	fprintf (log_f, "\n");
    if ( dycol )
	fprintf (log_f, "        y-error = column #%d\n\n", dycol);
    else
	fprintf (log_f, "        y-errors assumed equally distributed\n\n");
    fprintf (log_f, "function used for fitting: %s\n", fitfunction);

    /* read in parameters */

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

    if ( !strnicmp (paramsfile, "via ", 4) ) {   /* no values given */
	char *tmp2;

	fprintf (log_f, "use actual parameter values\n\n");
	paramsfile += 4;
	do {
	    tmp = strchr (paramsfile, ',');
	    if ( tmp != NULL )
		*tmp = '\0';
	    tmp2 = get_next_word (&paramsfile, &c);
	    if ( tmp2 == NULL )
		error_ex ("no parameter specified");
	    copy_max (par_name[num_params], tmp2, sizeof(fixstr));
	    a[num_params] = getdvar (par_name[num_params]);
	    num_params++;
	    paramsfile = tmp+1;
	}
	while ( tmp );
    }
    else {
	boolean fixed;
	double	tmp_par;

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

	fprintf (log_f, "take parameters from file: %s\n\n", paramsfile);
	if ( !(f = fopen (paramsfile, "r")) )
	    error_ex ("could not read parameter-file %s", paramsfile);
	while ( TRUE ) {
	    if ( !fgets (s = sstr, sizeof(sstr), f) )		/* EOF found */
		break;
	    if ( tmp = strstr (s, FIXED) ) {	      /* ignore fixed params */
		*tmp = '\0';
		fprintf (log_f, "FIXED:  %s\n", s);
		fixed = TRUE;
	    }
	    else
		fixed = FALSE;
            if ( tmp = strchr (s, '#') )
		*tmp = '\0';
	    if ( is_empty(s) )
		continue;
	    tmp = get_next_word (&s, &c);
	    if ( !is_variable (tmp)  ||  strlen(tmp) > MAX_VARLEN ) {
		fclose (f);
		error_ex ("syntax error in parameter file %s", paramsfile);
	    }
	    copy_max (par_name[num_params], tmp, sizeof(fixstr));
	    /* next must be '=' */
	    if ( c != '=' ) {
		tmp = strchr (s, '=');
		if ( tmp == NULL ) {
		    fclose (f);
		    error_ex ("syntax error in parameter file %s", paramsfile);
		}
		s = tmp+1;
	    }
	    tmp = get_next_word (&s, &c);
	    if ( !sscanf (tmp, "%lg", &tmp_par) ) {
		fclose (f);
		error_ex ("syntax error in parameter file %s", paramsfile);
	    }

	    /* 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);
		error_ex ("syntax error in parameter file %s", paramsfile);
	    }
	}
	fclose (f);
    }
    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;
    }

/* restore raw or cooked */
#if defined(AMIGA_SC_6_1)
    (void) rawcon (0);
#endif
}


/*****************************************************************
    replace for kbhit() under UNIXish systems
*****************************************************************/
#if !defined(MSDOS) && !defined(DOS386) && !defined(AMIGA_SC_6_1)
int	kbhit()
{
   int r;
   struct termio tmo, tmn;

   if (buff)
      return 1;

   ioctl(fileno(stdin),TCGETA,&tmo);
   tmn = tmo;
   tmn.c_iflag = tmn.c_lflag = 0;
   tmn.c_cc[VMIN]  = 0;      /* minimum number of characters to be read */
   tmn.c_cc[VTIME] = 1;      /* time-out in tenths of a second */
   ioctl(fileno(stdin),TCSETA,&tmn);

   r = read(fileno(stdin),&buff,1);
   ioctl(fileno(stdin),TCSETA,&tmo);
   return (r > 0);
}

int	getch()
{
   char c;
   struct termio tmo, tmn;

   if (buff) {
      c = buff;
      buff = 0;
      return (c & 127);
   }

   ioctl(fileno(stdin),TCGETA,&tmo);
   tmn = tmo;
   tmn.c_iflag = tmn.c_lflag = 0;
   tmn.c_cc[VMIN]  = 1;      /* minimum number of characters to be read */
   tmn.c_cc[VTIME] = 0;      /* time-out in tenths of a second */
   ioctl(fileno(stdin),TCSETA,&tmn);

   while (read(fileno(stdin),&c,1) != 1)
     ;	/* do nothing, just wait */
   ioctl(fileno(stdin),TCSETA,&tmo);
   return (c & 127);
}
#endif