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