2 static char *RCSid() { return RCSid("$Id: fit.c,v 1.56.2.4 2009/02/05 17:17:25 sfeam Exp $"); }
5 /* NOTICE: Change of Copyright Status
7 * The author of this module, Carsten Grammes, has expressed in
8 * personal email that he has no more interest in this code, and
9 * doesn't claim any copyright. He has agreed to put this module
10 * into the public domain.
12 * Lars Hecking 15-02-1999
16 * Nonlinear least squares fit according to the
17 * Marquardt-Levenberg-algorithm
19 * added as Patch to Gnuplot (v3.2 and higher)
22 * 930726: Recoding of the Unix-like raw console I/O routines by:
23 * Michele Marziani (marziani@ferrara.infn.it)
24 * drd: start unitialised variables at 1 rather than NEARLY_ZERO
25 * (fit is more likely to converge if started from 1 than 1e-30 ?)
27 * HBB (broeker@physik.rwth-aachen.de) : fit didn't calculate the errors
28 * in the 'physically correct' (:-) way, if a third data column containing
29 * the errors (or 'uncertainties') of the input data was given. I think
30 * I've fixed that, but I'm not sure I really understood the M-L-algo well
31 * enough to get it right. I deduced my change from the final steps of the
32 * equivalent algorithm for the linear case, which is much easier to
33 * understand. (I also made some minor, mostly cosmetic changes)
35 * HBB (again): added error checking for negative covar[i][i] values and
36 * for too many parameters being specified.
38 * drd: allow 3d fitting. Data value is now called fit_z internally,
39 * ie a 2d fit is z vs x, and a 3d fit is z vs x and y.
41 * Lars Hecking : review update command, for VMS in particular, where
42 * it is not necessary to rename the old file.
44 * HBB, 971023: lifted fixed limit on number of datapoints, and number
64 #if defined(VA_START) && defined(STDC_HEADERS)
65 static void Dblfn __PROTO((const char *fmt, ...));
67 static void Dblfn __PROTO(());
75 #if defined(MSDOS) || defined(DOS386) /* non-blocking IO stuff */
77 # ifndef _Windows /* WIN16 does define MSDOS .... */
81 #else /* !(MSDOS || DOS386) */
85 #endif /* !(MSDOS || DOS386) */
87 #if defined(ATARI) || defined(MTOS)
88 # define getchx() Crawcin()
89 static int kbhit(void);
93 OK, ML_ERROR, BETTER, WORSE
95 typedef enum marq_res marq_res_t;
101 #define INFINITY 1e30
102 #define NEARLY_ZERO 1e-30
104 /* create new variables with this value (was NEARLY_ZERO) */
105 #define INITIAL_VALUE 1.0
107 /* Relative change for derivatives */
110 #define MAX_DATA 2048
111 #define MAX_PARAMS 32
112 #define MAX_LAMBDA 1e20
113 #define MIN_LAMBDA 1e-20
114 #define LAMBDA_UP_FACTOR 10
115 #define LAMBDA_DOWN_FACTOR 10
117 #if defined(MSDOS) || defined(OS2) || defined(DOS386)
118 # define PLUSMINUS "\xF1" /* plusminus sign */
120 # define PLUSMINUS "+/-"
123 #define LASTFITCMDLENGTH 511
125 /* externally visible variables: */
127 /* saved copy of last 'fit' command -- for output by "save" */
130 /* log-file for fit command */
131 char *fitlogfile = NULL;
133 #ifdef GP_FIT_ERRVARS
134 TBOOLEAN fit_errorvariables = FALSE;
135 #endif /* GP_FIT_ERRVARS */
137 /* private variables: */
140 static int max_params;
142 static double epsilon = 1e-5; /* convergence limit */
143 static int maxiter = 0;
145 static char fit_script[256];
147 /* HBB/H.Harders 20020927: log file name now changeable from inside
149 static const char fitlogfile_default[] = "fit.log";
150 static const char GNUFITLOG[] = "FIT_LOG";
152 static const char *GP_FIXED = "# FIXED";
153 static const char *FITLIMIT = "FIT_LIMIT";
154 static const char *FITSTARTLAMBDA = "FIT_START_LAMBDA";
155 static const char *FITLAMBDAFACTOR = "FIT_LAMBDA_FACTOR";
156 static const char *FITMAXITER = "FIT_MAXITER"; /* HBB 970304: maxiter patch */
157 static const char *FITSCRIPT = "FIT_SCRIPT";
158 static const char *DEFAULT_CMD = "replot"; /* if no fitscript spec. */
159 static char last_fit_command[LASTFITCMDLENGTH+1] = "";
161 static FILE *log_f = NULL;
163 static int num_data, num_params;
165 static double *fit_x = 0, *fit_y = 0, *fit_z = 0, *err_data = 0, *a = 0;
166 static TBOOLEAN ctrlc_flag = FALSE;
167 static TBOOLEAN user_stop = FALSE;
169 static struct udft_entry func;
171 typedef char fixstr[MAX_ID_LEN+1];
173 static fixstr *par_name;
175 static double startup_lambda = 0, lambda_down_factor = LAMBDA_DOWN_FACTOR, lambda_up_factor = LAMBDA_UP_FACTOR;
177 /*****************************************************************
179 *****************************************************************/
181 static RETSIGTYPE ctrlc_handle __PROTO((int an_int));
182 static void ctrlc_setup __PROTO((void));
183 static marq_res_t marquardt __PROTO((double a[], double **alpha, double *chisq,
185 static TBOOLEAN analyze __PROTO((double a[], double **alpha, double beta[],
187 static void calculate __PROTO((double *zfunc, double **dzda, double a[]));
188 static void call_gnuplot __PROTO((double *par, double *data));
189 static TBOOLEAN fit_interrupt __PROTO((void));
190 static TBOOLEAN regress __PROTO((double a[]));
191 static void show_fit __PROTO((int i, double chisq, double last_chisq, double *a,
192 double lambda, FILE * device));
193 static void log_axis_restriction __PROTO((FILE *log_f, AXIS_INDEX axis));
194 static TBOOLEAN is_empty __PROTO((char *s));
195 static TBOOLEAN is_variable __PROTO((char *s));
196 static double getdvar __PROTO((const char *varname));
197 static int getivar __PROTO((const char *varname));
198 static void setvar __PROTO((char *varname, struct value data));
199 static char *get_next_word __PROTO((char **s, char *subst));
200 static double createdvar __PROTO((char *varname, double value));
201 static void splitpath __PROTO((char *s, char *p, char *f));
202 static void backup_file __PROTO((char *, const char *));
204 #ifdef GP_FIT_ERRVARS
205 static void setvarerr __PROTO((char *varname, double value));
208 /*****************************************************************
209 Small function to write the last fit command into a file
210 Arg: Pointer to the file; if NULL, nothing is written,
211 but only the size of the string is returned.
212 *****************************************************************/
215 wri_to_fil_last_fit_cmd(FILE *fp)
218 return strlen(last_fit_command);
220 return (size_t) fputs(last_fit_command, fp);
224 /*****************************************************************
225 This is called when a SIGINT occurs during fit
226 *****************************************************************/
229 ctrlc_handle(int an_int)
231 (void) an_int; /* avoid -Wunused warning */
232 /* reinstall signal handler (necessary on SysV) */
233 (void) signal(SIGINT, (sigfunc) ctrlc_handle);
238 /*****************************************************************
239 setup the ctrl_c signal handler
240 *****************************************************************/
245 * MSDOS defines signal(SIGINT) but doesn't handle it through
246 * real interrupts. So there remain cases in which a ctrl-c may
247 * be uncaught by signal. We must use kbhit() instead that really
248 * serves the keyboard interrupt (or write an own interrupt func
249 * which also generates #ifdefs)
251 * I hope that other OSes do it better, if not... add #ifdefs :-(
253 #if (defined(__EMX__) || !defined(MSDOS) && !defined(DOS386)) && !defined(ATARI) && !defined(MTOS)
254 (void) signal(SIGINT, (sigfunc) ctrlc_handle);
259 /*****************************************************************
260 getch that handles also function keys etc.
261 *****************************************************************/
262 #if defined(MSDOS) || defined(DOS386)
264 /* HBB 980317: added a prototype... */
265 int getchx __PROTO((void));
271 if (!c || c == 0xE0) {
279 /*****************************************************************
280 in case of fatal errors
281 *****************************************************************/
287 memcpy(fitbuf, " ", 9); /* start after GNUPLOT> */
288 sp = strchr(fitbuf, NUL);
289 while (*--sp == '\n');
290 strcpy(sp + 1, "\n\n"); /* terminate with exactly 2 newlines */
291 fputs(fitbuf, STANDARD);
293 fprintf(log_f, "BREAK: %s", fitbuf);
294 (void) fclose(log_f);
298 free_at(func.at); /* release perm. action table */
299 func.at = (struct at_type *) NULL;
301 /* restore original SIGINT function */
302 #if !defined(ATARI) && !defined(MTOS)
306 bail_to_command_line();
309 /* HBB 990829: removed the debug print routines */
310 /*****************************************************************
311 Marquardt's nonlinear least squares fit
312 *****************************************************************/
314 marquardt(double a[], double **C, double *chisq, double *lambda)
317 static double *da = 0, /* delta-step of the parameter */
318 *temp_a = 0, /* temptative new params set */
319 *d = 0, *tmp_d = 0, **tmp_C = 0, *residues = 0;
322 /* Initialization when lambda == -1 */
324 if (*lambda == -1) { /* Get first chi-square check */
325 TBOOLEAN analyze_ret;
327 temp_a = vec(num_params);
328 d = vec(num_data + num_params);
329 tmp_d = vec(num_data + num_params);
330 da = vec(num_params);
331 residues = vec(num_data + num_params);
332 tmp_C = matr(num_data + num_params, num_params);
334 analyze_ret = analyze(a, C, d, chisq);
336 /* Calculate a useful startup value for lambda, as given by Schwarz */
337 /* FIXME: this is doesn't turn out to be much better, really... */
338 if (startup_lambda != 0)
339 *lambda = startup_lambda;
342 for (i = 0; i < num_data; i++)
343 for (j = 0; j < num_params; j++)
344 *lambda += C[i][j] * C[i][j];
345 *lambda = sqrt(*lambda / num_data / num_params);
348 /* Fill in the lower square part of C (the diagonal is filled in on
349 each iteration, see below) */
350 for (i = 0; i < num_params; i++)
351 for (j = 0; j < i; j++)
352 C[num_data + i][j] = 0, C[num_data + j][i] = 0;
353 return analyze_ret ? OK : ML_ERROR;
355 /* once converged, free dynamic allocated vars */
366 /* Givens calculates in-place, so make working copies of C and d */
368 for (j = 0; j < num_data + num_params; j++)
369 memcpy(tmp_C[j], C[j], num_params * sizeof(double));
370 memcpy(tmp_d, d, num_data * sizeof(double));
372 /* fill in additional parts of tmp_C, tmp_d */
374 for (i = 0; i < num_params; i++) {
375 /* fill in low diag. of tmp_C ... */
376 tmp_C[num_data + i][i] = *lambda;
377 /* ... and low part of tmp_d */
378 tmp_d[num_data + i] = 0;
381 /* FIXME: residues[] isn't used at all. Why? Should it be used? */
383 Givens(tmp_C, tmp_d, da, residues, num_params + num_data, num_params, 1);
385 /* check if trial did ameliorate sum of squares */
387 for (j = 0; j < num_params; j++)
388 temp_a[j] = a[j] + da[j];
390 if (!analyze(temp_a, tmp_C, tmp_d, &tmp_chisq)) {
391 /* FIXME: will never be reached: always returns TRUE */
394 if (tmp_chisq < *chisq) { /* Success, accept new solution */
395 if (*lambda > MIN_LAMBDA) {
396 (void) putc('/', stderr);
397 *lambda /= lambda_down_factor;
400 for (j = 0; j < num_data; j++) {
401 memcpy(C[j], tmp_C[j], num_params * sizeof(double));
404 for (j = 0; j < num_params; j++)
407 } else { /* failure, increase lambda and return */
408 (void) putc('*', stderr);
409 *lambda *= lambda_up_factor;
415 /*****************************************************************
416 compute chi-square and numeric derivations
417 *****************************************************************/
418 /* used by marquardt to evaluate the linearized fitting matrix C and
419 * vector d, fills in only the top part of C and d I don't use a
420 * temporary array zfunc[] any more. Just use d[] instead. */
421 /* FIXME: in the new code, this function doesn't really do enough to
422 * be useful. Maybe it ought to be deleted, i.e. integrated with
425 analyze(double a[], double **C, double d[], double *chisq)
432 for (i = 0; i < num_data; i++) {
433 /* note: order reversed, as used by Schwarz */
434 d[i] = (d[i] - fit_z[i]) / err_data[i];
435 *chisq += d[i] * d[i];
436 for (j = 0; j < num_params; j++)
437 C[i][j] /= err_data[i];
439 /* FIXME: why return a value that is always TRUE ? */
444 /* To use the more exact, but slower two-side formula, activate the
446 /*#define TWO_SIDE_DIFFERENTIATION */
447 /*****************************************************************
448 compute function values and partial derivatives of chi-square
449 *****************************************************************/
451 calculate(double *zfunc, double **dzda, double a[])
455 double *tmp_high, *tmp_pars;
456 #ifdef TWO_SIDE_DIFFERENTIATION
460 tmp_high = vec(num_data); /* numeric derivations */
461 #ifdef TWO_SIDE_DIFFERENTIATION
462 tmp_low = vec(num_data);
464 tmp_pars = vec(num_params);
466 /* first function values */
468 call_gnuplot(a, zfunc);
470 /* then derivatives */
472 for (p = 0; p < num_params; p++)
474 for (p = 0; p < num_params; p++) {
475 tmp_a = fabs(a[p]) < NEARLY_ZERO ? NEARLY_ZERO : a[p];
476 tmp_pars[p] = tmp_a * (1 + DELTA);
477 call_gnuplot(tmp_pars, tmp_high);
478 #ifdef TWO_SIDE_DIFFERENTIATION
479 tmp_pars[p] = tmp_a * (1 - DELTA);
480 call_gnuplot(tmp_pars, tmp_low);
482 for (k = 0; k < num_data; k++)
483 #ifdef TWO_SIDE_DIFFERENTIATION
484 dzda[k][p] = (tmp_high[k] - tmp_low[k]) / (2 * tmp_a * DELTA);
486 dzda[k][p] = (tmp_high[k] - zfunc[k]) / (tmp_a * DELTA);
491 #ifdef TWO_SIDE_DIFFERENTIATION
499 /*****************************************************************
500 call internal gnuplot functions
501 *****************************************************************/
503 call_gnuplot(double *par, double *data)
508 /* set parameters first */
509 for (i = 0; i < num_params; i++) {
510 (void) Gcomplex(&v, par[i], 0.0);
511 setvar(par_name[i], v);
514 for (i = 0; i < num_data; i++) {
515 /* calculate fit-function value */
516 (void) Gcomplex(&func.dummy_values[0], fit_x[i], 0.0);
517 (void) Gcomplex(&func.dummy_values[1], fit_y[i], 0.0);
518 evaluate_at(func.at, &v);
520 Eex("Undefined value during function evaluation");
526 /*****************************************************************
527 handle user interrupts during fit
528 *****************************************************************/
533 fputs("\n\n(S)top fit, (C)ontinue, (E)xecute FIT_SCRIPT: ", STANDARD);
534 switch (getc(stdin)) {
539 fputs("Stop.", STANDARD);
545 fputs("Continue.", STANDARD);
554 tmp = (*fit_script) ? fit_script : DEFAULT_CMD;
555 fprintf(STANDARD, "executing: %s", tmp);
556 /* set parameters visible to gnuplot */
557 for (i = 0; i < num_params; i++) {
558 (void) Gcomplex(&v, a[i], 0.0);
559 setvar(par_name[i], v);
561 safe_strncpy(gp_input_line, tmp, gp_input_line_len);
570 /*****************************************************************
571 frame routine for the marquardt-fit
572 *****************************************************************/
576 double **covar, *dpar, **C, chisq, last_chisq, lambda;
579 struct udvt_entry *v; /* For exporting results to the user */
581 chisq = last_chisq = INFINITY;
582 C = matr(num_data + num_params, num_params);
583 lambda = -1; /* use sign as flag */
584 iter = 0; /* iteration counter */
586 /* ctrlc now serves as Hotkey */
589 /* Initialize internal variables and 1st chi-square check */
590 if ((res = marquardt(a, C, &chisq, &lambda)) == ML_ERROR)
591 Eex("FIT: error occurred during fit");
594 show_fit(iter, chisq, chisq, a, lambda, STANDARD);
595 show_fit(iter, chisq, chisq, a, lambda, log_f);
597 /* Reset flag describing fit result status */
598 v = add_udv_by_name("FIT_CONVERGED");
599 v->udv_undef = FALSE;
600 Ginteger(&v->udv_value, 0);
602 /* MAIN FIT LOOP: do the regression iteration */
604 /* HBB 981118: initialize new variable 'user_break' */
609 * MSDOS defines signal(SIGINT) but doesn't handle it through
610 * real interrupts. So there remain cases in which a ctrl-c may
611 * be uncaught by signal. We must use kbhit() instead that really
612 * serves the keyboard interrupt (or write an own interrupt func
613 * which also generates #ifdefs)
615 * I hope that other OSes do it better, if not... add #ifdefs :-(
616 * EMX does not have kbhit.
618 * HBB: I think this can be enabled for DJGPP V2. SIGINT is actually
619 * handled there, AFAIK.
621 #if ((defined(MSDOS) || defined(DOS386)) && !defined(__EMX__)) || defined(ATARI) || defined(MTOS)
631 show_fit(iter, chisq, last_chisq, a, lambda, STANDARD);
633 if (!fit_interrupt()) /* handle keys */
640 if ((res = marquardt(a, C, &chisq, &lambda)) == BETTER)
641 show_fit(iter, chisq, last_chisq, a, lambda, STANDARD);
642 } while ((res != ML_ERROR)
643 && (lambda < MAX_LAMBDA)
644 && ((maxiter == 0) || (iter <= maxiter))
646 || ((chisq > NEARLY_ZERO)
647 ? ((last_chisq - chisq) / chisq)
648 : (last_chisq - chisq)) > epsilon
654 #if !defined(ATARI) && !defined(MTOS)
655 /* restore original SIGINT function */
659 /* HBB 970304: the maxiter patch: */
660 if ((maxiter > 0) && (iter > maxiter)) {
661 Dblf2("\nMaximum iteration count (%d) reached. Fit stopped.\n", maxiter);
662 } else if (user_stop) {
663 Dblf2("\nThe fit was stopped by the user after %d iterations.\n", iter);
665 Dblf2("\nAfter %d iterations the fit converged.\n", iter);
666 v = add_udv_by_name("FIT_CONVERGED");
667 v->udv_undef = FALSE;
668 Ginteger(&v->udv_value, 1);
671 Dblf2("final sum of squares of residuals : %g\n", chisq);
672 if (chisq > NEARLY_ZERO) {
673 Dblf2("rel. change during last iteration : %g\n\n", (chisq - last_chisq) / chisq);
675 Dblf2("abs. change during last iteration : %g\n\n", (chisq - last_chisq));
679 Eex("FIT: error occurred during fit");
681 /* compute errors in the parameters */
683 #ifdef GP_FIT_ERRVARS
684 if (fit_errorvariables)
685 /* Set error variable to zero before doing this */
686 /* Thus making sure they are created */
687 for (i = 0; i < num_params; i++)
688 setvarerr(par_name[i], 0.0);
691 if (num_data == num_params) {
694 Dblf("\nExactly as many data points as there are parameters.\n");
695 Dblf("In this degenerate case, all errors are zero by definition.\n\n");
696 Dblf("Final set of parameters \n");
697 Dblf("======================= \n\n");
698 for (k = 0; k < num_params; k++)
699 Dblf3("%-15.15s = %-15g\n", par_name[k], a[k]);
700 } else if (chisq < NEARLY_ZERO) {
703 Dblf("\nHmmmm.... Sum of squared residuals is zero. Can't compute errors.\n\n");
704 Dblf("Final set of parameters \n");
705 Dblf("======================= \n\n");
706 for (k = 0; k < num_params; k++)
707 Dblf3("%-15.15s = %-15g\n", par_name[k], a[k]);
709 int ndf = num_data - num_params;
710 double stdfit = sqrt(chisq/ndf);
712 Dblf2("degrees of freedom (FIT_NDF) : %d\n", ndf);
713 Dblf2("rms of residuals (FIT_STDFIT) = sqrt(WSSR/ndf) : %g\n", stdfit);
714 Dblf2("variance of residuals (reduced chisquare) = WSSR/ndf : %g\n\n", chisq/ndf);
716 /* Export these to user-accessible variables */
717 v = add_udv_by_name("FIT_NDF");
718 v->udv_undef = FALSE;
719 Ginteger(&v->udv_value, ndf);
720 v = add_udv_by_name("FIT_STDFIT");
721 v->udv_undef = FALSE;
722 Gcomplex(&v->udv_value, stdfit, 0);
723 v = add_udv_by_name("FIT_WSSR");
724 v->udv_undef = FALSE;
725 Gcomplex(&v->udv_value, chisq, 0);
727 /* get covariance-, Korrelations- and Kurvature-Matrix */
728 /* and errors in the parameters */
730 /* compute covar[][] directly from C */
731 Givens(C, 0, 0, 0, num_data, num_params, 0);
733 /* Use lower square of C for covar */
734 covar = C + num_data;
735 Invert_RtR(C, covar, num_params);
737 /* calculate unscaled parameter errors in dpar[]: */
738 dpar = vec(num_params);
739 for (i = 0; i < num_params; i++) {
740 /* FIXME: can this still happen ? */
741 if (covar[i][i] <= 0.0) /* HBB: prevent floating point exception later on */
742 Eex("Calculation error: non-positive diagonal element in covar. matrix");
743 dpar[i] = sqrt(covar[i][i]);
746 /* transform covariances into correlations */
747 for (i = 0; i < num_params; i++) {
748 /* only lower triangle needs to be handled */
749 for (j = 0; j <= i; j++)
750 covar[i][j] /= dpar[i] * dpar[j];
753 /* scale parameter errors based on chisq */
754 chisq = sqrt(chisq / (num_data - num_params));
755 for (i = 0; i < num_params; i++)
758 Dblf("Final set of parameters Asymptotic Standard Error\n");
759 Dblf("======================= ==========================\n\n");
761 for (i = 0; i < num_params; i++) {
762 double temp = (fabs(a[i]) < NEARLY_ZERO)
764 : fabs(100.0 * dpar[i] / a[i]);
766 Dblf6("%-15.15s = %-15g %-3.3s %-12.4g (%.4g%%)\n",
767 par_name[i], a[i], PLUSMINUS, dpar[i], temp);
768 #ifdef GP_FIT_ERRVARS
769 if (fit_errorvariables)
770 setvarerr(par_name[i], dpar[i]);
774 Dblf("\n\ncorrelation matrix of the fit parameters:\n\n");
777 for (j = 0; j < num_params; j++)
778 Dblf2("%-6.6s ", par_name[j]);
781 for (i = 0; i < num_params; i++) {
782 Dblf2("%-15.15s", par_name[i]);
783 for (j = 0; j <= i; j++) {
784 /* Only print lower triangle of symmetric matrix */
785 Dblf2("%6.3f ", covar[i][j]);
793 /* HBB 990220: re-imported this snippet from older versions. Finally,
794 * some user noticed that it *is* necessary, after all. Not even
795 * Carsten Grammes himself remembered what it was for... :-(
796 * The thing is: the value of the last parameter is not reset to
797 * its original one after the derivatives have been calculated
799 /* restore last parameter's value (not done by calculate) */
802 Gcomplex(&val, a[num_params - 1], 0.0);
803 setvar(par_name[num_params - 1], val);
806 /* call destructor for allocated vars */
807 lambda = -2; /* flag value, meaning 'destruct!' */
808 (void) marquardt(a, C, &chisq, &lambda);
815 /*****************************************************************
816 display actual state of the fit
817 *****************************************************************/
821 double chisq, double last_chisq,
822 double *a, double lambda,
827 fprintf(device, "\n\n\
829 WSSR : %-15g delta(WSSR)/WSSR : %g\n\
830 delta(WSSR) : %-15g limit for stopping : %g\n\
831 lambda : %g\n\n%s parameter values\n\n",
832 i, chisq, chisq > NEARLY_ZERO ? (chisq - last_chisq) / chisq : 0.0,
833 chisq - last_chisq, epsilon, lambda,
834 (i > 0 ? "resultant" : "initial set of free"));
835 for (k = 0; k < num_params; k++)
836 fprintf(device, "%-15.15s = %g\n", par_name[k], a[k]);
841 /*****************************************************************
842 is_empty: check for valid string entries
843 *****************************************************************/
847 while (*s == ' ' || *s == '\t' || *s == '\n')
849 return (TBOOLEAN) (*s == '#' || *s == '\0');
853 /*****************************************************************
854 get next word of a multi-word string, advance pointer
855 *****************************************************************/
857 get_next_word(char **s, char *subst)
861 while (*tmp == ' ' || *tmp == '\t' || *tmp == '=')
863 if (*tmp == '\n' || *tmp == '\0') /* not found */
865 if ((*s = strpbrk(tmp, " =\t\n")) == NULL)
866 *s = tmp + strlen(tmp);
873 /*****************************************************************
874 check for variable identifiers
875 *****************************************************************/
880 if (!isalnum((unsigned char) *s) && *s != '_')
887 /*****************************************************************
889 *****************************************************************/
893 func.at = (struct at_type *) NULL; /* need to parse 1 time */
897 /*****************************************************************
898 Set a GNUPLOT user-defined variable
899 ******************************************************************/
902 setvar(char *varname, struct value data)
904 struct udvt_entry *udv_ptr = add_udv_by_name(varname);
905 udv_ptr->udv_value = data;
906 udv_ptr->udv_undef = FALSE;
909 #ifdef GP_FIT_ERRVARS
910 /*****************************************************************
911 Set a GNUPLOT user-defined variable for an error
912 variable: so take the parameter name, turn it
913 into an error parameter name (e.g. a to a_err)
915 ******************************************************************/
917 setvarerr(char *varname, double value)
919 struct value errval; /* This will hold the gnuplot value created from value*/
920 char* pErrValName; /* The name of the (new) error variable */
922 /* Create the variable name by appending _err */
923 pErrValName = gp_alloc(strlen(varname)+sizeof(char)+4, "");
925 sprintf(pErrValName,"%s_err",varname);
926 Gcomplex(&errval, value, 0.0);
927 setvar(pErrValName, errval);
930 #endif /* GP_FIT_ERRVARS */
932 /*****************************************************************
933 Read INTGR Variable value, return 0 if undefined or wrong type
934 *****************************************************************/
936 getivar(const char *varname)
938 struct udvt_entry *udv_ptr = first_udv;
941 if (!strcmp(varname, udv_ptr->udv_name))
942 return udv_ptr->udv_value.type == INTGR
943 ? udv_ptr->udv_value.v.int_val /* valid */
944 : 0; /* wrong type */
945 udv_ptr = udv_ptr->next_udv;
947 return 0; /* not in table */
951 /*****************************************************************
952 Read DOUBLE Variable value, return 0 if undefined or wrong type
953 I dont think it's a problem that it's an integer - div
954 *****************************************************************/
956 getdvar(const char *varname)
958 struct udvt_entry *udv_ptr = first_udv;
960 for (; udv_ptr; udv_ptr = udv_ptr->next_udv)
961 if (strcmp(varname, udv_ptr->udv_name) == 0)
962 return real(&(udv_ptr->udv_value));
964 /* get here => not found */
968 /*****************************************************************
970 - convert it from integer to real if necessary
971 - create it with value INITIAL_VALUE if not found or undefined
972 *****************************************************************/
974 createdvar(char *varname, double value)
976 struct udvt_entry *udv_ptr = first_udv;
978 for (; udv_ptr; udv_ptr = udv_ptr->next_udv)
979 if (strcmp(varname, udv_ptr->udv_name) == 0) {
980 if (udv_ptr->udv_undef) {
981 udv_ptr->udv_undef = 0;
982 (void) Gcomplex(&udv_ptr->udv_value, value, 0.0);
983 } else if (udv_ptr->udv_value.type == INTGR) {
984 (void) Gcomplex(&udv_ptr->udv_value, (double) udv_ptr->udv_value.v.int_val, 0.0);
986 return real(&(udv_ptr->udv_value));
988 /* get here => not found */
991 struct value tempval;
992 (void) Gcomplex(&tempval, value, 0.0);
993 setvar(varname, tempval);
1000 /*****************************************************************
1001 Split Identifier into path and filename
1002 *****************************************************************/
1004 splitpath(char *s, char *p, char *f)
1006 char *tmp = s + strlen(s) - 1;
1008 while (tmp >= s && *tmp != '\\' && *tmp != '/' && *tmp != ':')
1010 /* FIXME HBB 20010121: unsafe! Sizes of 'f' and 'p' are not known.
1011 * May write past buffer end. */
1013 memcpy(p, s, (size_t) (tmp - s + 1));
1014 p[tmp - s + 1] = NUL;
1018 /* argument: char *fn */
1019 #define VALID_FILENAME(fn) ((fn) != NULL && (*fn) != '\0')
1021 /*****************************************************************
1022 write the actual parameters to start parameter file
1023 *****************************************************************/
1025 update(char *pfile, char *npfile)
1027 char fnam[256], path[256], sstr[256], pname[64], tail[127], *s = sstr, *tmp, c;
1028 char ifilename[256], *ofilename;
1033 /* update pfile npfile:
1034 if npfile is a valid file name,
1035 take pfile as input file and
1036 npfile as output file
1038 if (VALID_FILENAME(npfile)) {
1039 safe_strncpy(ifilename, pfile, sizeof(ifilename));
1042 #ifdef BACKUP_FILESYSTEM
1043 /* filesystem will keep original as previous version */
1044 safe_strncpy(ifilename, pfile, sizeof(ifilename));
1046 backup_file(ifilename, pfile); /* will Eex if it fails */
1051 /* split into path and filename */
1052 splitpath(ifilename, path, fnam);
1053 if (!(of = loadpath_fopen(ifilename, "r")))
1054 Eex2("parameter file %s could not be read", ifilename);
1056 if (!(nf = fopen(ofilename, "w")))
1057 Eex2("new parameter file %s could not be created", ofilename);
1059 while (fgets(s = sstr, sizeof(sstr), of) != NULL) {
1062 fputs(s, nf); /* preserve comments */
1065 if ((tmp = strchr(s, '#')) != NULL) {
1066 safe_strncpy(tail, tmp, sizeof(tail));
1071 tmp = get_next_word(&s, &c);
1072 if (!is_variable(tmp) || strlen(tmp) > MAX_ID_LEN) {
1075 Eex2("syntax error in parameter file %s", fnam);
1077 safe_strncpy(pname, tmp, sizeof(pname));
1078 /* next must be '=' */
1080 tmp = strchr(s, '=');
1084 Eex2("syntax error in parameter file %s", fnam);
1088 tmp = get_next_word(&s, &c);
1089 if (!sscanf(tmp, "%lf", &pval)) {
1092 Eex2("syntax error in parameter file %s", fnam);
1094 if ((tmp = get_next_word(&s, &c)) != NULL) {
1097 Eex2("syntax error in parameter file %s", fnam);
1101 if ((pval = getdvar(pname)) == 0)
1102 pval = (double) getivar(pname);
1104 sprintf(sstr, "%g", pval);
1105 if (!strchr(sstr, '.') && !strchr(sstr, 'e'))
1106 strcat(sstr, ".0"); /* assure CMPLX-type */
1108 fprintf(nf, "%-15.15s = %-15.15s %s", pname, sstr, tail);
1111 if (fclose(nf) || fclose(of))
1112 Eex("I/O error during update");
1117 /*****************************************************************
1118 Backup a file by renaming it to something useful. Return
1119 the new name in tofile
1120 *****************************************************************/
1122 /* tofile must point to a char array[] or allocated data. See update() */
1125 backup_file(char *tofile, const char *fromfile)
1127 #if defined (WIN32) || defined(MSDOS) || defined(VMS)
1131 /* win32 needs to have two attempts at the rename, since it may
1132 * be running on win32s with msdos 8.3 names
1135 /* first attempt, for all o/s other than MSDOS */
1139 strcpy(tofile, fromfile);
1142 /* replace all dots with _, since we will be adding the only
1143 * dot allowed in VMS names
1145 while ((tmpn = strchr(tofile, '.')) != NULL)
1149 strcat(tofile, BACKUP_SUFFIX);
1151 if (rename(fromfile, tofile) == 0)
1152 return; /* hurrah */
1156 #if defined (WIN32) || defined(MSDOS)
1158 /* first attempt for msdos. Second attempt for win32s */
1160 /* Copy only the first 8 characters of the filename, to comply
1161 * with the restrictions of FAT filesystems. */
1162 safe_strncpy(tofile, fromfile, 8 + 1);
1164 while ((tmpn = strchr(tofile, '.')) != NULL)
1167 strcat(tofile, BACKUP_SUFFIX);
1169 if (rename(fromfile, tofile) == 0)
1170 return; /* success */
1172 #endif /* win32 || msdos */
1174 /* get here => rename failed. */
1175 Eex3("Could not rename file %s to %s", fromfile, tofile);
1178 /* A modified copy of save.c:save_range(), but this one reports
1179 * _current_ values, not the 'set' ones, by default */
1181 log_axis_restriction(FILE *log_f, AXIS_INDEX axis)
1184 AXIS *this_axis = axis_array + axis;
1186 fprintf(log_f, "\t%s range restricted to [", axis_defaults[axis].name);
1187 if (this_axis->autoscale & AUTOSCALE_MIN) {
1189 } else if (this_axis->is_timedata) {
1191 gstrftime(s, 80, this_axis->timefmt, this_axis->min);
1195 fprintf(log_f, "%#g", this_axis->min);
1198 fputs(" : ", log_f);
1199 if (this_axis->autoscale & AUTOSCALE_MAX) {
1201 } else if (this_axis->is_timedata) {
1203 gstrftime(s, 80, this_axis->timefmt, this_axis->max);
1207 fprintf(log_f, "%#g", this_axis->max);
1209 fputs("]\n", log_f);
1212 /*****************************************************************
1213 Interface to the classic gnuplot-software
1214 *****************************************************************/
1219 /* HBB 20000430: revised this completely, to make it more similar to
1220 * what plot3drequest() does */
1222 int dummy_x = -1, dummy_y = -1; /* eg fit [u=...] [v=...] */
1223 int zrange_token = -1;
1224 TBOOLEAN is_a_3d_fit;
1230 int token1, token2, token3;
1231 char *tmp, *file_name;
1235 /* first look for a restricted x fit range... */
1237 /* put stuff into arrays to simplify access */
1238 AXIS_INIT3D(FIRST_X_AXIS, 0, 0);
1239 AXIS_INIT3D(FIRST_Y_AXIS, 0, 0);
1240 AXIS_INIT3D(FIRST_Z_AXIS, 0, 1);
1242 /* use global default indices in axis.c to simplify access to
1243 * per-axis variables */
1244 x_axis = FIRST_X_AXIS;
1245 y_axis = FIRST_Y_AXIS;
1246 z_axis = FIRST_Z_AXIS;
1248 PARSE_NAMED_RANGE(FIRST_X_AXIS, dummy_x);
1249 PARSE_NAMED_RANGE(FIRST_Y_AXIS, dummy_y);
1250 /* HBB 980401: new: allow restricting the z range as well */
1251 if (equals(c_token, "["))
1252 zrange_token = c_token;
1253 PARSE_RANGE(FIRST_Z_AXIS);
1256 /* now compile the function */
1262 func.at = NULL; /* in case perm_at() does int_error */
1267 /* set the dummy variable names */
1270 copy_str(c_dummy_var[0], dummy_x, MAX_ID_LEN);
1272 strcpy(c_dummy_var[0], set_dummy_var[0]);
1275 copy_str(c_dummy_var[0], dummy_y, MAX_ID_LEN);
1277 strcpy(c_dummy_var[1], set_dummy_var[1]);
1279 func.at = perm_at();
1285 file_name = try_to_get_string();
1287 int_error(c_token, "missing filename");
1289 /* use datafile module to parse the datafile and qualifiers */
1290 df_set_plot_mode(MODE_QUERY); /* Does nothing except for binary datafiles */
1291 columns = df_open(file_name, 4); /* up to 4 using specs allowed */
1294 int_error(NO_CARET,"Can't read data file");
1296 int_error(c_token, "Need 2 to 4 using specs");
1297 is_a_3d_fit = (columns == 4);
1299 /* The following patch was made by Remko Scharroo, 25-Mar-1999
1300 * We need to check if one of the columns is time data, like
1301 * in plot2d and plot3d */
1303 if (X_AXIS.is_timedata) {
1305 int_error(c_token, "Need full using spec for x time data");
1307 if (Y_AXIS.is_timedata) {
1309 int_error(c_token, "Need using spec for y time data");
1311 df_axis[0] = FIRST_X_AXIS;
1312 df_axis[1] = (is_a_3d_fit) ? FIRST_Y_AXIS : FIRST_Z_AXIS;
1313 /* don't parse delta_z as times */
1314 df_axis[2] = (is_a_3d_fit) ? FIRST_Z_AXIS : -1;
1317 /* HBB 20000430: No need for such a check for the z axis. That
1318 * axis is only used iff columns==4, i.e. the detection method for
1319 * 3d fits requires a 4 column 'using', already. */
1320 /* End of patch by Remko Scharroo */
1323 /* HBB 980401: if this is a single-variable fit, we shouldn't have
1324 * allowed a variable name specifier for 'y': */
1325 if ((dummy_y >= 0) && !is_a_3d_fit)
1326 int_error(dummy_y, "Can't re-name 'y' in a one-variable fit");
1328 /* HBB 981210: two range specs mean different things, depending
1329 * on wether this is a 2D or 3D fit */
1331 if (zrange_token != -1)
1332 int_error(zrange_token, "Three range-specs not allowed in on-variable fit");
1334 /* 2D fit, 2 ranges: second range is for *z*, not y: */
1335 Z_AXIS.autoscale = Y_AXIS.autoscale;
1336 if (!(Y_AXIS.autoscale & AUTOSCALE_MIN))
1337 Z_AXIS.min = Y_AXIS.min;
1338 if (!(Y_AXIS.autoscale & AUTOSCALE_MAX))
1339 Z_AXIS.max = Y_AXIS.max;
1342 /* defer actually reading the data until we have parsed the rest
1347 tmpd = getdvar(FITLIMIT); /* get epsilon if given explicitly */
1348 if (tmpd < 1.0 && tmpd > 0.0)
1351 /* HBB 970304: maxiter patch */
1352 maxiter = getivar(FITMAXITER);
1354 /* get startup value for lambda, if given */
1355 tmpd = getdvar(FITSTARTLAMBDA);
1358 startup_lambda = tmpd;
1359 printf("Lambda Start value set: %g\n", startup_lambda);
1361 /* get lambda up/down factor, if given */
1362 tmpd = getdvar(FITLAMBDAFACTOR);
1365 lambda_up_factor = lambda_down_factor = tmpd;
1366 printf("Lambda scaling factors reset: %g\n", lambda_up_factor);
1369 if ((tmp = getenv(FITSCRIPT)) != NULL)
1370 safe_strncpy(fit_script, tmp, sizeof(fit_script));
1373 char *logfile = getfitlogfile();
1375 if (!log_f && !(log_f = fopen(logfile, "a")))
1376 Eex2("could not open log-file %s", logfile);
1382 fputs("\n\n*******************************************************************************\n", log_f);
1383 (void) time(&timer);
1384 fprintf(log_f, "%s\n\n", ctime(&timer));
1388 m_capture(&line, token2, token3 - 1);
1389 fprintf(log_f, "FIT: data read from %s\n", line);
1393 /* HBB 20040201: somebody insisted they need this ... */
1394 if (AUTOSCALE_BOTH != (X_AXIS.autoscale & AUTOSCALE_BOTH))
1395 log_axis_restriction(log_f, x_axis);
1396 if (AUTOSCALE_BOTH != (Y_AXIS.autoscale & AUTOSCALE_BOTH))
1397 log_axis_restriction(log_f, y_axis);
1399 max_data = MAX_DATA;
1400 fit_x = vec(max_data); /* start with max. value */
1401 fit_y = vec(max_data);
1402 fit_z = vec(max_data);
1404 /* first read in experimental data */
1406 err_data = vec(max_data);
1409 while ((i = df_readline(v, 4)) != DF_EOF) {
1410 if (num_data >= max_data) {
1411 /* increase max_data by factor of 1.5 */
1412 max_data = (max_data * 3) / 2;
1414 || !redim_vec(&fit_x, max_data)
1415 || !redim_vec(&fit_y, max_data)
1416 || !redim_vec(&fit_z, max_data)
1417 || !redim_vec(&err_data, max_data)
1419 /* Some of the reallocations went bad: */
1421 Eex2("Out of memory in fit: too many datapoints (%d)?", max_data);
1423 /* Just so we know that the routine is at work: */
1424 fprintf(STANDARD, "Max. number of data points scaled up to: %d\n",
1427 } /* if (need to extend storage space) */
1432 case DF_FIRST_BLANK:
1433 case DF_SECOND_BLANK:
1436 Eex2("bad data on line %d of datafile", df_line_number);
1438 case 1: /* only z provided */
1440 v[0] = (double) df_datum;
1442 case 2: /* x and z */
1446 /* only if the explicitly asked for 4 columns do we
1447 * do a 3d fit. (We can get here if they didn't
1448 * specify a using spec, and the file has 4 columns
1450 case 4: /* x, y, z, error */
1453 /* else fall through */
1454 case 3: /* x, z, error */
1455 v[3] = v[2]; /* error */
1456 v[2] = v[1]; /* z */
1461 /* skip this point if it is out of range */
1462 if (!(X_AXIS.autoscale & AUTOSCALE_MIN)
1463 && (v[0] < X_AXIS.min))
1465 if (!(X_AXIS.autoscale & AUTOSCALE_MAX)
1466 && (v[0] > X_AXIS.max))
1469 if (!(Y_AXIS.autoscale & AUTOSCALE_MIN)
1470 && (v[1] < Y_AXIS.min))
1472 if (!(Y_AXIS.autoscale & AUTOSCALE_MAX)
1473 && (v[1] > Y_AXIS.max))
1476 /* HBB 980401: check *z* range for all fits */
1477 if (!(Z_AXIS.autoscale & AUTOSCALE_MIN)
1478 && (v[2] < Z_AXIS.min))
1480 if (!(Z_AXIS.autoscale & AUTOSCALE_MAX)
1481 && (v[2] > Z_AXIS.max))
1484 fit_x[num_data] = v[0];
1485 fit_y[num_data] = v[1];
1486 fit_z[num_data] = v[2];
1488 /* we only use error if _explicitly_ asked for by a
1491 err_data[num_data++] = (columns > 2) ? v[3] : 1;
1496 Eex("No data to fit ");
1498 /* now resize fields to actual length: */
1499 redim_vec(&fit_x, num_data);
1500 redim_vec(&fit_y, num_data);
1501 redim_vec(&fit_z, num_data);
1502 redim_vec(&err_data, num_data);
1504 fprintf(log_f, " #datapoints = %d\n", num_data);
1507 fputs(" residuals are weighted equally (unit weight)\n\n", log_f);
1512 m_capture(&line, token1, token2 - 1);
1513 fprintf(log_f, "function used for fitting: %s\n", line);
1517 /* read in parameters */
1519 max_params = MAX_PARAMS; /* HBB 971023: make this resizeable */
1521 if (!equals(c_token++, "via"))
1522 int_error(c_token, "Need via and either parameter list or file");
1524 a = vec(max_params);
1525 par_name = (fixstr *) gp_alloc((max_params + 1) * sizeof(fixstr),
1529 if (isstringvalue(c_token)) { /* It's a parameter *file* */
1533 char sstr[MAX_LINE_LEN + 1];
1536 static char *viafile = NULL;
1537 free(viafile); /* Free previous name, if any */
1538 viafile = try_to_get_string(); /* Cannot fail since isstringvalue succeeded */
1539 fprintf(log_f, "fitted parameters and initial values from file: %s\n\n", viafile);
1540 if (!(f = loadpath_fopen(viafile, "r")))
1541 Eex2("could not read parameter-file \"%s\"", viafile);
1543 /* get parameters and values out of file and ignore fixed ones */
1546 if (!fgets(s = sstr, sizeof(sstr), f)) /* EOF found */
1548 if ((tmp = strstr(s, GP_FIXED)) != NULL) { /* ignore fixed params */
1550 fprintf(log_f, "FIXED: %s\n", s);
1554 if ((tmp = strchr(s, '#')) != NULL)
1558 tmp = get_next_word(&s, &c);
1559 if (!is_variable(tmp) || strlen(tmp) > MAX_ID_LEN) {
1561 Eex("syntax error in parameter file");
1563 safe_strncpy(par_name[num_params], tmp, sizeof(fixstr));
1564 /* next must be '=' */
1566 tmp = strchr(s, '=');
1569 Eex("syntax error in parameter file");
1573 tmp = get_next_word(&s, &c);
1574 if (sscanf(tmp, "%lf", &tmp_par) != 1) {
1576 Eex("syntax error in parameter file");
1578 /* make fixed params visible to GNUPLOT */
1580 struct value tempval;
1581 (void) Gcomplex(&tempval, tmp_par, 0.0);
1582 /* use parname as temp */
1583 setvar(par_name[num_params], tempval);
1585 if (num_params >= max_params) {
1586 max_params = (max_params * 3) / 2;
1588 || !redim_vec(&a, max_params)
1589 || !(par_name = gp_realloc(par_name, (max_params + 1) * sizeof(fixstr), "fit param resize"))
1592 Eex("Out of memory in fit: too many parameters?");
1595 a[num_params++] = tmp_par;
1598 if ((tmp = get_next_word(&s, &c)) != NULL) {
1600 Eex("syntax error in parameter file");
1606 /* not a string after via: it's a variable listing */
1608 fputs("fitted parameters initialized with current variable values\n\n", log_f);
1610 if (!isletter(c_token))
1611 Eex("no parameter specified");
1612 capture(par_name[num_params], c_token, c_token, (int) sizeof(par_name[0]));
1613 if (num_params >= max_params) {
1614 max_params = (max_params * 3) / 2;
1616 || !redim_vec(&a, max_params)
1617 || !(par_name = gp_realloc(par_name, (max_params + 1) * sizeof(fixstr), "fit param resize"))
1619 Eex("Out of memory in fit: too many parameters?");
1622 /* create variable if it doesn't exist */
1623 a[num_params] = createdvar(par_name[num_params], INITIAL_VALUE);
1625 } while (equals(++c_token, ",") && ++c_token);
1628 redim_vec(&a, num_params);
1629 par_name = (fixstr *) gp_realloc(par_name, (num_params + 1) * sizeof(fixstr), "fit param");
1631 if (num_data < num_params)
1632 Eex("Number of data points smaller than number of parameters");
1634 /* avoid parameters being equal to zero */
1635 for (i = 0; i < num_params; i++)
1640 if (num_params == 0)
1641 int_warn(NO_CARET, "No fittable parameters!\n");
1645 (void) fclose(log_f);
1654 free_at(func.at); /* release perm. action table */
1655 func.at = (struct at_type *) NULL;
1657 safe_strncpy(last_fit_command, gp_input_line, sizeof(last_fit_command));
1660 #if defined(ATARI) || defined(MTOS)
1665 struct timeval timeout;
1668 timeout.tv_usec = 0;
1669 rfds = (fd_set) (1L << fileno(stdin));
1670 return ((select(0, SELECT_TYPE_ARG234 &rfds, NULL, NULL, SELECT_TYPE_ARG5 &timeout) > 0) ? 1 : 0);
1675 * Print msg to stderr and log file
1677 #if defined(VA_START) && defined(STDC_HEADERS)
1679 Dblfn(const char *fmt, ...)
1682 Dblfn(const char *fmt, va_dcl)
1688 VA_START(args, fmt);
1689 # if defined(HAVE_VFPRINTF) || _LIBC
1690 vfprintf(STANDARD, fmt, args);
1692 VA_START(args, fmt);
1693 vfprintf(log_f, fmt, args);
1695 _doprnt(fmt, args, STANDARD);
1696 _doprnt(fmt, args, log_f);
1700 fprintf(STANDARD, fmt, a1, a2, a3, a4, a5, a6, a7, a8);
1701 fprintf(log_f, fmt, a1, a2, a3, a4, a5, a6, a7, a8);
1702 #endif /* VA_START */
1705 /* HBB/H.Harders NEW 20020927: make fit log filename exchangeable at
1706 * run-time, not only by setting an environment variable. */
1710 char *logfile = NULL;
1712 if (fitlogfile == NULL) {
1713 char *tmp = getenv(GNUFITLOG); /* open logfile */
1715 if (tmp != NULL && *tmp != '\0') {
1716 char *tmp2 = tmp + (strlen(tmp) - 1);
1718 /* if given log file name ends in path separator, treat it
1719 * as a directory to store the default "fit.log" in */
1720 if (*tmp2 == '/' || *tmp2 == '\\') {
1721 logfile = gp_alloc(strlen(tmp)
1722 + strlen(fitlogfile_default) + 1,
1724 strcpy(logfile, tmp);
1725 strcat(logfile, fitlogfile_default);
1727 logfile = gp_strdup(tmp);
1730 logfile = gp_strdup(fitlogfile_default);
1733 logfile = gp_strdup(fitlogfile);