Back to index

plt-scheme  4.2.1
fit.c
Go to the documentation of this file.
00001 /*  NOTICE: Change of Copyright Status
00002  *
00003  *  The author of this module, Carsten Grammes, has expressed in
00004  *  personal email that he has no more interest in this code, and
00005  *  doesn't claim any copyright. He has agreed to put this module
00006  *  into the public domain.
00007  *
00008  *  Lars Hecking  15-02-1999
00009  */
00010 
00011 /*
00012  *    Nonlinear least squares fit according to the
00013  *      Marquardt-Levenberg-algorithm
00014  *
00015  *      added as Patch to Gnuplot (v3.2 and higher)
00016  *      by Carsten Grammes
00017  *      Experimental Physics, University of Saarbruecken, Germany
00018  *
00019  *      Internet address: cagr@rz.uni-sb.de
00020  *
00021  *      Copyright of this module:  1993, 1998  Carsten Grammes
00022  *
00023  *      Permission to use, copy, and distribute this software and its
00024  *      documentation for any purpose with or without fee is hereby granted,
00025  *      provided that the above copyright notice appear in all copies and
00026  *      that both that copyright notice and this permission notice appear
00027  *      in supporting documentation.
00028  *
00029  *      This software is provided "as is" without express or implied warranty.
00030  *
00031  *      930726:     Recoding of the Unix-like raw console I/O routines by:
00032  *                  Michele Marziani (marziani@ferrara.infn.it)
00033  * drd: start unitialised variables at 1 rather than NEARLY_ZERO
00034  *  (fit is more likely to converge if started from 1 than 1e-30 ?)
00035  *
00036  * HBB (Broeker@physik.rwth-aachen.de) : fit didn't calculate the errors
00037  * in the 'physically correct' (:-) way, if a third data column containing
00038  * the errors (or 'uncertainties') of the input data was given. I think
00039  * I've fixed that, but I'm not sure I really understood the M-L-algo well
00040  * enough to get it right. I deduced my change from the final steps of the
00041  * equivalent algorithm for the linear case, which is much easier to
00042  * understand. (I also made some minor, mostly cosmetic changes)
00043  *
00044  * HBB (again): added error checking for negative covar[i][i] values and
00045  * for too many parameters being specified.
00046  *
00047  * drd: allow 3d fitting. Data value is now called fit_z internally,
00048  * ie a 2d fit is z vs x, and a 3d fit is z vs x and y.
00049  *
00050  * Lars Hecking : review update command, for VMS in particular, where
00051  * it is not necessary to rename the old file.
00052  *
00053  * HBB, 971023: lifted fixed limit on number of datapoints, and number
00054  * of parameters.
00055  */
00056 
00057 
00058 #define FIT_MAIN
00059 
00060 #define NULL 0
00061 
00062 //#include <scheme.h>
00063 
00064 
00065 
00066 #include "matrix.h"
00067 #include "fit.h"
00068 #include <math.h>
00069 
00070 /* #define STANDARD stderr */
00071 
00072 
00073 enum marq_res {
00074     OK, ERROR, BETTER, WORSE
00075 };
00076 typedef enum marq_res marq_res_t;
00077 
00078 #ifdef INFINITY
00079 # undef INFINITY
00080 #endif
00081 
00082 #define INFINITY    1e30
00083 #define NEARLY_ZERO 1e-30
00084 
00085 
00086 /* Relative change for derivatives */
00087 #define DELTA     0.001
00088 
00089 #define MAX_DATA    2048
00090 #define MAX_PARAMS  32
00091 #define MAX_LAMBDA  1e20
00092 #define MIN_LAMBDA  1e-20
00093 #define LAMBDA_UP_FACTOR 10
00094 #define LAMBDA_DOWN_FACTOR 10
00095 
00096 #define TBOOLEAN int
00097 
00098 # define PLUSMINUS   "+/-"
00099 
00100 
00101 
00102 /* HBB 971023: new, allow for dynamic adjustment of these: */
00103 static int max_data;
00104 static int max_params;
00105 
00106 static double epsilon = 1e-5;      /* convergence limit */
00107 static int maxiter = 0;            /* HBB 970304: maxiter patch */
00108 
00109 static char *FIXED = "# FIXED";
00110 static char *GNUFITLOG = "FIT_LOG";
00111 static char *FITLIMIT = "FIT_LIMIT";
00112 static char *FITSTARTLAMBDA = "FIT_START_LAMBDA";
00113 static char *FITLAMBDAFACTOR = "FIT_LAMBDA_FACTOR";
00114 static char *FITMAXITER = "FIT_MAXITER";  /* HBB 970304: maxiter patch */
00115 static char *FITSCRIPT = "FIT_SCRIPT";
00116 static char *DEFAULT_CMD = "replot";      /* if no fitscript spec. */
00117 
00118 
00119 static int num_data, num_params;
00120 static int columns;
00121 static double  *fit_x;
00122 static double    *fit_y;
00123 static double  *fit_z ;
00124 static double   *err_data;
00125 static double *a;
00126 
00127 
00128 /* static fixstr * par_name; */
00129 
00130 static double startup_lambda = 0;
00131 static double lambda_down_factor = LAMBDA_DOWN_FACTOR;
00132 static double lambda_up_factor = LAMBDA_UP_FACTOR;
00133 
00134 
00135 static void * current_fun;
00136 
00137 
00138 /*****************************************************************
00139                       internal vars to store results of fit
00140 *****************************************************************/
00141 
00142 double rms = 0;
00143 double varience = 0;
00144 double *asym_error;
00145 double *asym_error_percent;
00146 
00147 MZ_DLLEXPORT
00148 double get_rms()
00149 {return rms;}
00150 
00151 MZ_DLLEXPORT
00152 double get_varience()
00153 {return varience;}
00154 
00155 MZ_DLLEXPORT
00156 double * get_asym_error()
00157 {return asym_error;}
00158 
00159 MZ_DLLEXPORT
00160 double * get_asym_error_percent()
00161 {return asym_error_percent;}
00162 
00163 
00164 /*****************************************************************
00165                       internal Prototypes
00166 *****************************************************************/
00167 
00168 /*static void printmatrix __PROTO((double **C, int m, int n)); */
00169 static void print_matrix_and_vectors (double **C, double *d, double *r, int m, int n);
00170 static marq_res_t marquardt (double a[], double **alpha, double *chisq, 
00171                           double *lambda); 
00172 static TBOOLEAN analyze (double a[], double **alpha, double beta[], 
00173                       double *chisq);
00174 static void calculate (double *zfunc, double **dzda, double a[]);
00175 static void call_scheme (double *par, double *data);
00176 
00177 static TBOOLEAN regress (double a[]);
00178 //static void show_fit (int i, double chisq, double last_chisq, double *a, 
00179 //                         double lambda, FILE * device); 
00180 
00181 
00182 /*****************************************************************
00183     New utility routine: print a matrix (for debugging the alg.)
00184 *****************************************************************/
00185 static void printmatrix(C, m, n)
00186 double **C;
00187 int m, n;
00188 {
00189     int i, j;
00190 
00191     for (i = 0; i < m; i++) {
00192        for (j = 0; j < n - 1; j++);
00193          /* Dblf2("%.8g |", C[i][j]); */
00194        /* Dblf2("%.8g\n", C[i][j]); */
00195     }
00196     /* Dblf("\n"); */
00197 }
00198 
00199 /**************************************************************************
00200     Yet another debugging aid: print matrix, with diff. and residue vector
00201 **************************************************************************/
00202 static void print_matrix_and_vectors(C, d, r, m, n)
00203 double **C;
00204 double *d, *r;
00205 int m, n;
00206 {
00207     int i, j;
00208 
00209     for (i = 0; i < m; i++) {
00210        for (j = 0; j < n; j++);
00211          /* Dblf2("%8g ", C[i][j]); */
00212          /* Dblf3("| %8g | %8g\n", d[i], r[i]); */
00213     }
00214     /* Dblf("\n"); */
00215 }
00216 
00217 
00218 /*****************************************************************
00219     Marquardt's nonlinear least squares fit
00220 *****************************************************************/
00221 static marq_res_t marquardt(a, C, chisq, lambda)
00222 double a[];
00223 double **C;
00224 double *chisq;
00225 double *lambda;
00226 {
00227     int i, j;
00228     static double *da = 0,  /* delta-step of the parameter */
00229     *temp_a = 0,            /* temptative new params set   */
00230     *d = 0, *tmp_d = 0, **tmp_C = 0, *residues = 0;
00231     double tmp_chisq;
00232 
00233     /* Initialization when lambda == -1 */
00234 
00235     if (*lambda == -1) {    /* Get first chi-square check */
00236        TBOOLEAN analyze_ret;
00237 
00238        temp_a = vec(num_params);
00239        d = vec(num_data + num_params);
00240        tmp_d = vec(num_data + num_params);
00241        da = vec(num_params);
00242        residues = vec(num_data + num_params);
00243        tmp_C = matr(num_data + num_params, num_params);
00244 
00245        analyze_ret = analyze(a, C, d, chisq);
00246 
00247        /* Calculate a useful startup value for lambda, as given by Schwarz */
00248        /* FIXME: this is doesn't turn out to be much better, really... */
00249        if (startup_lambda != 0)
00250            *lambda = startup_lambda;
00251        else {
00252            *lambda = 0;
00253            for (i = 0; i < num_data; i++)
00254               for (j = 0; j < num_params; j++)
00255                   *lambda += C[i][j] * C[i][j];
00256            *lambda = sqrt(*lambda / num_data / num_params);
00257        }
00258 
00259        /* Fill in the lower square part of C (the diagonal is filled in on
00260           each iteration, see below) */
00261        for (i = 0; i < num_params; i++)
00262            for (j = 0; j < i; j++)
00263               C[num_data + i][j] = 0, C[num_data + j][i] = 0;
00264        /* printmatrix(C, num_data+num_params, num_params); */
00265        return analyze_ret ? OK : ERROR;
00266     }
00267     /* once converged, free dynamic allocated vars */
00268 
00269     if (*lambda == -2) {
00270        return OK;
00271     }
00272     /* Givens calculates in-place, so make working copies of C and d */
00273 
00274     for (j = 0; j < num_data + num_params; j++)
00275        memcpy(tmp_C[j], C[j], num_params * sizeof(double));
00276     memcpy(tmp_d, d, num_data * sizeof(double));
00277 
00278     /* fill in additional parts of tmp_C, tmp_d */
00279 
00280     for (i = 0; i < num_params; i++) {
00281        /* fill in low diag. of tmp_C ... */
00282        tmp_C[num_data + i][i] = *lambda;
00283        /* ... and low part of tmp_d */
00284        tmp_d[num_data + i] = 0;
00285     }
00286     /* printmatrix(tmp_C, num_data+num_params, num_params); */
00287 
00288     /* FIXME: residues[] isn't used at all. Why? Should it be used? */
00289 
00290     Givens(tmp_C, tmp_d, da, residues, num_params + num_data, num_params, 1);
00291     /*print_matrix_and_vectors (tmp_C, tmp_d, residues,
00292        num_params+num_data, num_params); */
00293 
00294     /* check if trial did ameliorate sum of squares */
00295 
00296     for (j = 0; j < num_params; j++)
00297        temp_a[j] = a[j] + da[j];
00298 
00299     if (!analyze(temp_a, tmp_C, tmp_d, &tmp_chisq)) {
00300        /* FIXME: will never be reached: always returns TRUE */
00301        return ERROR;
00302     }
00303 
00304     if (tmp_chisq < *chisq) {      /* Success, accept new solution */
00305        if (*lambda > MIN_LAMBDA) {
00306          /* (void) putc('/', stderr); */
00307           *lambda /= lambda_down_factor;
00308        }
00309        *chisq = tmp_chisq;
00310        for (j = 0; j < num_data; j++) {
00311            memcpy(C[j], tmp_C[j], num_params * sizeof(double));
00312            d[j] = tmp_d[j];
00313        }
00314        for (j = 0; j < num_params; j++)
00315            a[j] = temp_a[j];
00316        return BETTER;
00317     } else {                /* failure, increase lambda and return */
00318       /* (void) putc('*', stderr); */
00319       *lambda *= lambda_up_factor;
00320       return WORSE;
00321     }
00322 }
00323 
00324 
00325 /* FIXME: in the new code, this function doesn't really do enough to be
00326  * useful. Maybe it ought to be deleted, i.e. integrated with
00327  * calculate() ?
00328  */
00329 /*****************************************************************
00330     compute chi-square and numeric derivations
00331 *****************************************************************/
00332 static TBOOLEAN analyze(a, C, d, chisq)
00333 double a[];
00334 double **C;
00335 double d[];
00336 double *chisq;
00337 {
00338 /*
00339  *  used by marquardt to evaluate the linearized fitting matrix C
00340  *  and vector d, fills in only the top part of C and d
00341  *  I don't use a temporary array zfunc[] any more. Just use
00342  *  d[] instead.
00343  */
00344     int i, j;
00345 
00346     *chisq = 0;
00347     calculate(d, C, a);
00348 
00349     for (i = 0; i < num_data; i++) {
00350        /* note: order reversed, as used by Schwarz */
00351        d[i] = (d[i] - fit_z[i]) / err_data[i];
00352        *chisq += d[i] * d[i];
00353        for (j = 0; j < num_params; j++)
00354            C[i][j] /= err_data[i];
00355     }
00356     /* FIXME: why return a value that is always TRUE ? */
00357     return 1;
00358 }
00359 
00360 
00361 /* To use the more exact, but slower two-side formula, activate the
00362    following line: */
00363 
00364 #define TWO_SIDE_DIFFERENTIATION 
00365 
00366 /*****************************************************************
00367     compute function values and partial derivatives of chi-square
00368 *****************************************************************/
00369 static void calculate(zfunc, dzda, a)
00370 double *zfunc;
00371 double **dzda;
00372 double a[];
00373 {
00374     int k, p;
00375     double tmp_a;
00376     double *tmp_high, *tmp_pars;
00377 #ifdef TWO_SIDE_DIFFERENTIATION
00378     double *tmp_low;
00379 #endif
00380 
00381     tmp_high = vec(num_data);      /* numeric derivations */
00382 #ifdef TWO_SIDE_DIFFERENTIATION
00383     tmp_low = vec(num_data);
00384 #endif
00385     tmp_pars = vec(num_params);
00386 
00387     /* first function values */
00388 
00389     call_scheme(a, zfunc);
00390 
00391     /* then derivatives */
00392 
00393     for (p = 0; p < num_params; p++)
00394        tmp_pars[p] = a[p];
00395     for (p = 0; p < num_params; p++) {
00396        tmp_a = fabs(a[p]) < NEARLY_ZERO ? NEARLY_ZERO : a[p];
00397        tmp_pars[p] = tmp_a * (1 + DELTA);
00398        call_scheme(tmp_pars, tmp_high);
00399 #ifdef TWO_SIDE_DIFFERENTIATION
00400        tmp_pars[p] = tmp_a * (1 - DELTA);
00401        call_scheme(tmp_pars, tmp_low);
00402 #endif
00403        for (k = 0; k < num_data; k++)
00404 #ifdef TWO_SIDE_DIFFERENTIATION
00405            dzda[k][p] = (tmp_high[k] - tmp_low[k]) / (2 * tmp_a * DELTA);
00406 #else
00407            dzda[k][p] = (tmp_high[k] - zfunc[k]) / (tmp_a * DELTA);
00408 #endif
00409        tmp_pars[p] = a[p];
00410     }
00411 
00412 }
00413 
00414 
00415 /*****************************************************************
00416     evaluate the scheme function
00417 *****************************************************************/
00418 static void call_scheme(par, data)
00419 double *par;
00420 double *data;
00421 {
00422   int rators = 2 + num_params;
00423   double * rands =
00424     (double *) malloc(rators * sizeof(double));
00425 
00426   int i;
00427 
00428   /* set up the constant params */
00429   for(i = 0 ; i< num_params; i++) {
00430     rands[i+2] = par[i];
00431   }
00432 
00433   /* now calculate the function at the existing points */
00434   for (i = 0; i < num_data; i++) {
00435     rands[0] = fit_x[i];
00436     rands[1] = fit_y[i];
00437 
00438     data[i] = ((double (*) (int, double *) )current_fun) // ouch!
00439                (rators, rands);
00440   }
00441 
00442   free(rands);
00443 
00444 }
00445 
00446 /* /\***************************************************************** */
00447 /*     evaluate the scheme function */
00448 /* *****************************************************************\/ */
00449 /* static void call_scheme(par, data) */
00450 /* double *par; */
00451 /* double *data; */
00452 /* { */
00453 /*   int rators = 2 + num_params; */
00454 /*   Scheme_Object ** rands =  */
00455 /*     scheme_malloc(rators * sizeof(Scheme_Object)); */
00456 
00457 /*   int i; */
00458 
00459 /*   /\* set up the constant params *\/ */
00460 /*   for(i = 0 ; i< num_params; i++) { */
00461 /*     rands[i+2] = scheme_make_double(par[i]); */
00462 /*   } */
00463 
00464 /*   /\* now calculate the function at the existing points *\/ */
00465 /*     for (i = 0; i < num_data; i++) { */
00466 /*       rands[0] = scheme_make_double(fit_x[i]); */
00467 /*       rands[1] = scheme_make_double(fit_y[i]); */
00468 
00469 /*       data[i] = scheme_real_to_double(scheme_apply(current_fun, rators, rands)); */
00470 /*     } */
00471 /* } */
00472 
00473 /*****************************************************************
00474     Frame routine for the marquardt-fit
00475 *****************************************************************/
00476 static TBOOLEAN regress(a)
00477     double a[];
00478 {
00479   double **covar, *dpar, **C, chisq, last_chisq, lambda;
00480   int iter, i, j;
00481   marq_res_t res;
00482   
00483   chisq = last_chisq = INFINITY;
00484   C = matr(num_data + num_params, num_params);
00485   lambda = -1;          /* use sign as flag */
00486   iter = 0;             /* iteration counter  */
00487 
00488   /* Initialize internal variables and 1st chi-square check */
00489 
00490   if ((res = marquardt(a, C, &chisq, &lambda)) == ERROR)
00491     return 0; /* an error occurded */
00492 
00493   res = BETTER;
00494 
00495   /* show_fit(iter, chisq, chisq, a, lambda, STANDARD); */
00496 
00497   /* MAIN FIT LOOP: do the regression iteration */
00498 
00499   do {
00500     if (res == BETTER) {
00501       iter++;
00502       last_chisq = chisq;
00503     }
00504     if ((res = marquardt(a, C, &chisq, &lambda)) == BETTER)
00505       {};
00506     /* show_fit(iter, chisq, last_chisq, a, lambda, STANDARD); */
00507   } while ((res != ERROR)
00508           && (lambda < MAX_LAMBDA)
00509           && ((maxiter == 0) || (iter <= maxiter))
00510           && (res == WORSE
00511               || ((chisq > NEARLY_ZERO)
00512                  ? ((last_chisq - chisq) / chisq)
00513                  : (last_chisq - chisq)) > epsilon
00514               )
00515           );
00516 
00517     /* fit done */
00518 
00519   /* save all the info that was otherwise printed out */
00520 
00521   rms = sqrt(chisq / (num_data - num_params));
00522   varience = chisq / (num_data - num_params);
00523   asym_error = malloc (num_params * sizeof (double));  
00524   asym_error_percent = malloc (num_params * sizeof (double)) ;
00525 
00526   /* don't know what the following code does... */
00527 
00528   /* compute covar[][] directly from C */ 
00529   Givens(C, 0, 0, 0, num_data, num_params, 0); 
00530   covar = C + num_data; 
00531   Invert_RtR(C, covar, num_params);
00532   
00533   dpar = vec(num_params); 
00534   for (i = 0; i < num_params; i++) {
00535     /* FIXME: can this still happen ? */
00536       if (covar[i][i] <= 0.0)      /* HBB: prevent floating point exception later on */
00537        return 0;  /* Eex("Calculation error: non-positive diagonal element in covar. matrix"); */
00538     dpar[i] = sqrt(covar[i][i]);
00539   }  
00540 
00541   /* transform covariances into correlations */
00542   for (i = 0; i < num_params; i++) {
00543     /* only lower triangle needs to be handled */
00544       for (j = 0; j <= i; j++)
00545       covar[i][j] /= dpar[i] * dpar[j];
00546   }
00547 
00548   /* scale parameter errors based on chisq */
00549   chisq = sqrt(chisq / (num_data - num_params));
00550   for (i = 0; i < num_params; i++)
00551     dpar[i] *= chisq;
00552 
00553   for(i = 0; i< num_params; i++)
00554     {
00555       double temp = 
00556        (fabs(a[i]) < NEARLY_ZERO) ? 0.0 : fabs(100.0 * dpar[i] / a[i]);
00557       asym_error[i] = dpar[i];
00558       asym_error_percent[i] = temp;
00559     }
00560 
00561   return 1;
00562 
00563 
00564     /******** CRAP LEFT OVER FROM GNUPLOT ***********/
00565 
00566     /* HBB 970304: the maxiter patch: */
00567     /*
00568     if ((maxiter > 0) && (iter > maxiter)) {
00569        Dblf2("\nMaximum iteration count (%d) reached. Fit stopped.\n", maxiter);
00570     } else  {
00571        Dblf2("\nAfter %d iterations the fit converged.\n", iter);
00572     }
00573 
00574     Dblf2("final sum of squares of residuals : %g\n", chisq);
00575     if (chisq > NEARLY_ZERO) {
00576        Dblf2("rel. change during last iteration : %g\n\n", (chisq - last_chisq) / chisq);
00577     } else {
00578        Dblf2("abs. change during last iteration : %g\n\n", (chisq - last_chisq));
00579     }
00580 
00581     if (res == ERROR)
00582       // Eex("FIT: error occurred during fit");
00583     */
00584     /* compute errors in the parameters */
00585 
00586   /*   if (num_data == num_params) { */
00587 /*     int i; */
00588 
00589 /*     Dblf("\nExactly as many data points as there are parameters.\n"); */
00590 /*     Dblf("In this degenerate case, all errors are zero by definition.\n\n"); */
00591 /*     Dblf("Final set of parameters \n"); */
00592 /*     Dblf("======================= \n\n"); */
00593 /*     for (i = 0; i < num_params; i++) */
00594 /*         Dblf3("%-15.15s = %-15g\n", par_name[i], a[i]); */
00595 /*     } else if (chisq < NEARLY_ZERO) { */
00596 /*     int i; */
00597 
00598 /*     Dblf("\nHmmmm.... Sum of squared residuals is zero. Can't compute errors.\n\n"); */
00599 /*     Dblf("Final set of parameters \n"); */
00600 /*     Dblf("======================= \n\n"); */
00601 /*     for (i = 0; i < num_params; i++) */
00602 /*         Dblf3("%-15.15s = %-15g\n", par_name[i], a[i]); */
00603 /*     } else { */
00604 /*     Dblf2("degrees of freedom (ndf) : %d\n",  num_data - num_params); */
00605 /*     Dblf2("rms of residuals      (stdfit) = sqrt(WSSR/ndf)      : %g\n", sqrt(chisq / (num_data - num_params))); */
00606 /*     Dblf2("variance of residuals (reduced chisquare) = WSSR/ndf : %g\n\n", chisq / (num_data - num_params)); */
00607  
00608 /*     /\* get covariance-, Korrelations- and Kurvature-Matrix *\/ */
00609 /*     /\* and errors in the parameters                     *\/ */
00610 
00611 /*     /\* compute covar[][] directly from C *\/ */
00612 /*     Givens(C, 0, 0, 0, num_data, num_params, 0); */
00613 /*     /\*printmatrix(C, num_params, num_params); *\/ */
00614 
00615 /*     /\* Use lower square of C for covar *\/ */
00616 /*     covar = C + num_data; */
00617 /*     Invert_RtR(C, covar, num_params); */
00618 /*     /\*printmatrix(covar, num_params, num_params); *\/ */
00619 
00620 /*     /\* calculate unscaled parameter errors in dpar[]: *\/ */
00621 /*     dpar = vec(num_params); */
00622 /*     for (i = 0; i < num_params; i++) { */
00623 /*         /\* FIXME: can this still happen ? *\/ */
00624 /*         if (covar[i][i] <= 0.0) /\* HBB: prevent floating point exception later on *\/ */
00625 /*            Eex("Calculation error: non-positive diagonal element in covar. matrix"); */
00626 /*         dpar[i] = sqrt(covar[i][i]); */
00627 /*     } */
00628 
00629 /*     /\* transform covariances into correlations *\/ */
00630 /*     for (i = 0; i < num_params; i++) { */
00631 /*         /\* only lower triangle needs to be handled *\/ */
00632 /*         for (j = 0; j <= i; j++) */
00633 /*            covar[i][j] /= dpar[i] * dpar[j]; */
00634 /*     } */
00635 
00636 /*     /\* scale parameter errors based on chisq *\/ */
00637 /*     chisq = sqrt(chisq / (num_data - num_params)); */
00638 /*     for (i = 0; i < num_params; i++) */
00639 /*         dpar[i] *= chisq; */
00640 
00641 /*     Dblf("Final set of parameters            Asymptotic Standard Error\n"); */
00642 /*     Dblf("=======================            ==========================\n\n"); */
00643 
00644 /*     for (i = 0; i < num_params; i++) { */
00645 /*         double temp = */
00646 /*         (fabs(a[i]) < NEARLY_ZERO) ? 0.0 : fabs(100.0 * dpar[i] / a[i]); */
00647 /*         Dblf6("%-15.15s = %-15g  %-3.3s %-12.4g (%.4g%%)\n", */
00648 /*              par_name[i], a[i], PLUSMINUS, dpar[i], temp); */
00649 /*     } */
00650 
00651 /*     Dblf("\n\ncorrelation matrix of the fit parameters:\n\n"); */
00652 /*     Dblf("               "); */
00653 
00654 /*     for (j = 0; j < num_params; j++) */
00655 /*         Dblf2("%-6.6s ", par_name[j]); */
00656 
00657 /*     Dblf("\n"); */
00658 /*     for (i = 0; i < num_params; i++) { */
00659 /*         Dblf2("%-15.15s", par_name[i]); */
00660 /*         for (j = 0; j <= i; j++) { */
00661 /*            /\* Only print lower triangle of symmetric matrix *\/ */
00662 /*            Dblf2("%6.3f ", covar[i][j]); */
00663 /*         } */
00664 /*         Dblf("\n"); */
00665 /*     } */
00666 
00667 /*     free(dpar); */
00668 /*     } */
00669 
00670     return 1;
00671 }
00672 
00673 
00674 /*****************************************************************
00675     display actual state of the fit
00676 *****************************************************************/
00677 /* static void show_fit(i, chisq, last_chisq, a, lambda, device) */
00678 /* int i; */
00679 /* double chisq; */
00680 /* double last_chisq; */
00681 /* double *a; */
00682 /* double lambda; */
00683 /* FILE *device; */
00684 //{
00685   /*
00686     int k;
00687 
00688     fprintf(device, "\n\n\
00689 Iteration %d\n\
00690 WSSR        : %-15g   delta(WSSR)/WSSR   : %g\n\
00691 delta(WSSR) : %-15g   limit for stopping : %g\n\
00692 lambda   : %g\n\n%s parameter values\n\n",
00693       i, chisq, chisq > NEARLY_ZERO ? (chisq - last_chisq) / chisq : 0.0,
00694            chisq - last_chisq, epsilon, lambda,
00695            (i > 0 ? "resultant" : "initial set of free"));
00696     for (k = 0; k < num_params; k++)
00697        fprintf(device, "%-15.15s = %g\n", par_name[k], a[k]);
00698   */
00699 //}
00700 
00701 
00702 
00703 
00704 
00705 
00706 /*****************************************************************
00707     Interface to scheme
00708 *****************************************************************/
00709 MZ_DLLEXPORT
00710 double * do_fit(void * function,
00711               int n_values,
00712               double * x_values,
00713               double * y_values,
00714               double * z_values,
00715               double * errors,
00716               int n_parameters,
00717               double * parameters) {
00718 
00719   /* reset lambda and other parameters if desired */
00720   int i;
00721   current_fun = function;
00722 
00723   num_data = n_values;
00724   fit_x = x_values;
00725   fit_y = y_values;
00726   fit_z = z_values; /* value is stored in z */
00727   err_data = errors;
00728 
00729   a = parameters;
00730   num_params = n_parameters;
00731 
00732   /* redim_vec(&a, num_params); */
00733   /* par_name = (fixstr *) gp_realloc(par_name, (num_params + 1) * sizeof(fixstr), "fit param"); */
00734 
00735   /* avoid parameters being equal to zero */
00736   for (i = 0; i < num_params; i++) {
00737     if (a[i] == 0) {
00738       a[i] = NEARLY_ZERO; 
00739     }
00740   }
00741   
00742   if(regress(a)) {
00743     gc_cleanup();
00744     return a;
00745   }
00746   else { /* something went wrong */
00747     gc_cleanup(); 
00748     return NULL;
00749   }
00750 }