SDDS ToolKit Programs and Libraries for C and Python
•All Classes Files Functions Variables Macros Pages
sddssplinefit.c
Go to the documentation of this file.
1/**
2 * @file sddssplinefit.c
3 * @brief Performs nth-order spline least squares fitting for SDDS files.
4 *
5 * @details
6 * This program fits splines to data contained in SDDS (Self Describing Data Sets) files. It allows for
7 * various configurations such as specifying the order of the spline, the number of coefficients or breakpoints,
8 * handling of sigma values, and more. The program processes SDDS files and outputs the fitted results
9 * along with optional diagnostics or evaluation.
10 *
11 * @section Usage
12 * ```
13 * sddssplinefit <inputfile> <outputfile>
14 * [-pipe=[input][,output]]
15 * -independent=<xName>
16 * -dependent=<yName1-wildcard>[,<yName2-wildcard>...]
17 * [-sigmaIndependent=<xSigma>]
18 * [-sigmaDependent=<ySigmaFormatString>]
19 * [-order=<number>]
20 * [-coefficients=<number>]
21 * [-breakpoints=<number>]
22 * [-xOffset=<value>]
23 * [-xFactor=<value>]
24 * [-sigmas=<value>,{absolute|fractional}]
25 * [-modifySigmas]
26 * [-generateSigmas[={keepLargest,keepSmallest}]]
27 * [-sparse=<interval>]
28 * [-range=<lower>,<upper>[,fitOnly]]
29 * [-normalize[=<termNumber>]]
30 * [-verbose]
31 * [-evaluate=<filename>[,begin=<value>][,end=<value>][,number=<integer>][,derivatives=<order>][,basis]]
32 * [-infoFile=<filename>]
33 * [-copyParameters]
34 * ```
35 *
36 * @section Options
37 * | Required | Description |
38 * |---------------------------------------|---------------------------------------------------------------------------------------|
39 * | `-independent` | Specifies the name of the independent data column to use. |
40 * | `-dependent` | Specifies the names of dependent data columns to use, supporting wildcards. |
41 *
42 * | Optional | Description |
43 * |---------------------------------------|---------------------------------------------------------------------------------------|
44 * | `-pipe` | Uses pipes for input/output data flow. |
45 * | `-sigmaIndependent` | Specifies the name of the independent sigma column. |
46 * | `-sigmaDependent` | Specifies a format string for dependent sigma column names (e.g., `%sSigma`). |
47 * | `-order` | Specifies the order of the spline. Default is 4. |
48 * | `-coefficients` | Defines the number of coefficients. |
49 * | `-breakpoints` | Sets the number of breakpoints. |
50 * | `-xOffset` | Specifies the desired value of x to fit about. |
51 * | `-xFactor` | Specifies the factor by which to multiply x values before fitting. |
52 * | `-sigmas` | Specifies absolute or fractional sigma for all points. |
53 * | `-modifySigmas` | Modifies the y sigmas using the x sigmas and an initial fit. |
54 * | `-generateSigmas` | Generates y sigmas from the RMS deviation of an initial fit. |
55 * | `-sparse` | Specifies an interval at which to sample data. |
56 * | `-range` | Defines the range of the independent variable for fitting and evaluation. |
57 * | `-normalize` | Normalizes so that the specified term is unity. |
58 * | `-verbose` | Enables verbose output for additional information. |
59 * | `-evaluate` | Evaluates the spline fit and optionally computes derivatives or basis functions. |
60 * | `-infoFile` | Specifies a file to output fit information. |
61 * | `-copyParameters` | Copies parameters from the input file to the output file. |
62 *
63 * @subsection Incompatibilities
64 * - `-coefficients` is incompatible with:
65 * - `-breakpoints`
66 * - For `-modifySigmas`:
67 * - Requires `-sigmaIndependent`
68 * - For `-evaluate`:
69 * - Parameters must include valid range and number of evaluation points.
70 *
71 * @copyright
72 * - (c) 2002 The University of Chicago, as Operator of Argonne National Laboratory.
73 * - (c) 2002 The Regents of the University of California, as Operator of Los Alamos National Laboratory.
74 *
75 * @license
76 * This file is distributed under the terms of the Software License Agreement
77 * found in the file LICENSE included with this distribution.
78 *
79 * @author
80 * L. Emery, M. Borland, R. Soliday
81 */
82
83#include "mdb.h"
84#include "SDDS.h"
85#include "SDDSutils.h"
86#include "scan.h"
87
88#include <gsl/gsl_bspline.h>
89#include <gsl/gsl_multifit.h>
90#include <gsl/gsl_rng.h>
91#include <gsl/gsl_randist.h>
92#include <gsl/gsl_statistics.h>
93#include <gsl/gsl_version.h>
94
95void makeSubstitutions(char *buffer1, char *buffer2, char *template, char *nameRoot, char *symbolRoot, char *xName, char *xSymbol);
96char *changeInformation(SDDS_DATASET *SDDSout, char *name, char *nameRoot, char *symbolRoot, char *xName, char *xSymbol, char **template, char *newUnits);
97char **makeCoefficientUnits(SDDS_DATASET *SDDSout, char *xName, char *yName, long int order, long coeffs);
98void initializeOutputFile(SDDS_DATASET *SDDSout, SDDS_DATASET *SDDSoutInfo, char *output, char *outputInfo,
99 SDDS_DATASET *SDDSin, char *input, char *xName, char **yNames, char *xSigmaName,
100 char **ySigmaNames, long sigmasValid, long int order, long coeffs, long breaks,
101 long numCols, long copyParameters);
102void checkInputFile(SDDS_DATASET *SDDSin, char *xName, char **yNames, char *xSigmaName, char **ySigmaNames, long numYNames);
103long coefficient_index(long int order, long coeffs, long order_of_interest);
104
105char **ResolveColumnNames(SDDS_DATASET *SDDSin, char **wildcardList, long length, int32_t *numYNames);
106char **GenerateYSigmaNames(char *controlString, char **yNames, long numYNames);
107void RemoveElementFromStringArray(char **array, long index, long length);
108char **RemoveNonNumericColumnsFromNameArray(SDDS_DATASET *SDDSin, char **columns, int32_t *numColumns);
109
110int print_matrix(FILE *f, const gsl_matrix *m);
111
112/* Enumeration for option types */
113enum option_type {
114 CLO_DEPENDENT,
115 CLO_ORDER,
116 CLO_COEFFICIENTS,
117 CLO_BREAKPOINTS,
118 CLO_REVISEORDERS,
119 CLO_XXXXX,
120 CLO_MODIFYSIGMAS,
121 CLO_SIGMAS,
122 CLO_GENERATESIGMAS,
123 CLO_RANGE,
124 CLO_SPARSE,
125 CLO_NORMALIZE,
126 CLO_XFACTOR,
127 CLO_XOFFSET,
128 CLO_VERBOSE,
129 CLO_PIPE,
130 CLO_EVALUATE,
131 CLO_INDEPENDENT,
132 CLO_SIGMAINDEPENDENT,
133 CLO_SIGMADEPENDENT,
134 CLO_INFOFILE,
135 CLO_COPYPARAMETERS,
136 N_OPTIONS
137};
138
139char *option[N_OPTIONS] = {
140 "dependent",
141 "order",
142 "coefficients",
143 "breakpoints",
144 "reviseorders",
145 "splinebasis",
146 "modifysigmas",
147 "sigmas",
148 "generatesigmas",
149 "range",
150 "sparse",
151 "normalize",
152 "xfactor",
153 "xoffset",
154 "verbose",
155 "pipe",
156 "evaluate",
157 "independent",
158 "sigmaindependent",
159 "sigmadependent",
160 "infofile",
161 "copyparameters",
162};
163
164char *USAGE =
165 "sddssplinefit [<inputfile>] [<outputfile>]\n"
166 " [-pipe=[input][,output]]\n"
167 " -independent=<xName>\n"
168 " -dependent=<yName1-wildcard>[,<yName2-wildcard>...]\n"
169 " [-sigmaIndependent=<xSigma>]\n"
170 " [-sigmaDependent=<ySigmaFormatString>]\n"
171 " [-order=<number>]\n"
172 " [-coefficients=<number>]\n"
173 " [-breakpoints=<number>]\n"
174 " [-xOffset=<value>]\n"
175 " [-xFactor=<value>]\n"
176 " [-sigmas=<value>,{absolute|fractional}]\n"
177 " [-modifySigmas]\n"
178 " [-generateSigmas[={keepLargest,keepSmallest}]]\n"
179 " [-sparse=<interval>]\n"
180 " [-range=<lower>,<upper>[,fitOnly]]\n"
181 " [-normalize[=<termNumber>]]\n"
182 " [-verbose]\n"
183 " [-evaluate=<filename>[,begin=<value>][,end=<value>][,number=<integer>][,derivatives=<order>][,basis]]\n"
184 " [-infoFile=<filename>]\n"
185 " [-copyParameters]\n"
186 "Program by Louis Emery, started with Michael Borland polynomial fit program (" __DATE__ " " __TIME__ ", SVN revision: " SVN_VERSION ")\n";
187
188static char *additional_help = "\n\
189sddssplinefit performs spline fits of the form y = SUM(i){ A[i] * B(x-x_offset, i)}, where B(x,i) is the ith basis\n\
190spline function evaluated at x. Internally, sddssplinefit computes the A[i] coefficients, writes the fitted y values to the output file,\n\
191and estimates the errors in these values.\n";
192
193static char *additional_help2 = "\n\
194 -independent Specify the name of the independent data column to use.\n\
195 -dependent Specify the names of dependent data columns to use, supporting wildcards and separated by commas.\n\
196 -sigmaIndependent Specify the name of the independent sigma values column.\n\
197 -sigmaDependent Specify a printf-style control string to generate dependent sigma column names from independent variable names (e.g., %sSigma).\n\
198 -order Define the order of the spline. Default is 4.\n\
199 -coefficients Set the number of coefficients. Specify either coefficients or breakpoints, not both.\n\
200 -breakpoints Set the number of breakpoints. Condition enforced: breakpoints = coefficients + 2 - order.\n\
201 -xOffset Define the desired value of x to fit about.\n\
202 -xFactor Define the factor to multiply x values by before fitting.\n\
203 -sigmas Specify absolute or fractional sigma for all points.\n\
204 -modifySigmas Modify the y sigmas using the x sigmas and an initial fit.\n\
205 -generateSigmas Generate y sigmas from the RMS deviation of an initial fit.\n\
206 Optionally keep the sigmas from the data if larger/smaller than the RMS deviation.\n\
207 -sparse Specify an integer interval at which to sample data.\n\
208 -range Define the range of the independent variable over which to perform the fit and evaluation.\n\
209 If 'fitOnly' is given, the fit is compared to data over the original range.\n\
210 -normalize Normalize so that the specified term is unity.\n\
211 -verbose Enable verbose output for additional information.\n\
212 -evaluate Evaluate the spline fit and optionally compute derivatives and provide basis functions.\n\
213 -infoFile Specify a file to output fit information.\n\
214 -copyParameters Copy parameters from the input file to the output file.\n\n";
215
216#define ABSOLUTE_SIGMAS 0
217#define FRACTIONAL_SIGMAS 1
218#define N_SIGMAS_OPTIONS 2
219char *sigmas_options[N_SIGMAS_OPTIONS] = {"absolute", "fractional"};
220
221#define FLGS_GENERATESIGMAS 1
222#define FLGS_KEEPLARGEST 2
223#define FLGS_KEEPSMALLEST 4
224
225#define REVPOW_ACTIVE 0x0001
226#define REVPOW_VERBOSE 0x0002
227/* SDDS indices for output page */
228static long iOffset = -1, iOffsetO = -1, iFactor = -1, iFactorO = -1;
229static long *iChiSq = NULL, *iChiSqO = NULL, *iRmsResidual = NULL, *iRmsResidualO = NULL, *iSigLevel = NULL, *iSigLevelO = NULL;
230static long *iFitIsValid = NULL, *iFitIsValidO = NULL;
231
232/* These column indices are globals, so we don't have to pass them to functions */
233static long ix = -1, ixSigma = -1;
234static long *iy = NULL, *iySigma = NULL;
235static long *iFit = NULL, *iResidual = NULL;
236
237static long *iCoefficient = NULL, *iCoefficientSigma = NULL, *iCoefficientUnits = NULL;
238
239static char *xSymbol, **ySymbols;
240
241#define EVAL_BEGIN_GIVEN 0x0001U
242#define EVAL_END_GIVEN 0x0002U
243#define EVAL_NUMBER_GIVEN 0x0004U
244#define EVAL_DERIVATIVES 0x0008U
245#define EVAL_PROVIDEBASIS 0x0010U
246
247#define MAX_Y_SIGMA_NAME_SIZE 1024
248
249typedef struct
250{
251 char *file;
252 long initialized;
253 int64_t number, nderiv;
254 unsigned long flags;
255 double begin, end;
256 SDDS_DATASET dataset;
257 char ***yDerivName, ***yDerivUnits; /* a matrix of names */
258 long *iSpline; /* index of spline columns for dataset */
259 gsl_bspline_workspace *bw; /* eases the passing of parameters in function arguments */
260} EVAL_PARAMETERS;
261void setupEvaluationFile(EVAL_PARAMETERS *evalParameters, char *xName, char **yName, long yNames, SDDS_DATASET *SDDSin);
262void makeEvaluationTable(EVAL_PARAMETERS *evalParameters, double *x, int64_t points, gsl_vector *B,
263 gsl_matrix *cov, gsl_vector *c, char *xName, char **yName, long yNames, long iYName, long int order);
264
265int main(int argc, char **argv) {
266 double **y = NULL, **yFit = NULL, **sy = NULL, **diff = NULL;
267 double *x = NULL, *sx = NULL;
268 double xOffset, xScaleFactor;
269 double *xOrig = NULL, **yOrig = NULL, **yFitOrig = NULL, *sxOrig = NULL, **syOrig = NULL, **sy0 = NULL;
270 long coeffs, breaks, normTerm, ySigmasValid;
271 long orderGiven, coeffsGiven, breaksGiven;
272 int64_t i, j, points, pointsOrig;
273 double sigmas;
274 long sigmasMode, sparseInterval;
275 double *chi = NULL, xLow, xHigh, *rmsResidual = NULL;
276 char *xName = NULL, *yName = NULL, **yNames = NULL, *xSigmaName = NULL;
277 char **ySigmaNames = NULL, *ySigmaControlString = NULL;
278 char *input = NULL, *output = NULL;
279 SDDS_DATASET SDDSin, SDDSout, SDDSoutInfo;
280 long *isFit = NULL, iArg, modifySigmas;
281 long generateSigmas, verbose, ignoreSigmas;
282 long outputInitialized, copyParameters = 0;
283 long int order;
284 SCANNED_ARG *s_arg;
285 double xMin, xMax, revpowThreshold;
286 double rms_average(double *d_x, int64_t d_n);
287 char *infoFile = NULL;
288 unsigned long pipeFlags, reviseOrders;
289 EVAL_PARAMETERS evalParameters;
290 long rangeFitOnly = 0;
291
292 long colIndex;
293 long cloDependentIndex = -1, numDependentItems;
294 int32_t numYNames;
295
296 /* spline memory for working on a single column */
297 gsl_bspline_workspace *bw;
298 gsl_vector *B;
299 gsl_vector *c, *yGsl, *wGsl;
300 gsl_matrix *X, *cov;
301 gsl_multifit_linear_workspace *mw;
302 long degreesOfFreedom;
303 double totalSumSquare, Rsq;
304
305 SDDS_RegisterProgramName(argv[0]);
306 argc = scanargs(&s_arg, argc, argv);
307 if (argc < 2 || argc > (3 + N_OPTIONS)) {
308 fprintf(stderr, "usage: %s\n", USAGE);
309 fprintf(stderr, "%s%s", additional_help, additional_help2);
310 exit(EXIT_FAILURE);
311 }
312
313 input = output = NULL;
314 xName = yName = xSigmaName = ySigmaControlString = NULL;
315 yNames = ySigmaNames = NULL;
316 numDependentItems = 0;
317 modifySigmas = reviseOrders = 0;
318 /* 8 independent coefficients, a default that I think is reasonable */
319 coeffs = 8;
320 xMin = xMax = 0;
321 generateSigmas = 0;
322 sigmasMode = -1;
323 sigmas = 1;
324 sparseInterval = 1;
325 verbose = ignoreSigmas = 0;
326 normTerm = -1;
327 xOffset = 0;
328 xScaleFactor = 1;
329 pipeFlags = 0;
330 evalParameters.file = NULL;
331 infoFile = NULL;
332 orderGiven = 0;
333 coeffsGiven = 0;
334 breaksGiven = 0;
335
336 for (iArg = 1; iArg < argc; iArg++) {
337 if (s_arg[iArg].arg_type == OPTION) {
338 switch (match_string(s_arg[iArg].list[0], option, N_OPTIONS, 0)) {
339 case CLO_MODIFYSIGMAS:
340 modifySigmas = 1;
341 break;
342 case CLO_ORDER:
343 if (s_arg[iArg].n_items != 2 || sscanf(s_arg[iArg].list[1], "%ld", &order) != 1)
344 SDDS_Bomb("invalid -order syntax");
345 orderGiven = 1;
346 break;
347 case CLO_COEFFICIENTS:
348 if (s_arg[iArg].n_items != 2 || sscanf(s_arg[iArg].list[1], "%ld", &coeffs) != 1)
349 SDDS_Bomb("invalid -coefficients syntax");
350 coeffsGiven = 1;
351 break;
352 case CLO_BREAKPOINTS:
353 if (s_arg[iArg].n_items != 2 || sscanf(s_arg[iArg].list[1], "%ld", &breaks) != 1)
354 SDDS_Bomb("invalid -breakpoints syntax");
355 breaksGiven = 1;
356 break;
357 case CLO_RANGE:
358 rangeFitOnly = 0;
359 if ((s_arg[iArg].n_items != 3 && s_arg[iArg].n_items != 4) ||
360 1 != sscanf(s_arg[iArg].list[1], "%lf", &xMin) ||
361 1 != sscanf(s_arg[iArg].list[2], "%lf", &xMax) ||
362 xMin >= xMax)
363 SDDS_Bomb("incorrect -range syntax");
364 if (s_arg[iArg].n_items == 4) {
365 if (strncmp(str_tolower(s_arg[iArg].list[3]), "fitonly", strlen(s_arg[iArg].list[3])) == 0) {
366 rangeFitOnly = 1;
367 } else
368 SDDS_Bomb("incorrect -range syntax");
369 }
370 break;
371 case CLO_GENERATESIGMAS:
372 generateSigmas = FLGS_GENERATESIGMAS;
373 if (s_arg[iArg].n_items > 1) {
374 if (s_arg[iArg].n_items != 2)
375 SDDS_Bomb("incorrect -generateSigmas syntax");
376 if (strncmp(s_arg[iArg].list[1], "keepsmallest", strlen(s_arg[iArg].list[1])) == 0)
377 generateSigmas |= FLGS_KEEPSMALLEST;
378 if (strncmp(s_arg[iArg].list[1], "keeplargest", strlen(s_arg[iArg].list[1])) == 0)
379 generateSigmas |= FLGS_KEEPLARGEST;
380 if ((generateSigmas & FLGS_KEEPSMALLEST) && (generateSigmas & FLGS_KEEPLARGEST))
381 SDDS_Bomb("ambiguous -generateSigmas syntax");
382 }
383 break;
384 case CLO_XOFFSET:
385 if (s_arg[iArg].n_items != 2 || sscanf(s_arg[iArg].list[1], "%lf", &xOffset) != 1)
386 SDDS_Bomb("invalid -xOffset syntax");
387 break;
388 case CLO_SIGMAS:
389 if (s_arg[iArg].n_items != 3)
390 SDDS_Bomb("incorrect -sigmas syntax");
391 if (sscanf(s_arg[iArg].list[1], "%lf", &sigmas) != 1)
392 SDDS_Bomb("couldn't scan value for -sigmas");
393 if ((sigmasMode = match_string(s_arg[iArg].list[2], sigmas_options, N_SIGMAS_OPTIONS, 0)) < 0)
394 SDDS_Bomb("unrecognized -sigmas mode");
395 break;
396 case CLO_SPARSE:
397 if (s_arg[iArg].n_items != 2)
398 SDDS_Bomb("incorrect -sparse syntax");
399 if (sscanf(s_arg[iArg].list[1], "%ld", &sparseInterval) != 1)
400 SDDS_Bomb("couldn't scan value for -sparse");
401 if (sparseInterval < 1)
402 SDDS_Bomb("invalid -sparse value");
403 break;
404 case CLO_VERBOSE:
405 verbose = 1;
406 break;
407 case CLO_NORMALIZE:
408 normTerm = 0;
409 if (s_arg[iArg].n_items > 2 ||
410 (s_arg[iArg].n_items == 2 && sscanf(s_arg[iArg].list[1], "%ld", &normTerm) != 1) ||
411 normTerm < 0)
412 SDDS_Bomb("invalid -normalize syntax");
413 break;
414 case CLO_REVISEORDERS:
415 revpowThreshold = 0.1;
416 s_arg[iArg].n_items -= 1;
417 if (!scanItemList(&reviseOrders, s_arg[iArg].list + 1, &s_arg[iArg].n_items, 0,
418 "threshold", SDDS_DOUBLE, &revpowThreshold, 1, 0,
419 "verbose", -1, NULL, 1, REVPOW_VERBOSE, NULL))
420 SDDS_Bomb("invalid -reviseOrders syntax");
421 s_arg[iArg].n_items += 1;
422 reviseOrders |= REVPOW_ACTIVE;
423 revpowThreshold = fabs(revpowThreshold);
424 break;
425 case CLO_XFACTOR:
426 if (s_arg[iArg].n_items != 2 ||
427 sscanf(s_arg[iArg].list[1], "%lf", &xScaleFactor) != 1 || xScaleFactor == 0)
428 SDDS_Bomb("invalid -xFactor syntax");
429 break;
430 case CLO_INDEPENDENT:
431 if (s_arg[iArg].n_items != 2)
432 SDDS_Bomb("invalid -independent syntax");
433 xName = s_arg[iArg].list[1];
434 break;
435 case CLO_DEPENDENT:
436 numDependentItems = s_arg[iArg].n_items - 1;
437 cloDependentIndex = iArg;
438 if (numDependentItems < 1)
439 SDDS_Bomb("invalid -dependent syntax");
440 break;
441 case CLO_SIGMAINDEPENDENT:
442 if (s_arg[iArg].n_items != 2)
443 SDDS_Bomb("invalid -sigmaIndependent syntax");
444 xSigmaName = s_arg[iArg].list[1];
445 break;
446 case CLO_SIGMADEPENDENT:
447 if (s_arg[iArg].n_items != 2)
448 SDDS_Bomb("invalid -sigmaDependent syntax");
449 ySigmaControlString = s_arg[iArg].list[1];
450 break;
451 case CLO_PIPE:
452 if (!processPipeOption(s_arg[iArg].list + 1, s_arg[iArg].n_items - 1, &pipeFlags))
453 SDDS_Bomb("invalid -pipe syntax");
454 break;
455 case CLO_INFOFILE:
456 if (s_arg[iArg].n_items != 2)
457 SDDS_Bomb("invalid -infoFile syntax");
458 infoFile = s_arg[iArg].list[1];
459 break;
460 case CLO_EVALUATE:
461 if (s_arg[iArg].n_items < 2)
462 SDDS_Bomb("invalid -evaluate syntax");
463 evalParameters.file = s_arg[iArg].list[1];
464 s_arg[iArg].n_items -= 2;
465 s_arg[iArg].list += 2;
466 evalParameters.begin = 0;
467 evalParameters.end = 0;
468 evalParameters.nderiv = 0;
469 evalParameters.number = 0;
470 if (!scanItemList(&evalParameters.flags, s_arg[iArg].list, &s_arg[iArg].n_items, 0,
471 "begin", SDDS_DOUBLE, &evalParameters.begin, 1, EVAL_BEGIN_GIVEN,
472 "end", SDDS_DOUBLE, &evalParameters.end, 1, EVAL_END_GIVEN,
473 "derivatives", SDDS_LONG64, &evalParameters.nderiv, 1, EVAL_DERIVATIVES,
474 "basis", -1, NULL, 0, EVAL_PROVIDEBASIS,
475 "number", SDDS_LONG64, &evalParameters.number, 1, EVAL_NUMBER_GIVEN, NULL))
476 SDDS_Bomb("invalid -evaluate syntax");
477 s_arg[iArg].n_items += 2;
478 s_arg[iArg].list -= 2;
479 break;
480 case CLO_COPYPARAMETERS:
481 copyParameters = 1;
482 break;
483 default:
484 bomb("unknown switch", USAGE);
485 break;
486 }
487 } else {
488 if (input == NULL)
489 input = s_arg[iArg].list[0];
490 else if (output == NULL)
491 output = s_arg[iArg].list[0];
492 else
493 SDDS_Bomb("too many filenames");
494 }
495 }
496 /* your basic spline is order 4 (continuous second derivative) */
497 if (!orderGiven)
498 order = 4;
499 if (!breaksGiven)
500 breaks = coeffs + 2 - order;
501 if (!coeffsGiven)
502 coeffs = breaks - 2 + order;
503 if (breaksGiven && coeffsGiven)
504 SDDS_Bomb("You must specify only one of breakpoints or coefficients");
505
506 processFilenames("sddssplinefit", &input, &output, pipeFlags, 0, NULL);
507
508 if (!xName || !numDependentItems)
509 SDDS_Bomb("you must specify a column name for x and y");
510 if (modifySigmas && !xSigmaName)
511 SDDS_Bomb("you must specify x sigmas with -modifySigmas");
512 if (generateSigmas) {
513 if (modifySigmas)
514 SDDS_Bomb("you can't specify both -generateSigmas and -modifySigmas");
515 }
516 if (ySigmaControlString) {
517 if (sigmasMode != -1)
518 SDDS_Bomb("you can't specify both -sigmas and a y sigma name");
519 }
520 ySigmasValid = 0;
521 if (sigmasMode != -1 || generateSigmas || ySigmaControlString || modifySigmas)
522 ySigmasValid = 1;
523
524 if (!SDDS_InitializeInput(&SDDSin, input))
525 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
526 outputInitialized = 0;
527 yNames = ResolveColumnNames(&SDDSin, s_arg[cloDependentIndex].list + 1, numDependentItems, &numYNames);
528 if (ySigmaControlString != NULL)
529 ySigmaNames = GenerateYSigmaNames(ySigmaControlString, yNames, numYNames);
530
531 checkInputFile(&SDDSin, xName, yNames, xSigmaName, ySigmaNames, numYNames);
532 sy0 = tmalloc(sizeof(double *) * numYNames);
533 y = tmalloc(sizeof(double *) * numYNames);
534 yFit = tmalloc(sizeof(double *) * numYNames); /* this array of arrays is not in sddsmpfit */
535 sy = tmalloc(sizeof(double *) * numYNames);
536 isFit = tmalloc(sizeof(long) * numYNames);
537 chi = tmalloc(sizeof(double) * numYNames);
538 iCoefficient = tmalloc(sizeof(long) * numYNames);
539 iCoefficientSigma = tmalloc(sizeof(long) * numYNames);
540 iCoefficientUnits = tmalloc(sizeof(long) * numYNames);
541
542 while (SDDS_ReadPage(&SDDSin) > 0) {
543 if ((points = SDDS_CountRowsOfInterest(&SDDSin)) < coeffs) {
544 /* probably should emit an empty page here */
545 continue;
546 }
547 if (verbose)
548 fprintf(stdout, "number of points %" PRId64 "\n", points);
549 if (!(x = SDDS_GetColumnInDoubles(&SDDSin, xName))) {
550 fprintf(stderr, "error: unable to read column %s\n", xName);
551 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
552 }
553 for (i = 0; i < numYNames; i++) {
554 if (!(y[i] = SDDS_GetColumnInDoubles(&SDDSin, yNames[i]))) {
555 fprintf(stderr, "error: unable to read column %s\n", yNames[i]);
556 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
557 }
558 }
559 sx = NULL;
560 if (xSigmaName && !(sx = SDDS_GetColumnInDoubles(&SDDSin, xSigmaName))) {
561 fprintf(stderr, "error: unable to read column %s\n", xSigmaName);
562 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
563 }
564 for (colIndex = 0; colIndex < numYNames; colIndex++) {
565 sy0[colIndex] = tmalloc(sizeof(double) * points);
566 yFit[colIndex] = tmalloc(sizeof(double) * points);
567 }
568
569 if (ySigmaNames) {
570 for (i = 0; i < numYNames; i++) {
571 if (!(sy0[i] = SDDS_GetColumnInDoubles(&SDDSin, ySigmaNames[i]))) {
572 fprintf(stderr, "error: unable to read column %s\n", ySigmaNames[i]);
573 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
574 }
575 }
576 }
577
578 if (xMin != xMax || sparseInterval != 1) {
579 xOrig = tmalloc(sizeof(*xOrig) * points);
580 yOrig = tmalloc(sizeof(*yOrig) * numYNames);
581 yFitOrig = tmalloc(sizeof(*yFitOrig) * numYNames);
582 for (colIndex = 0; colIndex < numYNames; colIndex++) {
583 if (verbose)
584 fprintf(stdout, "Setting up a separate array for range or sparsing for column %s because of range option ...\n", yNames[colIndex]);
585 yOrig[colIndex] = tmalloc(sizeof(double) * points);
586 yFitOrig[colIndex] = tmalloc(sizeof(double) * points);
587 }
588 if (sx)
589 sxOrig = tmalloc(sizeof(*sxOrig) * points);
590 if (ySigmasValid) {
591 syOrig = tmalloc(sizeof(*syOrig) * numYNames);
592 for (colIndex = 0; colIndex < numYNames; colIndex++)
593 syOrig[colIndex] = tmalloc(sizeof(double) * points);
594 }
595 pointsOrig = points;
596 for (i = j = 0; i < points; i++) {
597 xOrig[i] = x[i];
598 if (sx)
599 sxOrig[i] = sx[i];
600 for (colIndex = 0; colIndex < numYNames; colIndex++) {
601 yOrig[colIndex][i] = y[colIndex][i];
602 if (ySigmasValid)
603 syOrig[colIndex][i] = sy0[colIndex][i];
604 }
605 }
606 if (xMin != xMax) {
607 for (i = j = 0; i < points; i++) {
608 if (xOrig[i] <= xMax && xOrig[i] >= xMin) {
609 x[j] = xOrig[i];
610 for (colIndex = 0; colIndex < numYNames; colIndex++) {
611 y[colIndex][j] = yOrig[colIndex][i];
612 if (ySigmasValid)
613 sy0[colIndex][j] = syOrig[colIndex][i];
614 }
615 if (sx)
616 sx[j] = sxOrig[i];
617 j++;
618 }
619 }
620 points = j;
621 }
622 if (sparseInterval != 1) {
623 for (i = j = 0; i < points; i++) {
624 if (i % sparseInterval == 0) {
625 x[j] = x[i];
626 for (colIndex = 0; colIndex < numYNames; colIndex++) {
627 y[colIndex][j] = y[colIndex][i];
628 if (ySigmasValid)
629 sy0[colIndex][j] = sy0[colIndex][i];
630 }
631 if (sx)
632 sx[j] = sx[i];
633 j++;
634 }
635 }
636 points = j;
637 }
638 } else {
639 /* normal processing, no ranges or sparsing */
640 xOrig = x;
641 yOrig = y;
642 sxOrig = sx;
643 syOrig = sy0;
644 pointsOrig = points;
645 }
646
647 find_min_max(&xLow, &xHigh, x, points);
648 if (verbose)
649 fprintf(stdout, "Range: xLow %lf; xHigh %lf; points %" PRId64 "\n", xLow, xHigh, points);
650 if (sigmasMode == ABSOLUTE_SIGMAS) {
651 for (colIndex = 0; colIndex < numYNames; colIndex++) {
652 for (i = 0; i < points; i++)
653 sy0[colIndex][i] = sigmas;
654 if (sy0[colIndex] != syOrig[colIndex])
655 for (i = 0; i < pointsOrig; i++)
656 syOrig[colIndex][i] = sigmas;
657 }
658 } else if (sigmasMode == FRACTIONAL_SIGMAS) {
659 for (colIndex = 0; colIndex < numYNames; colIndex++) {
660 for (i = 0; i < points; i++)
661 sy0[colIndex][i] = sigmas * fabs(y[colIndex][i]);
662 if (sy0[colIndex] != syOrig[colIndex])
663 for (i = 0; i < pointsOrig; i++)
664 syOrig[colIndex][i] = fabs(yOrig[colIndex][i]) * sigmas;
665 }
666 }
667
668 if (!ySigmasValid || generateSigmas)
669 for (colIndex = 0; colIndex < numYNames; colIndex++) {
670 for (i = 0; i < points; i++)
671 sy0[colIndex][i] = 1;
672 }
673 else
674 for (i = 0; i < points; i++)
675 for (colIndex = 0; colIndex < numYNames; colIndex++) {
676 if (sy0[colIndex][i] == 0)
677 SDDS_Bomb("y sigma = 0 for one or more points.");
678 }
679
680 diff = tmalloc(sizeof(*diff) * numYNames);
681 sy = tmalloc(sizeof(*sy) * numYNames);
682 for (colIndex = 0; colIndex < numYNames; colIndex++) {
683 diff[colIndex] = tmalloc(sizeof(double) * points);
684 sy[colIndex] = tmalloc(sizeof(double) * points);
685 }
686
687 for (i = 0; i < points; i++) {
688 for (colIndex = 0; colIndex < numYNames; colIndex++)
689 sy[colIndex][i] = sy0[colIndex][i];
690 }
691
692 /* this seems the places when things really start */
693
694 /* allocate a cubic bspline workspace (k = 4) */
695 /* k is order of spline; cubic infers k=4; breaks are number of splines */
696 /* each will be reused for each of the columns. No need to make them an array */
697 bw = gsl_bspline_alloc(order, breaks);
698 B = gsl_vector_alloc(coeffs); /* coeffs are the number of linear,
699 quadratic and cubic coeff on the
700 splines, say for k=4 */
701 X = gsl_matrix_alloc(points, coeffs);
702 c = gsl_vector_alloc(coeffs);
703 yGsl = gsl_vector_alloc(points);
704 wGsl = gsl_vector_alloc(points);
705 cov = gsl_matrix_alloc(coeffs, coeffs);
706 mw = gsl_multifit_linear_alloc(points, coeffs);
707 degreesOfFreedom = points - coeffs;
708 if (verbose)
709 fprintf(stdout, "Order %ld\ncoefficients %ld\nbreak points %ld\n", order, coeffs, breaks);
710 if (generateSigmas || modifySigmas)
711 fprintf(stderr, "generate sigmas or modify sigmas are not a feature in spline fitting.\n");
712
713 if (reviseOrders & REVPOW_ACTIVE)
714 fprintf(stderr, "revise orders is not a feature in spline fitting.\n");
715
716 if (!outputInitialized) {
717 initializeOutputFile(&SDDSout, &SDDSoutInfo, output, infoFile, &SDDSin, input, xName, yNames,
718 xSigmaName, ySigmaNames, ySigmasValid, order, coeffs, breaks, numYNames, copyParameters);
719 // free(output);
720 outputInitialized = 1;
721 /* we also want to setup the evaluation file only once */
722 if (evalParameters.file) {
723 /* check nderiv against order, repair if necessary */
724 if (evalParameters.nderiv >= order) {
725 evalParameters.nderiv = order - 1;
726 if (verbose)
727 fprintf(stderr, "Spline derivative order reduced to %" PRId64 " (i.e. order - 1)\n", evalParameters.nderiv);
728 }
729 evalParameters.bw = bw;
730 setupEvaluationFile(&evalParameters, xName, yNames, numYNames, &SDDSin);
731 }
732 }
733
734 rmsResidual = tmalloc(sizeof(double) * numYNames);
735
736 if (outputInitialized) {
737 if (!SDDS_StartPage(&SDDSout, rangeFitOnly ? pointsOrig : points) ||
738 (infoFile && !SDDS_StartPage(&SDDSoutInfo, coeffs)))
739 bomb("A", NULL);
740 if (copyParameters) {
741 if (!SDDS_CopyParameters(&SDDSout, &SDDSin))
742 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
743 if (infoFile && !SDDS_CopyParameters(&SDDSoutInfo, &SDDSin))
744 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
745 }
746 /* info file writing has been removed below for now. See sddsmpfit.c to retrieve the original */
747 if (!SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_INDEX, rangeFitOnly ? xOrig : x, rangeFitOnly ? pointsOrig : points, ix))
748 bomb("B", NULL);
749 for (colIndex = 0; colIndex < numYNames; colIndex++) {
750 /* do fit now for each column */
751 double Bj;
752 for (i = 0; i < points; i++) {
753 gsl_vector_set(yGsl, i, y[colIndex][i]);
754 /* if there was no sigmaY data given then sy = 1 will be used */
755 gsl_vector_set(wGsl, i, 1.0 / ipower(sy[colIndex][i], 2));
756 }
757 /* use uniform breakpoints on [low, high] */
758 /* what if we don't want uniform breakpoints? */
759 gsl_bspline_knots_uniform(xLow, xHigh, bw);
760 /* alternative is gsl_bspline_knots ( breakpts, bw) where
761 breakpts is the gsl vector of breakpoints that is
762 supplied by the user */
763
764 /* construct the fit matrix X */
765 for (i = 0; i < points; ++i) {
766 /* compute B_j(xi) for all j */
767 /* B is a vector */
768 gsl_bspline_eval(x[i], B, bw);
769 /* fill in row i of X */
770 for (j = 0; j < coeffs; ++j) {
771 Bj = gsl_vector_get(B, j);
772 /* X is some sort of large matrix points x coeffs, e.g. 100 x 8 */
773 gsl_matrix_set(X, i, j, Bj);
774 }
775 }
776 /* show the matrix X */
777 if (verbose == 2) {
778 fprintf(stderr, "X matrix %s:\n", yNames[colIndex]);
779 print_matrix(stderr, X);
780 }
781 /* do the fit */
782 gsl_multifit_wlinear(X, wGsl, yGsl, c, cov, &chi[colIndex], mw);
783 if (verbose)
784 fprintf(stdout, "conventionally-defined chi = sum( sqr(diff) * weight): %e\n", chi[colIndex]);
785 /* c is the answer */
786 if (verbose == 2) {
787 fprintf(stderr, "Covariance matrix for %s:\n", yNames[colIndex]);
788 print_matrix(stderr, cov);
789 }
790 /* weighted total sum of squares */
791 totalSumSquare = gsl_stats_wtss(wGsl->data, 1, yGsl->data, 1, yGsl->size);
792 Rsq = 1.0 - chi[colIndex] / totalSumSquare;
793 if (verbose)
794 fprintf(stdout, "(reduced) chisq/dof = %e, Rsq = %f\n", chi[colIndex] / degreesOfFreedom, Rsq);
795
796 double y_err; /* standard deviation of model at x[i]. Is it useful? */
797 for (i = 0; i < points; i++) {
798 gsl_bspline_eval(x[i], B, bw);
799 /* y = B c is the evaluated function. B must be interesting looking. */
800 gsl_multifit_linear_est(B, c, cov, &yFit[colIndex][i], &y_err);
801 }
802 if (rangeFitOnly) {
803 for (i = 0; i < pointsOrig; i++) {
804 diff[colIndex][i] = yOrig[colIndex][i] - yFitOrig[colIndex][i];
805 }
806 rmsResidual[colIndex] = rms_average(diff[colIndex], points);
807 /* the index iFit refers to the fit data column corresponding to colIndex */
808 if (!SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_INDEX, yOrig[colIndex], pointsOrig, iy[colIndex]) ||
809 !SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_INDEX, yFitOrig[colIndex], pointsOrig, iFit[colIndex]) ||
810 !SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_INDEX, diff[colIndex], pointsOrig, iResidual[colIndex]))
811 bomb("C", NULL);
812 } else {
813 for (i = 0; i < points; i++) {
814 diff[colIndex][i] = y[colIndex][i] - yFit[colIndex][i];
815 }
816 rmsResidual[colIndex] = rms_average(diff[colIndex], points);
817 /* the index iFit refers to the fit data column corresponding to colIndex */
818 if (!SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_INDEX, y[colIndex], points, iy[colIndex]) ||
819 !SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_INDEX, yFit[colIndex], points, iFit[colIndex]) ||
820 !SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_INDEX, diff[colIndex], points, iResidual[colIndex]))
821 bomb("C", NULL);
822 }
823 if (infoFile)
824 if (!SDDS_SetColumnFromDoubles(&SDDSoutInfo, SDDS_SET_BY_INDEX, c->data, coeffs, iCoefficient[colIndex]))
825 SDDS_PrintErrors(stderr, SDDS_EXIT_PrintErrors | SDDS_VERBOSE_PrintErrors);
826 if (evalParameters.file) {
827 evalParameters.bw = bw;
828 makeEvaluationTable(&evalParameters, x, points, B, cov, c, xName, yNames, numYNames, colIndex, order);
829 }
830 }
831
832 if (ixSigma != -1 &&
833 !SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_INDEX, rangeFitOnly ? sxOrig : sx, rangeFitOnly ? pointsOrig : points, ixSigma))
834 bomb("E", NULL);
835 if (infoFile) {
836 {
837 int32_t *Indices;
838 Indices = malloc(sizeof(*Indices) * coeffs);
839 for (i = 0; i < coeffs; i++)
840 Indices[i] = i;
841 if (!SDDS_SetColumn(&SDDSoutInfo, SDDS_SET_BY_NAME, Indices, coeffs, "Index"))
842 SDDS_PrintErrors(stderr, SDDS_EXIT_PrintErrors | SDDS_VERBOSE_PrintErrors);
843 }
844 }
845 for (colIndex = 0; colIndex < numYNames; colIndex++) {
846 if (ySigmasValid && iySigma[colIndex] != -1 &&
847 !SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_INDEX, rangeFitOnly ? syOrig[colIndex] : sy[colIndex],
848 rangeFitOnly ? pointsOrig : points, iySigma[colIndex]))
849 bomb("F", NULL);
850
851 /* info file has been removed for now. Splines have orders
852 but not used the same way as for regular least
853 squares of polynomials */
854 if (infoFile) {
855 if (!SDDS_SetParameters(&SDDSoutInfo, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE,
856 iRmsResidual[colIndex], rmsResidual[colIndex],
857 iChiSq[colIndex], (chi[colIndex] / degreesOfFreedom),
858 iSigLevel[colIndex], ChiSqrSigLevel(chi[colIndex], points - coeffs),
859 iOffset, xOffset, iFactor, xScaleFactor, iFitIsValid[colIndex], isFit[colIndex] ? 'y' : 'n', -1))
860 bomb("O", NULL);
861 }
862
863 /* writing the results for each page */
864 if (!SDDS_SetParameters(&SDDSout, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE,
865 iRmsResidualO[colIndex], rmsResidual[colIndex],
866 iChiSqO[colIndex], (chi[colIndex] / degreesOfFreedom),
867 iSigLevelO[colIndex], ChiSqrSigLevel(chi[colIndex], points - coeffs),
868 iOffsetO, xOffset, iFactorO, xScaleFactor, iFitIsValidO[colIndex], isFit[colIndex] ? 'y' : 'n', -1))
869 bomb("O", NULL);
870 }
871 if (!SDDS_WritePage(&SDDSout) || (infoFile && !SDDS_WritePage(&SDDSoutInfo)))
872 bomb("O", NULL);
873 }
874
875 /* this is the end of the page, the wrap-up before going to the next page. */
876 if (xOrig != x)
877 free(xOrig);
878 if (sxOrig != sx)
879 free(sxOrig);
880 free(x);
881 free(sx);
882 for (colIndex = 0; colIndex < numYNames; colIndex++) {
883 free(diff[colIndex]);
884 free(sy[colIndex]);
885 if (yOrig[colIndex] != y[colIndex])
886 free(yOrig[colIndex]);
887 if (syOrig && sy0 && syOrig[colIndex] != sy0[colIndex])
888 free(syOrig[colIndex]);
889 free(y[colIndex]);
890 if (sy0 && sy0[colIndex])
891 free(sy0[colIndex]);
892 if (yFit && yFit[colIndex])
893 free(yFit[colIndex]);
894 }
895 gsl_bspline_free(bw);
896 gsl_vector_free(B);
897 gsl_matrix_free(X);
898 gsl_vector_free(yGsl);
899 gsl_vector_free(wGsl);
900 // gsl_vector_alloc(c);
901 gsl_matrix_free(cov);
902 gsl_multifit_linear_free(mw);
903 }
904
905 if (yFit)
906 free(yFit);
907 if (outputInitialized) {
908 if (!SDDS_Terminate(&SDDSout)) {
909 SDDS_PrintErrors(stderr, SDDS_EXIT_PrintErrors | SDDS_VERBOSE_PrintErrors);
910 }
911 if (infoFile) {
912 if (!SDDS_Terminate(&SDDSoutInfo)) {
913 SDDS_PrintErrors(stderr, SDDS_EXIT_PrintErrors | SDDS_VERBOSE_PrintErrors);
914 }
915 }
916 if (evalParameters.file) {
917 if (!SDDS_Terminate(&evalParameters.dataset)) {
918 SDDS_PrintErrors(stderr, SDDS_EXIT_PrintErrors | SDDS_VERBOSE_PrintErrors);
919 }
920 }
921 }
922 if (!SDDS_Terminate(&SDDSin)) {
923 SDDS_PrintErrors(stderr, SDDS_EXIT_PrintErrors | SDDS_VERBOSE_PrintErrors);
924 }
925 free_scanargs(&s_arg, argc);
926
927 return EXIT_SUCCESS;
928}
929
930void RemoveElementFromStringArray(char **array, long index, long length) {
931 long lh;
932
933 for (lh = index; lh < length - 1; lh++)
934 array[lh] = array[lh + 1];
935}
936
937char **RemoveNonNumericColumnsFromNameArray(SDDS_DATASET *SDDSin, char **columns, int32_t *numColumns) {
938 long i, numNumericColumns = *numColumns;
939
940 for (i = 0; i < *numColumns; i++) {
941 if (SDDS_CheckColumn(SDDSin, columns[i], NULL, SDDS_ANY_NUMERIC_TYPE, NULL)) {
942 printf("Removing %s because not a numeric type.", columns[i]);
943 RemoveElementFromStringArray(columns, i, *numColumns);
944 numNumericColumns--;
945 }
946 }
947
948 *numColumns = numNumericColumns;
949 return (columns);
950}
951
952char **ResolveColumnNames(SDDS_DATASET *SDDSin, char **wildcardList, long length, int32_t *numYNames) {
953 char **result;
954 long i;
955
956 /* initially set the columns of interest to none, to make SDDS_OR work below */
957 SDDS_SetColumnsOfInterest(SDDSin, SDDS_MATCH_STRING, "", SDDS_AND);
958 for (i = 0; i < length; i++) {
959 SDDS_SetColumnsOfInterest(SDDSin, SDDS_MATCH_STRING, wildcardList[i], SDDS_OR);
960 }
961
962 if (!(result = SDDS_GetColumnNames(SDDSin, numYNames)) || *numYNames == 0)
963 bomb("Error matching columns in ResolveColumnNames: No matches.", NULL);
964
965 result = RemoveNonNumericColumnsFromNameArray(SDDSin, result, numYNames);
966 return (result);
967}
968
969char **GenerateYSigmaNames(char *controlString, char **yNames, long numYNames) {
970 long i, nameLength;
971 char **result, sigmaName[MAX_Y_SIGMA_NAME_SIZE];
972
973 result = tmalloc(sizeof(char *) * numYNames);
974 for (i = 0; i < numYNames; i++) {
975 sprintf(sigmaName, controlString, yNames[i]);
976 nameLength = strlen(sigmaName);
977 result[i] = tmalloc(sizeof(char) * (nameLength + 1));
978 strcpy(result[i], sigmaName);
979 }
980 return (result);
981}
982
983double rms_average(double *x, int64_t n) {
984 double sum2;
985 int64_t i;
986
987 for (i = sum2 = 0; i < n; i++)
988 sum2 += sqr(x[i]);
989
990 return (sqrt(sum2 / n));
991}
992
993void checkInputFile(SDDS_DATASET *SDDSin, char *xName, char **yNames, char *xSigmaName, char **ySigmaNames, long numYNames) {
994 char *ptr = NULL;
995 long i;
996
997 if (!(ptr = SDDS_FindColumn(SDDSin, FIND_NUMERIC_TYPE, xName, NULL)))
998 SDDS_Bomb("x column doesn't exist or is nonnumeric");
999 free(ptr);
1000
1001 /* y columns don't need to be checked because located using SDDS_SetColumnsOfInterest */
1002
1003 ptr = NULL;
1004 if (xSigmaName && !(ptr = SDDS_FindColumn(SDDSin, FIND_NUMERIC_TYPE, xSigmaName, NULL)))
1005 SDDS_Bomb("x sigma column doesn't exist or is nonnumeric");
1006 if (ptr)
1007 free(ptr);
1008
1009 if (ySigmaNames) {
1010 for (i = 0; i < numYNames; i++) {
1011 ptr = NULL;
1012 if (!(ptr = SDDS_FindColumn(SDDSin, FIND_NUMERIC_TYPE, ySigmaNames[i], NULL)))
1013 SDDS_Bomb("y sigma column doesn't exist or is nonnumeric");
1014 if (ptr)
1015 free(ptr);
1016 }
1017 }
1018}
1019
1020void initializeOutputFile(SDDS_DATASET *SDDSout, SDDS_DATASET *SDDSoutInfo, char *output, char *outputInfo,
1021 SDDS_DATASET *SDDSin, char *input, char *xName, char **yNames, char *xSigmaName,
1022 char **ySigmaNames, long sigmasValid, long int order, long coeffs, long breakpoints,
1023 long numCols, long copyParameters) {
1024 char buffer[SDDS_MAXLINE], buffer1[SDDS_MAXLINE], buffer2[SDDS_MAXLINE], buffer3[SDDS_MAXLINE];
1025 char *xUnits, *yUnits;
1026 // char ***coefUnits;
1027 long colIndex;
1028
1029 /* all array names followed by an 'O' contain the index of the parameter in the main output file; others refer to
1030 parameters in the infoFile */
1031 // coefUnits = tmalloc(sizeof(char **) * numCols);
1032 ySymbols = tmalloc(sizeof(char *) * numCols);
1033 iChiSq = tmalloc(sizeof(long) * numCols);
1034 iChiSqO = tmalloc(sizeof(long) * numCols);
1035 iRmsResidual = tmalloc(sizeof(long) * numCols);
1036 iRmsResidualO = tmalloc(sizeof(long) * numCols);
1037 iSigLevel = tmalloc(sizeof(long) * numCols);
1038 iSigLevelO = tmalloc(sizeof(long) * numCols);
1039 iFitIsValid = tmalloc(sizeof(long) * numCols);
1040 iFitIsValidO = tmalloc(sizeof(long) * numCols);
1041 iy = tmalloc(sizeof(long) * numCols);
1042 iySigma = tmalloc(sizeof(long) * numCols);
1043 iFit = tmalloc(sizeof(long) * numCols);
1044 iResidual = tmalloc(sizeof(long) * numCols);
1045
1046 for (colIndex = 0; colIndex < numCols; colIndex++) {
1047 ySymbols[colIndex] = NULL;
1048 // coefUnits[colIndex] = tmalloc(sizeof(char *) * coeffs);
1049 iChiSq[colIndex] = -1;
1050 iChiSqO[colIndex] = -1;
1051 iRmsResidual[colIndex] = -1;
1052 iRmsResidualO[colIndex] = -1;
1053 iSigLevel[colIndex] = -1;
1054 iSigLevelO[colIndex] = -1;
1055 iFitIsValid[colIndex] = -1;
1056 iFitIsValidO[colIndex] = -1;
1057 iy[colIndex] = -1;
1058 iySigma[colIndex] = -1;
1059 iFit[colIndex] = -1;
1060 iResidual[colIndex] = -1;
1061 }
1062
1063 if (!SDDS_InitializeOutput(SDDSout, SDDS_BINARY, 0, NULL, "sddssplinefit output: fitted data", output) ||
1064 !SDDS_TransferColumnDefinition(SDDSout, SDDSin, xName, NULL) ||
1065 SDDS_GetColumnInformation(SDDSout, "symbol", &xSymbol, SDDS_GET_BY_NAME, xName) != SDDS_STRING ||
1066 (xSigmaName && !SDDS_TransferColumnDefinition(SDDSout, SDDSin, xSigmaName, NULL)))
1067 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1068
1069 for (colIndex = 0; colIndex < numCols; colIndex++) {
1070 if (!SDDS_TransferColumnDefinition(SDDSout, SDDSin, yNames[colIndex], NULL) ||
1071 SDDS_GetColumnInformation(SDDSout, "symbol", &ySymbols[colIndex], SDDS_GET_BY_NAME, yNames[colIndex]) != SDDS_STRING ||
1072 (ySigmaNames && !SDDS_TransferColumnDefinition(SDDSout, SDDSin, ySigmaNames[colIndex], NULL)))
1073 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1074 }
1075 if (!xSymbol || SDDS_StringIsBlank(xSymbol))
1076 xSymbol = xName;
1077 for (colIndex = 0; colIndex < numCols; colIndex++)
1078 if (!ySymbols[colIndex] || SDDS_StringIsBlank(ySymbols[colIndex]))
1079 ySymbols[colIndex] = yNames[colIndex];
1080 ix = SDDS_GetColumnIndex(SDDSout, xName);
1081 for (colIndex = 0; colIndex < numCols; colIndex++) {
1082 iy[colIndex] = SDDS_GetColumnIndex(SDDSout, yNames[colIndex]);
1083 if (ySigmaNames)
1084 iySigma[colIndex] = SDDS_GetColumnIndex(SDDSout, ySigmaNames[colIndex]);
1085 }
1086 if (xSigmaName)
1087 ixSigma = SDDS_GetColumnIndex(SDDSout, xSigmaName);
1088 if (SDDS_NumberOfErrors())
1089 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1090
1091 for (colIndex = 0; colIndex < numCols; colIndex++) {
1092 sprintf(buffer, "%sFit", yNames[colIndex]);
1093 sprintf(buffer1, "Fit[%s]", ySymbols[colIndex]);
1094 if (!SDDS_TransferColumnDefinition(SDDSout, SDDSin, yNames[colIndex], buffer) ||
1095 !SDDS_ChangeColumnInformation(SDDSout, "symbol", buffer1, SDDS_SET_BY_NAME, buffer))
1096 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1097 if ((iFit[colIndex] = SDDS_GetColumnIndex(SDDSout, buffer)) < 0)
1098 SDDS_Bomb("unable to get index of just-defined fit output column");
1099
1100 sprintf(buffer, "%sResidual", yNames[colIndex]);
1101 sprintf(buffer1, "Residual[%s]", ySymbols[colIndex]);
1102 if (!SDDS_TransferColumnDefinition(SDDSout, SDDSin, yNames[colIndex], buffer) ||
1103 !SDDS_ChangeColumnInformation(SDDSout, "symbol", buffer1, SDDS_SET_BY_NAME, buffer))
1104 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1105 if (!(iResidual[colIndex] = SDDS_GetColumnIndex(SDDSout, buffer)))
1106 SDDS_Bomb("unable to get index of just-defined residual output column");
1107
1108 if (sigmasValid && !ySigmaNames) {
1109 sprintf(buffer, "%sSigma", yNames[colIndex]);
1110 if (!SDDS_TransferColumnDefinition(SDDSout, SDDSin, yNames[colIndex], buffer))
1111 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1112 iySigma[colIndex] = SDDS_GetColumnIndex(SDDSout, buffer);
1113 if (ySymbols[colIndex] && !SDDS_StringIsBlank(ySymbols[colIndex])) {
1114 sprintf(buffer1, "Sigma[%s]", ySymbols[colIndex]);
1115 if (!SDDS_ChangeColumnInformation(SDDSout, "symbol", buffer1, SDDS_SET_BY_NAME, buffer))
1116 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1117 }
1118 }
1119 }
1120
1121 if (outputInfo && !SDDS_InitializeOutput(SDDSoutInfo, SDDS_BINARY, 0, NULL, "sddsspline output: fit information", outputInfo))
1122 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1123
1124 if (outputInfo) {
1125
1126 if (SDDS_DefineParameter(SDDSoutInfo, "Order", NULL, NULL, "Order of term in fit", NULL, SDDS_LONG, 0) < 0)
1127 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1128
1129 if (SDDS_DefineParameter(SDDSoutInfo, "Coefficients", NULL, NULL, "Number of Coefficients", NULL, SDDS_LONG, 0) < 0)
1130 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1131
1132 if (SDDS_DefineParameter(SDDSoutInfo, "Breakpoints", NULL, NULL, "Number of breakpoints", NULL, SDDS_LONG, 0) < 0)
1133 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1134
1135 if (SDDS_GetColumnInformation(SDDSout, "units", &xUnits, SDDS_GET_BY_NAME, xName) != SDDS_STRING)
1136 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1137 sprintf(buffer, "%sOffset", xName);
1138 sprintf(buffer1, "Offset of %s for fit", xName);
1139 if ((iOffset = SDDS_DefineParameter(SDDSoutInfo, buffer, NULL, xUnits, buffer1, NULL, SDDS_DOUBLE, NULL)) < 0)
1140 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1141 sprintf(buffer, "%sScale", xName);
1142 sprintf(buffer1, "Scale factor of %s for fit", xName);
1143 if ((iFactor = SDDS_DefineParameter(SDDSoutInfo, buffer, NULL, xUnits, buffer1, NULL, SDDS_DOUBLE, NULL)) < 0)
1144 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1145
1146 if (SDDS_DefineColumn(SDDSoutInfo, "Index", NULL, NULL, "Index of spline coefficients", NULL, SDDS_LONG, 0) < 0)
1147 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1148 for (colIndex = 0; colIndex < numCols; colIndex++) {
1149 if (SDDS_GetColumnInformation(SDDSout, "units", &yUnits, SDDS_GET_BY_NAME, yNames[colIndex]) != SDDS_STRING)
1150 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1151
1152 sprintf(buffer1, "%sCoefficient", yNames[colIndex]);
1153 sprintf(buffer2, "%sCoefficientSigma", yNames[colIndex]);
1154
1155 if (SDDS_DefineColumn(SDDSoutInfo, buffer1, NULL, yUnits, "Coefficient of spline fit", NULL, SDDS_DOUBLE, 0) < 0 ||
1156 (sigmasValid && SDDS_DefineColumn(SDDSoutInfo, buffer2, "$gs$r$ba$n", "[CoefficientUnits]", "sigma of coefficient of term in fit", NULL, SDDS_DOUBLE, 0) < 0))
1157 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1158
1159 iCoefficient[colIndex] = SDDS_GetColumnIndex(SDDSoutInfo, buffer1); /* this index is used for setting values of that column */
1160 iCoefficientSigma[colIndex] = SDDS_GetColumnIndex(SDDSoutInfo, buffer2);
1161
1162 sprintf(buffer1, "%sReducedChiSquared", yNames[colIndex]);
1163 sprintf(buffer2, "%sRmsResidual", yNames[colIndex]);
1164 sprintf(buffer3, "%sSignificanceLevel", yNames[colIndex]);
1165
1166 if ((iChiSq[colIndex] = SDDS_DefineParameter(SDDSoutInfo, buffer1, "$gh$r$a2$n/(N-M)", NULL,
1167 "Reduced chi-squared of fit",
1168 NULL, SDDS_DOUBLE, NULL)) < 0 ||
1169 SDDS_GetColumnInformation(SDDSout, "units", &yUnits, SDDS_GET_BY_NAME, yNames[colIndex]) != SDDS_STRING ||
1170 (iRmsResidual[colIndex] =
1171 SDDS_DefineParameter(SDDSoutInfo, buffer2, "$gs$r$bres$n", yUnits, "RMS residual of fit", NULL, SDDS_DOUBLE, NULL)) < 0 ||
1172 (iSigLevel[colIndex] = SDDS_DefineParameter(SDDSoutInfo, buffer3, NULL, NULL, "Probability that data is from fit function", NULL, SDDS_DOUBLE, NULL)) < 0)
1173 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1174 if (yUnits)
1175 free(yUnits);
1176
1177 sprintf(buffer, "%sFitIsValid", yNames[colIndex]);
1178 if ((iFitIsValid[colIndex] = SDDS_DefineParameter(SDDSoutInfo, buffer, NULL, NULL, NULL, NULL, SDDS_CHARACTER, NULL)) < 0)
1179 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1180 }
1181 }
1182 if (SDDS_DefineParameter1(SDDSout, "Order", NULL, NULL, "Order of splines", NULL, SDDS_LONG, &order) < 0)
1183 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1184
1185 if (SDDS_DefineParameter1(SDDSout, "Coefficients", NULL, NULL, "Number of coeffs in fit", NULL, SDDS_LONG, &coeffs) < 0)
1186 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1187
1188 if (SDDS_DefineParameter1(SDDSout, "Breakpoints", NULL, NULL, "Number of break points in fit", NULL, SDDS_LONG, &breakpoints) < 0)
1189 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1190
1191 if (SDDS_GetColumnInformation(SDDSout, "units", &xUnits, SDDS_GET_BY_NAME, xName) != SDDS_STRING)
1192 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1193 sprintf(buffer, "%sOffset", xName);
1194 sprintf(buffer1, "Offset of %s for fit", xName);
1195 if ((iOffsetO = SDDS_DefineParameter(SDDSout, buffer, NULL, xUnits, buffer1, NULL, SDDS_DOUBLE, NULL)) < 0)
1196 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1197 sprintf(buffer, "%sScale", xName);
1198 sprintf(buffer1, "Scale factor of %s for fit", xName);
1199 if ((iFactorO = SDDS_DefineParameter(SDDSout, buffer, NULL, xUnits, buffer1, NULL, SDDS_DOUBLE, NULL)) < 0)
1200 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1201
1202 for (colIndex = 0; colIndex < numCols; colIndex++) {
1203
1204 sprintf(buffer1, "%sReducedChiSquared", yNames[colIndex]);
1205 sprintf(buffer2, "%sRmsResidual", yNames[colIndex]);
1206 sprintf(buffer3, "%sSignificanceLevel", yNames[colIndex]);
1207
1208 if ((iChiSqO[colIndex] = SDDS_DefineParameter(SDDSout, buffer1, "$gh$r$a2$n/(N-M)", NULL,
1209 "Reduced chi-squared of fit",
1210 NULL, SDDS_DOUBLE, NULL)) < 0 ||
1211 SDDS_GetColumnInformation(SDDSout, "units", &yUnits, SDDS_GET_BY_NAME, yNames[colIndex]) != SDDS_STRING ||
1212 (iRmsResidualO[colIndex] =
1213 SDDS_DefineParameter(SDDSout, buffer2, "$gs$r$bres$n", yUnits, "RMS residual of fit", NULL, SDDS_DOUBLE, NULL)) < 0 ||
1214 (iSigLevelO[colIndex] = SDDS_DefineParameter(SDDSout, buffer3, NULL, NULL, "Probability that data is from fit function", NULL, SDDS_DOUBLE, NULL)) < 0)
1215 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1216 if (yUnits)
1217 free(yUnits);
1218
1219 sprintf(buffer, "%sFitIsValid", yNames[colIndex]);
1220 if ((iFitIsValidO[colIndex] = SDDS_DefineParameter(SDDSout, buffer, NULL, NULL, NULL, NULL, SDDS_CHARACTER, NULL)) < 0)
1221 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1222 }
1223
1224 if (copyParameters) {
1225 if (!SDDS_TransferAllParameterDefinitions(SDDSout, SDDSin, 0))
1226 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1227 if (outputInfo && !SDDS_TransferAllParameterDefinitions(SDDSoutInfo, SDDSin, 0))
1228 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1229 }
1230
1231 if ((outputInfo && !SDDS_WriteLayout(SDDSoutInfo)) || !SDDS_WriteLayout(SDDSout))
1232 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1233
1234 // return coefUnits;
1235 return;
1236}
1237
1238void setupEvaluationFile(EVAL_PARAMETERS * evalParameters, char *xName, char **yName, long yNames, SDDS_DATASET *SDDSin) {
1239 long iYName, i, iCoeff, coeffs;
1240 char *xSymbol, *ySymbol;
1241 char *mainTemplateFirstDeriv[3] = {"%yNameDeriv", "Derivative w.r.t. %xSymbol of %ySymbol", "d[%ySymbol]/d[%xSymbol]"};
1242 char *mainTemplate[3];
1243 char buffer[1024];
1244 char ***yDerivName, ***yDerivUnits;
1245 long nderiv, derivOrder;
1246 long *iSpline;
1247 gsl_bspline_workspace *bw;
1248 SDDS_DATASET *SDDSout;
1249#if GSL_MAJOR_VERSION == 1
1250 gsl_bspline_deriv_workspace *bdw;
1251#endif
1252 SDDSout = &evalParameters->dataset;
1253 bw = evalParameters->bw;
1254 coeffs = gsl_bspline_ncoeffs(bw);
1255
1256 if (!SDDS_InitializeOutput(SDDSout, SDDS_BINARY, 0, NULL, "sddsspline output: evaluation of spline fits", evalParameters->file) ||
1257 !SDDS_TransferColumnDefinition(SDDSout, SDDSin, xName, NULL))
1258 SDDS_Bomb("Problem setting up evaluation file");
1259 if (SDDS_GetColumnInformation(SDDSout, "symbol", &xSymbol, SDDS_GET_BY_NAME, xName) != SDDS_STRING) {
1260 /* fprintf(stderr, "error: problem getting symbol for column %s\n", xName);*/
1261 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1262 }
1263 if (!xSymbol)
1264 SDDS_CopyString(&xSymbol, xName);
1265
1266 if (evalParameters->flags & EVAL_PROVIDEBASIS) {
1267 iSpline = tmalloc(sizeof(*iSpline) * coeffs); /* dataset index of spline column names */
1268 evalParameters->iSpline = iSpline;
1269 for (iCoeff = 0; iCoeff < coeffs; iCoeff++) {
1270 sprintf(buffer, "B%04ld", iCoeff);
1271 if ((iSpline[iCoeff] = SDDS_DefineColumn(SDDSout, buffer, NULL, NULL, NULL, NULL, SDDS_DOUBLE, 0)) < 0)
1272 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1273 }
1274 }
1275
1276 if (evalParameters->flags & EVAL_DERIVATIVES) {
1277 nderiv = evalParameters->nderiv;
1278 yDerivName = tmalloc(sizeof(*yDerivName) * (nderiv + 1));
1279 evalParameters->yDerivName = yDerivName;
1280 yDerivUnits = tmalloc(sizeof(*yDerivUnits) * (nderiv + 1));
1281 evalParameters->yDerivUnits = yDerivUnits;
1282 for (derivOrder = 0; derivOrder <= nderiv; derivOrder++) {
1283 yDerivName[derivOrder] = tmalloc(sizeof(**yDerivName) * yNames);
1284 yDerivUnits[derivOrder] = tmalloc(sizeof(**yDerivUnits) * yNames);
1285 }
1286 for (iYName = 0; iYName < yNames; iYName++) {
1287 if (!SDDS_TransferColumnDefinition(SDDSout, SDDSin, yName[iYName], NULL))
1288 SDDS_Bomb("Problem setting up evaluation file");
1289 yDerivName[0][iYName] = yName[iYName];
1290 /* first derivative is the first one */
1291 for (derivOrder = 1; derivOrder <= nderiv; derivOrder++) {
1292 for (i = 0; i < 3; i++) {
1293 if (derivOrder != 1) {
1294 switch (i) {
1295 case 0:
1296 /* name */
1297 sprintf(buffer, "%%yNameDeriv%ld", derivOrder);
1298 break;
1299 case 1:
1300 /* description */
1301 sprintf(buffer, "Derivative %ld w.r.t. %%xSymbol of %%ySymbol", derivOrder);
1302 break;
1303 case 2:
1304 /* symbol */
1305 sprintf(buffer, "d$a%ld$n[%%ySymbol]/d[%%xSymbol]$a%ld$n", derivOrder, derivOrder);
1306 break;
1307 }
1308 cp_str(&mainTemplate[i], buffer);
1309 } else {
1310 mainTemplate[i] = mainTemplateFirstDeriv[i];
1311 }
1312 }
1313 if (SDDS_GetColumnInformation(SDDSout, "symbol", &ySymbol, SDDS_GET_BY_NAME, yName[iYName]) != SDDS_STRING) {
1314 fprintf(stderr, "error: problem getting symbol for column %s\n", yName[iYName]);
1315 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1316 }
1317 if (!ySymbol || SDDS_StringIsBlank(ySymbol))
1318 SDDS_CopyString(&ySymbol, yName[iYName]);
1319 /* I'm using the function changeInformation from sddsderiv which requires an existing column of some kind*/
1320 if (SDDS_DefineColumn(SDDSout, "placeholderName", NULL, NULL, NULL, NULL, SDDS_DOUBLE, 0) < 0)
1321 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1322 yDerivUnits[derivOrder][iYName] = (char *)divideColumnUnits(SDDSout, yDerivName[derivOrder - 1][iYName], xName);
1323 yDerivName[derivOrder][iYName] = (char *)changeInformation(SDDSout, "placeholderName", yDerivName[0][iYName],
1324 ySymbol, xName, xSymbol, mainTemplate, yDerivUnits[derivOrder][iYName]);
1325 }
1326 }
1327 if (!SDDS_WriteLayout(SDDSout))
1328 SDDS_Bomb("Problem setting up evaluation file with derivatives");
1329 } else {
1330 for (iYName = 0; iYName < yNames; iYName++) {
1331 if (!SDDS_TransferColumnDefinition(SDDSout, SDDSin, yName[iYName], NULL))
1332 SDDS_Bomb("Problem setting up evaluation file");
1333 }
1334 if (!SDDS_WriteLayout(SDDSout))
1335 SDDS_Bomb("Problem setting up evaluation file");
1336 }
1337}
1338
1339/* called for each column, for each page */
1340void makeEvaluationTable(EVAL_PARAMETERS * evalParameters, double *x, int64_t points,
1341 gsl_vector *B, gsl_matrix *cov, gsl_vector *c,
1342 char *xName, char **yName, long yNames, long iYName, long int order) {
1343 static double *xEval = NULL, *yEval = NULL;
1344 static int64_t maxEvals = 0;
1345 double delta;
1346 double yerr;
1347 int64_t i;
1348 char ***yDerivName;
1349 double xi;
1350 gsl_matrix *dB = NULL;
1351 long nderiv, coeffs, iCoeff, derivOrder;
1352 double **yDeriv = NULL;
1353 double **Bspline;
1354 long *iSpline;
1355 gsl_bspline_workspace *bw;
1356 SDDS_DATASET *SDDSout;
1357#if GSL_MAJOR_VERSION == 1
1358 gsl_bspline_deriv_workspace *bdw;
1359#endif
1360 SDDSout = &evalParameters->dataset;
1361 yDerivName = evalParameters->yDerivName;
1362 iSpline = evalParameters->iSpline;
1363 bw = evalParameters->bw;
1364 coeffs = gsl_bspline_ncoeffs(bw);
1365
1366 if (!(evalParameters->flags & EVAL_BEGIN_GIVEN) || !(evalParameters->flags & EVAL_END_GIVEN)) {
1367 double min, max;
1368 find_min_max(&min, &max, x, points);
1369 if (!(evalParameters->flags & EVAL_BEGIN_GIVEN))
1370 evalParameters->begin = min;
1371 if (!(evalParameters->flags & EVAL_END_GIVEN))
1372 evalParameters->end = max;
1373 }
1374 if (!(evalParameters->flags & EVAL_NUMBER_GIVEN))
1375 evalParameters->number = points;
1376 if (evalParameters->number > 1)
1377 delta = (evalParameters->end - evalParameters->begin) / (evalParameters->number - 1);
1378 else
1379 delta = 0;
1380
1381 if (!xEval || maxEvals < evalParameters->number) {
1382 if (!(xEval = (double *)SDDS_Realloc(xEval, sizeof(*xEval) * evalParameters->number)) ||
1383 !(yEval = (double *)SDDS_Realloc(yEval, sizeof(*yEval) * evalParameters->number)))
1384 SDDS_Bomb("allocation failure");
1385 maxEvals = evalParameters->number;
1386 }
1387
1388 Bspline = tmalloc(sizeof(*Bspline) * coeffs);
1389 /* the allocation and calculation of Bspline is done only on the first column */
1390 if ((iYName == 0) && (evalParameters->flags & EVAL_PROVIDEBASIS)) {
1391 for (iCoeff = 0; iCoeff < coeffs; iCoeff++) {
1392 Bspline[iCoeff] = tmalloc(sizeof(**Bspline) * evalParameters->number);
1393 }
1394 }
1395
1396 if (!(evalParameters->flags & EVAL_DERIVATIVES)) {
1397 for (i = 0; i < evalParameters->number; i++) {
1398 xi = evalParameters->begin + i * delta;
1399 xEval[i] = xi;
1400 gsl_bspline_eval(xi, B, bw);
1401 gsl_multifit_linear_est(B, c, cov, &yEval[i], &yerr);
1402 if ((iYName == 0) && (evalParameters->flags & EVAL_PROVIDEBASIS)) {
1403 for (iCoeff = 0; iCoeff < coeffs; iCoeff++) {
1404 Bspline[iCoeff][i] = gsl_vector_get(B, iCoeff);
1405 }
1406 }
1407 }
1408 /* Oh, this filters so that StartPage and WritePage are run only once. Tricksty */
1409 /* On the other hand the x values is repeatedly assigned the values until the calls stop */
1410 if ((iYName == 0 &&
1411 !SDDS_StartPage(SDDSout, evalParameters->number)) ||
1412 !SDDS_SetColumnFromDoubles(SDDSout, SDDS_SET_BY_NAME, xEval, evalParameters->number, xName) ||
1413 !SDDS_SetColumnFromDoubles(SDDSout, SDDS_SET_BY_NAME, yEval, evalParameters->number, yName[iYName]))
1414 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1415 if ((iYName == 0) && (evalParameters->flags & EVAL_PROVIDEBASIS)) {
1416 for (iCoeff = 0; iCoeff < coeffs; iCoeff++)
1417 if (!SDDS_SetColumnFromDoubles(SDDSout, SDDS_SET_BY_INDEX, Bspline[iCoeff], evalParameters->number, iSpline[iCoeff]))
1418 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1419 }
1420 if ((iYName == yNames - 1 && !SDDS_WritePage(SDDSout)))
1421 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1422 } else {
1423 nderiv = evalParameters->nderiv;
1424 yDeriv = tmalloc(sizeof(double *) * (nderiv + 1));
1425 for (derivOrder = 0; derivOrder <= nderiv; derivOrder++) {
1426 yDeriv[derivOrder] = tmalloc(sizeof(double) * evalParameters->number);
1427 }
1428 dB = gsl_matrix_alloc(coeffs, nderiv + 1);
1429 for (i = 0; i < evalParameters->number; i++) {
1430 xi = evalParameters->begin + i * delta;
1431 xEval[i] = xi;
1432 gsl_bspline_eval(xi, B, bw);
1433 gsl_multifit_linear_est(B, c, cov, &yEval[i], &yerr);
1434#if GSL_MAJOR_VERSION == 1
1435 bdw = gsl_bspline_deriv_alloc(order);
1436 gsl_bspline_deriv_eval(xi, nderiv, dB, bw, bdw);
1437 gsl_bspline_deriv_free(bdw);
1438#else
1439 gsl_bspline_deriv_eval(xi, nderiv, dB, bw);
1440#endif
1441 for (derivOrder = 0; derivOrder <= nderiv; derivOrder++) {
1442 yDeriv[derivOrder][i] = 0;
1443 for (iCoeff = 0; iCoeff < coeffs; iCoeff++) {
1444 yDeriv[derivOrder][i] += gsl_vector_get(c, iCoeff) * gsl_matrix_get(dB, iCoeff, derivOrder);
1445 }
1446 }
1447 if ((iYName == 0) && (evalParameters->flags & EVAL_PROVIDEBASIS)) {
1448 for (iCoeff = 0; iCoeff < coeffs; iCoeff++) {
1449 Bspline[iCoeff][i] = gsl_vector_get(B, iCoeff);
1450 }
1451 }
1452 }
1453 gsl_matrix_free(dB);
1454 /* Oh, this filters so that StartPage and WritePage are run only once. Tricksty */
1455 /* On the other hand the x values is repeatedly assigned the values until the calls stop */
1456 if ((iYName == 0 &&
1457 !SDDS_StartPage(SDDSout, evalParameters->number)) ||
1458 !SDDS_SetColumnFromDoubles(SDDSout, SDDS_SET_BY_NAME, xEval, evalParameters->number, xName) ||
1459 !SDDS_SetColumnFromDoubles(SDDSout, SDDS_SET_BY_NAME, yEval, evalParameters->number, yName[iYName]))
1460 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1461 /* set derivatives */ /* don't need the 0th derivative */
1462 for (derivOrder = 1; derivOrder <= nderiv; derivOrder++) {
1463 if (!SDDS_SetColumnFromDoubles(SDDSout, SDDS_SET_BY_NAME, yDeriv[derivOrder], evalParameters->number, yDerivName[derivOrder][iYName]))
1464 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1465 }
1466 if ((iYName == 0) && (evalParameters->flags & EVAL_PROVIDEBASIS)) {
1467 for (iCoeff = 0; iCoeff < coeffs; iCoeff++)
1468 if (!SDDS_SetColumnFromDoubles(SDDSout, SDDS_SET_BY_INDEX, Bspline[iCoeff], evalParameters->number, iSpline[iCoeff]))
1469 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1470 }
1471 if ((iYName == yNames - 1 && !SDDS_WritePage(SDDSout)))
1472 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1473 if (yDeriv) {
1474 for (derivOrder = 0; derivOrder <= nderiv; derivOrder++) {
1475 free(yDeriv[derivOrder]);
1476 }
1477 free(yDeriv);
1478 }
1479 }
1480 if ((iYName == 0) && (evalParameters->flags & EVAL_PROVIDEBASIS)) {
1481 for (iCoeff = 0; iCoeff < coeffs; iCoeff++) {
1482 free(Bspline[iCoeff]);
1483 }
1484 }
1485 if (iYName == 0)
1486 free(Bspline);
1487}
1488
1489char *changeInformation(SDDS_DATASET * SDDSout, char *name, char *nameRoot, char *symbolRoot, char *xName, char *xSymbol, char **template, char *newUnits) {
1490 char buffer1[SDDS_MAXLINE], buffer2[SDDS_MAXLINE], *ptr;
1491
1492 if (!SDDS_ChangeColumnInformation(SDDSout, "units", newUnits, SDDS_PASS_BY_VALUE | SDDS_SET_BY_NAME, name))
1493 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1494
1495 makeSubstitutions(buffer1, buffer2, template[2], nameRoot, symbolRoot, xName, xSymbol);
1496 if (!SDDS_ChangeColumnInformation(SDDSout, "symbol", buffer2, SDDS_PASS_BY_VALUE | SDDS_SET_BY_NAME, name))
1497 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1498
1499 makeSubstitutions(buffer1, buffer2, template[1], nameRoot, symbolRoot, xName, xSymbol);
1500 if (!SDDS_ChangeColumnInformation(SDDSout, "description", buffer2, SDDS_PASS_BY_VALUE | SDDS_SET_BY_NAME, name))
1501 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1502
1503 makeSubstitutions(buffer1, buffer2, template[0], nameRoot, symbolRoot, xName, xSymbol);
1504 if (!SDDS_ChangeColumnInformation(SDDSout, "name", buffer2, SDDS_PASS_BY_VALUE | SDDS_SET_BY_NAME, name))
1505 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1506 SDDS_CopyString(&ptr, buffer2);
1507 return ptr;
1508}
1509
1510void makeSubstitutions(char *buffer1, char *buffer2, char *template, char *nameRoot, char *symbolRoot, char *xName, char *xSymbol) {
1511 strcpy(buffer2, template);
1512 replace_string(buffer1, buffer2, "%ySymbol", symbolRoot);
1513 replace_string(buffer2, buffer1, "%xSymbol", xSymbol);
1514 replace_string(buffer1, buffer2, "%yName", nameRoot);
1515 replace_string(buffer2, buffer1, "%xName", xName);
1516 strcpy(buffer1, buffer2);
1517}
1518
1519int print_matrix(FILE * f, const gsl_matrix *m) {
1520 int status, n = 0;
1521 size_t i, j;
1522 for (i = 0; i < m->size1; i++) {
1523 for (j = 0; j < m->size2; j++) {
1524 if ((status = fprintf(f, "%10.6lf ", gsl_matrix_get(m, i, j))) < 0)
1525 return -1;
1526 n += status;
1527 }
1528
1529 if ((status = fprintf(f, "\n")) < 0)
1530 return -1;
1531 n += status;
1532 }
1533
1534 return n;
1535}
SDDS.h
SDDS (Self Describing Data Set) Data Types Definitions and Function Prototypes.
SDDS_CopyParameters
int32_t SDDS_CopyParameters(SDDS_DATASET *SDDS_target, SDDS_DATASET *SDDS_source)
Definition SDDS_copy.c:286
SDDS_StartPage
int32_t SDDS_StartPage(SDDS_DATASET *SDDS_dataset, int64_t expected_n_rows)
Definition SDDS_dataprep.c:64
SDDS_SetParameters
int32_t SDDS_SetParameters(SDDS_DATASET *SDDS_dataset, int32_t mode,...)
Definition SDDS_dataprep.c:369
SDDS_SetColumnFromDoubles
int32_t SDDS_SetColumnFromDoubles(SDDS_DATASET *SDDS_dataset, int32_t mode, double *data, int64_t rows,...)
Sets the values for a single data column using double-precision floating-point numbers.
Definition SDDS_dataprep.c:1642
SDDS_SetColumn
int32_t SDDS_SetColumn(SDDS_DATASET *SDDS_dataset, int32_t mode, void *data, int64_t rows,...)
Sets the values for one data column in the current data table of an SDDS dataset.
Definition SDDS_dataprep.c:1536
SDDS_CountRowsOfInterest
int64_t SDDS_CountRowsOfInterest(SDDS_DATASET *SDDS_dataset)
Counts the number of rows marked as "of interest" in the current data table.
Definition SDDS_extract.c:364
SDDS_SetColumnsOfInterest
int32_t SDDS_SetColumnsOfInterest(SDDS_DATASET *SDDS_dataset, int32_t mode,...)
Sets the acceptance flags for columns based on specified naming criteria.
Definition SDDS_extract.c:432
SDDS_GetColumnInDoubles
double * SDDS_GetColumnInDoubles(SDDS_DATASET *SDDS_dataset, char *column_name)
Retrieves the data of a specified numerical column as an array of doubles, considering only rows mark...
Definition SDDS_extract.c:872
SDDS_ChangeColumnInformation
int32_t SDDS_ChangeColumnInformation(SDDS_DATASET *SDDS_dataset, char *field_name, void *memory, int32_t mode,...)
Modifies a specific field in a column definition within the SDDS dataset.
Definition SDDS_info.c:364
SDDS_GetColumnInformation
int32_t SDDS_GetColumnInformation(SDDS_DATASET *SDDS_dataset, char *field_name, void *memory, int32_t mode,...)
Retrieves information about a specified column in the SDDS dataset.
Definition SDDS_info.c:41
SDDS_InitializeInput
int32_t SDDS_InitializeInput(SDDS_DATASET *SDDS_dataset, char *filename)
Definition SDDS_input.c:49
SDDS_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_DefineParameter1
int32_t SDDS_DefineParameter1(SDDS_DATASET *SDDS_dataset, const char *name, const char *symbol, const char *units, const char *description, const char *format_string, int32_t type, void *fixed_value)
Defines a data parameter with a fixed numerical value.
Definition SDDS_output.c:1407
SDDS_InitializeOutput
int32_t SDDS_InitializeOutput(SDDS_DATASET *SDDS_dataset, int32_t data_mode, int32_t lines_per_row, const char *description, const char *contents, const char *filename)
Initializes the SDDS output dataset.
Definition SDDS_output.c:601
SDDS_WritePage
int32_t SDDS_WritePage(SDDS_DATASET *SDDS_dataset)
Writes the current data table to the output file.
Definition SDDS_output.c:1228
SDDS_DefineColumn
int32_t SDDS_DefineColumn(SDDS_DATASET *SDDS_dataset, const char *name, const char *symbol, const char *units, const char *description, const char *format_string, int32_t type, int32_t field_length)
Defines a data column within the SDDS dataset.
Definition SDDS_output.c:1715
SDDS_WriteLayout
int32_t SDDS_WriteLayout(SDDS_DATASET *SDDS_dataset)
Writes the SDDS layout header to the output file.
Definition SDDS_output.c:899
SDDS_DefineParameter
int32_t SDDS_DefineParameter(SDDS_DATASET *SDDS_dataset, const char *name, const char *symbol, const char *units, const char *description, const char *format_string, int32_t type, char *fixed_value)
Defines a data parameter with a fixed string value.
Definition SDDS_output.c:1472
SDDS_TransferColumnDefinition
int32_t SDDS_TransferColumnDefinition(SDDS_DATASET *target, SDDS_DATASET *source, char *name, char *newName)
Transfers a column definition from a source dataset to a target dataset.
Definition SDDS_transfer.c:35
SDDS_TransferAllParameterDefinitions
int32_t SDDS_TransferAllParameterDefinitions(SDDS_DATASET *SDDS_target, SDDS_DATASET *SDDS_source, uint32_t mode)
Transfers all parameter definitions from a source dataset to a target dataset.
Definition SDDS_transfer.c:343
SDDS_GetColumnIndex
int32_t SDDS_GetColumnIndex(SDDS_DATASET *SDDS_dataset, char *name)
Retrieves the index of a named column in the SDDS dataset.
Definition SDDS_utils.c:1309
SDDS_CheckColumn
int32_t SDDS_CheckColumn(SDDS_DATASET *SDDS_dataset, char *name, char *units, int32_t type, FILE *fp_message)
Checks if a column exists in the SDDS dataset with the specified name, units, and type.
Definition SDDS_utils.c:4562
SDDS_GetColumnNames
char ** SDDS_GetColumnNames(SDDS_DATASET *SDDS_dataset, int32_t *number)
Retrieves the names of all columns in the SDDS dataset.
Definition SDDS_utils.c:2460
SDDS_PrintErrors
void SDDS_PrintErrors(FILE *fp, int32_t mode)
Prints recorded error messages to a specified file stream.
Definition SDDS_utils.c:432
SDDS_RegisterProgramName
void SDDS_RegisterProgramName(const char *name)
Registers the executable program name for use in error messages.
Definition SDDS_utils.c:288
SDDS_NumberOfErrors
int32_t SDDS_NumberOfErrors()
Retrieves the number of errors recorded by SDDS library routines.
Definition SDDS_utils.c:304
SDDS_StringIsBlank
int32_t SDDS_StringIsBlank(char *s)
Checks if a string is blank (contains only whitespace characters).
Definition SDDS_utils.c:2404
SDDS_FindColumn
char * SDDS_FindColumn(SDDS_DATASET *SDDS_dataset, int32_t mode,...)
Finds the first column in the SDDS dataset that matches the specified criteria.
Definition SDDS_utils.c:4312
SDDS_Bomb
void SDDS_Bomb(char *message)
Terminates the program after printing an error message and recorded errors.
Definition SDDS_utils.c:342
SDDS_CopyString
int32_t SDDS_CopyString(char **target, const char *source)
Copies a source string to a target string with memory allocation.
Definition SDDS_utils.c:856
SDDS_Realloc
void * SDDS_Realloc(void *old_ptr, size_t new_size)
Reallocates memory to a new size.
Definition SDDS_utils.c:677
SDDS_STRING
#define SDDS_STRING
Identifier for the string data type.
Definition SDDStypes.h:85
SDDS_LONG
#define SDDS_LONG
Identifier for the signed 32-bit integer data type.
Definition SDDStypes.h:61
SDDS_CHARACTER
#define SDDS_CHARACTER
Identifier for the character data type.
Definition SDDStypes.h:91
SDDS_ANY_NUMERIC_TYPE
#define SDDS_ANY_NUMERIC_TYPE
Special identifier used by SDDS_Check*() routines to accept any numeric type.
Definition SDDStypes.h:157
SDDS_DOUBLE
#define SDDS_DOUBLE
Identifier for the double data type.
Definition SDDStypes.h:37
SDDS_LONG64
#define SDDS_LONG64
Identifier for the signed 64-bit integer data type.
Definition SDDStypes.h:49
SDDSutils.h
Utility functions for SDDS dataset manipulation and string array operations.
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
cp_str
char * cp_str(char **s, char *t)
Copies a string, allocating memory for storage.
Definition cp_str.c:28
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
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
replace_string
int replace_string(char *t, char *s, char *orig, char *repl)
Replace all occurrences of one string with another string.
Definition replace_string.c:29
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
free_scanargs
void free_scanargs(SCANNED_ARG **scanned, int argc)
Definition scanargs.c:584
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
ipower
double ipower(double x, long n)
Evaluate a power function x^n.
Definition lsfBasisFns.c:113
ChiSqrSigLevel
double ChiSqrSigLevel(double ChiSquared0, long nu)
Computes the probability that a chi-squared variable exceeds a given value.
Definition sigLevel.c:64
str_tolower
char * str_tolower(char *s)
Convert a string to lower case.
Definition str_tolower.c:27
EVAL_PARAMETERS
Definition sddsmpfit.c:286
SDDS_DATASET
Definition SDDS.h:308
Generated by doxygen 1.12.0