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

Detailed Description

Performs interpolation on SDDS formatted data.

The program reads SDDS data files, performs interpolation based on user-specified options, and writes the results back in SDDS format. It supports various interpolation methods, handles monotonicity enforcement, and manages out-of-range conditions.

Usage

sddsinterp [<inputfile>] [<outputfile>]
[-pipe=[input][,output]]
-columns=<independent-quantity>,<dependent-name>[,...]
[-exclude=<name>[,...]]
{
-atValues=<values-list> |
-sequence=<points>[,<start>,<end>] |
-equispaced=<spacing>[,<start>,<end>] |
-fillIn |
-fileValues=<valuesfile>[,column=<column-name>][,parallelPages]]
[-interpShort=-1|-2]
[-order=<number>]
[-printout[=bare][,stdout]]
[-forceMonotonic[={increasing|decreasing}]]
[-belowRange={value=<value>|skip|saturate|extrapolate|wrap}[,{abort|warn}]]
[-aboveRange={value=<value>|skip|saturate|extrapolate|wrap}[,{abort|warn}]]
[-majorOrder=row|column]

Options

Required Description
-columns Specify the independent and dependent columns.
Option Description
-pipe Use pipe for input and/or output.
-exclude Exclude specified columns from processing.
-atValues Interpolate at the specified list of values.
-sequence Generate a sequence of interpolation points.
-equispaced Generate equispaced interpolation points.
-fillIn Automatically fill in interpolation points based on data.
-fileValues Use values from a file for interpolation.
-interpShort Interpolate short columns with specified order (-1 or -2).
-order Set the interpolation order (default is 1).
-printout Output interpolated data in bare format and/or to stdout.
-forceMonotonic Enforce monotonicity in the data.
-belowRange Handle values below the interpolation range.
-aboveRange Handle values above the interpolation range.
-majorOrder Set the major order of data storage (row or column).

Incompatibilities

  • One and only one these options must be specified:
    • -atValues, -sequence, -equispaced, -fillIn, or -fileValues
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.
Authors
M. Borland, C. Saunders, R. Soliday, H. Shang

Definition in file sddsinterp.c.

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

Go to the source code of this file.

Functions

double * makeSequence (int64_t *points, double start, double end, double spacing, double *data, int64_t dataRows)
 
double * makeFillInSequence (double *x, int64_t dataRows, int64_t *nPointsReturn)
 
long checkMonotonicity (double *indepValue, int64_t rows)
 
int64_t forceMonotonicity (double *indepValue, double **y, long columns, int64_t rows, unsigned long direction)
 
int64_t combineDuplicatePoints (double *x, double **y, long columns, int64_t rows, double tolerance)
 
int main (int argc, char **argv)
 

Function Documentation

◆ checkMonotonicity()

long checkMonotonicity ( double * indepValue,
int64_t rows )

Definition at line 664 of file sddsinterp.c.

664 {
665 if (rows == 1)
666 return 1;
667 if (indepValue[rows - 1] > indepValue[0]) {
668 while (--rows > 0)
669 if (indepValue[rows] <= indepValue[rows - 1])
670 return 0;
671 return 1;
672 } else {
673 while (--rows > 0)
674 if (indepValue[rows] >= indepValue[rows - 1])
675 return 0;
676 return -1;
677 }
678}

◆ combineDuplicatePoints()

int64_t combineDuplicatePoints ( double * x,
double ** y,
long columns,
int64_t rows,
double tolerance )

Definition at line 680 of file sddsinterp.c.

680 {
681 double xmin, xmax;
682 long column;
683 int64_t i, j;
684
685 find_min_max(&xmin, &xmax, x, rows);
686 if (xmin == xmax)
687 SDDS_Bomb("interpolation data is invalid--no range in independent variable");
688 tolerance *= xmax - xmin;
689 for (i = 1; i < rows; i++) {
690 if (fabs(x[i] - x[i - 1]) <= tolerance) {
691 x[i - 1] = (x[i] + x[i - 1]) / 2;
692 for (column = 0; column < columns; column++)
693 y[column][i - 1] = (y[column][i] + y[column][i - 1]) / 2;
694 for (j = i + 1; j < rows; j++) {
695 x[j - 1] = x[j];
696 for (column = 0; column < columns; column++)
697 y[column][j - 1] = y[column][j];
698 }
699 rows--;
700 i--;
701 }
702 }
703 return rows;
704}
SDDS_Bomb
void SDDS_Bomb(char *message)
Terminates the program after printing an error message and recorded errors.
Definition SDDS_utils.c:342
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

◆ forceMonotonicity()

int64_t forceMonotonicity ( double * indepValue,
double ** y,
long columns,
int64_t rows,
unsigned long direction )

Definition at line 706 of file sddsinterp.c.

706 {
707 int64_t i, j;
708 long column, direction;
709 char *keep;
710 double reference;
711
712 if (rows < 2)
713 return rows;
714
715 if (mode & FORCE_INCREASING)
716 direction = 1;
717 else if (mode & FORCE_DECREASING)
718 direction = -1;
719 else
720 direction = x[1] > x[0] ? 1 : -1;
721
722 keep = tmalloc(sizeof(*keep) * rows);
723 reference = x[0];
724 keep[0] = 1; // Always keep the first point
725 for (i = 1; i < rows; i++) {
726 if (direction * (x[i] - reference) > 0) {
727 reference = x[i];
728 keep[i] = 1;
729 } else {
730 keep[i] = 0;
731 }
732 }
733 for (i = j = 1; i < rows; i++) {
734 if (keep[i]) {
735 if (i != j) {
736 for (column = 0; column < columns; column++)
737 y[column][j] = y[column][i];
738 x[j] = x[i];
739 }
740 j++;
741 }
742 }
743 rows = j;
744 free(keep);
745 return rows;
746}
tmalloc
void * tmalloc(uint64_t size_of_block)
Allocates a memory block of the specified size with zero initialization.
Definition array.c:59

◆ main()

int main ( int argc,
char ** argv )

Definition at line 175 of file sddsinterp.c.

175 {
176 int iArg;
177 char *indepQuantity, **depenQuantity, *fileValuesQuantity, *fileValuesFile, **exclude;
178 long depenQuantities, monotonicity, excludes;
179 char *input, *output;
180 long i, readCode, order, valuesReadCode, fillIn;
181 int64_t rows, row, interpPoints, j;
182 long sequencePoints, combineDuplicates, branch;
183 int32_t *rowFlag;
184 double sequenceStart, sequenceEnd;
185 double sequenceSpacing;
186 unsigned long flags, interpCode, printFlags, forceMonotonic;
187 SCANNED_ARG *scanned;
188 SDDS_DATASET SDDSin, SDDSout, SDDSvalues;
189 OUTRANGE_CONTROL aboveRange, belowRange;
190 double *atValue;
191 long atValues, doNotRead, parallelPages;
192 double *indepValue, **depenValue, *interpPoint, **outputData;
193 unsigned long pipeFlags, majorOrderFlag;
194 FILE *fpPrint;
195 short interpShort = 0, interpShortOrder = -1, *shortValue = NULL, columnMajorOrder = -1;
196 long nextPos;
197
198 SDDS_RegisterProgramName(argv[0]);
199 argc = scanargs(&scanned, argc, argv);
200 if (argc < 3 || argc > (3 + N_OPTIONS))
201 bomb(NULL, USAGE);
202
203 atValue = NULL;
204 atValues = fillIn = 0;
205 output = input = NULL;
206 combineDuplicates = branch = sequencePoints = parallelPages = 0;
207 indepQuantity = NULL;
208 depenQuantity = exclude = NULL;
209 depenQuantities = excludes = 0;
210 aboveRange.flags = belowRange.flags = OUTRANGE_SATURATE;
211 order = 1;
212 readCode = interpPoints = 0;
213 fileValuesFile = fileValuesQuantity = NULL;
214 sequenceStart = sequenceEnd = sequenceSpacing = 0;
215 printFlags = pipeFlags = 0;
216 forceMonotonic = 0;
217 indepValue = interpPoint = NULL;
218 depenValue = outputData = NULL;
219
220 for (iArg = 1; iArg < argc; iArg++) {
221 if (scanned[iArg].arg_type == OPTION) {
222 /* process options here */
223 switch (match_string(scanned[iArg].list[0], option, N_OPTIONS, 0)) {
224 case CLO_MAJOR_ORDER:
225 majorOrderFlag = 0;
226 scanned[iArg].n_items--;
227 if (scanned[iArg].n_items > 0 &&
228 (!scanItemList(&majorOrderFlag, scanned[iArg].list + 1, &scanned[iArg].n_items, 0,
229 "row", -1, NULL, 0, SDDS_ROW_MAJOR_ORDER,
230 "column", -1, NULL, 0, SDDS_COLUMN_MAJOR_ORDER, NULL)))
231 SDDS_Bomb("invalid -majorOrder syntax/values");
232 if (majorOrderFlag & SDDS_COLUMN_MAJOR_ORDER)
233 columnMajorOrder = 1;
234 else if (majorOrderFlag & SDDS_ROW_MAJOR_ORDER)
235 columnMajorOrder = 0;
236 break;
237 case CLO_ORDER:
238 if (scanned[iArg].n_items != 2 ||
239 sscanf(scanned[iArg].list[1], "%ld", &order) != 1 || order < 0)
240 SDDS_Bomb("invalid -order syntax/value");
241 break;
242 case CLO_ATVALUES:
243 if (scanned[iArg].n_items < 2)
244 SDDS_Bomb("invalid -atValues syntax");
245 if (atValue)
246 SDDS_Bomb("give -atValues only once");
247 atValue = tmalloc(sizeof(*atValue) * (atValues = scanned[iArg].n_items - 1));
248 for (i = 0; i < atValues; i++)
249 if (sscanf(scanned[iArg].list[i + 1], "%lf", &atValue[i]) != 1)
250 SDDS_Bomb("invalid -atValues value");
251 break;
252 case CLO_INTERP_SHORT:
253 if (scanned[iArg].n_items == 2) {
254 if (sscanf(scanned[iArg].list[1], "%hd", &interpShortOrder) != 1 ||
255 (interpShortOrder != -1 && interpShortOrder != -2))
256 SDDS_Bomb("invalid -interpShort value; must be -1 or -2");
257 }
258 interpShort = 1;
259 break;
260 case CLO_SEQUENCE:
261 if ((scanned[iArg].n_items != 2 && scanned[iArg].n_items != 4) ||
262 sscanf(scanned[iArg].list[1], "%ld", &sequencePoints) != 1 ||
263 sequencePoints < 2)
264 SDDS_Bomb("invalid -sequence syntax/value");
265 if (scanned[iArg].n_items == 4 &&
266 (sscanf(scanned[iArg].list[2], "%lf", &sequenceStart) != 1 ||
267 sscanf(scanned[iArg].list[3], "%lf", &sequenceEnd) != 1))
268 SDDS_Bomb("invalid -sequence syntax/value");
269 if (sequenceSpacing)
270 SDDS_Bomb("give only one of -sequence and -equispaced");
271 break;
272 case CLO_EQUISPACED:
273 if ((scanned[iArg].n_items != 2 && scanned[iArg].n_items != 4) ||
274 sscanf(scanned[iArg].list[1], "%lf", &sequenceSpacing) != 1 ||
275 sequenceSpacing <= 0)
276 SDDS_Bomb("invalid -equispaced syntax/value");
277 if (scanned[iArg].n_items == 4 &&
278 (sscanf(scanned[iArg].list[2], "%lf", &sequenceStart) != 1 ||
279 sscanf(scanned[iArg].list[3], "%lf", &sequenceEnd) != 1))
280 SDDS_Bomb("invalid -equispaced syntax/values");
281 if (sequencePoints)
282 SDDS_Bomb("give only one of -sequence and -equispaced");
283 break;
284 case CLO_COLUMNS:
285 if (indepQuantity)
286 SDDS_Bomb("only one -columns option may be given");
287 if (scanned[iArg].n_items < 2)
288 SDDS_Bomb("invalid -columns syntax");
289 indepQuantity = scanned[iArg].list[1];
290 if (scanned[iArg].n_items >= 2) {
291 depenQuantity = tmalloc(sizeof(*depenQuantity) * (depenQuantities = scanned[iArg].n_items - 2));
292 for (i = 0; i < depenQuantities; i++)
293 depenQuantity[i] = scanned[iArg].list[i + 2];
294 }
295 break;
296 case CLO_PRINTOUT:
297 if ((scanned[iArg].n_items -= 1) >= 1) {
298 if (!scanItemList(&printFlags, scanned[iArg].list + 1, &scanned[iArg].n_items, 0,
299 "bare", -1, NULL, 0, BARE_PRINTOUT,
300 "stdout", -1, NULL, 0, STDOUT_PRINTOUT, NULL))
301 SDDS_Bomb("invalid -printout syntax");
302 }
303 if (!(printFlags & BARE_PRINTOUT))
304 printFlags |= NORMAL_PRINTOUT;
305 break;
306 case CLO_FILEVALUES:
307 if (scanned[iArg].n_items < 2)
308 SDDS_Bomb("invalid -fileValues syntax");
309 fileValuesFile = scanned[iArg].list[1];
310 scanned[iArg].n_items -= 2;
311 if (!scanItemList(&flags, scanned[iArg].list + 2, &scanned[iArg].n_items, 0,
312 "column", SDDS_STRING, &fileValuesQuantity, 1, 0,
313 "parallelpages", -1, NULL, 0, FILEVALUES_PARALLEL_PAGES, NULL))
314 SDDS_Bomb("invalid -fileValues syntax");
315 if (flags & FILEVALUES_PARALLEL_PAGES)
316 parallelPages = 1;
317 break;
318 case CLO_COMBINEDUPLICATES:
319 SDDS_Bomb("-combineDuplicates option not implemented yet--send email to borland@aps.anl.gov");
320 combineDuplicates = 1;
321 break;
322 case CLO_BRANCH:
323 SDDS_Bomb("-branch option not implemented yet--send email to borland@aps.anl.gov");
324 if (scanned[iArg].n_items != 2 || sscanf(scanned[iArg].list[1], "%ld", &branch) != 1 || branch < 1)
325 SDDS_Bomb("invalid -branch syntax/value");
326 break;
327 case CLO_BELOWRANGE:
328 if ((scanned[iArg].n_items -= 1) < 1 ||
329 !scanItemList(&belowRange.flags, scanned[iArg].list + 1, &scanned[iArg].n_items, 0,
330 "value", SDDS_DOUBLE, &belowRange.value, 1, OUTRANGE_VALUE,
331 "skip", -1, NULL, 0, OUTRANGE_SKIP,
332 "saturate", -1, NULL, 0, OUTRANGE_SATURATE,
333 "extrapolate", -1, NULL, 0, OUTRANGE_EXTRAPOLATE,
334 "wrap", -1, NULL, 0, OUTRANGE_WRAP,
335 "abort", -1, NULL, 0, OUTRANGE_ABORT,
336 "warn", -1, NULL, 0, OUTRANGE_WARN, NULL))
337 SDDS_Bomb("invalid -belowRange syntax/value");
338 if ((i = bitsSet(belowRange.flags & (OUTRANGE_VALUE | OUTRANGE_SKIP | OUTRANGE_SATURATE | OUTRANGE_EXTRAPOLATE | OUTRANGE_WRAP | OUTRANGE_ABORT))) > 1)
339 SDDS_Bomb("incompatible keywords given for -belowRange");
340 if (i != 1)
341 belowRange.flags |= OUTRANGE_SATURATE;
342 break;
343 case CLO_ABOVERANGE:
344 if ((scanned[iArg].n_items -= 1) < 1 ||
345 !scanItemList(&aboveRange.flags, scanned[iArg].list + 1, &scanned[iArg].n_items, 0,
346 "value", SDDS_DOUBLE, &aboveRange.value, 1, OUTRANGE_VALUE,
347 "skip", -1, NULL, 0, OUTRANGE_SKIP,
348 "saturate", -1, NULL, 0, OUTRANGE_SATURATE,
349 "extrapolate", -1, NULL, 0, OUTRANGE_EXTRAPOLATE,
350 "wrap", -1, NULL, 0, OUTRANGE_WRAP,
351 "abort", -1, NULL, 0, OUTRANGE_ABORT,
352 "warn", -1, NULL, 0, OUTRANGE_WARN, NULL))
353 SDDS_Bomb("invalid -aboveRange syntax/value");
354 if ((i = bitsSet(aboveRange.flags & (OUTRANGE_VALUE | OUTRANGE_SKIP | OUTRANGE_SATURATE | OUTRANGE_EXTRAPOLATE | OUTRANGE_WRAP | OUTRANGE_ABORT))) > 1)
355 SDDS_Bomb("incompatible keywords given for -aboveRange");
356 if (i != 1)
357 aboveRange.flags |= OUTRANGE_SATURATE;
358 break;
359 case CLO_PIPE:
360 if (!processPipeOption(scanned[iArg].list + 1, scanned[iArg].n_items - 1, &pipeFlags))
361 SDDS_Bomb("invalid -pipe syntax");
362 break;
363 case CLO_EXCLUDE:
364 if (scanned[iArg].n_items < 2)
365 SDDS_Bomb("invalid -exclude syntax");
366 moveToStringArray(&exclude, &excludes, scanned[iArg].list + 1, scanned[iArg].n_items - 1);
367 break;
368 case CLO_FORCEMONOTONIC:
369 if ((scanned[iArg].n_items -= 1) > 0) {
370 if (!scanItemList(&forceMonotonic, scanned[iArg].list + 1, &scanned[iArg].n_items, 0,
371 "increasing", -1, NULL, 0, FORCE_INCREASING,
372 "decreasing", -1, NULL, 0, FORCE_DECREASING, NULL) ||
373 bitsSet(forceMonotonic) != 1)
374 SDDS_Bomb("invalid -forceMonotonic syntax/value");
375 } else
376 forceMonotonic = FORCE_MONOTONIC;
377 break;
378 case CLO_FILLIN:
379 fillIn = 1;
380 break;
381 default:
382 fprintf(stderr, "error: unknown/ambiguous option: %s\n", scanned[iArg].list[0]);
383 exit(EXIT_FAILURE);
384 break;
385 }
386 } else {
387 if (!input)
388 input = scanned[iArg].list[0];
389 else if (!output)
390 output = scanned[iArg].list[0];
391 else
392 SDDS_Bomb("too many filenames seen");
393 }
394 }
395
396 processFilenames("sddsinterp", &input, &output, pipeFlags, 0, NULL);
397
398 fpPrint = stderr;
399 if (printFlags & STDOUT_PRINTOUT)
400 fpPrint = stdout;
401
402 if (!indepQuantity)
403 SDDS_Bomb("supply the independent quantity name with the -columns option");
404
405 if ((atValues ? 1 : 0) + (fileValuesFile ? 1 : 0) + (sequencePoints ? 1 : 0) + fillIn + (sequenceSpacing > 0 ? 1 : 0) != 1)
406 SDDS_Bomb("you must give one and only one of -atValues, -fileValues, -sequence, -equispaced, and -fillIn");
407
408 if (!SDDS_InitializeInput(&SDDSin, input))
409 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
410
411 excludes = appendToStringArray(&exclude, excludes, indepQuantity);
412 if (!depenQuantities)
413 depenQuantities = appendToStringArray(&depenQuantity, depenQuantities, "*");
414
415 if ((depenQuantities = expandColumnPairNames(&SDDSin, &depenQuantity, NULL, depenQuantities, exclude, excludes, FIND_NUMERIC_TYPE, 0)) <= 0) {
416 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
417 SDDS_Bomb("no dependent quantities selected for interpolation");
418 }
419
420 if (!SDDS_InitializeOutput(&SDDSout, SDDS_BINARY, 0, NULL, "sddsinterp output", output))
421 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
422 if (!SDDS_TransferColumnDefinition(&SDDSout, &SDDSin, indepQuantity, NULL))
423 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
424 if (columnMajorOrder != -1)
425 SDDSout.layout.data_mode.column_major = columnMajorOrder;
426 else
427 SDDSout.layout.data_mode.column_major = SDDSin.layout.data_mode.column_major;
428 if (fileValuesFile && !SDDS_InitializeInput(&SDDSvalues, fileValuesFile))
429 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
430
431 if (SDDS_DefineParameter(&SDDSout, "InterpDataPage", NULL, NULL,
432 "Page of interpolation data file used to create this page", NULL, SDDS_LONG, NULL) < 0 ||
433 SDDS_DefineParameter(&SDDSout, "InterpPointsPage", NULL, NULL, "Page of interpolation points file used to create this page", NULL, SDDS_LONG, NULL) < 0)
434 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
435 for (i = 0; i < depenQuantities; i++)
436 if (!SDDS_TransferColumnDefinition(&SDDSout, &SDDSin, depenQuantity[i], NULL)) {
437 fprintf(stderr, "problem creating interpolated-output column %s\n", depenQuantity[i]);
438 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
439 }
440 if (!SDDS_TransferAllParameterDefinitions(&SDDSout, &SDDSin, SDDS_TRANSFER_KEEPOLD) ||
441 !SDDS_WriteLayout(&SDDSout))
442 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
443
444 doNotRead = 0;
445 interpPoint = NULL;
446 outputData = tmalloc(sizeof(*outputData) * (depenQuantities));
447 depenValue = tmalloc(sizeof(*depenValue) * (depenQuantities));
448 rowFlag = NULL;
449 valuesReadCode = 0;
450
451 while (doNotRead || (readCode = SDDS_ReadPage(&SDDSin)) > 0) {
452 rows = SDDS_CountRowsOfInterest(&SDDSin);
453 if (!(indepValue = SDDS_GetColumnInDoubles(&SDDSin, indepQuantity)))
454 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
455 if (atValues) {
456 interpPoint = atValue;
457 interpPoints = atValues;
458 } else if (fileValuesFile) {
459 if (interpPoint)
460 free(interpPoint);
461 if ((valuesReadCode = SDDS_ReadPage(&SDDSvalues)) == 0)
462 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
463 else if (valuesReadCode == -1) {
464 if (parallelPages) {
465 fprintf(stderr, "warning: file %s ends before file %s\n", fileValuesFile, input);
466 break;
467 } else {
468 /* "rewind" the values file */
469 if (!SDDS_Terminate(&SDDSvalues) ||
470 !SDDS_InitializeInput(&SDDSvalues, fileValuesFile))
471 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
472 if ((valuesReadCode = SDDS_ReadPage(&SDDSvalues)) < 1) {
473 fprintf(stderr, "error: unable to (re)read file %s\n", fileValuesFile);
474 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
475 exit(EXIT_FAILURE);
476 }
477 /* read the next page of the interpolation data file */
478 if ((readCode = SDDS_ReadPage(&SDDSin)) < 1) {
479 if (readCode == -1)
480 break;
481 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
482 }
483 rows = SDDS_CountRowsOfInterest(&SDDSin);
484 if (indepValue)
485 free(indepValue);
486 if (!(indepValue = SDDS_GetColumnInDoubles(&SDDSin, indepQuantity)))
487 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
488 }
489 }
490
491 if (!parallelPages)
492 doNotRead = 1;
493 interpPoints = SDDS_CountRowsOfInterest(&SDDSvalues);
494 interpPoint = SDDS_GetColumnInDoubles(&SDDSvalues, fileValuesQuantity);
495 if (SDDS_NumberOfErrors())
496 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
497 } else if (sequencePoints || sequenceSpacing) {
498 if (interpPoint)
499 free(interpPoint);
500 interpPoints = sequencePoints;
501 if (!(interpPoint = makeSequence(&interpPoints, sequenceStart, sequenceEnd, sequenceSpacing, indepValue, rows)))
502 exit(EXIT_FAILURE);
503 } else {
504 /* fillIn interpolation */
505 if (interpPoint)
506 free(interpPoint);
507 if (!(interpPoint = makeFillInSequence(indepValue, rows, &interpPoints)))
508 exit(EXIT_FAILURE);
509 }
510
511 for (i = 0; i < depenQuantities; i++)
512 outputData[i] = tmalloc(sizeof(*outputData[i]) * interpPoints);
513 rowFlag = trealloc(rowFlag, sizeof(*rowFlag) * interpPoints);
514 for (j = 0; j < interpPoints; j++)
515 rowFlag[j] = 1;
516 for (i = 0; i < depenQuantities; i++) {
517 if (!(depenValue[i] = SDDS_GetColumnInDoubles(&SDDSin, depenQuantity[i])))
518 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
519 }
520 if (forceMonotonic)
521 rows = forceMonotonicity(indepValue, depenValue, depenQuantities, rows, forceMonotonic);
522 else if (combineDuplicates)
523 rows = combineDuplicatePoints(indepValue, depenValue, depenQuantities, rows, 0.0);
524 if ((monotonicity = checkMonotonicity(indepValue, rows)) == 0)
525 SDDS_Bomb("independent data values do not change monotonically or repeated independent values exist");
526 if (interpShort)
527 shortValue = tmalloc(sizeof(*shortValue) * rows);
528
529 for (i = 0; i < depenQuantities; i++) {
530 if (interpShort) {
531 for (row = 0; row < rows; row++) {
532 shortValue[row] = (short)depenValue[i][row];
533 }
534 }
535 for (j = 0; j < interpPoints; j++) {
536 if (!interpShort) {
537 outputData[i][j] = interpolate(depenValue[i], indepValue, rows, interpPoint[j], &belowRange, &aboveRange, order, &interpCode, monotonicity);
538 } else {
539 outputData[i][j] = (double)interp_short(shortValue, indepValue, rows, interpPoint[j], 0, interpShortOrder, &interpCode, &nextPos);
540 }
541 if (interpCode) {
542 if (interpCode & OUTRANGE_ABORT) {
543 fprintf(stderr, "error: value %e is out of range for column %s\n", interpPoint[j], depenQuantity[i]);
544 exit(EXIT_FAILURE);
545 }
546 if (interpCode & OUTRANGE_WARN)
547 fprintf(stderr, "warning: value %e is out of range for column %s\n", interpPoint[j], depenQuantity[i]);
548 if (interpCode & OUTRANGE_SKIP)
549 rowFlag[j] = 0;
550 }
551 }
552 }
553 if (interpShort)
554 free(shortValue);
555 if (!SDDS_StartPage(&SDDSout, interpPoints) ||
556 !SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_NAME, interpPoint, interpPoints, indepQuantity))
557 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
558 if (!SDDS_SetParameters(&SDDSout, SDDS_BY_NAME | SDDS_PASS_BY_VALUE, "InterpDataPage", readCode, "InterpPointsPage", valuesReadCode, NULL) ||
559 !SDDS_CopyParameters(&SDDSout, &SDDSin))
560 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
561 for (i = 0; i < depenQuantities; i++)
562 if (!SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_NAME, outputData[i], interpPoints, depenQuantity[i]))
563 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
564 if (!SDDS_AssertRowFlags(&SDDSout, SDDS_FLAG_ARRAY, rowFlag, rows) || !SDDS_WritePage(&SDDSout))
565 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
566 if (printFlags & BARE_PRINTOUT) {
567 for (j = 0; j < interpPoints; j++)
568 if (rowFlag[j]) {
569 fprintf(fpPrint, "%21.15e ", interpPoint[j]);
570 for (i = 0; i < depenQuantities; i++)
571 fprintf(fpPrint, "%21.15e ", outputData[i][j]);
572 fputc('\n', fpPrint);
573 }
574 } else if (printFlags & NORMAL_PRINTOUT) {
575 for (j = 0; j < interpPoints; j++)
576 if (rowFlag[j]) {
577 fprintf(fpPrint, "%s=%21.15e ", indepQuantity, interpPoint[j]);
578 for (i = 0; i < depenQuantities; i++)
579 fprintf(fpPrint, "%s=%21.15e ", depenQuantity[i], outputData[i][j]);
580 fputc('\n', fpPrint);
581 }
582 }
583 if (indepValue)
584 free(indepValue);
585 indepValue = NULL;
586 for (i = 0; i < depenQuantities; i++) {
587 if (outputData[i])
588 free(outputData[i]);
589 outputData[i] = NULL;
590 if (depenValue[i])
591 free(depenValue[i]);
592 depenValue[i] = NULL;
593 }
594 if (fileValuesFile) {
595 if (interpPoint)
596 free(interpPoint);
597 interpPoint = NULL;
598 }
599 if (rowFlag)
600 free(rowFlag);
601 rowFlag = NULL;
602 }
603
604 if (!SDDS_Terminate(&SDDSin)) {
605 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
606 exit(EXIT_FAILURE);
607 }
608 if (!SDDS_Terminate(&SDDSout)) {
609 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
610 exit(EXIT_FAILURE);
611 }
612 if (fileValuesFile) {
613 if (!SDDS_Terminate(&SDDSvalues)) {
614 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
615 exit(EXIT_FAILURE);
616 }
617 }
618
619 return EXIT_SUCCESS;
620}
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_AssertRowFlags
int32_t SDDS_AssertRowFlags(SDDS_DATASET *SDDS_dataset, uint32_t mode,...)
Sets acceptance flags for rows based on specified criteria.
Definition SDDS_extract.c:148
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_Terminate
int32_t SDDS_Terminate(SDDS_DATASET *SDDS_dataset)
Definition SDDS_input.c:1448
SDDS_ReadPage
int32_t SDDS_ReadPage(SDDS_DATASET *SDDS_dataset)
Definition SDDS_input.c:1006
SDDS_InitializeOutput
int32_t SDDS_InitializeOutput(SDDS_DATASET *SDDS_dataset, int32_t data_mode, int32_t lines_per_row, const char *description, const char *contents, const char *filename)
Initializes the SDDS output dataset.
Definition SDDS_output.c:601
SDDS_WritePage
int32_t SDDS_WritePage(SDDS_DATASET *SDDS_dataset)
Writes the current data table to the output file.
Definition SDDS_output.c:1228
SDDS_WriteLayout
int32_t SDDS_WriteLayout(SDDS_DATASET *SDDS_dataset)
Writes the SDDS layout header to the output file.
Definition SDDS_output.c:899
SDDS_DefineParameter
int32_t SDDS_DefineParameter(SDDS_DATASET *SDDS_dataset, const char *name, const char *symbol, const char *units, const char *description, const char *format_string, int32_t type, char *fixed_value)
Defines a data parameter with a fixed string value.
Definition SDDS_output.c:1472
SDDS_TransferColumnDefinition
int32_t SDDS_TransferColumnDefinition(SDDS_DATASET *target, SDDS_DATASET *source, char *name, char *newName)
Transfers a column definition from a source dataset to a target dataset.
Definition SDDS_transfer.c:35
SDDS_TransferAllParameterDefinitions
int32_t SDDS_TransferAllParameterDefinitions(SDDS_DATASET *SDDS_target, SDDS_DATASET *SDDS_source, uint32_t mode)
Transfers all parameter definitions from a source dataset to a target dataset.
Definition SDDS_transfer.c:343
SDDS_PrintErrors
void SDDS_PrintErrors(FILE *fp, int32_t mode)
Prints recorded error messages to a specified file stream.
Definition SDDS_utils.c:432
SDDS_RegisterProgramName
void SDDS_RegisterProgramName(const char *name)
Registers the executable program name for use in error messages.
Definition SDDS_utils.c:288
SDDS_NumberOfErrors
int32_t SDDS_NumberOfErrors()
Retrieves the number of errors recorded by SDDS library routines.
Definition SDDS_utils.c:304
SDDS_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_DOUBLE
#define SDDS_DOUBLE
Identifier for the double data type.
Definition SDDStypes.h:37
trealloc
void * trealloc(void *old_ptr, uint64_t size_of_block)
Reallocates a memory block to a new size.
Definition array.c:181
bitsSet
long bitsSet(unsigned long data)
Counts the number of set bits (1s) in the given data.
Definition binary.c:52
bomb
void bomb(char *error, char *usage)
Reports error messages to the terminal and aborts the program.
Definition bomb.c:26
interp_short
short interp_short(short *f, double *x, int64_t n, double xo, long warnings, short order, unsigned long *returnCode, long *next_start_pos)
Performs interpolation for short integer data types.
Definition interp.c:281
interpolate
double interpolate(double *f, double *x, int64_t n, double xo, OUTRANGE_CONTROL *belowRange, OUTRANGE_CONTROL *aboveRange, long order, unsigned long *returnCode, long M)
Performs interpolation with range control options.
Definition interp.c:160
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
SDDS_DATASET
Definition SDDS.h:308

◆ makeFillInSequence()

double * makeFillInSequence ( double * x,
int64_t dataRows,
int64_t * nPointsReturn )

Definition at line 622 of file sddsinterp.c.

622 {
623 int64_t i;
624 double dxMin, dx, xMin, xMax;
625
626 dxMin = DBL_MAX;
627 xMin = xMax = x[0];
628 for (i = 1; i < dataRows; i++) {
629 if ((dx = fabs(x[i] - x[i - 1])) < dxMin && dx > 0)
630 dxMin = dx;
631 if (x[i] < xMin)
632 xMin = x[i];
633 if (x[i] > xMax)
634 xMax = x[i];
635 }
636 *nPointsReturn = (int64_t)((xMax - xMin) / dxMin + 1);
637 return makeSequence(nPointsReturn, xMin, xMax, 0.0, x, dataRows);
638}

◆ makeSequence()

double * makeSequence ( int64_t * points,
double start,
double end,
double spacing,
double * data,
int64_t dataRows )

Definition at line 640 of file sddsinterp.c.

640 {
641 int64_t i;
642 double delta, *sequence;
643 if (start == end) {
644 if (!find_min_max(&start, &end, data, dataRows))
645 return NULL;
646 }
647 if (*points > 1) {
648 delta = (end - start) / (*points - 1);
649 } else if (spacing >= 0) {
650 delta = spacing;
651 if (delta == 0)
652 return NULL;
653 *points = (int64_t)((end - start) / delta + 1.5);
654 } else {
655 *points = 1;
656 delta = 0;
657 }
658 sequence = tmalloc(sizeof(*sequence) * (*points));
659 for (i = 0; i < *points; i++)
660 sequence[i] = start + delta * i;
661 return sequence;
662}
Generated by doxygen 1.12.0