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