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