#include "fit.h"

#include "alloc.h"

#include "command.h"

#include "datablock.h"

#include "datafile.h"

#include "eval.h"

#include "gplocale.h"

#include "gp_time.h"

#include "matrix.h"

#include "misc.h"

#include "plot.h"

#include "setshow.h"

#include "scanner.h" /* For legal_identifier() */

#include "specfun.h"

#include "util.h"

#include <signal.h>

#if defined(VA_START) && defined(STDC_HEADERS)

static void Dblfn(const char *fmt, ...);

static void Dblfn();

#define Dblf Dblfn

#define Dblf2 Dblfn

#define Dblf3 Dblfn

#define Dblf5 Dblfn

#define Dblf6 Dblfn

#if defined(MSDOS) /* non-blocking IO stuff */

# include <io.h>

# include <conio.h>

# include <dos.h>

#ifdef _WIN32

# include <fcntl.h>

# include "win/winmain.h"

#ifdef INFINITY

# undef INFINITY

#define INFINITY 1e30

#define NEARLY_ZERO 1e-30

#define INITIAL_VALUE 1.0

#define DELTA 0.001

#define MAX_DATA 2048

#define MAX_PARAMS 32

#define MAX_LAMBDA 1e20

#define MIN_LAMBDA 1e-20

#define LAMBDA_UP_FACTOR 10

#define LAMBDA_DOWN_FACTOR 10

#if defined(MSDOS) || defined(OS2)

# define PLUSMINUS "\xF1" /* plusminus sign */

# define PLUSMINUS "+/-"

#define LASTFITCMDLENGTH 511

#define STANDARD stderr

#define BACKUP_SUFFIX ".old"

#define SQR(x) ((x) * (x))

enum marq_res {

OK, ML_ERROR, BETTER, WORSE

};

typedef enum marq_res marq_res_t;

JMP_BUF *fit_env = NULL;

char *fitlogfile = NULL;

TBOOLEAN fit_suppress_log = FALSE;

TBOOLEAN fit_errorvariables = TRUE;

TBOOLEAN fit_covarvariables = FALSE;

verbosity_level fit_verbosity = BRIEF;

TBOOLEAN fit_errorscaling = TRUE;

TBOOLEAN fit_prescale = TRUE;

char *fit_script = NULL;

int fit_wrap = 0;

TBOOLEAN fit_v4compatible = FALSE;

char *last_fit_command = NULL;

const char * FITLIMIT = "FIT_LIMIT";

const char * FITSTARTLAMBDA = "FIT_START_LAMBDA";

const char * FITLAMBDAFACTOR = "FIT_LAMBDA_FACTOR";

const char * FITMAXITER = "FIT_MAXITER";

static double epsilon = DEF_FIT_LIMIT; /* relative convergence limit */

double epsilon_abs = 0.0; /* default to zero non-relative limit */

int maxiter = 0;

static double startup_lambda = 0;

static double lambda_down_factor = LAMBDA_DOWN_FACTOR;

static double lambda_up_factor = LAMBDA_UP_FACTOR;

static const char fitlogfile_default[] = "fit.log";

static const char GNUFITLOG[] = "FIT_LOG";

static FILE *log_f = NULL;

static FILE *via_f = NULL;

static TBOOLEAN fit_show_lambda = TRUE;

static const char *GP_FIXED = "# FIXED";

static const char *FITSCRIPT = "FIT_SCRIPT";

static const char *DEFAULT_CMD = "replot"; /* if no fitscript spec. */

static int num_data;

static int num_params;

static int num_indep; /* # independent variables in fit function */

static int num_errors; /* # error columns */

static TBOOLEAN err_cols[MAX_NUM_VAR+1]; /* TRUE if variable has an associated error */

static int columns; /* # values read from data file for each point */

static double *fit_x = 0; /* all independent variable values,

static double *fit_z = 0; /* dependent data values */

static double *err_data = 0; /* standard deviations of indep. and dependent data */

static double *a = 0; /* array of fitting parameters */

static double **regress_C = 0; /* global copy of C matrix in regress */

static void (* regress_cleanup)(void) = NULL; /* memory cleanup function callback */

static TBOOLEAN user_stop = FALSE;

static double *scale_params = 0; /* scaling values for parameters */

static struct udft_entry func;

static fixstr *par_name;

static t_value **par_udv; /* array of pointers to the "via" variables */

static fixstr *last_par_name = NULL;

static int last_num_params = 0;

static char *last_dummy_var[MAX_NUM_VAR];

static udvt_entry *fit_dummy_udvs[MAX_NUM_VAR];

#if !defined(_WIN32) || defined(WGP_CONSOLE)

static RETSIGTYPE ctrlc_handle(int an_int);

static void ctrlc_setup(void);

static void fit_main(void);

static marq_res_t marquardt(double a[], double **alpha, double *chisq, double *lambda);

static void analyze(double a[], double **alpha, double beta[],

double *chisq, double **deriv);

static void calculate(double *zfunc, double **dzda, double a[]);

static void calc_derivatives(const double *par, double *data, double **deriv);

static TBOOLEAN fit_interrupt(void);

static TBOOLEAN regress(double a[]);

static void regress_init(void);

static void regress_finalize(int iter, double chisq, double last_chisq, double lambda, double **covar);

static void fit_show(int i, double chisq, double last_chisq, double *a,

double lambda, FILE * device);

static void fit_show_brief(int iter, double chisq, double last_chisq, double *parms,

double lambda, FILE * device);

static void show_results(double chisq, double last_chisq, double* a, double* dpar, double** corel);

static void log_axis_restriction(FILE *log_f, int param,

double min, double max, int autoscale, char *name);

static void print_function_definitions(struct at_type *at, FILE * device);

static TBOOLEAN is_empty(char *s);

static intgr_t getivar(const char *varname);

static double getdvar(const char *varname);

static double createdvar(char *varname, double value);

static void setvar(char *varname, double value);

static void setvarerr(char *varname, double value);

static void setvarcovar(char *varname1, char *varname2, double value);

static char *get_next_word(char **s, char *subst);

fit_command()

{

static JMP_BUF fit_jumppoint;

if (evaluate_inside_functionblock && inside_plot_command)

int_error(NO_CARET, "fit command not possible in this context");

inside_plot_command = TRUE;

/* Set up an exception handler for errors that occur during "fit".

fit_env = &fit_jumppoint;

if (SETJMP(*fit_env,1)) {

fit_env = NULL;

fprintf(stderr, "*** FIT ERROR ***\n");

free(last_fit_command);

last_fit_command = NULL;

while (!END_OF_COMMAND)

c_token++;

Ginteger( &(add_udv_by_name("FIT_ERROR")->udv_value), 1);

inside_plot_command = FALSE;

return;

}

fit_main();

fit_env = NULL;

Ginteger( &(add_udv_by_name("FIT_ERROR")->udv_value), 0);

inside_plot_command = FALSE;

}

#if !defined(_WIN32) || defined(WGP_CONSOLE)

static RETSIGTYPE

ctrlc_handle(int an_int)

{

(void) an_int; /* avoid -Wunused warning */

/* reinstall signal handler (necessary on SysV) */

(void) signal(SIGINT, (sigfunc) ctrlc_handle);

ctrlc_flag = TRUE;

}

static void

ctrlc_setup()

{

#if (defined(__EMX__) || !defined(MSDOS)) && (!defined(_WIN32) || defined(WGP_CONSOLE))

(void) signal(SIGINT, (sigfunc) ctrlc_handle);

}

#if defined(MSDOS)

int getchx(void);

getchx()

{

int c = getch();

if (!c || c == 0xE0) {

c <<= 8;

c |= getch();

}

return c;

}

error_ex(int t_num, const char *str, ...)

{

char buf[128];

va_list args;

va_start(args, str);

vsnprintf(buf, sizeof(buf), str, args);

va_end(args);

/* cleanup - free memory */

if (log_f) {

fprintf(log_f, "BREAK: %s", buf);

fclose(log_f);

log_f = NULL;

}

if (via_f) {

fclose(via_f);

via_f = NULL;

}

free(fit_x);

free(fit_z);

free(err_data);

free(a);

a = fit_x = fit_z = err_data = NULL;

if (func.at) {

free_at(func.at); /* release perm. action table */

func.at = (struct at_type *) NULL;

}

if (regress_cleanup != NULL)

(*regress_cleanup)();

/* the datafile may still be open */

df_close();

/* restore original SIGINT function */

interrupt_setup();

/* Usually int_error() returns to the top level command parser

int_error(t_num, buf);

}

static marq_res_t

marquardt(double a[], double **C, double *chisq, double *lambda)

{

int i, j;

static double *da = 0, /* delta-step of the parameter */

*temp_a = 0, /* temptative new params set */

*d = 0, *tmp_d = 0, **tmp_C = 0, *residues = 0, **deriv = 0;

double tmp_chisq;

/* Initialization when lambda == -1 */

if (*lambda == -1) { /* Get first chi-square check */

temp_a = vec(num_params);

d = vec(num_data + num_params);

tmp_d = vec(num_data + num_params);

da = vec(num_params);

residues = vec(num_data + num_params);

tmp_C = matr(num_data + num_params, num_params);

deriv = NULL;

if (num_errors > 1)

deriv = matr(num_errors - 1, num_data);

analyze(a, C, d, chisq, deriv);

/* Calculate a useful startup value for lambda, as given by Schwarz */

if (startup_lambda != 0)

*lambda = startup_lambda;

else {

*lambda = 0;

for (i = 0; i < num_data; i++)

for (j = 0; j < num_params; j++)

*lambda += C[i][j] * C[i][j];

*lambda = sqrt(*lambda / num_data / num_params);

}

/* Fill in the lower square part of C (the diagonal is filled in on

for (i = 0; i < num_params; i++)

for (j = 0; j < i; j++)

C[num_data + i][j] = 0, C[num_data + j][i] = 0;

return OK;

}

/* once converged, free allocated vars */

if (*lambda == -2) {

free(d);

free(tmp_d);

free(da);

free(temp_a);

free(residues);

free_matr(tmp_C);

free_matr(deriv);

/* may be called more than once */

d = tmp_d = da = temp_a = residues = (double *) NULL;

tmp_C = deriv = (double **) NULL;

return OK;

}

/* Givens calculates in-place, so make working copies of C and d */

for (j = 0; j < num_data + num_params; j++)

memcpy(tmp_C[j], C[j], num_params * sizeof(double));

memcpy(tmp_d, d, num_data * sizeof(double));

/* fill in additional parts of tmp_C, tmp_d */

for (i = 0; i < num_params; i++) {

/* fill in low diag. of tmp_C ... */

tmp_C[num_data + i][i] = *lambda;

/* ... and low part of tmp_d */

tmp_d[num_data + i] = 0;

}

Givens(tmp_C, tmp_d, da, num_params + num_data, num_params);

/* check if trial did ameliorate sum of squares */

for (j = 0; j < num_params; j++)

temp_a[j] = a[j] + da[j];

analyze(temp_a, tmp_C, tmp_d, &tmp_chisq, deriv);

/* tsm patchset 230: Changed < to <= in next line */

/* so finding exact minimum stops iteration instead of just increasing lambda. */

/* Disadvantage is that if lambda is large enough so that chisq doesn't change */

/* is taken as success. */

if (tmp_chisq <= *chisq) { /* Success, accept new solution */

if (*lambda > MIN_LAMBDA) {

if (fit_verbosity == VERBOSE)

putc('/', stderr);

*lambda /= lambda_down_factor;

}

/* update chisq, C, d, a */

*chisq = tmp_chisq;

for (j = 0; j < num_data; j++) {

memcpy(C[j], tmp_C[j], num_params * sizeof(double));

d[j] = tmp_d[j];

}

for (j = 0; j < num_params; j++)

a[j] = temp_a[j];

return BETTER;

} else { /* failure, increase lambda and return */

*lambda *= lambda_up_factor;

if (fit_verbosity == VERBOSE)

(void) putc('*', stderr);

else if (fit_verbosity == BRIEF) /* one-line report even if chisq increases */

fit_show_brief(-1, tmp_chisq, *chisq, temp_a, *lambda, STANDARD);

return WORSE;

}

}

static double

effective_error(double **deriv, int i)

{

double tot_err;

int j, k;

if (num_errors <= 1) /* z-errors or equal weights */

tot_err = err_data[i];

else {

/* "Effective variance" according to

tot_err = SQR(err_data[i * num_errors + (num_errors - 1)]);

for (j = 0, k = 0; j < num_indep; j++) {

if (err_cols[j]) {

tot_err += SQR(deriv[k][i] * err_data[i * num_errors + k]);

k++;

}

}

tot_err = sqrt(tot_err);

}

return tot_err;

}

static void

analyze(double a[], double **C, double d[], double *chisq, double ** deriv)

{

int i, j;

calculate(d, C, a);

/* derivatives in indep. variables are required for

if (num_errors > 1)

calc_derivatives(a, d, deriv);

for (i = 0; i < num_data; i++) {

double err = effective_error(deriv, i);

/* note: order reversed, as used by Schwarz */

d[i] = (d[i] - fit_z[i]) / err;

for (j = 0; j < num_params; j++)

C[i][j] /= err;

}

*chisq = sumsq_vec(num_data, d);

}

static void

calculate(double *zfunc, double **dzda, double a[])

{

int k, p;

double tmp_a;

double *tmp_high, *tmp_pars;

#ifdef TWO_SIDE_DIFFERENTIATION

double *tmp_low;

tmp_high = vec(num_data); /* numeric derivations */

#ifdef TWO_SIDE_DIFFERENTIATION

tmp_low = vec(num_data);

tmp_pars = vec(num_params);

/* first function values */

call_gnuplot(a, zfunc);

/* then derivatives in parameters */

for (p = 0; p < num_params; p++)

tmp_pars[p] = a[p];

for (p = 0; p < num_params; p++) {

tmp_a = fabs(a[p]) < NEARLY_ZERO ? NEARLY_ZERO : a[p];

tmp_pars[p] = tmp_a * (1 + DELTA);

call_gnuplot(tmp_pars, tmp_high);

#ifdef TWO_SIDE_DIFFERENTIATION

tmp_pars[p] = tmp_a * (1 - DELTA);

call_gnuplot(tmp_pars, tmp_low);

for (k = 0; k < num_data; k++)

#ifdef TWO_SIDE_DIFFERENTIATION

dzda[k][p] = (tmp_high[k] - tmp_low[k]) / (2 * tmp_a * DELTA);

dzda[k][p] = (tmp_high[k] - zfunc[k]) / (tmp_a * DELTA);

tmp_pars[p] = a[p];

}

#ifdef TWO_SIDE_DIFFERENTIATION

free(tmp_low);

free(tmp_high);

free(tmp_pars);

}

call_gnuplot(const double *par, double *data)

{

int i, j;

struct value v;

/* set parameters first */

for (i = 0; i < num_params; i++)

Gcomplex(par_udv[i], par[i] * scale_params[i], 0.0);

for (i = 0; i < num_data; i++) {

/* calculate fit-function value */

/* initialize extra dummy variables from the corresponding

for (j = 0; j < MAX_NUM_VAR; j++) {

double dummy_value;

struct udvt_entry *udv = fit_dummy_udvs[j];

if (!udv)

int_error(NO_CARET, "Internal error: lost a dummy parameter!");

if (udv->udv_value.type == CMPLX || udv->udv_value.type == INTGR)

dummy_value = real(&(udv->udv_value));

else

dummy_value = 0.0;

Gcomplex(&func.dummy_values[j], dummy_value, 0.0);

}

/* set actual dummy variables from file data */

for (j = 0; j < num_indep; j++)

Gcomplex(&func.dummy_values[j],

fit_x[i * num_indep + j], 0.0);

evaluate_at(func.at, &v);

if (undefined || isnan(real(&v))) {

/* Print useful info on undefined-function error. */

Dblf("\nCurrent data point\n");

Dblf("=========================\n");

Dblf3("%-15s = %i out of %i\n", "#", i + 1, num_data);

for (j = 0; j < num_indep; j++)

Dblf3("%-15.15s = %-15g\n", c_dummy_var[j], par[j] * scale_params[j]);

Dblf3("%-15.15s = %-15g\n", "z", fit_z[i]);

Dblf("\nCurrent set of parameters\n");

Dblf("=========================\n");

for (j = 0; j < num_params; j++)

Dblf3("%-15.15s = %-15g\n", par_name[j], par[j] * scale_params[j]);

Dblf("\n");

if (undefined) {

Eex("Undefined value during function evaluation");

} else {

Eex("Function evaluation yields NaN (\"not a number\")");

}

}

data[i] = real(&v);

}

}

static void

calc_derivatives(const double *par, double *data, double **deriv)

{

int i, j, k, m;

struct value v;

double h;

/* set parameters first */

for (i = 0; i < num_params; i++)

Gcomplex(par_udv[i], par[i] * scale_params[i], 0.0);

for (i = 0; i < num_data; i++) { /* loop over data points */

for (j = 0, m = 0; j < num_indep; j++) { /* loop over indep. variables */

double tmp_high;

double tmp_x;

#ifdef TWO_SIDE_DIFFERENTIATION

double tmp_low;

/* only calculate derivatives if necessary */

if (!err_cols[j])

continue;

/* set actual dummy variables from file data */

for (k = 0; k < num_indep; k++) {

if (j != k)

Gcomplex(&func.dummy_values[k],

fit_x[i * num_indep + k], 0.0);

}

tmp_x = fit_x[i * num_indep + j];

/* optimal step size */

h = GPMAX(DELTA * fabs(tmp_x), 8*1e-8*(fabs(tmp_x) + 1e-8));

Gcomplex(&func.dummy_values[j], tmp_x + h, 0.0);

evaluate_at(func.at, &v);

tmp_high = real(&v);

#ifdef TWO_SIDE_DIFFERENTIATION

Gcomplex(&func.dummy_values[j], tmp_x - h, 0.0);

evaluate_at(func.at, &v);

tmp_low = real(&v);

deriv[m][i] = (tmp_high - tmp_low) / (2 * h);

deriv[m][i] = (tmp_high - data[i]) / h;

m++;

}

}

}

static TBOOLEAN

fit_interrupt()

{

while (TRUE) {

fputs("\n\n(S)top fit, (C)ontinue, (E)xecute FIT_SCRIPT: ", STANDARD);

#ifdef _WIN32

WinRaiseConsole();

switch (getchar()) {

case EOF:

case 's':

case 'S':

fputs("Stop.\n", STANDARD);

user_stop = TRUE;

return FALSE;

case 'c':

case 'C':

fputs("Continue.\n", STANDARD);

return TRUE;

case 'e':

case 'E':{

int i;

const char *tmp = getfitscript();

fprintf(STANDARD, "executing: %s\n", tmp);

/* FIXME: Shouldn't we also set FIT_STDFIT etc? */

/* set parameters visible to gnuplot */

for (i = 0; i < num_params; i++)

Gcomplex(par_udv[i], a[i] * scale_params[i], 0.0);

do_string(tmp);

}

}

}

return TRUE;

}

const char *

getfitscript(void)

{

char *tmp;

if (fit_script != NULL)

return fit_script;

if ((tmp = getenv(FITSCRIPT)) != NULL)

return tmp;

else

return DEFAULT_CMD;

}

static void

regress_init(void)

{

struct udvt_entry *v; /* For exporting results to the user */

/* Reset flag describing fit result status */

v = add_udv_by_name("FIT_CONVERGED");

Ginteger(&v->udv_value, 0);

/* Ctrl-C now serves as Hotkey */

ctrlc_setup();

/* HBB 981118: initialize new variable 'user_break' */

user_stop = FALSE;

}

static void

regress_finalize(int iter, double chisq, double last_chisq, double lambda, double **covar)

{

int i, j;

struct udvt_entry *v; /* For exporting results to the user */

int ndf;

int niter;

double stdfit;

double pvalue;

double *dpar;

double **corel = NULL;

TBOOLEAN covar_invalid = FALSE;

/* restore original SIGINT function */

interrupt_setup();

/* tsm patchset 230: final progress report labels to console */

if (fit_verbosity == BRIEF)

fit_show_brief(-2, chisq, chisq, a, lambda, STANDARD);

/* tsm patchset 230: final progress report to log file */

if (!fit_suppress_log) {

if (fit_verbosity == VERBOSE)

fit_show(iter, chisq, last_chisq, a, lambda, log_f);

else

fit_show_brief(iter, chisq, last_chisq, a, lambda, log_f);

}

/* test covariance matrix */

if (covar != NULL) {

for (i = 0; i < num_params; i++) {

/* diagonal elements must be larger than zero */

if (covar[i][i] <= 0.0) {

/* Not a fatal error, but prevent floating point exception later on */

Dblf2("Calculation error: non-positive diagonal element in covar. matrix of parameter '%s'.\n", par_name[i]);

covar_invalid = TRUE;

}

}

}

/* HBB 970304: the maxiter patch: */

if ((maxiter > 0) && (iter > maxiter)) {

Dblf2("\nMaximum iteration count (%d) reached. Fit stopped.\n", maxiter);

} else if (user_stop) {

Dblf2("\nThe fit was stopped by the user after %d iterations.\n", iter);

} else if (lambda >= MAX_LAMBDA) {

Dblf2("\nThe maximum lambda = %e was exceeded. Fit stopped.\n", MAX_LAMBDA);

} else if (covar_invalid) {

Dblf2("\nThe covariance matrix is invalid. Fit did not converge properly.\n");

} else {

Dblf2("\nAfter %d iterations the fit converged.\n", iter);

v = add_udv_by_name("FIT_CONVERGED");

Ginteger(&v->udv_value, 1);

}

/* fit results */

ndf = num_data - num_params;

stdfit = sqrt(chisq / ndf);

pvalue = 1. - chisq_cdf(ndf, chisq);

niter = iter;

/* Export these to user-accessible variables */

v = add_udv_by_name("FIT_NDF");

Ginteger(&v->udv_value, ndf);

v = add_udv_by_name("FIT_STDFIT");

Gcomplex(&v->udv_value, stdfit, 0);

v = add_udv_by_name("FIT_WSSR");

Gcomplex(&v->udv_value, chisq, 0);

v = add_udv_by_name("FIT_P");

Gcomplex(&v->udv_value, pvalue, 0);

v = add_udv_by_name("FIT_NITER");

Ginteger(&v->udv_value, niter);

/* Save final parameters. Depending on the backend and

for (i = 0; i < num_params; i++)

Gcomplex(par_udv[i], a[i] * scale_params[i], 0.0);

/* Set error and covariance variables to zero,

if (fit_errorvariables) {

for (i = 0; i < num_params; i++)

setvarerr(par_name[i], 0.0);

}

if (fit_covarvariables) {

/* first, remove all previous covariance variables */

del_udv_by_name("FIT_COV_", TRUE);

for (i = 0; i < num_params; i++) {

for (j = 0; j < i; j++) {

setvarcovar(par_name[i], par_name[j], 0.0);

setvarcovar(par_name[j], par_name[i], 0.0);

}

setvarcovar(par_name[i], par_name[i], 0.0);

}

}

/* calculate unscaled parameter errors in dpar[]: */

dpar = vec(num_params);

if ((covar != NULL) && !covar_invalid) {

/* calculate unscaled parameter errors in dpar[]: */

for (i = 0; i < num_params; i++) {

dpar[i] = sqrt(covar[i][i]);

}

/* transform covariances into correlations */

corel = matr(num_params, num_params);

for (i = 0; i < num_params; i++) {

/* only lower triangle needs to be handled */

for (j = 0; j < i; j++)

corel[i][j] = covar[i][j] / (dpar[i] * dpar[j]);

corel[i][i] = 1.;

}

} else {

/* set all errors to zero if covariance matrix invalid or unavailable */

for (i = 0; i < num_params; i++)

dpar[i] = 0.0;

}

if (fit_errorscaling || (num_errors == 0)) {

/* scale parameter errors based on chisq */

double temp = sqrt(chisq / (num_data - num_params));

for (i = 0; i < num_params; i++)

dpar[i] *= temp;

}

/* Save user error variables. */

if (fit_errorvariables) {

for (i = 0; i < num_params; i++)

setvarerr(par_name[i], dpar[i] * scale_params[i]);

}

/* fill covariance variables if needed */

if (fit_covarvariables && (covar != NULL) && !covar_invalid) {

double scale =

(fit_errorscaling || (num_errors == 0)) ?

(chisq / (num_data - num_params)) : 1.0;

for (i = 0; i < num_params; i++) {

/* only lower triangle needs to be handled */

for (j = 0; j <= i; j++) {

double temp = scale * scale_params[i] * scale_params[j];

setvarcovar(par_name[i], par_name[j], covar[i][j] * temp);

setvarcovar(par_name[j], par_name[i], covar[i][j] * temp);

}

}