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