Performs a Gaussian fit on input data. More...

#include "mdb.h"
#include "SDDS.h"
#include "scan.h"

Macros
#define	GUESS_BASELINE_GIVEN 0x0001

#define	FIX_BASELINE_GIVEN (0x0001 << 4)

#define	GUESS_HEIGHT_GIVEN 0x0002

#define	FIX_HEIGHT_GIVEN (0x0002 << 4)

#define	GUESS_MEAN_GIVEN 0x0004

#define	FIX_MEAN_GIVEN (0x0004 << 4)

#define	GUESS_SIGMA_GIVEN 0x0008

#define	FIX_SIGMA_GIVEN (0x0008 << 4)

#define	BASELINE_INDEX 0

#define	HEIGHT_INDEX 1

#define	MEAN_INDEX 2

#define	SIGMA_INDEX 3

Enumerations
enum	option_type { SET_FITRANGE , SET_GUESSES , SET_VERBOSITY , SET_COLUMNS , SET_TOLERANCE , SET_FULLOUTPUT , SET_STEPSIZE , SET_LIMITS , SET_PIPE , SET_FIXVALUE , SET_MAJOR_ORDER , N_OPTIONS }

Functions
void	report (double res, double *a, long pass, long n_eval, long n_dimen)

void	setupOutputFile (SDDS_DATASET OutputTable, long xIndex, long yIndex, long syIndex, long fitIndex, long residualIndex, long fullOutput, char output, SDDS_DATASET InputTable, char xName, char yName, char *syName, short columnMajorOrder)

long	computeStartingPoint (double a, double da, double x, double y, int64_t n, unsigned long guessFlags, double sigmaGuess, double meanGuess, double baselineGuess, double heightGuess, double stepSize)

double	fitFunction (double a, long invalid)

int64_t	makeFilteredCopy (double xFit, double yFit, double *syFit, double x, double y, double sy, int64_t n, double lower, double upper)

int	main (int argc, char **argv)

Variables
char *	option [N_OPTIONS]

static char *	USAGE

static double *	xDataFit = NULL

static double *	yDataFit = NULL

static double *	syDataFit = NULL

static int64_t	nDataFit = 0

Detailed Description

Performs a Gaussian fit on input data.

This program fits the data to the Gaussian function:

y(n) = baseline + height * exp(-0.5 * (x(n) - mean)^2 / sigma^2)

The program supports various options for specifying fit ranges, initial guesses, verbosity levels, output formats, and more. It processes input data from a file or via a pipe and outputs the fitted results to a file or pipe.

Usage

To execute the Gaussian fitting program, use the following command structure:

sddsgfit [<inputfile>] [<outputfile>] [options]

Options

The program accepts the following command-line options:

-pipe=[input][,output] -columns=<x-name>,<y-name>[,ySigma=<sy-name>] -fitRange=<lower>|<parameter-name>,<upper>|<parameter-name> -fullOutput -verbosity=<integer> -stepSize=<factor> -tolerance=

-guesses=[baseline=

|<parameter-name>][,mean=

|<parameter-name>][,height=

|<parameter-name>][,sigma=

|<parameter-name>] -fixValue=[baseline=

|<parameter-name>][,mean=

|<parameter-name>][,height=

|<parameter-name>][,sigma=

|<parameter-name>] -limits=[evaluations=<number>][,passes=<number>] -majorOrder=row|column

Copyright

(c) 2002 The University of Chicago, as Operator of Argonne National Laboratory.
(c) 2002 The Regents of the University of California, as Operator of Los Alamos National Laboratory.

License: This file is distributed under the terms of the Software License Agreement found in the file LICENSE included with this distribution.

Author: M. Borland, C. Saunders, R. Soliday, H. Shang

Definition in file sddsgfit.c.

Macro Definition Documentation

◆ BASELINE_INDEX

#define BASELINE_INDEX 0

Definition at line 123 of file sddsgfit.c.

◆ FIX_BASELINE_GIVEN

#define FIX_BASELINE_GIVEN (0x0001 << 4)

Definition at line 115 of file sddsgfit.c.

◆ FIX_HEIGHT_GIVEN

#define FIX_HEIGHT_GIVEN (0x0002 << 4)

Definition at line 117 of file sddsgfit.c.

◆ FIX_MEAN_GIVEN

#define FIX_MEAN_GIVEN (0x0004 << 4)

Definition at line 119 of file sddsgfit.c.

◆ FIX_SIGMA_GIVEN

#define FIX_SIGMA_GIVEN (0x0008 << 4)

Definition at line 121 of file sddsgfit.c.

◆ GUESS_BASELINE_GIVEN

#define GUESS_BASELINE_GIVEN 0x0001

Definition at line 114 of file sddsgfit.c.

◆ GUESS_HEIGHT_GIVEN

#define GUESS_HEIGHT_GIVEN 0x0002

Definition at line 116 of file sddsgfit.c.

◆ GUESS_MEAN_GIVEN

#define GUESS_MEAN_GIVEN 0x0004

Definition at line 118 of file sddsgfit.c.

◆ GUESS_SIGMA_GIVEN

#define GUESS_SIGMA_GIVEN 0x0008

Definition at line 120 of file sddsgfit.c.

◆ HEIGHT_INDEX

#define HEIGHT_INDEX 1

Definition at line 124 of file sddsgfit.c.

◆ MEAN_INDEX

#define MEAN_INDEX 2

Definition at line 125 of file sddsgfit.c.

◆ SIGMA_INDEX

#define SIGMA_INDEX 3

Definition at line 126 of file sddsgfit.c.

Enumeration Type Documentation

◆ option_type

enum option_type

Definition at line 53 of file sddsgfit.c.

                 {
  SET_FITRANGE,
  SET_GUESSES,
  SET_VERBOSITY,
  SET_COLUMNS,
  SET_TOLERANCE,
  SET_FULLOUTPUT,
  SET_STEPSIZE,
  SET_LIMITS,
  SET_PIPE,
  SET_FIXVALUE,
  SET_MAJOR_ORDER,
  N_OPTIONS
};

Function Documentation

◆ computeStartingPoint()

long computeStartingPoint	(	double *	a,
		double *	da,
		double *	x,
		double *	y,
		int64_t	n,
		unsigned long	guessFlags,
		double	sigmaGuess,
		double	meanGuess,
		double	baselineGuess,
		double	heightGuess,
		double	stepSize )

Definition at line 696 of file sddsgfit.c.

                                           {
  double xhalf, dhalf, ymin, ymax, xcenter, tmp, xmax, xmin;
  int64_t i;
 
  if (n < 5) {
    return 0;
  }
 
  /* first find maximum y value and corresponding x value, plus max x */
  xcenter = 0;
  ymax = xmax = -DBL_MAX;
  ymin = xmin = DBL_MAX;
  for (i = 0; i < n; i++) {
    if (xmax < fabs(x[i])) {
      xmax = fabs(x[i]);
    }
    if (xmin > fabs(x[i])) { /* Corrected comparison from < to > */
      xmin = fabs(x[i]);
    }
    if (ymax < y[i]) {
      ymax = y[i];
      xcenter = x[i];
    }
    if (ymin > y[i]) {
      ymin = y[i];
    }
  }
 
  /*  now find approximate half-max point */
  xhalf = 0;
  dhalf = DBL_MAX;
  for (i = 0; i < n; i++) {
    tmp = fabs((fabs(y[i] - ymax) / (ymax - ymin)) - 0.5);
    if (tmp < dhalf) {
      xhalf = x[i];
      dhalf = tmp;
    }
  }
  if (dhalf != DBL_MAX) {
    a[SIGMA_INDEX] = fabs(xhalf - xcenter) / 1.177; /* starting sigma */
  } else {
    a[SIGMA_INDEX] = xmax - xmin;
  }
  a[MEAN_INDEX] = xcenter;       /* starting mean */
  a[BASELINE_INDEX] = ymin;      /* starting baseline */
  a[HEIGHT_INDEX] = ymax - ymin; /* starting height */
 
  if (guessFlags & (GUESS_SIGMA_GIVEN + FIX_SIGMA_GIVEN)) {
    a[SIGMA_INDEX] = sigmaGuess;
  }
  if (guessFlags & (GUESS_MEAN_GIVEN + FIX_MEAN_GIVEN)) {
    a[MEAN_INDEX] = meanGuess;
  }
  if (guessFlags & (GUESS_BASELINE_GIVEN + FIX_BASELINE_GIVEN)) {
    a[BASELINE_INDEX] = baselineGuess;
  }
  if (guessFlags & (GUESS_HEIGHT_GIVEN + FIX_HEIGHT_GIVEN)) {
    a[HEIGHT_INDEX] = heightGuess;
  }
 
  /* step sizes */
  for (i = 0; i < 4; i++) {
    if (!(da[i] = a[i] * stepSize))
      da[i] = stepSize;
  }
 
  return 1;
}

◆ fitFunction()

double fitFunction	(	double *	a,
		long *	invalid )

Definition at line 659 of file sddsgfit.c.

                                             {
  double sum, tmp, mean, sigma, base, norm;
  int64_t i;
 
  *invalid = 0;
  sigma = a[SIGMA_INDEX];
  mean = a[MEAN_INDEX];
  base = a[BASELINE_INDEX];
  norm = a[HEIGHT_INDEX];
 
  if (!syDataFit) {
    for (i = sum = 0; i < nDataFit; i++) {
      tmp = (xDataFit[i] - mean) / sigma;
      tmp = yDataFit[i] - base - norm * exp(-sqr(tmp) / 2);
      sum += sqr(tmp);
    }
    return (sum / nDataFit);
  } else {
    for (i = sum = 0; i < nDataFit; i++) {
      tmp = (xDataFit[i] - mean) / sigma;
      tmp = (yDataFit[i] - base - norm * exp(-sqr(tmp) / 2)) / syDataFit[i];
      sum += sqr(tmp);
    }
    return (sum / nDataFit);
  }
}

◆ main()

int main	(	int	argc,
		char **	argv )

Definition at line 128 of file sddsgfit.c.

                                {
  double *xData = NULL, *yData = NULL, *syData = NULL;
  int64_t nData = 0;
  SDDS_DATASET InputTable, OutputTable;
  SCANNED_ARG *s_arg;
  long i_arg;
  char *input, *output, *xName, *yName, *syName;
  long xIndex, yIndex, fitIndex, residualIndex, syIndex_out, retval;
  double *fitData, *residualData, rmsResidual, chiSqr, sigLevel;
  unsigned long guessFlags, dummyFlags, pipeFlags, majorOrderFlag;
  double sigmaGuess, meanGuess, baselineGuess, heightGuess;
  double tolerance, stepSize;
  double a[4], da[4], lower, upper, result;
  double aLow[4], aHigh[4];
  short disable[4] = {0, 0, 0, 0};
  int32_t nEvalMax = 5000, nPassMax = 100;
  long nEval, verbosity, fullOutput = 0;
  int64_t i;
  short columnMajorOrder = -1;
  char *lowerPar, *upperPar, *baselineGuessPar, *sigmaGuessPar, *meanGuessPar, *heightGuessPar;
 
  SDDS_RegisterProgramName(argv[0]);
  argc = scanargs(&s_arg, argc, argv);
  if (argc < 2 || argc > (2 + N_OPTIONS)) {
    bomb(NULL, USAGE);
  }
 
  for (i = 0; i < 4; i++) {
    aLow[i] = -(aHigh[i] = DBL_MAX);
  }
  aLow[SIGMA_INDEX] = 0;
  input = output = NULL;
  stepSize = 1e-2;
  tolerance = 1e-8;
  verbosity = 0;
  guessFlags = sigmaGuess = heightGuess = baselineGuess = meanGuess = pipeFlags = 0;
  xName = yName = syName = NULL;
  lower = upper = 0;
  lowerPar = upperPar = sigmaGuessPar = heightGuessPar = baselineGuessPar = meanGuessPar = NULL;
 
  for (i_arg = 1; i_arg < argc; i_arg++) {
    if (s_arg[i_arg].arg_type == OPTION) {
      switch (match_string(s_arg[i_arg].list[0], option, N_OPTIONS, 0)) {
      case SET_MAJOR_ORDER:
        majorOrderFlag = 0;
        s_arg[i_arg].n_items--;
        if (s_arg[i_arg].n_items > 0 && (!scanItemList(&majorOrderFlag, s_arg[i_arg].list + 1,
                                                       &s_arg[i_arg].n_items, 0, "row", -1, NULL, 0, SDDS_ROW_MAJOR_ORDER,
                                                       "column", -1, NULL, 0, SDDS_COLUMN_MAJOR_ORDER, NULL))) {
          SDDS_Bomb("invalid -majorOrder syntax/values");
        }
        if (majorOrderFlag & SDDS_COLUMN_MAJOR_ORDER) {
          columnMajorOrder = 1;
        } else if (majorOrderFlag & SDDS_ROW_MAJOR_ORDER) {
          columnMajorOrder = 0;
        }
        break;
 
      case SET_FITRANGE:
        if (s_arg[i_arg].n_items != 3) {
          SDDS_Bomb("incorrect -fitRange syntax");
        }
        if (s_arg[i_arg].list[1][0] == '@') {
          cp_str(&lowerPar, s_arg[i_arg].list[1] + 1);
        } else if (sscanf(s_arg[i_arg].list[1], "%lf", &lower) != 1) {
          SDDS_Bomb("invalid fitRange lower value provided");
        }
        if (s_arg[i_arg].list[2][0] == '@') {
          cp_str(&upperPar, s_arg[i_arg].list[2] + 1);
        } else if (sscanf(s_arg[i_arg].list[2], "%lf", &upper) != 1) {
          SDDS_Bomb("invalid fitRange upper value provided");
        }
        break;
 
      case SET_TOLERANCE:
        if (s_arg[i_arg].n_items != 2 ||
            sscanf(s_arg[i_arg].list[1], "%lf", &tolerance) != 1 ||
            tolerance <= 0) {
          SDDS_Bomb("incorrect -tolerance syntax");
        }
        break;
 
      case SET_STEPSIZE:
        if (s_arg[i_arg].n_items != 2 ||
            sscanf(s_arg[i_arg].list[1], "%lf", &stepSize) != 1 ||
            stepSize <= 0) {
          SDDS_Bomb("incorrect -stepSize syntax");
        }
        break;
 
      case SET_VERBOSITY:
        if (s_arg[i_arg].n_items != 2 ||
            sscanf(s_arg[i_arg].list[1], "%ld", &verbosity) != 1) {
          SDDS_Bomb("incorrect -verbosity syntax");
        }
        break;
 
      case SET_GUESSES:
        if (s_arg[i_arg].n_items < 2) {
          SDDS_Bomb("incorrect -guesses syntax");
        }
        s_arg[i_arg].n_items -= 1;
        dummyFlags = guessFlags;
        if (!scanItemList(&guessFlags, s_arg[i_arg].list + 1, &s_arg[i_arg].n_items, 0,
                          "baseline", SDDS_STRING, &baselineGuessPar, 1, GUESS_BASELINE_GIVEN,
                          "height", SDDS_STRING, &heightGuessPar, 1, GUESS_HEIGHT_GIVEN,
                          "mean", SDDS_STRING, &meanGuessPar, 1, GUESS_MEAN_GIVEN,
                          "sigma", SDDS_STRING, &sigmaGuessPar, 1, GUESS_SIGMA_GIVEN, NULL)) {
          SDDS_Bomb("invalid -guesses syntax");
        }
        if (baselineGuessPar) {
          if (baselineGuessPar[0] == '@') {
            baselineGuessPar++;
          } else {
            if (sscanf(baselineGuessPar, "%lf", &baselineGuess) != 1) {
              SDDS_Bomb("Invalid baseline guess value provided.");
            }
            free(baselineGuessPar);
            baselineGuessPar = NULL;
          }
        }
        if (heightGuessPar) {
          if (heightGuessPar[0] == '@') {
            heightGuessPar++;
          } else {
            if (sscanf(heightGuessPar, "%lf", &heightGuess) != 1) {
              SDDS_Bomb("Invalid height guess value provided.");
            }
            free(heightGuessPar);
            heightGuessPar = NULL;
          }
        }
        if (meanGuessPar) {
          if (meanGuessPar[0] == '@') {
            meanGuessPar++;
          } else {
            if (sscanf(meanGuessPar, "%lf", &meanGuess) != 1) {
              SDDS_Bomb("Invalid mean guess value provided.");
            }
            free(meanGuessPar);
            meanGuessPar = NULL;
          }
        }
        if (sigmaGuessPar) {
          if (sigmaGuessPar[0] == '@') {
            sigmaGuessPar++;
          } else {
            if (sscanf(sigmaGuessPar, "%lf", &sigmaGuess) != 1) {
              SDDS_Bomb("Invalid sigma guess value provided.");
            }
            free(sigmaGuessPar);
            sigmaGuessPar = NULL;
          }
        }
        if ((dummyFlags >> 4) & guessFlags) {
          SDDS_Bomb("can't have -fixedValue and -guesses for the same item");
        }
        guessFlags |= dummyFlags;
        break;
 
      case SET_FIXVALUE:
        if (s_arg[i_arg].n_items < 2) {
          SDDS_Bomb("incorrect -fixValue syntax");
        }
        s_arg[i_arg].n_items -= 1;
        dummyFlags = guessFlags;
        if (!scanItemList(&guessFlags, s_arg[i_arg].list + 1, &s_arg[i_arg].n_items, 0,
                          "baseline", SDDS_STRING, &baselineGuessPar, 1, FIX_BASELINE_GIVEN,
                          "height", SDDS_STRING, &heightGuessPar, 1, FIX_HEIGHT_GIVEN,
                          "mean", SDDS_STRING, &meanGuessPar, 1, FIX_MEAN_GIVEN,
                          "sigma", SDDS_STRING, &sigmaGuessPar, 1, FIX_SIGMA_GIVEN, NULL)) {
          SDDS_Bomb("invalid -fixValue syntax");
        }
        if (dummyFlags & (guessFlags >> 4)) {
          SDDS_Bomb("can't have -fixValue and -guesses for the same item");
        }
        guessFlags |= dummyFlags;
        if (baselineGuessPar) {
          if (baselineGuessPar[0] == '@') {
            baselineGuessPar++;
          } else {
            if (sscanf(baselineGuessPar, "%lf", &baselineGuess) != 1) {
              SDDS_Bomb("Invalid baseline guess value provided.");
            }
            free(baselineGuessPar);
            baselineGuessPar = NULL;
          }
        }
        if (heightGuessPar) {
          if (heightGuessPar[0] == '@') {
            heightGuessPar++;
          } else {
            if (sscanf(heightGuessPar, "%lf", &heightGuess) != 1) {
              SDDS_Bomb("Invalid height guess value provided.");
            }
            free(heightGuessPar);
            heightGuessPar = NULL;
          }
        }
        if (meanGuessPar) {
          if (meanGuessPar[0] == '@') {
            meanGuessPar++;
          } else {
            if (sscanf(meanGuessPar, "%lf", &meanGuess) != 1) {
              SDDS_Bomb("Invalid mean guess value provided.");
            }
            free(meanGuessPar);
            meanGuessPar = NULL;
          }
        }
        if (sigmaGuessPar) {
          if (sigmaGuessPar[0] == '@') {
            sigmaGuessPar++;
          } else {
            if (sscanf(sigmaGuessPar, "%lf", &sigmaGuess) != 1) {
              SDDS_Bomb("Invalid sigma guess value provided.");
            }
            free(sigmaGuessPar);
            sigmaGuessPar = NULL;
          }
        }
        break;
 
      case SET_COLUMNS:
        if (s_arg[i_arg].n_items != 3 && s_arg[i_arg].n_items != 4) {
          SDDS_Bomb("invalid -columns syntax");
        }
        xName = s_arg[i_arg].list[1];
        yName = s_arg[i_arg].list[2];
        s_arg[i_arg].n_items -= 3;
        if (!scanItemList(&dummyFlags, s_arg[i_arg].list + 3, &s_arg[i_arg].n_items, 0,
                          "ysigma", SDDS_STRING, &syName, 1, 0, NULL)) {
          SDDS_Bomb("invalid -columns syntax");
        }
        break;
 
      case SET_FULLOUTPUT:
        fullOutput = 1;
        break;
 
      case SET_LIMITS:
        if (s_arg[i_arg].n_items < 2) {
          SDDS_Bomb("incorrect -limits syntax");
        }
        s_arg[i_arg].n_items -= 1;
        if (!scanItemList(&dummyFlags, s_arg[i_arg].list + 1, &s_arg[i_arg].n_items, 0,
                          "evaluations", SDDS_LONG, &nEvalMax, 1, 0,
                          "passes", SDDS_LONG, &nPassMax, 1, 0, NULL) ||
            nEvalMax <= 0 || nPassMax <= 0) {
          SDDS_Bomb("invalid -limits syntax");
        }
        break;
 
      case SET_PIPE:
        if (!processPipeOption(s_arg[i_arg].list + 1, s_arg[i_arg].n_items - 1, &pipeFlags)) {
          SDDS_Bomb("invalid -pipe syntax");
        }
        break;
 
      default:
        fprintf(stderr, "error: unknown/ambiguous option: %s\n", s_arg[i_arg].list[0]);
        exit(EXIT_FAILURE);
        break;
      }
    } else {
      if (input == NULL) {
        input = s_arg[i_arg].list[0];
      } else if (output == NULL) {
        output = s_arg[i_arg].list[0];
      } else {
        SDDS_Bomb("too many filenames");
      }
    }
  }
 
  processFilenames("sddsgfit", &input, &output, pipeFlags, 0, NULL);
 
  for (i = 0; i < 4; i++) {
    if ((guessFlags >> 4) & (1 << i)) {
      disable[i] = 1;
    }
  }
 
  if (!xName || !yName) {
    SDDS_Bomb("-columns option must be given");
  }
 
  if (!SDDS_InitializeInput(&InputTable, input)) {
    SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
  }
  if (!SDDS_FindColumn(&InputTable, FIND_NUMERIC_TYPE, xName, NULL) ||
      !SDDS_FindColumn(&InputTable, FIND_NUMERIC_TYPE, yName, NULL) ||
      (syName && !SDDS_FindColumn(&InputTable, FIND_NUMERIC_TYPE, syName, NULL))) {
    SDDS_Bomb("one or more of the given data columns is nonexistent or nonnumeric");
  }
 
  setupOutputFile(&OutputTable, &xIndex, &yIndex, &syIndex_out, &fitIndex, &residualIndex,
                  fullOutput, output, &InputTable, xName, yName, syName, columnMajorOrder);
 
  while ((retval = SDDS_ReadPage(&InputTable)) > 0) {
    xData = yData = syData = NULL;
    fitData = residualData = NULL;
    if (!(xData = SDDS_GetColumnInDoubles(&InputTable, xName)) ||
        !(yData = SDDS_GetColumnInDoubles(&InputTable, yName))) {
      SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
    }
    if (syName && !(syData = SDDS_GetColumnInDoubles(&InputTable, syName))) {
      SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
    }
    if (lowerPar && !SDDS_GetParameterAsDouble(&InputTable, lowerPar, &lower)) {
      SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
    }
    if (upperPar && !SDDS_GetParameterAsDouble(&InputTable, upperPar, &upper)) {
      SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
    }
    if (baselineGuessPar && !SDDS_GetParameterAsDouble(&InputTable, baselineGuessPar, &baselineGuess)) {
      SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
    }
    if (heightGuessPar && !SDDS_GetParameterAsDouble(&InputTable, heightGuessPar, &heightGuess)) {
      SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
    }
    if (meanGuessPar && !SDDS_GetParameterAsDouble(&InputTable, meanGuessPar, &meanGuess)) {
      SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
    }
    if (sigmaGuessPar && !SDDS_GetParameterAsDouble(&InputTable, sigmaGuessPar, &sigmaGuess)) {
      SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
    }
 
    if ((nData = SDDS_CountRowsOfInterest(&InputTable)) < 5) {
      continue;
    }
    if (lower < upper) {
      if ((nDataFit = makeFilteredCopy(&xDataFit, &yDataFit, &syDataFit, xData, yData, syData, nData, lower, upper)) < 5) {
        continue;
      }
    } else {
      xDataFit = xData;
      yDataFit = yData;
      syDataFit = syData;
      nDataFit = nData;
    }
 
    if (!computeStartingPoint(a, da, xDataFit, yDataFit, nDataFit, guessFlags, sigmaGuess, meanGuess,
                              baselineGuess, heightGuess, stepSize)) {
      fprintf(stderr, "error: couldn't compute starting point for page %ld--skipping\n", retval);
      continue;
    }
    if (verbosity > 2) {
      fprintf(stderr, "starting values:  sigma=%.6e  mean=%.6e  baseline=%.6e  height=%.6e\n",
              a[SIGMA_INDEX], a[MEAN_INDEX], a[BASELINE_INDEX], a[HEIGHT_INDEX]);
    }
    if (verbosity > 3) {
      fprintf(stderr, "starting steps:   sigma=%.6e  mean=%.6e  baseline=%.6e  height=%.6e\n",
              da[SIGMA_INDEX], da[MEAN_INDEX], da[BASELINE_INDEX], da[HEIGHT_INDEX]);
    }
 
    nEval = simplexMin(&result, a, da, aLow, aHigh, disable, 4, -DBL_MAX, tolerance,
                       fitFunction, (verbosity > 0 ? report : NULL),
                       nEvalMax, nPassMax, 12, 3, 1.0, 0);
    if (xData != xDataFit) {
      free(xDataFit);
      free(yDataFit);
      if (syDataFit) {
        free(syDataFit);
      }
    }
 
    if (verbosity > 3) {
      fprintf(stderr, "%ld evaluations of fit function required, giving result %e\n", nEval, result);
    }
 
    fitData = trealloc(fitData, sizeof(*fitData) * nData);
    residualData = trealloc(residualData, sizeof(*residualData) * nData);
    for (i = 0, result = 0; i < nData; i++) {
      fitData[i] = a[BASELINE_INDEX] + a[HEIGHT_INDEX] * exp(-ipow((xData[i] - a[MEAN_INDEX]) / a[SIGMA_INDEX], 2) / 2);
      residualData[i] = yData[i] - fitData[i];
      result += sqr(residualData[i]);
    }
    rmsResidual = sqrt(result / nData);
    if (syData) {
      for (i = 0, chiSqr = 0; i < nData; i++) {
        chiSqr += sqr(residualData[i] / syData[i]);
      }
    } else {
      double sy2 = result / (nData - 4);
      for (i = 0, chiSqr = 0; i < nData; i++) {
        chiSqr += sqr(residualData[i]) / sy2;
      }
    }
    sigLevel = ChiSqrSigLevel(chiSqr, nData - 4);
    if (verbosity > 0) {
      fprintf(stderr, "sigma: %.15e\nmean: %.15e\nbaseline: %.15e\nheight: %.15e\n",
              a[SIGMA_INDEX], a[MEAN_INDEX], a[BASELINE_INDEX], a[HEIGHT_INDEX]);
    }
    if (verbosity > 1) {
      if (syData) {
        fprintf(stderr, "Significance level: %.5e\n", sigLevel);
      }
      fprintf(stderr, "RMS deviation: %.15e\n", rmsResidual);
    }
 
    if (!SDDS_StartPage(&OutputTable, nData) ||
        !SDDS_CopyParameters(&OutputTable, &InputTable) ||
        !SDDS_SetColumn(&OutputTable, SDDS_SET_BY_INDEX, xData, nData, xIndex) ||
        !SDDS_SetColumn(&OutputTable, SDDS_SET_BY_INDEX, fitData, nData, fitIndex) ||
        !SDDS_SetParameters(&OutputTable, SDDS_PASS_BY_VALUE | SDDS_SET_BY_NAME,
                            "gfitSigma", a[SIGMA_INDEX],
                            "gfitMean", a[MEAN_INDEX],
                            "gfitBaseline", a[BASELINE_INDEX],
                            "gfitHeight", a[HEIGHT_INDEX],
                            "gfitRmsResidual", rmsResidual,
                            "gfitSigLevel", sigLevel, NULL) ||
        (fullOutput &&
         (!SDDS_SetColumn(&OutputTable, SDDS_SET_BY_INDEX, yData, nData, yIndex) ||
          !SDDS_SetColumn(&OutputTable, SDDS_SET_BY_INDEX, residualData, nData, residualIndex) ||
          (syName && !SDDS_SetColumn(&OutputTable, SDDS_SET_BY_INDEX, syData, nData, syIndex_out)))) ||
        !SDDS_WritePage(&OutputTable)) {
      SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
    }
 
    if (xData) {
      free(xData);
    }
    if (yData) {
      free(yData);
    }
    if (syData) {
      free(syData);
    }
    if (fitData) {
      free(fitData);
    }
    if (residualData) {
      free(residualData);
    }
  }
 
  if (!SDDS_Terminate(&InputTable) || !SDDS_Terminate(&OutputTable)) {
    SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
    exit(EXIT_FAILURE);
  }
  if (lowerPar) {
    free(lowerPar);
  }
  if (upperPar) {
    free(upperPar);
  }
  if (baselineGuessPar) {
    free(baselineGuessPar);
  }
  if (heightGuessPar) {
    free(heightGuessPar);
  }
  if (sigmaGuessPar) {
    free(sigmaGuessPar);
  }
  free_scanargs(&s_arg, argc);
  return EXIT_SUCCESS;
}

◆ makeFilteredCopy()

int64_t makeFilteredCopy	(	double **	xFit,
		double **	yFit,
		double **	syFit,
		double *	x,
		double *	y,
		double *	sy,
		int64_t	n,
		double	lower,
		double	upper )

Definition at line 767 of file sddsgfit.c.

                                                                            {
  int64_t i, j;
 
  if (!(*xFit = (double *)malloc(sizeof(**xFit) * n)) ||
      !(*yFit = (double *)malloc(sizeof(**yFit) * n)) ||
      (sy && !(*syFit = (double *)malloc(sizeof(**syFit) * n)))) {
    return 0;
  }
 
  for (i = j = 0; i < n; i++) {
    if (x[i] < lower || x[i] > upper) {
      continue;
    }
    (*xFit)[j] = x[i];
    (*yFit)[j] = y[i];
    if (sy) {
      (*syFit)[j] = sy[i];
    }
    j++;
  }
  return j;
}

◆ report()

void report	(	double	res,
		double *	a,
		long	pass,
		long	n_eval,
		long	n_dimen )

Definition at line 686 of file sddsgfit.c.

                                                                      {
  long i;
 
  fprintf(stderr, "pass %ld, after %ld evaluations: result = %.16e\na = ", pass, nEval, y);
  for (i = 0; i < n_dimen; i++) {
    fprintf(stderr, "%.8e ", x[i]);
  }
  fputc('\n', stderr);
}

◆ setupOutputFile()

void setupOutputFile	(	SDDS_DATASET *	OutputTable,
		long *	xIndex,
		long *	yIndex,
		long *	syIndex,
		long *	fitIndex,
		long *	residualIndex,
		long	fullOutput,
		char *	output,
		SDDS_DATASET *	InputTable,
		char *	xName,
		char *	yName,
		char *	syName,
		short	columnMajorOrder )

Definition at line 588 of file sddsgfit.c.

                                                                                     {
  char *name, *yUnits, *description, *xUnits;
  int32_t typeValue = SDDS_DOUBLE;
  static char *residualNamePart = "Residual";
  static char *residualDescriptionPart = "Residual of Gaussian fit to ";
 
  if (!SDDS_InitializeOutput(OutputTable, SDDS_BINARY, 0, NULL, "sddsgfit output", output) ||
      !SDDS_TransferColumnDefinition(OutputTable, InputTable, xName, NULL) ||
      !SDDS_ChangeColumnInformation(OutputTable, "type", &typeValue, SDDS_BY_NAME, xName) ||
      (*xIndex = SDDS_GetColumnIndex(OutputTable, xName)) < 0 ||
      !SDDS_GetColumnInformation(InputTable, "units", &xUnits, SDDS_BY_NAME, xName) ||
      !SDDS_GetColumnInformation(InputTable, "units", &yUnits, SDDS_BY_NAME, yName)) {
    SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
  }
 
  if (columnMajorOrder != -1) {
    OutputTable->layout.data_mode.column_major = columnMajorOrder;
  } else {
    OutputTable->layout.data_mode.column_major = InputTable->layout.data_mode.column_major;
  }
 
  name = tmalloc(sizeof(*name) * (strlen(yName) + strlen(residualNamePart) + 1));
  description = tmalloc(sizeof(*description) * (strlen(yName) + strlen(residualDescriptionPart) + 1));
 
  if (fullOutput) {
    if (!SDDS_TransferColumnDefinition(OutputTable, InputTable, yName, NULL) ||
        !SDDS_ChangeColumnInformation(OutputTable, "type", &typeValue, SDDS_BY_NAME, yName) ||
        (*yIndex = SDDS_GetColumnIndex(OutputTable, yName)) < 0) {
      SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
    }
    if (syName && (!SDDS_TransferColumnDefinition(OutputTable, InputTable, syName, NULL) ||
                   !SDDS_ChangeColumnInformation(OutputTable, "type", &typeValue, SDDS_BY_NAME, syName) ||
                   (*syIndex = SDDS_GetColumnIndex(OutputTable, syName)) < 0)) {
      SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
    }
    sprintf(name, "%s%s", yName, residualNamePart);
    sprintf(description, "%s%s", yName, residualDescriptionPart);
    if ((*residualIndex = SDDS_DefineColumn(OutputTable, name, NULL, yUnits, description, NULL, SDDS_DOUBLE, 0)) < 0) {
      SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
    }
  }
 
  sprintf(name, "%sFit", yName);
  sprintf(description, "Gaussian fit to %s", yName);
  if ((*fitIndex = SDDS_DefineColumn(OutputTable, name, NULL, yUnits, description, NULL, SDDS_DOUBLE, 0)) < 0) {
    SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
  }
 
  if (SDDS_DefineParameter(OutputTable, "gfitBaseline", NULL, yUnits, "Baseline from Gaussian fit",
                           NULL, SDDS_DOUBLE, 0) < 0 ||
      SDDS_DefineParameter(OutputTable, "gfitHeight", NULL, yUnits, "Height from Gaussian fit",
                           NULL, SDDS_DOUBLE, 0) < 0 ||
      SDDS_DefineParameter(OutputTable, "gfitMean", "$gm$r", xUnits, "Mean from Gaussian fit",
                           NULL, SDDS_DOUBLE, 0) < 0 ||
      SDDS_DefineParameter(OutputTable, "gfitSigma", "$gs$r", xUnits, "Sigma from Gaussian fit",
                           NULL, SDDS_DOUBLE, 0) < 0 ||
      SDDS_DefineParameter(OutputTable, "gfitRmsResidual", NULL, yUnits, "RMS residual from Gaussian fit",
                           NULL, SDDS_DOUBLE, 0) < 0 ||
      SDDS_DefineParameter(OutputTable, "gfitSigLevel", NULL, NULL, "Significance level from chi-squared test",
                           NULL, SDDS_DOUBLE, 0) < 0 ||
      !SDDS_TransferAllParameterDefinitions(OutputTable, InputTable, SDDS_TRANSFER_KEEPOLD) ||
      !SDDS_WriteLayout(OutputTable)) {
    SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
  }
 
  free(name);
  free(description);
}

Variable Documentation

◆ nDataFit

int64_t nDataFit = 0

static

Definition at line 112 of file sddsgfit.c.

◆ option

char* option[N_OPTIONS]

Initial value:

= {
  "fitrange",
  "guesses",
  "verbosity",
  "columns",
  "tolerance",
  "fulloutput",
  "stepsize",
  "limits",
  "pipe",
  "fixvalue",
  "majorOrder",
}

Definition at line 68 of file sddsgfit.c.

                          {
  "fitrange",
  "guesses",
  "verbosity",
  "columns",
  "tolerance",
  "fulloutput",
  "stepsize",
  "limits",
  "pipe",
  "fixvalue",
  "majorOrder",
};

◆ syDataFit

double * syDataFit = NULL

static

Definition at line 111 of file sddsgfit.c.

◆ USAGE

char* USAGE

static

Initial value:

=
  "Usage: sddsgfit [<inputfile>] [<outputfile>] [-pipe=[input][,output]]\n"
  "  -columns=<x-name>,<y-name>[,ySigma=<sy-name>]\n"
  "  -fitRange=<lower>|@<parameter-name>,<upper>|@<parameter-name>\n"
  "  -fullOutput\n"
  "  -verbosity=<integer>\n"
  "  -stepSize=<factor>\n"
  "  -tolerance=<value>\n"
  "  -guesses=[baseline=<value>|@<parameter-name>][,mean=<value>|@<parameter-name>]"
  "[,height=<value>|@<parameter-name>][,sigma=<value>|@<parameter-name>]\n"
  "  -fixValue=[baseline=<value>|@<parameter-name>][,mean=<value>|@<parameter-name>]"
  "[,height=<value>|@<parameter-name>][,sigma=<value>|@<parameter-name>]\n"
  "  -limits=[evaluations=<number>][,passes=<number>]\n"
  "  -majorOrder=row|column\n\n"
  "Performs a Gaussian fit of the form:\n"
  "  y = <baseline> + <height> * exp(-0.5 * (x - <mean>)^2 / <sigma>^2)\n"
  "Program by Michael Borland. (" __DATE__ " " __TIME__ ", SVN revision: " SVN_VERSION ")\n"

Definition at line 82 of file sddsgfit.c.

◆ xDataFit

double* xDataFit = NULL

static

Definition at line 111 of file sddsgfit.c.

◆ yDataFit

double * yDataFit = NULL

static

Definition at line 111 of file sddsgfit.c.

Macros

Enumerations

Functions

Variables

Detailed Description

Usage

Options

Macro Definition Documentation

◆ BASELINE_INDEX

◆ FIX_BASELINE_GIVEN

◆ FIX_HEIGHT_GIVEN

◆ FIX_MEAN_GIVEN

◆ FIX_SIGMA_GIVEN

◆ GUESS_BASELINE_GIVEN

◆ GUESS_HEIGHT_GIVEN

◆ GUESS_MEAN_GIVEN

◆ GUESS_SIGMA_GIVEN

◆ HEIGHT_INDEX

◆ MEAN_INDEX

◆ SIGMA_INDEX

Enumeration Type Documentation

◆ option_type

Function Documentation

◆ computeStartingPoint()

◆ fitFunction()

◆ main()

◆ makeFilteredCopy()

◆ report()

◆ setupOutputFile()

Variable Documentation

◆ nDataFit

◆ option

◆ syDataFit

◆ USAGE

◆ xDataFit

◆ yDataFit