SDDSlib
Loading...
Searching...
No Matches
sddspfit.c
Go to the documentation of this file.
1/**
2 * @file sddspfit.c
3 * @brief Performs nth order polynomial least squares fitting for SDDS files.
4 *
5 * @details
6 * sddspfit fits data to the form:
7 *
8 * \f[
9 * y = \sum_{i}{ A[i] \cdot P(x - x_{\text{offset}}, i) }
10 * \f]
11 *
12 * where \‍( P(x, i) \‍) is the ith basis function evaluated at \‍( x \‍).
13 * By default, \‍( P(x, i) = x^i \‍), but Chebyshev T polynomials can also be used as the basis functions.
14 *
15 * The program outputs the coefficients \‍( A[i] \‍) and their estimated errors.
16 *
17 * ### Usage
18 *
19 * \code{.sh}
20 * sddspfit [<inputfile>] [<outputfile>] [-pipe=[input][,output]]
21 * -columns=<xname>,<yname>[,xSigma=<name>][,ySigma=<name>]
22 * [ {-terms=<number> [-symmetry={none|odd|even}] | -orders=<number>[,...]} ]
23 * [-reviseOrders [=threshold=<chiValue>] [,verbose] [,complete=<chiThreshold>] [,goodEnough=<chiValue>]]
24 * [-chebyshev [=convert]]
25 * [-xOffset=<value>] [-autoOffset] [-xFactor=<value>]
26 * [-sigmas=<value>,{absolute|fractional}]
27 * [-modifySigmas] [-generateSigmas[={keepLargest|keepSmallest}]]
28 * [-sparse=<interval>] [-range=<lower>,<upper>[,fitOnly]]
29 * [-normalize[=<termNumber>]] [-verbose]
30 * [-evaluate=<filename>[,begin=<value>] [,end=<value>] [,number=<integer>]
31 * [,valuesFile=<filename>,valuesColumn=<string>[,reusePage]]]
32 * [-fitLabelFormat=<sprintf-string>] [-copyParameters] [-majorOrder={row|column}]
33 * \endcode
34 *
35 * @copyright
36 * - (c) 2002 The University of Chicago, as Operator of Argonne National Laboratory.
37 * - (c) 2002 The Regents of the University of California, as Operator of Los Alamos National Laboratory.
38 *
39 * @license
40 * This file is distributed under the terms of the Software License Agreement
41 * found in the file LICENSE included with this distribution.
42 *
43 * @author M. Borland, C. Saunders, R. Soliday, H. Shang
44 */
45
46#include "mdb.h"
47#include "SDDS.h"
48#include "scan.h"
49#include "../../2d_interpolate/nn/nan.h"
50
51void print_coefs(FILE *fprec, double x_offset, double x_scale, long chebyshev,
52 double *coef, double *s_coef, int32_t *order, long n_terms,
53 double chi, long norm_term, char *prepend);
54char **makeCoefficientUnits(SDDS_DATASET *SDDSout, char *xName, char *yName,
55 int32_t *order, long terms);
56long setCoefficientData(SDDS_DATASET *SDDSout, double *coef, double *coefSigma,
57 char **coefUnits, int32_t *order, long terms, long chebyshev,
58 char *fitLabelFormat, char *rpnSeqBuffer);
59char **initializeOutputFile(SDDS_DATASET *SDDSout, char *output,
60 SDDS_DATASET *SDDSin, char *input, char *xName,
61 char *yName, char *xSigmaName, char *ySigmaName,
62 long sigmasValid, int32_t *order, long terms,
63 long chebyshev, long copyParameters);
64void checkInputFile(SDDS_DATASET *SDDSin, char *xName, char *yName,
65 char *xSigmaName, char *ySigmaName);
66long coefficient_index(int32_t *order, long terms, long order_of_interest);
67void makeFitLabel(char *buffer, long bufsize, char *fitLabelFormat,
68 double *coef, int32_t *order, long terms, long chebyshev);
69void createRpnSequence(char *buffer, long bufsize, double *coef, int32_t *order,
70 long terms);
71
72double tcheby(double x, long n);
73double dtcheby(double x, long n);
74double ipower(double x, long n);
75double dipower(double x, long n);
76long reviseFitOrders(double *x, double *y, double *sy, int64_t points,
77 long terms, int32_t *order, double *coef,
78 double *coefSigma, double *diff,
79 double (*basis_fn)(double xa, long ordera),
80 unsigned long reviseOrders, double xOffset,
81 double xScaleFactor, long normTerm, long ySigmasValid,
82 long chebyshev, double revpowThreshold,
83 double revpowCompleteThres, double goodEnoughChi);
84long reviseFitOrders1(double *x, double *y, double *sy, int64_t points,
85 long terms, int32_t *order, double *coef,
86 double *coefSigma, double *diff,
87 double (*basis_fn)(double xa, long ordera),
88 unsigned long reviseOrders, double xOffset,
89 double xScaleFactor, long normTerm, long ySigmasValid,
90 long chebyshev, double revpowThreshold,
91 double goodEnoughChi);
92void compareOriginalToFit(double *x, double *y, double **residual,
93 int64_t points, double *rmsResidual, double *coef,
94 int32_t *order, long terms);
95
96typedef struct {
97 long nTerms; /* Maximum power is nTerms-1 */
98 double *coef; /* Coefficients for Chebyshev T polynomials */
99} CHEBYSHEV_COEF;
100
101CHEBYSHEV_COEF *makeChebyshevCoefficients(long maxOrder, long *nPoly);
102void convertFromChebyshev(long termsT, int32_t *orderT, double *coefT, double *coefSigmaT,
103 long *termsOrdinaryRet, int32_t **orderOrdinaryRet, double **coefOrdinaryRet, double **coefSigmaOrdinaryRet);
104
105/* Enumeration for option types */
106enum option_type {
107 CLO_COLUMNS,
108 CLO_ORDERS,
109 CLO_TERMS,
110 CLO_SYMMETRY,
111 CLO_REVISEORDERS,
112 CLO_CHEBYSHEV,
113 CLO_MODIFYSIGMAS,
114 CLO_SIGMAS,
115 CLO_GENERATESIGMAS,
116 CLO_RANGE,
117 CLO_SPARSE,
118 CLO_NORMALIZE,
119 CLO_XFACTOR,
120 CLO_XOFFSET,
121 CLO_VERBOSE,
122 CLO_FITLABELFORMAT,
123 CLO_PIPE,
124 CLO_EVALUATE,
125 CLO_AUTOOFFSET,
126 CLO_COPY_PARAMETERS,
127 CLO_MAJOR_ORDER,
128 N_OPTIONS
129};
130
131char *option[N_OPTIONS] = {
132 "columns",
133 "orders",
134 "terms",
135 "symmetry",
136 "reviseorders",
137 "chebyshev",
138 "modifysigmas",
139 "sigmas",
140 "generatesigmas",
141 "range",
142 "sparse",
143 "normalize",
144 "xfactor",
145 "xoffset",
146 "verbose",
147 "fitlabelformat",
148 "pipe",
149 "evaluate",
150 "autooffset",
151 "copyparameters",
152 "majorOrder",
153};
154
155char *USAGE =
156 "sddspfit [<inputfile>] [<outputfile>] [-pipe=[input][,output]]\n\
157 -columns=<xname>,<yname>[,xSigma=<name>][,ySigma=<name>]\n\
158 [ {-terms=<number> [-symmetry={none|odd|even}] | -orders=<number>[,...]} ]\n\
159 [-reviseOrders [=threshold=<chiValue>] [,verbose] [,complete=<chiThreshold>] [,goodEnough=<chiValue>]]\n\
160 [-chebyshev [=convert]]\n\
161 [-xOffset=<value>] [-autoOffset] [-xFactor=<value>]\n\
162 [-sigmas=<value>,{absolute|fractional}] \n\
163 [-modifySigmas] [-generateSigmas[={keepLargest|keepSmallest}]]\n\
164 [-sparse=<interval>] [-range=<lower>,<upper>[,fitOnly]]\n\
165 [-normalize[=<termNumber>]] [-verbose]\n\
166 [-evaluate=<filename>[,begin=<value>] [,end=<value>] [,number=<integer>] \n\
167 [,valuesFile=<filename>,valuesColumn=<string>[,reusePage]]]\n\
168 [-fitLabelFormat=<sprintf-string>] [-copyParameters] [-majorOrder={row|column}]\n\n\
169Program by Michael Borland. (" __DATE__ " " __TIME__ ", SVN revision: " SVN_VERSION ")\n";
170
171static char *additional_help1 = "\n\
172sddspfit fits data to the form y = SUM(i){ A[i] * P(x - x_offset, i) }, where P(x, i) is the ith basis\n\
173function evaluated at x. By default, P(x, i) = x^i. Chebyshev T polynomials can also be selected as the basis functions.\n\n\
174 -columns Specify names of data columns to use.\n\
175 -terms Number of terms desired in fit.\n\
176 -symmetry Symmetry of desired fit about x_offset.\n\
177 -orders Orders (P[i]) to use in fitting.\n\
178 -reviseOrders Modify the orders used in the fit to eliminate poorly-determined coefficients based on fitting\n\
179 of the first data page. The algorithm adds one order at a time, terminating when the reduced\n\
180 chi-squared is less than the \"goodEnough\" value (default: 1) or when the new term does not improve\n\
181 the reduced chi-squared by more than the threshold value (default: 0.1). It next tries removing terms one at a time.\n\
182 Finally, if the resulting best reduced chi-squared is greater than the threshold given with the \"complete\" option,\n\
183 it also tries all possible combinations of allowed terms.\n\
184 -chebyshev Use Chebyshev T polynomials (xOffset is set automatically).\n\
185 Giving the `convert` option causes the fit to be written out in terms of ordinary polynomials.\n\
186 -majorOrder Specify output file in row or column major order.\n\
187 -xOffset Desired value of x to fit about.\n";
188
189static char *additional_help2 =
190 " -autoOffset Automatically offset x values by the mean x value for fitting.\n\
191 Helpful if x values are very large in magnitude.\n\
192 -xFactor Desired factor to multiply x values by before fitting.\n\
193 -sigmas Specify absolute or fractional sigma for all points.\n\
194 -modifySigmas Modify the y sigmas using the x sigmas and an initial fit.\n\
195 -generateSigmas Generate y sigmas from the RMS deviation from an initial fit.\n\
196 Optionally keep the sigmas from the data if larger/smaller than RMS deviation.\n\
197 -sparse Specify integer interval at which to sample data.\n\
198 -range Specify range of independent variable over which to perform fit and evaluation.\n\
199 If 'fitOnly' is given, then fit is compared to data over the original range.\n\
200 -normalize Normalize so that the specified term is unity.\n\
201 -verbose Generates extra output that may be useful.\n\
202 -evaluate Specify evaluation of fit over a selected range of equispaced points,\n\
203 or at values listed in a file.\n\
204 -copyParameters If given, program copies all parameters from the input file into the main output file.\n\
205 By default, no parameters are copied.\n\n";
206
207#define NO_SYMMETRY 0
208#define EVEN_SYMMETRY 1
209#define ODD_SYMMETRY 2
210#define N_SYMMETRY_OPTIONS 3
211char *symmetry_options[N_SYMMETRY_OPTIONS] = {"none", "even", "odd"};
212
213#define ABSOLUTE_SIGMAS 0
214#define FRACTIONAL_SIGMAS 1
215#define N_SIGMAS_OPTIONS 2
216char *sigmas_options[N_SIGMAS_OPTIONS] = {"absolute", "fractional"};
217
218#define FLGS_GENERATESIGMAS 1
219#define FLGS_KEEPLARGEST 2
220#define FLGS_KEEPSMALLEST 4
221
222#define REVPOW_ACTIVE 0x0001
223#define REVPOW_VERBOSE 0x0002
224#define REVPOW_COMPLETE 0x0004
225
226/* SDDS indices for output page */
227static long iIntercept = -1, iInterceptSigma = -1;
228static long iSlope = -1, iSlopeSigma = -1;
229static long iCurvature = -1, iCurvatureSigma = -1;
230static long *iTerm = NULL, *iTermSig = NULL;
231static long iOffset = -1, iFactor = -1;
232static long iChiSq = -1, iRmsResidual = -1, iSigLevel = -1;
233static long iFitIsValid = -1, iFitLabel = -1, iTerms = -1;
234static long iRpnSequence;
235
236static long ix = -1, iy = -1, ixSigma = -1, iySigma = -1;
237static long iFit = -1, iResidual = -1;
238
239static char *xSymbol, *ySymbol;
240
241#define EVAL_BEGIN_GIVEN 0x0001U
242#define EVAL_END_GIVEN 0x0002U
243#define EVAL_NUMBER_GIVEN 0x0004U
244#define EVAL_VALUESFILE_GIVEN 0x0008U
245#define EVAL_VALUESCOLUMN_GIVEN 0x0010U
246#define EVAL_REUSE_PAGE_GIVEN 0x0020U
247
248typedef struct {
249 char *file;
250 short initialized;
251 int64_t number;
252 unsigned long flags;
253 double begin, end;
254 SDDS_DATASET dataset;
255 short inputInitialized;
256 char *valuesFile, *valuesColumn;
257 SDDS_DATASET valuesDataset;
258} EVAL_PARAMETERS;
259void makeEvaluationTable(EVAL_PARAMETERS *evalParameters, double *x, int64_t n,
260 double *coef, int32_t *order, long terms,
261 SDDS_DATASET *SDDSin, char *xName, char *yName);
262
263static double (*basis_fn)(double xa, long ordera);
264static double (*basis_dfn)(double xa, long ordera);
265
266int main(int argc, char **argv) {
267 double *x = NULL, *y = NULL, *sy = NULL, *sx = NULL, *diff = NULL, xOffset,
268 xScaleFactor;
269 double *xOrig = NULL, *yOrig = NULL, *sxOrig, *syOrig, *sy0 = NULL;
270 long terms, normTerm, ySigmasValid;
271 int64_t i, j, points, pointsOrig;
272 long symmetry, chebyshev, autoOffset, copyParameters = 0;
273 double sigmas;
274 long sigmasMode, sparseInterval;
275 unsigned long flags;
276 double *coef, *coefSigma;
277 double chi, xLow, xHigh, rmsResidual;
278 char *xName, *yName, *xSigmaName, *ySigmaName;
279 char *input, *output, **coefUnits;
280 SDDS_DATASET SDDSin, SDDSout;
281 long isFit, iArg, modifySigmas;
282 long generateSigmas, verbose, ignoreSigmas;
283 long npages = 0, invalid = 0;
284 int32_t *order;
285 SCANNED_ARG *s_arg;
286 double xMin, xMax, revpowThreshold, revpowCompleteThres, goodEnoughChi;
287 long rangeFitOnly = 0;
288 double rms_average(double *d_x, int64_t d_n);
289 char *fitLabelFormat = "%g";
290 static char rpnSeqBuffer[SDDS_MAXLINE];
291 unsigned long pipeFlags, reviseOrders, majorOrderFlag;
292 EVAL_PARAMETERS evalParameters;
293 short columnMajorOrder = -1;
294
295 sxOrig = syOrig = NULL;
296 rmsResidual = 0;
297
298 SDDS_RegisterProgramName(argv[0]);
299 argc = scanargs(&s_arg, argc, argv);
300 if (argc < 2 || argc > (3 + N_OPTIONS)) {
301 fprintf(stderr, "usage: %s%s%s\n", USAGE, additional_help1,
302 additional_help2);
303 exit(EXIT_FAILURE);
304 }
305
306 input = output = NULL;
307 xName = yName = xSigmaName = ySigmaName = NULL;
308 modifySigmas = reviseOrders = chebyshev = 0;
309 order = NULL;
310 symmetry = NO_SYMMETRY;
311 xMin = xMax = 0;
312 autoOffset = 0;
313 generateSigmas = 0;
314 sigmasMode = -1;
315 sigmas = 1;
316 sparseInterval = 1;
317 terms = 2;
318 verbose = ignoreSigmas = 0;
319 normTerm = -1;
320 xOffset = 0;
321 xScaleFactor = 1;
322 coefUnits = NULL;
323 basis_fn = ipower;
324 basis_dfn = dipower;
325 pipeFlags = 0;
326 evalParameters.file = evalParameters.valuesFile = evalParameters.valuesColumn = NULL;
327 evalParameters.initialized = evalParameters.inputInitialized = 0;
328
329 for (iArg = 1; iArg < argc; iArg++) {
330 if (s_arg[iArg].arg_type == OPTION) {
331 switch (match_string(s_arg[iArg].list[0], option, N_OPTIONS, 0)) {
332 case CLO_MAJOR_ORDER:
333 majorOrderFlag = 0;
334 s_arg[iArg].n_items--;
335 if (s_arg[iArg].n_items > 0 &&
336 (!scanItemList(&majorOrderFlag, s_arg[iArg].list + 1,
337 &s_arg[iArg].n_items, 0, "row", -1, NULL, 0,
338 SDDS_ROW_MAJOR_ORDER, "column", -1, NULL, 0,
339 SDDS_COLUMN_MAJOR_ORDER, NULL)))
340 SDDS_Bomb("invalid -majorOrder syntax/values");
341 if (majorOrderFlag & SDDS_COLUMN_MAJOR_ORDER)
342 columnMajorOrder = 1;
343 else if (majorOrderFlag & SDDS_ROW_MAJOR_ORDER)
344 columnMajorOrder = 0;
345 break;
346 case CLO_MODIFYSIGMAS:
347 modifySigmas = 1;
348 break;
349 case CLO_AUTOOFFSET:
350 autoOffset = 1;
351 break;
352 case CLO_ORDERS:
353 if (s_arg[iArg].n_items < 2)
354 SDDS_Bomb("invalid -orders syntax");
355 order = tmalloc(sizeof(*order) * (terms = s_arg[iArg].n_items - 1));
356 for (i = 1; i < s_arg[iArg].n_items; i++) {
357 if (sscanf(s_arg[iArg].list[i], "%" SCNd32, order + i - 1) != 1)
358 SDDS_Bomb("unable to scan order from -orders list");
359 }
360 break;
361 case CLO_RANGE:
362 rangeFitOnly = 0;
363 if ((s_arg[iArg].n_items != 3 && s_arg[iArg].n_items != 4) ||
364 1 != sscanf(s_arg[iArg].list[1], "%lf", &xMin) ||
365 1 != sscanf(s_arg[iArg].list[2], "%lf", &xMax) || xMin >= xMax)
366 SDDS_Bomb("incorrect -range syntax");
367 if (s_arg[iArg].n_items == 4) {
368 if (strncmp(str_tolower(s_arg[iArg].list[3]), "fitonly",
369 strlen(s_arg[iArg].list[3])) == 0) {
370 rangeFitOnly = 1;
371 } else
372 SDDS_Bomb("incorrect -range syntax");
373 }
374 break;
375 case CLO_GENERATESIGMAS:
376 generateSigmas = FLGS_GENERATESIGMAS;
377 if (s_arg[iArg].n_items > 1) {
378 if (s_arg[iArg].n_items != 2)
379 SDDS_Bomb("incorrect -generateSigmas syntax");
380 if (strncmp(s_arg[iArg].list[1], "keepsmallest",
381 strlen(s_arg[iArg].list[1])) == 0)
382 generateSigmas |= FLGS_KEEPSMALLEST;
383 if (strncmp(s_arg[iArg].list[1], "keeplargest",
384 strlen(s_arg[iArg].list[1])) == 0)
385 generateSigmas |= FLGS_KEEPLARGEST;
386 if ((generateSigmas & FLGS_KEEPSMALLEST) &&
387 (generateSigmas & FLGS_KEEPLARGEST))
388 SDDS_Bomb("ambiguous -generateSigmas syntax");
389 }
390 break;
391 case CLO_TERMS:
392 if (s_arg[iArg].n_items != 2 ||
393 sscanf(s_arg[iArg].list[1], "%ld", &terms) != 1)
394 SDDS_Bomb("invalid -terms syntax");
395 break;
396 case CLO_XOFFSET:
397 if (s_arg[iArg].n_items != 2 ||
398 sscanf(s_arg[iArg].list[1], "%lf", &xOffset) != 1)
399 SDDS_Bomb("invalid -xOffset syntax");
400 break;
401 case CLO_SYMMETRY:
402 if (s_arg[iArg].n_items == 2) {
403 if ((symmetry = match_string(s_arg[iArg].list[1], symmetry_options,
404 N_SYMMETRY_OPTIONS, 0)) < 0)
405 SDDS_Bomb("unknown option used with -symmetry");
406 } else
407 SDDS_Bomb("incorrect -symmetry syntax");
408 break;
409 case CLO_SIGMAS:
410 if (s_arg[iArg].n_items != 3)
411 SDDS_Bomb("incorrect -sigmas syntax");
412 if (sscanf(s_arg[iArg].list[1], "%lf", &sigmas) != 1)
413 SDDS_Bomb("couldn't scan value for -sigmas");
414 if ((sigmasMode = match_string(s_arg[iArg].list[2], sigmas_options,
415 N_SIGMAS_OPTIONS, 0)) < 0)
416 SDDS_Bomb("unrecognized -sigmas mode");
417 break;
418 case CLO_SPARSE:
419 if (s_arg[iArg].n_items != 2)
420 SDDS_Bomb("incorrect -sparse syntax");
421 if (sscanf(s_arg[iArg].list[1], "%ld", &sparseInterval) != 1)
422 SDDS_Bomb("couldn't scan value for -sparse");
423 if (sparseInterval < 1)
424 SDDS_Bomb("invalid -sparse value");
425 break;
426 case CLO_VERBOSE:
427 verbose = 1;
428 break;
429 case CLO_NORMALIZE:
430 normTerm = 0;
431 if (s_arg[iArg].n_items > 2 ||
432 (s_arg[iArg].n_items == 2 &&
433 sscanf(s_arg[iArg].list[1], "%ld", &normTerm) != 1) ||
434 normTerm < 0)
435 SDDS_Bomb("invalid -normalize syntax");
436 break;
437 case CLO_REVISEORDERS:
438 revpowThreshold = 0.1;
439 revpowCompleteThres = 10;
440 goodEnoughChi = 1;
441 s_arg[iArg].n_items -= 1;
442 if (!scanItemList(&reviseOrders, s_arg[iArg].list + 1,
443 &s_arg[iArg].n_items, 0,
444 "threshold", SDDS_DOUBLE, &revpowThreshold, 1, 0,
445 "complete", SDDS_DOUBLE, &revpowCompleteThres, 1, REVPOW_COMPLETE,
446 "goodenough", SDDS_DOUBLE, &goodEnoughChi, 1, 0,
447 "verbose", -1, NULL, 1, REVPOW_VERBOSE, NULL) ||
448 revpowThreshold < 0 || revpowCompleteThres < 0 || goodEnoughChi < 0)
449 SDDS_Bomb("invalid -reviseOrders syntax");
450 reviseOrders |= REVPOW_ACTIVE;
451 break;
452 case CLO_CHEBYSHEV:
453 if (s_arg[iArg].n_items > 2 ||
454 (s_arg[iArg].n_items == 2 &&
455 strncmp(s_arg[iArg].list[1], "convert",
456 strlen(s_arg[iArg].list[1])) != 0))
457 SDDS_Bomb("invalid -chebyshev syntax");
458 chebyshev = s_arg[iArg].n_items; /* 1: use chebyshev polynomials; 2: also convert to ordinary form */
459 basis_fn = tcheby;
460 basis_dfn = dtcheby;
461 break;
462 case CLO_XFACTOR:
463 if (s_arg[iArg].n_items != 2 ||
464 sscanf(s_arg[iArg].list[1], "%lf", &xScaleFactor) != 1 ||
465 xScaleFactor == 0)
466 SDDS_Bomb("invalid -xFactor syntax");
467 break;
468 case CLO_COLUMNS:
469 if (s_arg[iArg].n_items < 3 || s_arg[iArg].n_items > 5)
470 SDDS_Bomb("invalid -columns syntax");
471 xName = s_arg[iArg].list[1];
472 yName = s_arg[iArg].list[2];
473 s_arg[iArg].n_items -= 3;
474 if (!scanItemList(&flags, s_arg[iArg].list + 3, &s_arg[iArg].n_items, 0,
475 "xsigma", SDDS_STRING, &xSigmaName, 1, 0, "ysigma",
476 SDDS_STRING, &ySigmaName, 1, 0, NULL))
477 SDDS_Bomb("invalid -columns syntax");
478 break;
479 case CLO_FITLABELFORMAT:
480 if (s_arg[iArg].n_items != 2)
481 SDDS_Bomb("invalid -fitLabelFormat syntax");
482 fitLabelFormat = s_arg[iArg].list[1];
483 break;
484 case CLO_PIPE:
485 if (!processPipeOption(s_arg[iArg].list + 1, s_arg[iArg].n_items - 1,
486 &pipeFlags))
487 SDDS_Bomb("invalid -pipe syntax");
488 break;
489 case CLO_EVALUATE:
490 if (s_arg[iArg].n_items < 2)
491 SDDS_Bomb("invalid -evaluate syntax");
492 evalParameters.file = s_arg[iArg].list[1];
493 s_arg[iArg].n_items -= 2;
494 s_arg[iArg].list += 2;
495 if (!scanItemList(&evalParameters.flags, s_arg[iArg].list,
496 &s_arg[iArg].n_items, 0,
497 "begin", SDDS_DOUBLE, &evalParameters.begin, 1, EVAL_BEGIN_GIVEN,
498 "end", SDDS_DOUBLE, &evalParameters.end, 1, EVAL_END_GIVEN,
499 "number", SDDS_LONG64, &evalParameters.number, 1, EVAL_NUMBER_GIVEN,
500 "valuesfile", SDDS_STRING, &evalParameters.valuesFile, 1, EVAL_VALUESFILE_GIVEN,
501 "valuescolumn", SDDS_STRING, &evalParameters.valuesColumn, 1, EVAL_VALUESCOLUMN_GIVEN,
502 "reusepage", 0, NULL, 0, EVAL_REUSE_PAGE_GIVEN,
503 NULL))
504 SDDS_Bomb("invalid -evaluate syntax");
505 if (evalParameters.flags & EVAL_VALUESFILE_GIVEN || evalParameters.flags & EVAL_VALUESCOLUMN_GIVEN) {
506 if (evalParameters.flags & (EVAL_BEGIN_GIVEN | EVAL_END_GIVEN | EVAL_NUMBER_GIVEN))
507 SDDS_Bomb("invalid -evaluate syntax: given begin/end/number or valuesFile/valuesColumn, not a mixture.");
508 if (!(evalParameters.flags & EVAL_VALUESFILE_GIVEN && evalParameters.flags & EVAL_VALUESCOLUMN_GIVEN))
509 SDDS_Bomb("invalid -evaluate syntax: give both valuesFile and valuesColumn, not just one");
510 }
511 evalParameters.initialized = 0;
512 break;
513 case CLO_COPY_PARAMETERS:
514 copyParameters = 1;
515 break;
516 default:
517 bomb("unknown switch", USAGE);
518 break;
519 }
520 } else {
521 if (input == NULL)
522 input = s_arg[iArg].list[0];
523 else if (output == NULL)
524 output = s_arg[iArg].list[0];
525 else
526 SDDS_Bomb("too many filenames");
527 }
528 }
529
530 processFilenames("sddspfit", &input, &output, pipeFlags, 0, NULL);
531
532 if (symmetry && order)
533 SDDS_Bomb("can't specify both -symmetry and -orders");
534 if (chebyshev && order)
535 SDDS_Bomb("can't specify both -chebyshev and -orders");
536 if (chebyshev && symmetry)
537 SDDS_Bomb("can't specify both -chebyshev and -symmetry");
538 if (!xName || !yName)
539 SDDS_Bomb("you must specify a column name for x and y");
540
541 if (modifySigmas && !xSigmaName)
542 SDDS_Bomb("you must specify x sigmas with -modifySigmas");
543 if (generateSigmas) {
544 if (modifySigmas)
545 SDDS_Bomb("you can't specify both -generateSigmas and -modifySigmas");
546 }
547 if (ySigmaName) {
548 if (sigmasMode != -1)
549 SDDS_Bomb("you can't specify both -sigmas and a y sigma name");
550 }
551 ySigmasValid = 0;
552 if (sigmasMode != -1 || generateSigmas || ySigmaName || modifySigmas)
553 ySigmasValid = 1;
554
555 if (normTerm >= 0 && normTerm >= terms)
556 SDDS_Bomb("can't normalize to that term--not that many terms");
557 if (reviseOrders && !(sigmasMode != -1 || generateSigmas || ySigmaName))
558 SDDS_Bomb(
559 "can't use -reviseOrders unless a y sigma or -generateSigmas is given");
560
561 if (symmetry == EVEN_SYMMETRY) {
562 order = tmalloc(sizeof(*order) * terms);
563 for (i = 0; i < terms; i++)
564 order[i] = 2 * i;
565 } else if (symmetry == ODD_SYMMETRY) {
566 order = tmalloc(sizeof(*order) * terms);
567 for (i = 0; i < terms; i++)
568 order[i] = 2 * i + 1;
569 } else if (!order) {
570 order = tmalloc(sizeof(*order) * terms);
571 for (i = 0; i < terms; i++)
572 order[i] = i;
573 }
574 coef = tmalloc(sizeof(*coef) * terms);
575 coefSigma = tmalloc(sizeof(*coefSigma) * terms);
576 iTerm = tmalloc(sizeof(*iTerm) * terms);
577 iTermSig = tmalloc(sizeof(*iTermSig) * terms);
578
579 if (!SDDS_InitializeInput(&SDDSin, input))
580 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
581 checkInputFile(&SDDSin, xName, yName, xSigmaName, ySigmaName);
582 coefUnits = initializeOutputFile(&SDDSout, output, &SDDSin, input, xName,
583 yName, xSigmaName, ySigmaName, ySigmasValid,
584 order, terms, chebyshev, copyParameters);
585 if (columnMajorOrder != -1)
586 SDDSout.layout.data_mode.column_major = columnMajorOrder;
587 else
588 SDDSout.layout.data_mode.column_major =
589 SDDSin.layout.data_mode.column_major;
590 while (SDDS_ReadPage(&SDDSin) > 0) {
591 npages++;
592 invalid = 0;
593 if ((points = SDDS_CountRowsOfInterest(&SDDSin)) < terms) {
594 pointsOrig = 0;
595 invalid = 1;
596 isFit = 0;
597 } else {
598 if (!(x = SDDS_GetColumnInDoubles(&SDDSin, xName))) {
599 fprintf(stderr, "error: unable to read column %s\n", xName);
600 SDDS_PrintErrors(stderr,
601 SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
602 }
603 if (!(y = SDDS_GetColumnInDoubles(&SDDSin, yName))) {
604 fprintf(stderr, "error: unable to read column %s\n", yName);
605 SDDS_PrintErrors(stderr,
606 SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
607 }
608 sx = NULL;
609 if (xSigmaName && !(sx = SDDS_GetColumnInDoubles(&SDDSin, xSigmaName))) {
610 fprintf(stderr, "error: unable to read column %s\n", xSigmaName);
611 SDDS_PrintErrors(stderr,
612 SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
613 }
614 sy0 = NULL;
615 if (ySigmaName && !(sy0 = SDDS_GetColumnInDoubles(&SDDSin, ySigmaName))) {
616 fprintf(stderr, "error: unable to read column %s\n", ySigmaName);
617 SDDS_PrintErrors(stderr,
618 SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
619 }
620 if (!sy0)
621 sy0 = tmalloc(sizeof(*sy0) * points);
622
623 if (xMin != xMax || sparseInterval != 1) {
624 xOrig = tmalloc(sizeof(*xOrig) * points);
625 yOrig = tmalloc(sizeof(*yOrig) * points);
626 if (sx)
627 sxOrig = tmalloc(sizeof(*sxOrig) * points);
628 if (ySigmasValid)
629 syOrig = tmalloc(sizeof(*syOrig) * points);
630 pointsOrig = points;
631 for (i = j = 0; i < points; i++) {
632 xOrig[i] = x[i];
633 yOrig[i] = y[i];
634 if (ySigmasValid)
635 syOrig[i] = sy0[i];
636 if (sx)
637 sxOrig[i] = sx[i];
638 }
639 if (xMin != xMax) {
640 for (i = j = 0; i < points; i++) {
641 if (xOrig[i] <= xMax && xOrig[i] >= xMin) {
642 x[j] = xOrig[i];
643 y[j] = yOrig[i];
644 if (ySigmasValid)
645 sy0[j] = syOrig[i];
646 if (sx)
647 sx[j] = sxOrig[i];
648 j++;
649 }
650 }
651 points = j;
652 }
653 if (sparseInterval != 1) {
654 for (i = j = 0; i < points; i++) {
655 if (i % sparseInterval == 0) {
656 x[j] = x[i];
657 y[j] = y[i];
658 if (ySigmasValid)
659 sy0[j] = sy0[i];
660 if (sx)
661 sx[j] = sx[i];
662 j++;
663 }
664 }
665 points = j;
666 }
667 } else {
668 xOrig = x;
669 yOrig = y;
670 sxOrig = sx;
671 syOrig = sy0;
672 pointsOrig = points;
673 }
674
675 find_min_max(&xLow, &xHigh, x, points);
676
677 if (sigmasMode == ABSOLUTE_SIGMAS) {
678 for (i = 0; i < points; i++)
679 sy0[i] = sigmas;
680 if (sy0 != syOrig)
681 for (i = 0; i < pointsOrig; i++)
682 syOrig[i] = sigmas;
683 } else if (sigmasMode == FRACTIONAL_SIGMAS) {
684 for (i = 0; i < points; i++)
685 sy0[i] = sigmas * fabs(y[i]);
686 if (sy0 != syOrig)
687 for (i = 0; i < pointsOrig; i++)
688 syOrig[i] = fabs(yOrig[i]) * sigmas;
689 }
690
691 if (!ySigmasValid || generateSigmas)
692 for (i = 0; i < points; i++)
693 sy0[i] = 1;
694 else
695 for (i = 0; i < points; i++)
696 if (sy0[i] == 0)
697 SDDS_Bomb("y sigma = 0 for one or more points.");
698
699 diff = tmalloc(sizeof(*x) * points);
700 sy = tmalloc(sizeof(*sy) * points);
701 for (i = 0; i < points; i++)
702 sy[i] = sy0[i];
703
704 if (autoOffset && !compute_average(&xOffset, x, points))
705 xOffset = 0;
706
707 set_argument_offset(xOffset);
708 set_argument_scale(xScaleFactor);
709 if (chebyshev) {
710 if (xOffset) {
711 /* User has provided an offset, adjust scale factor to match range of data */
712 xScaleFactor = MAX(fabs(xHigh - xOffset), fabs(xLow - xOffset));
713 } else {
714 /* Compute the offset and scale factor to match range of data */
715 xOffset = (xHigh + xLow) / 2;
716 xScaleFactor = (xHigh - xLow) / 2;
717 }
718 set_argument_offset(xOffset);
719 set_argument_scale(xScaleFactor);
720 }
721
722 if (generateSigmas || modifySigmas) {
723 /* do an initial fit */
724 isFit = lsfg(x, y, sy, points, terms, order, coef, coefSigma, &chi,
725 diff, basis_fn);
726 if (!isFit)
727 SDDS_Bomb("initial fit failed.");
728 if (verbose) {
729 fputs("initial_fit:", stdout);
730 print_coefs(stdout, xOffset, xScaleFactor, chebyshev, coef, NULL,
731 order, terms, chi, normTerm, "");
732 fprintf(stdout, "unweighted rms deviation from fit: %21.15e\n",
733 rms_average(diff, points));
734 }
735 if (modifySigmas) {
736 if (!ySigmasValid) {
737 for (i = 0; i < points; i++)
738 sy[i] =
739 fabs(eval_sum(basis_dfn, coef, order, terms, x[i]) * sx[i]);
740 } else
741 for (i = 0; i < points; i++) {
742 sy[i] = sqrt(
743 sqr(sy0[i]) +
744 sqr(eval_sum(basis_dfn, coef, order, terms, x[i]) * sx[i]));
745 }
746 }
747 if (generateSigmas) {
748 double sigma;
749 for (i = sigma = 0; i < points; i++) {
750 sigma += sqr(diff[i]);
751 }
752 sigma = sqrt(sigma / (points - terms));
753 for (i = 0; i < points; i++) {
754 if (generateSigmas & FLGS_KEEPSMALLEST) {
755 if (sigma < sy[i])
756 sy[i] = sigma;
757 } else if (generateSigmas & FLGS_KEEPLARGEST) {
758 if (sigma > sy[i])
759 sy[i] = sigma;
760 } else {
761 sy[i] = sigma;
762 }
763 }
764 for (i = 0; i < pointsOrig; i++) {
765 if (generateSigmas & FLGS_KEEPSMALLEST) {
766 if (sigma < sy0[i])
767 sy0[i] = sigma;
768 } else if (generateSigmas & FLGS_KEEPLARGEST) {
769 if (sigma > sy0[i])
770 sy0[i] = sigma;
771 } else {
772 sy0[i] = sigma;
773 }
774 }
775 }
776 }
777
778 if (reviseOrders & REVPOW_ACTIVE) {
779 terms = reviseFitOrders(
780 x, y, sy, points, terms, order, coef, coefSigma, diff, basis_fn,
781 reviseOrders, xOffset, xScaleFactor, normTerm, ySigmasValid,
782 chebyshev, revpowThreshold, revpowCompleteThres, goodEnoughChi);
783 reviseOrders = 0;
784 }
785
786 isFit = lsfg(x, y, sy, points, terms, order, coef, coefSigma, &chi, diff,
787 basis_fn);
788 if (isFit) {
789 rmsResidual = rms_average(diff, points);
790 if (verbose) {
791 print_coefs(stdout, xOffset, xScaleFactor, chebyshev, coef,
792 (ySigmasValid ? coefSigma : NULL), order, terms, chi,
793 normTerm, "");
794 fprintf(stdout, "unweighted rms deviation from fit: %21.15e\n",
795 rmsResidual);
796 }
797 } else if (verbose)
798 fprintf(stdout, "fit failed.\n");
799
800 if (evalParameters.file)
801 makeEvaluationTable(&evalParameters, x, points, coef, order, terms,
802 &SDDSin, xName, yName);
803 }
804
805 if (!SDDS_StartPage(&SDDSout, rangeFitOnly ? pointsOrig : points))
806 SDDS_PrintErrors(stderr,
807 SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
808 rpnSeqBuffer[0] = 0;
809 if (!invalid) {
810 setCoefficientData(&SDDSout, coef, (ySigmasValid ? coefSigma : NULL),
811 coefUnits, order, terms, chebyshev, fitLabelFormat,
812 rpnSeqBuffer);
813 if (rangeFitOnly) {
814 double *residual;
815 compareOriginalToFit(xOrig, yOrig, &residual, pointsOrig, &rmsResidual,
816 coef, order, terms);
817
818 if (!SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_INDEX, xOrig,
819 pointsOrig, ix) ||
820 !SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_INDEX, yOrig,
821 pointsOrig, iy) ||
822 !SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_INDEX, residual,
823 pointsOrig, iResidual))
824 SDDS_PrintErrors(stderr,
825 SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
826 for (i = 0; i < pointsOrig; i++)
827 residual[i] = yOrig[i] - residual[i]; /* computes fit values from
828 residual and y */
829 if (!SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_INDEX, residual,
830 pointsOrig, iFit))
831 SDDS_PrintErrors(stderr,
832 SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
833
834 if (ixSigma != -1 &&
835 !SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_INDEX, sxOrig,
836 pointsOrig, ixSigma))
837 SDDS_PrintErrors(stderr,
838 SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
839 if (ySigmasValid && iySigma != -1 &&
840 !SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_INDEX, syOrig,
841 pointsOrig, iySigma))
842 SDDS_PrintErrors(stderr,
843 SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
844 free(residual);
845 } else {
846 for (i = 0; i < points; i++)
847 diff[i] =
848 -diff[i]; /* convert from (Fit-y) to (y-Fit) to get residual */
849 if (!SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_INDEX, x, points,
850 ix) ||
851 !SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_INDEX, y, points,
852 iy) ||
853 !SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_INDEX, diff,
854 points, iResidual))
855 SDDS_PrintErrors(stderr,
856 SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
857 for (i = 0; i < points; i++)
858 diff[i] =
859 y[i] - diff[i]; /* computes fit values from residual and y */
860 if (!SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_INDEX, diff,
861 points, iFit))
862 SDDS_PrintErrors(stderr,
863 SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
864
865 if (ixSigma != -1 &&
866 !SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_INDEX, sx, points,
867 ixSigma))
868 SDDS_PrintErrors(stderr,
869 SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
870 if (ySigmasValid && iySigma != -1 &&
871 !SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_INDEX, sy, points,
872 iySigma))
873 SDDS_PrintErrors(stderr,
874 SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
875 }
876 }
877
878 if (copyParameters && !SDDS_CopyParameters(&SDDSout, &SDDSin))
879 SDDS_PrintErrors(stderr,
880 SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
881 if (!SDDS_SetParameters(
882 &SDDSout, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE, iRpnSequence,
883 invalid ? "" : rpnSeqBuffer, iRmsResidual,
884 invalid ? NaN : rmsResidual, iChiSq, invalid ? NaN : chi, iTerms,
885 terms, iSigLevel,
886 invalid ? NaN : ChiSqrSigLevel(chi, points - terms), iOffset,
887 invalid ? NaN : xOffset, iFactor, invalid ? NaN : xScaleFactor,
888 iFitIsValid, isFit ? 'y' : 'n', -1) ||
889 !SDDS_WritePage(&SDDSout))
890 SDDS_PrintErrors(stderr,
891 SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
892 if (!invalid) {
893 free(diff);
894 free(sy);
895 if (xOrig != x)
896 free(xOrig);
897 if (yOrig != y)
898 free(yOrig);
899 if (syOrig != sy0)
900 free(syOrig);
901 if (sxOrig != sx)
902 free(sxOrig);
903 free(x);
904 free(y);
905 free(sx);
906 free(sy0);
907 }
908 }
909 if (!SDDS_Terminate(&SDDSin) || !SDDS_Terminate(&SDDSout)) {
910 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
911 exit(EXIT_FAILURE);
912 }
913 if (evalParameters.initialized && !SDDS_Terminate(&(evalParameters.dataset)))
914 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
915 /* free_scanargs(&s_arg, argc); */
916 free(coef);
917 free(coefSigma);
918 if (coefUnits)
919 free(coefUnits);
920 if (order)
921 free(order);
922
923 return (EXIT_SUCCESS);
924}
925
926void print_coefs(FILE *fpo, double xOffset, double xScaleFactor, long chebyshev,
927 double *coef, double *coefSigma, int32_t *order, long terms,
928 double chi, long normTerm, char *prepend) {
929 long i;
930
931 if (chebyshev)
932 fprintf(fpo, "%s%ld-term Chebyshev T polynomial least-squares fit about "
933 "x=%21.15e, scaled by %21.15e:\n",
934 prepend, terms, xOffset, xScaleFactor);
935 else
936 fprintf(fpo, "%s%ld-term polynomial least-squares fit about x=%21.15e:\n",
937 prepend, terms, xOffset);
938 if (normTerm >= 0 && terms > normTerm) {
939 if (coef[normTerm] != 0)
940 fprintf(fpo, "%s coefficients are normalized with factor %21.15e to "
941 "make a[%ld]==1\n",
942 prepend, coef[normTerm], (order ? order[normTerm] : normTerm));
943 else {
944 fprintf(fpo, "%s can't normalize coefficients as requested: a[%ld]==0\n",
945 prepend, (order ? order[normTerm] : normTerm));
946 normTerm = -1;
947 }
948 } else
949 normTerm = -1;
950
951 for (i = 0; i < terms; i++) {
952 fprintf(fpo, "%sa[%ld] = %21.15e ", prepend, (order ? order[i] : i),
953 (normTerm < 0 ? coef[i] : coef[i] / coef[normTerm]));
954 if (coefSigma)
955 fprintf(
956 fpo, "+/- %21.15e\n",
957 (normTerm < 0 ? coefSigma[i] : coefSigma[i] / fabs(coef[normTerm])));
958 else
959 fputc('\n', fpo);
960 }
961 if (coefSigma)
962 fprintf(fpo, "%sreduced chi-squared = %21.15e\n", prepend, chi);
963}
964
965void makeFitLabel(char *buffer, long bufsize, char *fitLabelFormat,
966 double *coef, int32_t *order, long terms, long chebyshev) {
967 long i;
968 static char buffer1[SDDS_MAXLINE], buffer2[SDDS_MAXLINE];
969
970 sprintf(buffer, "%s = ", ySymbol);
971 for (i = 0; i < terms; i++) {
972 buffer1[0] = 0;
973 if (coef[i] == 0)
974 continue;
975 if (order[i] == 0) {
976 if (coef[i] > 0) {
977 strcat(buffer1, "+");
978 sprintf(buffer1 + 1, fitLabelFormat, coef[i]);
979 } else
980 sprintf(buffer1, fitLabelFormat, coef[i]);
981 } else {
982 if (coef[i] > 0) {
983 strcat(buffer1, "+");
984 sprintf(buffer1 + 1, fitLabelFormat, coef[i]);
985 } else
986 sprintf(buffer1, fitLabelFormat, coef[i]);
987 if (order[i] >= 1) {
988 strcat(buffer1, "*");
989 if (chebyshev != 1) {
990 strcat(buffer1, xSymbol);
991 if (order[i] > 1) {
992 sprintf(buffer2, "$a%d$n", order[i]);
993 strcat(buffer1, buffer2);
994 }
995 } else {
996 sprintf(buffer2, "T$b%d$n(%s)", order[i], xSymbol);
997 strcat(buffer1, buffer2);
998 }
999 }
1000 }
1001 if ((long)(strlen(buffer) + strlen(buffer1)) > (long)(0.95 * bufsize)) {
1002 fprintf(stderr, "buffer overflow making fit label!\n");
1003 return;
1004 }
1005 strcat(buffer, buffer1);
1006 }
1007}
1008
1009double rms_average(double *x, int64_t n) {
1010 double sum2;
1011 int64_t i;
1012
1013 for (i = sum2 = 0; i < n; i++)
1014 sum2 += sqr(x[i]);
1015
1016 return (sqrt(sum2 / n));
1017}
1018
1019long coefficient_index(int32_t *order, long terms, long order_of_interest) {
1020 long i;
1021 for (i = 0; i < terms; i++)
1022 if (order[i] == order_of_interest)
1023 return (i);
1024 return (-1);
1025}
1026
1027void checkInputFile(SDDS_DATASET *SDDSin, char *xName, char *yName,
1028 char *xSigmaName, char *ySigmaName) {
1029 char *ptr = NULL;
1030 if (!(ptr = SDDS_FindColumn(SDDSin, FIND_NUMERIC_TYPE, xName, NULL)))
1031 SDDS_Bomb("x column doesn't exist or is nonnumeric");
1032 free(ptr);
1033 ptr = NULL;
1034 if (!(ptr = SDDS_FindColumn(SDDSin, FIND_NUMERIC_TYPE, yName, NULL)))
1035 SDDS_Bomb("y column doesn't exist or is nonnumeric");
1036 free(ptr);
1037 ptr = NULL;
1038 if (xSigmaName &&
1039 !(ptr = SDDS_FindColumn(SDDSin, FIND_NUMERIC_TYPE, xSigmaName, NULL)))
1040 SDDS_Bomb("x sigma column doesn't exist or is nonnumeric");
1041 if (ptr)
1042 free(ptr);
1043 ptr = NULL;
1044 if (ySigmaName &&
1045 !(ptr = SDDS_FindColumn(SDDSin, FIND_NUMERIC_TYPE, ySigmaName, NULL)))
1046 SDDS_Bomb("y sigma column doesn't exist or is nonnumeric");
1047 if (ptr)
1048 free(ptr);
1049 ptr = NULL;
1050}
1051
1052char **initializeOutputFile(SDDS_DATASET *SDDSout, char *output,
1053 SDDS_DATASET *SDDSin, char *input, char *xName,
1054 char *yName, char *xSigmaName, char *ySigmaName,
1055 long sigmasValid, int32_t *order, long terms,
1056 long chebyshev, long copyParameters) {
1057 char buffer[SDDS_MAXLINE], buffer1[SDDS_MAXLINE], *xUnits, *yUnits,
1058 **coefUnits;
1059 long i;
1060
1061 if (!SDDS_InitializeOutput(SDDSout, SDDS_BINARY, 0, NULL, "sddspfit output",
1062 output) ||
1063 !SDDS_TransferColumnDefinition(SDDSout, SDDSin, xName, NULL) ||
1064 !SDDS_TransferColumnDefinition(SDDSout, SDDSin, yName, NULL) ||
1065 SDDS_GetColumnInformation(SDDSout, "symbol", &xSymbol, SDDS_GET_BY_NAME,
1066 xName) != SDDS_STRING ||
1067 SDDS_GetColumnInformation(SDDSout, "symbol", &ySymbol, SDDS_GET_BY_NAME,
1068 yName) != SDDS_STRING ||
1069 (xSigmaName &&
1070 !SDDS_TransferColumnDefinition(SDDSout, SDDSin, xSigmaName, NULL)) ||
1071 (ySigmaName &&
1072 !SDDS_TransferColumnDefinition(SDDSout, SDDSin, ySigmaName, NULL)))
1073 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1074 if (!xSymbol || SDDS_StringIsBlank(xSymbol))
1075 xSymbol = xName;
1076 if (!ySymbol || SDDS_StringIsBlank(ySymbol))
1077 ySymbol = yName;
1078 ix = SDDS_GetColumnIndex(SDDSout, xName);
1079 iy = SDDS_GetColumnIndex(SDDSout, yName);
1080 if (ySigmaName)
1081 iySigma = SDDS_GetColumnIndex(SDDSout, ySigmaName);
1082 if (xSigmaName)
1083 ixSigma = SDDS_GetColumnIndex(SDDSout, xSigmaName);
1084 if (SDDS_NumberOfErrors())
1085 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1086
1087 sprintf(buffer, "%sFit", yName);
1088 sprintf(buffer1, "Fit[%s]", ySymbol);
1089 if (!SDDS_TransferColumnDefinition(SDDSout, SDDSin, yName, buffer) ||
1090 !SDDS_ChangeColumnInformation(SDDSout, "symbol", buffer1,
1091 SDDS_SET_BY_NAME, buffer))
1092 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1093 if ((iFit = SDDS_GetColumnIndex(SDDSout, buffer)) < 0)
1094 SDDS_Bomb("unable to get index of just-defined fit output column");
1095
1096 sprintf(buffer, "%sResidual", yName);
1097 sprintf(buffer1, "Residual[%s]", ySymbol);
1098 if (!SDDS_TransferColumnDefinition(SDDSout, SDDSin, yName, buffer) ||
1099 !SDDS_ChangeColumnInformation(SDDSout, "symbol", buffer1,
1100 SDDS_SET_BY_NAME, buffer))
1101 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1102
1103 if ((iResidual = SDDS_GetColumnIndex(SDDSout, buffer)) < 0)
1104 SDDS_Bomb("unable to get index of just-defined residual output column");
1105
1106 if (sigmasValid && !ySigmaName) {
1107 sprintf(buffer, "%sSigma", yName);
1108 if (!SDDS_TransferColumnDefinition(SDDSout, SDDSin, yName, buffer))
1109 SDDS_PrintErrors(stderr,
1110 SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1111 iySigma = SDDS_GetColumnIndex(SDDSout, buffer);
1112 if (ySymbol && !SDDS_StringIsBlank(ySymbol)) {
1113 sprintf(buffer1, "Sigma[%s]", ySymbol);
1114 if (!SDDS_ChangeColumnInformation(SDDSout, "symbol", buffer1,
1115 SDDS_SET_BY_NAME, buffer))
1116 SDDS_PrintErrors(stderr,
1117 SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1118 }
1119 }
1120
1121 if (!(coefUnits = makeCoefficientUnits(SDDSout, xName, yName, order, terms)))
1122 SDDS_Bomb("unable to make coefficient units");
1123
1124 if (SDDS_DefineArray(SDDSout, "Order", NULL, NULL, "Order of term in fit",
1125 NULL, SDDS_LONG, 0, 1, "FitResults") < 0 ||
1126 SDDS_DefineArray(SDDSout, "Coefficient", "a", "[CoefficientUnits]",
1127 "Coefficient of term in fit", NULL, SDDS_DOUBLE, 0, 1,
1128 "FitResults") < 0 ||
1129 (sigmasValid &&
1130 SDDS_DefineArray(SDDSout, "CoefficientSigma", "$gs$r$ba$n",
1131 "[CoefficientUnits]",
1132 "Sigma of coefficient of term in fit", NULL,
1133 SDDS_DOUBLE, 0, 1, "FitResults") < 0) ||
1134 SDDS_DefineArray(SDDSout, "CoefficientUnits", NULL, NULL, NULL, NULL,
1135 SDDS_STRING, 0, 1, "FitResults") < 0)
1136 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1137
1138 if (SDDS_DefineParameter(SDDSout, "Basis", NULL, NULL,
1139 "Function basis for fit", NULL, SDDS_STRING,
1140 chebyshev ? (chebyshev == 1 ? "Chebyshev T polynomials" : "Converted Chebyshev T polynomials")
1141 : "ordinary polynomials") < 0 ||
1142 (iChiSq = SDDS_DefineParameter(
1143 SDDSout, "ReducedChiSquared", "$gh$r$a2$n/(N-M)", NULL,
1144 "Reduced chi-squared of fit", NULL, SDDS_DOUBLE, NULL)) < 0 ||
1145 SDDS_GetColumnInformation(SDDSout, "units", &yUnits, SDDS_GET_BY_NAME,
1146 yName) != SDDS_STRING ||
1147 (iRmsResidual = SDDS_DefineParameter(
1148 SDDSout, "RmsResidual", "$gs$r$bres$n", yUnits,
1149 "RMS residual of fit", NULL, SDDS_DOUBLE, NULL)) < 0 ||
1150 (iSigLevel =
1151 SDDS_DefineParameter(SDDSout, "SignificanceLevel", NULL, NULL,
1152 "Probability that data is from fit function",
1153 NULL, SDDS_DOUBLE, NULL)) < 0 ||
1154 (iRpnSequence = SDDS_DefineParameter(SDDSout, "RpnSequence", NULL, NULL,
1155 "Rpn sequence to evaluate the fit",
1156 NULL, SDDS_STRING, NULL)) < 0) {
1157 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
1158 exit(EXIT_FAILURE);
1159 }
1160 if (yUnits)
1161 free(yUnits);
1162
1163 if (SDDS_GetColumnInformation(SDDSout, "units", &xUnits, SDDS_GET_BY_NAME,
1164 xName) != SDDS_STRING)
1165 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1166 sprintf(buffer, "%sOffset", xName);
1167 sprintf(buffer1, "Offset of %s for fit", xName);
1168 if ((iOffset = SDDS_DefineParameter(SDDSout, buffer, NULL, xUnits, buffer1,
1169 NULL, SDDS_DOUBLE, NULL)) < 0)
1170 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1171 sprintf(buffer, "%sScale", xName);
1172 sprintf(buffer1, "Scale factor of %s for fit", xName);
1173 if ((iFactor = SDDS_DefineParameter(SDDSout, buffer, NULL, xUnits, buffer1,
1174 NULL, SDDS_DOUBLE, NULL)) < 0)
1175 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1176
1177 if ((iFitIsValid = SDDS_DefineParameter(SDDSout, "FitIsValid", NULL, NULL,
1178 NULL, NULL, SDDS_CHARACTER, NULL)) <
1179 0)
1180 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1181
1182 if ((iTerms = SDDS_DefineParameter(SDDSout, "Terms", NULL, NULL,
1183 "Number of terms in fit", NULL, SDDS_LONG,
1184 NULL)) < 0)
1185 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1186
1187 iFitLabel = SDDS_DefineParameter(SDDSout, "sddspfitLabel", NULL, NULL, NULL,
1188 NULL, SDDS_STRING, NULL);
1189 if (!chebyshev) {
1190 if ((i = coefficient_index(order, terms, 0)) >= 0) {
1191 iIntercept =
1192 SDDS_DefineParameter(SDDSout, "Intercept", "Intercept", coefUnits[i],
1193 "Intercept of fit", NULL, SDDS_DOUBLE, NULL);
1194 if (sigmasValid)
1195 iInterceptSigma = SDDS_DefineParameter(
1196 SDDSout, "InterceptSigma", "InterceptSigma", coefUnits[i],
1197 "Sigma of intercept of fit", NULL, SDDS_DOUBLE, NULL);
1198 }
1199 if ((i = coefficient_index(order, terms, 1)) >= 0) {
1200 iSlope = SDDS_DefineParameter(SDDSout, "Slope", "Slope", coefUnits[i],
1201 "Slope of fit", NULL, SDDS_DOUBLE, NULL);
1202 if (sigmasValid)
1203 iSlopeSigma = SDDS_DefineParameter(SDDSout, "SlopeSigma", "SlopeSigma",
1204 coefUnits[i], "Sigma of slope of fit",
1205 NULL, SDDS_DOUBLE, NULL);
1206 }
1207 if ((i = coefficient_index(order, terms, 2)) >= 0) {
1208 iCurvature =
1209 SDDS_DefineParameter(SDDSout, "Curvature", "Curvature", coefUnits[i],
1210 "Curvature of fit", NULL, SDDS_DOUBLE, NULL);
1211 if (sigmasValid)
1212 iCurvatureSigma = SDDS_DefineParameter(
1213 SDDSout, "CurvatureSigma", "CurvatureSigma", coefUnits[i],
1214 "Sigma of curvature of fit", NULL, SDDS_DOUBLE, NULL);
1215 }
1216 }
1217
1218 for (i = 0; i < terms; i++) {
1219 char s[100];
1220 sprintf(s, "Coefficient%02ld", (long)order[i]);
1221 iTerm[i] = SDDS_DefineParameter(SDDSout, s, s, coefUnits[i], NULL, NULL,
1222 SDDS_DOUBLE, NULL);
1223 }
1224 for (i = 0; i < terms; i++) {
1225 char s[100];
1226 if (sigmasValid) {
1227 sprintf(s, "Coefficient%02ldSigma", (long)order[i]);
1228 iTermSig[i] = SDDS_DefineParameter(SDDSout, s, s, coefUnits[i], NULL,
1229 NULL, SDDS_DOUBLE, NULL);
1230 } else {
1231 iTermSig[i] = -1;
1232 }
1233 }
1234
1235 if (SDDS_NumberOfErrors())
1236 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1237
1238 if (copyParameters &&
1239 !SDDS_TransferAllParameterDefinitions(SDDSout, SDDSin, 0))
1240 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1241
1242 if (!SDDS_WriteLayout(SDDSout))
1243 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1244
1245 return coefUnits;
1246}
1247
1248long setCoefficientData(SDDS_DATASET *SDDSout, double *coef, double *coefSigma,
1249 char **coefUnits, int32_t *order, long terms, long chebyshev,
1250 char *fitLabelFormat, char *rpnSeqBuffer) {
1251 long termIndex, i;
1252 long invalid = 0;
1253 static char fitLabelBuffer[SDDS_MAXLINE];
1254
1255 if (chebyshev != 2) {
1256 createRpnSequence(rpnSeqBuffer, SDDS_MAXLINE, coef, order, terms);
1257 if (!SDDS_SetArrayVararg(SDDSout, "Order", SDDS_POINTER_ARRAY, order,
1258 terms) ||
1259 !SDDS_SetArrayVararg(SDDSout, "Coefficient", SDDS_POINTER_ARRAY, coef,
1260 terms) ||
1261 (coefSigma &&
1262 !SDDS_SetArrayVararg(SDDSout, "CoefficientSigma", SDDS_POINTER_ARRAY,
1263 coefSigma, terms)) ||
1264 !SDDS_SetArrayVararg(SDDSout, "CoefficientUnits", SDDS_POINTER_ARRAY,
1265 coefUnits, terms))
1266 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1267
1268 termIndex = coefficient_index(order, terms, 0);
1269
1270 if (iIntercept != -1 &&
1271 !SDDS_SetParameters(SDDSout, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE,
1272 iIntercept, invalid ? NaN : coef[termIndex], -1))
1273 SDDS_PrintErrors(stderr,
1274 SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1275 if (iInterceptSigma != -1 &&
1276 !SDDS_SetParameters(SDDSout, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE,
1277 iInterceptSigma,
1278 invalid ? NaN : coefSigma[termIndex], -1))
1279 SDDS_PrintErrors(stderr,
1280 SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1281 if (!invalid)
1282 termIndex = coefficient_index(order, terms, 1);
1283 if (iSlope != -1 &&
1284 !SDDS_SetParameters(SDDSout, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE,
1285 iSlope, invalid ? NaN : coef[termIndex], -1))
1286 SDDS_PrintErrors(stderr,
1287 SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1288 if (iSlopeSigma != -1 &&
1289 !SDDS_SetParameters(SDDSout, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE,
1290 iSlopeSigma, invalid ? NaN : coefSigma[termIndex],
1291 -1))
1292 SDDS_PrintErrors(stderr,
1293 SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1294 if (!invalid)
1295 termIndex = coefficient_index(order, terms, 2);
1296 if (iCurvature != -1 &&
1297 !SDDS_SetParameters(SDDSout, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE,
1298 iCurvature, invalid ? NaN : coef[termIndex], -1))
1299 SDDS_PrintErrors(stderr,
1300 SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1301 if (iCurvatureSigma != -1 &&
1302 !SDDS_SetParameters(SDDSout, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE,
1303 iCurvatureSigma,
1304 invalid ? NaN : coefSigma[termIndex], -1))
1305 SDDS_PrintErrors(stderr,
1306 SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1307 if (iFitLabel != -1 && !invalid) {
1308 makeFitLabel(fitLabelBuffer, SDDS_MAXLINE, fitLabelFormat, coef, order,
1309 terms, chebyshev);
1310 if (!SDDS_SetParameters(SDDSout, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE,
1311 iFitLabel, fitLabelBuffer, -1))
1312 SDDS_PrintErrors(stderr,
1313 SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1314 }
1315 for (i = 0; i < terms; i++) {
1316 if (!SDDS_SetParameters(SDDSout, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE,
1317 iTerm[i], coef[i], -1))
1318 SDDS_PrintErrors(stderr,
1319 SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1320 if (iTermSig[i] != -1)
1321 if (!SDDS_SetParameters(SDDSout,
1322 SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE,
1323 iTermSig[i], coefSigma[i], -1))
1324 SDDS_PrintErrors(stderr,
1325 SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1326 }
1327 } else {
1328 long termsC;
1329 int32_t *orderC;
1330 double *coefC, *coefSigmaC;
1331 convertFromChebyshev(terms, order, coef, coefSigma, &termsC, &orderC, &coefC, &coefSigmaC);
1332 setCoefficientData(SDDSout, coefC, coefSigmaC, coefUnits, orderC, termsC, 0, fitLabelFormat,
1333 rpnSeqBuffer);
1334 }
1335
1336 return 1;
1337}
1338
1339char **makeCoefficientUnits(SDDS_DATASET *SDDSout, char *xName, char *yName,
1340 int32_t *order, long terms) {
1341 char *xUnits, *yUnits, buffer[SDDS_MAXLINE];
1342 char **coefUnits = NULL;
1343 long i;
1344
1345 if (!SDDS_GetColumnInformation(SDDSout, "units", &xUnits, SDDS_GET_BY_NAME,
1346 xName) ||
1347 !SDDS_GetColumnInformation(SDDSout, "units", &yUnits, SDDS_GET_BY_NAME,
1348 yName))
1349 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1350
1351 coefUnits = tmalloc(sizeof(*coefUnits) * terms);
1352 if (!xUnits || SDDS_StringIsBlank(xUnits)) {
1353 if (!yUnits || SDDS_StringIsBlank(yUnits))
1354 SDDS_CopyString(&yUnits, "");
1355 for (i = 0; i < terms; i++)
1356 coefUnits[i] = yUnits;
1357 } else {
1358 if (!yUnits || SDDS_StringIsBlank(yUnits))
1359 SDDS_CopyString(&yUnits, "1");
1360 for (i = 0; i < terms; i++) {
1361 if (order[i] == 0) {
1362 if (strcmp(yUnits, "1") != 0)
1363 SDDS_CopyString(coefUnits + i, yUnits);
1364 else
1365 SDDS_CopyString(coefUnits + i, "");
1366 } else if (strcmp(xUnits, yUnits) == 0) {
1367 if (order[i] > 1)
1368 sprintf(buffer, "1/%s$a%d$n", xUnits, order[i] - 1);
1369 else
1370 strcpy(buffer, "");
1371 SDDS_CopyString(coefUnits + i, buffer);
1372 } else {
1373 if (order[i] > 1)
1374 sprintf(buffer, "%s/%s$a%d$n", yUnits, xUnits, order[i]);
1375 else
1376 sprintf(buffer, "%s/%s", yUnits, xUnits);
1377 SDDS_CopyString(coefUnits + i, buffer);
1378 }
1379 }
1380 }
1381 return coefUnits;
1382}
1383
1384void compareOriginalToFit(double *x, double *y, double **residual,
1385 int64_t points, double *rmsResidual, double *coef,
1386 int32_t *order, long terms) {
1387 int64_t i;
1388 double residualSum2, fit;
1389
1390 *residual = tmalloc(sizeof(**residual) * points);
1391
1392 residualSum2 = 0;
1393 for (i = 0; i < points; i++) {
1394 fit = eval_sum(basis_fn, coef, order, terms, x[i]);
1395 (*residual)[i] = y[i] - fit;
1396 residualSum2 += sqr((*residual)[i]);
1397 }
1398 *rmsResidual = sqrt(residualSum2 / points);
1399}
1400
1401void makeEvaluationTable(EVAL_PARAMETERS *evalParameters, double *x,
1402 int64_t points, double *coef, int32_t *order,
1403 long terms, SDDS_DATASET *SDDSin, char *xName,
1404 char *yName) {
1405 double *xEval, *yEval, delta;
1406 int64_t i;
1407 yEval = NULL;
1408 if (!evalParameters->initialized) {
1409 if (!SDDS_InitializeOutput(&evalParameters->dataset, SDDS_BINARY, 0, NULL,
1410 "sddspfit evaluation table",
1411 evalParameters->file) ||
1412 !SDDS_TransferColumnDefinition(&evalParameters->dataset, SDDSin, xName,
1413 NULL) ||
1414 !SDDS_TransferColumnDefinition(&evalParameters->dataset, SDDSin, yName,
1415 NULL) ||
1416 !SDDS_WriteLayout(&evalParameters->dataset))
1417 SDDS_PrintErrors(stderr,
1418 SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1419 evalParameters->initialized = 1;
1420 }
1421
1422 if (evalParameters->flags & EVAL_VALUESFILE_GIVEN) {
1423 if (!evalParameters->inputInitialized) {
1424 if (!SDDS_InitializeInput(&(evalParameters->valuesDataset), evalParameters->valuesFile)) {
1425 fprintf(stderr, "error: unable to initialize %s\n", evalParameters->valuesFile);
1426 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1427 }
1428 if (!SDDS_ReadPage(&(evalParameters->valuesDataset))) {
1429 fprintf(stderr, "error: unable to read page from %s\n", evalParameters->valuesFile);
1430 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1431 }
1432 evalParameters->inputInitialized = 1;
1433 } else {
1434 if (!(evalParameters->flags & EVAL_REUSE_PAGE_GIVEN) &&
1435 !SDDS_ReadPage(&(evalParameters->valuesDataset))) {
1436 fprintf(stderr, "error: unable to read page from %s\n", evalParameters->valuesFile);
1437 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1438 }
1439 }
1440 if (!(xEval = SDDS_GetColumnInDoubles(&(evalParameters->valuesDataset), evalParameters->valuesColumn))) {
1441 fprintf(stderr, "error: unable to read column %s from %s\n", evalParameters->valuesColumn, evalParameters->valuesFile);
1442 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1443 }
1444 evalParameters->number = SDDS_CountRowsOfInterest(&(evalParameters->valuesDataset));
1445 } else {
1446 if (!(evalParameters->flags & EVAL_BEGIN_GIVEN) ||
1447 !(evalParameters->flags & EVAL_END_GIVEN)) {
1448 double min, max;
1449 find_min_max(&min, &max, x, points);
1450 if (!(evalParameters->flags & EVAL_BEGIN_GIVEN))
1451 evalParameters->begin = min;
1452 if (!(evalParameters->flags & EVAL_END_GIVEN))
1453 evalParameters->end = max;
1454 }
1455 if (!(evalParameters->flags & EVAL_NUMBER_GIVEN))
1456 evalParameters->number = points;
1457 if (evalParameters->number > 1)
1458 delta = (evalParameters->end - evalParameters->begin) /
1459 (evalParameters->number - 1);
1460 else
1461 delta = 0;
1462
1463 if (!(xEval = (double *)malloc(sizeof(*xEval) * evalParameters->number)))
1464 SDDS_Bomb("allocation failure");
1465
1466 for (i = 0; i < evalParameters->number; i++)
1467 xEval[i] = evalParameters->begin + i * delta;
1468 }
1469
1470 if (!(yEval = (double *)malloc(sizeof(*yEval) * evalParameters->number)))
1471 SDDS_Bomb("allocation failure");
1472 for (i = 0; i < evalParameters->number; i++)
1473 yEval[i] = eval_sum(basis_fn, coef, order, terms, xEval[i]);
1474
1475 if (!SDDS_StartPage(&evalParameters->dataset, evalParameters->number) ||
1476 !SDDS_SetColumnFromDoubles(&evalParameters->dataset, SDDS_SET_BY_NAME,
1477 xEval, evalParameters->number, xName) ||
1478 !SDDS_SetColumnFromDoubles(&evalParameters->dataset, SDDS_SET_BY_NAME,
1479 yEval, evalParameters->number, yName) ||
1480 !SDDS_WritePage(&evalParameters->dataset))
1481 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1482 free(xEval);
1483 free(yEval);
1484}
1485
1486long reviseFitOrders(double *x, double *y, double *sy, int64_t points,
1487 long terms, int32_t *order, double *coef,
1488 double *coefSigma, double *diff,
1489 double (*basis_fn)(double xa, long ordera),
1490 unsigned long reviseOrders, double xOffset,
1491 double xScaleFactor, long normTerm, long ySigmasValid,
1492 long chebyshev, double revpowThreshold,
1493 double revpowCompleteThreshold, double goodEnoughChi) {
1494 double bestChi, chi;
1495 long bestTerms, newTerms, newBest, *termUsed;
1496 int32_t *newOrder, *bestOrder;
1497 long i, ip, j;
1498 long origTerms, *origOrder;
1499
1500 bestOrder = tmalloc(sizeof(*bestOrder) * terms);
1501 newOrder = tmalloc(sizeof(*newOrder) * terms);
1502 termUsed = tmalloc(sizeof(*termUsed) * terms);
1503 origOrder = tmalloc(sizeof(*origOrder) * terms);
1504 origTerms = terms;
1505 for (i = 0; i < terms; i++)
1506 origOrder[i] = order[i];
1507 qsort((void *)order, terms, sizeof(*order), long_cmpasc);
1508 bestOrder[0] = newOrder[0] = order[0];
1509 termUsed[0] = 1;
1510 newTerms = bestTerms = 1;
1511 /* do a fit */
1512 if (!lsfg(x, y, sy, points, newTerms, newOrder, coef, coefSigma, &bestChi,
1513 diff, basis_fn))
1514 SDDS_Bomb("revise-orders fit failed.");
1515 if (reviseOrders & REVPOW_VERBOSE) {
1516 fputs("fit to revise orders:", stdout);
1517 print_coefs(stdout, xOffset, xScaleFactor, chebyshev, coef,
1518 (ySigmasValid ? coefSigma : NULL), bestOrder, bestTerms,
1519 bestChi, normTerm, "");
1520 fprintf(stdout, "unweighted rms deviation from fit: %21.15e\n",
1521 rms_average(diff, points));
1522 }
1523
1524 do {
1525 newBest = 0;
1526 newTerms = newTerms + 1;
1527 for (ip = 1; ip < terms; ip++) {
1528 if (termUsed[ip])
1529 continue;
1530 newOrder[newTerms - 1] = order[ip];
1531 if (!lsfg(x, y, sy, points, newTerms, newOrder, coef, coefSigma, &chi,
1532 diff, basis_fn))
1533 SDDS_Bomb("revise-orders fit failed.");
1534 if (reviseOrders & REVPOW_VERBOSE) {
1535 fputs("trial fit:", stdout);
1536 print_coefs(stdout, xOffset, xScaleFactor, chebyshev, coef,
1537 (ySigmasValid ? coefSigma : NULL), newOrder, newTerms,
1538 chi, normTerm, "");
1539 fprintf(stdout, "unweighted rms deviation from fit: %21.15e\n",
1540 rms_average(diff, points));
1541 }
1542 if ((bestChi > goodEnoughChi && chi < bestChi) ||
1543 (chi + revpowThreshold < bestChi && newTerms < bestTerms)) {
1544 bestChi = chi;
1545 bestTerms = newTerms;
1546 newBest = 1;
1547 termUsed[ip] = 1;
1548 for (i = 0; i < newTerms; i++)
1549 bestOrder[i] = newOrder[i];
1550 if (reviseOrders & REVPOW_VERBOSE) {
1551 fputs("new best fit:", stdout);
1552 print_coefs(stdout, xOffset, xScaleFactor, chebyshev, coef,
1553 (ySigmasValid ? coefSigma : NULL), bestOrder, bestTerms,
1554 bestChi, normTerm, "");
1555 fprintf(stdout, "unweighted rms deviation from fit: %21.15e\n",
1556 rms_average(diff, points));
1557 }
1558 break;
1559 }
1560 }
1561 } while (newBest && bestChi > goodEnoughChi);
1562
1563 terms = bestTerms;
1564 for (ip = 0; ip < terms; ip++)
1565 order[ip] = bestOrder[ip];
1566
1567 if (newBest) {
1568 do {
1569 newBest = 0;
1570 for (ip = 0; ip < terms; ip++) {
1571 for (i = j = 0; i < terms; i++) {
1572 if (i != ip)
1573 newOrder[j++] = order[i];
1574 }
1575 newTerms = terms - 1;
1576 if (!lsfg(x, y, sy, points, newTerms, newOrder, coef, coefSigma, &chi, diff,
1577 basis_fn))
1578 SDDS_Bomb("revise-orders fit failed.");
1579 if ((bestChi > goodEnoughChi && chi < goodEnoughChi) ||
1580 (chi + revpowThreshold < bestChi && newTerms < terms)) {
1581 bestChi = chi;
1582 terms = newTerms;
1583 newBest = 1;
1584 for (i = 0; i < newTerms; i++)
1585 order[i] = newOrder[i];
1586 if (reviseOrders & REVPOW_VERBOSE) {
1587 fputs("new best fit:", stdout);
1588 print_coefs(stdout, xOffset, xScaleFactor, chebyshev, coef,
1589 (ySigmasValid ? coefSigma : NULL), order, terms,
1590 bestChi, normTerm, "");
1591 fprintf(stdout, "unweighted rms deviation from fit: %21.15e\n",
1592 rms_average(diff, points));
1593 }
1594 break;
1595 }
1596 }
1597 } while (newBest && terms > 1 && bestChi > goodEnoughChi);
1598 }
1599
1600 free(bestOrder);
1601 free(termUsed);
1602 free(newOrder);
1603
1604 if ((reviseOrders & REVPOW_COMPLETE) && bestChi > revpowCompleteThreshold) {
1605 terms = origTerms;
1606 for (i = 0; i < origTerms; i++)
1607 order[i] = origOrder[i];
1608 if (reviseOrders & REVPOW_VERBOSE)
1609 fprintf(stdout, "Result unsatisfactory---attempting complete trials\n");
1610 return reviseFitOrders1(x, y, sy, points, terms, order, coef, coefSigma,
1611 diff, basis_fn, reviseOrders, xOffset, xScaleFactor,
1612 normTerm, ySigmasValid, chebyshev, revpowThreshold,
1613 goodEnoughChi);
1614 }
1615
1616 free(origOrder);
1617 return terms;
1618}
1619
1620long reviseFitOrders1(double *x, double *y, double *sy, int64_t points,
1621 long terms, int32_t *order, double *coef,
1622 double *coefSigma, double *diff,
1623 double (*basis_fn)(double xa, long ordera),
1624 unsigned long reviseOrders, double xOffset,
1625 double xScaleFactor, long normTerm, long ySigmasValid,
1626 long chebyshev, double revpowThreshold,
1627 double goodEnoughChi) {
1628 double bestChi, chi;
1629 long bestTerms, newTerms;
1630 int32_t *newOrder = NULL, *bestOrder;
1631 long i, ip, j;
1632 long *counter = NULL, *counterLim = NULL;
1633
1634 if (!(bestOrder = malloc(sizeof(*bestOrder) * terms)) ||
1635 !(newOrder = malloc(sizeof(*newOrder) * terms)) ||
1636 !(counter = calloc(sizeof(*counter), terms)) ||
1637 !(counterLim = calloc(sizeof(*counterLim), terms))) {
1638 fprintf(stderr, "Error: memory allocation failure (%ld terms)\n", terms);
1639 SDDS_Bomb(NULL);
1640 }
1641 for (i = 0; i < terms; i++)
1642 counterLim[i] = 2;
1643 qsort((void *)order, terms, sizeof(*order), long_cmpasc);
1644 /* do a fit */
1645 if (!lsfg(x, y, sy, points, 2, order, coef, coefSigma, &bestChi, diff,
1646 basis_fn))
1647 SDDS_Bomb("revise-orders fit failed.");
1648 for (i = 0; i < 2; i++)
1649 bestOrder[i] = order[i];
1650 bestTerms = 2;
1651 if (reviseOrders & REVPOW_VERBOSE) {
1652 fputs("starting fit to revise orders:", stdout);
1653 print_coefs(stdout, xOffset, xScaleFactor, chebyshev, coef,
1654 (ySigmasValid ? coefSigma : NULL), order, 1, bestChi, normTerm,
1655 "");
1656 fprintf(stdout, "unweighted rms deviation from fit: %21.15e\n",
1657 rms_average(diff, points));
1658 }
1659
1660 newTerms = 1;
1661 while (advance_counter(counter, counterLim, terms) >= 0) {
1662 for (i = j = 0; i < terms; i++) {
1663 if (counter[i])
1664 newOrder[j++] = order[i];
1665 }
1666 newTerms = j;
1667 if (!lsfg(x, y, sy, points, newTerms, newOrder, coef, coefSigma, &chi, diff,
1668 basis_fn))
1669 SDDS_Bomb("revise-orders fit failed.");
1670 if ((chi < goodEnoughChi && newTerms < bestTerms) ||
1671 (bestChi > goodEnoughChi && chi < bestChi)) {
1672 bestChi = chi;
1673 bestTerms = newTerms;
1674 for (i = 0; i < newTerms; i++)
1675 bestOrder[i] = newOrder[i];
1676 if (reviseOrders & REVPOW_VERBOSE) {
1677 fputs("new best fit:", stdout);
1678 print_coefs(stdout, xOffset, xScaleFactor, chebyshev, coef,
1679 (ySigmasValid ? coefSigma : NULL), bestOrder, bestTerms,
1680 bestChi, normTerm, "");
1681 fprintf(stdout, "unweighted rms deviation from fit: %21.15e\n",
1682 rms_average(diff, points));
1683 }
1684 }
1685 }
1686
1687 terms = bestTerms;
1688 for (ip = 0; ip < terms; ip++)
1689 order[ip] = bestOrder[ip];
1690
1691 free(bestOrder);
1692 free(newOrder);
1693 free(counter);
1694 free(counterLim);
1695 return terms;
1696}
1697
1698void createRpnSequence(char *buffer, long bufsize, double *coef, int32_t *order,
1699 long terms) {
1700 long i, j, maxOrder;
1701 static char buffer1[SDDS_MAXLINE];
1702 double coef1;
1703 double offset;
1704
1705 buffer[0] = 0;
1706 maxOrder = 0;
1707 for (i = 0; i < terms; i++) {
1708 if (maxOrder < order[i])
1709 maxOrder = order[i];
1710 }
1711 if (maxOrder == 0) {
1712 snprintf(buffer, SDDS_MAXLINE, "%.15e", coef[0]);
1713 return;
1714 }
1715 offset = get_argument_offset();
1716 if (offset != 0)
1717 snprintf(buffer, SDDS_MAXLINE, "%le - ", offset);
1718 for (i = 2; i <= maxOrder; i++) {
1719 strcat(buffer, "= ");
1720 if ((strlen(buffer) + 2) > bufsize) {
1721 fprintf(stderr, "buffer overflow making rpn expression!\n");
1722 return;
1723 }
1724 }
1725 for (i = maxOrder; i >= 0; i--) {
1726 for (j = 0; j < terms; j++) {
1727 if (order[j] == i)
1728 break;
1729 }
1730 if (j == terms)
1731 coef1 = 0;
1732 else
1733 coef1 = coef[j];
1734 if (i == maxOrder)
1735 snprintf(buffer1, SDDS_MAXLINE, "%.15g * ", coef1);
1736 else if (i == 0 && order[j] == 0) {
1737 if (coef1 != 0)
1738 snprintf(buffer1, SDDS_MAXLINE, "%.15g + ", coef1);
1739 else
1740 continue;
1741 } else {
1742 if (coef1 == 0)
1743 strcpy(buffer1, "* ");
1744 else
1745 snprintf(buffer1, SDDS_MAXLINE, "%.15g + * ", coef1);
1746 }
1747 if ((strlen(buffer) + strlen(buffer1)) >= bufsize) {
1748 fprintf(stderr, "buffer overflow making rpn expression!\n");
1749 return;
1750 }
1751 strcat(buffer, buffer1);
1752 }
1753}
1754
1755/* Make array of structures giving coefficients for Chebyshev T polynomials,
1756 * which allows converting to ordinary polynomials.
1757 */
1758
1759CHEBYSHEV_COEF *makeChebyshevCoefficients(long maxOrder, long *nPoly) {
1760 CHEBYSHEV_COEF *coef;
1761 long i, j;
1762
1763 if (maxOrder < 2)
1764 *nPoly = 2;
1765 else
1766 *nPoly = maxOrder + 1;
1767
1768 coef = tmalloc(sizeof(*coef) * (*nPoly));
1769
1770 coef[0].nTerms = 1;
1771 coef[0].coef = tmalloc(sizeof(*(coef[0].coef)) * coef[0].nTerms);
1772 coef[0].coef[0] = 1;
1773
1774 coef[1].nTerms = 2;
1775 coef[1].coef = tmalloc(sizeof(*(coef[1].coef)) * coef[1].nTerms);
1776 coef[1].coef[0] = 0;
1777 coef[1].coef[1] = 1;
1778
1779 for (i = 2; i < *nPoly; i++) {
1780 coef[i].nTerms = coef[i - 1].nTerms + 1;
1781 coef[i].coef = calloc(coef[i].nTerms, sizeof(*(coef[i].coef)));
1782 for (j = 0; j < coef[i - 2].nTerms; j++)
1783 coef[i].coef[j] = -coef[i - 2].coef[j];
1784 for (j = 0; j < coef[i - 1].nTerms; j++)
1785 coef[i].coef[j + 1] += 2 * coef[i - 1].coef[j];
1786 }
1787 /*
1788 for (i = 0; i < *nPoly; i++) {
1789 printf("T%ld: ", i);
1790 for (j = 0; j < coef[i].nTerms; j++) {
1791 if (coef[i].coef[j])
1792 printf("%c%lg*x^%ld ", coef[i].coef[j] < 0 ? '-' : '+', fabs(coef[i].coef[j]), j);
1793 }
1794 printf("\n");
1795 }
1796 */
1797 return coef;
1798}
1799
1800void convertFromChebyshev(long termsT, int32_t *orderT, double *coefT, double *coefSigmaT,
1801 long *termsOrdinaryRet, int32_t **orderOrdinaryRet, double **coefOrdinaryRet, double **coefSigmaOrdinaryRet) {
1802 long i, maxOrder;
1803 long termsOrdinary;
1804 int32_t *orderOrdinary;
1805 double *coefOrdinary, *coefSigmaOrdinary, scale;
1806 static CHEBYSHEV_COEF *chebyCoef = NULL;
1807 static long nChebyCoef = 0, chebyMaxOrder = 0;
1808
1809 maxOrder = 0;
1810 for (i = 0; i < termsT; i++)
1811 if (orderT[i] > maxOrder)
1812 maxOrder = orderT[i];
1813
1814 termsOrdinary = maxOrder + 1;
1815 orderOrdinary = tmalloc(sizeof(*orderOrdinary) * termsOrdinary);
1816 coefOrdinary = calloc(termsOrdinary, sizeof(*coefOrdinary));
1817 if (coefSigmaT)
1818 coefSigmaOrdinary = calloc(termsOrdinary, sizeof(*coefSigmaOrdinary));
1819 else
1820 coefSigmaOrdinary = NULL;
1821
1822 if (chebyCoef == NULL || maxOrder > chebyMaxOrder) {
1823 if (chebyCoef) {
1824 for (i = 0; i < nChebyCoef; i++)
1825 free(chebyCoef[i].coef);
1826 free(chebyCoef);
1827 }
1828 chebyCoef = makeChebyshevCoefficients(chebyMaxOrder = maxOrder, &nChebyCoef);
1829 }
1830
1831 for (i = 0; i < termsT; i++) {
1832 long j;
1833 for (j = 0; j < chebyCoef[orderT[i]].nTerms; j++) {
1834 coefOrdinary[j] += coefT[i] * chebyCoef[i].coef[j];
1835 if (coefSigmaT)
1836 coefSigmaOrdinary[j] += sqr(coefSigmaT[i] * chebyCoef[i].coef[j]);
1837 }
1838 }
1839 scale = get_argument_scale();
1840 for (i = 0; i < termsOrdinary; i++) {
1841 if (coefSigmaT)
1842 coefSigmaOrdinary[i] = sqrt(coefSigmaOrdinary[i]) / ipow(scale, i);
1843 orderOrdinary[i] = i;
1844 coefOrdinary[i] /= ipow(scale, i);
1845 }
1846 *termsOrdinaryRet = termsOrdinary;
1847 *orderOrdinaryRet = orderOrdinary;
1848 *coefOrdinaryRet = coefOrdinary;
1849 *coefSigmaOrdinaryRet = coefSigmaOrdinary;
1850}
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_SetArrayVararg
int32_t SDDS_SetArrayVararg(SDDS_DATASET *SDDS_dataset, char *array_name, int32_t mode, void *data_pointer,...)
Sets the values of an array variable in the SDDS dataset using variable arguments for dimensions.
Definition SDDS_dataprep.c:1157
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_CountRowsOfInterest
int64_t SDDS_CountRowsOfInterest(SDDS_DATASET *SDDS_dataset)
Counts the number of rows marked as "of interest" in the current data table.
Definition SDDS_extract.c:364
SDDS_GetColumnInDoubles
double * SDDS_GetColumnInDoubles(SDDS_DATASET *SDDS_dataset, char *column_name)
Retrieves the data of a specified numerical column as an array of doubles, considering only rows mark...
Definition SDDS_extract.c:872
SDDS_ChangeColumnInformation
int32_t SDDS_ChangeColumnInformation(SDDS_DATASET *SDDS_dataset, char *field_name, void *memory, int32_t mode,...)
Modifies a specific field in a column definition within the SDDS dataset.
Definition SDDS_info.c:364
SDDS_GetColumnInformation
int32_t SDDS_GetColumnInformation(SDDS_DATASET *SDDS_dataset, char *field_name, void *memory, int32_t mode,...)
Retrieves information about a specified column in the SDDS dataset.
Definition SDDS_info.c:41
SDDS_InitializeInput
int32_t SDDS_InitializeInput(SDDS_DATASET *SDDS_dataset, char *filename)
Definition SDDS_input.c:49
SDDS_Terminate
int32_t SDDS_Terminate(SDDS_DATASET *SDDS_dataset)
Definition SDDS_input.c:1448
SDDS_ReadPage
int32_t SDDS_ReadPage(SDDS_DATASET *SDDS_dataset)
Definition SDDS_input.c:1006
SDDS_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:595
SDDS_DefineArray
int32_t SDDS_DefineArray(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, int32_t dimensions, const char *group_name)
Defines a data array within the SDDS dataset.
Definition SDDS_output.c:1592
SDDS_WritePage
int32_t SDDS_WritePage(SDDS_DATASET *SDDS_dataset)
Writes the current data table to the output file.
Definition SDDS_output.c:1222
SDDS_WriteLayout
int32_t SDDS_WriteLayout(SDDS_DATASET *SDDS_dataset)
Writes the SDDS layout header to the output file.
Definition SDDS_output.c:893
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:1466
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_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_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_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
advance_counter
long advance_counter(long *counter, long *max_count, long n_indices)
Advances the counter array based on maximum counts.
Definition counter.c:51
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
ipow
double ipow(const double x, const int64_t p)
Compute x raised to the power p (x^p).
Definition ipow.c:33
get_argument_offset
double get_argument_offset()
Get the current argument offset applied before function evaluations.
Definition lsfBasisFns.c:53
get_argument_scale
double get_argument_scale()
Get the current argument scale factor used before function evaluations.
Definition lsfBasisFns.c:60
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
compute_average
long compute_average(double *value, double *data, int64_t n)
Computes the average of an array of doubles.
Definition median.c:144
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
set_argument_offset
void set_argument_offset(double offset)
Set the offset applied to the input argument of basis functions.
Definition lsfBasisFns.c:33
dtcheby
double dtcheby(double x, long n)
Evaluate the derivative of the Chebyshev polynomial T_n(x).
Definition lsfBasisFns.c:92
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
long_cmpasc
int long_cmpasc(const void *a, const void *b)
Compare two long integers in ascending order.
Definition sortfunctions.c:128
str_tolower
char * str_tolower(char *s)
Convert a string to lower case.
Definition str_tolower.c:27
CHEBYSHEV_COEF
Definition sddspfit.c:96
EVAL_PARAMETERS
Definition sddsmpfit.c:247
SDDS_DATASET
Definition SDDS.h:303
Generated by doxygen 1.12.0