version 1.1, 2000/01/09 17:00:50 |
version 1.1.1.3, 2003/09/15 07:09:24 |
|
|
static char *RCSid = "$Id$"; |
static char *RCSid = "$Id$"; |
#endif |
#endif |
|
|
|
/* NOTICE: Change of Copyright Status |
|
* |
|
* The author of this module, Carsten Grammes, has expressed in |
|
* personal email that he has no more interest in this code, and |
|
* doesn't claim any copyright. He has agreed to put this module |
|
* into the public domain. |
|
* |
|
* Lars Hecking 15-02-1999 |
|
*/ |
|
|
/* |
/* |
* Nonlinear least squares fit according to the |
* Nonlinear least squares fit according to the |
* Marquardt-Levenberg-algorithm |
* Marquardt-Levenberg-algorithm |
Line 89 extern char c_dummy_var[MAX_NUM_VAR][MAX_ID_LEN+1]; |
|
Line 99 extern char c_dummy_var[MAX_NUM_VAR][MAX_ID_LEN+1]; |
|
extern int c_token; |
extern int c_token; |
extern int df_datum, df_line_number; |
extern int df_datum, df_line_number; |
|
|
|
/* following 2 external arrays are needed to use time data */ |
|
extern int datatype[]; |
|
extern int df_timecol[]; |
|
|
enum marq_res { |
enum marq_res { |
OK, ERROR, BETTER, WORSE |
OK, ERROR, BETTER, WORSE |
Line 606 static TBOOLEAN fit_interrupt() |
|
Line 619 static TBOOLEAN fit_interrupt() |
|
frame routine for the marquardt-fit |
frame routine for the marquardt-fit |
*****************************************************************/ |
*****************************************************************/ |
static TBOOLEAN regress(a) |
static TBOOLEAN regress(a) |
double a[]; |
double a[]; |
{ |
{ |
double **covar, *dpar, **C, chisq, last_chisq, lambda; |
double **covar, *dpar, **C, chisq, last_chisq, lambda; |
int iter, i, j; |
int iter, i, j; |
|
|
|
|
/* Initialize internal variables and 1st chi-square check */ |
/* Initialize internal variables and 1st chi-square check */ |
if ((res = marquardt(a, C, &chisq, &lambda)) == ERROR) |
if ((res = marquardt(a, C, &chisq, &lambda)) == ERROR) |
Eex("FIT: error occured during fit"); |
Eex("FIT: error occurred during fit"); |
res = BETTER; |
res = BETTER; |
|
|
show_fit(iter, chisq, chisq, a, lambda, STANDARD); |
show_fit(iter, chisq, chisq, a, lambda, STANDARD); |
|
|
} |
} |
|
|
if (res == ERROR) |
if (res == ERROR) |
Eex("FIT: error occured during fit"); |
Eex("FIT: error occurred during fit"); |
|
|
/* compute errors in the parameters */ |
/* compute errors in the parameters */ |
|
|
|
|
for (i = 0; i < num_params; i++) |
for (i = 0; i < num_params; i++) |
Dblf3("%-15.15s = %-15g\n", par_name[i], a[i]); |
Dblf3("%-15.15s = %-15g\n", par_name[i], a[i]); |
} else { |
} else { |
Dblf2("degrees of freedom (ndf) : %d\n", num_data - num_params); |
Dblf2("degrees of freedom (ndf) : %d\n", num_data - num_params); |
Dblf2("rms of residuals (stdfit) = sqrt(WSSR/ndf) : %g\n", sqrt(chisq / (num_data - num_params))); |
Dblf2("rms of residuals (stdfit) = sqrt(WSSR/ndf) : %g\n", sqrt(chisq / (num_data - num_params))); |
Dblf2("variance of residuals (reduced chisquare) = WSSR/ndf : %g\n\n", chisq / (num_data - num_params)); |
Dblf2("variance of residuals (reduced chisquare) = WSSR/ndf : %g\n\n", chisq / (num_data - num_params)); |
|
|
/* get covariance-, Korrelations- and Kurvature-Matrix */ |
/* get covariance-, Korrelations- and Kurvature-Matrix */ |
/* and errors in the parameters */ |
/* and errors in the parameters */ |
|
|
free(dpar); |
free(dpar); |
} |
} |
|
|
|
/* HBB 990220: re-imported this snippet from older versions. Finally, |
|
* some user noticed that it *is* necessary, after all. Not even |
|
* Carsten Grammes himself remembered what it was for... :-( |
|
* The thing is: the value of the last parameter is not reset to |
|
* its original one after the derivatives have been calculated |
|
*/ |
|
/* restore last parameter's value (not done by calculate) */ |
|
{ |
|
struct value val; |
|
Gcomplex (&val, a[num_params-1], 0.0); |
|
setvar (par_name[num_params-1], val); |
|
} |
|
|
/* call destructor for allocated vars */ |
/* call destructor for allocated vars */ |
lambda = -2; /* flag value, meaning 'destruct!' */ |
lambda = -2; /* flag value, meaning 'destruct!' */ |
(void) marquardt(a, C, &chisq, &lambda); |
(void) marquardt(a, C, &chisq, &lambda); |
Line 903 struct value data; |
|
Line 929 struct value data; |
|
|
|
if (!udv_ptr) { /* generate new entry */ |
if (!udv_ptr) { /* generate new entry */ |
last->next_udv = udv_ptr = (struct udvt_entry *) |
last->next_udv = udv_ptr = (struct udvt_entry *) |
gp_alloc((unsigned int) sizeof(struct udvt_entry), "fit setvar"); |
gp_alloc(sizeof(struct udvt_entry), "fit setvar"); |
udv_ptr->next_udv = NULL; |
udv_ptr->next_udv = NULL; |
} |
} |
safe_strncpy(udv_ptr->udv_name, varname, sizeof(udv_ptr->udv_name)); |
safe_strncpy(udv_ptr->udv_name, varname, sizeof(udv_ptr->udv_name)); |
|
|
if (columns == 1) |
if (columns == 1) |
int_error("Need 2 to 4 using specs", c_token); |
int_error("Need 2 to 4 using specs", c_token); |
|
|
|
/* The following patch was made by Remko Scharroo, 25-Mar-1999 |
|
* We need to check if one of the columns is time data, like |
|
* in plot2d and plot3d */ |
|
|
|
if (datatype[FIRST_X_AXIS] == TIME) { |
|
if (columns < 2) |
|
int_error("Need full using spec for x time data", c_token); |
|
df_timecol[0] = 1; |
|
} |
|
if (datatype[FIRST_Y_AXIS] == TIME) { |
|
if (columns < 1) |
|
int_error("Need using spec for y time data", c_token); |
|
df_timecol[1] = 1; |
|
} |
|
/* HBB 990326: added this check. Just in case some wants to fit |
|
* time/date data depending on two other variables ... */ |
|
if (datatype[FIRST_Z_AXIS] == TIME) { |
|
if (columns < 4) |
|
int_error("Need full using spec for z time data", c_token); |
|
else |
|
df_timecol[2] = 1; |
|
} |
|
/* End of patch by Remko Scharroo */ |
|
|
/* HBB 980401: if this is a single-variable fit, we shouldn't have |
/* HBB 980401: if this is a single-variable fit, we shouldn't have |
* allowed a variable name specifier for 'y': */ |
* allowed a variable name specifier for 'y': */ |
if ((dummy_y >= 0) && (columns < 4)) |
if ((dummy_y >= 0) && (columns < 4)) |
|
|
|
|
/* HBB 981210: two range specs mean different things, depending |
/* HBB 981210: two range specs mean different things, depending |
* on wether this is a 2D or 3D fit */ |
* on wether this is a 2D or 3D fit */ |
if (columns<4) { |
if (columns < 4) { |
if (zrange_token != -1) |
if (zrange_token != -1) |
int_error("Three range-specs not allowed in on-variable fit", zrange_token); |
int_error("Three range-specs not allowed in on-variable fit", zrange_token); |
else { |
else { |
/* 2D fit, 2 ranges: second range is for *z*, not y: */ |
/* 2D fit, 2 ranges: second range is for *z*, not y: */ |
autorange_z = autorange_y; |
autorange_z = autorange_y; |
min_z = min_y; |
/* HBB 20010926: Bug fixed. Checks were the wrong way round */ |
max_z = max_y; |
if (! (autorange_y & 1)) |
} |
min_z = min_y; |
|
if (! (autorange_y & 2)) |
|
max_z = max_y; |
|
} |
} |
} |
|
|
/* defer actually reading the data until we have parsed the rest |
/* defer actually reading the data until we have parsed the rest |
|
|
TBOOLEAN fixed; |
TBOOLEAN fixed; |
double tmp_par; |
double tmp_par; |
char c, *s; |
char c, *s; |
char sstr[MAX_LINE_LEN]; |
char sstr[MAX_LINE_LEN+1]; |
FILE *f; |
FILE *f; |
|
|
quote_str(sstr, c_token++, MAX_LINE_LEN); |
quote_str(sstr, c_token++, MAX_LINE_LEN); |