sddsgfit.c File Reference

Detailed Description

A program for performing Gaussian fitting on input data.

This program fits input data to a Gaussian function using customizable options. It reads data from SDDS (Self Describing Data Set) files, processes input columns, and generates fitted results. Users can specify fit ranges, initial parameter guesses, verbosity levels, and output formats.

Usage

sddsgfit [<inputfile>] [<outputfile>]
         [-pipe=[input][,output]]
          -columns=<x-name>,<y-name>[,ySigma=<sy-name>] 
         [-fitRange=<lower>|@<parameter-name>,<upper>|@<parameter-name>] 
         [-fullOutput] 
         [-verbosity=<integer>] 
         [-stepSize=<factor>] 
         [-tolerance=<value>] 
         [-guesses=[baseline=<value>|@<parameter-name>][,mean=<value>|@<parameter-name>][,height=<value>|@<parameter-name>][,sigma=<value>|@<parameter-name>]] 
         [-fixValue=[baseline=<value>|@<parameter-name>][,mean=<value>|@<parameter-name>][,height=<value>|@<parameter-name>][,sigma=<value>|@<parameter-name>]] 
         [-limits=[evaluations=<number>][,passes=<number>]] 
         [-majorOrder=row|column]

Options

Required	Description
-columns	Specifies the column names for x and y data.

Optional	Description
-pipe	Specifies input/output via pipes.
-fitRange	Sets the range for fitting.
-fullOutput	Enables full output, including residuals.
-verbosity	Controls verbosity level (higher values give more detailed logs).
-stepSize	Sets the optimization step size.
-tolerance	Specifies the convergence tolerance.
-guesses	Provides initial parameter guesses (baseline, mean, height, sigma).
-fixValue	Fixes specific parameter values during fitting.
-limits	Sets limits for evaluations and passes.
-majorOrder	Defines the data order for output (row or column major).

Incompatibilities

• -fixValue is incompatible with -guesses for the same parameter.

The program fits data to the Gaussian function:

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

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.

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

Go to the source code of this file.

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)

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 711 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 674 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 143 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 782 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 701 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 603 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);
}