Back to index

plt-scheme  4.2.1
Defines | Enumerations | Functions | Variables
fit.c File Reference
#include "matrix.h"
#include "fit.h"
#include <math.h>

Go to the source code of this file.

Defines

#define FIT_MAIN
#define NULL   0
#define INFINITY   1e30
#define NEARLY_ZERO   1e-30
#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
#define TBOOLEAN   int
#define PLUSMINUS   "+/-"
#define TWO_SIDE_DIFFERENTIATION

Enumerations

enum  marq_res { OK, ERROR, BETTER, WORSE }

Functions

MZ_DLLEXPORT double get_rms ()
MZ_DLLEXPORT double get_varience ()
MZ_DLLEXPORT double * get_asym_error ()
MZ_DLLEXPORT double * get_asym_error_percent ()
static void print_matrix_and_vectors (double **C, double *d, double *r, int m, int n)
static marq_res_t marquardt (double a[], double **alpha, double *chisq, double *lambda)
static TBOOLEAN analyze (double a[], double **alpha, double beta[], double *chisq)
static void calculate (double *zfunc, double **dzda, double a[])
static void call_scheme (double *par, double *data)
static TBOOLEAN regress (double a[])
static void printmatrix (double **C, int m, int n)
static marq_res_t marquardt (a, double **C, double *chisq, double *lambda)
static TBOOLEAN analyze (a, double **C, d, double *chisq)
static void calculate (double *zfunc, double **dzda, a)
static TBOOLEAN regress (a)
MZ_DLLEXPORT double * do_fit (void *function, int n_values, double *x_values, double *y_values, double *z_values, double *errors, int n_parameters, double *parameters)

Variables

static int max_data
static int max_params
static double epsilon = 1e-5
static int maxiter = 0
static char * FIXED = "# FIXED"
static char * GNUFITLOG = "FIT_LOG"
static char * FITLIMIT = "FIT_LIMIT"
static char * FITSTARTLAMBDA = "FIT_START_LAMBDA"
static char * FITLAMBDAFACTOR = "FIT_LAMBDA_FACTOR"
static char * FITMAXITER = "FIT_MAXITER"
static char * FITSCRIPT = "FIT_SCRIPT"
static char * DEFAULT_CMD = "replot"
static int num_data
static int num_params
static int columns
static double * fit_x
static double * fit_y
static double * fit_z
static double * err_data
static double * a
static double startup_lambda = 0
static double lambda_down_factor = LAMBDA_DOWN_FACTOR
static double lambda_up_factor = LAMBDA_UP_FACTOR
static voidcurrent_fun
double rms = 0
double varience = 0
double * asym_error
double * asym_error_percent

Define Documentation

#define DELTA   0.001

Definition at line 87 of file fit.c.

#define FIT_MAIN

Definition at line 58 of file fit.c.

#define INFINITY   1e30

Definition at line 82 of file fit.c.

#define LAMBDA_DOWN_FACTOR   10

Definition at line 94 of file fit.c.

#define LAMBDA_UP_FACTOR   10

Definition at line 93 of file fit.c.

#define MAX_DATA   2048

Definition at line 89 of file fit.c.

#define MAX_LAMBDA   1e20

Definition at line 91 of file fit.c.

#define MAX_PARAMS   32

Definition at line 90 of file fit.c.

#define MIN_LAMBDA   1e-20

Definition at line 92 of file fit.c.

#define NEARLY_ZERO   1e-30

Definition at line 83 of file fit.c.

#define NULL   0

Definition at line 60 of file fit.c.

#define PLUSMINUS   "+/-"

Definition at line 98 of file fit.c.

#define TBOOLEAN   int

Definition at line 96 of file fit.c.

Definition at line 364 of file fit.c.


Enumeration Type Documentation

enum marq_res
Enumerator:
OK 
ERROR 
BETTER 
WORSE 

Definition at line 73 of file fit.c.

              {
    OK, ERROR, BETTER, WORSE
};

Function Documentation

static TBOOLEAN analyze ( double  a[],
double **  alpha,
double  beta[],
double *  chisq 
) [static]

Here is the caller graph for this function:

static TBOOLEAN analyze ( a  ,
double **  C,
,
double *  chisq 
) [static]

Definition at line 332 of file fit.c.

{
/*
 *  used by marquardt to evaluate the linearized fitting matrix C
 *  and vector d, fills in only the top part of C and d
 *  I don't use a temporary array zfunc[] any more. Just use
 *  d[] instead.
 */
    int i, j;

    *chisq = 0;
    calculate(d, C, a);

    for (i = 0; i < num_data; i++) {
       /* note: order reversed, as used by Schwarz */
       d[i] = (d[i] - fit_z[i]) / err_data[i];
       *chisq += d[i] * d[i];
       for (j = 0; j < num_params; j++)
           C[i][j] /= err_data[i];
    }
    /* FIXME: why return a value that is always TRUE ? */
    return 1;
}

Here is the call graph for this function:

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

Here is the caller graph for this function:

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

Definition at line 369 of file fit.c.

{
    int k, p;
    double tmp_a;
    double *tmp_high, *tmp_pars;
#ifdef TWO_SIDE_DIFFERENTIATION
    double *tmp_low;
#endif

    tmp_high = vec(num_data);      /* numeric derivations */
#ifdef TWO_SIDE_DIFFERENTIATION
    tmp_low = vec(num_data);
#endif
    tmp_pars = vec(num_params);

    /* first function values */

    call_scheme(a, zfunc);

    /* then derivatives */

    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_scheme(tmp_pars, tmp_high);
#ifdef TWO_SIDE_DIFFERENTIATION
       tmp_pars[p] = tmp_a * (1 - DELTA);
       call_scheme(tmp_pars, tmp_low);
#endif
       for (k = 0; k < num_data; k++)
#ifdef TWO_SIDE_DIFFERENTIATION
           dzda[k][p] = (tmp_high[k] - tmp_low[k]) / (2 * tmp_a * DELTA);
#else
           dzda[k][p] = (tmp_high[k] - zfunc[k]) / (tmp_a * DELTA);
#endif
       tmp_pars[p] = a[p];
    }

}

Here is the call graph for this function:

static void call_scheme ( double *  par,
double *  data 
) [static]

Definition at line 418 of file fit.c.

{
  int rators = 2 + num_params;
  double * rands =
    (double *) malloc(rators * sizeof(double));

  int i;

  /* set up the constant params */
  for(i = 0 ; i< num_params; i++) {
    rands[i+2] = par[i];
  }

  /* now calculate the function at the existing points */
  for (i = 0; i < num_data; i++) {
    rands[0] = fit_x[i];
    rands[1] = fit_y[i];

    data[i] = ((double (*) (int, double *) )current_fun) // ouch!
               (rators, rands);
  }

  free(rands);

}

Here is the caller graph for this function:

MZ_DLLEXPORT double* do_fit ( void function,
int  n_values,
double *  x_values,
double *  y_values,
double *  z_values,
double *  errors,
int  n_parameters,
double *  parameters 
)

Definition at line 710 of file fit.c.

                                   {

  /* reset lambda and other parameters if desired */
  int i;
  current_fun = function;

  num_data = n_values;
  fit_x = x_values;
  fit_y = y_values;
  fit_z = z_values; /* value is stored in z */
  err_data = errors;

  a = parameters;
  num_params = n_parameters;

  /* redim_vec(&a, num_params); */
  /* par_name = (fixstr *) gp_realloc(par_name, (num_params + 1) * sizeof(fixstr), "fit param"); */

  /* avoid parameters being equal to zero */
  for (i = 0; i < num_params; i++) {
    if (a[i] == 0) {
      a[i] = NEARLY_ZERO; 
    }
  }
  
  if(regress(a)) {
    gc_cleanup();
    return a;
  }
  else { /* something went wrong */
    gc_cleanup(); 
    return NULL;
  }
}

Here is the call graph for this function:

Definition at line 156 of file fit.c.

{return asym_error;}

Definition at line 160 of file fit.c.

MZ_DLLEXPORT double get_rms ( )

Definition at line 148 of file fit.c.

{return rms;}

Definition at line 152 of file fit.c.

{return varience;}
static marq_res_t marquardt ( double  a[],
double **  alpha,
double *  chisq,
double *  lambda 
) [static]

Here is the caller graph for this function:

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

Definition at line 221 of file fit.c.

{
    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;
    double tmp_chisq;

    /* Initialization when lambda == -1 */

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

       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);

       analyze_ret = analyze(a, C, d, chisq);

       /* Calculate a useful startup value for lambda, as given by Schwarz */
       /* FIXME: this is doesn't turn out to be much better, really... */
       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
          each iteration, see below) */
       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;
       /* printmatrix(C, num_data+num_params, num_params); */
       return analyze_ret ? OK : ERROR;
    }
    /* once converged, free dynamic allocated vars */

    if (*lambda == -2) {
       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;
    }
    /* printmatrix(tmp_C, num_data+num_params, num_params); */

    /* FIXME: residues[] isn't used at all. Why? Should it be used? */

    Givens(tmp_C, tmp_d, da, residues, num_params + num_data, num_params, 1);
    /*print_matrix_and_vectors (tmp_C, tmp_d, residues,
       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];

    if (!analyze(temp_a, tmp_C, tmp_d, &tmp_chisq)) {
       /* FIXME: will never be reached: always returns TRUE */
       return ERROR;
    }

    if (tmp_chisq < *chisq) {      /* Success, accept new solution */
       if (*lambda > MIN_LAMBDA) {
         /* (void) putc('/', stderr); */
          *lambda /= lambda_down_factor;
       }
       *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 */
      /* (void) putc('*', stderr); */
      *lambda *= lambda_up_factor;
      return WORSE;
    }
}

Here is the call graph for this function:

static void print_matrix_and_vectors ( double **  C,
double *  d,
double *  r,
int  m,
int  n 
) [static]

Definition at line 202 of file fit.c.

{
    int i, j;

    for (i = 0; i < m; i++) {
       for (j = 0; j < n; j++);
         /* Dblf2("%8g ", C[i][j]); */
         /* Dblf3("| %8g | %8g\n", d[i], r[i]); */
    }
    /* Dblf("\n"); */
}
static void printmatrix ( double **  C,
int  m,
int  n 
) [static]

Definition at line 185 of file fit.c.

{
    int i, j;

    for (i = 0; i < m; i++) {
       for (j = 0; j < n - 1; j++);
         /* Dblf2("%.8g |", C[i][j]); */
       /* Dblf2("%.8g\n", C[i][j]); */
    }
    /* Dblf("\n"); */
}
static TBOOLEAN regress ( double  a[]) [static]

Here is the caller graph for this function:

static TBOOLEAN regress ( a  ) [static]

Definition at line 476 of file fit.c.

{
  double **covar, *dpar, **C, chisq, last_chisq, lambda;
  int iter, i, j;
  marq_res_t res;
  
  chisq = last_chisq = INFINITY;
  C = matr(num_data + num_params, num_params);
  lambda = -1;          /* use sign as flag */
  iter = 0;             /* iteration counter  */

  /* Initialize internal variables and 1st chi-square check */

  if ((res = marquardt(a, C, &chisq, &lambda)) == ERROR)
    return 0; /* an error occurded */

  res = BETTER;

  /* show_fit(iter, chisq, chisq, a, lambda, STANDARD); */

  /* MAIN FIT LOOP: do the regression iteration */

  do {
    if (res == BETTER) {
      iter++;
      last_chisq = chisq;
    }
    if ((res = marquardt(a, C, &chisq, &lambda)) == BETTER)
      {};
    /* show_fit(iter, chisq, last_chisq, a, lambda, STANDARD); */
  } while ((res != ERROR)
          && (lambda < MAX_LAMBDA)
          && ((maxiter == 0) || (iter <= maxiter))
          && (res == WORSE
              || ((chisq > NEARLY_ZERO)
                 ? ((last_chisq - chisq) / chisq)
                 : (last_chisq - chisq)) > epsilon
              )
          );

    /* fit done */

  /* save all the info that was otherwise printed out */

  rms = sqrt(chisq / (num_data - num_params));
  varience = chisq / (num_data - num_params);
  asym_error = malloc (num_params * sizeof (double));  
  asym_error_percent = malloc (num_params * sizeof (double)) ;

  /* don't know what the following code does... */

  /* compute covar[][] directly from C */ 
  Givens(C, 0, 0, 0, num_data, num_params, 0); 
  covar = C + num_data; 
  Invert_RtR(C, covar, num_params);
  
  dpar = vec(num_params); 
  for (i = 0; i < num_params; i++) {
    /* FIXME: can this still happen ? */
      if (covar[i][i] <= 0.0)      /* HBB: prevent floating point exception later on */
       return 0;  /* Eex("Calculation error: non-positive diagonal element in covar. matrix"); */
    dpar[i] = sqrt(covar[i][i]);
  }  

  /* transform covariances into correlations */
  for (i = 0; i < num_params; i++) {
    /* only lower triangle needs to be handled */
      for (j = 0; j <= i; j++)
      covar[i][j] /= dpar[i] * dpar[j];
  }

  /* scale parameter errors based on chisq */
  chisq = sqrt(chisq / (num_data - num_params));
  for (i = 0; i < num_params; i++)
    dpar[i] *= chisq;

  for(i = 0; i< num_params; i++)
    {
      double temp = 
       (fabs(a[i]) < NEARLY_ZERO) ? 0.0 : fabs(100.0 * dpar[i] / a[i]);
      asym_error[i] = dpar[i];
      asym_error_percent[i] = temp;
    }

  return 1;


    /******** CRAP LEFT OVER FROM GNUPLOT ***********/

    /* HBB 970304: the maxiter patch: */
    /*
    if ((maxiter > 0) && (iter > maxiter)) {
       Dblf2("\nMaximum iteration count (%d) reached. Fit stopped.\n", maxiter);
    } else  {
       Dblf2("\nAfter %d iterations the fit converged.\n", iter);
    }

    Dblf2("final sum of squares of residuals : %g\n", chisq);
    if (chisq > NEARLY_ZERO) {
       Dblf2("rel. change during last iteration : %g\n\n", (chisq - last_chisq) / chisq);
    } else {
       Dblf2("abs. change during last iteration : %g\n\n", (chisq - last_chisq));
    }

    if (res == ERROR)
      // Eex("FIT: error occurred during fit");
    */
    /* compute errors in the parameters */

  /*   if (num_data == num_params) { */
/*     int i; */

/*     Dblf("\nExactly as many data points as there are parameters.\n"); */
/*     Dblf("In this degenerate case, all errors are zero by definition.\n\n"); */
/*     Dblf("Final set of parameters \n"); */
/*     Dblf("======================= \n\n"); */
/*     for (i = 0; i < num_params; i++) */
/*         Dblf3("%-15.15s = %-15g\n", par_name[i], a[i]); */
/*     } else if (chisq < NEARLY_ZERO) { */
/*     int i; */

/*     Dblf("\nHmmmm.... Sum of squared residuals is zero. Can't compute errors.\n\n"); */
/*     Dblf("Final set of parameters \n"); */
/*     Dblf("======================= \n\n"); */
/*     for (i = 0; i < num_params; i++) */
/*         Dblf3("%-15.15s = %-15g\n", par_name[i], a[i]); */
/*     } else { */
/*     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("variance of residuals (reduced chisquare) = WSSR/ndf : %g\n\n", chisq / (num_data - num_params)); */
 
/*     /\* get covariance-, Korrelations- and Kurvature-Matrix *\/ */
/*     /\* and errors in the parameters                     *\/ */

/*     /\* compute covar[][] directly from C *\/ */
/*     Givens(C, 0, 0, 0, num_data, num_params, 0); */
/*     /\*printmatrix(C, num_params, num_params); *\/ */

/*     /\* Use lower square of C for covar *\/ */
/*     covar = C + num_data; */
/*     Invert_RtR(C, covar, num_params); */
/*     /\*printmatrix(covar, num_params, num_params); *\/ */

/*     /\* calculate unscaled parameter errors in dpar[]: *\/ */
/*     dpar = vec(num_params); */
/*     for (i = 0; i < num_params; i++) { */
/*         /\* FIXME: can this still happen ? *\/ */
/*         if (covar[i][i] <= 0.0) /\* HBB: prevent floating point exception later on *\/ */
/*            Eex("Calculation error: non-positive diagonal element in covar. matrix"); */
/*         dpar[i] = sqrt(covar[i][i]); */
/*     } */

/*     /\* transform covariances into correlations *\/ */
/*     for (i = 0; i < num_params; i++) { */
/*         /\* only lower triangle needs to be handled *\/ */
/*         for (j = 0; j <= i; j++) */
/*            covar[i][j] /= dpar[i] * dpar[j]; */
/*     } */

/*     /\* scale parameter errors based on chisq *\/ */
/*     chisq = sqrt(chisq / (num_data - num_params)); */
/*     for (i = 0; i < num_params; i++) */
/*         dpar[i] *= chisq; */

/*     Dblf("Final set of parameters            Asymptotic Standard Error\n"); */
/*     Dblf("=======================            ==========================\n\n"); */

/*     for (i = 0; i < num_params; i++) { */
/*         double temp = */
/*         (fabs(a[i]) < NEARLY_ZERO) ? 0.0 : fabs(100.0 * dpar[i] / a[i]); */
/*         Dblf6("%-15.15s = %-15g  %-3.3s %-12.4g (%.4g%%)\n", */
/*              par_name[i], a[i], PLUSMINUS, dpar[i], temp); */
/*     } */

/*     Dblf("\n\ncorrelation matrix of the fit parameters:\n\n"); */
/*     Dblf("               "); */

/*     for (j = 0; j < num_params; j++) */
/*         Dblf2("%-6.6s ", par_name[j]); */

/*     Dblf("\n"); */
/*     for (i = 0; i < num_params; i++) { */
/*         Dblf2("%-15.15s", par_name[i]); */
/*         for (j = 0; j <= i; j++) { */
/*            /\* Only print lower triangle of symmetric matrix *\/ */
/*            Dblf2("%6.3f ", covar[i][j]); */
/*         } */
/*         Dblf("\n"); */
/*     } */

/*     free(dpar); */
/*     } */

    return 1;
}

Here is the call graph for this function:


Variable Documentation

double* a [static]

Definition at line 125 of file fit.c.

double* asym_error

Definition at line 144 of file fit.c.

Definition at line 145 of file fit.c.

int columns [static]

Definition at line 120 of file fit.c.

void* current_fun [static]

Definition at line 135 of file fit.c.

char* DEFAULT_CMD = "replot" [static]

Definition at line 116 of file fit.c.

double epsilon = 1e-5 [static]

Definition at line 106 of file fit.c.

double* err_data [static]

Definition at line 124 of file fit.c.

double* fit_x [static]

Definition at line 121 of file fit.c.

double* fit_y [static]

Definition at line 122 of file fit.c.

double* fit_z [static]

Definition at line 123 of file fit.c.

char* FITLAMBDAFACTOR = "FIT_LAMBDA_FACTOR" [static]

Definition at line 113 of file fit.c.

char* FITLIMIT = "FIT_LIMIT" [static]

Definition at line 111 of file fit.c.

char* FITMAXITER = "FIT_MAXITER" [static]

Definition at line 114 of file fit.c.

char* FITSCRIPT = "FIT_SCRIPT" [static]

Definition at line 115 of file fit.c.

char* FITSTARTLAMBDA = "FIT_START_LAMBDA" [static]

Definition at line 112 of file fit.c.

char* FIXED = "# FIXED" [static]

Definition at line 109 of file fit.c.

char* GNUFITLOG = "FIT_LOG" [static]

Definition at line 110 of file fit.c.

Definition at line 131 of file fit.c.

Definition at line 132 of file fit.c.

int max_data [static]

Definition at line 103 of file fit.c.

int max_params [static]

Definition at line 104 of file fit.c.

int maxiter = 0 [static]

Definition at line 107 of file fit.c.

int num_data [static]

Definition at line 119 of file fit.c.

int num_params [static]

Definition at line 119 of file fit.c.

double rms = 0

Definition at line 142 of file fit.c.

double startup_lambda = 0 [static]

Definition at line 130 of file fit.c.

double varience = 0

Definition at line 143 of file fit.c.