Annotation of OpenXM_contrib/gnuplot/fit.c, Revision 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