SDDS ToolKit Programs and Libraries for C and Python
•All Classes Files Functions Variables Macros Pages
Functions
sddslorentzianfit.c File Reference

Detailed Description

Perform a Lorentzian fit on data using the SDDS library.

This program reads data from an SDDS file, performs a Lorentzian fit, and writes the results to an SDDS file. The fitting function is:

\[ y(x) = baseline + height \cdot \frac{\gamma^2}{\gamma^2 + (x - center)^2} \]

Various options allow the user to customize the fit, including specifying fitting ranges, initial parameter guesses, and verbosity levels.

Usage

sddslorentzianfit [<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>][,center=<value>|@<parameter-name>][,height=<value>|@<parameter-name>][,gamma=<value>|@<parameter-name>]]
[-fixValue=[baseline=<value>|@<parameter-name>][,center=<value>|@<parameter-name>][,height=<value>|@<parameter-name>][,gamma=<value>|@<parameter-name>]]
[-limits=[evaluations=<number>][,passes=<number>]]
[-majorOrder=row|column]

Options

Required Description
-columns Specify x and y column names, and optionally y sigma.
Optional Description
-pipe Use standard input/output pipes.
-fitRange Define fitting range with bounds, which can be constants or parameters.
-fullOutput Enable detailed output including residuals.
-verbosity Set verbosity level.
-stepSize Step size factor for fitting algorithm.
-tolerance Tolerance for convergence.
-guesses Provide initial guesses for parameters (baseline, center, height, gamma).
-fixValue Fix specific parameters to provided values.
-limits Set limits on evaluations and passes.
-majorOrder Specify major order for data storage (row or column).

Incompatibilities

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

Requirements

  • For -fitRange:
    • Both lower and upper bounds are required.
  • For -limits:
    • Must specify evaluations or passes or both.
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

Definition in file sddslorentzianfit.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 gammaGuess, double centerGuess, 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 gammaGuess,
double centerGuess,
double baselineGuess,
double heightGuess,
double stepSize )

Definition at line 645 of file sddslorentzianfit.c.

647 {
648 double xhalf, dhalf, ymin, ymax, xcenter, tmp, xmax, xmin;
649 int64_t i;
650
651 if (n < 5)
652 return 0;
653
654 /* First find maximum y value and corresponding x value, plus max x */
655 xcenter = 0;
656 ymax = xmax = -DBL_MAX;
657 ymin = xmin = DBL_MAX;
658 for (i = 0; i < n; i++) {
659 if (xmax < (tmp = fabs(x[i])))
660 xmax = tmp;
661 if (xmin < tmp)
662 xmin = tmp;
663 if (ymax < y[i]) {
664 ymax = y[i];
665 xcenter = x[i];
666 }
667 if (ymin > y[i])
668 ymin = y[i];
669 }
670
671 /* Now find approximate half-max point */
672 xhalf = 0;
673 dhalf = DBL_MAX;
674 for (i = 0; i < n; i++) {
675 tmp = fabs(fabs(y[i] - ymax) / (ymax - ymin) - 0.5);
676 if (tmp < dhalf) {
677 xhalf = x[i];
678 dhalf = tmp;
679 }
680 }
681 if (dhalf != DBL_MAX)
682 a[GAMMA_INDEX] = fabs(xhalf - xcenter) / 1.177; /* Starting gamma */
683 else
684 a[GAMMA_INDEX] = xmax - xmin;
685 a[CENTER_INDEX] = xcenter; /* Starting center */
686 a[BASELINE_INDEX] = ymin; /* Starting baseline */
687 a[HEIGHT_INDEX] = ymax - ymin; /* Starting height */
688
689 if (guessFlags & (GUESS_GAMMA_GIVEN + FIX_GAMMA_GIVEN))
690 a[GAMMA_INDEX] = gammaGuess;
691 if (guessFlags & (GUESS_CENTER_GIVEN + FIX_CENTER_GIVEN))
692 a[CENTER_INDEX] = centerGuess;
693 if (guessFlags & (GUESS_BASELINE_GIVEN + FIX_BASELINE_GIVEN))
694 a[BASELINE_INDEX] = baselineGuess;
695 if (guessFlags & (GUESS_HEIGHT_GIVEN + FIX_HEIGHT_GIVEN))
696 a[HEIGHT_INDEX] = heightGuess;
697
698 /* Step sizes */
699 for (i = 0; i < 4; i++)
700 if (!(da[i] = a[i] * stepSize))
701 da[i] = stepSize;
702
703 return 1;
704}

◆ fitFunction()

double fitFunction ( double * a,
long * invalid )

Definition at line 609 of file sddslorentzianfit.c.

609 {
610 double sum, tmp, center, gamma, base, norm;
611 int64_t i;
612
613 *invalid = 0;
614 gamma = a[GAMMA_INDEX];
615 center = a[CENTER_INDEX];
616 base = a[BASELINE_INDEX];
617 norm = a[HEIGHT_INDEX];
618
619 if (!syDataFit) {
620 for (i = sum = 0; i < nDataFit; i++) {
621 tmp = (xDataFit[i] - center) / gamma;
622 tmp = norm / (1 + sqr(tmp)) + base;
623 sum += sqr(tmp - yDataFit[i]);
624 }
625 return (sum / nDataFit);
626 } else {
627 for (i = sum = 0; i < nDataFit; i++) {
628 tmp = (xDataFit[i] - center) / gamma;
629 tmp = norm / (1 + sqr(tmp)) + base;
630 sum += sqr((tmp - yDataFit[i]) / syDataFit[i]);
631 }
632 return (sum / nDataFit);
633 }
634}

◆ main()

int main ( int argc,
char ** argv )

Definition at line 171 of file sddslorentzianfit.c.

171 {
172 double *xData = NULL, *yData = NULL, *syData = NULL;
173 int64_t nData = 0;
174 SDDS_DATASET InputTable, OutputTable;
175 SCANNED_ARG *s_arg;
176 long i_arg;
177 char *input, *output, *xName, *yName, *syName;
178 long xIndex, yIndex, fitIndex, residualIndex, syIndex_col, retval;
179 double *fitData, *residualData, rmsResidual, chiSqr, sigLevel;
180 unsigned long guessFlags, dummyFlags, pipeFlags, majorOrderFlag;
181 double gammaGuess, centerGuess, baselineGuess, heightGuess;
182 double tolerance, stepSize;
183 double a[4], da[4], lower, upper, result;
184 double aLow[4], aHigh[4];
185 short disable[4] = {0, 0, 0, 0};
186 int32_t nEvalMax = 5000, nPassMax = 100;
187 long nEval, verbosity, fullOutput = 0;
188 int64_t i;
189 short columnMajorOrder = -1;
190 char *lowerPar, *upperPar, *baselineGuessPar, *gammaGuessPar, *centerGuessPar, *heightGuessPar;
191
192 SDDS_RegisterProgramName(argv[0]);
193 argc = scanargs(&s_arg, argc, argv);
194 if (argc < 2 || argc > (2 + N_OPTIONS)) {
195 fprintf(stderr, "%s", USAGE);
196 exit(EXIT_FAILURE);
197 }
198
199 for (i = 0; i < 4; i++)
200 aLow[i] = -(aHigh[i] = DBL_MAX);
201 aLow[GAMMA_INDEX] = 0;
202 input = output = NULL;
203 stepSize = 1e-2;
204 tolerance = 1e-8;
205 verbosity = 0;
206 guessFlags = gammaGuess = heightGuess = baselineGuess = centerGuess = pipeFlags = 0;
207 xName = yName = syName = NULL;
208 lower = upper = 0;
209 lowerPar = upperPar = gammaGuessPar = heightGuessPar = baselineGuessPar = centerGuessPar = NULL;
210 for (i_arg = 1; i_arg < argc; i_arg++) {
211 if (s_arg[i_arg].arg_type == OPTION) {
212 switch (match_string(s_arg[i_arg].list[0], option, N_OPTIONS, 0)) {
213 case SET_MAJOR_ORDER:
214 majorOrderFlag = 0;
215 s_arg[i_arg].n_items--;
216 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)))
217 SDDS_Bomb("invalid -majorOrder syntax/values");
218 if (majorOrderFlag & SDDS_COLUMN_MAJOR_ORDER)
219 columnMajorOrder = 1;
220 else if (majorOrderFlag & SDDS_ROW_MAJOR_ORDER)
221 columnMajorOrder = 0;
222 break;
223 case SET_FITRANGE:
224 if (s_arg[i_arg].n_items != 3)
225 SDDS_Bomb("incorrect -fitRange syntax");
226 if (s_arg[i_arg].list[1][0] == '@') {
227 cp_str(&lowerPar, s_arg[i_arg].list[1] + 1);
228 } else if (sscanf(s_arg[i_arg].list[1], "%lf", &lower) != 1)
229 SDDS_Bomb("invalid fitRange lower value provided");
230 if (s_arg[i_arg].list[2][0] == '@') {
231 cp_str(&upperPar, s_arg[i_arg].list[2] + 1);
232 } else if (sscanf(s_arg[i_arg].list[2], "%lf", &upper) != 1)
233 SDDS_Bomb("invalid fitRange upper value provided");
234 break;
235 case SET_TOLERANCE:
236 if (s_arg[i_arg].n_items != 2 || sscanf(s_arg[i_arg].list[1], "%lf", &tolerance) != 1 || tolerance <= 0)
237 SDDS_Bomb("incorrect -tolerance syntax");
238 break;
239 case SET_STEPSIZE:
240 if (s_arg[i_arg].n_items != 2 || sscanf(s_arg[i_arg].list[1], "%lf", &stepSize) != 1 || stepSize <= 0)
241 SDDS_Bomb("incorrect -stepSize syntax");
242 break;
243 case SET_VERBOSITY:
244 if (s_arg[i_arg].n_items != 2 || sscanf(s_arg[i_arg].list[1], "%ld", &verbosity) != 1)
245 SDDS_Bomb("incorrect -verbosity syntax");
246 break;
247 case SET_GUESSES:
248 if (s_arg[i_arg].n_items < 2)
249 SDDS_Bomb("incorrect -guesses syntax");
250 s_arg[i_arg].n_items -= 1;
251 dummyFlags = guessFlags;
252 if (!scanItemList(&guessFlags, s_arg[i_arg].list + 1, &s_arg[i_arg].n_items, 0,
253 "baseline", SDDS_STRING, &baselineGuessPar, 1, GUESS_BASELINE_GIVEN,
254 "height", SDDS_STRING, &heightGuessPar, 1, GUESS_HEIGHT_GIVEN,
255 "center", SDDS_STRING, &centerGuessPar, 1, GUESS_CENTER_GIVEN,
256 "gamma", SDDS_STRING, &gammaGuessPar, 1, GUESS_GAMMA_GIVEN, NULL))
257 SDDS_Bomb("invalid -guesses syntax");
258 if (baselineGuessPar) {
259 if (baselineGuessPar[0] == '@') {
260 baselineGuessPar++;
261 } else {
262 if (sscanf(baselineGuessPar, "%lf", &baselineGuess) != 1)
263 SDDS_Bomb("Invalid baseline guess value provided.");
264 free(baselineGuessPar);
265 baselineGuessPar = NULL;
266 }
267 }
268 if (heightGuessPar) {
269 if (heightGuessPar[0] == '@') {
270 heightGuessPar++;
271 } else {
272 if (sscanf(heightGuessPar, "%lf", &heightGuess) != 1)
273 SDDS_Bomb("Invalid height guess value provided.");
274 free(heightGuessPar);
275 heightGuessPar = NULL;
276 }
277 }
278 if (centerGuessPar) {
279 if (centerGuessPar[0] == '@') {
280 centerGuessPar++;
281 } else {
282 if (sscanf(centerGuessPar, "%lf", &centerGuess) != 1)
283 SDDS_Bomb("Invalid center guess value provided.");
284 free(centerGuessPar);
285 centerGuessPar = NULL;
286 }
287 }
288 if (gammaGuessPar) {
289 if (gammaGuessPar[0] == '@') {
290 gammaGuessPar++;
291 } else {
292 if (sscanf(gammaGuessPar, "%lf", &gammaGuess) != 1)
293 SDDS_Bomb("Invalid gamma guess value provided.");
294 free(gammaGuessPar);
295 gammaGuessPar = NULL;
296 }
297 }
298 if ((dummyFlags >> 4) & guessFlags)
299 SDDS_Bomb("can't have -fixValue and -guesses for the same item");
300 guessFlags |= dummyFlags;
301 break;
302 case SET_FIXVALUE:
303 if (s_arg[i_arg].n_items < 2)
304 SDDS_Bomb("incorrect -fixValue syntax");
305 s_arg[i_arg].n_items -= 1;
306 dummyFlags = guessFlags;
307 if (!scanItemList(&guessFlags, s_arg[i_arg].list + 1, &s_arg[i_arg].n_items, 0,
308 "baseline", SDDS_STRING, &baselineGuessPar, 1, FIX_BASELINE_GIVEN,
309 "height", SDDS_STRING, &heightGuessPar, 1, FIX_HEIGHT_GIVEN,
310 "center", SDDS_STRING, &centerGuessPar, 1, FIX_CENTER_GIVEN,
311 "gamma", SDDS_STRING, &gammaGuessPar, 1, FIX_GAMMA_GIVEN, NULL))
312 SDDS_Bomb("invalid -fixValue syntax");
313 if (dummyFlags & (guessFlags >> 4))
314 SDDS_Bomb("can't have -fixValue and -guesses for the same item");
315 guessFlags |= dummyFlags;
316 if (baselineGuessPar) {
317 if (baselineGuessPar[0] == '@') {
318 baselineGuessPar++;
319 } else {
320 if (sscanf(baselineGuessPar, "%lf", &baselineGuess) != 1)
321 SDDS_Bomb("Invalid baseline guess value provided.");
322 free(baselineGuessPar);
323 baselineGuessPar = NULL;
324 }
325 }
326 if (heightGuessPar) {
327 if (heightGuessPar[0] == '@') {
328 heightGuessPar++;
329 } else {
330 if (sscanf(heightGuessPar, "%lf", &heightGuess) != 1)
331 SDDS_Bomb("Invalid height guess value provided.");
332 free(heightGuessPar);
333 heightGuessPar = NULL;
334 }
335 }
336 if (centerGuessPar) {
337 if (centerGuessPar[0] == '@') {
338 centerGuessPar++;
339 } else {
340 if (sscanf(centerGuessPar, "%lf", &centerGuess) != 1)
341 SDDS_Bomb("Invalid center guess value provided.");
342 free(centerGuessPar);
343 centerGuessPar = NULL;
344 }
345 }
346 if (gammaGuessPar) {
347 if (gammaGuessPar[0] == '@') {
348 gammaGuessPar++;
349 } else {
350 if (sscanf(gammaGuessPar, "%lf", &gammaGuess) != 1)
351 SDDS_Bomb("Invalid gamma guess value provided.");
352 free(gammaGuessPar);
353 gammaGuessPar = NULL;
354 }
355 }
356 break;
357 case SET_COLUMNS:
358 if (s_arg[i_arg].n_items != 3 && s_arg[i_arg].n_items != 4)
359 SDDS_Bomb("invalid -columns syntax");
360 xName = s_arg[i_arg].list[1];
361 yName = s_arg[i_arg].list[2];
362 s_arg[i_arg].n_items -= 3;
363 if (!scanItemList(&dummyFlags, s_arg[i_arg].list + 3, &s_arg[i_arg].n_items, 0,
364 "ysigma", SDDS_STRING, &syName, 1, 0, NULL))
365 SDDS_Bomb("invalid -columns syntax");
366 break;
367 case SET_FULLOUTPUT:
368 fullOutput = 1;
369 break;
370 case SET_LIMITS:
371 if (s_arg[i_arg].n_items < 2)
372 SDDS_Bomb("incorrect -limits syntax");
373 s_arg[i_arg].n_items -= 1;
374 if (!scanItemList(&dummyFlags, s_arg[i_arg].list + 1, &s_arg[i_arg].n_items, 0,
375 "evaluations", SDDS_LONG, &nEvalMax, 1, 0,
376 "passes", SDDS_LONG, &nPassMax, 1, 0, NULL) ||
377 nEvalMax <= 0 || nPassMax <= 0)
378 SDDS_Bomb("invalid -limits syntax");
379 break;
380 case SET_PIPE:
381 if (!processPipeOption(s_arg[i_arg].list + 1, s_arg[i_arg].n_items - 1, &pipeFlags))
382 SDDS_Bomb("invalid -pipe syntax");
383 break;
384 default:
385 fprintf(stderr, "Error: Unknown or ambiguous option: %s\n", s_arg[i_arg].list[0]);
386 fprintf(stderr, "%s", USAGE);
387 exit(EXIT_FAILURE);
388 break;
389 }
390 } else {
391 if (input == NULL)
392 input = s_arg[i_arg].list[0];
393 else if (output == NULL)
394 output = s_arg[i_arg].list[0];
395 else
396 SDDS_Bomb("too many filenames");
397 }
398 }
399
400 processFilenames("sddslorentzianfit", &input, &output, pipeFlags, 0, NULL);
401
402 for (i = 0; i < 4; i++) {
403 if ((guessFlags >> 4) & (1 << i)) {
404 disable[i] = 1;
405 }
406 }
407
408 if (!xName || !yName) {
409 fprintf(stderr, "Error: -columns option must be specified.\n");
410 fprintf(stderr, "%s", USAGE);
411 exit(EXIT_FAILURE);
412 }
413
414 if (!SDDS_InitializeInput(&InputTable, input))
415 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
416 if (!SDDS_FindColumn(&InputTable, FIND_NUMERIC_TYPE, xName, NULL) ||
417 !SDDS_FindColumn(&InputTable, FIND_NUMERIC_TYPE, yName, NULL) ||
418 (syName && !SDDS_FindColumn(&InputTable, FIND_NUMERIC_TYPE, syName, NULL)))
419 SDDS_Bomb("One or more of the specified data columns do not exist or are non-numeric.");
420
421 setupOutputFile(&OutputTable, &xIndex, &yIndex, &syIndex_col, &fitIndex, &residualIndex,
422 fullOutput, output, &InputTable, xName, yName, syName, columnMajorOrder);
423
424 while ((retval = SDDS_ReadPage(&InputTable)) > 0) {
425 xData = yData = syData = NULL;
426 fitData = residualData = NULL;
427 if (!(xData = SDDS_GetColumnInDoubles(&InputTable, xName)) ||
428 !(yData = SDDS_GetColumnInDoubles(&InputTable, yName)))
429 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
430 if (syName && !(syData = SDDS_GetColumnInDoubles(&InputTable, syName)))
431 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
432 if (lowerPar && !SDDS_GetParameterAsDouble(&InputTable, lowerPar, &lower))
433 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
434 if (upperPar && !SDDS_GetParameterAsDouble(&InputTable, upperPar, &upper))
435 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
436 if (baselineGuessPar && !SDDS_GetParameterAsDouble(&InputTable, baselineGuessPar, &baselineGuess))
437 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
438 if (heightGuessPar && !SDDS_GetParameterAsDouble(&InputTable, heightGuessPar, &heightGuess))
439 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
440 if (centerGuessPar && !SDDS_GetParameterAsDouble(&InputTable, centerGuessPar, &centerGuess))
441 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
442 if (gammaGuessPar && !SDDS_GetParameterAsDouble(&InputTable, gammaGuessPar, &gammaGuess))
443 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
444
445 if ((nData = SDDS_CountRowsOfInterest(&InputTable)) < 5)
446 continue;
447 if (lower < upper) {
448 if ((nDataFit = makeFilteredCopy(&xDataFit, &yDataFit, &syDataFit, xData, yData, syData, nData, lower, upper)) < 5)
449 continue;
450 } else {
451 xDataFit = xData;
452 yDataFit = yData;
453 syDataFit = syData;
454 nDataFit = nData;
455 }
456
457 if (!computeStartingPoint(a, da, xDataFit, yDataFit, nDataFit, guessFlags, gammaGuess, centerGuess, baselineGuess, heightGuess, stepSize)) {
458 fprintf(stderr, "Error: Couldn't compute starting point for page %ld--skipping\n", retval);
459 continue;
460 }
461 if (verbosity > 2)
462 fprintf(stderr, "Starting values: gamma=%.6e center=%.6e baseline=%.6e height=%.6e\n", a[GAMMA_INDEX], a[CENTER_INDEX], a[BASELINE_INDEX], a[HEIGHT_INDEX]);
463 if (verbosity > 3)
464 fprintf(stderr, "Starting steps: gamma=%.6e center=%.6e baseline=%.6e height=%.6e\n", da[GAMMA_INDEX], da[CENTER_INDEX], da[BASELINE_INDEX], da[HEIGHT_INDEX]);
465
466 nEval = simplexMin(&result, a, da, aLow, aHigh, disable, 4, -DBL_MAX, tolerance, fitFunction, (verbosity > 0 ? report : NULL), nEvalMax, nPassMax, 12, 3, 1.0, 0);
467 if (xData != xDataFit) {
468 free(xDataFit);
469 free(yDataFit);
470 if (syDataFit)
471 free(syDataFit);
472 }
473
474 if (verbosity > 3)
475 fprintf(stderr, "%ld evaluations of fit function required, giving result %e\n", nEval, result);
476
477 fitData = trealloc(fitData, sizeof(*fitData) * nData);
478 residualData = trealloc(residualData, sizeof(*residualData) * nData);
479 for (i = result = 0; i < nData; i++) {
480 fitData[i] = a[BASELINE_INDEX] + a[HEIGHT_INDEX] / (1 + sqr((xDataFit[i] - a[CENTER_INDEX]) / a[GAMMA_INDEX]));
481 residualData[i] = yData[i] - fitData[i];
482 result += sqr(residualData[i]);
483 }
484 rmsResidual = sqrt(result / nData);
485 if (syData) {
486 for (i = chiSqr = 0; i < nData; i++)
487 chiSqr += sqr(residualData[i] / syData[i]);
488 } else {
489 double sy2;
490 sy2 = result / (nData - 4);
491 for (i = chiSqr = 0; i < nData; i++)
492 chiSqr += sqr(residualData[i]) / sy2;
493 }
494 sigLevel = ChiSqrSigLevel(chiSqr, nData - 4);
495 if (verbosity > 0) {
496 fprintf(stderr, "gamma: %.15e\ncenter: %.15e\nbaseline: %.15e\nheight: %.15e\n", a[GAMMA_INDEX], a[CENTER_INDEX], a[BASELINE_INDEX], a[HEIGHT_INDEX]);
497 }
498 if (verbosity > 1) {
499 if (syData)
500 fprintf(stderr, "Significance level: %.5e\n", sigLevel);
501 fprintf(stderr, "RMS deviation: %.15e\n", rmsResidual);
502 }
503
504 if (!SDDS_StartPage(&OutputTable, nData) ||
505 !SDDS_CopyParameters(&OutputTable, &InputTable) ||
506 !SDDS_SetColumn(&OutputTable, SDDS_SET_BY_INDEX, xData, nData, xIndex) ||
507 !SDDS_SetColumn(&OutputTable, SDDS_SET_BY_INDEX, fitData, nData, fitIndex) ||
508 !SDDS_SetParameters(&OutputTable, SDDS_PASS_BY_VALUE | SDDS_SET_BY_NAME,
509 "lorentzianfitGamma", a[GAMMA_INDEX],
510 "lorentzianfitCenter", a[CENTER_INDEX],
511 "lorentzianfitBaseline", a[BASELINE_INDEX],
512 "lorentzianfitHeight", a[HEIGHT_INDEX],
513 "lorentzianfitRmsResidual", rmsResidual,
514 "lorentzianfitSigLevel", sigLevel, NULL) ||
515 (fullOutput &&
516 (!SDDS_SetColumn(&OutputTable, SDDS_SET_BY_INDEX, yData, nData, yIndex) ||
517 !SDDS_SetColumn(&OutputTable, SDDS_SET_BY_INDEX, residualData, nData, residualIndex) ||
518 (syName && !SDDS_SetColumn(&OutputTable, SDDS_SET_BY_INDEX, syData, nData, syIndex_col)))) ||
519 !SDDS_WritePage(&OutputTable))
520 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
521 if (xData)
522 free(xData);
523 if (yData)
524 free(yData);
525 if (syData)
526 free(syData);
527 if (fitData)
528 free(fitData);
529 if (residualData)
530 free(residualData);
531 }
532 if (!SDDS_Terminate(&InputTable) || !SDDS_Terminate(&OutputTable)) {
533 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
534 exit(EXIT_FAILURE);
535 }
536 if (lowerPar)
537 free(lowerPar);
538 if (upperPar)
539 free(upperPar);
540 if (baselineGuessPar)
541 free(baselineGuessPar);
542 if (heightGuessPar)
543 free(heightGuessPar);
544 if (gammaGuessPar)
545 free(gammaGuessPar);
546 free_scanargs(&s_arg, argc);
547 return EXIT_SUCCESS;
548}
SDDS_CopyParameters
int32_t SDDS_CopyParameters(SDDS_DATASET *SDDS_target, SDDS_DATASET *SDDS_source)
Definition SDDS_copy.c:286
SDDS_StartPage
int32_t SDDS_StartPage(SDDS_DATASET *SDDS_dataset, int64_t expected_n_rows)
Definition SDDS_dataprep.c:64
SDDS_SetParameters
int32_t SDDS_SetParameters(SDDS_DATASET *SDDS_dataset, int32_t mode,...)
Definition SDDS_dataprep.c:369
SDDS_SetColumn
int32_t SDDS_SetColumn(SDDS_DATASET *SDDS_dataset, int32_t mode, void *data, int64_t rows,...)
Sets the values for one data column in the current data table of an SDDS dataset.
Definition SDDS_dataprep.c:1536
SDDS_GetParameterAsDouble
double * SDDS_GetParameterAsDouble(SDDS_DATASET *SDDS_dataset, char *parameter_name, double *memory)
Retrieves the value of a specified parameter as a double from the current data table of an SDDS datas...
Definition SDDS_extract.c:2782
SDDS_CountRowsOfInterest
int64_t SDDS_CountRowsOfInterest(SDDS_DATASET *SDDS_dataset)
Counts the number of rows marked as "of interest" in the current data table.
Definition SDDS_extract.c:364
SDDS_GetColumnInDoubles
double * SDDS_GetColumnInDoubles(SDDS_DATASET *SDDS_dataset, char *column_name)
Retrieves the data of a specified numerical column as an array of doubles, considering only rows mark...
Definition SDDS_extract.c:872
SDDS_InitializeInput
int32_t SDDS_InitializeInput(SDDS_DATASET *SDDS_dataset, char *filename)
Definition SDDS_input.c:49
SDDS_Terminate
int32_t SDDS_Terminate(SDDS_DATASET *SDDS_dataset)
Definition SDDS_input.c:1448
SDDS_ReadPage
int32_t SDDS_ReadPage(SDDS_DATASET *SDDS_dataset)
Definition SDDS_input.c:1006
SDDS_WritePage
int32_t SDDS_WritePage(SDDS_DATASET *SDDS_dataset)
Writes the current data table to the output file.
Definition SDDS_output.c:1228
SDDS_PrintErrors
void SDDS_PrintErrors(FILE *fp, int32_t mode)
Prints recorded error messages to a specified file stream.
Definition SDDS_utils.c:432
SDDS_RegisterProgramName
void SDDS_RegisterProgramName(const char *name)
Registers the executable program name for use in error messages.
Definition SDDS_utils.c:288
SDDS_FindColumn
char * SDDS_FindColumn(SDDS_DATASET *SDDS_dataset, int32_t mode,...)
Finds the first column in the SDDS dataset that matches the specified criteria.
Definition SDDS_utils.c:4312
SDDS_Bomb
void SDDS_Bomb(char *message)
Terminates the program after printing an error message and recorded errors.
Definition SDDS_utils.c:342
SDDS_STRING
#define SDDS_STRING
Identifier for the string data type.
Definition SDDStypes.h:85
SDDS_LONG
#define SDDS_LONG
Identifier for the signed 32-bit integer data type.
Definition SDDStypes.h:61
trealloc
void * trealloc(void *old_ptr, uint64_t size_of_block)
Reallocates a memory block to a new size.
Definition array.c:181
cp_str
char * cp_str(char **s, char *t)
Copies a string, allocating memory for storage.
Definition cp_str.c:28
match_string
long match_string(char *string, char **option, long n_options, long mode)
Matches a given string against an array of option strings based on specified modes.
Definition match_string.c:34
scanargs
int scanargs(SCANNED_ARG **scanned, int argc, char **argv)
Definition scanargs.c:36
processPipeOption
long processPipeOption(char **item, long items, unsigned long *flags)
Definition scanargs.c:356
processFilenames
void processFilenames(char *programName, char **input, char **output, unsigned long pipeFlags, long noWarnings, long *tmpOutputUsed)
Definition scanargs.c:390
free_scanargs
void free_scanargs(SCANNED_ARG **scanned, int argc)
Definition scanargs.c:584
scanItemList
long scanItemList(unsigned long *flags, char **item, long *items, unsigned long mode,...)
Scans a list of items and assigns values based on provided keywords and types.
Definition scanitemlist.c:41
ChiSqrSigLevel
double ChiSqrSigLevel(double ChiSquared0, long nu)
Computes the probability that a chi-squared variable exceeds a given value.
Definition sigLevel.c:64
simplexMin
long simplexMin(double *yReturn, double *xGuess, double *dxGuess, double *xLowerLimit, double *xUpperLimit, short *disable, long dimensions, double target, double tolerance, double(*func)(double *x, long *invalid), void(*report)(double ymin, double *xmin, long pass, long evals, long dims), long maxEvaluations, long maxPasses, long maxDivisions, double divisorFactor, double passRangeFactor, unsigned long flags)
Top-level convenience function for simplex-based minimization.
Definition simplex.c:472
SDDS_DATASET
Definition SDDS.h:308

◆ 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 706 of file sddslorentzianfit.c.

707 {
708 int64_t i, j;
709
710 if (!(*xFit = (double *)malloc(sizeof(**xFit) * n)) ||
711 !(*yFit = (double *)malloc(sizeof(**yFit) * n)) ||
712 (sy && !(*syFit = (double *)malloc(sizeof(**syFit) * n))))
713 return 0;
714
715 for (i = j = 0; i < n; i++) {
716 if (x[i] < lower || x[i] > upper)
717 continue;
718 (*xFit)[j] = x[i];
719 (*yFit)[j] = y[i];
720 if (sy)
721 (*syFit)[j] = sy[i];
722 j++;
723 }
724 return j;
725}

◆ report()

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

Definition at line 636 of file sddslorentzianfit.c.

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

◆ 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 550 of file sddslorentzianfit.c.

552 {
553 char *name, *yUnits, *description, *xUnits;
554 int32_t typeValue = SDDS_DOUBLE;
555 static char *residualNamePart = "Residual";
556 static char *residualDescriptionPart = "Residual of Lorentzian fit to ";
557
558 if (!SDDS_InitializeOutput(OutputTable, SDDS_BINARY, 0, NULL, "sddslorentzianfit output", output) ||
559 !SDDS_TransferColumnDefinition(OutputTable, InputTable, xName, NULL) ||
560 !SDDS_ChangeColumnInformation(OutputTable, "type", &typeValue, SDDS_BY_NAME, xName) ||
561 (*xIndex = SDDS_GetColumnIndex(OutputTable, xName)) < 0 ||
562 !SDDS_GetColumnInformation(InputTable, "units", &xUnits, SDDS_BY_NAME, xName) ||
563 !SDDS_GetColumnInformation(InputTable, "units", &yUnits, SDDS_BY_NAME, yName))
564 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
565 if (columnMajorOrder != -1)
566 OutputTable->layout.data_mode.column_major = columnMajorOrder;
567 else
568 OutputTable->layout.data_mode.column_major = InputTable->layout.data_mode.column_major;
569 name = tmalloc(sizeof(*name) * (strlen(yName) + strlen(residualNamePart) + 1));
570 description = tmalloc(sizeof(*description) * (strlen(yName) + strlen(residualDescriptionPart) + 1));
571
572 if (fullOutput) {
573 if (!SDDS_TransferColumnDefinition(OutputTable, InputTable, yName, NULL) ||
574 !SDDS_ChangeColumnInformation(OutputTable, "type", &typeValue, SDDS_BY_NAME, yName) ||
575 (*yIndex = SDDS_GetColumnIndex(OutputTable, yName)) < 0)
576 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
577 if (syName && (!SDDS_TransferColumnDefinition(OutputTable, InputTable, syName, NULL) ||
578 !SDDS_ChangeColumnInformation(OutputTable, "type", &typeValue, SDDS_BY_NAME, syName) ||
579 (*syIndex = SDDS_GetColumnIndex(OutputTable, syName)) < 0))
580 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
581 sprintf(name, "%s%s", yName, residualNamePart);
582 sprintf(description, "%s%s", yName, residualDescriptionPart);
583 if ((*residualIndex = SDDS_DefineColumn(OutputTable, name, NULL, yUnits, description, NULL, SDDS_DOUBLE, 0)) < 0)
584 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
585 }
586
587 sprintf(name, "%sFit", yName);
588 sprintf(description, "Lorentzian fit to %s", yName);
589 if ((*fitIndex = SDDS_DefineColumn(OutputTable, name, NULL, yUnits, description, NULL, SDDS_DOUBLE, 0)) < 0)
590 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
591
592 if (SDDS_DefineParameter(OutputTable, "lorentzianfitBaseline", NULL, yUnits, "Baseline from Lorentzian fit",
593 NULL, SDDS_DOUBLE, 0) < 0 ||
594 SDDS_DefineParameter(OutputTable, "lorentzianfitHeight", NULL, yUnits, "Height from Lorentzian fit",
595 NULL, SDDS_DOUBLE, 0) < 0 ||
596 SDDS_DefineParameter(OutputTable, "lorentzianfitCenter", NULL, xUnits, "Center from Lorentzian fit",
597 NULL, SDDS_DOUBLE, 0) < 0 ||
598 SDDS_DefineParameter(OutputTable, "lorentzianfitGamma", NULL, xUnits, "Gamma from Lorentzian fit",
599 NULL, SDDS_DOUBLE, 0) < 0 ||
600 SDDS_DefineParameter(OutputTable, "lorentzianfitRmsResidual", NULL, yUnits, "RMS residual from Lorentzian fit",
601 NULL, SDDS_DOUBLE, 0) < 0 ||
602 SDDS_DefineParameter(OutputTable, "lorentzianfitSigLevel", NULL, NULL, "Significance level from chi-squared test",
603 NULL, SDDS_DOUBLE, 0) < 0 ||
604 !SDDS_TransferAllParameterDefinitions(OutputTable, InputTable, SDDS_TRANSFER_KEEPOLD) ||
605 !SDDS_WriteLayout(OutputTable))
606 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
607}
SDDS_ChangeColumnInformation
int32_t SDDS_ChangeColumnInformation(SDDS_DATASET *SDDS_dataset, char *field_name, void *memory, int32_t mode,...)
Modifies a specific field in a column definition within the SDDS dataset.
Definition SDDS_info.c:364
SDDS_GetColumnInformation
int32_t SDDS_GetColumnInformation(SDDS_DATASET *SDDS_dataset, char *field_name, void *memory, int32_t mode,...)
Retrieves information about a specified column in the SDDS dataset.
Definition SDDS_info.c:41
SDDS_InitializeOutput
int32_t SDDS_InitializeOutput(SDDS_DATASET *SDDS_dataset, int32_t data_mode, int32_t lines_per_row, const char *description, const char *contents, const char *filename)
Initializes the SDDS output dataset.
Definition SDDS_output.c:601
SDDS_DefineColumn
int32_t SDDS_DefineColumn(SDDS_DATASET *SDDS_dataset, const char *name, const char *symbol, const char *units, const char *description, const char *format_string, int32_t type, int32_t field_length)
Defines a data column within the SDDS dataset.
Definition SDDS_output.c:1715
SDDS_WriteLayout
int32_t SDDS_WriteLayout(SDDS_DATASET *SDDS_dataset)
Writes the SDDS layout header to the output file.
Definition SDDS_output.c:899
SDDS_DefineParameter
int32_t SDDS_DefineParameter(SDDS_DATASET *SDDS_dataset, const char *name, const char *symbol, const char *units, const char *description, const char *format_string, int32_t type, char *fixed_value)
Defines a data parameter with a fixed string value.
Definition SDDS_output.c:1472
SDDS_TransferColumnDefinition
int32_t SDDS_TransferColumnDefinition(SDDS_DATASET *target, SDDS_DATASET *source, char *name, char *newName)
Transfers a column definition from a source dataset to a target dataset.
Definition SDDS_transfer.c:35
SDDS_TransferAllParameterDefinitions
int32_t SDDS_TransferAllParameterDefinitions(SDDS_DATASET *SDDS_target, SDDS_DATASET *SDDS_source, uint32_t mode)
Transfers all parameter definitions from a source dataset to a target dataset.
Definition SDDS_transfer.c:343
SDDS_GetColumnIndex
int32_t SDDS_GetColumnIndex(SDDS_DATASET *SDDS_dataset, char *name)
Retrieves the index of a named column in the SDDS dataset.
Definition SDDS_utils.c:1309
SDDS_DOUBLE
#define SDDS_DOUBLE
Identifier for the double data type.
Definition SDDStypes.h:37
tmalloc
void * tmalloc(uint64_t size_of_block)
Allocates a memory block of the specified size with zero initialization.
Definition array.c:59
Generated by doxygen 1.12.0