SDDS ToolKit Programs and Libraries for C and Python
Loading...
Searching...
No Matches
sddssinefit.c
Go to the documentation of this file.
1/**
2 * @file sddssinefit.c
3 * @brief Performs a sinusoidal fit on input data.
4 *
5 * @details
6 * This program reads input data from an SDDS (Self Describing Data Sets) file,
7 * fits the data to one of two sinusoidal models:
8 *
9 * @f[
10 * y(n) = a_0 + a_1 \sin(2\pi a_2 x(n) + a_3)
11 * @f]
12 *
13 * or
14 *
15 * @f[
16 * y(n) = a_0 + a_1 \sin(2\pi a_2 x(n) + a_3) + a_5 x(n)
17 * @f]
18 *
19 * or
20 *
21 * @f[
22 * y(n) = a_0 + a_1 \sin(2\pi a_2 x(n) + a_3)*exp(x(n)*a_5) + a_6 x(n)
23 * @f]
24 *
25 * Based on user-provided parameters and options, the program performs fitting and outputs
26 * the fitted data along with residuals to an output SDDS file.
27 *
28 * @section Usage
29 * ```
30 * sddssinefit [<inputfile>] [<outputfile>]
31 * [-pipe=<input>[,<output>]]
32 * [-fulloutput]
33 * -columns=<x-name>,<y-name>
34 * [-tolerance=<value>]
35 * [-limits=evaluations=<number>,passes=<number>]
36 * [-verbosity=<integer>]
37 * [-guess=constant=<constant>,factor=<factor>,frequency=<freq>,phase=<phase>,slope=<slope>,rate=<value>]
38 * [-lockFrequency]
39 * [-addSlope]
40 * [-addExponential={grow|decay}]
41 * [-majorOrder=row|column]
42 * ```
43 *
44 * @section Options
45 * | Required | Description |
46 * |---------------------------------------|---------------------------------------------------------------------------------------|
47 * | `-columns` | Specifies the x and y data column names. |
48 *
49 * | Optional | Description |
50 * |---------------------------------------|---------------------------------------------------------------------------------------|
51 * | `-pipe` | Use standard input/output for data streams. |
52 * | `-fulloutput` | Includes full output with residuals. |
53 * | `-tolerance` | Sets the tolerance for the fitting algorithm (default: 1e-6). |
54 * | `-limits=evaluations` | Sets maximum number of evaluations and passes (default: 5000 evaluations, 25 passes).|
55 * | `-verbosity` | Sets verbosity level (default: 0). |
56 * | `-guess` | Provides initial guesses for fit parameters. |
57 * | `-lockFrequency` | Locks the frequency parameter during fitting. |
58 * | `-addSlope` | Includes a slope term in the fit. |
59 * | `-addExponential` | Includes decaying expontential factor in the fit. |
60 * | `-majorOrder` | Specifies the major order for data processing. |
61 *
62 * @subsection Incompatibilities
63 * - `-lockFrequency` cannot be used with `-guess=frequency=<freq>`.
64 *
65 * @subsection Requirements
66 * - For `-guess`, at least one of the guess parameters (`constant`, `factor`, `frequency`, `phase`, `rate`, `slope`) is required.
67 *
68 * @copyright
69 * - (c) 2002 The University of Chicago, as Operator of Argonne National Laboratory.
70 * - (c) 2002 The Regents of the University of California, as Operator of Los Alamos National Laboratory.
71 *
72 * @license
73 * This file is distributed under the terms of the Software License Agreement
74 * found in the file LICENSE included with this distribution.
75 *
76 * @author
77 * M. Borland, C. Saunders, R. Soliday, L. Emery, H. Shang, N. Kuklev
78 */
79
80#include "mdb.h"
81#include "SDDS.h"
82#include "scan.h"
83
84#define EXPONENTIAL_DECAY 0
85#define EXPONENTIAL_GROW 1
86#define EXPONENTIAL_OPTIONS 2
87static char *expotentialOptions[EXPONENTIAL_OPTIONS] = {
88 "decay", "grow"
89};
90
91/* Enumeration for option types */
92enum option_type {
93 SET_TOLERANCE,
94 SET_VERBOSITY,
95 SET_CLUE,
96 SET_GUESS,
97 SET_COLUMNS,
98 SET_FULLOUTPUT,
99 SET_LIMITS,
100 SET_PIPE,
101 SET_MAJOR_ORDER,
102 SET_LOCK_FREQ,
103 SET_ADD_SLOPE,
104 SET_ADD_EXPONENTIAL,
105 N_OPTIONS
106};
107
108char *option[N_OPTIONS] = {
109 "tolerance",
110 "verbosity",
111 "clue",
112 "guess",
113 "columns",
114 "fulloutput",
115 "limits",
116 "pipe",
117 "majorOrder",
118 "lockFrequency",
119 "addSlope",
120 "addExponential",
121};
122
123static char *USAGE =
124 "sddssinefit [<inputfile>] [<outputfile>] \n"
125 " [-pipe=<input>[,<output>]]\n"
126 " [-fulloutput]\n"
127 " [-columns=<x-name>,<y-name>]\n"
128 " [-tolerance=<value>]\n"
129 " [-limits=evaluations=<number>,passes=<number>]\n"
130 " [-verbosity=<integer>]\n"
131 " [-guess=constant=<constant>,factor=<factor>,frequency=<freq>,phase=<phase>,slope=<slope>,rate=<value>]\n"
132 " [-lockFrequency]\n"
133 " [-addSlope] [-addExponential={grow|decay}\n"
134 " [-majorOrder=row|column]\n\n"
135 "Description:\n"
136 " Performs a sinusoidal fit of the form:\n"
137 " y = <constant> + <factor>*sin(2*PI*<freq>*x + <phase>)\n"
138 " y = <constant> + <factor>*sin(2*PI*<freq>*x + <phase>) + <slope>*x\n\n"
139 " or\n"
140 " y = <constant> + <factor>*sin(2*PI*<freq>*x + <phase>)*exp(<rate>*x) + <slope>*x\n\n"
141 "Options:\n"
142 " <inputfile> : Path to the input SDDS file.\n"
143 " <outputfile> : Path to the output SDDS file.\n"
144 " -pipe=<input>,<output> : Use standard input/output for data streams.\n"
145 " -fulloutput : Include full output with residuals.\n"
146 " -columns=<x-name>,<y-name> : Specify the names of the x and y data columns.\n"
147 " -tolerance=<value> : Set the tolerance for the fitting algorithm (default: 1e-6).\n"
148 " -limits=evaluations=<n>,passes=<m> : Set maximum number of evaluations and passes (default: 5000 evaluations, 25 passes).\n"
149 " -verbosity=<integer> : Set verbosity level (default: 0).\n"
150 " -guess=constant=<c>,factor=<f>,frequency=<freq>,phase=<p>,slope=<s>,rate=<r> : Provide initial guesses for fit parameters.\n"
151 " -lockFrequency : Lock the frequency parameter during fitting.\n"
152 " -addSlope : Include a slope term in the fit.\n"
153 " -addExponential : Include exponential decay of the sinusoid.\n"
154 " -majorOrder=row|column : Specify the major order for data processing.\n\n"
155 "Author:\n"
156 " Michael Borland\n"
157 " (" __DATE__ " " __TIME__ ", SVN revision: " SVN_VERSION ")\n";
158
159static double *xData, *yData;
160static int64_t nData;
161static double yMin, yMax, xMin, xMax;
162static double fit[6];
163static short disable[6] = {0,0,0,0,1,1};
164static short addSlope, addExponential;
165static double *fitData, *residualData, rmsResidual;
166
167double fitFunction(double *a, long *invalid);
168void report(double res, double *a, long pass, long n_eval, long n_dimen);
169void setupOutputFile(SDDS_DATASET *OutputTable, int32_t *xIndex, int32_t *yIndex, int32_t *fitIndex,
170 int32_t *residualIndex, char *output, long fullOutput, SDDS_DATASET *InputTable,
171 char *xName, char *yName, short columnMajorOrder, short addSlope, short addExponential);
172char *makeInverseUnits(char *units);
173
174long verbosity;
175
176#define GUESS_CONSTANT_GIVEN 0x0001
177#define GUESS_FACTOR_GIVEN 0x0002
178#define GUESS_FREQ_GIVEN 0x0004
179#define GUESS_PHASE_GIVEN 0x0008
180#define GUESS_SLOPE_GIVEN 0x0010
181#define GUESS_RATE_GIVEN 0x0020
182
183int main(int argc, char **argv) {
184 double tolerance, result;
185 int32_t nEvalMax = 5000, nPassMax = 25;
186 double a[6], da[6];
187 double alo[6], ahi[6];
188 int64_t zeroes;
189 SDDS_DATASET InputTable, OutputTable;
190 SCANNED_ARG *s_arg;
191 long i_arg, fullOutput;
192 int64_t i;
193 char *input, *output, *xName, *yName;
194 int32_t xIndex, yIndex, fitIndex, residualIndex;
195 long retval;
196 unsigned long guessFlags, dummyFlags, pipeFlags;
197 double constantGuess, factorGuess, freqGuess, phaseGuess, slopeGuess, rateGuess;
198 double firstZero, lastZero;
199 unsigned long simplexFlags = 0, majorOrderFlag;
200 short columnMajorOrder = -1;
201 short lockFreq = 0;
202
203 SDDS_RegisterProgramName(argv[0]);
204 argc = scanargs(&s_arg, argc, argv);
205 if (argc < 2 || argc > (2 + N_OPTIONS))
206 bomb(NULL, USAGE);
207
208 input = output = NULL;
209 tolerance = 1e-6;
210 verbosity = fullOutput = 0;
211 xName = yName = NULL;
212 guessFlags = 0;
213 pipeFlags = 0;
214 constantGuess = factorGuess = freqGuess = phaseGuess = slopeGuess = rateGuess = 0;
215
216 for (i_arg = 1; i_arg < argc; i_arg++) {
217 if (s_arg[i_arg].arg_type == OPTION) {
218 switch (match_string(s_arg[i_arg].list[0], option, N_OPTIONS, 0)) {
219 case SET_MAJOR_ORDER:
220 majorOrderFlag = 0;
221 s_arg[i_arg].n_items--;
222 if (s_arg[i_arg].n_items > 0 &&
223 (!scanItemList(&majorOrderFlag, s_arg[i_arg].list + 1, &s_arg[i_arg].n_items, 0,
224 "row", -1, NULL, 0, SDDS_ROW_MAJOR_ORDER,
225 "column", -1, NULL, 0, SDDS_COLUMN_MAJOR_ORDER, NULL)))
226 SDDS_Bomb("invalid -majorOrder syntax/values");
227 if (majorOrderFlag & SDDS_COLUMN_MAJOR_ORDER)
228 columnMajorOrder = 1;
229 else if (majorOrderFlag & SDDS_ROW_MAJOR_ORDER)
230 columnMajorOrder = 0;
231 break;
232 case SET_TOLERANCE:
233 if (s_arg[i_arg].n_items != 2 || sscanf(s_arg[i_arg].list[1], "%lf", &tolerance) != 1)
234 SDDS_Bomb("incorrect -tolerance syntax");
235 break;
236 case SET_VERBOSITY:
237 if (s_arg[i_arg].n_items != 2 || sscanf(s_arg[i_arg].list[1], "%ld", &verbosity) != 1)
238 SDDS_Bomb("incorrect -verbosity syntax");
239 break;
240 case SET_GUESS:
241 if (s_arg[i_arg].n_items < 2)
242 SDDS_Bomb("incorrect -guess syntax");
243 s_arg[i_arg].n_items -= 1;
244 if (!scanItemList(&guessFlags, s_arg[i_arg].list + 1, &s_arg[i_arg].n_items, 0,
245 "constant", SDDS_DOUBLE, &constantGuess, 1, GUESS_CONSTANT_GIVEN,
246 "factor", SDDS_DOUBLE, &factorGuess, 1, GUESS_FACTOR_GIVEN,
247 "frequency", SDDS_DOUBLE, &freqGuess, 1, GUESS_FREQ_GIVEN,
248 "phase", SDDS_DOUBLE, &phaseGuess, 1, GUESS_PHASE_GIVEN,
249 "slope", SDDS_DOUBLE, &slopeGuess, 1, GUESS_SLOPE_GIVEN,
250 "rate", SDDS_DOUBLE, &rateGuess, 1, GUESS_RATE_GIVEN, NULL,
251 NULL))
252 SDDS_Bomb("invalid -guess syntax");
253 break;
254 case SET_COLUMNS:
255 if (s_arg[i_arg].n_items != 3)
256 SDDS_Bomb("invalid -columns syntax");
257 xName = s_arg[i_arg].list[1];
258 yName = s_arg[i_arg].list[2];
259 break;
260 case SET_FULLOUTPUT:
261 fullOutput = 1;
262 break;
263 case SET_LOCK_FREQ:
264 lockFreq = 1;
265 break;
266 case SET_ADD_SLOPE:
267 addSlope = 1;
268 break;
269 case SET_ADD_EXPONENTIAL:
270 if (s_arg[i_arg].n_items!=2)
271 SDDS_Bomb("incorrect -addExponential syntax");
272 switch (match_string(s_arg[i_arg].list[1], expotentialOptions, EXPONENTIAL_OPTIONS, 0)) {
273 case EXPONENTIAL_DECAY:
274 addExponential = -1;
275 break;
276 case EXPONENTIAL_GROW:
277 addExponential = 1;
278 break;
279 default:
280 SDDS_Bomb("incorrect -addExponential syntax");
281 break;
282 }
283 break;
284 case SET_LIMITS:
285 if (s_arg[i_arg].n_items < 2)
286 SDDS_Bomb("incorrect -limits syntax");
287 s_arg[i_arg].n_items -= 1;
288 if (!scanItemList(&dummyFlags, s_arg[i_arg].list + 1, &s_arg[i_arg].n_items, 0,
289 "evaluations", SDDS_LONG, &nEvalMax, 1, 0,
290 "passes", SDDS_LONG, &nPassMax, 1, 0, NULL) ||
291 nEvalMax <= 0 || nPassMax <= 0)
292 SDDS_Bomb("invalid -limits syntax");
293 break;
294 case SET_PIPE:
295 if (!processPipeOption(s_arg[i_arg].list + 1, s_arg[i_arg].n_items - 1, &pipeFlags))
296 SDDS_Bomb("invalid -pipe syntax");
297 break;
298 default:
299 fprintf(stderr, "Error: Unknown or ambiguous option: %s\n", s_arg[i_arg].list[0]);
300 exit(EXIT_FAILURE);
301 break;
302 }
303 } else {
304 if (input == NULL)
305 input = s_arg[i_arg].list[0];
306 else if (output == NULL)
307 output = s_arg[i_arg].list[0];
308 else
309 SDDS_Bomb("Too many filenames provided.");
310 }
311 }
312
313 if (addSlope)
314 disable[4] = 0;
315 if (addExponential)
316 disable[5] = 0;
317
318 if ((guessFlags & GUESS_SLOPE_GIVEN) && !addSlope)
319 SDDS_Bomb("-guess=slope given but -addSlope not given");
320
321 if ((guessFlags & GUESS_RATE_GIVEN) && !addExponential)
322 SDDS_Bomb("-guess=rate given but -addExponential not given");
323
324 processFilenames("sddssinefit", &input, &output, pipeFlags, 0, NULL);
325
326 if (!xName || !yName)
327 SDDS_Bomb("-columns option must be specified.");
328
329 if (!SDDS_InitializeInput(&InputTable, input) ||
330 SDDS_GetColumnIndex(&InputTable, xName) < 0 ||
331 SDDS_GetColumnIndex(&InputTable, yName) < 0)
332 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
333
334 setupOutputFile(&OutputTable, &xIndex, &yIndex, &fitIndex, &residualIndex, output, fullOutput,
335 &InputTable, xName, yName, columnMajorOrder, addSlope, addExponential);
336
337 fitData = residualData = NULL;
338
339 alo[0] = -(ahi[0] = DBL_MAX);
340 alo[1] = alo[2] = 0;
341 ahi[1] = ahi[2] = DBL_MAX;
342 alo[3] = -(ahi[3] = PIx2);
343 if (addSlope)
344 alo[4] = -(ahi[4] = DBL_MAX);
345 if (addExponential)
346 alo[5] = -(ahi[5] = DBL_MAX);
347
348 firstZero = lastZero = 0;
349 if (verbosity>5) {
350 simplexFlags = SIMPLEX_VERBOSE_LEVEL1;
351 }
352 while ((retval = SDDS_ReadPage(&InputTable)) > 0) {
353 if (!(xData = SDDS_GetColumnInDoubles(&InputTable, xName)) ||
354 !(yData = SDDS_GetColumnInDoubles(&InputTable, yName)))
355 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
356 if ((nData = SDDS_CountRowsOfInterest(&InputTable)) < 4)
357 continue;
358
359 fitData = SDDS_Realloc(fitData, sizeof(*fitData) * nData);
360 residualData = SDDS_Realloc(residualData, sizeof(*residualData) * nData);
361
362 find_min_max(&yMin, &yMax, yData, nData);
363 find_min_max(&xMin, &xMax, xData, nData);
364 zeroes = 0;
365 for (i = 1; i < nData; i++)
366 if (yData[i] * yData[i - 1] <= 0) {
367 i++;
368 if (!zeroes)
369 firstZero = (xData[i] + xData[i - 1]) / 2;
370 else
371 lastZero = (xData[i] + xData[i - 1]) / 2;
372 zeroes++;
373 if (zeroes>5)
374 break;
375 }
376 a[0] = (yMin + yMax) / 2;
377 a[1] = (yMax - yMin) / 2;
378 if (!zeroes)
379 a[2] = 2 / fabs(xMax - xMin);
380 else
381 a[2] = zeroes / (2 * fabs(lastZero - firstZero));
382 a[3] = 0;
383 a[4] = 0;
384 a[5] = 0;
385 if (guessFlags & GUESS_CONSTANT_GIVEN)
386 a[0] = constantGuess;
387 if (guessFlags & GUESS_FACTOR_GIVEN)
388 a[1] = factorGuess;
389 if (guessFlags & GUESS_FREQ_GIVEN)
390 a[2] = freqGuess;
391 if (guessFlags & GUESS_PHASE_GIVEN)
392 a[3] = phaseGuess;
393 if (guessFlags & GUESS_SLOPE_GIVEN)
394 a[4] = slopeGuess;
395 if (guessFlags & GUESS_RATE_GIVEN)
396 a[5] = rateGuess;
397
398 alo[1] = a[1] / 2;
399 if (!(da[0] = a[0] * 0.1))
400 da[0] = 0.01;
401 if (!(da[1] = a[1] * 0.1))
402 da[1] = 0.01;
403 da[2] = a[2] * 0.25;
404 da[3] = 0.01;
405 if (addSlope)
406 da[4] = 0.01;
407 if (addExponential) {
408 if (a[5])
409 da[5] = a[5]*0.01;
410 else
411 da[5] = 0.001;
412 }
413
414 if (lockFreq) {
415 alo[2] = ahi[2] = a[2];
416 da[2] = 0;
417 disable[2] = 1;
418 }
419
420 simplexMin(&result, a, da, alo, ahi, disable, 6, -DBL_MAX, tolerance, fitFunction,
421 (verbosity > 0 ? report : NULL), nEvalMax, nPassMax, 12, 3, 1.0, simplexFlags);
422
423 for (i = result = 0; i < nData; i++)
424 result += sqr(residualData[i]);
425 rmsResidual = sqrt(result / nData);
426 if (verbosity > 1) {
427 fprintf(stderr, "RMS deviation: %.15e\n", rmsResidual);
428 fprintf(stderr, "(RMS deviation)/(largest value): %.15e\n", rmsResidual / MAX(fabs(yMin), fabs(yMax)));
429 }
430 if (verbosity > 0) {
431 if (addSlope) {
432 fprintf(stderr, "Coefficients of fit to the form y = a0 + a1*sin(2*PI*a2*x + a3) + a4*x, a = \n");
433 for (i = 0; i < 5; i++)
434 fprintf(stderr, "%.8e ", a[i]);
435 } else {
436 fprintf(stderr, "Coefficients of fit to the form y = a0 + a1*sin(2*PI*a2*x + a3), a = \n");
437 for (i = 0; i < 4; i++)
438 fprintf(stderr, "%.8e ", a[i]);
439 }
440 fprintf(stderr, "\n");
441 }
442
443 if (!SDDS_StartPage(&OutputTable, nData) ||
444 !SDDS_SetColumn(&OutputTable, SDDS_SET_BY_INDEX, xData, nData, xIndex) ||
445 !SDDS_SetColumn(&OutputTable, SDDS_SET_BY_INDEX, fitData, nData, fitIndex) ||
446 !SDDS_SetParameters(&OutputTable, SDDS_PASS_BY_VALUE | SDDS_SET_BY_NAME,
447 "sinefitConstant", a[0],
448 "sinefitFactor", a[1],
449 "sinefitFrequency", a[2],
450 "sinefitPhase", a[3],
451 "sinefitRmsResidual", rmsResidual,
452 NULL) ||
453 (addSlope &&
454 (!SDDS_SetParameters(&OutputTable, SDDS_PASS_BY_VALUE | SDDS_SET_BY_NAME,
455 "sinefitSlope", a[4], NULL))) ||
456 (addExponential &&
457 (!SDDS_SetParameters(&OutputTable, SDDS_PASS_BY_VALUE | SDDS_SET_BY_NAME,
458 "sinefitRate", a[5], NULL))) ||
459 (fullOutput && (!SDDS_SetColumn(&OutputTable, SDDS_SET_BY_INDEX, yData, nData, yIndex) ||
460 !SDDS_SetColumn(&OutputTable, SDDS_SET_BY_INDEX, residualData, nData, residualIndex))) ||
461 !SDDS_WritePage(&OutputTable))
462 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
463 }
464
465 if (SDDS_Terminate(&InputTable) != 1) {
466 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
467 return EXIT_FAILURE;
468 }
469 if (SDDS_Terminate(&OutputTable) != 1) {
470 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
471 return EXIT_FAILURE;
472 }
473
474 return EXIT_SUCCESS;
475}
476
477void setupOutputFile(SDDS_DATASET *OutputTable, int32_t *xIndex, int32_t *yIndex, int32_t *fitIndex,
478 int32_t *residualIndex, char *output, long fullOutput, SDDS_DATASET *InputTable,
479 char *xName, char *yName, short columnMajorOrder, short addSlope, short addExponential) {
480 char *name, *yUnits, *description, *xUnits, *inverse_xUnits;
481 int32_t typeValue = SDDS_DOUBLE;
482 static char *residualNamePart = "Residual";
483 static char *residualDescriptionPart = "Residual of sinusoidal fit to ";
484
485 if (!SDDS_InitializeOutput(OutputTable, SDDS_BINARY, 0, NULL, "sddssinefit output", output) ||
486 !SDDS_TransferColumnDefinition(OutputTable, InputTable, xName, NULL) ||
487 !SDDS_ChangeColumnInformation(OutputTable, "type", &typeValue, SDDS_BY_NAME, xName) ||
488 (*xIndex = SDDS_GetColumnIndex(OutputTable, xName)) < 0 ||
489 !SDDS_GetColumnInformation(InputTable, "units", &xUnits, SDDS_BY_NAME, xName) ||
490 !SDDS_GetColumnInformation(InputTable, "units", &yUnits, SDDS_BY_NAME, yName)) {
491 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
492 exit(EXIT_FAILURE);
493 }
494 if (columnMajorOrder != -1)
495 OutputTable->layout.data_mode.column_major = columnMajorOrder;
496 else
497 OutputTable->layout.data_mode.column_major = InputTable->layout.data_mode.column_major;
498
499 name = tmalloc(sizeof(*name) * (strlen(yName) + strlen(residualNamePart) + 1));
500 description = tmalloc(sizeof(*description) * (strlen(yName) + strlen(residualDescriptionPart) + 1));
501
502 if (fullOutput) {
503 if (!SDDS_TransferColumnDefinition(OutputTable, InputTable, yName, NULL) ||
504 !SDDS_ChangeColumnInformation(OutputTable, "type", &typeValue, SDDS_BY_NAME, yName) ||
505 (*yIndex = SDDS_GetColumnIndex(OutputTable, yName)) < 0)
506 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
507 sprintf(name, "%s%s", yName, residualNamePart);
508 sprintf(description, "%s%s", yName, residualDescriptionPart);
509 if ((*residualIndex = SDDS_DefineColumn(OutputTable, name, NULL, yUnits, description, NULL, SDDS_DOUBLE, 0)) < 0)
510 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
511 }
512
513 sprintf(name, "%sFit", yName);
514 sprintf(description, "Sinusoidal fit to %s", yName);
515 if ((*fitIndex = SDDS_DefineColumn(OutputTable, name, NULL, yUnits, description, NULL, SDDS_DOUBLE, 0)) < 0)
516 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
517
518 inverse_xUnits = makeInverseUnits(xUnits);
519
520 if (SDDS_DefineParameter(OutputTable, "sinefitConstant", NULL, yUnits, "Constant term from sinusoidal fit",
521 NULL, SDDS_DOUBLE, 0) < 0 ||
522 SDDS_DefineParameter(OutputTable, "sinefitFactor", NULL, yUnits, "Factor from sinusoidal fit",
523 NULL, SDDS_DOUBLE, 0) < 0 ||
524 SDDS_DefineParameter(OutputTable, "sinefitFrequency", NULL, inverse_xUnits, "Frequency from sinusoidal fit",
525 NULL, SDDS_DOUBLE, 0) < 0 ||
526 SDDS_DefineParameter(OutputTable, "sinefitPhase", NULL, xUnits, "Phase from sinusoidal fit",
527 NULL, SDDS_DOUBLE, 0) < 0 ||
528 SDDS_DefineParameter(OutputTable, "sinefitRmsResidual", NULL, yUnits, "RMS residual from sinusoidal fit",
529 NULL, SDDS_DOUBLE, 0) < 0)
530 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
531
532 if (addSlope) {
533 if (SDDS_DefineParameter(OutputTable, "sinefitSlope", NULL, yUnits, "Slope term added to sinusoidal fit",
534 NULL, SDDS_DOUBLE, 0) < 0)
535 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
536 }
537 if (addExponential) {
538 if (SDDS_DefineParameter(OutputTable, "sinefitRate", NULL, inverse_xUnits, "Exponential decay rate for sinusoidal fit",
539 NULL, SDDS_DOUBLE, 0) < 0)
540 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
541 }
542 if (!SDDS_WriteLayout(OutputTable))
543 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
544
545 free(name);
546 free(description);
547}
548
549char *makeInverseUnits(char *units) {
550 char *inverseUnits;
551
552 if (!units || SDDS_StringIsBlank(units))
553 return NULL;
554 inverseUnits = tmalloc(sizeof(*inverseUnits) * (strlen(units) + 5));
555
556 if (strncmp(units, "1/(", 3) == 0 && units[strlen(units) - 1] == ')') {
557 /* Special case of "1/(<unit>)" */
558 strcpy(inverseUnits, units + 3);
559 inverseUnits[strlen(inverseUnits) - 1] = '\0';
560 } else if (!strchr(units, ' ')) {
561 /* Special case of units string without spaces */
562 sprintf(inverseUnits, "1/%s", units);
563 } else {
564 /* General case */
565 sprintf(inverseUnits, "1/(%s)", units);
566 }
567
568 return inverseUnits;
569}
570
571double fitFunction(double *a, long *invalid) {
572 int64_t i;
573 double chi, min_chi, fitValue;
574
575 min_chi = DBL_MAX;
576
577 *invalid = 0;
578 if (verbosity>5) {
579 fprintf(stderr, "Trial: a = %e, %e, %e, %e, %e, %e\n", a[0], a[1], a[2], a[3], a[4], a[5]);
580 }
581 for (i = chi = 0; i < nData; i++) {
582 fitValue = a[1] * sin(PIx2 * a[2] * xData[i] + a[3]);
583 if (addExponential)
584 fitValue *= exp(xData[i]*a[5]);
585 if (addSlope)
586 fitValue += a[4] * xData[i];
587 fitValue += a[0];
588 fitData[i] = fitValue;
589 residualData[i] = yData[i] - fitValue;
590 chi += sqr(residualData[i]);
591 }
592 if (isnan(chi) || isinf(chi))
593 *invalid = 1;
594 if (verbosity > 5)
595 fprintf(stderr, " --> chi = %e, invalid = %ld\n", chi, *invalid);
596 if (min_chi > chi) {
597 min_chi = chi;
598 fit[0] = a[0];
599 fit[1] = a[1];
600 fit[2] = a[2];
601 fit[3] = a[3];
602 fit[4] = a[4];
603 fit[5] = a[5];
604 if (verbosity > 2)
605 fprintf(stderr, "New best chi = %e: a = %e, %e, %e, %e, %e, %e\n",
606 chi, fit[0], fit[1], fit[2], fit[3], fit[4], fit[5]);
607 }
608 return chi;
609}
610
611void report(double y, double *x, long pass, long nEval, long n_dimen) {
612 long i;
613
614 fprintf(stderr, "Pass %ld, after %ld evaluations: result = %.16e\n", pass, nEval, y);
615 fprintf(stderr, "a = ");
616 for (i = 0; i < 6; i++)
617 fprintf(stderr, "%.8e ", x[i]);
618 fputc('\n', stderr);
619}
SDDS.h
SDDS (Self Describing Data Set) Data Types Definitions and Function Prototypes.
SDDS_StartPage
int32_t SDDS_StartPage(SDDS_DATASET *SDDS_dataset, int64_t expected_n_rows)
Definition SDDS_dataprep.c:70
SDDS_SetParameters
int32_t SDDS_SetParameters(SDDS_DATASET *SDDS_dataset, int32_t mode,...)
Definition SDDS_dataprep.c:375
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:1542
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_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_InitializeInput
int32_t SDDS_InitializeInput(SDDS_DATASET *SDDS_dataset, char *filename)
Definition SDDS_input.c:50
SDDS_Terminate
int32_t SDDS_Terminate(SDDS_DATASET *SDDS_dataset)
Definition SDDS_input.c:1471
SDDS_ReadPage
int32_t SDDS_ReadPage(SDDS_DATASET *SDDS_dataset)
Definition SDDS_input.c:1007
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:602
SDDS_WritePage
int32_t SDDS_WritePage(SDDS_DATASET *SDDS_dataset)
Writes the current data table to the output file.
Definition SDDS_output.c:1229
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:1716
SDDS_WriteLayout
int32_t SDDS_WriteLayout(SDDS_DATASET *SDDS_dataset)
Writes the SDDS layout header to the output file.
Definition SDDS_output.c:900
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:1473
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_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:1375
SDDS_PrintErrors
void SDDS_PrintErrors(FILE *fp, int32_t mode)
Prints recorded error messages to a specified file stream.
Definition SDDS_utils.c:474
SDDS_RegisterProgramName
void SDDS_RegisterProgramName(const char *name)
Registers the executable program name for use in error messages.
Definition SDDS_utils.c:318
SDDS_StringIsBlank
int32_t SDDS_StringIsBlank(char *s)
Checks if a string is blank (contains only whitespace characters).
Definition SDDS_utils.c:2470
SDDS_Bomb
void SDDS_Bomb(char *message)
Terminates the program after printing an error message and recorded errors.
Definition SDDS_utils.c:380
SDDS_Realloc
void * SDDS_Realloc(void *old_ptr, size_t new_size)
Reallocates memory to a new size.
Definition SDDS_utils.c:743
SDDS_LONG
#define SDDS_LONG
Identifier for the signed 32-bit integer data type.
Definition SDDStypes.h:61
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:65
bomb
void bomb(char *error, char *usage)
Reports error messages to the terminal and aborts the program.
Definition bomb.c:26
find_min_max
int find_min_max(double *min, double *max, double *list, int64_t n)
Finds the minimum and maximum values in a list of doubles.
Definition findMinMax.c:33
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:357
processFilenames
void processFilenames(char *programName, char **input, char **output, unsigned long pipeFlags, long noWarnings, long *tmpOutputUsed)
Definition scanargs.c:391
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
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:502
SDDS_DATASET
Definition SDDS.h:308
Generated by doxygen 1.12.0