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