SDDS ToolKit Programs and Libraries for C and Python
•All Classes Files Functions Variables Macros Pages
sddsmpfit.c
Go to the documentation of this file.
1/**
2 * @file sddsmpfit.c
3 * @brief Performs polynomial least-squares fitting on SDDS files.
4 *
5 * @details
6 * This program reads SDDS (Self Describing Data Sets) files, fits the data using either ordinary or Chebyshev T polynomials,
7 * and outputs the results including fitted data, residuals, and fit parameters with their uncertainties.
8 * Additional features include normalization, sigma modification, and fit evaluation.
9 *
10 * The fitting model is:
11 * \f[
12 * y = \sum_{i=0}^{N-1} A[i] \cdot P(x - x_{offset}, i)
13 * \f]
14 * where \‍( P(x, i) \‍) is the ith basis function, typically \‍( x^i \‍), or a Chebyshev polynomial.
15 *
16 * @section Usage
17 * ```
18 * sddsmpfit [<inputfile>] [<outputfile>]
19 * [-pipe=[input][,output]]
20 * -independent=<xName>
21 * -dependent=<yname1-wildcard>[,<yname2-wildcard>...]
22 * [-sigmaIndependent=<xSigma>]
23 * [-sigmaDependent=<ySigmaFormatString>]
24 * {
25 * -terms=<number> [-symmetry={none|odd|even}] |
26 * -orders=<number>[,<number>...]
27 * }
28 * [-reviseOrders[=threshold=<value>][,verbose]]
29 * [-chebyshev[=convert]]
30 * [-xOffset=<value>]
31 * [-xFactor=<value>]
32 * [-sigmas=<value>,{absolute|fractional}]
33 * [-minimumSigma=<value>]
34 * [-modifySigmas]
35 * [-generateSigmas={keepLargest|keepSmallest}]
36 * [-sparse=<interval>]
37 * [-range=<lower>,<upper>[,fitOnly]]
38 * [-normalize[=<termNumber>]]
39 * [-verbose]
40 * [-evaluate=<filename>[,begin=<value>][,end=<value>][,number=<integer>]]
41 * [-fitLabelFormat=<sprintf-string>]
42 * [-infoFile=<filename>]
43 * [-copyParameters]
44 * ```
45 *
46 * @section Options
47 * | Required | Description |
48 * |---------------------------------------|---------------------------------------------------------|
49 * | `-independent` | Specify the independent data column name. |
50 * | `-dependent` | Specify dependent data columns, using wildcards if needed. |
51 *
52 * | Optional | Description |
53 * |-------------------------------------|--------------------------------------------------------------------------------------|
54 * | `-pipe` | Use standard input and/or output. |
55 * | `-terms` | Number of terms to include in the fit. |
56 * | `-symmetry` | Symmetry of the fit about x_offset. |
57 * | `-orders` | Specify the polynomial orders to include in the fit. |
58 * | `-reviseOrders` | Revise the orders to eliminate poorly-determined coefficients. |
59 * | `-chebyshev` | Use Chebyshev T polynomials. Specify `convert` to convert back to ordinary polynomials.|
60 * | `-xOffset` | Set the x-offset for the fit. |
61 * | `-xFactor` | Set a scaling factor for x values. |
62 * | `-sigmas` | Set sigma values as absolute or fractional. |
63 * | `-minimumSigma` | Set a minimum sigma value; smaller values are replaced. |
64 * | `-modifySigmas` | Modify y-sigmas using x-sigmas and the initial fit. |
65 * | `-generateSigmas` | Generate y-sigmas based on RMS deviation, optionally keeping the largest/smallest sigmas.|
66 * | `-sparse` | Sample data at specified intervals. |
67 * | `-range` | Fit and evaluate within the specified range. |
68 * | `-normalize` | Normalize coefficients so that a specific term equals 1. |
69 * | `-evaluate` | Evaluate the fit over a specified range. |
70 * | `-fitLabelFormat` | Specify a format string for fit labels. |
71 * | `-infoFile` | Output fit coefficients and statistics to the specified information file. |
72 * | `-copyParameters` | Copy parameters from the input to the output SDDS file. |
73 * | `-verbose` | Enable verbose output for debugging or additional information. |
74 *
75 * @subsection Incompatibilities
76 * - `-terms` is incompatible with `-orders`.
77 * - `-generateSigmas` is incompatible with `-modifySigmas`.
78 * - `-sigmas` cannot be used with `-sigmaDependent`.
79 * - `-reviseOrders` requires `-generateSigmas`, `-sigmas`, or `-sigmaDependent`.
80 *
81 * @copyright
82 * - (c) 2002 The University of Chicago, as Operator of Argonne National Laboratory.
83 * - (c) 2002 The Regents of the University of California, as Operator of Los Alamos National Laboratory.
84 *
85 * @license
86 * This file is distributed under the terms of the Software License Agreement
87 * found in the file LICENSE included with this distribution.
88 *
89 * @author
90 * M. Borland, Brad Dolin, R. Soliday, H. Shang
91 */
92
93#include "mdb.h"
94#include "SDDS.h"
95#include "scan.h"
96
97void print_coefs(FILE *fprec, double x_offset, double x_scale, long chebyshev, double *coef, double *s_coef,
98 int32_t *order, long n_terms, double chi, long norm_term, char *prepend);
99char **makeCoefficientUnits(SDDS_DATASET *SDDSout, char *xName, char *yName, int32_t *order, long terms);
100long setCoefficientData(SDDS_DATASET *SDDSout, double *coef, double *coefSigma, char **coefUnits, long *order, long terms);
101char ***initializeOutputFile(SDDS_DATASET *SDDSout, SDDS_DATASET *SDDSoutInfo, char *output, char *outputInfo,
102 SDDS_DATASET *SDDSin, char *input, char *xName, char **yNames, char *xSigmaName,
103 char **ySigmaNames, long sigmasValid, int32_t *order, long terms, long isChebyshev,
104 long numCols, long copyParameters);
105void checkInputFile(SDDS_DATASET *SDDSin, char *xName, char **yNames, char *xSigmaName, char **ySigmaNames, long numYNames);
106long coefficient_index(int32_t *order, long terms, long order_of_interest);
107void makeFitLabel(char *buffer, long bufsize, char *fitLabelFormat, double *coef, int32_t *order, long terms, long colIndex);
108
109char **ResolveColumnNames(SDDS_DATASET *SDDSin, char **wildcardList, long length, int32_t *numYNames);
110char **GenerateYSigmaNames(char *controlString, char **yNames, long numYNames);
111void RemoveElementFromStringArray(char **array, long index, long length);
112char **RemoveNonNumericColumnsFromNameArray(SDDS_DATASET *SDDSin, char **columns, int32_t *numColumns);
113void compareOriginalToFit(double *x, double *y, double **residual, int64_t points, double *rmsResidual, double *coef, int32_t *order, long terms);
114
115void set_argument_offset(double offset);
116void set_argument_scale(double scale);
117double tcheby(double x, long n);
118double dtcheby(double x, long n);
119double ipower(double x, long n);
120double dipower(double x, long n);
121
122/* Enumeration for option types */
123enum option_type {
124 CLO_DEPENDENT,
125 CLO_ORDERS,
126 CLO_TERMS,
127 CLO_SYMMETRY,
128 CLO_REVISEORDERS,
129 CLO_CHEBYSHEV,
130 CLO_MODIFYSIGMAS,
131 CLO_SIGMAS,
132 CLO_GENERATESIGMAS,
133 CLO_RANGE,
134 CLO_SPARSE,
135 CLO_NORMALIZE,
136 CLO_XFACTOR,
137 CLO_XOFFSET,
138 CLO_VERBOSE,
139 CLO_FITLABELFORMAT,
140 CLO_PIPE,
141 CLO_EVALUATE,
142 CLO_INDEPENDENT,
143 CLO_SIGMAINDEPENDENT,
144 CLO_SIGMADEPENDENT,
145 CLO_INFOFILE,
146 CLO_COPYPARAMETERS,
147 CLO_MINSIGMA,
148 N_OPTIONS
149};
150
151char *option[N_OPTIONS] = {
152 "dependent",
153 "orders",
154 "terms",
155 "symmetry",
156 "reviseorders",
157 "chebyshev",
158 "modifysigmas",
159 "sigmas",
160 "generatesigmas",
161 "range",
162 "sparse",
163 "normalize",
164 "xfactor",
165 "xoffset",
166 "verbose",
167 "fitlabelformat",
168 "pipe",
169 "evaluate",
170 "independent",
171 "sigmaindependent",
172 "sigmadependent",
173 "infofile",
174 "copyparameters",
175 "minimumsigma",
176};
177
178char *USAGE =
179 "sddsmpfit [<inputfile>] [<outputfile>]\n"
180 " [-pipe=[input][,output]]\n"
181 " -independent=<xName>\n"
182 " -dependent=<yname1-wildcard>[,<yname2-wildcard>...]\n"
183 " [-sigmaIndependent=<xSigma>]\n"
184 " [-sigmaDependent=<ySigmaFormatString>]\n"
185 " {\n"
186 " -terms=<number> [-symmetry={none|odd|even}] | \n"
187 " -orders=<number>[,<number>...] \n"
188 " }\n"
189 " [-reviseOrders[=threshold=<value>][,verbose]]\n"
190 " [-chebyshev[=convert]]\n"
191 " [-xOffset=<value>] \n"
192 " [-xFactor=<value>]\n"
193 " [-sigmas=<value>,{absolute|fractional}] \n"
194 " [-minimumSigma=<value>]\n"
195 " [-modifySigmas] \n"
196 " [-generateSigmas={keepLargest|keepSmallest}]\n"
197 " [-sparse=<interval>] \n"
198 " [-range=<lower>,<upper>[,fitOnly]]\n"
199 " [-normalize[=<termNumber>]] \n"
200 " [-verbose]\n"
201 " [-evaluate=<filename>[,begin=<value>][,end=<value>][,number=<integer>]]\n"
202 " [-fitLabelFormat=<sprintf-string>] \n"
203 " [-infoFile=<filename>]\n"
204 " [-copyParameters]\n"
205 "Program by Michael Borland, revised by Brad Dolin.\n"
206 "Compiled on " __DATE__ " " __TIME__ ", SVN revision: " SVN_VERSION "\n";
207
208static char *additional_help = "\n\
209sddsmpfit does fits to the form y = SUM(i){ A[i] *P(x-x_offset, i)}, where P(x,i) is the ith basis\n\
210function evaluated at x. sddsmpfit returns the A[i] and estimates of the errors in the values.\n\
211By default P(x,i) = x^i. One can also select Chebyshev T polynomials as the basis functions.\n\n\
212-independent specify name of independent data column to use.\n\
213-dependent specify names of dependent data columns to use, using wildcards,\n\
214 separated by commas.\n\
215-sigmaIndependent specify name of independent sigma values to use\n\
216-sigmaDependent specify names of dependent sigma values to use by specifying a printf-style control\n\
217 string to generate the names from the independent variable names (e.g., %sSigma)\n\
218-terms number of terms desired in fit.\n\
219-symmetry symmetry of desired fit about x_offset.\n\
220-orders orders (P[i]) to use in fitting.\n\
221-reviseOrders the orders used in the fit are modified from the specified ones\n\
222 in order eliminate poorly-determined coefficients, based on fitting\n\
223 of the first data page.\n";
224static char *additional_help2 = "-chebyshev use Chebyshev T polynomials (xOffset is set automatically).\n\
225 Giving the `convert' option causes the fit to be written out in\n\
226 terms of ordinary polynomials.\n\
227-xOffset desired value of x to fit about.\n\
228-xFactor desired factor to multiply x values by before fitting.\n\
229-sigmas specify absolute or fractional sigma for all points.\n\
230-minimumSigma specify minimum sigma value. If the value is less than this\n\
231 it is replaced by this value.\n\
232-modifySigmas modify the y sigmas using the x sigmas and an initial fit.\n\
233-generateSigmas generate y sigmas from the rms deviation from an initial fit.\n\
234 optionally keep the sigmas from the data if larger/smaller than rms\n\
235 deviation.\n\
236-sparse specify integer interval at which to sample data.\n\
237-range specify range of independent variable over which to perform fit and evaluation.\n\
238 If 'fitOnly' is given, then fit is compared to data over the original range.\n\
239-normalize normalize so that specified term is unity.\n\
240-evaluate specify evaluation of fit over a selected range of\n\
241 equispaced points.\n\
242-infoFile specify name of optional information file containing coefficients and fit statistics.\n\
243-copyParameters specify that parameters from input should be copied to output.\n\
244-verbose generates extra output that may be useful.\n\n";
245
246#define NO_SYMMETRY 0
247#define EVEN_SYMMETRY 1
248#define ODD_SYMMETRY 2
249#define N_SYMMETRY_OPTIONS 3
250char *symmetry_options[N_SYMMETRY_OPTIONS] = {"none", "even", "odd"};
251
252#define ABSOLUTE_SIGMAS 0
253#define FRACTIONAL_SIGMAS 1
254#define N_SIGMAS_OPTIONS 2
255char *sigmas_options[N_SIGMAS_OPTIONS] = {"absolute", "fractional"};
256
257#define FLGS_GENERATESIGMAS 1
258#define FLGS_KEEPLARGEST 2
259#define FLGS_KEEPSMALLEST 4
260
261#define REVPOW_ACTIVE 0x0001
262#define REVPOW_VERBOSE 0x0002
263/* SDDS indices for output page */
264static long *iIntercept = NULL, *iInterceptO = NULL, *iInterceptSigma = NULL, *iInterceptSigmaO = NULL;
265static long *iSlope = NULL, *iSlopeO = NULL, *iSlopeSigma = NULL, *iSlopeSigmaO = NULL;
266static long *iCurvature = NULL, *iCurvatureO = NULL, *iCurvatureSigma = NULL, *iCurvatureSigmaO = NULL;
267static long iOffset = -1, iOffsetO = -1, iFactor = -1, iFactorO = -1;
268static long *iChiSq = NULL, *iChiSqO = NULL, *iRmsResidual = NULL, *iRmsResidualO = NULL, *iSigLevel = NULL, *iSigLevelO = NULL;
269static long *iFitIsValid = NULL, *iFitIsValidO = NULL, *iFitLabel = NULL, *iFitLabelO = NULL, iTerms = -1, iTermsO = -1;
270
271static long ix = -1, ixSigma = -1;
272static long *iy = NULL, *iySigma = NULL;
273static long *iFit = NULL, *iResidual = NULL;
274
275static long iOrder = -1, *iCoefficient = NULL, *iCoefficientSigma = NULL, *iCoefficientUnits = NULL;
276
277static char *xSymbol, **ySymbols;
278
279#define EVAL_BEGIN_GIVEN 0x0001U
280#define EVAL_END_GIVEN 0x0002U
281#define EVAL_NUMBER_GIVEN 0x0004U
282
283#define MAX_Y_SIGMA_NAME_SIZE 1024
284
285typedef struct
286{
287 char *file;
288 long initialized;
289 int64_t number;
290 unsigned long flags;
291 double begin, end;
292 SDDS_DATASET dataset;
293} EVAL_PARAMETERS;
294void setupEvaluationFile(EVAL_PARAMETERS *evalParameters, char *xName, char **yName, long yNames, SDDS_DATASET *SDDSin);
295void makeEvaluationTable(EVAL_PARAMETERS *evalParameters, double *x, int64_t points, double *coef, int32_t *order,
296 long terms, char *xName, char **yName, long yNames, long iYName);
297
298static double (*basis_fn)(double xa, long ordera);
299static double (*basis_dfn)(double xa, long ordera);
300
301int main(int argc, char **argv) {
302 double **y = NULL, **sy = NULL, **diff = NULL;
303 double *x = NULL, *sx = NULL;
304 double xOffset, xScaleFactor;
305 double *xOrig = NULL, **yOrig = NULL, *sxOrig = NULL, **syOrig = NULL, **sy0 = NULL;
306 long terms, normTerm, ip, ySigmasValid;
307 int64_t i, j, points, pointsOrig;
308 long symmetry, chebyshev, termsGiven;
309 double sigmas, minimumSigma;
310 long sigmasMode, sparseInterval;
311 double **coef = NULL, **coefSigma = NULL;
312 double *chi = NULL, xLow, xHigh, *rmsResidual = NULL;
313 char *xName = NULL, *yName = NULL, **yNames = NULL, *xSigmaName = NULL;
314 char **ySigmaNames = NULL, *ySigmaControlString = NULL;
315 char *input = NULL, *output = NULL;
316 SDDS_DATASET SDDSin, SDDSout, SDDSoutInfo;
317 long *isFit = NULL, iArg, modifySigmas, termIndex;
318 long generateSigmas, verbose, ignoreSigmas;
319 long outputInitialized, copyParameters = 0;
320 int32_t *order = NULL;
321 SCANNED_ARG *s_arg;
322 double xMin, xMax, revpowThreshold;
323 double rms_average(double *d_x, int64_t d_n);
324 char *infoFile = NULL;
325 char *fitLabelFormat = "%g";
326 static char fitLabelBuffer[SDDS_MAXLINE];
327 unsigned long pipeFlags, reviseOrders;
328 EVAL_PARAMETERS evalParameters;
329 long rangeFitOnly = 0;
330
331 long colIndex;
332 long cloDependentIndex = -1, numDependentItems;
333 int32_t numYNames;
334
335 SDDS_RegisterProgramName(argv[0]);
336 argc = scanargs(&s_arg, argc, argv);
337 if (argc < 2 || argc > (3 + N_OPTIONS)) {
338 fprintf(stderr, "usage: %s\n", USAGE);
339 fprintf(stderr, "%s%s", additional_help, additional_help2);
340 exit(EXIT_FAILURE);
341 }
342
343 input = output = NULL;
344 xName = yName = xSigmaName = ySigmaControlString = NULL;
345 yNames = ySigmaNames = NULL;
346 numDependentItems = 0;
347 modifySigmas = reviseOrders = chebyshev = 0;
348 order = NULL;
349 symmetry = NO_SYMMETRY;
350 xMin = xMax = 0;
351 generateSigmas = 0;
352 sigmasMode = -1;
353 sigmas = 1;
354 minimumSigma = 0;
355 sparseInterval = 1;
356 terms = 2;
357 verbose = ignoreSigmas = 0;
358 normTerm = -1;
359 xOffset = 0;
360 xScaleFactor = 1;
361 basis_fn = ipower;
362 basis_dfn = dipower;
363 pipeFlags = 0;
364 evalParameters.file = NULL;
365 infoFile = NULL;
366 termsGiven = 0;
367
368 for (iArg = 1; iArg < argc; iArg++) {
369 if (s_arg[iArg].arg_type == OPTION) {
370 switch (match_string(s_arg[iArg].list[0], option, N_OPTIONS, 0)) {
371 case CLO_MODIFYSIGMAS:
372 modifySigmas = 1;
373 break;
374 case CLO_ORDERS:
375 if (termsGiven)
376 SDDS_Bomb("give -order or -terms, not both");
377 if (s_arg[iArg].n_items < 2)
378 SDDS_Bomb("invalid -orders syntax");
379 order = tmalloc(sizeof(*order) * (terms = s_arg[iArg].n_items - 1));
380 for (i = 1; i < s_arg[iArg].n_items; i++) {
381 if (sscanf(s_arg[iArg].list[i], "%" SCNd32, order + i - 1) != 1)
382 SDDS_Bomb("unable to scan order from -orders list");
383 }
384 break;
385 case CLO_RANGE:
386 rangeFitOnly = 0;
387 if ((s_arg[iArg].n_items != 3 && s_arg[iArg].n_items != 4) || 1 != sscanf(s_arg[iArg].list[1], "%lf", &xMin) || 1 != sscanf(s_arg[iArg].list[2], "%lf", &xMax) || xMin >= xMax)
388 SDDS_Bomb("incorrect -range syntax");
389 if (s_arg[iArg].n_items == 4) {
390 if (strncmp(str_tolower(s_arg[iArg].list[3]), "fitonly", strlen(s_arg[iArg].list[3])) == 0) {
391 rangeFitOnly = 1;
392 } else
393 SDDS_Bomb("incorrect -range syntax");
394 }
395 break;
396 case CLO_GENERATESIGMAS:
397 generateSigmas = FLGS_GENERATESIGMAS;
398 if (s_arg[iArg].n_items > 1) {
399 if (s_arg[iArg].n_items != 2)
400 SDDS_Bomb("incorrect -generateSigmas synax");
401 if (strncmp(s_arg[iArg].list[1], "keepsmallest", strlen(s_arg[iArg].list[1])) == 0)
402 generateSigmas |= FLGS_KEEPSMALLEST;
403 if (strncmp(s_arg[iArg].list[1], "keeplargest", strlen(s_arg[iArg].list[1])) == 0)
404 generateSigmas |= FLGS_KEEPLARGEST;
405 if ((generateSigmas & FLGS_KEEPSMALLEST) && (generateSigmas & FLGS_KEEPLARGEST))
406 SDDS_Bomb("ambiguous -generateSigmas synax");
407 }
408 break;
409 case CLO_TERMS:
410 if (order)
411 SDDS_Bomb("give -order or -terms, not both");
412 if (s_arg[iArg].n_items != 2 || sscanf(s_arg[iArg].list[1], "%ld", &terms) != 1)
413 SDDS_Bomb("invalid -terms syntax");
414 termsGiven = 1;
415 break;
416 case CLO_XOFFSET:
417 if (s_arg[iArg].n_items != 2 || sscanf(s_arg[iArg].list[1], "%lf", &xOffset) != 1)
418 SDDS_Bomb("invalid -xOffset syntax");
419 break;
420 case CLO_SYMMETRY:
421 if (s_arg[iArg].n_items == 2) {
422 if ((symmetry = match_string(s_arg[iArg].list[1], symmetry_options, N_SYMMETRY_OPTIONS, 0)) < 0)
423 SDDS_Bomb("unknown option used with -symmetry");
424 } else
425 SDDS_Bomb("incorrect -symmetry syntax");
426 break;
427 case CLO_SIGMAS:
428 if (s_arg[iArg].n_items != 3)
429 SDDS_Bomb("incorrect -sigmas syntax");
430 if (sscanf(s_arg[iArg].list[1], "%lf", &sigmas) != 1)
431 SDDS_Bomb("couldn't scan value for -sigmas");
432 if ((sigmasMode = match_string(s_arg[iArg].list[2], sigmas_options, N_SIGMAS_OPTIONS, 0)) < 0)
433 SDDS_Bomb("unrecognized -sigmas mode");
434 break;
435 case CLO_MINSIGMA:
436 if (s_arg[iArg].n_items != 2)
437 SDDS_Bomb("incorrect -minimumSigma syntax");
438 if (sscanf(s_arg[iArg].list[1], "%lf", &minimumSigma) != 1)
439 SDDS_Bomb("couldn't scan value for -minimumSigma");
440 break;
441 case CLO_SPARSE:
442 if (s_arg[iArg].n_items != 2)
443 SDDS_Bomb("incorrect -sparse syntax");
444 if (sscanf(s_arg[iArg].list[1], "%ld", &sparseInterval) != 1)
445 SDDS_Bomb("couldn't scan value for -sparse");
446 if (sparseInterval < 1)
447 SDDS_Bomb("invalid -sparse value");
448 break;
449 case CLO_VERBOSE:
450 verbose = 1;
451 break;
452 case CLO_NORMALIZE:
453 normTerm = 0;
454 if (s_arg[iArg].n_items > 2 ||
455 (s_arg[iArg].n_items == 2 && sscanf(s_arg[iArg].list[1], "%ld", &normTerm) != 1) ||
456 normTerm < 0)
457 SDDS_Bomb("invalid -normalize syntax");
458 break;
459 case CLO_REVISEORDERS:
460 revpowThreshold = 0.1;
461 s_arg[iArg].n_items -= 1;
462 if (!scanItemList(&reviseOrders, s_arg[iArg].list + 1, &s_arg[iArg].n_items, 0,
463 "threshold", SDDS_DOUBLE, &revpowThreshold, 1, 0,
464 "verbose", -1, NULL, 1, REVPOW_VERBOSE, NULL))
465 SDDS_Bomb("invalid -reviseOrders syntax");
466 reviseOrders |= REVPOW_ACTIVE;
467 revpowThreshold = fabs(revpowThreshold);
468 break;
469 case CLO_CHEBYSHEV:
470 if (s_arg[iArg].n_items > 2 ||
471 (s_arg[iArg].n_items == 2 && strncmp(s_arg[iArg].list[1], "convert", strlen(s_arg[iArg].list[1])) != 0))
472 SDDS_Bomb("invalid -chebyshev syntax");
473 chebyshev = s_arg[iArg].n_items;
474 basis_fn = tcheby;
475 basis_dfn = dtcheby;
476 break;
477 case CLO_XFACTOR:
478 if (s_arg[iArg].n_items != 2 ||
479 sscanf(s_arg[iArg].list[1], "%lf", &xScaleFactor) != 1 || xScaleFactor == 0)
480 SDDS_Bomb("invalid -xFactor syntax");
481 break;
482 case CLO_INDEPENDENT:
483 if (s_arg[iArg].n_items != 2)
484 SDDS_Bomb("invalid -independent syntax");
485 xName = s_arg[iArg].list[1];
486 break;
487 case CLO_DEPENDENT:
488 numDependentItems = s_arg[iArg].n_items - 1;
489 cloDependentIndex = iArg;
490 if (numDependentItems < 1)
491 SDDS_Bomb("invalid -dependent syntax");
492 break;
493 case CLO_SIGMAINDEPENDENT:
494 if (s_arg[iArg].n_items != 2)
495 SDDS_Bomb("invalid -sigmaIndependent syntax");
496 xSigmaName = s_arg[iArg].list[1];
497 break;
498 case CLO_SIGMADEPENDENT:
499 if (s_arg[iArg].n_items != 2)
500 SDDS_Bomb("invalid -sigmaDependent syntax");
501 ySigmaControlString = s_arg[iArg].list[1];
502 break;
503 case CLO_FITLABELFORMAT:
504 if (s_arg[iArg].n_items != 2)
505 SDDS_Bomb("invalid -fitLabelFormat syntax");
506 fitLabelFormat = s_arg[iArg].list[1];
507 break;
508 case CLO_PIPE:
509 if (!processPipeOption(s_arg[iArg].list + 1, s_arg[iArg].n_items - 1, &pipeFlags))
510 SDDS_Bomb("invalid -pipe syntax");
511 break;
512 case CLO_INFOFILE:
513 if (s_arg[iArg].n_items != 2)
514 SDDS_Bomb("invalid -infoFile syntax");
515 infoFile = s_arg[iArg].list[1];
516 break;
517 case CLO_EVALUATE:
518 if (s_arg[iArg].n_items < 2)
519 SDDS_Bomb("invalid -evaluate syntax");
520 evalParameters.file = s_arg[iArg].list[1];
521 s_arg[iArg].n_items -= 2;
522 s_arg[iArg].list += 2;
523 if (!scanItemList(&evalParameters.flags, s_arg[iArg].list, &s_arg[iArg].n_items, 0,
524 "begin", SDDS_DOUBLE, &evalParameters.begin, 1, EVAL_BEGIN_GIVEN,
525 "end", SDDS_DOUBLE, &evalParameters.end, 1, EVAL_END_GIVEN,
526 "number", SDDS_LONG64, &evalParameters.number, 1, EVAL_NUMBER_GIVEN, NULL))
527 SDDS_Bomb("invalid -evaluate syntax");
528 break;
529 case CLO_COPYPARAMETERS:
530 copyParameters = 1;
531 break;
532 default:
533 bomb("unknown switch", USAGE);
534 break;
535 }
536 } else {
537 if (input == NULL)
538 input = s_arg[iArg].list[0];
539 else if (output == NULL)
540 output = s_arg[iArg].list[0];
541 else
542 SDDS_Bomb("too many filenames");
543 }
544 }
545
546 processFilenames("sddsmpfit", &input, &output, pipeFlags, 0, NULL);
547
548 if (symmetry && order)
549 SDDS_Bomb("can't specify both -symmetry and -orders");
550 if (!xName || !numDependentItems)
551 SDDS_Bomb("you must specify a column name for x and y");
552 if (modifySigmas && !xSigmaName)
553 SDDS_Bomb("you must specify x sigmas with -modifySigmas");
554 if (generateSigmas) {
555 if (modifySigmas)
556 SDDS_Bomb("you can't specify both -generateSigmas and -modifySigmas");
557 }
558 if (ySigmaControlString) {
559 if (sigmasMode != -1)
560 SDDS_Bomb("you can't specify both -sigmas and a y sigma name");
561 }
562 ySigmasValid = 0;
563 if (sigmasMode != -1 || generateSigmas || ySigmaControlString || modifySigmas)
564 ySigmasValid = 1;
565
566 if (normTerm >= 0 && normTerm >= terms)
567 SDDS_Bomb("can't normalize to that term--not that many terms");
568 if (reviseOrders && !(sigmasMode != -1 || generateSigmas || ySigmaNames))
569 SDDS_Bomb("can't use -reviseOrders unless a y sigma or -generateSigmas is given");
570
571 if (symmetry == EVEN_SYMMETRY) {
572 order = tmalloc(sizeof(*order) * terms);
573 for (i = 0; i < terms; i++)
574 order[i] = 2 * i;
575 } else if (symmetry == ODD_SYMMETRY) {
576 order = tmalloc(sizeof(*order) * terms);
577 for (i = 0; i < terms; i++)
578 order[i] = 2 * i + 1;
579 } else if (!order) {
580 order = tmalloc(sizeof(*order) * terms);
581 for (i = 0; i < terms; i++)
582 order[i] = i;
583 }
584
585 if (!SDDS_InitializeInput(&SDDSin, input))
586 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
587 outputInitialized = 0;
588 yNames = ResolveColumnNames(&SDDSin, s_arg[cloDependentIndex].list + 1, numDependentItems, &numYNames);
589 if (ySigmaControlString != NULL)
590 ySigmaNames = GenerateYSigmaNames(ySigmaControlString, yNames, numYNames);
591
592 checkInputFile(&SDDSin, xName, yNames, xSigmaName, ySigmaNames, numYNames);
593 sy0 = tmalloc(sizeof(double *) * numYNames);
594 y = tmalloc(sizeof(double *) * numYNames);
595 sy = tmalloc(sizeof(double *) * numYNames);
596 isFit = tmalloc(sizeof(long) * numYNames);
597 chi = tmalloc(sizeof(double) * numYNames);
598 coef = tmalloc(sizeof(double *) * numYNames);
599 coefSigma = tmalloc(sizeof(double *) * numYNames);
600 for (colIndex = 0; colIndex < numYNames; colIndex++) {
601 coef[colIndex] = tmalloc(sizeof(double) * terms);
602 coefSigma[colIndex] = tmalloc(sizeof(double) * terms);
603 }
604 iCoefficient = tmalloc(sizeof(long) * numYNames);
605 iCoefficientSigma = tmalloc(sizeof(long) * numYNames);
606 iCoefficientUnits = tmalloc(sizeof(long) * numYNames);
607
608 while (SDDS_ReadPage(&SDDSin) > 0) {
609 if ((points = SDDS_CountRowsOfInterest(&SDDSin)) < terms) {
610 /* probably should emit an empty page here */
611 continue;
612 }
613 if (!(x = SDDS_GetColumnInDoubles(&SDDSin, xName))) {
614 fprintf(stderr, "error: unable to read column %s\n", xName);
615 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
616 }
617 for (i = 0; i < numYNames; i++) {
618 if (!(y[i] = SDDS_GetColumnInDoubles(&SDDSin, yNames[i]))) {
619 fprintf(stderr, "error: unable to read column %s\n", yNames[i]);
620 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
621 }
622 }
623 sx = NULL;
624 if (xSigmaName && !(sx = SDDS_GetColumnInDoubles(&SDDSin, xSigmaName))) {
625 fprintf(stderr, "error: unable to read column %s\n", xSigmaName);
626 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
627 }
628 for (colIndex = 0; colIndex < numYNames; colIndex++)
629 sy0[colIndex] = tmalloc(sizeof(double) * points);
630 if (ySigmaNames) {
631 for (i = 0; i < numYNames; i++) {
632 if (!(sy0[i] = SDDS_GetColumnInDoubles(&SDDSin, ySigmaNames[i]))) {
633 fprintf(stderr, "error: unable to read column %s\n", ySigmaNames[i]);
634 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
635 }
636 }
637 }
638
639 if (minimumSigma > 0) {
640 int64_t j;
641 for (i = 0; i < numYNames; i++) {
642 for (j = 0; j < points; j++)
643 if (sy0[i][j] < minimumSigma)
644 sy0[i][j] = minimumSigma;
645 }
646 }
647
648 if (xMin != xMax || sparseInterval != 1) {
649 xOrig = tmalloc(sizeof(*xOrig) * points);
650 yOrig = tmalloc(sizeof(*yOrig) * numYNames);
651 for (colIndex = 0; colIndex < numYNames; colIndex++)
652 yOrig[colIndex] = tmalloc(sizeof(double) * points);
653 if (sx)
654 sxOrig = tmalloc(sizeof(*sxOrig) * points);
655 if (ySigmasValid) {
656 syOrig = tmalloc(sizeof(*syOrig) * numYNames);
657 for (colIndex = 0; colIndex < numYNames; colIndex++)
658 syOrig[colIndex] = tmalloc(sizeof(double) * points);
659 }
660 pointsOrig = points;
661 for (i = j = 0; i < points; i++) {
662 xOrig[i] = x[i];
663 if (sx)
664 sxOrig[i] = sx[i];
665 for (colIndex = 0; colIndex < numYNames; colIndex++) {
666 yOrig[colIndex][i] = y[colIndex][i];
667 if (ySigmasValid)
668 syOrig[colIndex][i] = sy0[colIndex][i];
669 }
670 }
671 if (xMin != xMax) {
672 for (i = j = 0; i < points; i++) {
673 if (xOrig[i] <= xMax && xOrig[i] >= xMin) {
674 x[j] = xOrig[i];
675 for (colIndex = 0; colIndex < numYNames; colIndex++) {
676 y[colIndex][j] = yOrig[colIndex][i];
677 if (ySigmasValid)
678 sy0[colIndex][j] = syOrig[colIndex][i];
679 }
680 if (sx)
681 sx[j] = sxOrig[i];
682 j++;
683 }
684 }
685 points = j;
686 }
687 if (sparseInterval != 1) {
688 for (i = j = 0; i < points; i++) {
689 if (i % sparseInterval == 0) {
690 x[j] = x[i];
691 for (colIndex = 0; colIndex < numYNames; colIndex++) {
692 y[colIndex][j] = y[colIndex][i];
693 if (ySigmasValid)
694 sy0[colIndex][j] = sy0[colIndex][i];
695 }
696 if (sx)
697 sx[j] = sx[i];
698 j++;
699 }
700 }
701 points = j;
702 }
703 } else {
704 xOrig = x;
705 yOrig = y;
706 sxOrig = sx;
707 syOrig = sy0;
708 pointsOrig = points;
709 }
710
711 find_min_max(&xLow, &xHigh, x, points);
712
713 if (sigmasMode == ABSOLUTE_SIGMAS) {
714 for (colIndex = 0; colIndex < numYNames; colIndex++) {
715 for (i = 0; i < points; i++)
716 sy0[colIndex][i] = sigmas;
717 if (sy0[colIndex] != syOrig[colIndex])
718 for (i = 0; i < pointsOrig; i++)
719 syOrig[colIndex][i] = sigmas;
720 }
721 } else if (sigmasMode == FRACTIONAL_SIGMAS) {
722 for (colIndex = 0; colIndex < numYNames; colIndex++) {
723 for (i = 0; i < points; i++)
724 sy0[colIndex][i] = sigmas * fabs(y[colIndex][i]);
725 if (sy0[colIndex] != syOrig[colIndex])
726 for (i = 0; i < pointsOrig; i++)
727 syOrig[colIndex][i] = fabs(yOrig[colIndex][i]) * sigmas;
728 }
729 }
730
731 for (i = 0; i < numYNames; i++) {
732 if (minimumSigma > 0) {
733 int64_t j;
734 for (j = 0; j < points; j++)
735 if (sy0[i][j] < minimumSigma)
736 sy0[i][j] = minimumSigma;
737 }
738 }
739
740 if (!ySigmasValid || generateSigmas)
741 for (colIndex = 0; colIndex < numYNames; colIndex++) {
742 for (i = 0; i < points; i++)
743 sy0[colIndex][i] = 1;
744 }
745 else
746 for (i = 0; i < points; i++)
747 for (colIndex = 0; colIndex < numYNames; colIndex++) {
748 if (sy0[colIndex][i] == 0)
749 SDDS_Bomb("y sigma = 0 for one or more points.");
750 }
751
752 diff = tmalloc(sizeof(*diff) * numYNames);
753 sy = tmalloc(sizeof(*sy) * numYNames);
754 for (colIndex = 0; colIndex < numYNames; colIndex++) {
755 diff[colIndex] = tmalloc(sizeof(double) * points);
756 sy[colIndex] = tmalloc(sizeof(double) * points);
757 }
758
759 for (i = 0; i < points; i++) {
760 for (colIndex = 0; colIndex < numYNames; colIndex++)
761 sy[colIndex][i] = sy0[colIndex][i];
762 }
763
764 set_argument_offset(xOffset);
765 set_argument_scale(xScaleFactor);
766 if (chebyshev) {
767 xOffset = (xHigh + xLow) / 2;
768 set_argument_offset(xOffset);
769 set_argument_scale(xScaleFactor = (xHigh - xLow) / 2);
770 }
771
772 if (generateSigmas || modifySigmas) {
773 /* do an initial fit */
774 for (colIndex = 0; colIndex < numYNames; colIndex++) {
775 isFit[colIndex] = lsfg(x, y[colIndex], sy[colIndex], points, terms, order, coef[colIndex], coefSigma[colIndex], &chi[colIndex], diff[colIndex], basis_fn);
776 if (!isFit[colIndex]) {
777 fprintf(stderr, "Column %s: ", yNames[colIndex]);
778 SDDS_Bomb("initial fit failed.");
779 }
780 if (verbose) {
781 fprintf(stderr, "Column %s: ", yNames[colIndex]);
782 fputs("initial_fit:", stderr);
783 print_coefs(stderr, xOffset, xScaleFactor, chebyshev, coef[colIndex], NULL, order, terms, chi[colIndex], normTerm, "");
784 fprintf(stderr, "unweighted rms deviation from fit: %21.15le\n", rms_average(diff[colIndex], points));
785 }
786 if (modifySigmas) {
787 if (!ySigmasValid) {
788 for (i = 0; i < points; i++)
789 sy[colIndex][i] = fabs(eval_sum(basis_dfn, coef[colIndex], order, terms, x[i]) * sx[i]);
790 } else
791 for (i = 0; i < points; i++) {
792 sy[colIndex][i] = sqrt(sqr(sy0[colIndex][i]) + sqr(eval_sum(basis_dfn, coef[colIndex], order, terms, x[i]) * sx[i]));
793 }
794 }
795 if (generateSigmas) {
796 double sigma;
797 for (i = sigma = 0; i < points; i++) {
798 sigma += sqr(diff[colIndex][i]);
799 }
800 sigma = sqrt(sigma / (points - terms));
801 for (i = 0; i < points; i++) {
802 if (generateSigmas & FLGS_KEEPSMALLEST) {
803 if (sigma < sy[colIndex][i])
804 sy[colIndex][i] = sigma;
805 } else if (generateSigmas & FLGS_KEEPLARGEST) {
806 if (sigma > sy[colIndex][i])
807 sy[colIndex][i] = sigma;
808 } else {
809 sy[colIndex][i] = sigma;
810 }
811 }
812 for (i = 0; i < pointsOrig; i++) {
813 if (generateSigmas & FLGS_KEEPSMALLEST) {
814 if (sigma < sy0[colIndex][i])
815 sy0[colIndex][i] = sigma;
816 } else if (generateSigmas & FLGS_KEEPLARGEST) {
817 if (sigma > sy0[colIndex][i])
818 sy0[colIndex][i] = sigma;
819 } else {
820 sy0[colIndex][i] = sigma;
821 }
822 }
823 }
824 }
825 }
826
827 if (reviseOrders & REVPOW_ACTIVE) {
828 double bestChi;
829 long bestTerms, newBest;
830 int32_t *bestOrder;
831
832 bestTerms = terms;
833 bestOrder = tmalloc(sizeof(*bestOrder) * bestTerms);
834 for (ip = 0; ip < terms; ip++)
835 bestOrder[ip] = order[ip];
836 /* do a fit */
837 for (colIndex = 0; colIndex < numYNames; colIndex++) {
838 isFit[colIndex] = lsfg(x, y[colIndex], sy[colIndex], points, bestTerms, bestOrder, coef[colIndex], coefSigma[colIndex], &bestChi, diff[colIndex], basis_fn);
839 if (!isFit[colIndex]) {
840 fprintf(stderr, "Column %s: ", yNames[colIndex]);
841 SDDS_Bomb("revise-orders fit failed.");
842 if (reviseOrders & REVPOW_VERBOSE) {
843 fprintf(stderr, "Column %s: ", yNames[colIndex]);
844 fputs("fit to revise orders:", stderr);
845 print_coefs(stderr, xOffset, xScaleFactor, chebyshev, coef[colIndex], (ySigmasValid ? coefSigma[colIndex] : NULL), bestOrder, bestTerms, bestChi, normTerm, "");
846 fprintf(stderr, "unweighted rms deviation from fit: %21.15le\n", rms_average(diff[colIndex], points));
847 }
848 }
849
850 do {
851 newBest = 0;
852 terms = bestTerms - 1;
853 for (ip = bestTerms - 1; ip >= 0; ip--) {
854 for (i = j = 0; i < bestTerms; i++)
855 if (i != ip)
856 order[j++] = bestOrder[i];
857 isFit[colIndex] = lsfg(x, y[colIndex], sy[colIndex], points, terms, order, coef[colIndex], coefSigma[colIndex], &chi[colIndex], diff[colIndex], basis_fn);
858 if (!isFit[colIndex]) {
859 fprintf(stderr, "Column %s: ", yNames[colIndex]);
860 SDDS_Bomb("revise-orders fit failed.");
861 }
862 if (reviseOrders & REVPOW_VERBOSE) {
863 fprintf(stderr, "Column %s: ", yNames[colIndex]);
864 fputs("new trial fit:", stderr);
865 print_coefs(stderr, xOffset, xScaleFactor, chebyshev, coef[colIndex], (ySigmasValid ? coefSigma[colIndex] : NULL), order, terms, chi[colIndex], normTerm, "");
866 fprintf(stderr, "unweighted rms deviation from fit: %21.15le\n", rms_average(diff[colIndex], points));
867 }
868 if (chi[colIndex] - bestChi < revpowThreshold) {
869 bestChi = chi[colIndex];
870 bestTerms = terms;
871 newBest = 1;
872 for (i = 0; i < terms; i++)
873 bestOrder[i] = order[i];
874 if (reviseOrders & REVPOW_VERBOSE) {
875 fputs("new best fit:", stderr);
876 print_coefs(stderr, xOffset, xScaleFactor, chebyshev, coef[colIndex], (ySigmasValid ? coefSigma[colIndex] : NULL), bestOrder, bestTerms, bestChi, normTerm, "");
877 fprintf(stderr, "unweighted rms deviation from fit: %21.15le\n", rms_average(diff[colIndex], points));
878 }
879 break;
880 }
881 }
882 if (bestTerms == 1)
883 break;
884 } while (newBest);
885 terms = bestTerms;
886 for (ip = 0; ip < terms; ip++)
887 order[ip] = bestOrder[ip];
888 free(bestOrder);
889 reviseOrders = 0;
890 }
891 }
892
893 if (!outputInitialized) {
894 initializeOutputFile(&SDDSout, &SDDSoutInfo, output, infoFile, &SDDSin, input, xName, yNames, xSigmaName, ySigmaNames, ySigmasValid, order, terms, chebyshev, numYNames, copyParameters);
895 free(output);
896 outputInitialized = 1;
897 }
898 if (evalParameters.file)
899 setupEvaluationFile(&evalParameters, xName, yNames, numYNames, &SDDSin);
900
901 rmsResidual = tmalloc(sizeof(double) * numYNames);
902 for (colIndex = 0; colIndex < numYNames; colIndex++) {
903 isFit[colIndex] = lsfg(x, y[colIndex], sy[colIndex], points, terms, order, coef[colIndex], coefSigma[colIndex], &chi[colIndex], diff[colIndex], basis_fn);
904 if (isFit[colIndex]) {
905 rmsResidual[colIndex] = rms_average(diff[colIndex], points);
906 if (verbose) {
907 fprintf(stderr, "Column: %s\n", yNames[colIndex]);
908 print_coefs(stderr, xOffset, xScaleFactor, chebyshev, coef[colIndex], (ySigmasValid ? coefSigma[colIndex] : NULL), order, terms, chi[colIndex], normTerm, "");
909 fprintf(stderr, "unweighted rms deviation from fit: %21.15le\n", rmsResidual[colIndex]);
910 }
911 } else if (verbose)
912 fprintf(stderr, "fit failed for %s.\n", yNames[colIndex]);
913
914 if (evalParameters.file)
915 makeEvaluationTable(&evalParameters, x, points, coef[colIndex], order, terms, xName, yNames, numYNames, colIndex);
916 }
917
918 if (outputInitialized) {
919 if (!SDDS_StartPage(&SDDSout, rangeFitOnly ? pointsOrig : points) ||
920 (infoFile && !SDDS_StartPage(&SDDSoutInfo, terms)))
921 bomb("A", NULL);
922 if (copyParameters) {
923 if (!SDDS_CopyParameters(&SDDSout, &SDDSin))
924 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
925 if (infoFile && !SDDS_CopyParameters(&SDDSoutInfo, &SDDSin))
926 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
927 }
928 if (!SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_INDEX, rangeFitOnly ? xOrig : x, rangeFitOnly ? pointsOrig : points, ix) ||
929 (infoFile && !SDDS_SetColumnFromLongs(&SDDSoutInfo, SDDS_SET_BY_INDEX, order, terms, iOrder)))
930 bomb("B", NULL);
931 for (colIndex = 0; colIndex < numYNames; colIndex++) {
932 if (rangeFitOnly) {
933 double *residual, rmsResidual0;
934 compareOriginalToFit(xOrig, yOrig[colIndex], &residual, pointsOrig, &rmsResidual0, coef[colIndex], order, terms);
935 if (!SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_INDEX, yOrig[colIndex], pointsOrig, iy[colIndex]) ||
936 !SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_INDEX, residual, pointsOrig, iResidual[colIndex]))
937 bomb("C", NULL);
938 for (i = 0; i < pointsOrig; i++)
939 residual[i] = yOrig[colIndex][i] - residual[i]; /* computes fit values from residual and y */
940 if (!SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_INDEX, residual, pointsOrig, iFit[colIndex]))
941 bomb("D", NULL);
942 free(residual);
943 } else {
944 for (i = 0; i < points; i++)
945 diff[colIndex][i] = -diff[colIndex][i]; /* convert from (Fit-y) to (y-Fit) to get residual */
946 if (!SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_INDEX, y[colIndex], points, iy[colIndex]) ||
947 !SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_INDEX, diff[colIndex], points, iResidual[colIndex]))
948 bomb("C", NULL);
949 for (i = 0; i < points; i++)
950 diff[colIndex][i] = y[colIndex][i] - diff[colIndex][i]; /* computes fit values from residual and y */
951 if (!SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_INDEX, diff[colIndex], points, iFit[colIndex]))
952 bomb("D", NULL);
953 }
954 }
955 if (ixSigma != -1 &&
956 !SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_INDEX, rangeFitOnly ? sxOrig : sx, rangeFitOnly ? pointsOrig : points, ixSigma))
957 bomb("E", NULL);
958 for (colIndex = 0; colIndex < numYNames; colIndex++) {
959 if (ySigmasValid && iySigma[colIndex] != -1 &&
960 !SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_INDEX, rangeFitOnly ? syOrig[colIndex] : sy[colIndex], rangeFitOnly ? pointsOrig : points, iySigma[colIndex]))
961 bomb("F", NULL);
962
963 if (infoFile) {
964 termIndex = coefficient_index(order, terms, 0);
965 if (iIntercept[colIndex] != -1 &&
966 !SDDS_SetParameters(&SDDSoutInfo, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE, iIntercept[colIndex], coef[colIndex][termIndex], -1))
967 bomb("G", NULL);
968 if (iInterceptSigma[colIndex] != -1 &&
969 !SDDS_SetParameters(&SDDSoutInfo, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE, iInterceptSigma[colIndex], coefSigma[colIndex][termIndex], -1))
970 bomb("H", NULL);
971
972 termIndex = coefficient_index(order, terms, 1);
973 if (iSlope[colIndex] != -1 &&
974 !SDDS_SetParameters(&SDDSoutInfo, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE, iSlope[colIndex], coef[colIndex][termIndex], -1))
975 bomb("I", NULL);
976 if (iSlopeSigma[colIndex] != -1 &&
977 !SDDS_SetParameters(&SDDSoutInfo, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE, iSlopeSigma[colIndex], coefSigma[colIndex][termIndex], -1))
978 bomb("J", NULL);
979
980 termIndex = coefficient_index(order, terms, 2);
981 if (iCurvature[colIndex] != -1 &&
982 !SDDS_SetParameters(&SDDSoutInfo, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE, iCurvature[colIndex], coef[colIndex][termIndex], -1))
983 bomb("K", NULL);
984 if (iCurvatureSigma[colIndex] != -1 &&
985 !SDDS_SetParameters(&SDDSoutInfo, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE, iCurvatureSigma[colIndex], coefSigma[colIndex][termIndex], -1))
986 bomb("L", NULL);
987 if (iFitLabel[colIndex] != -1) {
988 makeFitLabel(fitLabelBuffer, SDDS_MAXLINE, fitLabelFormat, coef[colIndex], order, terms, colIndex);
989 if (!SDDS_SetParameters(&SDDSoutInfo, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE, iFitLabel[colIndex], fitLabelBuffer, -1))
990 bomb("M", NULL);
991 }
992 if (!SDDS_SetColumnFromDoubles(&SDDSoutInfo, SDDS_SET_BY_INDEX, coef[colIndex], terms, iCoefficient[colIndex]) ||
993 (ySigmasValid &&
994 !SDDS_SetColumnFromDoubles(&SDDSoutInfo, SDDS_SET_BY_INDEX, coefSigma[colIndex], terms, iCoefficientSigma[colIndex])))
995 bomb("N", NULL);
996 if (!SDDS_SetParameters(&SDDSoutInfo, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE,
997 iRmsResidual[colIndex], rmsResidual[colIndex], iChiSq[colIndex], chi[colIndex],
998 iTerms, terms, iSigLevel[colIndex], ChiSqrSigLevel(chi[colIndex], points - terms),
999 iOffset, xOffset, iFactor, xScaleFactor, iFitIsValid[colIndex], isFit[colIndex] ? 'y' : 'n', -1))
1000 bomb("O", NULL);
1001 }
1002
1003 termIndex = coefficient_index(order, terms, 0);
1004 if (iInterceptO[colIndex] != -1 &&
1005 !SDDS_SetParameters(&SDDSout, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE, iInterceptO[colIndex], coef[colIndex][termIndex], -1))
1006 bomb("G", NULL);
1007 if (iInterceptSigmaO[colIndex] != -1 &&
1008 !SDDS_SetParameters(&SDDSout, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE, iInterceptSigmaO[colIndex], coefSigma[colIndex][termIndex], -1))
1009 bomb("H", NULL);
1010
1011 termIndex = coefficient_index(order, terms, 1);
1012 if (iSlopeO[colIndex] != -1 &&
1013 !SDDS_SetParameters(&SDDSout, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE, iSlopeO[colIndex], coef[colIndex][termIndex], -1))
1014 bomb("I", NULL);
1015 if (iSlopeSigmaO[colIndex] != -1 &&
1016 !SDDS_SetParameters(&SDDSout, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE, iSlopeSigmaO[colIndex], coefSigma[colIndex][termIndex], -1))
1017 bomb("J", NULL);
1018
1019 termIndex = coefficient_index(order, terms, 2);
1020 if (iCurvatureO[colIndex] != -1 &&
1021 !SDDS_SetParameters(&SDDSout, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE, iCurvatureO[colIndex], coef[colIndex][termIndex], -1))
1022 bomb("K", NULL);
1023 if (iCurvatureSigmaO[colIndex] != -1 &&
1024 !SDDS_SetParameters(&SDDSout, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE, iCurvatureSigmaO[colIndex], coefSigma[colIndex][termIndex], -1))
1025 bomb("L", NULL);
1026 if (iFitLabelO[colIndex] != -1) {
1027 makeFitLabel(fitLabelBuffer, SDDS_MAXLINE, fitLabelFormat, coef[colIndex], order, terms, colIndex);
1028 if (!SDDS_SetParameters(&SDDSout, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE, iFitLabelO[colIndex], fitLabelBuffer, -1))
1029 bomb("M", NULL);
1030 }
1031 if (!SDDS_SetParameters(&SDDSout, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE,
1032 iRmsResidualO[colIndex], rmsResidual[colIndex], iChiSqO[colIndex], chi[colIndex],
1033 iTermsO, terms, iSigLevelO[colIndex], ChiSqrSigLevel(chi[colIndex], points - terms),
1034 iOffsetO, xOffset, iFactorO, xScaleFactor, iFitIsValidO[colIndex], isFit[colIndex] ? 'y' : 'n', -1))
1035 bomb("O", NULL);
1036 }
1037 if (!SDDS_WritePage(&SDDSout) || (infoFile && !SDDS_WritePage(&SDDSoutInfo)))
1038 bomb("O", NULL);
1039 }
1040 if (xOrig != x)
1041 free(xOrig);
1042 if (sxOrig != sx)
1043 free(sxOrig);
1044 free(x);
1045 free(sx);
1046 for (colIndex = 0; colIndex < numYNames; colIndex++) {
1047 free(diff[colIndex]);
1048 free(sy[colIndex]);
1049 if (yOrig[colIndex] != y[colIndex])
1050 free(yOrig[colIndex]);
1051 if (syOrig && sy0 && syOrig[colIndex] != sy0[colIndex])
1052 free(syOrig[colIndex]);
1053 free(y[colIndex]);
1054 if (sy0 && sy0[colIndex])
1055 free(sy0[colIndex]);
1056 }
1057 }
1058 return (EXIT_SUCCESS);
1059}
1060
1061void print_coefs(FILE * fpo, double xOffset, double xScaleFactor, long chebyshev, double *coef, double *coefSigma,
1062 int32_t *order, long terms, double chi, long normTerm, char *prepend) {
1063 long i;
1064
1065 if (chebyshev)
1066 fprintf(fpo, "%s%ld-term Chebyshev T polynomial least-squares fit about x=%21.15le, scaled by %21.15le:\n", prepend, terms, xOffset, xScaleFactor);
1067 else
1068 fprintf(fpo, "%s%ld-term polynomial least-squares fit about x=%21.15le:\n", prepend, terms, xOffset);
1069 if (normTerm >= 0 && terms > normTerm) {
1070 if (coef[normTerm] != 0)
1071 fprintf(fpo, "%s coefficients are normalized with factor %21.15le to make a[%ld]==1\n", prepend, coef[normTerm], (order ? order[normTerm] : normTerm));
1072 else {
1073 fprintf(fpo, "%s can't normalize coefficients as requested: a[%ld]==0\n", prepend, (order ? order[normTerm] : normTerm));
1074 normTerm = -1;
1075 }
1076 } else
1077 normTerm = -1;
1078
1079 for (i = 0; i < terms; i++) {
1080 fprintf(fpo, "%sa[%ld] = %21.15le ", prepend, (order ? order[i] : i), (normTerm < 0 ? coef[i] : coef[i] / coef[normTerm]));
1081 if (coefSigma)
1082 fprintf(fpo, "+/- %21.15le\n", (normTerm < 0 ? coefSigma[i] : coefSigma[i] / fabs(coef[normTerm])));
1083 else
1084 fputc('\n', fpo);
1085 }
1086 if (coefSigma)
1087 fprintf(fpo, "%sreduced chi-squared = %21.15le\n", prepend, chi);
1088}
1089
1090void RemoveElementFromStringArray(char **array, long index, long length) {
1091 long lh;
1092
1093 for (lh = index; lh < length - 1; lh++)
1094 array[lh] = array[lh + 1];
1095}
1096
1097char **RemoveNonNumericColumnsFromNameArray(SDDS_DATASET * SDDSin, char **columns, int32_t *numColumns) {
1098 long i, numNumericColumns = *numColumns;
1099
1100 for (i = 0; i < *numColumns; i++) {
1101 if (SDDS_CheckColumn(SDDSin, columns[i], NULL, SDDS_ANY_NUMERIC_TYPE, NULL)) {
1102 printf("Removing %s because not a numeric type.\n", columns[i]);
1103 RemoveElementFromStringArray(columns, i, *numColumns);
1104 numNumericColumns--;
1105 }
1106 }
1107
1108 *numColumns = numNumericColumns;
1109 return (columns);
1110}
1111
1112char **ResolveColumnNames(SDDS_DATASET * SDDSin, char **wildcardList, long length, int32_t *numYNames) {
1113 char **result;
1114 long i;
1115
1116 /* initially set the columns of interest to none, to make SDDS_OR work below */
1117 SDDS_SetColumnsOfInterest(SDDSin, SDDS_MATCH_STRING, "", SDDS_AND);
1118 for (i = 0; i < length; i++) {
1119 SDDS_SetColumnsOfInterest(SDDSin, SDDS_MATCH_STRING, wildcardList[i], SDDS_OR);
1120 }
1121
1122 if (!(result = SDDS_GetColumnNames(SDDSin, numYNames)) || *numYNames == 0)
1123 bomb("Error matching columns in ResolveColumnNames: No matches.", NULL);
1124
1125 result = RemoveNonNumericColumnsFromNameArray(SDDSin, result, numYNames);
1126 return (result);
1127}
1128
1129char **GenerateYSigmaNames(char *controlString, char **yNames, long numYNames) {
1130 long i, nameLength;
1131 char **result, sigmaName[MAX_Y_SIGMA_NAME_SIZE];
1132
1133 result = tmalloc(sizeof(char *) * numYNames);
1134 for (i = 0; i < numYNames; i++) {
1135 sprintf(sigmaName, controlString, yNames[i]);
1136 nameLength = strlen(sigmaName);
1137 result[i] = tmalloc(sizeof(char) * (nameLength + 1));
1138 strcpy(result[i], sigmaName);
1139 }
1140 return (result);
1141}
1142
1143void makeFitLabel(char *buffer, long bufsize, char *fitLabelFormat, double *coef, int32_t *order, long terms, long colIndex) {
1144 long i;
1145 static char buffer1[SDDS_MAXLINE], buffer2[SDDS_MAXLINE];
1146
1147 sprintf(buffer, "%s = ", ySymbols[colIndex]);
1148 for (i = 0; i < terms; i++) {
1149 if (order[i] == 0)
1150 sprintf(buffer1, fitLabelFormat, coef[i]);
1151 else {
1152 if (coef[i] >= 0) {
1153 strcpy(buffer1, " +");
1154 sprintf(buffer1 + 2, fitLabelFormat, coef[i]);
1155 } else
1156 sprintf(buffer1, fitLabelFormat, coef[i]);
1157 strcat(buffer1, "*");
1158 strcat(buffer1, xSymbol);
1159 if (order[i] > 1) {
1160 sprintf(buffer2, "$a%" PRId32 "$n", order[i]);
1161 strcat(buffer1, buffer2);
1162 }
1163 }
1164 if ((long)(strlen(buffer) + strlen(buffer1)) > (long)(0.95 * bufsize)) {
1165 fprintf(stderr, "buffer overflow making fit label!\n");
1166 return;
1167 }
1168 strcat(buffer, buffer1);
1169 }
1170}
1171
1172double rms_average(double *x, int64_t n) {
1173 double sum2;
1174 int64_t i;
1175
1176 for (i = sum2 = 0; i < n; i++)
1177 sum2 += sqr(x[i]);
1178
1179 return (sqrt(sum2 / n));
1180}
1181
1182long coefficient_index(int32_t * order, long terms, long order_of_interest) {
1183 long i;
1184 for (i = 0; i < terms; i++)
1185 if (order[i] == order_of_interest)
1186 return (i);
1187 return (-1);
1188}
1189
1190void checkInputFile(SDDS_DATASET * SDDSin, char *xName, char **yNames, char *xSigmaName, char **ySigmaNames, long numYNames) {
1191 char *ptr = NULL;
1192 long i;
1193
1194 if (!(ptr = SDDS_FindColumn(SDDSin, FIND_NUMERIC_TYPE, xName, NULL)))
1195 SDDS_Bomb("x column doesn't exist or is nonnumeric");
1196 free(ptr);
1197
1198 /* y columns don't need to be checked because located using SDDS_SetColumnsOfInterest */
1199
1200 ptr = NULL;
1201 if (xSigmaName && !(ptr = SDDS_FindColumn(SDDSin, FIND_NUMERIC_TYPE, xSigmaName, NULL)))
1202 SDDS_Bomb("x sigma column doesn't exist or is nonnumeric");
1203 if (ptr)
1204 free(ptr);
1205
1206 if (ySigmaNames) {
1207 for (i = 0; i < numYNames; i++) {
1208 ptr = NULL;
1209 if (!(ptr = SDDS_FindColumn(SDDSin, FIND_NUMERIC_TYPE, ySigmaNames[i], NULL)))
1210 SDDS_Bomb("y sigma column doesn't exist or is nonnumeric");
1211 if (ptr)
1212 free(ptr);
1213 }
1214 }
1215}
1216
1217char ***initializeOutputFile(SDDS_DATASET * SDDSout, SDDS_DATASET * SDDSoutInfo, char *output, char *outputInfo,
1218 SDDS_DATASET *SDDSin, char *input, char *xName, char **yNames, char *xSigmaName,
1219 char **ySigmaNames, long sigmasValid, int32_t *order, long terms, long isChebyshev,
1220 long numCols, long copyParameters) {
1221 char buffer[SDDS_MAXLINE], buffer1[SDDS_MAXLINE], buffer2[SDDS_MAXLINE], buffer3[SDDS_MAXLINE];
1222 char *xUnits, *yUnits, ***coefUnits;
1223 long i, colIndex;
1224
1225 /* all array names followed by an 'O' contain the index of the parameter in the main output file; others refer to
1226 parameters in the infoFile */
1227 coefUnits = tmalloc(sizeof(char **) * numCols);
1228 ySymbols = tmalloc(sizeof(char *) * numCols);
1229 iIntercept = tmalloc(sizeof(long) * numCols);
1230 iInterceptO = tmalloc(sizeof(long) * numCols);
1231 iInterceptSigma = tmalloc(sizeof(long) * numCols);
1232 iInterceptSigmaO = tmalloc(sizeof(long) * numCols);
1233 iSlope = tmalloc(sizeof(long) * numCols);
1234 iSlopeO = tmalloc(sizeof(long) * numCols);
1235 iSlopeSigma = tmalloc(sizeof(long) * numCols);
1236 iSlopeSigmaO = tmalloc(sizeof(long) * numCols);
1237 iCurvature = tmalloc(sizeof(long) * numCols);
1238 iCurvatureO = tmalloc(sizeof(long) * numCols);
1239 iCurvatureSigma = tmalloc(sizeof(long) * numCols);
1240 iCurvatureSigmaO = tmalloc(sizeof(long) * numCols);
1241 iChiSq = tmalloc(sizeof(long) * numCols);
1242 iChiSqO = tmalloc(sizeof(long) * numCols);
1243 iRmsResidual = tmalloc(sizeof(long) * numCols);
1244 iRmsResidualO = tmalloc(sizeof(long) * numCols);
1245 iSigLevel = tmalloc(sizeof(long) * numCols);
1246 iSigLevelO = tmalloc(sizeof(long) * numCols);
1247 iFitIsValid = tmalloc(sizeof(long) * numCols);
1248 iFitIsValidO = tmalloc(sizeof(long) * numCols);
1249 iFitLabel = tmalloc(sizeof(long) * numCols);
1250 iFitLabelO = tmalloc(sizeof(long) * numCols);
1251 iy = tmalloc(sizeof(long) * numCols);
1252 iySigma = tmalloc(sizeof(long) * numCols);
1253 iFit = tmalloc(sizeof(long) * numCols);
1254 iResidual = tmalloc(sizeof(long) * numCols);
1255
1256 for (colIndex = 0; colIndex < numCols; colIndex++) {
1257 ySymbols[colIndex] = NULL;
1258 coefUnits[colIndex] = tmalloc(sizeof(char *) * terms);
1259 iInterceptSigma[colIndex] = -1;
1260 iInterceptSigmaO[colIndex] = -1;
1261 iIntercept[colIndex] = -1;
1262 iInterceptO[colIndex] = -1;
1263 iInterceptSigma[colIndex] = -1;
1264 iInterceptSigmaO[colIndex] = -1;
1265 iSlope[colIndex] = -1;
1266 iSlopeO[colIndex] = -1;
1267 iSlopeSigma[colIndex] = -1;
1268 iSlopeSigmaO[colIndex] = -1;
1269 iCurvature[colIndex] = -1;
1270 iCurvatureO[colIndex] = -1;
1271 iCurvatureSigma[colIndex] = -1;
1272 iCurvatureSigmaO[colIndex] = -1;
1273 iChiSq[colIndex] = -1;
1274 iChiSqO[colIndex] = -1;
1275 iRmsResidual[colIndex] = -1;
1276 iRmsResidualO[colIndex] = -1;
1277 iSigLevel[colIndex] = -1;
1278 iSigLevelO[colIndex] = -1;
1279 iFitIsValid[colIndex] = -1;
1280 iFitIsValidO[colIndex] = -1;
1281 iFitLabel[colIndex] = -1;
1282 iFitLabelO[colIndex] = -1;
1283 iy[colIndex] = -1;
1284 iySigma[colIndex] = -1;
1285 iFit[colIndex] = -1;
1286 iResidual[colIndex] = -1;
1287 }
1288
1289 if (!SDDS_InitializeOutput(SDDSout, SDDS_BINARY, 0, NULL, "sddsmpfit output: fitted data", output) ||
1290 !SDDS_TransferColumnDefinition(SDDSout, SDDSin, xName, NULL) ||
1291 SDDS_GetColumnInformation(SDDSout, "symbol", &xSymbol, SDDS_GET_BY_NAME, xName) != SDDS_STRING ||
1292 (xSigmaName && !SDDS_TransferColumnDefinition(SDDSout, SDDSin, xSigmaName, NULL)))
1293 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1294
1295 for (colIndex = 0; colIndex < numCols; colIndex++) {
1296 if (!SDDS_TransferColumnDefinition(SDDSout, SDDSin, yNames[colIndex], NULL) ||
1297 SDDS_GetColumnInformation(SDDSout, "symbol", &ySymbols[colIndex], SDDS_GET_BY_NAME, yNames[colIndex]) != SDDS_STRING ||
1298 (ySigmaNames && !SDDS_TransferColumnDefinition(SDDSout, SDDSin, ySigmaNames[colIndex], NULL)))
1299 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1300 }
1301 if (!xSymbol || SDDS_StringIsBlank(xSymbol))
1302 xSymbol = xName;
1303 for (colIndex = 0; colIndex < numCols; colIndex++)
1304 if (!ySymbols[colIndex] || SDDS_StringIsBlank(ySymbols[colIndex]))
1305 ySymbols[colIndex] = yNames[colIndex];
1306 ix = SDDS_GetColumnIndex(SDDSout, xName);
1307 for (colIndex = 0; colIndex < numCols; colIndex++) {
1308 iy[colIndex] = SDDS_GetColumnIndex(SDDSout, yNames[colIndex]);
1309 if (ySigmaNames)
1310 iySigma[colIndex] = SDDS_GetColumnIndex(SDDSout, ySigmaNames[colIndex]);
1311 }
1312 if (xSigmaName)
1313 ixSigma = SDDS_GetColumnIndex(SDDSout, xSigmaName);
1314 if (SDDS_NumberOfErrors())
1315 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1316
1317 for (colIndex = 0; colIndex < numCols; colIndex++) {
1318 sprintf(buffer, "%sFit", yNames[colIndex]);
1319 sprintf(buffer1, "Fit[%s]", ySymbols[colIndex]);
1320 if (!SDDS_TransferColumnDefinition(SDDSout, SDDSin, yNames[colIndex], buffer) ||
1321 !SDDS_ChangeColumnInformation(SDDSout, "symbol", buffer1, SDDS_SET_BY_NAME, buffer))
1322 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1323 if ((iFit[colIndex] = SDDS_GetColumnIndex(SDDSout, buffer)) < 0)
1324 SDDS_Bomb("unable to get index of just-defined fit output column");
1325
1326 sprintf(buffer, "%sResidual", yNames[colIndex]);
1327 sprintf(buffer1, "Residual[%s]", ySymbols[colIndex]);
1328 if (!SDDS_TransferColumnDefinition(SDDSout, SDDSin, yNames[colIndex], buffer) ||
1329 !SDDS_ChangeColumnInformation(SDDSout, "symbol", buffer1, SDDS_SET_BY_NAME, buffer))
1330 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1331 if (!(iResidual[colIndex] = SDDS_GetColumnIndex(SDDSout, buffer)))
1332 SDDS_Bomb("unable to get index of just-defined residual output column");
1333
1334 if (sigmasValid && !ySigmaNames) {
1335 sprintf(buffer, "%sSigma", yNames[colIndex]);
1336 if (!SDDS_TransferColumnDefinition(SDDSout, SDDSin, yNames[colIndex], buffer))
1337 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1338 iySigma[colIndex] = SDDS_GetColumnIndex(SDDSout, buffer);
1339 if (ySymbols[colIndex] && !SDDS_StringIsBlank(ySymbols[colIndex])) {
1340 sprintf(buffer1, "Sigma[%s]", ySymbols[colIndex]);
1341 if (!SDDS_ChangeColumnInformation(SDDSout, "symbol", buffer1, SDDS_SET_BY_NAME, buffer))
1342 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1343 }
1344 }
1345
1346 if (!(coefUnits[colIndex] = makeCoefficientUnits(SDDSout, xName, yNames[colIndex], order, terms)))
1347 SDDS_Bomb("unable to make coefficient units");
1348 }
1349
1350 if (outputInfo && !SDDS_InitializeOutput(SDDSoutInfo, SDDS_BINARY, 0, NULL, "sddsmpfit output: fit information", outputInfo))
1351 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1352
1353 if (outputInfo) {
1354 if ((SDDS_DefineColumn(SDDSoutInfo, "Order", NULL, NULL, "Order of term in fit", NULL, SDDS_LONG, 0) < 0) ||
1355 SDDS_DefineParameter(SDDSoutInfo, "Basis", NULL, NULL, "Function basis for fit", NULL, SDDS_STRING, isChebyshev ? "Chebyshev T polynomials" : "ordinary polynomials") < 0)
1356 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1357
1358 if ((iTerms = SDDS_DefineParameter(SDDSoutInfo, "Terms", NULL, NULL, "Number of terms in fit", NULL, SDDS_LONG, NULL)) < 0)
1359 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1360
1361 if (SDDS_GetColumnInformation(SDDSout, "units", &xUnits, SDDS_GET_BY_NAME, xName) != SDDS_STRING)
1362 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1363 sprintf(buffer, "%sOffset", xName);
1364 sprintf(buffer1, "Offset of %s for fit", xName);
1365 if ((iOffset = SDDS_DefineParameter(SDDSoutInfo, buffer, NULL, xUnits, buffer1, NULL, SDDS_DOUBLE, NULL)) < 0)
1366 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1367 sprintf(buffer, "%sScale", xName);
1368 sprintf(buffer1, "Scale factor of %s for fit", xName);
1369 if ((iFactor = SDDS_DefineParameter(SDDSoutInfo, buffer, NULL, xUnits, buffer1, NULL, SDDS_DOUBLE, NULL)) < 0)
1370 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1371
1372 for (colIndex = 0; colIndex < numCols; colIndex++) {
1373
1374 sprintf(buffer1, "%sCoefficient", yNames[colIndex]);
1375 sprintf(buffer2, "%sCoefficientSigma", yNames[colIndex]);
1376 sprintf(buffer3, "%sCoefficientUnits", yNames[colIndex]);
1377
1378 if (SDDS_DefineColumn(SDDSoutInfo, buffer1, NULL, "[CoefficientUnits]", "Coefficient of term in fit", NULL, SDDS_DOUBLE, 0) < 0 ||
1379 (sigmasValid && SDDS_DefineColumn(SDDSoutInfo, buffer2, "$gs$r$ba$n", "[CoefficientUnits]", "sigma of coefficient of term in fit", NULL, SDDS_DOUBLE, 0) < 0))
1380 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1381
1382 iOrder = SDDS_GetColumnIndex(SDDSoutInfo, "Order");
1383 iCoefficient[colIndex] = SDDS_GetColumnIndex(SDDSoutInfo, buffer1);
1384 iCoefficientSigma[colIndex] = SDDS_GetColumnIndex(SDDSoutInfo, buffer2);
1385 iCoefficientUnits[colIndex] = SDDS_GetColumnIndex(SDDSoutInfo, buffer3);
1386
1387 sprintf(buffer1, "%sReducedChiSquared", yNames[colIndex]);
1388 sprintf(buffer2, "%sRmsResidual", yNames[colIndex]);
1389 sprintf(buffer3, "%sSignificanceLevel", yNames[colIndex]);
1390
1391 if ((iChiSq[colIndex] = SDDS_DefineParameter(SDDSoutInfo, buffer1, "$gh$r$a2$n/(N-M)", NULL,
1392 "Reduced chi-squared of fit",
1393 NULL, SDDS_DOUBLE, NULL)) < 0 ||
1394 SDDS_GetColumnInformation(SDDSout, "units", &yUnits, SDDS_GET_BY_NAME, yNames[colIndex]) != SDDS_STRING ||
1395 (iRmsResidual[colIndex] =
1396 SDDS_DefineParameter(SDDSoutInfo, buffer2, "$gs$r$bres$n", yUnits, "RMS residual of fit", NULL, SDDS_DOUBLE, NULL)) < 0 ||
1397 (iSigLevel[colIndex] = SDDS_DefineParameter(SDDSoutInfo, buffer3, NULL, NULL, "Probability that data is from fit function", NULL, SDDS_DOUBLE, NULL)) < 0)
1398 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1399 if (yUnits)
1400 free(yUnits);
1401
1402 sprintf(buffer, "%sFitIsValid", yNames[colIndex]);
1403 if ((iFitIsValid[colIndex] = SDDS_DefineParameter(SDDSoutInfo, buffer, NULL, NULL, NULL, NULL, SDDS_CHARACTER, NULL)) < 0)
1404 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1405
1406 if (!isChebyshev) {
1407
1408 sprintf(buffer, "%sSddsmpfitlabel", yNames[colIndex]);
1409 iFitLabel[colIndex] = SDDS_DefineParameter(SDDSoutInfo, buffer, NULL, NULL, NULL, NULL, SDDS_STRING, NULL);
1410 if ((i = coefficient_index(order, terms, 0)) >= 0) {
1411 sprintf(buffer, "%sIntercept", yNames[colIndex]);
1412 iIntercept[colIndex] = SDDS_DefineParameter(SDDSoutInfo, buffer, buffer, coefUnits[colIndex][i], "Intercept of fit", NULL, SDDS_DOUBLE, NULL);
1413 sprintf(buffer, "%sInterceptSigma", yNames[colIndex]);
1414 if (sigmasValid)
1415 iInterceptSigma[colIndex] = SDDS_DefineParameter(SDDSoutInfo, buffer, buffer, coefUnits[colIndex][i], "Sigma of intercept of fit", NULL, SDDS_DOUBLE, NULL);
1416 }
1417 sprintf(buffer, "%sSlope", yNames[colIndex]);
1418 if ((i = coefficient_index(order, terms, 1)) >= 0) {
1419 iSlope[colIndex] = SDDS_DefineParameter(SDDSoutInfo, buffer, buffer, coefUnits[colIndex][i], "Slope of fit", NULL, SDDS_DOUBLE, NULL);
1420 if (sigmasValid) {
1421 sprintf(buffer, "%sSlopeSigma", yNames[colIndex]);
1422 iSlopeSigma[colIndex] = SDDS_DefineParameter(SDDSoutInfo, buffer, buffer, coefUnits[colIndex][i], "Sigma of slope of fit", NULL, SDDS_DOUBLE, NULL);
1423 }
1424 }
1425 if ((i = coefficient_index(order, terms, 2)) >= 0) {
1426 sprintf(buffer, "%sCurvature", yNames[colIndex]);
1427 iCurvature[colIndex] = SDDS_DefineParameter(SDDSoutInfo, buffer, buffer, coefUnits[colIndex][i], "Curvature of fit", NULL, SDDS_DOUBLE, NULL);
1428 if (sigmasValid) {
1429 sprintf(buffer, "%sCurvatureSigma", yNames[colIndex]);
1430 iCurvatureSigma[colIndex] = SDDS_DefineParameter(SDDSoutInfo, buffer, buffer, coefUnits[colIndex][i], "Sigma of curvature of fit", NULL, SDDS_DOUBLE, NULL);
1431 }
1432 }
1433 if (SDDS_NumberOfErrors())
1434 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1435 }
1436 }
1437 }
1438
1439 if (SDDS_DefineParameter(SDDSout, "Basis", NULL, NULL, "Function basis for fit", NULL, SDDS_STRING, isChebyshev ? "Chebyshev T polynomials" : "ordinary polynomials") < 0)
1440 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1441
1442 if ((iTermsO = SDDS_DefineParameter(SDDSout, "Terms", NULL, NULL, "Number of terms in fit", NULL, SDDS_LONG, NULL)) < 0)
1443 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1444
1445 if (SDDS_GetColumnInformation(SDDSout, "units", &xUnits, SDDS_GET_BY_NAME, xName) != SDDS_STRING)
1446 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1447 sprintf(buffer, "%sOffset", xName);
1448 sprintf(buffer1, "Offset of %s for fit", xName);
1449 if ((iOffsetO = SDDS_DefineParameter(SDDSout, buffer, NULL, xUnits, buffer1, NULL, SDDS_DOUBLE, NULL)) < 0)
1450 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1451 sprintf(buffer, "%sScale", xName);
1452 sprintf(buffer1, "Scale factor of %s for fit", xName);
1453 if ((iFactorO = SDDS_DefineParameter(SDDSout, buffer, NULL, xUnits, buffer1, NULL, SDDS_DOUBLE, NULL)) < 0)
1454 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1455
1456 for (colIndex = 0; colIndex < numCols; colIndex++) {
1457
1458 sprintf(buffer1, "%sReducedChiSquared", yNames[colIndex]);
1459 sprintf(buffer2, "%sRmsResidual", yNames[colIndex]);
1460 sprintf(buffer3, "%sSignificanceLevel", yNames[colIndex]);
1461
1462 if ((iChiSqO[colIndex] = SDDS_DefineParameter(SDDSout, buffer1, "$gh$r$a2$n/(N-M)", NULL,
1463 "Reduced chi-squared of fit",
1464 NULL, SDDS_DOUBLE, NULL)) < 0 ||
1465 SDDS_GetColumnInformation(SDDSout, "units", &yUnits, SDDS_GET_BY_NAME, yNames[colIndex]) != SDDS_STRING ||
1466 (iRmsResidualO[colIndex] =
1467 SDDS_DefineParameter(SDDSout, buffer2, "$gs$r$bres$n", yUnits, "RMS residual of fit", NULL, SDDS_DOUBLE, NULL)) < 0 ||
1468 (iSigLevelO[colIndex] = SDDS_DefineParameter(SDDSout, buffer3, NULL, NULL, "Probability that data is from fit function", NULL, SDDS_DOUBLE, NULL)) < 0)
1469 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1470 if (yUnits)
1471 free(yUnits);
1472
1473 sprintf(buffer, "%sFitIsValid", yNames[colIndex]);
1474 if ((iFitIsValidO[colIndex] = SDDS_DefineParameter(SDDSout, buffer, NULL, NULL, NULL, NULL, SDDS_CHARACTER, NULL)) < 0)
1475 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1476
1477 if (!isChebyshev) {
1478 sprintf(buffer, "%sSddsmpfitlabel", yNames[colIndex]);
1479 iFitLabelO[colIndex] = SDDS_DefineParameter(SDDSout, buffer, NULL, NULL, NULL, NULL, SDDS_STRING, NULL);
1480 if ((i = coefficient_index(order, terms, 0)) >= 0) {
1481 sprintf(buffer, "%sIntercept", yNames[colIndex]);
1482 iInterceptO[colIndex] = SDDS_DefineParameter(SDDSout, buffer, buffer, coefUnits[colIndex][i], "Intercept of fit", NULL, SDDS_DOUBLE, NULL);
1483 sprintf(buffer, "%sInterceptSigma", yNames[colIndex]);
1484 if (sigmasValid)
1485 iInterceptSigmaO[colIndex] = SDDS_DefineParameter(SDDSout, buffer, buffer, coefUnits[colIndex][i], "Sigma of intercept of fit", NULL, SDDS_DOUBLE, NULL);
1486 }
1487 sprintf(buffer, "%sSlope", yNames[colIndex]);
1488 if ((i = coefficient_index(order, terms, 1)) >= 0) {
1489 iSlopeO[colIndex] = SDDS_DefineParameter(SDDSout, buffer, buffer, coefUnits[colIndex][i], "Slope of fit", NULL, SDDS_DOUBLE, NULL);
1490 if (sigmasValid) {
1491 sprintf(buffer, "%sSlopeSigma", yNames[colIndex]);
1492 iSlopeSigmaO[colIndex] = SDDS_DefineParameter(SDDSout, buffer, buffer, coefUnits[colIndex][i], "Sigma of slope of fit", NULL, SDDS_DOUBLE, NULL);
1493 }
1494 }
1495 if ((i = coefficient_index(order, terms, 2)) >= 0) {
1496 sprintf(buffer, "%sCurvature", yNames[colIndex]);
1497 iCurvatureO[colIndex] = SDDS_DefineParameter(SDDSout, buffer, buffer, coefUnits[colIndex][i], "Curvature of fit", NULL, SDDS_DOUBLE, NULL);
1498 if (sigmasValid) {
1499 sprintf(buffer, "%sCurvatureSigma", yNames[colIndex]);
1500 iCurvatureSigmaO[colIndex] = SDDS_DefineParameter(SDDSout, buffer, buffer, coefUnits[colIndex][i], "Sigma of curvature of fit", NULL, SDDS_DOUBLE, NULL);
1501 }
1502 }
1503 if (SDDS_NumberOfErrors())
1504 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1505 }
1506 }
1507
1508 if (copyParameters) {
1509 if (!SDDS_TransferAllParameterDefinitions(SDDSout, SDDSin, 0))
1510 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1511 if (outputInfo && !SDDS_TransferAllParameterDefinitions(SDDSoutInfo, SDDSin, 0))
1512 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1513 }
1514
1515 if ((outputInfo && !SDDS_WriteLayout(SDDSoutInfo)) || !SDDS_WriteLayout(SDDSout))
1516 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1517
1518 return coefUnits;
1519}
1520
1521char **makeCoefficientUnits(SDDS_DATASET * SDDSout, char *xName, char *yName, int32_t *order, long terms) {
1522 char *xUnits, *yUnits, buffer[SDDS_MAXLINE];
1523 char **coefUnits = NULL;
1524 long i;
1525
1526 if (!SDDS_GetColumnInformation(SDDSout, "units", &xUnits, SDDS_GET_BY_NAME, xName) ||
1527 !SDDS_GetColumnInformation(SDDSout, "units", &yUnits, SDDS_GET_BY_NAME, yName))
1528 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1529
1530 coefUnits = tmalloc(sizeof(*coefUnits) * terms);
1531 if (!xUnits || SDDS_StringIsBlank(xUnits)) {
1532 if (!yUnits || SDDS_StringIsBlank(yUnits))
1533 SDDS_CopyString(&yUnits, "");
1534 for (i = 0; i < terms; i++)
1535 coefUnits[i] = yUnits;
1536 } else {
1537 if (!yUnits || SDDS_StringIsBlank(yUnits))
1538 SDDS_CopyString(&yUnits, "1");
1539 for (i = 0; i < terms; i++) {
1540 if (order[i] == 0) {
1541 if (strcmp(yUnits, "1") != 0)
1542 SDDS_CopyString(coefUnits + i, yUnits);
1543 else
1544 SDDS_CopyString(coefUnits + i, "");
1545 } else if (strcmp(xUnits, yUnits) == 0) {
1546 if (order[i] > 1)
1547 sprintf(buffer, "1/%s$a%" PRId32 "$n", xUnits, order[i] - 1);
1548 else
1549 strcpy(buffer, "");
1550 SDDS_CopyString(coefUnits + i, buffer);
1551 } else {
1552 if (order[i] > 1)
1553 sprintf(buffer, "%s/%s$a%" PRId32 "$n", yUnits, xUnits, order[i]);
1554 else
1555 sprintf(buffer, "%s/%s", yUnits, xUnits);
1556 SDDS_CopyString(coefUnits + i, buffer);
1557 }
1558 }
1559 }
1560 return coefUnits;
1561}
1562
1563void setupEvaluationFile(EVAL_PARAMETERS * evalParameters, char *xName, char **yName, long yNames, SDDS_DATASET *SDDSin) {
1564 long i;
1565 SDDS_DATASET *SDDSout;
1566 SDDSout = &evalParameters->dataset;
1567 if (!SDDS_InitializeOutput(SDDSout, SDDS_BINARY, 0, NULL, "sddsmpfit output: evaluation of fits", evalParameters->file) ||
1568 !SDDS_TransferColumnDefinition(SDDSout, SDDSin, xName, NULL))
1569 SDDS_Bomb("Problem setting up evaluation file");
1570 for (i = 0; i < yNames; i++)
1571 if (!SDDS_TransferColumnDefinition(SDDSout, SDDSin, yName[i], NULL))
1572 SDDS_Bomb("Problem setting up evaluation file");
1573 if (!SDDS_WriteLayout(SDDSout))
1574 SDDS_Bomb("Problem setting up evaluation file");
1575}
1576
1577void makeEvaluationTable(EVAL_PARAMETERS * evalParameters, double *x, int64_t points, double *coef, int32_t *order, long terms,
1578 char *xName, char **yName, long yNames, long iYName) {
1579 static double *xEval = NULL, *yEval = NULL;
1580 static int64_t maxEvals = 0;
1581 double delta;
1582 int64_t i;
1583
1584 if (!(evalParameters->flags & EVAL_BEGIN_GIVEN) || !(evalParameters->flags & EVAL_END_GIVEN)) {
1585 double min, max;
1586 find_min_max(&min, &max, x, points);
1587 if (!(evalParameters->flags & EVAL_BEGIN_GIVEN))
1588 evalParameters->begin = min;
1589 if (!(evalParameters->flags & EVAL_END_GIVEN))
1590 evalParameters->end = max;
1591 }
1592 if (!(evalParameters->flags & EVAL_NUMBER_GIVEN))
1593 evalParameters->number = points;
1594 if (evalParameters->number > 1)
1595 delta = (evalParameters->end - evalParameters->begin) / (evalParameters->number - 1);
1596 else
1597 delta = 0;
1598
1599 if (!xEval || maxEvals < evalParameters->number) {
1600 if (!(xEval = (double *)SDDS_Realloc(xEval, sizeof(*xEval) * evalParameters->number)) ||
1601 !(yEval = (double *)SDDS_Realloc(yEval, sizeof(*yEval) * evalParameters->number)))
1602 SDDS_Bomb("allocation failure");
1603 maxEvals = evalParameters->number;
1604 }
1605
1606 for (i = 0; i < evalParameters->number; i++) {
1607 xEval[i] = evalParameters->begin + i * delta;
1608 yEval[i] = eval_sum(basis_fn, coef, order, terms, xEval[i]);
1609 }
1610
1611 if ((iYName == 0 &&
1612 !SDDS_StartPage(&evalParameters->dataset, evalParameters->number)) ||
1613 !SDDS_SetColumnFromDoubles(&evalParameters->dataset, SDDS_SET_BY_NAME, xEval, evalParameters->number, xName) ||
1614 !SDDS_SetColumnFromDoubles(&evalParameters->dataset, SDDS_SET_BY_NAME, yEval, evalParameters->number, yName[iYName]) ||
1615 (iYName == yNames - 1 && !SDDS_WritePage(&evalParameters->dataset)))
1616 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1617}
1618
1619void compareOriginalToFit(double *x, double *y, double **residual, int64_t points, double *rmsResidual, double *coef, int32_t *order, long terms) {
1620 int64_t i;
1621 double residualSum2, fit;
1622
1623 *residual = tmalloc(sizeof(**residual) * points);
1624
1625 residualSum2 = 0;
1626 for (i = 0; i < points; i++) {
1627 fit = eval_sum(basis_fn, coef, order, terms, x[i]);
1628 (*residual)[i] = y[i] - fit;
1629 residualSum2 += sqr((*residual)[i]);
1630 }
1631 *rmsResidual = sqrt(residualSum2 / points);
1632}
SDDS.h
SDDS (Self Describing Data Set) Data Types Definitions and Function Prototypes.
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_SetColumnFromDoubles
int32_t SDDS_SetColumnFromDoubles(SDDS_DATASET *SDDS_dataset, int32_t mode, double *data, int64_t rows,...)
Sets the values for a single data column using double-precision floating-point numbers.
Definition SDDS_dataprep.c:1642
SDDS_SetColumnFromLongs
int32_t SDDS_SetColumnFromLongs(SDDS_DATASET *SDDS_dataset, int32_t mode, int32_t *data, int64_t rows,...)
Sets the values for a single data column using long integer numbers.
Definition SDDS_dataprep.c:2045
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_SetColumnsOfInterest
int32_t SDDS_SetColumnsOfInterest(SDDS_DATASET *SDDS_dataset, int32_t mode,...)
Sets the acceptance flags for columns based on specified naming criteria.
Definition SDDS_extract.c:432
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:49
SDDS_ReadPage
int32_t SDDS_ReadPage(SDDS_DATASET *SDDS_dataset)
Definition SDDS_input.c:1006
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_WritePage
int32_t SDDS_WritePage(SDDS_DATASET *SDDS_dataset)
Writes the current data table to the output file.
Definition SDDS_output.c:1228
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_CheckColumn
int32_t SDDS_CheckColumn(SDDS_DATASET *SDDS_dataset, char *name, char *units, int32_t type, FILE *fp_message)
Checks if a column exists in the SDDS dataset with the specified name, units, and type.
Definition SDDS_utils.c:4562
SDDS_GetColumnNames
char ** SDDS_GetColumnNames(SDDS_DATASET *SDDS_dataset, int32_t *number)
Retrieves the names of all columns in the SDDS dataset.
Definition SDDS_utils.c:2460
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_NumberOfErrors
int32_t SDDS_NumberOfErrors()
Retrieves the number of errors recorded by SDDS library routines.
Definition SDDS_utils.c:304
SDDS_StringIsBlank
int32_t SDDS_StringIsBlank(char *s)
Checks if a string is blank (contains only whitespace characters).
Definition SDDS_utils.c:2404
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_CopyString
int32_t SDDS_CopyString(char **target, const char *source)
Copies a source string to a target string with memory allocation.
Definition SDDS_utils.c:856
SDDS_Realloc
void * SDDS_Realloc(void *old_ptr, size_t new_size)
Reallocates memory to a new size.
Definition SDDS_utils.c:677
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
SDDS_CHARACTER
#define SDDS_CHARACTER
Identifier for the character data type.
Definition SDDStypes.h:91
SDDS_ANY_NUMERIC_TYPE
#define SDDS_ANY_NUMERIC_TYPE
Special identifier used by SDDS_Check*() routines to accept any numeric type.
Definition SDDStypes.h:157
SDDS_DOUBLE
#define SDDS_DOUBLE
Identifier for the double data type.
Definition SDDStypes.h:37
SDDS_LONG64
#define SDDS_LONG64
Identifier for the signed 64-bit integer data type.
Definition SDDStypes.h:49
tmalloc
void * tmalloc(uint64_t size_of_block)
Allocates a memory block of the specified size with zero initialization.
Definition array.c:59
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
eval_sum
double eval_sum(double(*fn)(double x, long ord), double *coef, int32_t *order, long n_coefs, double x0)
Evaluate a sum of basis functions.
Definition lsfBasisFns.c:145
lsfg
long lsfg(double *xd, double *yd, double *sy, long n_pts, long n_terms, int32_t *order, double *coef, double *s_coef, double *chi, double *diff, double(*fn)(double x, long ord))
Computes generalized least squares fits using a function passed by the caller.
Definition lsfg.c:30
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
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
set_argument_scale
void set_argument_scale(double scale)
Set the scale factor applied to the input argument of basis functions.
Definition lsfBasisFns.c:43
dtcheby
double dtcheby(double x, long n)
Evaluate the derivative of the Chebyshev polynomial T_n(x).
Definition lsfBasisFns.c:92
set_argument_offset
void set_argument_offset(double offset)
Set the offset applied to the input argument of basis functions.
Definition lsfBasisFns.c:33
ipower
double ipower(double x, long n)
Evaluate a power function x^n.
Definition lsfBasisFns.c:113
dipower
double dipower(double x, long n)
Evaluate the derivative of x^n.
Definition lsfBasisFns.c:127
tcheby
double tcheby(double x, long n)
Evaluate the Chebyshev polynomial of the first kind T_n(x).
Definition lsfBasisFns.c:72
ChiSqrSigLevel
double ChiSqrSigLevel(double ChiSquared0, long nu)
Computes the probability that a chi-squared variable exceeds a given value.
Definition sigLevel.c:64
str_tolower
char * str_tolower(char *s)
Convert a string to lower case.
Definition str_tolower.c:27
EVAL_PARAMETERS
Definition sddsmpfit.c:286
SDDS_DATASET
Definition SDDS.h:308
Generated by doxygen 1.12.0