SDDS ToolKit Programs and Libraries for C and Python
•All Classes Files Functions Variables Macros Pages
Functions
sddsmpfit.c File Reference

Detailed Description

Performs polynomial least-squares fitting on SDDS files.

This program reads SDDS (Self Describing Data Sets) files, fits the data using either ordinary or Chebyshev T polynomials, and outputs the results including fitted data, residuals, and fit parameters with their uncertainties. Additional features include normalization, sigma modification, and fit evaluation.

The fitting model is:

\[ y = \sum_{i=0}^{N-1} A[i] \cdot P(x - x_{offset}, i) \]

where ( P(x, i) ) is the ith basis function, typically ( x^i ), or a Chebyshev polynomial.

Usage

sddsmpfit [<inputfile>] [<outputfile>]
[-pipe=[input][,output]]
-independent=<xName>
-dependent=<yname1-wildcard>[,<yname2-wildcard>...]
[-sigmaIndependent=<xSigma>]
[-sigmaDependent=<ySigmaFormatString>]
{
-terms=<number> [-symmetry={none|odd|even}] |
-orders=<number>[,<number>...]
}
[-reviseOrders[=threshold=<value>][,verbose]]
[-chebyshev[=convert]]
[-xOffset=<value>]
[-xFactor=<value>]
[-sigmas=<value>,{absolute|fractional}]
[-minimumSigma=<value>]
[-modifySigmas]
[-generateSigmas={keepLargest|keepSmallest}]
[-sparse=<interval>]
[-range=<lower>,<upper>[,fitOnly]]
[-normalize[=<termNumber>]]
[-verbose]
[-evaluate=<filename>[,begin=<value>][,end=<value>][,number=<integer>]]
[-fitLabelFormat=<sprintf-string>]
[-infoFile=<filename>]
[-copyParameters]

Options

Required Description
-independent Specify the independent data column name.
-dependent Specify dependent data columns, using wildcards if needed.
Optional Description
-pipe Use standard input and/or output.
-terms Number of terms to include in the fit.
-symmetry Symmetry of the fit about x_offset.
-orders Specify the polynomial orders to include in the fit.
-reviseOrders Revise the orders to eliminate poorly-determined coefficients.
-chebyshev Use Chebyshev T polynomials. Specify convert to convert back to ordinary polynomials.
-xOffset Set the x-offset for the fit.
-xFactor Set a scaling factor for x values.
-sigmas Set sigma values as absolute or fractional.
-minimumSigma Set a minimum sigma value; smaller values are replaced.
-modifySigmas Modify y-sigmas using x-sigmas and the initial fit.
-generateSigmas Generate y-sigmas based on RMS deviation, optionally keeping the largest/smallest sigmas.
-sparse Sample data at specified intervals.
-range Fit and evaluate within the specified range.
-normalize Normalize coefficients so that a specific term equals 1.
-evaluate Evaluate the fit over a specified range.
-fitLabelFormat Specify a format string for fit labels.
-infoFile Output fit coefficients and statistics to the specified information file.
-copyParameters Copy parameters from the input to the output SDDS file.
-verbose Enable verbose output for debugging or additional information.

Incompatibilities

  • -terms is incompatible with -orders.
  • -generateSigmas is incompatible with -modifySigmas.
  • -sigmas cannot be used with -sigmaDependent.
  • -reviseOrders requires -generateSigmas, -sigmas, or -sigmaDependent.
Copyright
  • (c) 2002 The University of Chicago, as Operator of Argonne National Laboratory.
  • (c) 2002 The Regents of the University of California, as Operator of Los Alamos National Laboratory.
License
This file is distributed under the terms of the Software License Agreement found in the file LICENSE included with this distribution.
Author
M. Borland, Brad Dolin, R. Soliday, H. Shang

Definition in file sddsmpfit.c.

#include "mdb.h"
#include "SDDS.h"
#include "scan.h"

Go to the source code of this file.

Functions

void print_coefs (FILE *fprec, double x_offset, double x_scale, long chebyshev, double *coef, double *s_coef, int32_t *order, long n_terms, double chi, long norm_term, char *prepend)
 
char ** makeCoefficientUnits (SDDS_DATASET *SDDSout, char *xName, char *yName, int32_t *order, long terms)
 
long setCoefficientData (SDDS_DATASET *SDDSout, double *coef, double *coefSigma, char **coefUnits, long *order, long terms)
 
char *** initializeOutputFile (SDDS_DATASET *SDDSout, SDDS_DATASET *SDDSoutInfo, char *output, char *outputInfo, SDDS_DATASET *SDDSin, char *input, char *xName, char **yNames, char *xSigmaName, char **ySigmaNames, long sigmasValid, int32_t *order, long terms, long isChebyshev, long numCols, long copyParameters)
 
void checkInputFile (SDDS_DATASET *SDDSin, char *xName, char **yNames, char *xSigmaName, char **ySigmaNames, long numYNames)
 
long coefficient_index (int32_t *order, long terms, long order_of_interest)
 
void makeFitLabel (char *buffer, long bufsize, char *fitLabelFormat, double *coef, int32_t *order, long terms, long colIndex)
 
char ** ResolveColumnNames (SDDS_DATASET *SDDSin, char **wildcardList, long length, int32_t *numYNames)
 
char ** GenerateYSigmaNames (char *controlString, char **yNames, long numYNames)
 
void RemoveElementFromStringArray (char **array, long index, long length)
 
char ** RemoveNonNumericColumnsFromNameArray (SDDS_DATASET *SDDSin, char **columns, int32_t *numColumns)
 
void compareOriginalToFit (double *x, double *y, double **residual, int64_t points, double *rmsResidual, double *coef, int32_t *order, long terms)
 
void set_argument_offset (double offset)
 Set the offset applied to the input argument of basis functions.
 
void set_argument_scale (double scale)
 Set the scale factor applied to the input argument of basis functions.
 
double tcheby (double x, long n)
 Evaluate the Chebyshev polynomial of the first kind T_n(x).
 
double dtcheby (double x, long n)
 Evaluate the derivative of the Chebyshev polynomial T_n(x).
 
double ipower (double x, long n)
 Evaluate a power function x^n.
 
double dipower (double x, long n)
 Evaluate the derivative of x^n.
 
void setupEvaluationFile (EVAL_PARAMETERS *evalParameters, char *xName, char **yName, long yNames, SDDS_DATASET *SDDSin)
 
void makeEvaluationTable (EVAL_PARAMETERS *evalParameters, double *x, int64_t points, double *coef, int32_t *order, long terms, char *xName, char **yName, long yNames, long iYName)
 
int main (int argc, char **argv)
 
double rms_average (double *x, int64_t n)
 

Function Documentation

◆ checkInputFile()

void checkInputFile ( SDDS_DATASET * SDDSin,
char * xName,
char ** yNames,
char * xSigmaName,
char ** ySigmaNames,
long numYNames )

Definition at line 1190 of file sddsmpfit.c.

1190 {
1191 char *ptr = NULL;
1192 long i;
1193
1194 if (!(ptr = SDDS_FindColumn(SDDSin, FIND_NUMERIC_TYPE, xName, NULL)))
1195 SDDS_Bomb("x column doesn't exist or is nonnumeric");
1196 free(ptr);
1197
1198 /* y columns don't need to be checked because located using SDDS_SetColumnsOfInterest */
1199
1200 ptr = NULL;
1201 if (xSigmaName && !(ptr = SDDS_FindColumn(SDDSin, FIND_NUMERIC_TYPE, xSigmaName, NULL)))
1202 SDDS_Bomb("x sigma column doesn't exist or is nonnumeric");
1203 if (ptr)
1204 free(ptr);
1205
1206 if (ySigmaNames) {
1207 for (i = 0; i < numYNames; i++) {
1208 ptr = NULL;
1209 if (!(ptr = SDDS_FindColumn(SDDSin, FIND_NUMERIC_TYPE, ySigmaNames[i], NULL)))
1210 SDDS_Bomb("y sigma column doesn't exist or is nonnumeric");
1211 if (ptr)
1212 free(ptr);
1213 }
1214 }
1215}
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

◆ coefficient_index()

long coefficient_index ( int32_t * order,
long terms,
long order_of_interest )

Definition at line 1182 of file sddsmpfit.c.

1182 {
1183 long i;
1184 for (i = 0; i < terms; i++)
1185 if (order[i] == order_of_interest)
1186 return (i);
1187 return (-1);
1188}

◆ compareOriginalToFit()

void compareOriginalToFit ( double * x,
double * y,
double ** residual,
int64_t points,
double * rmsResidual,
double * coef,
int32_t * order,
long terms )

Definition at line 1619 of file sddsmpfit.c.

1619 {
1620 int64_t i;
1621 double residualSum2, fit;
1622
1623 *residual = tmalloc(sizeof(**residual) * points);
1624
1625 residualSum2 = 0;
1626 for (i = 0; i < points; i++) {
1627 fit = eval_sum(basis_fn, coef, order, terms, x[i]);
1628 (*residual)[i] = y[i] - fit;
1629 residualSum2 += sqr((*residual)[i]);
1630 }
1631 *rmsResidual = sqrt(residualSum2 / points);
1632}
tmalloc
void * tmalloc(uint64_t size_of_block)
Allocates a memory block of the specified size with zero initialization.
Definition array.c:59
eval_sum
double eval_sum(double(*fn)(double x, long ord), double *coef, int32_t *order, long n_coefs, double x0)
Evaluate a sum of basis functions.
Definition lsfBasisFns.c:145

◆ dipower()

double dipower ( double x,
long n )

Evaluate the derivative of x^n.

This function returns d/dx [ (x - x_offset)/x_scale ]^n = n * ((x - x_offset)/x_scale)^(n-1) / x_scale.

Parameters
xThe point at which to evaluate the derivative.
nThe exponent.
Returns
The derivative of the power function at the given x.

Definition at line 127 of file lsfBasisFns.c.

127 {
128 x = (x - x_offset) / x_scale;
129 return (n * ipow(x, n - 1));
130}
ipow
double ipow(const double x, const int64_t p)
Compute x raised to the power p (x^p).
Definition ipow.c:33

◆ dtcheby()

double dtcheby ( double x,
long n )

Evaluate the derivative of the Chebyshev polynomial T_n(x).

This function returns d/dx [T_n((x - x_offset)/x_scale)]. If x is out of the domain [-1,1], it is clipped to ±1 before evaluation.

Parameters
xThe point at which to evaluate the derivative of T_n.
nThe order of the Chebyshev polynomial.
Returns
The derivative dT_n/dx at the given x.

Definition at line 92 of file lsfBasisFns.c.

92 {
93 x = (x - x_offset) / x_scale;
94 if (x > 1 || x < -1) {
95 /* fprintf(stderr, "warning: argument %e is out of range for tcheby()\n",
96 * x); */
97 x = SIGN(x);
98 }
99 if (x != 1 && x != -1)
100 return (n * sin(n * acos(x)) / sqrt(1 - sqr(x)));
101 return (1.0 * n * n);
102}

◆ GenerateYSigmaNames()

char ** GenerateYSigmaNames ( char * controlString,
char ** yNames,
long numYNames )

Definition at line 1129 of file sddsmpfit.c.

1129 {
1130 long i, nameLength;
1131 char **result, sigmaName[MAX_Y_SIGMA_NAME_SIZE];
1132
1133 result = tmalloc(sizeof(char *) * numYNames);
1134 for (i = 0; i < numYNames; i++) {
1135 sprintf(sigmaName, controlString, yNames[i]);
1136 nameLength = strlen(sigmaName);
1137 result[i] = tmalloc(sizeof(char) * (nameLength + 1));
1138 strcpy(result[i], sigmaName);
1139 }
1140 return (result);
1141}

◆ initializeOutputFile()

char *** initializeOutputFile ( SDDS_DATASET * SDDSout,
SDDS_DATASET * SDDSoutInfo,
char * output,
char * outputInfo,
SDDS_DATASET * SDDSin,
char * input,
char * xName,
char ** yNames,
char * xSigmaName,
char ** ySigmaNames,
long sigmasValid,
int32_t * order,
long terms,
long isChebyshev,
long numCols,
long copyParameters )

Definition at line 1217 of file sddsmpfit.c.

1220 {
1221 char buffer[SDDS_MAXLINE], buffer1[SDDS_MAXLINE], buffer2[SDDS_MAXLINE], buffer3[SDDS_MAXLINE];
1222 char *xUnits, *yUnits, ***coefUnits;
1223 long i, colIndex;
1224
1225 /* all array names followed by an 'O' contain the index of the parameter in the main output file; others refer to
1226 parameters in the infoFile */
1227 coefUnits = tmalloc(sizeof(char **) * numCols);
1228 ySymbols = tmalloc(sizeof(char *) * numCols);
1229 iIntercept = tmalloc(sizeof(long) * numCols);
1230 iInterceptO = tmalloc(sizeof(long) * numCols);
1231 iInterceptSigma = tmalloc(sizeof(long) * numCols);
1232 iInterceptSigmaO = tmalloc(sizeof(long) * numCols);
1233 iSlope = tmalloc(sizeof(long) * numCols);
1234 iSlopeO = tmalloc(sizeof(long) * numCols);
1235 iSlopeSigma = tmalloc(sizeof(long) * numCols);
1236 iSlopeSigmaO = tmalloc(sizeof(long) * numCols);
1237 iCurvature = tmalloc(sizeof(long) * numCols);
1238 iCurvatureO = tmalloc(sizeof(long) * numCols);
1239 iCurvatureSigma = tmalloc(sizeof(long) * numCols);
1240 iCurvatureSigmaO = tmalloc(sizeof(long) * numCols);
1241 iChiSq = tmalloc(sizeof(long) * numCols);
1242 iChiSqO = tmalloc(sizeof(long) * numCols);
1243 iRmsResidual = tmalloc(sizeof(long) * numCols);
1244 iRmsResidualO = tmalloc(sizeof(long) * numCols);
1245 iSigLevel = tmalloc(sizeof(long) * numCols);
1246 iSigLevelO = tmalloc(sizeof(long) * numCols);
1247 iFitIsValid = tmalloc(sizeof(long) * numCols);
1248 iFitIsValidO = tmalloc(sizeof(long) * numCols);
1249 iFitLabel = tmalloc(sizeof(long) * numCols);
1250 iFitLabelO = tmalloc(sizeof(long) * numCols);
1251 iy = tmalloc(sizeof(long) * numCols);
1252 iySigma = tmalloc(sizeof(long) * numCols);
1253 iFit = tmalloc(sizeof(long) * numCols);
1254 iResidual = tmalloc(sizeof(long) * numCols);
1255
1256 for (colIndex = 0; colIndex < numCols; colIndex++) {
1257 ySymbols[colIndex] = NULL;
1258 coefUnits[colIndex] = tmalloc(sizeof(char *) * terms);
1259 iInterceptSigma[colIndex] = -1;
1260 iInterceptSigmaO[colIndex] = -1;
1261 iIntercept[colIndex] = -1;
1262 iInterceptO[colIndex] = -1;
1263 iInterceptSigma[colIndex] = -1;
1264 iInterceptSigmaO[colIndex] = -1;
1265 iSlope[colIndex] = -1;
1266 iSlopeO[colIndex] = -1;
1267 iSlopeSigma[colIndex] = -1;
1268 iSlopeSigmaO[colIndex] = -1;
1269 iCurvature[colIndex] = -1;
1270 iCurvatureO[colIndex] = -1;
1271 iCurvatureSigma[colIndex] = -1;
1272 iCurvatureSigmaO[colIndex] = -1;
1273 iChiSq[colIndex] = -1;
1274 iChiSqO[colIndex] = -1;
1275 iRmsResidual[colIndex] = -1;
1276 iRmsResidualO[colIndex] = -1;
1277 iSigLevel[colIndex] = -1;
1278 iSigLevelO[colIndex] = -1;
1279 iFitIsValid[colIndex] = -1;
1280 iFitIsValidO[colIndex] = -1;
1281 iFitLabel[colIndex] = -1;
1282 iFitLabelO[colIndex] = -1;
1283 iy[colIndex] = -1;
1284 iySigma[colIndex] = -1;
1285 iFit[colIndex] = -1;
1286 iResidual[colIndex] = -1;
1287 }
1288
1289 if (!SDDS_InitializeOutput(SDDSout, SDDS_BINARY, 0, NULL, "sddsmpfit output: fitted data", output) ||
1290 !SDDS_TransferColumnDefinition(SDDSout, SDDSin, xName, NULL) ||
1291 SDDS_GetColumnInformation(SDDSout, "symbol", &xSymbol, SDDS_GET_BY_NAME, xName) != SDDS_STRING ||
1292 (xSigmaName && !SDDS_TransferColumnDefinition(SDDSout, SDDSin, xSigmaName, NULL)))
1293 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1294
1295 for (colIndex = 0; colIndex < numCols; colIndex++) {
1296 if (!SDDS_TransferColumnDefinition(SDDSout, SDDSin, yNames[colIndex], NULL) ||
1297 SDDS_GetColumnInformation(SDDSout, "symbol", &ySymbols[colIndex], SDDS_GET_BY_NAME, yNames[colIndex]) != SDDS_STRING ||
1298 (ySigmaNames && !SDDS_TransferColumnDefinition(SDDSout, SDDSin, ySigmaNames[colIndex], NULL)))
1299 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1300 }
1301 if (!xSymbol || SDDS_StringIsBlank(xSymbol))
1302 xSymbol = xName;
1303 for (colIndex = 0; colIndex < numCols; colIndex++)
1304 if (!ySymbols[colIndex] || SDDS_StringIsBlank(ySymbols[colIndex]))
1305 ySymbols[colIndex] = yNames[colIndex];
1306 ix = SDDS_GetColumnIndex(SDDSout, xName);
1307 for (colIndex = 0; colIndex < numCols; colIndex++) {
1308 iy[colIndex] = SDDS_GetColumnIndex(SDDSout, yNames[colIndex]);
1309 if (ySigmaNames)
1310 iySigma[colIndex] = SDDS_GetColumnIndex(SDDSout, ySigmaNames[colIndex]);
1311 }
1312 if (xSigmaName)
1313 ixSigma = SDDS_GetColumnIndex(SDDSout, xSigmaName);
1314 if (SDDS_NumberOfErrors())
1315 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1316
1317 for (colIndex = 0; colIndex < numCols; colIndex++) {
1318 sprintf(buffer, "%sFit", yNames[colIndex]);
1319 sprintf(buffer1, "Fit[%s]", ySymbols[colIndex]);
1320 if (!SDDS_TransferColumnDefinition(SDDSout, SDDSin, yNames[colIndex], buffer) ||
1321 !SDDS_ChangeColumnInformation(SDDSout, "symbol", buffer1, SDDS_SET_BY_NAME, buffer))
1322 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1323 if ((iFit[colIndex] = SDDS_GetColumnIndex(SDDSout, buffer)) < 0)
1324 SDDS_Bomb("unable to get index of just-defined fit output column");
1325
1326 sprintf(buffer, "%sResidual", yNames[colIndex]);
1327 sprintf(buffer1, "Residual[%s]", ySymbols[colIndex]);
1328 if (!SDDS_TransferColumnDefinition(SDDSout, SDDSin, yNames[colIndex], buffer) ||
1329 !SDDS_ChangeColumnInformation(SDDSout, "symbol", buffer1, SDDS_SET_BY_NAME, buffer))
1330 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1331 if (!(iResidual[colIndex] = SDDS_GetColumnIndex(SDDSout, buffer)))
1332 SDDS_Bomb("unable to get index of just-defined residual output column");
1333
1334 if (sigmasValid && !ySigmaNames) {
1335 sprintf(buffer, "%sSigma", yNames[colIndex]);
1336 if (!SDDS_TransferColumnDefinition(SDDSout, SDDSin, yNames[colIndex], buffer))
1337 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1338 iySigma[colIndex] = SDDS_GetColumnIndex(SDDSout, buffer);
1339 if (ySymbols[colIndex] && !SDDS_StringIsBlank(ySymbols[colIndex])) {
1340 sprintf(buffer1, "Sigma[%s]", ySymbols[colIndex]);
1341 if (!SDDS_ChangeColumnInformation(SDDSout, "symbol", buffer1, SDDS_SET_BY_NAME, buffer))
1342 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1343 }
1344 }
1345
1346 if (!(coefUnits[colIndex] = makeCoefficientUnits(SDDSout, xName, yNames[colIndex], order, terms)))
1347 SDDS_Bomb("unable to make coefficient units");
1348 }
1349
1350 if (outputInfo && !SDDS_InitializeOutput(SDDSoutInfo, SDDS_BINARY, 0, NULL, "sddsmpfit output: fit information", outputInfo))
1351 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1352
1353 if (outputInfo) {
1354 if ((SDDS_DefineColumn(SDDSoutInfo, "Order", NULL, NULL, "Order of term in fit", NULL, SDDS_LONG, 0) < 0) ||
1355 SDDS_DefineParameter(SDDSoutInfo, "Basis", NULL, NULL, "Function basis for fit", NULL, SDDS_STRING, isChebyshev ? "Chebyshev T polynomials" : "ordinary polynomials") < 0)
1356 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1357
1358 if ((iTerms = SDDS_DefineParameter(SDDSoutInfo, "Terms", NULL, NULL, "Number of terms in fit", NULL, SDDS_LONG, NULL)) < 0)
1359 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1360
1361 if (SDDS_GetColumnInformation(SDDSout, "units", &xUnits, SDDS_GET_BY_NAME, xName) != SDDS_STRING)
1362 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1363 sprintf(buffer, "%sOffset", xName);
1364 sprintf(buffer1, "Offset of %s for fit", xName);
1365 if ((iOffset = SDDS_DefineParameter(SDDSoutInfo, buffer, NULL, xUnits, buffer1, NULL, SDDS_DOUBLE, NULL)) < 0)
1366 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1367 sprintf(buffer, "%sScale", xName);
1368 sprintf(buffer1, "Scale factor of %s for fit", xName);
1369 if ((iFactor = SDDS_DefineParameter(SDDSoutInfo, buffer, NULL, xUnits, buffer1, NULL, SDDS_DOUBLE, NULL)) < 0)
1370 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1371
1372 for (colIndex = 0; colIndex < numCols; colIndex++) {
1373
1374 sprintf(buffer1, "%sCoefficient", yNames[colIndex]);
1375 sprintf(buffer2, "%sCoefficientSigma", yNames[colIndex]);
1376 sprintf(buffer3, "%sCoefficientUnits", yNames[colIndex]);
1377
1378 if (SDDS_DefineColumn(SDDSoutInfo, buffer1, NULL, "[CoefficientUnits]", "Coefficient of term in fit", NULL, SDDS_DOUBLE, 0) < 0 ||
1379 (sigmasValid && SDDS_DefineColumn(SDDSoutInfo, buffer2, "$gs$r$ba$n", "[CoefficientUnits]", "sigma of coefficient of term in fit", NULL, SDDS_DOUBLE, 0) < 0))
1380 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1381
1382 iOrder = SDDS_GetColumnIndex(SDDSoutInfo, "Order");
1383 iCoefficient[colIndex] = SDDS_GetColumnIndex(SDDSoutInfo, buffer1);
1384 iCoefficientSigma[colIndex] = SDDS_GetColumnIndex(SDDSoutInfo, buffer2);
1385 iCoefficientUnits[colIndex] = SDDS_GetColumnIndex(SDDSoutInfo, buffer3);
1386
1387 sprintf(buffer1, "%sReducedChiSquared", yNames[colIndex]);
1388 sprintf(buffer2, "%sRmsResidual", yNames[colIndex]);
1389 sprintf(buffer3, "%sSignificanceLevel", yNames[colIndex]);
1390
1391 if ((iChiSq[colIndex] = SDDS_DefineParameter(SDDSoutInfo, buffer1, "$gh$r$a2$n/(N-M)", NULL,
1392 "Reduced chi-squared of fit",
1393 NULL, SDDS_DOUBLE, NULL)) < 0 ||
1394 SDDS_GetColumnInformation(SDDSout, "units", &yUnits, SDDS_GET_BY_NAME, yNames[colIndex]) != SDDS_STRING ||
1395 (iRmsResidual[colIndex] =
1396 SDDS_DefineParameter(SDDSoutInfo, buffer2, "$gs$r$bres$n", yUnits, "RMS residual of fit", NULL, SDDS_DOUBLE, NULL)) < 0 ||
1397 (iSigLevel[colIndex] = SDDS_DefineParameter(SDDSoutInfo, buffer3, NULL, NULL, "Probability that data is from fit function", NULL, SDDS_DOUBLE, NULL)) < 0)
1398 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1399 if (yUnits)
1400 free(yUnits);
1401
1402 sprintf(buffer, "%sFitIsValid", yNames[colIndex]);
1403 if ((iFitIsValid[colIndex] = SDDS_DefineParameter(SDDSoutInfo, buffer, NULL, NULL, NULL, NULL, SDDS_CHARACTER, NULL)) < 0)
1404 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1405
1406 if (!isChebyshev) {
1407
1408 sprintf(buffer, "%sSddsmpfitlabel", yNames[colIndex]);
1409 iFitLabel[colIndex] = SDDS_DefineParameter(SDDSoutInfo, buffer, NULL, NULL, NULL, NULL, SDDS_STRING, NULL);
1410 if ((i = coefficient_index(order, terms, 0)) >= 0) {
1411 sprintf(buffer, "%sIntercept", yNames[colIndex]);
1412 iIntercept[colIndex] = SDDS_DefineParameter(SDDSoutInfo, buffer, buffer, coefUnits[colIndex][i], "Intercept of fit", NULL, SDDS_DOUBLE, NULL);
1413 sprintf(buffer, "%sInterceptSigma", yNames[colIndex]);
1414 if (sigmasValid)
1415 iInterceptSigma[colIndex] = SDDS_DefineParameter(SDDSoutInfo, buffer, buffer, coefUnits[colIndex][i], "Sigma of intercept of fit", NULL, SDDS_DOUBLE, NULL);
1416 }
1417 sprintf(buffer, "%sSlope", yNames[colIndex]);
1418 if ((i = coefficient_index(order, terms, 1)) >= 0) {
1419 iSlope[colIndex] = SDDS_DefineParameter(SDDSoutInfo, buffer, buffer, coefUnits[colIndex][i], "Slope of fit", NULL, SDDS_DOUBLE, NULL);
1420 if (sigmasValid) {
1421 sprintf(buffer, "%sSlopeSigma", yNames[colIndex]);
1422 iSlopeSigma[colIndex] = SDDS_DefineParameter(SDDSoutInfo, buffer, buffer, coefUnits[colIndex][i], "Sigma of slope of fit", NULL, SDDS_DOUBLE, NULL);
1423 }
1424 }
1425 if ((i = coefficient_index(order, terms, 2)) >= 0) {
1426 sprintf(buffer, "%sCurvature", yNames[colIndex]);
1427 iCurvature[colIndex] = SDDS_DefineParameter(SDDSoutInfo, buffer, buffer, coefUnits[colIndex][i], "Curvature of fit", NULL, SDDS_DOUBLE, NULL);
1428 if (sigmasValid) {
1429 sprintf(buffer, "%sCurvatureSigma", yNames[colIndex]);
1430 iCurvatureSigma[colIndex] = SDDS_DefineParameter(SDDSoutInfo, buffer, buffer, coefUnits[colIndex][i], "Sigma of curvature of fit", NULL, SDDS_DOUBLE, NULL);
1431 }
1432 }
1433 if (SDDS_NumberOfErrors())
1434 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1435 }
1436 }
1437 }
1438
1439 if (SDDS_DefineParameter(SDDSout, "Basis", NULL, NULL, "Function basis for fit", NULL, SDDS_STRING, isChebyshev ? "Chebyshev T polynomials" : "ordinary polynomials") < 0)
1440 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1441
1442 if ((iTermsO = SDDS_DefineParameter(SDDSout, "Terms", NULL, NULL, "Number of terms in fit", NULL, SDDS_LONG, NULL)) < 0)
1443 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1444
1445 if (SDDS_GetColumnInformation(SDDSout, "units", &xUnits, SDDS_GET_BY_NAME, xName) != SDDS_STRING)
1446 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1447 sprintf(buffer, "%sOffset", xName);
1448 sprintf(buffer1, "Offset of %s for fit", xName);
1449 if ((iOffsetO = SDDS_DefineParameter(SDDSout, buffer, NULL, xUnits, buffer1, NULL, SDDS_DOUBLE, NULL)) < 0)
1450 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1451 sprintf(buffer, "%sScale", xName);
1452 sprintf(buffer1, "Scale factor of %s for fit", xName);
1453 if ((iFactorO = SDDS_DefineParameter(SDDSout, buffer, NULL, xUnits, buffer1, NULL, SDDS_DOUBLE, NULL)) < 0)
1454 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1455
1456 for (colIndex = 0; colIndex < numCols; colIndex++) {
1457
1458 sprintf(buffer1, "%sReducedChiSquared", yNames[colIndex]);
1459 sprintf(buffer2, "%sRmsResidual", yNames[colIndex]);
1460 sprintf(buffer3, "%sSignificanceLevel", yNames[colIndex]);
1461
1462 if ((iChiSqO[colIndex] = SDDS_DefineParameter(SDDSout, buffer1, "$gh$r$a2$n/(N-M)", NULL,
1463 "Reduced chi-squared of fit",
1464 NULL, SDDS_DOUBLE, NULL)) < 0 ||
1465 SDDS_GetColumnInformation(SDDSout, "units", &yUnits, SDDS_GET_BY_NAME, yNames[colIndex]) != SDDS_STRING ||
1466 (iRmsResidualO[colIndex] =
1467 SDDS_DefineParameter(SDDSout, buffer2, "$gs$r$bres$n", yUnits, "RMS residual of fit", NULL, SDDS_DOUBLE, NULL)) < 0 ||
1468 (iSigLevelO[colIndex] = SDDS_DefineParameter(SDDSout, buffer3, NULL, NULL, "Probability that data is from fit function", NULL, SDDS_DOUBLE, NULL)) < 0)
1469 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1470 if (yUnits)
1471 free(yUnits);
1472
1473 sprintf(buffer, "%sFitIsValid", yNames[colIndex]);
1474 if ((iFitIsValidO[colIndex] = SDDS_DefineParameter(SDDSout, buffer, NULL, NULL, NULL, NULL, SDDS_CHARACTER, NULL)) < 0)
1475 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1476
1477 if (!isChebyshev) {
1478 sprintf(buffer, "%sSddsmpfitlabel", yNames[colIndex]);
1479 iFitLabelO[colIndex] = SDDS_DefineParameter(SDDSout, buffer, NULL, NULL, NULL, NULL, SDDS_STRING, NULL);
1480 if ((i = coefficient_index(order, terms, 0)) >= 0) {
1481 sprintf(buffer, "%sIntercept", yNames[colIndex]);
1482 iInterceptO[colIndex] = SDDS_DefineParameter(SDDSout, buffer, buffer, coefUnits[colIndex][i], "Intercept of fit", NULL, SDDS_DOUBLE, NULL);
1483 sprintf(buffer, "%sInterceptSigma", yNames[colIndex]);
1484 if (sigmasValid)
1485 iInterceptSigmaO[colIndex] = SDDS_DefineParameter(SDDSout, buffer, buffer, coefUnits[colIndex][i], "Sigma of intercept of fit", NULL, SDDS_DOUBLE, NULL);
1486 }
1487 sprintf(buffer, "%sSlope", yNames[colIndex]);
1488 if ((i = coefficient_index(order, terms, 1)) >= 0) {
1489 iSlopeO[colIndex] = SDDS_DefineParameter(SDDSout, buffer, buffer, coefUnits[colIndex][i], "Slope of fit", NULL, SDDS_DOUBLE, NULL);
1490 if (sigmasValid) {
1491 sprintf(buffer, "%sSlopeSigma", yNames[colIndex]);
1492 iSlopeSigmaO[colIndex] = SDDS_DefineParameter(SDDSout, buffer, buffer, coefUnits[colIndex][i], "Sigma of slope of fit", NULL, SDDS_DOUBLE, NULL);
1493 }
1494 }
1495 if ((i = coefficient_index(order, terms, 2)) >= 0) {
1496 sprintf(buffer, "%sCurvature", yNames[colIndex]);
1497 iCurvatureO[colIndex] = SDDS_DefineParameter(SDDSout, buffer, buffer, coefUnits[colIndex][i], "Curvature of fit", NULL, SDDS_DOUBLE, NULL);
1498 if (sigmasValid) {
1499 sprintf(buffer, "%sCurvatureSigma", yNames[colIndex]);
1500 iCurvatureSigmaO[colIndex] = SDDS_DefineParameter(SDDSout, buffer, buffer, coefUnits[colIndex][i], "Sigma of curvature of fit", NULL, SDDS_DOUBLE, NULL);
1501 }
1502 }
1503 if (SDDS_NumberOfErrors())
1504 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1505 }
1506 }
1507
1508 if (copyParameters) {
1509 if (!SDDS_TransferAllParameterDefinitions(SDDSout, SDDSin, 0))
1510 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1511 if (outputInfo && !SDDS_TransferAllParameterDefinitions(SDDSoutInfo, SDDSin, 0))
1512 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1513 }
1514
1515 if ((outputInfo && !SDDS_WriteLayout(SDDSoutInfo)) || !SDDS_WriteLayout(SDDSout))
1516 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1517
1518 return coefUnits;
1519}
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_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_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_PrintErrors
void SDDS_PrintErrors(FILE *fp, int32_t mode)
Prints recorded error messages to a specified file stream.
Definition SDDS_utils.c:432
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_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

◆ ipower()

double ipower ( double x,
long n )

Evaluate a power function x^n.

This function returns ( (x - x_offset)/x_scale )^n.

Parameters
xThe point at which to evaluate the power.
nThe exponent.
Returns
The value of ((x - x_offset)/x_scale)^n.

Definition at line 113 of file lsfBasisFns.c.

113 {
114 x = (x - x_offset) / x_scale;
115 return (ipow(x, n));
116}

◆ main()

int main ( int argc,
char ** argv )

Definition at line 301 of file sddsmpfit.c.

301 {
302 double **y = NULL, **sy = NULL, **diff = NULL;
303 double *x = NULL, *sx = NULL;
304 double xOffset, xScaleFactor;
305 double *xOrig = NULL, **yOrig = NULL, *sxOrig = NULL, **syOrig = NULL, **sy0 = NULL;
306 long terms, normTerm, ip, ySigmasValid;
307 int64_t i, j, points, pointsOrig;
308 long symmetry, chebyshev, termsGiven;
309 double sigmas, minimumSigma;
310 long sigmasMode, sparseInterval;
311 double **coef = NULL, **coefSigma = NULL;
312 double *chi = NULL, xLow, xHigh, *rmsResidual = NULL;
313 char *xName = NULL, *yName = NULL, **yNames = NULL, *xSigmaName = NULL;
314 char **ySigmaNames = NULL, *ySigmaControlString = NULL;
315 char *input = NULL, *output = NULL;
316 SDDS_DATASET SDDSin, SDDSout, SDDSoutInfo;
317 long *isFit = NULL, iArg, modifySigmas, termIndex;
318 long generateSigmas, verbose, ignoreSigmas;
319 long outputInitialized, copyParameters = 0;
320 int32_t *order = NULL;
321 SCANNED_ARG *s_arg;
322 double xMin, xMax, revpowThreshold;
323 double rms_average(double *d_x, int64_t d_n);
324 char *infoFile = NULL;
325 char *fitLabelFormat = "%g";
326 static char fitLabelBuffer[SDDS_MAXLINE];
327 unsigned long pipeFlags, reviseOrders;
328 EVAL_PARAMETERS evalParameters;
329 long rangeFitOnly = 0;
330
331 long colIndex;
332 long cloDependentIndex = -1, numDependentItems;
333 int32_t numYNames;
334
335 SDDS_RegisterProgramName(argv[0]);
336 argc = scanargs(&s_arg, argc, argv);
337 if (argc < 2 || argc > (3 + N_OPTIONS)) {
338 fprintf(stderr, "usage: %s\n", USAGE);
339 fprintf(stderr, "%s%s", additional_help, additional_help2);
340 exit(EXIT_FAILURE);
341 }
342
343 input = output = NULL;
344 xName = yName = xSigmaName = ySigmaControlString = NULL;
345 yNames = ySigmaNames = NULL;
346 numDependentItems = 0;
347 modifySigmas = reviseOrders = chebyshev = 0;
348 order = NULL;
349 symmetry = NO_SYMMETRY;
350 xMin = xMax = 0;
351 generateSigmas = 0;
352 sigmasMode = -1;
353 sigmas = 1;
354 minimumSigma = 0;
355 sparseInterval = 1;
356 terms = 2;
357 verbose = ignoreSigmas = 0;
358 normTerm = -1;
359 xOffset = 0;
360 xScaleFactor = 1;
361 basis_fn = ipower;
362 basis_dfn = dipower;
363 pipeFlags = 0;
364 evalParameters.file = NULL;
365 infoFile = NULL;
366 termsGiven = 0;
367
368 for (iArg = 1; iArg < argc; iArg++) {
369 if (s_arg[iArg].arg_type == OPTION) {
370 switch (match_string(s_arg[iArg].list[0], option, N_OPTIONS, 0)) {
371 case CLO_MODIFYSIGMAS:
372 modifySigmas = 1;
373 break;
374 case CLO_ORDERS:
375 if (termsGiven)
376 SDDS_Bomb("give -order or -terms, not both");
377 if (s_arg[iArg].n_items < 2)
378 SDDS_Bomb("invalid -orders syntax");
379 order = tmalloc(sizeof(*order) * (terms = s_arg[iArg].n_items - 1));
380 for (i = 1; i < s_arg[iArg].n_items; i++) {
381 if (sscanf(s_arg[iArg].list[i], "%" SCNd32, order + i - 1) != 1)
382 SDDS_Bomb("unable to scan order from -orders list");
383 }
384 break;
385 case CLO_RANGE:
386 rangeFitOnly = 0;
387 if ((s_arg[iArg].n_items != 3 && s_arg[iArg].n_items != 4) || 1 != sscanf(s_arg[iArg].list[1], "%lf", &xMin) || 1 != sscanf(s_arg[iArg].list[2], "%lf", &xMax) || xMin >= xMax)
388 SDDS_Bomb("incorrect -range syntax");
389 if (s_arg[iArg].n_items == 4) {
390 if (strncmp(str_tolower(s_arg[iArg].list[3]), "fitonly", strlen(s_arg[iArg].list[3])) == 0) {
391 rangeFitOnly = 1;
392 } else
393 SDDS_Bomb("incorrect -range syntax");
394 }
395 break;
396 case CLO_GENERATESIGMAS:
397 generateSigmas = FLGS_GENERATESIGMAS;
398 if (s_arg[iArg].n_items > 1) {
399 if (s_arg[iArg].n_items != 2)
400 SDDS_Bomb("incorrect -generateSigmas synax");
401 if (strncmp(s_arg[iArg].list[1], "keepsmallest", strlen(s_arg[iArg].list[1])) == 0)
402 generateSigmas |= FLGS_KEEPSMALLEST;
403 if (strncmp(s_arg[iArg].list[1], "keeplargest", strlen(s_arg[iArg].list[1])) == 0)
404 generateSigmas |= FLGS_KEEPLARGEST;
405 if ((generateSigmas & FLGS_KEEPSMALLEST) && (generateSigmas & FLGS_KEEPLARGEST))
406 SDDS_Bomb("ambiguous -generateSigmas synax");
407 }
408 break;
409 case CLO_TERMS:
410 if (order)
411 SDDS_Bomb("give -order or -terms, not both");
412 if (s_arg[iArg].n_items != 2 || sscanf(s_arg[iArg].list[1], "%ld", &terms) != 1)
413 SDDS_Bomb("invalid -terms syntax");
414 termsGiven = 1;
415 break;
416 case CLO_XOFFSET:
417 if (s_arg[iArg].n_items != 2 || sscanf(s_arg[iArg].list[1], "%lf", &xOffset) != 1)
418 SDDS_Bomb("invalid -xOffset syntax");
419 break;
420 case CLO_SYMMETRY:
421 if (s_arg[iArg].n_items == 2) {
422 if ((symmetry = match_string(s_arg[iArg].list[1], symmetry_options, N_SYMMETRY_OPTIONS, 0)) < 0)
423 SDDS_Bomb("unknown option used with -symmetry");
424 } else
425 SDDS_Bomb("incorrect -symmetry syntax");
426 break;
427 case CLO_SIGMAS:
428 if (s_arg[iArg].n_items != 3)
429 SDDS_Bomb("incorrect -sigmas syntax");
430 if (sscanf(s_arg[iArg].list[1], "%lf", &sigmas) != 1)
431 SDDS_Bomb("couldn't scan value for -sigmas");
432 if ((sigmasMode = match_string(s_arg[iArg].list[2], sigmas_options, N_SIGMAS_OPTIONS, 0)) < 0)
433 SDDS_Bomb("unrecognized -sigmas mode");
434 break;
435 case CLO_MINSIGMA:
436 if (s_arg[iArg].n_items != 2)
437 SDDS_Bomb("incorrect -minimumSigma syntax");
438 if (sscanf(s_arg[iArg].list[1], "%lf", &minimumSigma) != 1)
439 SDDS_Bomb("couldn't scan value for -minimumSigma");
440 break;
441 case CLO_SPARSE:
442 if (s_arg[iArg].n_items != 2)
443 SDDS_Bomb("incorrect -sparse syntax");
444 if (sscanf(s_arg[iArg].list[1], "%ld", &sparseInterval) != 1)
445 SDDS_Bomb("couldn't scan value for -sparse");
446 if (sparseInterval < 1)
447 SDDS_Bomb("invalid -sparse value");
448 break;
449 case CLO_VERBOSE:
450 verbose = 1;
451 break;
452 case CLO_NORMALIZE:
453 normTerm = 0;
454 if (s_arg[iArg].n_items > 2 ||
455 (s_arg[iArg].n_items == 2 && sscanf(s_arg[iArg].list[1], "%ld", &normTerm) != 1) ||
456 normTerm < 0)
457 SDDS_Bomb("invalid -normalize syntax");
458 break;
459 case CLO_REVISEORDERS:
460 revpowThreshold = 0.1;
461 s_arg[iArg].n_items -= 1;
462 if (!scanItemList(&reviseOrders, s_arg[iArg].list + 1, &s_arg[iArg].n_items, 0,
463 "threshold", SDDS_DOUBLE, &revpowThreshold, 1, 0,
464 "verbose", -1, NULL, 1, REVPOW_VERBOSE, NULL))
465 SDDS_Bomb("invalid -reviseOrders syntax");
466 reviseOrders |= REVPOW_ACTIVE;
467 revpowThreshold = fabs(revpowThreshold);
468 break;
469 case CLO_CHEBYSHEV:
470 if (s_arg[iArg].n_items > 2 ||
471 (s_arg[iArg].n_items == 2 && strncmp(s_arg[iArg].list[1], "convert", strlen(s_arg[iArg].list[1])) != 0))
472 SDDS_Bomb("invalid -chebyshev syntax");
473 chebyshev = s_arg[iArg].n_items;
474 basis_fn = tcheby;
475 basis_dfn = dtcheby;
476 break;
477 case CLO_XFACTOR:
478 if (s_arg[iArg].n_items != 2 ||
479 sscanf(s_arg[iArg].list[1], "%lf", &xScaleFactor) != 1 || xScaleFactor == 0)
480 SDDS_Bomb("invalid -xFactor syntax");
481 break;
482 case CLO_INDEPENDENT:
483 if (s_arg[iArg].n_items != 2)
484 SDDS_Bomb("invalid -independent syntax");
485 xName = s_arg[iArg].list[1];
486 break;
487 case CLO_DEPENDENT:
488 numDependentItems = s_arg[iArg].n_items - 1;
489 cloDependentIndex = iArg;
490 if (numDependentItems < 1)
491 SDDS_Bomb("invalid -dependent syntax");
492 break;
493 case CLO_SIGMAINDEPENDENT:
494 if (s_arg[iArg].n_items != 2)
495 SDDS_Bomb("invalid -sigmaIndependent syntax");
496 xSigmaName = s_arg[iArg].list[1];
497 break;
498 case CLO_SIGMADEPENDENT:
499 if (s_arg[iArg].n_items != 2)
500 SDDS_Bomb("invalid -sigmaDependent syntax");
501 ySigmaControlString = s_arg[iArg].list[1];
502 break;
503 case CLO_FITLABELFORMAT:
504 if (s_arg[iArg].n_items != 2)
505 SDDS_Bomb("invalid -fitLabelFormat syntax");
506 fitLabelFormat = s_arg[iArg].list[1];
507 break;
508 case CLO_PIPE:
509 if (!processPipeOption(s_arg[iArg].list + 1, s_arg[iArg].n_items - 1, &pipeFlags))
510 SDDS_Bomb("invalid -pipe syntax");
511 break;
512 case CLO_INFOFILE:
513 if (s_arg[iArg].n_items != 2)
514 SDDS_Bomb("invalid -infoFile syntax");
515 infoFile = s_arg[iArg].list[1];
516 break;
517 case CLO_EVALUATE:
518 if (s_arg[iArg].n_items < 2)
519 SDDS_Bomb("invalid -evaluate syntax");
520 evalParameters.file = s_arg[iArg].list[1];
521 s_arg[iArg].n_items -= 2;
522 s_arg[iArg].list += 2;
523 if (!scanItemList(&evalParameters.flags, s_arg[iArg].list, &s_arg[iArg].n_items, 0,
524 "begin", SDDS_DOUBLE, &evalParameters.begin, 1, EVAL_BEGIN_GIVEN,
525 "end", SDDS_DOUBLE, &evalParameters.end, 1, EVAL_END_GIVEN,
526 "number", SDDS_LONG64, &evalParameters.number, 1, EVAL_NUMBER_GIVEN, NULL))
527 SDDS_Bomb("invalid -evaluate syntax");
528 break;
529 case CLO_COPYPARAMETERS:
530 copyParameters = 1;
531 break;
532 default:
533 bomb("unknown switch", USAGE);
534 break;
535 }
536 } else {
537 if (input == NULL)
538 input = s_arg[iArg].list[0];
539 else if (output == NULL)
540 output = s_arg[iArg].list[0];
541 else
542 SDDS_Bomb("too many filenames");
543 }
544 }
545
546 processFilenames("sddsmpfit", &input, &output, pipeFlags, 0, NULL);
547
548 if (symmetry && order)
549 SDDS_Bomb("can't specify both -symmetry and -orders");
550 if (!xName || !numDependentItems)
551 SDDS_Bomb("you must specify a column name for x and y");
552 if (modifySigmas && !xSigmaName)
553 SDDS_Bomb("you must specify x sigmas with -modifySigmas");
554 if (generateSigmas) {
555 if (modifySigmas)
556 SDDS_Bomb("you can't specify both -generateSigmas and -modifySigmas");
557 }
558 if (ySigmaControlString) {
559 if (sigmasMode != -1)
560 SDDS_Bomb("you can't specify both -sigmas and a y sigma name");
561 }
562 ySigmasValid = 0;
563 if (sigmasMode != -1 || generateSigmas || ySigmaControlString || modifySigmas)
564 ySigmasValid = 1;
565
566 if (normTerm >= 0 && normTerm >= terms)
567 SDDS_Bomb("can't normalize to that term--not that many terms");
568 if (reviseOrders && !(sigmasMode != -1 || generateSigmas || ySigmaNames))
569 SDDS_Bomb("can't use -reviseOrders unless a y sigma or -generateSigmas is given");
570
571 if (symmetry == EVEN_SYMMETRY) {
572 order = tmalloc(sizeof(*order) * terms);
573 for (i = 0; i < terms; i++)
574 order[i] = 2 * i;
575 } else if (symmetry == ODD_SYMMETRY) {
576 order = tmalloc(sizeof(*order) * terms);
577 for (i = 0; i < terms; i++)
578 order[i] = 2 * i + 1;
579 } else if (!order) {
580 order = tmalloc(sizeof(*order) * terms);
581 for (i = 0; i < terms; i++)
582 order[i] = i;
583 }
584
585 if (!SDDS_InitializeInput(&SDDSin, input))
586 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
587 outputInitialized = 0;
588 yNames = ResolveColumnNames(&SDDSin, s_arg[cloDependentIndex].list + 1, numDependentItems, &numYNames);
589 if (ySigmaControlString != NULL)
590 ySigmaNames = GenerateYSigmaNames(ySigmaControlString, yNames, numYNames);
591
592 checkInputFile(&SDDSin, xName, yNames, xSigmaName, ySigmaNames, numYNames);
593 sy0 = tmalloc(sizeof(double *) * numYNames);
594 y = tmalloc(sizeof(double *) * numYNames);
595 sy = tmalloc(sizeof(double *) * numYNames);
596 isFit = tmalloc(sizeof(long) * numYNames);
597 chi = tmalloc(sizeof(double) * numYNames);
598 coef = tmalloc(sizeof(double *) * numYNames);
599 coefSigma = tmalloc(sizeof(double *) * numYNames);
600 for (colIndex = 0; colIndex < numYNames; colIndex++) {
601 coef[colIndex] = tmalloc(sizeof(double) * terms);
602 coefSigma[colIndex] = tmalloc(sizeof(double) * terms);
603 }
604 iCoefficient = tmalloc(sizeof(long) * numYNames);
605 iCoefficientSigma = tmalloc(sizeof(long) * numYNames);
606 iCoefficientUnits = tmalloc(sizeof(long) * numYNames);
607
608 while (SDDS_ReadPage(&SDDSin) > 0) {
609 if ((points = SDDS_CountRowsOfInterest(&SDDSin)) < terms) {
610 /* probably should emit an empty page here */
611 continue;
612 }
613 if (!(x = SDDS_GetColumnInDoubles(&SDDSin, xName))) {
614 fprintf(stderr, "error: unable to read column %s\n", xName);
615 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
616 }
617 for (i = 0; i < numYNames; i++) {
618 if (!(y[i] = SDDS_GetColumnInDoubles(&SDDSin, yNames[i]))) {
619 fprintf(stderr, "error: unable to read column %s\n", yNames[i]);
620 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
621 }
622 }
623 sx = NULL;
624 if (xSigmaName && !(sx = SDDS_GetColumnInDoubles(&SDDSin, xSigmaName))) {
625 fprintf(stderr, "error: unable to read column %s\n", xSigmaName);
626 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
627 }
628 for (colIndex = 0; colIndex < numYNames; colIndex++)
629 sy0[colIndex] = tmalloc(sizeof(double) * points);
630 if (ySigmaNames) {
631 for (i = 0; i < numYNames; i++) {
632 if (!(sy0[i] = SDDS_GetColumnInDoubles(&SDDSin, ySigmaNames[i]))) {
633 fprintf(stderr, "error: unable to read column %s\n", ySigmaNames[i]);
634 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
635 }
636 }
637 }
638
639 if (minimumSigma > 0) {
640 int64_t j;
641 for (i = 0; i < numYNames; i++) {
642 for (j = 0; j < points; j++)
643 if (sy0[i][j] < minimumSigma)
644 sy0[i][j] = minimumSigma;
645 }
646 }
647
648 if (xMin != xMax || sparseInterval != 1) {
649 xOrig = tmalloc(sizeof(*xOrig) * points);
650 yOrig = tmalloc(sizeof(*yOrig) * numYNames);
651 for (colIndex = 0; colIndex < numYNames; colIndex++)
652 yOrig[colIndex] = tmalloc(sizeof(double) * points);
653 if (sx)
654 sxOrig = tmalloc(sizeof(*sxOrig) * points);
655 if (ySigmasValid) {
656 syOrig = tmalloc(sizeof(*syOrig) * numYNames);
657 for (colIndex = 0; colIndex < numYNames; colIndex++)
658 syOrig[colIndex] = tmalloc(sizeof(double) * points);
659 }
660 pointsOrig = points;
661 for (i = j = 0; i < points; i++) {
662 xOrig[i] = x[i];
663 if (sx)
664 sxOrig[i] = sx[i];
665 for (colIndex = 0; colIndex < numYNames; colIndex++) {
666 yOrig[colIndex][i] = y[colIndex][i];
667 if (ySigmasValid)
668 syOrig[colIndex][i] = sy0[colIndex][i];
669 }
670 }
671 if (xMin != xMax) {
672 for (i = j = 0; i < points; i++) {
673 if (xOrig[i] <= xMax && xOrig[i] >= xMin) {
674 x[j] = xOrig[i];
675 for (colIndex = 0; colIndex < numYNames; colIndex++) {
676 y[colIndex][j] = yOrig[colIndex][i];
677 if (ySigmasValid)
678 sy0[colIndex][j] = syOrig[colIndex][i];
679 }
680 if (sx)
681 sx[j] = sxOrig[i];
682 j++;
683 }
684 }
685 points = j;
686 }
687 if (sparseInterval != 1) {
688 for (i = j = 0; i < points; i++) {
689 if (i % sparseInterval == 0) {
690 x[j] = x[i];
691 for (colIndex = 0; colIndex < numYNames; colIndex++) {
692 y[colIndex][j] = y[colIndex][i];
693 if (ySigmasValid)
694 sy0[colIndex][j] = sy0[colIndex][i];
695 }
696 if (sx)
697 sx[j] = sx[i];
698 j++;
699 }
700 }
701 points = j;
702 }
703 } else {
704 xOrig = x;
705 yOrig = y;
706 sxOrig = sx;
707 syOrig = sy0;
708 pointsOrig = points;
709 }
710
711 find_min_max(&xLow, &xHigh, x, points);
712
713 if (sigmasMode == ABSOLUTE_SIGMAS) {
714 for (colIndex = 0; colIndex < numYNames; colIndex++) {
715 for (i = 0; i < points; i++)
716 sy0[colIndex][i] = sigmas;
717 if (sy0[colIndex] != syOrig[colIndex])
718 for (i = 0; i < pointsOrig; i++)
719 syOrig[colIndex][i] = sigmas;
720 }
721 } else if (sigmasMode == FRACTIONAL_SIGMAS) {
722 for (colIndex = 0; colIndex < numYNames; colIndex++) {
723 for (i = 0; i < points; i++)
724 sy0[colIndex][i] = sigmas * fabs(y[colIndex][i]);
725 if (sy0[colIndex] != syOrig[colIndex])
726 for (i = 0; i < pointsOrig; i++)
727 syOrig[colIndex][i] = fabs(yOrig[colIndex][i]) * sigmas;
728 }
729 }
730
731 for (i = 0; i < numYNames; i++) {
732 if (minimumSigma > 0) {
733 int64_t j;
734 for (j = 0; j < points; j++)
735 if (sy0[i][j] < minimumSigma)
736 sy0[i][j] = minimumSigma;
737 }
738 }
739
740 if (!ySigmasValid || generateSigmas)
741 for (colIndex = 0; colIndex < numYNames; colIndex++) {
742 for (i = 0; i < points; i++)
743 sy0[colIndex][i] = 1;
744 }
745 else
746 for (i = 0; i < points; i++)
747 for (colIndex = 0; colIndex < numYNames; colIndex++) {
748 if (sy0[colIndex][i] == 0)
749 SDDS_Bomb("y sigma = 0 for one or more points.");
750 }
751
752 diff = tmalloc(sizeof(*diff) * numYNames);
753 sy = tmalloc(sizeof(*sy) * numYNames);
754 for (colIndex = 0; colIndex < numYNames; colIndex++) {
755 diff[colIndex] = tmalloc(sizeof(double) * points);
756 sy[colIndex] = tmalloc(sizeof(double) * points);
757 }
758
759 for (i = 0; i < points; i++) {
760 for (colIndex = 0; colIndex < numYNames; colIndex++)
761 sy[colIndex][i] = sy0[colIndex][i];
762 }
763
764 set_argument_offset(xOffset);
765 set_argument_scale(xScaleFactor);
766 if (chebyshev) {
767 xOffset = (xHigh + xLow) / 2;
768 set_argument_offset(xOffset);
769 set_argument_scale(xScaleFactor = (xHigh - xLow) / 2);
770 }
771
772 if (generateSigmas || modifySigmas) {
773 /* do an initial fit */
774 for (colIndex = 0; colIndex < numYNames; colIndex++) {
775 isFit[colIndex] = lsfg(x, y[colIndex], sy[colIndex], points, terms, order, coef[colIndex], coefSigma[colIndex], &chi[colIndex], diff[colIndex], basis_fn);
776 if (!isFit[colIndex]) {
777 fprintf(stderr, "Column %s: ", yNames[colIndex]);
778 SDDS_Bomb("initial fit failed.");
779 }
780 if (verbose) {
781 fprintf(stderr, "Column %s: ", yNames[colIndex]);
782 fputs("initial_fit:", stderr);
783 print_coefs(stderr, xOffset, xScaleFactor, chebyshev, coef[colIndex], NULL, order, terms, chi[colIndex], normTerm, "");
784 fprintf(stderr, "unweighted rms deviation from fit: %21.15le\n", rms_average(diff[colIndex], points));
785 }
786 if (modifySigmas) {
787 if (!ySigmasValid) {
788 for (i = 0; i < points; i++)
789 sy[colIndex][i] = fabs(eval_sum(basis_dfn, coef[colIndex], order, terms, x[i]) * sx[i]);
790 } else
791 for (i = 0; i < points; i++) {
792 sy[colIndex][i] = sqrt(sqr(sy0[colIndex][i]) + sqr(eval_sum(basis_dfn, coef[colIndex], order, terms, x[i]) * sx[i]));
793 }
794 }
795 if (generateSigmas) {
796 double sigma;
797 for (i = sigma = 0; i < points; i++) {
798 sigma += sqr(diff[colIndex][i]);
799 }
800 sigma = sqrt(sigma / (points - terms));
801 for (i = 0; i < points; i++) {
802 if (generateSigmas & FLGS_KEEPSMALLEST) {
803 if (sigma < sy[colIndex][i])
804 sy[colIndex][i] = sigma;
805 } else if (generateSigmas & FLGS_KEEPLARGEST) {
806 if (sigma > sy[colIndex][i])
807 sy[colIndex][i] = sigma;
808 } else {
809 sy[colIndex][i] = sigma;
810 }
811 }
812 for (i = 0; i < pointsOrig; i++) {
813 if (generateSigmas & FLGS_KEEPSMALLEST) {
814 if (sigma < sy0[colIndex][i])
815 sy0[colIndex][i] = sigma;
816 } else if (generateSigmas & FLGS_KEEPLARGEST) {
817 if (sigma > sy0[colIndex][i])
818 sy0[colIndex][i] = sigma;
819 } else {
820 sy0[colIndex][i] = sigma;
821 }
822 }
823 }
824 }
825 }
826
827 if (reviseOrders & REVPOW_ACTIVE) {
828 double bestChi;
829 long bestTerms, newBest;
830 int32_t *bestOrder;
831
832 bestTerms = terms;
833 bestOrder = tmalloc(sizeof(*bestOrder) * bestTerms);
834 for (ip = 0; ip < terms; ip++)
835 bestOrder[ip] = order[ip];
836 /* do a fit */
837 for (colIndex = 0; colIndex < numYNames; colIndex++) {
838 isFit[colIndex] = lsfg(x, y[colIndex], sy[colIndex], points, bestTerms, bestOrder, coef[colIndex], coefSigma[colIndex], &bestChi, diff[colIndex], basis_fn);
839 if (!isFit[colIndex]) {
840 fprintf(stderr, "Column %s: ", yNames[colIndex]);
841 SDDS_Bomb("revise-orders fit failed.");
842 if (reviseOrders & REVPOW_VERBOSE) {
843 fprintf(stderr, "Column %s: ", yNames[colIndex]);
844 fputs("fit to revise orders:", stderr);
845 print_coefs(stderr, xOffset, xScaleFactor, chebyshev, coef[colIndex], (ySigmasValid ? coefSigma[colIndex] : NULL), bestOrder, bestTerms, bestChi, normTerm, "");
846 fprintf(stderr, "unweighted rms deviation from fit: %21.15le\n", rms_average(diff[colIndex], points));
847 }
848 }
849
850 do {
851 newBest = 0;
852 terms = bestTerms - 1;
853 for (ip = bestTerms - 1; ip >= 0; ip--) {
854 for (i = j = 0; i < bestTerms; i++)
855 if (i != ip)
856 order[j++] = bestOrder[i];
857 isFit[colIndex] = lsfg(x, y[colIndex], sy[colIndex], points, terms, order, coef[colIndex], coefSigma[colIndex], &chi[colIndex], diff[colIndex], basis_fn);
858 if (!isFit[colIndex]) {
859 fprintf(stderr, "Column %s: ", yNames[colIndex]);
860 SDDS_Bomb("revise-orders fit failed.");
861 }
862 if (reviseOrders & REVPOW_VERBOSE) {
863 fprintf(stderr, "Column %s: ", yNames[colIndex]);
864 fputs("new trial fit:", stderr);
865 print_coefs(stderr, xOffset, xScaleFactor, chebyshev, coef[colIndex], (ySigmasValid ? coefSigma[colIndex] : NULL), order, terms, chi[colIndex], normTerm, "");
866 fprintf(stderr, "unweighted rms deviation from fit: %21.15le\n", rms_average(diff[colIndex], points));
867 }
868 if (chi[colIndex] - bestChi < revpowThreshold) {
869 bestChi = chi[colIndex];
870 bestTerms = terms;
871 newBest = 1;
872 for (i = 0; i < terms; i++)
873 bestOrder[i] = order[i];
874 if (reviseOrders & REVPOW_VERBOSE) {
875 fputs("new best fit:", stderr);
876 print_coefs(stderr, xOffset, xScaleFactor, chebyshev, coef[colIndex], (ySigmasValid ? coefSigma[colIndex] : NULL), bestOrder, bestTerms, bestChi, normTerm, "");
877 fprintf(stderr, "unweighted rms deviation from fit: %21.15le\n", rms_average(diff[colIndex], points));
878 }
879 break;
880 }
881 }
882 if (bestTerms == 1)
883 break;
884 } while (newBest);
885 terms = bestTerms;
886 for (ip = 0; ip < terms; ip++)
887 order[ip] = bestOrder[ip];
888 free(bestOrder);
889 reviseOrders = 0;
890 }
891 }
892
893 if (!outputInitialized) {
894 initializeOutputFile(&SDDSout, &SDDSoutInfo, output, infoFile, &SDDSin, input, xName, yNames, xSigmaName, ySigmaNames, ySigmasValid, order, terms, chebyshev, numYNames, copyParameters);
895 free(output);
896 outputInitialized = 1;
897 }
898 if (evalParameters.file)
899 setupEvaluationFile(&evalParameters, xName, yNames, numYNames, &SDDSin);
900
901 rmsResidual = tmalloc(sizeof(double) * numYNames);
902 for (colIndex = 0; colIndex < numYNames; colIndex++) {
903 isFit[colIndex] = lsfg(x, y[colIndex], sy[colIndex], points, terms, order, coef[colIndex], coefSigma[colIndex], &chi[colIndex], diff[colIndex], basis_fn);
904 if (isFit[colIndex]) {
905 rmsResidual[colIndex] = rms_average(diff[colIndex], points);
906 if (verbose) {
907 fprintf(stderr, "Column: %s\n", yNames[colIndex]);
908 print_coefs(stderr, xOffset, xScaleFactor, chebyshev, coef[colIndex], (ySigmasValid ? coefSigma[colIndex] : NULL), order, terms, chi[colIndex], normTerm, "");
909 fprintf(stderr, "unweighted rms deviation from fit: %21.15le\n", rmsResidual[colIndex]);
910 }
911 } else if (verbose)
912 fprintf(stderr, "fit failed for %s.\n", yNames[colIndex]);
913
914 if (evalParameters.file)
915 makeEvaluationTable(&evalParameters, x, points, coef[colIndex], order, terms, xName, yNames, numYNames, colIndex);
916 }
917
918 if (outputInitialized) {
919 if (!SDDS_StartPage(&SDDSout, rangeFitOnly ? pointsOrig : points) ||
920 (infoFile && !SDDS_StartPage(&SDDSoutInfo, terms)))
921 bomb("A", NULL);
922 if (copyParameters) {
923 if (!SDDS_CopyParameters(&SDDSout, &SDDSin))
924 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
925 if (infoFile && !SDDS_CopyParameters(&SDDSoutInfo, &SDDSin))
926 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
927 }
928 if (!SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_INDEX, rangeFitOnly ? xOrig : x, rangeFitOnly ? pointsOrig : points, ix) ||
929 (infoFile && !SDDS_SetColumnFromLongs(&SDDSoutInfo, SDDS_SET_BY_INDEX, order, terms, iOrder)))
930 bomb("B", NULL);
931 for (colIndex = 0; colIndex < numYNames; colIndex++) {
932 if (rangeFitOnly) {
933 double *residual, rmsResidual0;
934 compareOriginalToFit(xOrig, yOrig[colIndex], &residual, pointsOrig, &rmsResidual0, coef[colIndex], order, terms);
935 if (!SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_INDEX, yOrig[colIndex], pointsOrig, iy[colIndex]) ||
936 !SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_INDEX, residual, pointsOrig, iResidual[colIndex]))
937 bomb("C", NULL);
938 for (i = 0; i < pointsOrig; i++)
939 residual[i] = yOrig[colIndex][i] - residual[i]; /* computes fit values from residual and y */
940 if (!SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_INDEX, residual, pointsOrig, iFit[colIndex]))
941 bomb("D", NULL);
942 free(residual);
943 } else {
944 for (i = 0; i < points; i++)
945 diff[colIndex][i] = -diff[colIndex][i]; /* convert from (Fit-y) to (y-Fit) to get residual */
946 if (!SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_INDEX, y[colIndex], points, iy[colIndex]) ||
947 !SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_INDEX, diff[colIndex], points, iResidual[colIndex]))
948 bomb("C", NULL);
949 for (i = 0; i < points; i++)
950 diff[colIndex][i] = y[colIndex][i] - diff[colIndex][i]; /* computes fit values from residual and y */
951 if (!SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_INDEX, diff[colIndex], points, iFit[colIndex]))
952 bomb("D", NULL);
953 }
954 }
955 if (ixSigma != -1 &&
956 !SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_INDEX, rangeFitOnly ? sxOrig : sx, rangeFitOnly ? pointsOrig : points, ixSigma))
957 bomb("E", NULL);
958 for (colIndex = 0; colIndex < numYNames; colIndex++) {
959 if (ySigmasValid && iySigma[colIndex] != -1 &&
960 !SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_INDEX, rangeFitOnly ? syOrig[colIndex] : sy[colIndex], rangeFitOnly ? pointsOrig : points, iySigma[colIndex]))
961 bomb("F", NULL);
962
963 if (infoFile) {
964 termIndex = coefficient_index(order, terms, 0);
965 if (iIntercept[colIndex] != -1 &&
966 !SDDS_SetParameters(&SDDSoutInfo, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE, iIntercept[colIndex], coef[colIndex][termIndex], -1))
967 bomb("G", NULL);
968 if (iInterceptSigma[colIndex] != -1 &&
969 !SDDS_SetParameters(&SDDSoutInfo, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE, iInterceptSigma[colIndex], coefSigma[colIndex][termIndex], -1))
970 bomb("H", NULL);
971
972 termIndex = coefficient_index(order, terms, 1);
973 if (iSlope[colIndex] != -1 &&
974 !SDDS_SetParameters(&SDDSoutInfo, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE, iSlope[colIndex], coef[colIndex][termIndex], -1))
975 bomb("I", NULL);
976 if (iSlopeSigma[colIndex] != -1 &&
977 !SDDS_SetParameters(&SDDSoutInfo, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE, iSlopeSigma[colIndex], coefSigma[colIndex][termIndex], -1))
978 bomb("J", NULL);
979
980 termIndex = coefficient_index(order, terms, 2);
981 if (iCurvature[colIndex] != -1 &&
982 !SDDS_SetParameters(&SDDSoutInfo, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE, iCurvature[colIndex], coef[colIndex][termIndex], -1))
983 bomb("K", NULL);
984 if (iCurvatureSigma[colIndex] != -1 &&
985 !SDDS_SetParameters(&SDDSoutInfo, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE, iCurvatureSigma[colIndex], coefSigma[colIndex][termIndex], -1))
986 bomb("L", NULL);
987 if (iFitLabel[colIndex] != -1) {
988 makeFitLabel(fitLabelBuffer, SDDS_MAXLINE, fitLabelFormat, coef[colIndex], order, terms, colIndex);
989 if (!SDDS_SetParameters(&SDDSoutInfo, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE, iFitLabel[colIndex], fitLabelBuffer, -1))
990 bomb("M", NULL);
991 }
992 if (!SDDS_SetColumnFromDoubles(&SDDSoutInfo, SDDS_SET_BY_INDEX, coef[colIndex], terms, iCoefficient[colIndex]) ||
993 (ySigmasValid &&
994 !SDDS_SetColumnFromDoubles(&SDDSoutInfo, SDDS_SET_BY_INDEX, coefSigma[colIndex], terms, iCoefficientSigma[colIndex])))
995 bomb("N", NULL);
996 if (!SDDS_SetParameters(&SDDSoutInfo, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE,
997 iRmsResidual[colIndex], rmsResidual[colIndex], iChiSq[colIndex], chi[colIndex],
998 iTerms, terms, iSigLevel[colIndex], ChiSqrSigLevel(chi[colIndex], points - terms),
999 iOffset, xOffset, iFactor, xScaleFactor, iFitIsValid[colIndex], isFit[colIndex] ? 'y' : 'n', -1))
1000 bomb("O", NULL);
1001 }
1002
1003 termIndex = coefficient_index(order, terms, 0);
1004 if (iInterceptO[colIndex] != -1 &&
1005 !SDDS_SetParameters(&SDDSout, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE, iInterceptO[colIndex], coef[colIndex][termIndex], -1))
1006 bomb("G", NULL);
1007 if (iInterceptSigmaO[colIndex] != -1 &&
1008 !SDDS_SetParameters(&SDDSout, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE, iInterceptSigmaO[colIndex], coefSigma[colIndex][termIndex], -1))
1009 bomb("H", NULL);
1010
1011 termIndex = coefficient_index(order, terms, 1);
1012 if (iSlopeO[colIndex] != -1 &&
1013 !SDDS_SetParameters(&SDDSout, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE, iSlopeO[colIndex], coef[colIndex][termIndex], -1))
1014 bomb("I", NULL);
1015 if (iSlopeSigmaO[colIndex] != -1 &&
1016 !SDDS_SetParameters(&SDDSout, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE, iSlopeSigmaO[colIndex], coefSigma[colIndex][termIndex], -1))
1017 bomb("J", NULL);
1018
1019 termIndex = coefficient_index(order, terms, 2);
1020 if (iCurvatureO[colIndex] != -1 &&
1021 !SDDS_SetParameters(&SDDSout, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE, iCurvatureO[colIndex], coef[colIndex][termIndex], -1))
1022 bomb("K", NULL);
1023 if (iCurvatureSigmaO[colIndex] != -1 &&
1024 !SDDS_SetParameters(&SDDSout, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE, iCurvatureSigmaO[colIndex], coefSigma[colIndex][termIndex], -1))
1025 bomb("L", NULL);
1026 if (iFitLabelO[colIndex] != -1) {
1027 makeFitLabel(fitLabelBuffer, SDDS_MAXLINE, fitLabelFormat, coef[colIndex], order, terms, colIndex);
1028 if (!SDDS_SetParameters(&SDDSout, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE, iFitLabelO[colIndex], fitLabelBuffer, -1))
1029 bomb("M", NULL);
1030 }
1031 if (!SDDS_SetParameters(&SDDSout, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE,
1032 iRmsResidualO[colIndex], rmsResidual[colIndex], iChiSqO[colIndex], chi[colIndex],
1033 iTermsO, terms, iSigLevelO[colIndex], ChiSqrSigLevel(chi[colIndex], points - terms),
1034 iOffsetO, xOffset, iFactorO, xScaleFactor, iFitIsValidO[colIndex], isFit[colIndex] ? 'y' : 'n', -1))
1035 bomb("O", NULL);
1036 }
1037 if (!SDDS_WritePage(&SDDSout) || (infoFile && !SDDS_WritePage(&SDDSoutInfo)))
1038 bomb("O", NULL);
1039 }
1040 if (xOrig != x)
1041 free(xOrig);
1042 if (sxOrig != sx)
1043 free(sxOrig);
1044 free(x);
1045 free(sx);
1046 for (colIndex = 0; colIndex < numYNames; colIndex++) {
1047 free(diff[colIndex]);
1048 free(sy[colIndex]);
1049 if (yOrig[colIndex] != y[colIndex])
1050 free(yOrig[colIndex]);
1051 if (syOrig && sy0 && syOrig[colIndex] != sy0[colIndex])
1052 free(syOrig[colIndex]);
1053 free(y[colIndex]);
1054 if (sy0 && sy0[colIndex])
1055 free(sy0[colIndex]);
1056 }
1057 }
1058 return (EXIT_SUCCESS);
1059}
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_SetColumnFromLongs
int32_t SDDS_SetColumnFromLongs(SDDS_DATASET *SDDS_dataset, int32_t mode, int32_t *data, int64_t rows,...)
Sets the values for a single data column using long integer numbers.
Definition SDDS_dataprep.c:2045
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_InitializeInput
int32_t SDDS_InitializeInput(SDDS_DATASET *SDDS_dataset, char *filename)
Definition SDDS_input.c:49
SDDS_ReadPage
int32_t SDDS_ReadPage(SDDS_DATASET *SDDS_dataset)
Definition SDDS_input.c:1006
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_RegisterProgramName
void SDDS_RegisterProgramName(const char *name)
Registers the executable program name for use in error messages.
Definition SDDS_utils.c:288
SDDS_LONG64
#define SDDS_LONG64
Identifier for the signed 64-bit integer data type.
Definition SDDStypes.h:49
bomb
void bomb(char *error, char *usage)
Reports error messages to the terminal and aborts the program.
Definition bomb.c:26
find_min_max
int find_min_max(double *min, double *max, double *list, int64_t n)
Finds the minimum and maximum values in a list of doubles.
Definition findMinMax.c:33
lsfg
long lsfg(double *xd, double *yd, double *sy, long n_pts, long n_terms, int32_t *order, double *coef, double *s_coef, double *chi, double *diff, double(*fn)(double x, long ord))
Computes generalized least squares fits using a function passed by the caller.
Definition lsfg.c:30
match_string
long match_string(char *string, char **option, long n_options, long mode)
Matches a given string against an array of option strings based on specified modes.
Definition match_string.c:34
scanargs
int scanargs(SCANNED_ARG **scanned, int argc, char **argv)
Definition scanargs.c:36
processPipeOption
long processPipeOption(char **item, long items, unsigned long *flags)
Definition scanargs.c:356
processFilenames
void processFilenames(char *programName, char **input, char **output, unsigned long pipeFlags, long noWarnings, long *tmpOutputUsed)
Definition scanargs.c:390
scanItemList
long scanItemList(unsigned long *flags, char **item, long *items, unsigned long mode,...)
Scans a list of items and assigns values based on provided keywords and types.
Definition scanitemlist.c:41
set_argument_scale
void set_argument_scale(double scale)
Set the scale factor applied to the input argument of basis functions.
Definition lsfBasisFns.c:43
dtcheby
double dtcheby(double x, long n)
Evaluate the derivative of the Chebyshev polynomial T_n(x).
Definition lsfBasisFns.c:92
set_argument_offset
void set_argument_offset(double offset)
Set the offset applied to the input argument of basis functions.
Definition lsfBasisFns.c:33
ipower
double ipower(double x, long n)
Evaluate a power function x^n.
Definition lsfBasisFns.c:113
dipower
double dipower(double x, long n)
Evaluate the derivative of x^n.
Definition lsfBasisFns.c:127
tcheby
double tcheby(double x, long n)
Evaluate the Chebyshev polynomial of the first kind T_n(x).
Definition lsfBasisFns.c:72
ChiSqrSigLevel
double ChiSqrSigLevel(double ChiSquared0, long nu)
Computes the probability that a chi-squared variable exceeds a given value.
Definition sigLevel.c:64
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

◆ makeCoefficientUnits()

char ** makeCoefficientUnits ( SDDS_DATASET * SDDSout,
char * xName,
char * yName,
int32_t * order,
long terms )

Definition at line 1521 of file sddsmpfit.c.

1521 {
1522 char *xUnits, *yUnits, buffer[SDDS_MAXLINE];
1523 char **coefUnits = NULL;
1524 long i;
1525
1526 if (!SDDS_GetColumnInformation(SDDSout, "units", &xUnits, SDDS_GET_BY_NAME, xName) ||
1527 !SDDS_GetColumnInformation(SDDSout, "units", &yUnits, SDDS_GET_BY_NAME, yName))
1528 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1529
1530 coefUnits = tmalloc(sizeof(*coefUnits) * terms);
1531 if (!xUnits || SDDS_StringIsBlank(xUnits)) {
1532 if (!yUnits || SDDS_StringIsBlank(yUnits))
1533 SDDS_CopyString(&yUnits, "");
1534 for (i = 0; i < terms; i++)
1535 coefUnits[i] = yUnits;
1536 } else {
1537 if (!yUnits || SDDS_StringIsBlank(yUnits))
1538 SDDS_CopyString(&yUnits, "1");
1539 for (i = 0; i < terms; i++) {
1540 if (order[i] == 0) {
1541 if (strcmp(yUnits, "1") != 0)
1542 SDDS_CopyString(coefUnits + i, yUnits);
1543 else
1544 SDDS_CopyString(coefUnits + i, "");
1545 } else if (strcmp(xUnits, yUnits) == 0) {
1546 if (order[i] > 1)
1547 sprintf(buffer, "1/%s$a%" PRId32 "$n", xUnits, order[i] - 1);
1548 else
1549 strcpy(buffer, "");
1550 SDDS_CopyString(coefUnits + i, buffer);
1551 } else {
1552 if (order[i] > 1)
1553 sprintf(buffer, "%s/%s$a%" PRId32 "$n", yUnits, xUnits, order[i]);
1554 else
1555 sprintf(buffer, "%s/%s", yUnits, xUnits);
1556 SDDS_CopyString(coefUnits + i, buffer);
1557 }
1558 }
1559 }
1560 return coefUnits;
1561}
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

◆ makeEvaluationTable()

void makeEvaluationTable ( EVAL_PARAMETERS * evalParameters,
double * x,
int64_t points,
double * coef,
int32_t * order,
long terms,
char * xName,
char ** yName,
long yNames,
long iYName )

Definition at line 1577 of file sddsmpfit.c.

1578 {
1579 static double *xEval = NULL, *yEval = NULL;
1580 static int64_t maxEvals = 0;
1581 double delta;
1582 int64_t i;
1583
1584 if (!(evalParameters->flags & EVAL_BEGIN_GIVEN) || !(evalParameters->flags & EVAL_END_GIVEN)) {
1585 double min, max;
1586 find_min_max(&min, &max, x, points);
1587 if (!(evalParameters->flags & EVAL_BEGIN_GIVEN))
1588 evalParameters->begin = min;
1589 if (!(evalParameters->flags & EVAL_END_GIVEN))
1590 evalParameters->end = max;
1591 }
1592 if (!(evalParameters->flags & EVAL_NUMBER_GIVEN))
1593 evalParameters->number = points;
1594 if (evalParameters->number > 1)
1595 delta = (evalParameters->end - evalParameters->begin) / (evalParameters->number - 1);
1596 else
1597 delta = 0;
1598
1599 if (!xEval || maxEvals < evalParameters->number) {
1600 if (!(xEval = (double *)SDDS_Realloc(xEval, sizeof(*xEval) * evalParameters->number)) ||
1601 !(yEval = (double *)SDDS_Realloc(yEval, sizeof(*yEval) * evalParameters->number)))
1602 SDDS_Bomb("allocation failure");
1603 maxEvals = evalParameters->number;
1604 }
1605
1606 for (i = 0; i < evalParameters->number; i++) {
1607 xEval[i] = evalParameters->begin + i * delta;
1608 yEval[i] = eval_sum(basis_fn, coef, order, terms, xEval[i]);
1609 }
1610
1611 if ((iYName == 0 &&
1612 !SDDS_StartPage(&evalParameters->dataset, evalParameters->number)) ||
1613 !SDDS_SetColumnFromDoubles(&evalParameters->dataset, SDDS_SET_BY_NAME, xEval, evalParameters->number, xName) ||
1614 !SDDS_SetColumnFromDoubles(&evalParameters->dataset, SDDS_SET_BY_NAME, yEval, evalParameters->number, yName[iYName]) ||
1615 (iYName == yNames - 1 && !SDDS_WritePage(&evalParameters->dataset)))
1616 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1617}
SDDS_Realloc
void * SDDS_Realloc(void *old_ptr, size_t new_size)
Reallocates memory to a new size.
Definition SDDS_utils.c:677

◆ makeFitLabel()

void makeFitLabel ( char * buffer,
long bufsize,
char * fitLabelFormat,
double * coef,
int32_t * order,
long terms,
long colIndex )

Definition at line 1143 of file sddsmpfit.c.

1143 {
1144 long i;
1145 static char buffer1[SDDS_MAXLINE], buffer2[SDDS_MAXLINE];
1146
1147 sprintf(buffer, "%s = ", ySymbols[colIndex]);
1148 for (i = 0; i < terms; i++) {
1149 if (order[i] == 0)
1150 sprintf(buffer1, fitLabelFormat, coef[i]);
1151 else {
1152 if (coef[i] >= 0) {
1153 strcpy(buffer1, " +");
1154 sprintf(buffer1 + 2, fitLabelFormat, coef[i]);
1155 } else
1156 sprintf(buffer1, fitLabelFormat, coef[i]);
1157 strcat(buffer1, "*");
1158 strcat(buffer1, xSymbol);
1159 if (order[i] > 1) {
1160 sprintf(buffer2, "$a%" PRId32 "$n", order[i]);
1161 strcat(buffer1, buffer2);
1162 }
1163 }
1164 if ((long)(strlen(buffer) + strlen(buffer1)) > (long)(0.95 * bufsize)) {
1165 fprintf(stderr, "buffer overflow making fit label!\n");
1166 return;
1167 }
1168 strcat(buffer, buffer1);
1169 }
1170}

◆ print_coefs()

void print_coefs ( FILE * fprec,
double x_offset,
double x_scale,
long chebyshev,
double * coef,
double * s_coef,
int32_t * order,
long n_terms,
double chi,
long norm_term,
char * prepend )

Definition at line 1061 of file sddsmpfit.c.

1062 {
1063 long i;
1064
1065 if (chebyshev)
1066 fprintf(fpo, "%s%ld-term Chebyshev T polynomial least-squares fit about x=%21.15le, scaled by %21.15le:\n", prepend, terms, xOffset, xScaleFactor);
1067 else
1068 fprintf(fpo, "%s%ld-term polynomial least-squares fit about x=%21.15le:\n", prepend, terms, xOffset);
1069 if (normTerm >= 0 && terms > normTerm) {
1070 if (coef[normTerm] != 0)
1071 fprintf(fpo, "%s coefficients are normalized with factor %21.15le to make a[%ld]==1\n", prepend, coef[normTerm], (order ? order[normTerm] : normTerm));
1072 else {
1073 fprintf(fpo, "%s can't normalize coefficients as requested: a[%ld]==0\n", prepend, (order ? order[normTerm] : normTerm));
1074 normTerm = -1;
1075 }
1076 } else
1077 normTerm = -1;
1078
1079 for (i = 0; i < terms; i++) {
1080 fprintf(fpo, "%sa[%ld] = %21.15le ", prepend, (order ? order[i] : i), (normTerm < 0 ? coef[i] : coef[i] / coef[normTerm]));
1081 if (coefSigma)
1082 fprintf(fpo, "+/- %21.15le\n", (normTerm < 0 ? coefSigma[i] : coefSigma[i] / fabs(coef[normTerm])));
1083 else
1084 fputc('\n', fpo);
1085 }
1086 if (coefSigma)
1087 fprintf(fpo, "%sreduced chi-squared = %21.15le\n", prepend, chi);
1088}

◆ RemoveElementFromStringArray()

void RemoveElementFromStringArray ( char ** array,
long index,
long length )

Definition at line 1090 of file sddsmpfit.c.

1090 {
1091 long lh;
1092
1093 for (lh = index; lh < length - 1; lh++)
1094 array[lh] = array[lh + 1];
1095}

◆ RemoveNonNumericColumnsFromNameArray()

char ** RemoveNonNumericColumnsFromNameArray ( SDDS_DATASET * SDDSin,
char ** columns,
int32_t * numColumns )

Definition at line 1097 of file sddsmpfit.c.

1097 {
1098 long i, numNumericColumns = *numColumns;
1099
1100 for (i = 0; i < *numColumns; i++) {
1101 if (SDDS_CheckColumn(SDDSin, columns[i], NULL, SDDS_ANY_NUMERIC_TYPE, NULL)) {
1102 printf("Removing %s because not a numeric type.\n", columns[i]);
1103 RemoveElementFromStringArray(columns, i, *numColumns);
1104 numNumericColumns--;
1105 }
1106 }
1107
1108 *numColumns = numNumericColumns;
1109 return (columns);
1110}
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_ANY_NUMERIC_TYPE
#define SDDS_ANY_NUMERIC_TYPE
Special identifier used by SDDS_Check*() routines to accept any numeric type.
Definition SDDStypes.h:157

◆ ResolveColumnNames()

char ** ResolveColumnNames ( SDDS_DATASET * SDDSin,
char ** wildcardList,
long length,
int32_t * numYNames )

Definition at line 1112 of file sddsmpfit.c.

1112 {
1113 char **result;
1114 long i;
1115
1116 /* initially set the columns of interest to none, to make SDDS_OR work below */
1117 SDDS_SetColumnsOfInterest(SDDSin, SDDS_MATCH_STRING, "", SDDS_AND);
1118 for (i = 0; i < length; i++) {
1119 SDDS_SetColumnsOfInterest(SDDSin, SDDS_MATCH_STRING, wildcardList[i], SDDS_OR);
1120 }
1121
1122 if (!(result = SDDS_GetColumnNames(SDDSin, numYNames)) || *numYNames == 0)
1123 bomb("Error matching columns in ResolveColumnNames: No matches.", NULL);
1124
1125 result = RemoveNonNumericColumnsFromNameArray(SDDSin, result, numYNames);
1126 return (result);
1127}
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_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

◆ rms_average()

double rms_average ( double * x,
int64_t n )

Definition at line 1172 of file sddsmpfit.c.

1172 {
1173 double sum2;
1174 int64_t i;
1175
1176 for (i = sum2 = 0; i < n; i++)
1177 sum2 += sqr(x[i]);
1178
1179 return (sqrt(sum2 / n));
1180}

◆ set_argument_offset()

void set_argument_offset ( double offset)

Set the offset applied to the input argument of basis functions.

This function updates the global variable used to shift the input argument before evaluating polynomial or Chebyshev functions.

Parameters
offsetThe new offset to apply to the argument.

Definition at line 33 of file lsfBasisFns.c.

33{ x_offset = offset; }

◆ set_argument_scale()

void set_argument_scale ( double scale)

Set the scale factor applied to the input argument of basis functions.

This function updates the global variable used to scale the input argument before evaluating polynomial or Chebyshev functions. It ensures the scale factor is not zero.

Parameters
scaleThe new scale factor for the argument.

Definition at line 43 of file lsfBasisFns.c.

43 {
44 if (!(x_scale = scale))
45 bomb("argument scale factor is zero", NULL);
46}

◆ setupEvaluationFile()

void setupEvaluationFile ( EVAL_PARAMETERS * evalParameters,
char * xName,
char ** yName,
long yNames,
SDDS_DATASET * SDDSin )

Definition at line 1563 of file sddsmpfit.c.

1563 {
1564 long i;
1565 SDDS_DATASET *SDDSout;
1566 SDDSout = &evalParameters->dataset;
1567 if (!SDDS_InitializeOutput(SDDSout, SDDS_BINARY, 0, NULL, "sddsmpfit output: evaluation of fits", evalParameters->file) ||
1568 !SDDS_TransferColumnDefinition(SDDSout, SDDSin, xName, NULL))
1569 SDDS_Bomb("Problem setting up evaluation file");
1570 for (i = 0; i < yNames; i++)
1571 if (!SDDS_TransferColumnDefinition(SDDSout, SDDSin, yName[i], NULL))
1572 SDDS_Bomb("Problem setting up evaluation file");
1573 if (!SDDS_WriteLayout(SDDSout))
1574 SDDS_Bomb("Problem setting up evaluation file");
1575}

◆ tcheby()

double tcheby ( double x,
long n )

Evaluate the Chebyshev polynomial of the first kind T_n(x).

Given x and an order n, this function returns T_n((x - x_offset) / x_scale). If x is out of the domain [-1,1], it is clipped to ±1 before evaluation.

Parameters
xThe point at which to evaluate the Chebyshev polynomial.
nThe order of the Chebyshev polynomial.
Returns
The value of T_n(x).

Definition at line 72 of file lsfBasisFns.c.

72 {
73 x = (x - x_offset) / x_scale;
74 if (x > 1 || x < -1) {
75 /* fprintf(stderr, "warning: argument %e is out of range for tcheby()\n",
76 * x); */
77 x = SIGN(x);
78 }
79 return (cos(n * acos(x)));
80}
Generated by doxygen 1.12.0