SDDSlib
Loading...
Searching...
No Matches
Classes | Macros | Enumerations | Functions | Variables
sddsmpfit.c File Reference

Performs nth-order polynomial least squares fitting for SDDS files. More...

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

Go to the source code of this file.

Classes

struct  EVAL_PARAMETERS
 

Macros

#define NO_SYMMETRY   0
 
#define EVEN_SYMMETRY   1
 
#define ODD_SYMMETRY   2
 
#define N_SYMMETRY_OPTIONS   3
 
#define ABSOLUTE_SIGMAS   0
 
#define FRACTIONAL_SIGMAS   1
 
#define N_SIGMAS_OPTIONS   2
 
#define FLGS_GENERATESIGMAS   1
 
#define FLGS_KEEPLARGEST   2
 
#define FLGS_KEEPSMALLEST   4
 
#define REVPOW_ACTIVE   0x0001
 
#define REVPOW_VERBOSE   0x0002
 
#define EVAL_BEGIN_GIVEN   0x0001U
 
#define EVAL_END_GIVEN   0x0002U
 
#define EVAL_NUMBER_GIVEN   0x0004U
 
#define MAX_Y_SIGMA_NAME_SIZE   1024
 

Enumerations

enum  option_type {
  CLO_DEPENDENT , CLO_ORDERS , CLO_TERMS , CLO_SYMMETRY ,
  CLO_REVISEORDERS , CLO_CHEBYSHEV , CLO_MODIFYSIGMAS , CLO_SIGMAS ,
  CLO_GENERATESIGMAS , CLO_RANGE , CLO_SPARSE , CLO_NORMALIZE ,
  CLO_XFACTOR , CLO_XOFFSET , CLO_VERBOSE , CLO_FITLABELFORMAT ,
  CLO_PIPE , CLO_EVALUATE , CLO_INDEPENDENT , CLO_SIGMAINDEPENDENT ,
  CLO_SIGMADEPENDENT , CLO_INFOFILE , CLO_COPYPARAMETERS , CLO_MINSIGMA ,
  N_OPTIONS
}
 

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)
 

Variables

char * option [N_OPTIONS]
 
char * USAGE
 
static char * additional_help
 
static char * additional_help2
 
char * symmetry_options [N_SYMMETRY_OPTIONS] = {"none", "even", "odd"}
 
char * sigmas_options [N_SIGMAS_OPTIONS] = {"absolute", "fractional"}
 
static long * iIntercept = NULL
 
static long * iInterceptO = NULL
 
static long * iInterceptSigma = NULL
 
static long * iInterceptSigmaO = NULL
 
static long * iSlope = NULL
 
static long * iSlopeO = NULL
 
static long * iSlopeSigma = NULL
 
static long * iSlopeSigmaO = NULL
 
static long * iCurvature = NULL
 
static long * iCurvatureO = NULL
 
static long * iCurvatureSigma = NULL
 
static long * iCurvatureSigmaO = NULL
 
static long iOffset = -1
 
static long iOffsetO = -1
 
static long iFactor = -1
 
static long iFactorO = -1
 
static long * iChiSq = NULL
 
static long * iChiSqO = NULL
 
static long * iRmsResidual = NULL
 
static long * iRmsResidualO = NULL
 
static long * iSigLevel = NULL
 
static long * iSigLevelO = NULL
 
static long * iFitIsValid = NULL
 
static long * iFitIsValidO = NULL
 
static long * iFitLabel = NULL
 
static long * iFitLabelO = NULL
 
static long iTerms = -1
 
static long iTermsO = -1
 
static long ix = -1
 
static long ixSigma = -1
 
static long * iy = NULL
 
static long * iySigma = NULL
 
static long * iFit = NULL
 
static long * iResidual = NULL
 
static long iOrder = -1
 
static long * iCoefficient = NULL
 
static long * iCoefficientSigma = NULL
 
static long * iCoefficientUnits = NULL
 
static char * xSymbol
 
static char ** ySymbols
 
static double(* basis_fn )(double xa, long ordera)
 
static double(* basis_dfn )(double xa, long ordera)
 

Detailed Description

Performs nth-order polynomial least squares fitting for SDDS files.

This program fits data from SDDS (Self Describing Data Set) files using either ordinary polynomials or Chebyshev T polynomials. It outputs the fitted data, residuals, and fit parameters with their uncertainties. Additional functionalities include sigma modification, sigma generation, order revision, and fit evaluation.

The fitting model is of the form:

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

where ( P(x, i) ) is the ith basis function evaluated at ( x ). By default, ( P(x, i) = x^i ), but Chebyshev T polynomials can also be selected as the basis functions.

The program reads input SDDS files, performs the fitting, and writes the output to specified SDDS files. Various options allow customization of the fitting process, including normalization, range selection, sigma handling, and more.

@usage

sddsmpfit [-pipe=[input][,output]] [<inputfile>] [<outputfile>]
-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]
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.

Macro Definition Documentation

◆ ABSOLUTE_SIGMAS

#define ABSOLUTE_SIGMAS   0

Definition at line 213 of file sddsmpfit.c.

◆ EVAL_BEGIN_GIVEN

#define EVAL_BEGIN_GIVEN   0x0001U

Definition at line 240 of file sddsmpfit.c.

◆ EVAL_END_GIVEN

#define EVAL_END_GIVEN   0x0002U

Definition at line 241 of file sddsmpfit.c.

◆ EVAL_NUMBER_GIVEN

#define EVAL_NUMBER_GIVEN   0x0004U

Definition at line 242 of file sddsmpfit.c.

◆ EVEN_SYMMETRY

#define EVEN_SYMMETRY   1

Definition at line 208 of file sddsmpfit.c.

◆ FLGS_GENERATESIGMAS

#define FLGS_GENERATESIGMAS   1

Definition at line 218 of file sddsmpfit.c.

◆ FLGS_KEEPLARGEST

#define FLGS_KEEPLARGEST   2

Definition at line 219 of file sddsmpfit.c.

◆ FLGS_KEEPSMALLEST

#define FLGS_KEEPSMALLEST   4

Definition at line 220 of file sddsmpfit.c.

◆ FRACTIONAL_SIGMAS

#define FRACTIONAL_SIGMAS   1

Definition at line 214 of file sddsmpfit.c.

◆ MAX_Y_SIGMA_NAME_SIZE

#define MAX_Y_SIGMA_NAME_SIZE   1024

Definition at line 244 of file sddsmpfit.c.

◆ N_SIGMAS_OPTIONS

#define N_SIGMAS_OPTIONS   2

Definition at line 215 of file sddsmpfit.c.

◆ N_SYMMETRY_OPTIONS

#define N_SYMMETRY_OPTIONS   3

Definition at line 210 of file sddsmpfit.c.

◆ NO_SYMMETRY

#define NO_SYMMETRY   0

Definition at line 207 of file sddsmpfit.c.

◆ ODD_SYMMETRY

#define ODD_SYMMETRY   2

Definition at line 209 of file sddsmpfit.c.

◆ REVPOW_ACTIVE

#define REVPOW_ACTIVE   0x0001

Definition at line 222 of file sddsmpfit.c.

◆ REVPOW_VERBOSE

#define REVPOW_VERBOSE   0x0002

Definition at line 223 of file sddsmpfit.c.

Enumeration Type Documentation

◆ option_type

enum option_type

Definition at line 82 of file sddsmpfit.c.

82 {
83 CLO_DEPENDENT,
84 CLO_ORDERS,
85 CLO_TERMS,
86 CLO_SYMMETRY,
87 CLO_REVISEORDERS,
88 CLO_CHEBYSHEV,
89 CLO_MODIFYSIGMAS,
90 CLO_SIGMAS,
91 CLO_GENERATESIGMAS,
92 CLO_RANGE,
93 CLO_SPARSE,
94 CLO_NORMALIZE,
95 CLO_XFACTOR,
96 CLO_XOFFSET,
97 CLO_VERBOSE,
98 CLO_FITLABELFORMAT,
99 CLO_PIPE,
100 CLO_EVALUATE,
101 CLO_INDEPENDENT,
102 CLO_SIGMAINDEPENDENT,
103 CLO_SIGMADEPENDENT,
104 CLO_INFOFILE,
105 CLO_COPYPARAMETERS,
106 CLO_MINSIGMA,
107 N_OPTIONS
108};

Function Documentation

◆ checkInputFile()

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

Definition at line 1151 of file sddsmpfit.c.

1151 {
1152 char *ptr = NULL;
1153 long i;
1154
1155 if (!(ptr = SDDS_FindColumn(SDDSin, FIND_NUMERIC_TYPE, xName, NULL)))
1156 SDDS_Bomb("x column doesn't exist or is nonnumeric");
1157 free(ptr);
1158
1159 /* y columns don't need to be checked because located using SDDS_SetColumnsOfInterest */
1160
1161 ptr = NULL;
1162 if (xSigmaName && !(ptr = SDDS_FindColumn(SDDSin, FIND_NUMERIC_TYPE, xSigmaName, NULL)))
1163 SDDS_Bomb("x sigma column doesn't exist or is nonnumeric");
1164 if (ptr)
1165 free(ptr);
1166
1167 if (ySigmaNames) {
1168 for (i = 0; i < numYNames; i++) {
1169 ptr = NULL;
1170 if (!(ptr = SDDS_FindColumn(SDDSin, FIND_NUMERIC_TYPE, ySigmaNames[i], NULL)))
1171 SDDS_Bomb("y sigma column doesn't exist or is nonnumeric");
1172 if (ptr)
1173 free(ptr);
1174 }
1175 }
1176}
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 1143 of file sddsmpfit.c.

1143 {
1144 long i;
1145 for (i = 0; i < terms; i++)
1146 if (order[i] == order_of_interest)
1147 return (i);
1148 return (-1);
1149}

◆ compareOriginalToFit()

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

Definition at line 1580 of file sddsmpfit.c.

1580 {
1581 int64_t i;
1582 double residualSum2, fit;
1583
1584 *residual = tmalloc(sizeof(**residual) * points);
1585
1586 residualSum2 = 0;
1587 for (i = 0; i < points; i++) {
1588 fit = eval_sum(basis_fn, coef, order, terms, x[i]);
1589 (*residual)[i] = y[i] - fit;
1590 residualSum2 += sqr((*residual)[i]);
1591 }
1592 *rmsResidual = sqrt(residualSum2 / points);
1593}
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 1090 of file sddsmpfit.c.

1090 {
1091 long i, nameLength;
1092 char **result, sigmaName[MAX_Y_SIGMA_NAME_SIZE];
1093
1094 result = tmalloc(sizeof(char *) * numYNames);
1095 for (i = 0; i < numYNames; i++) {
1096 sprintf(sigmaName, controlString, yNames[i]);
1097 nameLength = strlen(sigmaName);
1098 result[i] = tmalloc(sizeof(char) * (nameLength + 1));
1099 strcpy(result[i], sigmaName);
1100 }
1101 return (result);
1102}

◆ 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 1178 of file sddsmpfit.c.

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

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

◆ makeCoefficientUnits()

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

Definition at line 1482 of file sddsmpfit.c.

1482 {
1483 char *xUnits, *yUnits, buffer[SDDS_MAXLINE];
1484 char **coefUnits = NULL;
1485 long i;
1486
1487 if (!SDDS_GetColumnInformation(SDDSout, "units", &xUnits, SDDS_GET_BY_NAME, xName) ||
1488 !SDDS_GetColumnInformation(SDDSout, "units", &yUnits, SDDS_GET_BY_NAME, yName))
1489 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1490
1491 coefUnits = tmalloc(sizeof(*coefUnits) * terms);
1492 if (!xUnits || SDDS_StringIsBlank(xUnits)) {
1493 if (!yUnits || SDDS_StringIsBlank(yUnits))
1494 SDDS_CopyString(&yUnits, "");
1495 for (i = 0; i < terms; i++)
1496 coefUnits[i] = yUnits;
1497 } else {
1498 if (!yUnits || SDDS_StringIsBlank(yUnits))
1499 SDDS_CopyString(&yUnits, "1");
1500 for (i = 0; i < terms; i++) {
1501 if (order[i] == 0) {
1502 if (strcmp(yUnits, "1") != 0)
1503 SDDS_CopyString(coefUnits + i, yUnits);
1504 else
1505 SDDS_CopyString(coefUnits + i, "");
1506 } else if (strcmp(xUnits, yUnits) == 0) {
1507 if (order[i] > 1)
1508 sprintf(buffer, "1/%s$a%" PRId32 "$n", xUnits, order[i] - 1);
1509 else
1510 strcpy(buffer, "");
1511 SDDS_CopyString(coefUnits + i, buffer);
1512 } else {
1513 if (order[i] > 1)
1514 sprintf(buffer, "%s/%s$a%" PRId32 "$n", yUnits, xUnits, order[i]);
1515 else
1516 sprintf(buffer, "%s/%s", yUnits, xUnits);
1517 SDDS_CopyString(coefUnits + i, buffer);
1518 }
1519 }
1520 }
1521 return coefUnits;
1522}
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 1538 of file sddsmpfit.c.

1539 {
1540 static double *xEval = NULL, *yEval = NULL;
1541 static int64_t maxEvals = 0;
1542 double delta;
1543 int64_t i;
1544
1545 if (!(evalParameters->flags & EVAL_BEGIN_GIVEN) || !(evalParameters->flags & EVAL_END_GIVEN)) {
1546 double min, max;
1547 find_min_max(&min, &max, x, points);
1548 if (!(evalParameters->flags & EVAL_BEGIN_GIVEN))
1549 evalParameters->begin = min;
1550 if (!(evalParameters->flags & EVAL_END_GIVEN))
1551 evalParameters->end = max;
1552 }
1553 if (!(evalParameters->flags & EVAL_NUMBER_GIVEN))
1554 evalParameters->number = points;
1555 if (evalParameters->number > 1)
1556 delta = (evalParameters->end - evalParameters->begin) / (evalParameters->number - 1);
1557 else
1558 delta = 0;
1559
1560 if (!xEval || maxEvals < evalParameters->number) {
1561 if (!(xEval = (double *)SDDS_Realloc(xEval, sizeof(*xEval) * evalParameters->number)) ||
1562 !(yEval = (double *)SDDS_Realloc(yEval, sizeof(*yEval) * evalParameters->number)))
1563 SDDS_Bomb("allocation failure");
1564 maxEvals = evalParameters->number;
1565 }
1566
1567 for (i = 0; i < evalParameters->number; i++) {
1568 xEval[i] = evalParameters->begin + i * delta;
1569 yEval[i] = eval_sum(basis_fn, coef, order, terms, xEval[i]);
1570 }
1571
1572 if ((iYName == 0 &&
1573 !SDDS_StartPage(&evalParameters->dataset, evalParameters->number)) ||
1574 !SDDS_SetColumnFromDoubles(&evalParameters->dataset, SDDS_SET_BY_NAME, xEval, evalParameters->number, xName) ||
1575 !SDDS_SetColumnFromDoubles(&evalParameters->dataset, SDDS_SET_BY_NAME, yEval, evalParameters->number, yName[iYName]) ||
1576 (iYName == yNames - 1 && !SDDS_WritePage(&evalParameters->dataset)))
1577 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1578}
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 1104 of file sddsmpfit.c.

1104 {
1105 long i;
1106 static char buffer1[SDDS_MAXLINE], buffer2[SDDS_MAXLINE];
1107
1108 sprintf(buffer, "%s = ", ySymbols[colIndex]);
1109 for (i = 0; i < terms; i++) {
1110 if (order[i] == 0)
1111 sprintf(buffer1, fitLabelFormat, coef[i]);
1112 else {
1113 if (coef[i] >= 0) {
1114 strcpy(buffer1, " +");
1115 sprintf(buffer1 + 2, fitLabelFormat, coef[i]);
1116 } else
1117 sprintf(buffer1, fitLabelFormat, coef[i]);
1118 strcat(buffer1, "*");
1119 strcat(buffer1, xSymbol);
1120 if (order[i] > 1) {
1121 sprintf(buffer2, "$a%" PRId32 "$n", order[i]);
1122 strcat(buffer1, buffer2);
1123 }
1124 }
1125 if ((long)(strlen(buffer) + strlen(buffer1)) > (long)(0.95 * bufsize)) {
1126 fprintf(stderr, "buffer overflow making fit label!\n");
1127 return;
1128 }
1129 strcat(buffer, buffer1);
1130 }
1131}

◆ 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 1022 of file sddsmpfit.c.

1023 {
1024 long i;
1025
1026 if (chebyshev)
1027 fprintf(fpo, "%s%ld-term Chebyshev T polynomial least-squares fit about x=%21.15le, scaled by %21.15le:\n", prepend, terms, xOffset, xScaleFactor);
1028 else
1029 fprintf(fpo, "%s%ld-term polynomial least-squares fit about x=%21.15le:\n", prepend, terms, xOffset);
1030 if (normTerm >= 0 && terms > normTerm) {
1031 if (coef[normTerm] != 0)
1032 fprintf(fpo, "%s coefficients are normalized with factor %21.15le to make a[%ld]==1\n", prepend, coef[normTerm], (order ? order[normTerm] : normTerm));
1033 else {
1034 fprintf(fpo, "%s can't normalize coefficients as requested: a[%ld]==0\n", prepend, (order ? order[normTerm] : normTerm));
1035 normTerm = -1;
1036 }
1037 } else
1038 normTerm = -1;
1039
1040 for (i = 0; i < terms; i++) {
1041 fprintf(fpo, "%sa[%ld] = %21.15le ", prepend, (order ? order[i] : i), (normTerm < 0 ? coef[i] : coef[i] / coef[normTerm]));
1042 if (coefSigma)
1043 fprintf(fpo, "+/- %21.15le\n", (normTerm < 0 ? coefSigma[i] : coefSigma[i] / fabs(coef[normTerm])));
1044 else
1045 fputc('\n', fpo);
1046 }
1047 if (coefSigma)
1048 fprintf(fpo, "%sreduced chi-squared = %21.15le\n", prepend, chi);
1049}

◆ RemoveElementFromStringArray()

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

Definition at line 1051 of file sddsmpfit.c.

1051 {
1052 long lh;
1053
1054 for (lh = index; lh < length - 1; lh++)
1055 array[lh] = array[lh + 1];
1056}

◆ RemoveNonNumericColumnsFromNameArray()

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

Definition at line 1058 of file sddsmpfit.c.

1058 {
1059 long i, numNumericColumns = *numColumns;
1060
1061 for (i = 0; i < *numColumns; i++) {
1062 if (SDDS_CheckColumn(SDDSin, columns[i], NULL, SDDS_ANY_NUMERIC_TYPE, NULL)) {
1063 printf("Removing %s because not a numeric type.\n", columns[i]);
1064 RemoveElementFromStringArray(columns, i, *numColumns);
1065 numNumericColumns--;
1066 }
1067 }
1068
1069 *numColumns = numNumericColumns;
1070 return (columns);
1071}
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 1073 of file sddsmpfit.c.

1073 {
1074 char **result;
1075 long i;
1076
1077 /* initially set the columns of interest to none, to make SDDS_OR work below */
1078 SDDS_SetColumnsOfInterest(SDDSin, SDDS_MATCH_STRING, "", SDDS_AND);
1079 for (i = 0; i < length; i++) {
1080 SDDS_SetColumnsOfInterest(SDDSin, SDDS_MATCH_STRING, wildcardList[i], SDDS_OR);
1081 }
1082
1083 if (!(result = SDDS_GetColumnNames(SDDSin, numYNames)) || *numYNames == 0)
1084 bomb("Error matching columns in ResolveColumnNames: No matches.", NULL);
1085
1086 result = RemoveNonNumericColumnsFromNameArray(SDDSin, result, numYNames);
1087 return (result);
1088}
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 1133 of file sddsmpfit.c.

1133 {
1134 double sum2;
1135 int64_t i;
1136
1137 for (i = sum2 = 0; i < n; i++)
1138 sum2 += sqr(x[i]);
1139
1140 return (sqrt(sum2 / n));
1141}

◆ 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 1524 of file sddsmpfit.c.

1524 {
1525 long i;
1526 SDDS_DATASET *SDDSout;
1527 SDDSout = &evalParameters->dataset;
1528 if (!SDDS_InitializeOutput(SDDSout, SDDS_BINARY, 0, NULL, "sddsmpfit output: evaluation of fits", evalParameters->file) ||
1529 !SDDS_TransferColumnDefinition(SDDSout, SDDSin, xName, NULL))
1530 SDDS_Bomb("Problem setting up evaluation file");
1531 for (i = 0; i < yNames; i++)
1532 if (!SDDS_TransferColumnDefinition(SDDSout, SDDSin, yName[i], NULL))
1533 SDDS_Bomb("Problem setting up evaluation file");
1534 if (!SDDS_WriteLayout(SDDSout))
1535 SDDS_Bomb("Problem setting up evaluation file");
1536}

◆ 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}

Variable Documentation

◆ additional_help

char* additional_help
static
Initial value:
= "\n\
sddsmpfit does fits to the form y = SUM(i){ A[i] *P(x-x_offset, i)}, where P(x,i) is the ith basis\n\
function evaluated at x. sddsmpfit returns the A[i] and estimates of the errors in the values.\n\
By default P(x,i) = x^i. One can also select Chebyshev T polynomials as the basis functions.\n\n\
-independent specify name of independent data column to use.\n\
-dependent specify names of dependent data columns to use, using wildcards,\n\
separated by commas.\n\
-sigmaIndependent specify name of independent sigma values to use\n\
-sigmaDependent specify names of dependent sigma values to use by specifying a printf-style control\n\
string to generate the names from the independent variable names (e.g., %sSigma)\n\
-terms number of terms desired in fit.\n\
-symmetry symmetry of desired fit about x_offset.\n\
-orders orders (P[i]) to use in fitting.\n\
-reviseOrders the orders used in the fit are modified from the specified ones\n\
in order eliminate poorly-determined coefficients, based on fitting\n\
of the first data page.\n"

Definition at line 169 of file sddsmpfit.c.

◆ additional_help2

char* additional_help2
static
Initial value:
= "-chebyshev use Chebyshev T polynomials (xOffset is set automatically).\n\
Giving the `convert' option causes the fit to be written out in\n\
terms of ordinary polynomials.\n\
-xOffset desired value of x to fit about.\n\
-xFactor desired factor to multiply x values by before fitting.\n\
-sigmas specify absolute or fractional sigma for all points.\n\
-minimumSigma specify minimum sigma value. If the value is less than this\n\
it is replaced by this value.\n\
-modifySigmas modify the y sigmas using the x sigmas and an initial fit.\n\
-generateSigmas generate y sigmas from the rms deviation from an initial fit.\n\
optionally keep the sigmas from the data if larger/smaller than rms\n\
deviation.\n\
-sparse specify integer interval at which to sample data.\n\
-range specify range of independent variable over which to perform fit and evaluation.\n\
If 'fitOnly' is given, then fit is compared to data over the original range.\n\
-normalize normalize so that specified term is unity.\n\
-evaluate specify evaluation of fit over a selected range of\n\
equispaced points.\n\
-infoFile specify name of optional information file containing coefficients and fit statistics.\n\
-copyParameters specify that parameters from input should be copied to output.\n\
-verbose generates extra output that may be useful.\n\n"

Definition at line 185 of file sddsmpfit.c.

◆ basis_dfn

double(* basis_dfn) (double xa, long ordera) ( double xa,
long ordera )
static

Definition at line 260 of file sddsmpfit.c.

◆ basis_fn

double(* basis_fn) (double xa, long ordera) ( double xa,
long ordera )
static

Definition at line 259 of file sddsmpfit.c.

◆ iChiSq

long* iChiSq = NULL
static

Definition at line 229 of file sddsmpfit.c.

◆ iChiSqO

long * iChiSqO = NULL
static

Definition at line 229 of file sddsmpfit.c.

◆ iCoefficient

long * iCoefficient = NULL
static

Definition at line 236 of file sddsmpfit.c.

◆ iCoefficientSigma

long * iCoefficientSigma = NULL
static

Definition at line 236 of file sddsmpfit.c.

◆ iCoefficientUnits

long * iCoefficientUnits = NULL
static

Definition at line 236 of file sddsmpfit.c.

◆ iCurvature

long* iCurvature = NULL
static

Definition at line 227 of file sddsmpfit.c.

◆ iCurvatureO

long * iCurvatureO = NULL
static

Definition at line 227 of file sddsmpfit.c.

◆ iCurvatureSigma

long * iCurvatureSigma = NULL
static

Definition at line 227 of file sddsmpfit.c.

◆ iCurvatureSigmaO

long * iCurvatureSigmaO = NULL
static

Definition at line 227 of file sddsmpfit.c.

◆ iFactor

long iFactor = -1
static

Definition at line 228 of file sddsmpfit.c.

◆ iFactorO

long iFactorO = -1
static

Definition at line 228 of file sddsmpfit.c.

◆ iFit

long* iFit = NULL
static

Definition at line 234 of file sddsmpfit.c.

◆ iFitIsValid

long* iFitIsValid = NULL
static

Definition at line 230 of file sddsmpfit.c.

◆ iFitIsValidO

long * iFitIsValidO = NULL
static

Definition at line 230 of file sddsmpfit.c.

◆ iFitLabel

long * iFitLabel = NULL
static

Definition at line 230 of file sddsmpfit.c.

◆ iFitLabelO

long * iFitLabelO = NULL
static

Definition at line 230 of file sddsmpfit.c.

◆ iIntercept

long* iIntercept = NULL
static

Definition at line 225 of file sddsmpfit.c.

◆ iInterceptO

long * iInterceptO = NULL
static

Definition at line 225 of file sddsmpfit.c.

◆ iInterceptSigma

long * iInterceptSigma = NULL
static

Definition at line 225 of file sddsmpfit.c.

◆ iInterceptSigmaO

long * iInterceptSigmaO = NULL
static

Definition at line 225 of file sddsmpfit.c.

◆ iOffset

long iOffset = -1
static

Definition at line 228 of file sddsmpfit.c.

◆ iOffsetO

long iOffsetO = -1
static

Definition at line 228 of file sddsmpfit.c.

◆ iOrder

long iOrder = -1
static

Definition at line 236 of file sddsmpfit.c.

◆ iResidual

long * iResidual = NULL
static

Definition at line 234 of file sddsmpfit.c.

◆ iRmsResidual

long * iRmsResidual = NULL
static

Definition at line 229 of file sddsmpfit.c.

◆ iRmsResidualO

long * iRmsResidualO = NULL
static

Definition at line 229 of file sddsmpfit.c.

◆ iSigLevel

long * iSigLevel = NULL
static

Definition at line 229 of file sddsmpfit.c.

◆ iSigLevelO

long * iSigLevelO = NULL
static

Definition at line 229 of file sddsmpfit.c.

◆ iSlope

long* iSlope = NULL
static

Definition at line 226 of file sddsmpfit.c.

◆ iSlopeO

long * iSlopeO = NULL
static

Definition at line 226 of file sddsmpfit.c.

◆ iSlopeSigma

long * iSlopeSigma = NULL
static

Definition at line 226 of file sddsmpfit.c.

◆ iSlopeSigmaO

long * iSlopeSigmaO = NULL
static

Definition at line 226 of file sddsmpfit.c.

◆ iTerms

long iTerms = -1
static

Definition at line 230 of file sddsmpfit.c.

◆ iTermsO

long iTermsO = -1
static

Definition at line 230 of file sddsmpfit.c.

◆ ix

long ix = -1
static

Definition at line 232 of file sddsmpfit.c.

◆ ixSigma

long ixSigma = -1
static

Definition at line 232 of file sddsmpfit.c.

◆ iy

long* iy = NULL
static

Definition at line 233 of file sddsmpfit.c.

◆ iySigma

long * iySigma = NULL
static

Definition at line 233 of file sddsmpfit.c.

◆ option

char* option[N_OPTIONS]
Initial value:
= {
"dependent",
"orders",
"terms",
"symmetry",
"reviseorders",
"chebyshev",
"modifysigmas",
"sigmas",
"generatesigmas",
"range",
"sparse",
"normalize",
"xfactor",
"xoffset",
"verbose",
"fitlabelformat",
"pipe",
"evaluate",
"independent",
"sigmaindependent",
"sigmadependent",
"infofile",
"copyparameters",
"minimumsigma",
}

Definition at line 110 of file sddsmpfit.c.

110 {
111 "dependent",
112 "orders",
113 "terms",
114 "symmetry",
115 "reviseorders",
116 "chebyshev",
117 "modifysigmas",
118 "sigmas",
119 "generatesigmas",
120 "range",
121 "sparse",
122 "normalize",
123 "xfactor",
124 "xoffset",
125 "verbose",
126 "fitlabelformat",
127 "pipe",
128 "evaluate",
129 "independent",
130 "sigmaindependent",
131 "sigmadependent",
132 "infofile",
133 "copyparameters",
134 "minimumsigma",
135};

◆ sigmas_options

char* sigmas_options[N_SIGMAS_OPTIONS] = {"absolute", "fractional"}

Definition at line 216 of file sddsmpfit.c.

216{"absolute", "fractional"};

◆ symmetry_options

char* symmetry_options[N_SYMMETRY_OPTIONS] = {"none", "even", "odd"}

Definition at line 211 of file sddsmpfit.c.

211{"none", "even", "odd"};

◆ USAGE

char* USAGE

Definition at line 137 of file sddsmpfit.c.

◆ xSymbol

char* xSymbol
static

Definition at line 238 of file sddsmpfit.c.

◆ ySymbols

char ** ySymbols
static

Definition at line 238 of file sddsmpfit.c.

Generated by doxygen 1.12.0