Annotation of OpenXM_contrib/gnuplot/fit.c, Revision 1.1.1.2
1.1 maekawa 1: #ifndef lint
1.1.1.2 ! maekawa 2: static char *RCSid = "$Id: fit.c,v 1.23.2.3 1999/10/01 18:22:14 lhecking Exp $";
1.1 maekawa 3: #endif
4:
5: /*
6: * Nonlinear least squares fit according to the
7: * Marquardt-Levenberg-algorithm
8: *
9: * added as Patch to Gnuplot (v3.2 and higher)
10: * by Carsten Grammes
11: * Experimental Physics, University of Saarbruecken, Germany
12: *
13: * Internet address: cagr@rz.uni-sb.de
14: *
15: * Copyright of this module: 1993, 1998 Carsten Grammes
16: *
17: * Permission to use, copy, and distribute this software and its
18: * documentation for any purpose with or without fee is hereby granted,
19: * provided that the above copyright notice appear in all copies and
20: * that both that copyright notice and this permission notice appear
21: * in supporting documentation.
22: *
23: * This software is provided "as is" without express or implied warranty.
24: *
25: * 930726: Recoding of the Unix-like raw console I/O routines by:
26: * Michele Marziani (marziani@ferrara.infn.it)
27: * drd: start unitialised variables at 1 rather than NEARLY_ZERO
28: * (fit is more likely to converge if started from 1 than 1e-30 ?)
29: *
30: * HBB (Broeker@physik.rwth-aachen.de) : fit didn't calculate the errors
31: * in the 'physically correct' (:-) way, if a third data column containing
32: * the errors (or 'uncertainties') of the input data was given. I think
33: * I've fixed that, but I'm not sure I really understood the M-L-algo well
34: * enough to get it right. I deduced my change from the final steps of the
35: * equivalent algorithm for the linear case, which is much easier to
36: * understand. (I also made some minor, mostly cosmetic changes)
37: *
38: * HBB (again): added error checking for negative covar[i][i] values and
39: * for too many parameters being specified.
40: *
41: * drd: allow 3d fitting. Data value is now called fit_z internally,
42: * ie a 2d fit is z vs x, and a 3d fit is z vs x and y.
43: *
44: * Lars Hecking : review update command, for VMS in particular, where
45: * it is not necessary to rename the old file.
46: *
47: * HBB, 971023: lifted fixed limit on number of datapoints, and number
48: * of parameters.
49: */
50:
51:
52: #define FIT_MAIN
53:
54: #include <signal.h>
55:
56: #include "plot.h"
57: #include "matrix.h"
58: #include "fit.h"
59: #include "setshow.h" /* for load_range */
60: #include "alloc.h"
61:
62: #if defined(MSDOS) || defined(DOS386) /* non-blocking IO stuff */
63: # include <io.h>
64: # ifndef _Windows /* WIN16 does define MSDOS .... */
65: # include <conio.h>
66: # endif
67: # include <dos.h>
68: #else /* !(MSDOS || DOS386) */
69: # ifndef VMS
70: # include <fcntl.h>
71: # endif /* !VMS */
72: #endif /* !(MSDOS || DOS386) */
73:
74: #if defined(ATARI) || defined(MTOS)
75: # define getchx() Crawcin()
76: static int kbhit(void);
77: #endif
78:
79: #define STANDARD stderr /* compatible with gnuplot philosophy */
80:
81: #define BACKUP_SUFFIX ".old"
82:
83:
84: /* access external global variables (ought to make a globals.h someday) */
85:
86: extern struct udft_entry *dummy_func;
87: extern char dummy_var[MAX_NUM_VAR][MAX_ID_LEN+1];
88: extern char c_dummy_var[MAX_NUM_VAR][MAX_ID_LEN+1];
89: extern int c_token;
90: extern int df_datum, df_line_number;
91:
1.1.1.2 ! maekawa 92: /* following 2 external arrays are needed to use time data */
! 93: extern int datatype[];
! 94: extern int df_timecol[];
1.1 maekawa 95:
96: enum marq_res {
97: OK, ERROR, BETTER, WORSE
98: };
99: typedef enum marq_res marq_res_t;
100:
101: #ifdef INFINITY
102: # undef INFINITY
103: #endif
104:
105: #define INFINITY 1e30
106: #define NEARLY_ZERO 1e-30
107:
108: /* create new variables with this value (was NEARLY_ZERO) */
109: #define INITIAL_VALUE 1.0
110:
111: /* Relative change for derivatives */
112: #define DELTA 0.001
113:
114: #define MAX_DATA 2048
115: #define MAX_PARAMS 32
116: #define MAX_LAMBDA 1e20
117: #define MIN_LAMBDA 1e-20
118: #define LAMBDA_UP_FACTOR 10
119: #define LAMBDA_DOWN_FACTOR 10
120:
121: #if defined(MSDOS) || defined(OS2) || defined(DOS386)
122: # define PLUSMINUS "\xF1" /* plusminus sign */
123: #else
124: # define PLUSMINUS "+/-"
125: #endif
126:
127: #define LASTFITCMDLENGTH 511
128:
129: /* HBB 971023: new, allow for dynamic adjustment of these: */
130: static int max_data;
131: static int max_params;
132:
133: static double epsilon = 1e-5; /* convergence limit */
134: static int maxiter = 0; /* HBB 970304: maxiter patch */
135:
136: static char fit_script[128];
137: static char logfile[128] = "fit.log";
138: static char *FIXED = "# FIXED";
139: static char *GNUFITLOG = "FIT_LOG";
140: static char *FITLIMIT = "FIT_LIMIT";
141: static char *FITSTARTLAMBDA = "FIT_START_LAMBDA";
142: static char *FITLAMBDAFACTOR = "FIT_LAMBDA_FACTOR";
143: static char *FITMAXITER = "FIT_MAXITER"; /* HBB 970304: maxiter patch */
144: static char *FITSCRIPT = "FIT_SCRIPT";
145: static char *DEFAULT_CMD = "replot"; /* if no fitscript spec. */
146: static char last_fit_command[LASTFITCMDLENGTH+1] = "";
147:
148:
149: static FILE *log_f = NULL;
150:
151: static int num_data, num_params;
152: static int columns;
153: static double *fit_x = 0, *fit_y = 0, *fit_z = 0, *err_data = 0, *a = 0;
154: static TBOOLEAN ctrlc_flag = FALSE;
155: static TBOOLEAN user_stop = FALSE;
156:
157: static struct udft_entry func;
158:
159: typedef char fixstr[MAX_ID_LEN+1];
160:
161: static fixstr *par_name;
162:
163: static double startup_lambda = 0, lambda_down_factor = LAMBDA_DOWN_FACTOR,
164: lambda_up_factor = LAMBDA_UP_FACTOR;
165:
166: /*****************************************************************
167: internal Prototypes
168: *****************************************************************/
169:
170: static RETSIGTYPE ctrlc_handle __PROTO((int an_int));
171: static void ctrlc_setup __PROTO((void));
172: static void printmatrix __PROTO((double **C, int m, int n));
173: static void print_matrix_and_vectors __PROTO((double **C, double *d, double *r, int m, int n));
174: static marq_res_t marquardt __PROTO((double a[], double **alpha, double *chisq,
175: double *lambda));
176: static TBOOLEAN analyze __PROTO((double a[], double **alpha, double beta[],
177: double *chisq));
178: static void calculate __PROTO((double *zfunc, double **dzda, double a[]));
179: static void call_gnuplot __PROTO((double *par, double *data));
180: static TBOOLEAN fit_interrupt __PROTO((void));
181: static TBOOLEAN regress __PROTO((double a[]));
182: static void show_fit __PROTO((int i, double chisq, double last_chisq, double *a,
183: double lambda, FILE * device));
184: static TBOOLEAN is_empty __PROTO((char *s));
185: static TBOOLEAN is_variable __PROTO((char *s));
186: static double getdvar __PROTO((char *varname));
187: static double createdvar __PROTO((char *varname, double value));
188: static void splitpath __PROTO((char *s, char *p, char *f));
189: static void backup_file __PROTO((char *, const char *));
190:
191:
192: /*****************************************************************
193: Small function to write the last fit command into a file
194: Arg: Pointer to the file; if NULL, nothing is written,
195: but only the size of the string is returned.
196: *****************************************************************/
197:
198: size_t wri_to_fil_last_fit_cmd(fp)
199: FILE *fp;
200: {
201: if (fp == NULL)
202: return strlen(last_fit_command);
203: else
204: return (size_t) fputs(last_fit_command, fp);
205: }
206:
207:
208: /*****************************************************************
209: This is called when a SIGINT occurs during fit
210: *****************************************************************/
211:
212: static RETSIGTYPE ctrlc_handle(an_int)
213: int an_int;
214: {
215: #ifdef OS2
216: (void) signal(an_int, SIG_ACK);
217: #else
218: /* reinstall signal handler (necessary on SysV) */
219: (void) signal(SIGINT, (sigfunc) ctrlc_handle);
220: #endif
221: ctrlc_flag = TRUE;
222: }
223:
224:
225: /*****************************************************************
226: setup the ctrl_c signal handler
227: *****************************************************************/
228: static void ctrlc_setup()
229: {
230: /*
231: * MSDOS defines signal(SIGINT) but doesn't handle it through
232: * real interrupts. So there remain cases in which a ctrl-c may
233: * be uncaught by signal. We must use kbhit() instead that really
234: * serves the keyboard interrupt (or write an own interrupt func
235: * which also generates #ifdefs)
236: *
237: * I hope that other OSes do it better, if not... add #ifdefs :-(
238: */
239: #if (defined(__EMX__) || !defined(MSDOS) && !defined(DOS386)) && !defined(ATARI) && !defined(MTOS)
240: (void) signal(SIGINT, (sigfunc) ctrlc_handle);
241: #endif
242: }
243:
244:
245: /*****************************************************************
246: getch that handles also function keys etc.
247: *****************************************************************/
248: #if defined(MSDOS) || defined(DOS386)
249:
250: /* HBB 980317: added a prototype... */
251: int getchx __PROTO((void));
252:
253: int getchx()
254: {
255: register int c = getch();
256: if (!c || c == 0xE0) {
257: c <<= 8;
258: c |= getch();
259: }
260: return c;
261: }
262: #endif
263:
264: /*****************************************************************
265: in case of fatal errors
266: *****************************************************************/
267: void error_ex()
268: {
269: char *sp;
270:
271: strncpy(fitbuf, " ", 9); /* start after GNUPLOT> */
272: sp = strchr(fitbuf, NUL);
273: while (*--sp == '\n')
274: ;
275: strcpy(sp + 1, "\n\n"); /* terminate with exactly 2 newlines */
276: fputs(fitbuf, STANDARD);
277: if (log_f) {
278: fprintf(log_f, "BREAK: %s", fitbuf);
279: (void) fclose(log_f);
280: log_f = NULL;
281: }
282: if (func.at) {
283: free(func.at); /* release perm. action table */
284: func.at = (struct at_type *) NULL;
285: }
286: /* restore original SIGINT function */
287: #if !defined(ATARI) && !defined(MTOS)
288: interrupt_setup();
289: #endif
290:
291: bail_to_command_line();
292: }
293:
294:
295: /*****************************************************************
296: New utility routine: print a matrix (for debugging the alg.)
297: *****************************************************************/
298: static void printmatrix(C, m, n)
299: double **C;
300: int m, n;
301: {
302: int i, j;
303:
304: for (i = 0; i < m; i++) {
305: for (j = 0; j < n - 1; j++)
306: Dblf2("%.8g |", C[i][j]);
307: Dblf2("%.8g\n", C[i][j]);
308: }
309: Dblf("\n");
310: }
311:
312: /**************************************************************************
313: Yet another debugging aid: print matrix, with diff. and residue vector
314: **************************************************************************/
315: static void print_matrix_and_vectors(C, d, r, m, n)
316: double **C;
317: double *d, *r;
318: int m, n;
319: {
320: int i, j;
321:
322: for (i = 0; i < m; i++) {
323: for (j = 0; j < n; j++)
324: Dblf2("%8g ", C[i][j]);
325: Dblf3("| %8g | %8g\n", d[i], r[i]);
326: }
327: Dblf("\n");
328: }
329:
330:
331: /*****************************************************************
332: Marquardt's nonlinear least squares fit
333: *****************************************************************/
334: static marq_res_t marquardt(a, C, chisq, lambda)
335: double a[];
336: double **C;
337: double *chisq;
338: double *lambda;
339: {
340: int i, j;
341: static double *da = 0, /* delta-step of the parameter */
342: *temp_a = 0, /* temptative new params set */
343: *d = 0, *tmp_d = 0, **tmp_C = 0, *residues = 0;
344: double tmp_chisq;
345:
346: /* Initialization when lambda == -1 */
347:
348: if (*lambda == -1) { /* Get first chi-square check */
349: TBOOLEAN analyze_ret;
350:
351: temp_a = vec(num_params);
352: d = vec(num_data + num_params);
353: tmp_d = vec(num_data + num_params);
354: da = vec(num_params);
355: residues = vec(num_data + num_params);
356: tmp_C = matr(num_data + num_params, num_params);
357:
358: analyze_ret = analyze(a, C, d, chisq);
359:
360: /* Calculate a useful startup value for lambda, as given by Schwarz */
361: /* FIXME: this is doesn't turn out to be much better, really... */
362: if (startup_lambda != 0)
363: *lambda = startup_lambda;
364: else {
365: *lambda = 0;
366: for (i = 0; i < num_data; i++)
367: for (j = 0; j < num_params; j++)
368: *lambda += C[i][j] * C[i][j];
369: *lambda = sqrt(*lambda / num_data / num_params);
370: }
371:
372: /* Fill in the lower square part of C (the diagonal is filled in on
373: each iteration, see below) */
374: for (i = 0; i < num_params; i++)
375: for (j = 0; j < i; j++)
376: C[num_data + i][j] = 0, C[num_data + j][i] = 0;
377: /*printmatrix(C, num_data+num_params, num_params); */
378: return analyze_ret ? OK : ERROR;
379: }
380: /* once converged, free dynamic allocated vars */
381:
382: if (*lambda == -2) {
383: free(d);
384: free(tmp_d);
385: free(da);
386: free(temp_a);
387: free(residues);
388: free_matr(tmp_C);
389: return OK;
390: }
391: /* Givens calculates in-place, so make working copies of C and d */
392:
393: for (j = 0; j < num_data + num_params; j++)
394: memcpy(tmp_C[j], C[j], num_params * sizeof(double));
395: memcpy(tmp_d, d, num_data * sizeof(double));
396:
397: /* fill in additional parts of tmp_C, tmp_d */
398:
399: for (i = 0; i < num_params; i++) {
400: /* fill in low diag. of tmp_C ... */
401: tmp_C[num_data + i][i] = *lambda;
402: /* ... and low part of tmp_d */
403: tmp_d[num_data + i] = 0;
404: }
405: /* printmatrix(tmp_C, num_data+num_params, num_params); */
406:
407: /* FIXME: residues[] isn't used at all. Why? Should it be used? */
408:
409: Givens(tmp_C, tmp_d, da, residues, num_params + num_data, num_params, 1);
410: /*print_matrix_and_vectors (tmp_C, tmp_d, residues,
411: num_params+num_data, num_params); */
412:
413: /* check if trial did ameliorate sum of squares */
414:
415: for (j = 0; j < num_params; j++)
416: temp_a[j] = a[j] + da[j];
417:
418: if (!analyze(temp_a, tmp_C, tmp_d, &tmp_chisq)) {
419: /* FIXME: will never be reached: always returns TRUE */
420: return ERROR;
421: }
422:
423: if (tmp_chisq < *chisq) { /* Success, accept new solution */
424: if (*lambda > MIN_LAMBDA) {
425: (void) putc('/', stderr);
426: *lambda /= lambda_down_factor;
427: }
428: *chisq = tmp_chisq;
429: for (j = 0; j < num_data; j++) {
430: memcpy(C[j], tmp_C[j], num_params * sizeof(double));
431: d[j] = tmp_d[j];
432: }
433: for (j = 0; j < num_params; j++)
434: a[j] = temp_a[j];
435: return BETTER;
436: } else { /* failure, increase lambda and return */
437: (void) putc('*', stderr);
438: *lambda *= lambda_up_factor;
439: return WORSE;
440: }
441: }
442:
443:
444: /* FIXME: in the new code, this function doesn't really do enough to be
445: * useful. Maybe it ought to be deleted, i.e. integrated with
446: * calculate() ?
447: */
448: /*****************************************************************
449: compute chi-square and numeric derivations
450: *****************************************************************/
451: static TBOOLEAN analyze(a, C, d, chisq)
452: double a[];
453: double **C;
454: double d[];
455: double *chisq;
456: {
457: /*
458: * used by marquardt to evaluate the linearized fitting matrix C
459: * and vector d, fills in only the top part of C and d
460: * I don't use a temporary array zfunc[] any more. Just use
461: * d[] instead.
462: */
463: int i, j;
464:
465: *chisq = 0;
466: calculate(d, C, a);
467:
468: for (i = 0; i < num_data; i++) {
469: /* note: order reversed, as used by Schwarz */
470: d[i] = (d[i] - fit_z[i]) / err_data[i];
471: *chisq += d[i] * d[i];
472: for (j = 0; j < num_params; j++)
473: C[i][j] /= err_data[i];
474: }
475: /* FIXME: why return a value that is always TRUE ? */
476: return TRUE;
477: }
478:
479:
480: /* To use the more exact, but slower two-side formula, activate the
481: following line: */
482: /*#define TWO_SIDE_DIFFERENTIATION */
483: /*****************************************************************
484: compute function values and partial derivatives of chi-square
485: *****************************************************************/
486: static void calculate(zfunc, dzda, a)
487: double *zfunc;
488: double **dzda;
489: double a[];
490: {
491: int k, p;
492: double tmp_a;
493: double *tmp_high, *tmp_pars;
494: #ifdef TWO_SIDE_DIFFERENTIATION
495: double *tmp_low;
496: #endif
497:
498: tmp_high = vec(num_data); /* numeric derivations */
499: #ifdef TWO_SIDE_DIFFERENTIATION
500: tmp_low = vec(num_data);
501: #endif
502: tmp_pars = vec(num_params);
503:
504: /* first function values */
505:
506: call_gnuplot(a, zfunc);
507:
508: /* then derivatives */
509:
510: for (p = 0; p < num_params; p++)
511: tmp_pars[p] = a[p];
512: for (p = 0; p < num_params; p++) {
513: tmp_a = fabs(a[p]) < NEARLY_ZERO ? NEARLY_ZERO : a[p];
514: tmp_pars[p] = tmp_a * (1 + DELTA);
515: call_gnuplot(tmp_pars, tmp_high);
516: #ifdef TWO_SIDE_DIFFERENTIATION
517: tmp_pars[p] = tmp_a * (1 - DELTA);
518: call_gnuplot(tmp_pars, tmp_low);
519: #endif
520: for (k = 0; k < num_data; k++)
521: #ifdef TWO_SIDE_DIFFERENTIATION
522: dzda[k][p] = (tmp_high[k] - tmp_low[k]) / (2 * tmp_a * DELTA);
523: #else
524: dzda[k][p] = (tmp_high[k] - zfunc[k]) / (tmp_a * DELTA);
525: #endif
526: tmp_pars[p] = a[p];
527: }
528:
529: #ifdef TWO_SIDE_DIFFERENTIATION
530: free(tmp_low);
531: #endif
532: free(tmp_high);
533: free(tmp_pars);
534: }
535:
536:
537: /*****************************************************************
538: call internal gnuplot functions
539: *****************************************************************/
540: static void call_gnuplot(par, data)
541: double *par;
542: double *data;
543: {
544: int i;
545: struct value v;
546:
547: /* set parameters first */
548: for (i = 0; i < num_params; i++) {
549: (void) Gcomplex(&v, par[i], 0.0);
550: setvar(par_name[i], v);
551: }
552:
553: for (i = 0; i < num_data; i++) {
554: /* calculate fit-function value */
555: (void) Gcomplex(&func.dummy_values[0], fit_x[i], 0.0);
556: (void) Gcomplex(&func.dummy_values[1], fit_y[i], 0.0);
557: evaluate_at(func.at, &v);
558: if (undefined)
559: Eex("Undefined value during function evaluation");
560: data[i] = real(&v);
561: }
562: }
563:
564:
565: /*****************************************************************
566: handle user interrupts during fit
567: *****************************************************************/
568: static TBOOLEAN fit_interrupt()
569: {
570: while (TRUE) {
571: fputs("\n\n(S)top fit, (C)ontinue, (E)xecute FIT_SCRIPT: ", STANDARD);
572: switch (getc(stdin)) {
573:
574: case EOF:
575: case 's':
576: case 'S':
577: fputs("Stop.", STANDARD);
578: user_stop = TRUE;
579: return FALSE;
580:
581: case 'c':
582: case 'C':
583: fputs("Continue.", STANDARD);
584: return TRUE;
585:
586: case 'e':
587: case 'E':{
588: int i;
589: struct value v;
590: char *tmp;
591:
592: tmp = (fit_script != 0 && *fit_script) ? fit_script : DEFAULT_CMD;
593: fprintf(STANDARD, "executing: %s", tmp);
594: /* set parameters visible to gnuplot */
595: for (i = 0; i < num_params; i++) {
596: (void) Gcomplex(&v, a[i], 0.0);
597: setvar(par_name[i], v);
598: }
599: sprintf(input_line, tmp);
600: (void) do_line();
601: }
602: }
603: }
604: return TRUE;
605: }
606:
607:
608: /*****************************************************************
609: frame routine for the marquardt-fit
610: *****************************************************************/
611: static TBOOLEAN regress(a)
612: double a[];
613: {
614: double **covar, *dpar, **C, chisq, last_chisq, lambda;
615: int iter, i, j;
616: marq_res_t res;
617:
618: chisq = last_chisq = INFINITY;
619: C = matr(num_data + num_params, num_params);
620: lambda = -1; /* use sign as flag */
621: iter = 0; /* iteration counter */
622:
623: /* ctrlc now serves as Hotkey */
624: ctrlc_setup();
625:
626: /* Initialize internal variables and 1st chi-square check */
627: if ((res = marquardt(a, C, &chisq, &lambda)) == ERROR)
628: Eex("FIT: error occured during fit");
629: res = BETTER;
630:
631: show_fit(iter, chisq, chisq, a, lambda, STANDARD);
632: show_fit(iter, chisq, chisq, a, lambda, log_f);
633:
634: /* MAIN FIT LOOP: do the regression iteration */
635:
636: /* HBB 981118: initialize new variable 'user_break' */
637: user_stop = FALSE;
638:
639: do {
640: /*
641: * MSDOS defines signal(SIGINT) but doesn't handle it through
642: * real interrupts. So there remain cases in which a ctrl-c may
643: * be uncaught by signal. We must use kbhit() instead that really
644: * serves the keyboard interrupt (or write an own interrupt func
645: * which also generates #ifdefs)
646: *
647: * I hope that other OSes do it better, if not... add #ifdefs :-(
648: * EMX does not have kbhit.
649: *
650: * HBB: I think this can be enabled for DJGPP V2. SIGINT is actually
651: * handled there, AFAIK.
652: */
653: #if ((defined(MSDOS) || defined(DOS386)) && !defined(__EMX__)) || defined(ATARI) || defined(MTOS)
654: if (kbhit()) {
655: do {
656: getchx();
657: } while (kbhit());
658: ctrlc_flag = TRUE;
659: }
660: #endif
661:
662: if (ctrlc_flag) {
663: show_fit(iter, chisq, last_chisq, a, lambda, STANDARD);
664: ctrlc_flag = FALSE;
665: if (!fit_interrupt()) /* handle keys */
666: break;
667: }
668: if (res == BETTER) {
669: iter++;
670: last_chisq = chisq;
671: }
672: if ((res = marquardt(a, C, &chisq, &lambda)) == BETTER)
673: show_fit(iter, chisq, last_chisq, a, lambda, STANDARD);
674: } while ((res != ERROR)
675: && (lambda < MAX_LAMBDA)
676: && ((maxiter == 0) || (iter <= maxiter))
677: && (res == WORSE
678: || ((chisq > NEARLY_ZERO)
679: ? ((last_chisq - chisq) / chisq)
680: : (last_chisq - chisq)) > epsilon
681: )
682: );
683:
684: /* fit done */
685:
686: #if !defined(ATARI) && !defined(MTOS)
687: /* restore original SIGINT function */
688: interrupt_setup();
689: #endif
690:
691: /* HBB 970304: the maxiter patch: */
692: if ((maxiter > 0) && (iter > maxiter)) {
693: Dblf2("\nMaximum iteration count (%d) reached. Fit stopped.\n", maxiter);
694: } else if (user_stop) {
695: Dblf2("\nThe fit was stopped by the user after %d iterations.\n", iter);
696: } else {
697: Dblf2("\nAfter %d iterations the fit converged.\n", iter);
698: }
699:
700: Dblf2("final sum of squares of residuals : %g\n", chisq);
701: if (chisq > NEARLY_ZERO) {
702: Dblf2("rel. change during last iteration : %g\n\n", (chisq - last_chisq) / chisq);
703: } else {
704: Dblf2("abs. change during last iteration : %g\n\n", (chisq - last_chisq));
705: }
706:
707: if (res == ERROR)
708: Eex("FIT: error occured during fit");
709:
710: /* compute errors in the parameters */
711:
712: if (num_data == num_params) {
713: int i;
714:
715: Dblf("\nExactly as many data points as there are parameters.\n");
716: Dblf("In this degenerate case, all errors are zero by definition.\n\n");
717: Dblf("Final set of parameters \n");
718: Dblf("======================= \n\n");
719: for (i = 0; i < num_params; i++)
720: Dblf3("%-15.15s = %-15g\n", par_name[i], a[i]);
721: } else if (chisq < NEARLY_ZERO) {
722: int i;
723:
724: Dblf("\nHmmmm.... Sum of squared residuals is zero. Can't compute errors.\n\n");
725: Dblf("Final set of parameters \n");
726: Dblf("======================= \n\n");
727: for (i = 0; i < num_params; i++)
728: Dblf3("%-15.15s = %-15g\n", par_name[i], a[i]);
729: } else {
1.1.1.2 ! maekawa 730: Dblf2("degrees of freedom (ndf) : %d\n", num_data - num_params);
! 731: Dblf2("rms of residuals (stdfit) = sqrt(WSSR/ndf) : %g\n", sqrt(chisq / (num_data - num_params)));
! 732: Dblf2("variance of residuals (reduced chisquare) = WSSR/ndf : %g\n\n", chisq / (num_data - num_params));
1.1 maekawa 733:
734: /* get covariance-, Korrelations- and Kurvature-Matrix */
735: /* and errors in the parameters */
736:
737: /* compute covar[][] directly from C */
738: Givens(C, 0, 0, 0, num_data, num_params, 0);
739: /*printmatrix(C, num_params, num_params); */
740:
741: /* Use lower square of C for covar */
742: covar = C + num_data;
743: Invert_RtR(C, covar, num_params);
744: /*printmatrix(covar, num_params, num_params); */
745:
746: /* calculate unscaled parameter errors in dpar[]: */
747: dpar = vec(num_params);
748: for (i = 0; i < num_params; i++) {
749: /* FIXME: can this still happen ? */
750: if (covar[i][i] <= 0.0) /* HBB: prevent floating point exception later on */
751: Eex("Calculation error: non-positive diagonal element in covar. matrix");
752: dpar[i] = sqrt(covar[i][i]);
753: }
754:
755: /* transform covariances into correlations */
756: for (i = 0; i < num_params; i++) {
757: /* only lower triangle needs to be handled */
758: for (j = 0; j <= i; j++)
759: covar[i][j] /= dpar[i] * dpar[j];
760: }
761:
762: /* scale parameter errors based on chisq */
763: chisq = sqrt(chisq / (num_data - num_params));
764: for (i = 0; i < num_params; i++)
765: dpar[i] *= chisq;
766:
767: Dblf("Final set of parameters Asymptotic Standard Error\n");
768: Dblf("======================= ==========================\n\n");
769:
770: for (i = 0; i < num_params; i++) {
771: double temp =
772: (fabs(a[i]) < NEARLY_ZERO) ? 0.0 : fabs(100.0 * dpar[i] / a[i]);
773: Dblf6("%-15.15s = %-15g %-3.3s %-12.4g (%.4g%%)\n",
774: par_name[i], a[i], PLUSMINUS, dpar[i], temp);
775: }
776:
777: Dblf("\n\ncorrelation matrix of the fit parameters:\n\n");
778: Dblf(" ");
779:
780: for (j = 0; j < num_params; j++)
781: Dblf2("%-6.6s ", par_name[j]);
782:
783: Dblf("\n");
784: for (i = 0; i < num_params; i++) {
785: Dblf2("%-15.15s", par_name[i]);
786: for (j = 0; j <= i; j++) {
787: /* Only print lower triangle of symmetric matrix */
788: Dblf2("%6.3f ", covar[i][j]);
789: }
790: Dblf("\n");
791: }
792:
793: free(dpar);
794: }
795:
1.1.1.2 ! maekawa 796: /* HBB 990220: re-imported this snippet from older versions. Finally,
! 797: * some user noticed that it *is* necessary, after all. Not even
! 798: * Carsten Grammes himself remembered what it was for... :-(
! 799: * The thing is: the value of the last parameter is not reset to
! 800: * its original one after the derivatives have been calculated
! 801: */
! 802: /* restore last parameter's value (not done by calculate) */
! 803: {
! 804: struct value val;
! 805: Gcomplex (&val, a[num_params-1], 0.0);
! 806: setvar (par_name[num_params-1], val);
! 807: }
! 808:
1.1 maekawa 809: /* call destructor for allocated vars */
810: lambda = -2; /* flag value, meaning 'destruct!' */
811: (void) marquardt(a, C, &chisq, &lambda);
812:
813: free_matr(C);
814: return TRUE;
815: }
816:
817:
818: /*****************************************************************
819: display actual state of the fit
820: *****************************************************************/
821: static void show_fit(i, chisq, last_chisq, a, lambda, device)
822: int i;
823: double chisq;
824: double last_chisq;
825: double *a;
826: double lambda;
827: FILE *device;
828: {
829: int k;
830:
831: fprintf(device, "\n\n\
832: Iteration %d\n\
833: WSSR : %-15g delta(WSSR)/WSSR : %g\n\
834: delta(WSSR) : %-15g limit for stopping : %g\n\
835: lambda : %g\n\n%s parameter values\n\n",
836: i, chisq, chisq > NEARLY_ZERO ? (chisq - last_chisq) / chisq : 0.0,
837: chisq - last_chisq, epsilon, lambda,
838: (i > 0 ? "resultant" : "initial set of free"));
839: for (k = 0; k < num_params; k++)
840: fprintf(device, "%-15.15s = %g\n", par_name[k], a[k]);
841: }
842:
843:
844:
845: /*****************************************************************
846: is_empty: check for valid string entries
847: *****************************************************************/
848: static TBOOLEAN is_empty(s)
849: char *s;
850: {
851: while (*s == ' ' || *s == '\t' || *s == '\n')
852: s++;
853: return (TBOOLEAN) (*s == '#' || *s == '\0');
854: }
855:
856:
857: /*****************************************************************
858: get next word of a multi-word string, advance pointer
859: *****************************************************************/
860: char *get_next_word(s, subst)
861: char **s;
862: char *subst;
863: {
864: char *tmp = *s;
865:
866: while (*tmp == ' ' || *tmp == '\t' || *tmp == '=')
867: tmp++;
868: if (*tmp == '\n' || *tmp == '\0') /* not found */
869: return NULL;
870: if ((*s = strpbrk(tmp, " =\t\n")) == NULL)
871: *s = tmp + strlen(tmp);
872: *subst = **s;
873: *(*s)++ = '\0';
874: return tmp;
875: }
876:
877:
878: /*****************************************************************
879: check for variable identifiers
880: *****************************************************************/
881: static TBOOLEAN is_variable(s)
882: char *s;
883: {
884: while (*s != '\0') {
885: if (!isalnum((int) *s) && *s != '_')
886: return FALSE;
887: s++;
888: }
889: return TRUE;
890: }
891:
892: /*****************************************************************
893: first time settings
894: *****************************************************************/
895: void init_fit()
896: {
897: func.at = (struct at_type *) NULL; /* need to parse 1 time */
898: }
899:
900:
901: /*****************************************************************
902: Set a GNUPLOT user-defined variable
903: ******************************************************************/
904:
905: void setvar(varname, data)
906: char *varname;
907: struct value data;
908: {
909: register struct udvt_entry *udv_ptr = first_udv, *last = first_udv;
910:
911: /* check if it's already in the table... */
912:
913: while (udv_ptr) {
914: last = udv_ptr;
915: if (!strcmp(varname, udv_ptr->udv_name))
916: break;
917: udv_ptr = udv_ptr->next_udv;
918: }
919:
920: if (!udv_ptr) { /* generate new entry */
921: last->next_udv = udv_ptr = (struct udvt_entry *)
1.1.1.2 ! maekawa 922: gp_alloc(sizeof(struct udvt_entry), "fit setvar");
1.1 maekawa 923: udv_ptr->next_udv = NULL;
924: }
925: safe_strncpy(udv_ptr->udv_name, varname, sizeof(udv_ptr->udv_name));
926: udv_ptr->udv_value = data;
927: udv_ptr->udv_undef = FALSE;
928: }
929:
930:
931:
932: /*****************************************************************
933: Read INTGR Variable value, return 0 if undefined or wrong type
934: *****************************************************************/
935: int getivar(varname)
936: char *varname;
937: {
938: register struct udvt_entry *udv_ptr = first_udv;
939:
940: while (udv_ptr) {
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;
946: }
947: return 0; /* not in table */
948: }
949:
950:
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: *****************************************************************/
955: static double getdvar(varname)
956: char *varname;
957: {
958: register struct udvt_entry *udv_ptr = first_udv;
959:
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));
963:
964: /* get here => not found */
965: return 0;
966: }
967:
968: /*****************************************************************
969: like getdvar, but
970: - convert it from integer to real if necessary
971: - create it with value INITIAL_VALUE if not found or undefined
972: *****************************************************************/
973: static double createdvar(varname, value)
974: char *varname;
975: double value;
976: {
977: register struct udvt_entry *udv_ptr = first_udv;
978:
979: for (; udv_ptr; udv_ptr = udv_ptr->next_udv)
980: if (strcmp(varname, udv_ptr->udv_name) == 0) {
981: if (udv_ptr->udv_undef) {
982: udv_ptr->udv_undef = 0;
983: (void) Gcomplex(&udv_ptr->udv_value, value, 0.0);
984: } else if (udv_ptr->udv_value.type == INTGR) {
985: (void) Gcomplex(&udv_ptr->udv_value, (double) udv_ptr->udv_value.v.int_val, 0.0);
986: }
987: return real(&(udv_ptr->udv_value));
988: }
989: /* get here => not found */
990:
991: {
992: struct value tempval;
993: (void) Gcomplex(&tempval, value, 0.0);
994: setvar(varname, tempval);
995: }
996:
997: return value;
998: }
999:
1000:
1001: /*****************************************************************
1002: Split Identifier into path and filename
1003: *****************************************************************/
1004: static void splitpath(s, p, f)
1005: char *s;
1006: char *p;
1007: char *f;
1008: {
1009: register char *tmp;
1010: tmp = s + strlen(s) - 1;
1011: while (*tmp != '\\' && *tmp != '/' && *tmp != ':' && tmp - s >= 0)
1012: tmp--;
1013: strcpy(f, tmp + 1);
1014: safe_strncpy(p, s, (size_t) (tmp - s + 1));
1015: p[tmp - s + 1] = NUL;
1016: }
1017:
1018:
1019: /*****************************************************************
1020: write the actual parameters to start parameter file
1021: *****************************************************************/
1022: void update(pfile, npfile)
1023: char *pfile, *npfile;
1024: {
1025: char fnam[256], path[256], sstr[256], pname[64], tail[127], *s = sstr,
1026: *tmp, c;
1027: char ifilename[256], *ofilename;
1028: FILE *of, *nf;
1029: double pval;
1030:
1031:
1032: /* update pfile npfile:
1033: if npfile is a valid file name,
1034: take pfile as input file and
1035: npfile as output file
1036: */
1037: if (VALID_FILENAME(npfile)) {
1038: safe_strncpy(ifilename, pfile, sizeof(ifilename));
1039: ofilename = npfile;
1040: } else {
1041: #ifdef BACKUP_FILESYSTEM
1042: /* filesystem will keep original as previous version */
1043: safe_strncpy(ifilename, pfile, sizeof(ifilename));
1044: #else
1045: backup_file(ifilename, pfile); /* will Eex if it fails */
1046: #endif
1047: ofilename = pfile;
1048: }
1049:
1050: /* split into path and filename */
1051: splitpath(ifilename, path, fnam);
1052: if (!(of = fopen(ifilename, "r")))
1053: Eex2("parameter file %s could not be read", ifilename);
1054:
1055: if (!(nf = fopen(ofilename, "w")))
1056: Eex2("new parameter file %s could not be created", ofilename);
1057:
1058: while (fgets(s = sstr, sizeof(sstr), of) != NULL) {
1059:
1060: if (is_empty(s)) {
1061: fputs(s, nf); /* preserve comments */
1062: continue;
1063: }
1064: if ((tmp = strchr(s, '#')) != NULL) {
1065: safe_strncpy(tail, tmp, sizeof(tail));
1066: *tmp = NUL;
1067: } else
1068: strcpy(tail, "\n");
1069:
1070: tmp = get_next_word(&s, &c);
1071: if (!is_variable(tmp) || strlen(tmp) > MAX_ID_LEN) {
1072: (void) fclose(nf);
1073: (void) fclose(of);
1074: Eex2("syntax error in parameter file %s", fnam);
1075: }
1076: safe_strncpy(pname, tmp, sizeof(pname));
1077: /* next must be '=' */
1078: if (c != '=') {
1079: tmp = strchr(s, '=');
1080: if (tmp == NULL) {
1081: (void) fclose(nf);
1082: (void) fclose(of);
1083: Eex2("syntax error in parameter file %s", fnam);
1084: }
1085: s = tmp + 1;
1086: }
1087: tmp = get_next_word(&s, &c);
1088: if (!sscanf(tmp, "%lf", &pval)) {
1089: (void) fclose(nf);
1090: (void) fclose(of);
1091: Eex2("syntax error in parameter file %s", fnam);
1092: }
1093: if ((tmp = get_next_word(&s, &c)) != NULL) {
1094: (void) fclose(nf);
1095: (void) fclose(of);
1096: Eex2("syntax error in parameter file %s", fnam);
1097: }
1098: /* now modify */
1099:
1100: if ((pval = getdvar(pname)) == 0)
1101: pval = (double) getivar(pname);
1102:
1103: sprintf(sstr, "%g", pval);
1104: if (!strchr(sstr, '.') && !strchr(sstr, 'e'))
1105: strcat(sstr, ".0"); /* assure CMPLX-type */
1106:
1107: fprintf(nf, "%-15.15s = %-15.15s %s", pname, sstr, tail);
1108: }
1109:
1110: if (fclose(nf) || fclose(of))
1111: Eex("I/O error during update");
1112:
1113: }
1114:
1115:
1116: /*****************************************************************
1117: Backup a file by renaming it to something useful. Return
1118: the new name in tofile
1119: *****************************************************************/
1120: static void backup_file(tofile, fromfile)
1121: char *tofile;
1122: const char *fromfile;
1123: /* tofile must point to a char array[] or allocated data. See update() */
1124: {
1125: #if defined (WIN32) || defined(MSDOS) || defined(VMS)
1126: char *tmpn;
1127: #endif
1128:
1129: /* win32 needs to have two attempts at the rename, since it may
1130: * be running on win32s with msdos 8.3 names
1131: */
1132:
1133: /* first attempt, for all o/s other than MSDOS */
1134:
1135: #ifndef MSDOS
1136:
1137: strcpy(tofile, fromfile);
1138:
1139: #ifdef VMS
1140: /* replace all dots with _, since we will be adding the only
1141: * dot allowed in VMS names
1142: */
1143: while ((tmpn = strchr(tofile, '.')) != NULL)
1144: *tmpn = '_';
1145: #endif /*VMS */
1146:
1147: strcat(tofile, BACKUP_SUFFIX);
1148:
1149: if (rename(fromfile, tofile) == 0)
1150: return; /* hurrah */
1151: #endif
1152:
1153:
1154: #if defined (WIN32) || defined(MSDOS)
1155:
1156: /* first attempt for msdos. Second attempt for win32s */
1157:
1158: /* beware : strncpy is very dangerous since it does not guarantee
1159: * to terminate the string. Copy up to 8 characters. If exactly
1160: * chars were copied, the string is not terminated. If the
1161: * source string was shorter than 8 chars, no harm is done
1162: * (here) by writing to offset 8.
1163: */
1164: safe_strncpy(tofile, fromfile, 8);
1165: /* tofile[8] = NUL; */
1166:
1167: while ((tmpn = strchr(tofile, '.')) != NULL)
1168: *tmpn = '_';
1169:
1170: strcat(tofile, BACKUP_SUFFIX);
1171:
1172: if (rename(fromfile, tofile) == 0)
1173: return; /* success */
1174:
1175: #endif /* win32 || msdos */
1176:
1177: /* get here => rename failed. */
1178: Eex3("Could not rename file %s to %s", fromfile, tofile);
1179: }
1180:
1181:
1182: /*****************************************************************
1183: Interface to the classic gnuplot-software
1184: *****************************************************************/
1185:
1186: void do_fit()
1187: {
1188: TBOOLEAN autorange_x = 3, autorange_y = 3; /* yes */
1189: /* HBB 980401: new: z range specification */
1190: TBOOLEAN autorange_z = 3;
1191: double min_x, max_x; /* range to fit */
1192: double min_y, max_y; /* range to fit */
1193: /* HBB 980401: new: z range specification */
1194: double min_z, max_z; /* range to fit */
1195: int dummy_x = -1, dummy_y = -1; /* eg fit [u=...] [v=...] */
1196: /* HBB 981210: memorize position of possible third [ : ] spec: */
1197: int zrange_token = -1;
1198:
1199: int i;
1200: double v[4];
1201: double tmpd;
1202: time_t timer;
1203: int token1, token2, token3;
1204: char *tmp;
1205:
1206: /* first look for a restricted x fit range... */
1207:
1208: if (equals(c_token, "[")) {
1209: c_token++;
1210: if (isletter(c_token)) {
1211: if (equals(c_token + 1, "=")) {
1212: dummy_x = c_token;
1213: c_token += 2;
1214: }
1215: /* else parse it as a xmin expression */
1216: }
1217: autorange_x = load_range(FIRST_X_AXIS, &min_x, &max_x, autorange_x);
1218: if (!equals(c_token, "]"))
1219: int_error("']' expected", c_token);
1220: c_token++;
1221: }
1222: /* ... and y */
1223:
1224: if (equals(c_token, "[")) {
1225: c_token++;
1226: if (isletter(c_token)) {
1227: if (equals(c_token + 1, "=")) {
1228: dummy_y = c_token;
1229: c_token += 2;
1230: }
1231: /* else parse it as a ymin expression */
1232: }
1233: autorange_y = load_range(FIRST_Y_AXIS, &min_y, &max_y, autorange_y);
1234: if (!equals(c_token, "]"))
1235: int_error("']' expected", c_token);
1236: c_token++;
1237: }
1238: /* HBB 980401: new: allow restricting the z range as well */
1239: /* ... and z */
1240:
1241: if (equals(c_token, "[")) {
1242: zrange_token = c_token++;
1243: if (isletter(c_token))
1244: /* do *not* allow the z range being given with a variable name */
1245: int_error("Can't re-name dependent variable", c_token);
1246: autorange_z = load_range(FIRST_Z_AXIS, &min_z, &max_z, autorange_z);
1247: if (!equals(c_token, "]"))
1248: int_error("']' expected", c_token);
1249: c_token++;
1250: #if 0 /* HBB 981210: move this to a later point */
1251: } else {
1252: /* Just in case I muck up things below: make sure that the z
1253: * range is the same as the y range, if it didn't get specified
1254: */
1255: autorange_z = autorange_y;
1256: min_z = min_y;
1257: max_z = max_y;
1258: #endif
1259: }
1260:
1261:
1262: /* now compile the function */
1263:
1264: token1 = c_token;
1265:
1266: if (func.at) {
1267: free(func.at);
1268: func.at = NULL; /* in case perm_at() does int_error */
1269: }
1270: dummy_func = &func;
1271:
1272:
1273: /* set the dummy variable names */
1274:
1275: if (dummy_x >= 0)
1276: copy_str(c_dummy_var[0], dummy_x, MAX_ID_LEN);
1277: else
1278: strcpy(c_dummy_var[0], dummy_var[0]);
1279:
1280: if (dummy_y >= 0)
1281: copy_str(c_dummy_var[0], dummy_y, MAX_ID_LEN);
1282: else
1283: strcpy(c_dummy_var[1], dummy_var[1]);
1284:
1285: func.at = perm_at();
1286: dummy_func = NULL;
1287:
1288: token2 = c_token;
1289:
1290: /* use datafile module to parse the datafile and qualifiers */
1291:
1292: columns = df_open(4); /* up to 4 using specs allowed */
1293: if (columns == 1)
1294: int_error("Need 2 to 4 using specs", c_token);
1295:
1.1.1.2 ! maekawa 1296: /* The following patch was made by Remko Scharroo, 25-Mar-1999
! 1297: * We need to check if one of the columns is time data, like
! 1298: * in plot2d and plot3d */
! 1299:
! 1300: if (datatype[FIRST_X_AXIS] == TIME) {
! 1301: if (columns < 2)
! 1302: int_error("Need full using spec for x time data", c_token);
! 1303: df_timecol[0] = 1;
! 1304: }
! 1305: if (datatype[FIRST_Y_AXIS] == TIME) {
! 1306: if (columns < 1)
! 1307: int_error("Need using spec for y time data", c_token);
! 1308: df_timecol[1] = 1;
! 1309: }
! 1310: /* HBB 990326: added this check. Just in case some wants to fit
! 1311: * time/date data depending on two other variables ... */
! 1312: if (datatype[FIRST_Z_AXIS] == TIME) {
! 1313: if (columns < 4)
! 1314: int_error("Need full using spec for z time data", c_token);
! 1315: else
! 1316: df_timecol[2] = 1;
! 1317: }
! 1318: /* End of patch by Remko Scharroo */
! 1319:
1.1 maekawa 1320: /* HBB 980401: if this is a single-variable fit, we shouldn't have
1321: * allowed a variable name specifier for 'y': */
1322: if ((dummy_y >= 0) && (columns < 4))
1323: int_error("Can't re-name 'y' in a one-variable fit", dummy_y);
1324:
1325: /* HBB 981210: two range specs mean different things, depending
1326: * on wether this is a 2D or 3D fit */
1.1.1.2 ! maekawa 1327: if (columns < 4) {
! 1328: if (zrange_token != -1)
! 1329: int_error("Three range-specs not allowed in on-variable fit", zrange_token);
! 1330: else {
! 1331: /* 2D fit, 2 ranges: second range is for *z*, not y: */
! 1332: autorange_z = autorange_y;
! 1333: if (autorange_y & 1)
! 1334: min_z = min_y;
! 1335: if (autorange_y & 2)
! 1336: max_z = max_y;
! 1337: }
1.1 maekawa 1338: }
1339:
1340: /* defer actually reading the data until we have parsed the rest
1341: * of the line */
1342:
1343: token3 = c_token;
1344:
1345: tmpd = getdvar(FITLIMIT); /* get epsilon if given explicitly */
1346: if (tmpd < 1.0 && tmpd > 0.0)
1347: epsilon = tmpd;
1348:
1349: /* HBB 970304: maxiter patch */
1350: maxiter = getivar(FITMAXITER);
1351:
1352: /* get startup value for lambda, if given */
1353: tmpd = getdvar(FITSTARTLAMBDA);
1354:
1355: if (tmpd > 0.0) {
1356: startup_lambda = tmpd;
1357: printf("Lambda Start value set: %g\n", startup_lambda);
1358: }
1359:
1360: /* get lambda up/down factor, if given */
1361: tmpd = getdvar(FITLAMBDAFACTOR);
1362:
1363: if (tmpd > 0.0) {
1364: lambda_up_factor = lambda_down_factor = tmpd;
1365: printf("Lambda scaling factors reset: %g\n", lambda_up_factor);
1366: }
1367: *fit_script = NUL;
1368: if ((tmp = getenv(FITSCRIPT)) != NULL)
1369: strcpy(fit_script, tmp);
1370:
1371: tmp = getenv(GNUFITLOG); /* open logfile */
1372: if (tmp != NULL) {
1373: char *tmp2 = &tmp[strlen(tmp) - 1];
1374: if (*tmp2 == '/' || *tmp2 == '\\')
1375: sprintf(logfile, "%s%s", tmp, logfile);
1376: else
1377: strcpy(logfile, tmp);
1378: }
1379: if (!log_f && /* div */ !(log_f = fopen(logfile, "a")))
1380: Eex2("could not open log-file %s", logfile);
1381: fputs("\n\n*******************************************************************************\n", log_f);
1382: (void) time(&timer);
1383: fprintf(log_f, "%s\n\n", ctime(&timer));
1384: {
1385: char *line = NULL;
1386:
1387: m_capture(&line, token2, token3 - 1);
1388: fprintf(log_f, "FIT: data read from %s\n", line);
1389: free(line);
1390: }
1391:
1392: max_data = MAX_DATA;
1393: fit_x = vec(max_data); /* start with max. value */
1394: fit_y = vec(max_data);
1395: fit_z = vec(max_data);
1396:
1397: /* first read in experimental data */
1398:
1399: err_data = vec(max_data);
1400: num_data = 0;
1401:
1402: while ((i = df_readline(v, 4)) != EOF) {
1403: if (num_data >= max_data) {
1404: /* increase max_data by factor of 1.5 */
1405: max_data = (max_data * 3) / 2;
1406: if (0
1407: || !redim_vec(&fit_x, max_data)
1408: || !redim_vec(&fit_y, max_data)
1409: || !redim_vec(&fit_z, max_data)
1410: || !redim_vec(&err_data, max_data)
1411: ) {
1412: /* Some of the reallocations went bad: */
1413: df_close();
1414: Eex2("Out of memory in fit: too many datapoints (%d)?", max_data);
1415: } else {
1416: /* Just so we know that the routine is at work: */
1417: fprintf(STANDARD, "Max. number of data points scaled up to: %d\n",
1418: max_data);
1419: }
1420: }
1421: switch (i) {
1422: case DF_UNDEFINED:
1423: case DF_FIRST_BLANK:
1424: case DF_SECOND_BLANK:
1425: continue;
1426: case 0:
1427: Eex2("bad data on line %d of datafile", df_line_number);
1428: break;
1429: case 1: /* only z provided */
1430: v[2] = v[0];
1431: v[0] = (double) df_datum;
1432: break;
1433: case 2: /* x and z */
1434: v[2] = v[1];
1435: break;
1436:
1437: /* only if the explicitly asked for 4 columns do we
1438: * do a 3d fit. (We can get here if they didn't
1439: * specify a using spec, and the file has 4 columns
1440: */
1441: case 4: /* x, y, z, error */
1442: if (columns == 4)
1443: break;
1444: /* else fall through */
1445: case 3: /* x, z, error */
1446: v[3] = v[2]; /* error */
1447: v[2] = v[1]; /* z */
1448: break;
1449:
1450: }
1451:
1452: /* skip this point if it is out of range */
1453: if (!(autorange_x & 1) && (v[0] < min_x))
1454: continue;
1455: if (!(autorange_x & 2) && (v[0] > max_x))
1456: continue;
1457: /* HBB 980401: yrange is only relevant for 3d-fits! */
1458: if (columns == 4) {
1459: if (!(autorange_y & 1) && (v[1] < min_y))
1460: continue;
1461: if (!(autorange_y & 2) && (v[1] > max_y))
1462: continue;
1463: }
1464: /* HBB 980401: check *z* range for all fits */
1465: if (!(autorange_z & 1) && (v[2] < min_z))
1466: continue;
1467: if (!(autorange_z & 2) && (v[2] > max_z))
1468: continue;
1469:
1470: fit_x[num_data] = v[0];
1471: fit_y[num_data] = v[1];
1472: fit_z[num_data] = v[2];
1473:
1474: /* we only use error if _explicitly_ asked for by a
1475: * using spec
1476: */
1477: err_data[num_data++] = (columns > 2) ? v[3] : 1;
1478: }
1479: df_close();
1480:
1481: if (num_data <= 1)
1482: Eex("No data to fit ");
1483:
1484: /* now resize fields to actual length: */
1485: redim_vec(&fit_x, num_data);
1486: redim_vec(&fit_y, num_data);
1487: redim_vec(&fit_z, num_data);
1488: redim_vec(&err_data, num_data);
1489:
1490: fprintf(log_f, " #datapoints = %d\n", num_data);
1491:
1492: if (columns < 3)
1493: fputs(" residuals are weighted equally (unit weight)\n\n", log_f);
1494:
1495: {
1496: char *line = NULL;
1497:
1498: m_capture(&line, token1, token2 - 1);
1499: fprintf(log_f, "function used for fitting: %s\n", line);
1500: free(line);
1501: }
1502:
1503: /* read in parameters */
1504:
1505: max_params = MAX_PARAMS; /* HBB 971023: make this resizeable */
1506:
1507: if (!equals(c_token, "via"))
1508: int_error("Need via and either parameter list or file", c_token);
1509:
1510: a = vec(max_params);
1511: par_name = (fixstr *) gp_alloc((max_params + 1) * sizeof(fixstr), "fit param");
1512: num_params = 0;
1513:
1514: if (isstring(++c_token)) { /* It's a parameter *file* */
1515: TBOOLEAN fixed;
1516: double tmp_par;
1517: char c, *s;
1.1.1.2 ! maekawa 1518: char sstr[MAX_LINE_LEN+1];
1.1 maekawa 1519: FILE *f;
1520:
1521: quote_str(sstr, c_token++, MAX_LINE_LEN);
1522:
1523: /* get parameters and values out of file and ignore fixed ones */
1524:
1525: fprintf(log_f, "fitted parameters and initial values from file: %s\n\n", sstr);
1526: if (!(f = fopen(sstr, "r")))
1527: Eex2("could not read parameter-file %s", sstr);
1528: while (TRUE) {
1529: if (!fgets(s = sstr, (int) sizeof(sstr), f)) /* EOF found */
1530: break;
1531: if ((tmp = strstr(s, FIXED)) != NULL) { /* ignore fixed params */
1532: *tmp = NUL;
1533: fprintf(log_f, "FIXED: %s\n", s);
1534: fixed = TRUE;
1535: } else
1536: fixed = FALSE;
1537: if ((tmp = strchr(s, '#')) != NULL)
1538: *tmp = NUL;
1539: if (is_empty(s))
1540: continue;
1541: tmp = get_next_word(&s, &c);
1542: if (!is_variable(tmp) || strlen(tmp) > MAX_ID_LEN) {
1543: (void) fclose(f);
1544: Eex("syntax error in parameter file");
1545: }
1546: safe_strncpy(par_name[num_params], tmp, sizeof(fixstr));
1547: /* next must be '=' */
1548: if (c != '=') {
1549: tmp = strchr(s, '=');
1550: if (tmp == NULL) {
1551: (void) fclose(f);
1552: Eex("syntax error in parameter file");
1553: }
1554: s = tmp + 1;
1555: }
1556: tmp = get_next_word(&s, &c);
1557: if (sscanf(tmp, "%lf", &tmp_par) != 1) {
1558: (void) fclose(f);
1559: Eex("syntax error in parameter file");
1560: }
1561: /* make fixed params visible to GNUPLOT */
1562: if (fixed) {
1563: struct value tempval;
1564: (void) Gcomplex(&tempval, tmp_par, 0.0);
1565: /* use parname as temp */
1566: setvar(par_name[num_params], tempval);
1567: } else {
1568: if (num_params >= max_params) {
1569: max_params = (max_params * 3) / 2;
1570: if (0
1571: || !redim_vec(&a, max_params)
1572: || !(par_name = gp_realloc(par_name, (max_params + 1) * sizeof(fixstr), "fit param resize"))
1573: ) {
1574: (void) fclose(f);
1575: Eex("Out of memory in fit: too many parameters?");
1576: }
1577: }
1578: a[num_params++] = tmp_par;
1579: }
1580:
1581: if ((tmp = get_next_word(&s, &c)) != NULL) {
1582: (void) fclose(f);
1583: Eex("syntax error in parameter file");
1584: }
1585: }
1586: (void) fclose(f);
1587:
1588: } else {
1589: /* not a string after via: it's a variable listing */
1590:
1591: fputs("fitted parameters initialized with current variable values\n\n", log_f);
1592: do {
1593: if (!isletter(c_token))
1594: Eex("no parameter specified");
1595: capture(par_name[num_params], c_token, c_token, (int) sizeof(par_name[0]));
1596: if (num_params >= max_params) {
1597: max_params = (max_params * 3) / 2;
1598: if (0
1599: || !redim_vec(&a, max_params)
1600: || !(par_name = gp_realloc(par_name, (max_params + 1) * sizeof(fixstr), "fit param resize"))
1601: ) {
1602: Eex("Out of memory in fit: too many parameters?");
1603: }
1604: }
1605: /* create variable if it doesn't exist */
1606: a[num_params] = createdvar(par_name[num_params], INITIAL_VALUE);
1607: ++num_params;
1608: } while (equals(++c_token, ",") && ++c_token);
1609: }
1610:
1611: redim_vec(&a, num_params);
1612: par_name = (fixstr *) gp_realloc(par_name, (num_params + 1) * sizeof(fixstr), "fit param");
1613:
1614: if (num_data < num_params)
1615: Eex("Number of data points smaller than number of parameters");
1616:
1617: /* avoid parameters being equal to zero */
1618: for (i = 0; i < num_params; i++)
1619: if (a[i] == 0)
1620: a[i] = NEARLY_ZERO;
1621:
1622: (void) regress(a);
1623:
1624: (void) fclose(log_f);
1625: log_f = NULL;
1626: free(fit_x);
1627: free(fit_y);
1628: free(fit_z);
1629: free(err_data);
1630: free(a);
1631: free(par_name);
1632: if (func.at) {
1633: free(func.at); /* release perm. action table */
1634: func.at = (struct at_type *) NULL;
1635: }
1636: safe_strncpy(last_fit_command, input_line, sizeof(last_fit_command));
1637: }
1638:
1639: #if defined(ATARI) || defined(MTOS)
1640: static int kbhit()
1641: {
1642: fd_set rfds;
1643: struct timeval timeout;
1644:
1645: timeout.tv_sec = 0;
1646: timeout.tv_usec = 0;
1647: rfds = (fd_set) (1L << fileno(stdin));
1648: return ((select(0, &rfds, NULL, NULL, &timeout) > 0) ? 1 : 0);
1649: }
1650: #endif
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to fit.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gnuplot