190 {
192 long stats;
194 long requests;
195 int64_t count;
196 SCANNED_ARG *scanned;
198 char *independent;
199 int32_t power;
200 int64_t pointsToStat;
201 long pointsToStat0, overlap;
202 int64_t rows, outputRowsMax, outputRows, outputRow;
203 long iArg, code, iStat, iColumn;
204 int64_t startRow, rowOffset;
205 char *input, *output, *windowColumn;
206 double *inputData, *outputData, topLimit, bottomLimit, *inputDataOffset, result, sum1, sum2, *newData;
207 double windowWidth, *windowData, slope, intercept, variance;
208 long windowIndex;
209 unsigned long pipeFlags, scanFlags, majorOrderFlag;
210 char s[100];
211 long lastRegion, region, windowRef, partialOk;
212 short columnMajorOrder = -1;
213
215 argc =
scanargs(&scanned, argc, argv);
216 if (argc < 2) {
217 bomb(
"too few arguments", USAGE);
218 }
219 newData = NULL;
220 result = 0;
221 input = output = NULL;
222 stat = NULL;
223 request = NULL;
224 stats = requests = pipeFlags = 0;
225 pointsToStat = -1;
226 partialOk = 0;
227 overlap = 1;
228 windowColumn = NULL;
229 windowData = NULL;
230
231 for (iArg = 1; iArg < argc; iArg++) {
232 scanFlags = 0;
233 if (scanned[iArg].arg_type == OPTION) {
234
235 switch (code =
match_string(scanned[iArg].list[0], option, N_OPTIONS, 0)) {
236 case SET_MAJOR_ORDER:
237 majorOrderFlag = 0;
238 scanned[iArg].n_items--;
239 if (scanned[iArg].n_items > 0 &&
240 (!
scanItemList(&majorOrderFlag, scanned[iArg].list + 1, &scanned[iArg].n_items, 0,
241 "row", -1, NULL, 0, SDDS_ROW_MAJOR_ORDER,
242 "column", -1, NULL, 0, SDDS_COLUMN_MAJOR_ORDER, NULL)))
243 SDDS_Bomb(
"invalid -majorOrder syntax/values");
244 if (majorOrderFlag & SDDS_COLUMN_MAJOR_ORDER)
245 columnMajorOrder = 1;
246 else if (majorOrderFlag & SDDS_ROW_MAJOR_ORDER)
247 columnMajorOrder = 0;
248 break;
249 case SET_POINTS:
250 if (scanned[iArg].n_items != 2 || sscanf(scanned[iArg].list[1], "%" SCNd64, &pointsToStat) != 1 ||
251 (pointsToStat <= 2 && pointsToStat != 0))
253 break;
254 case SET_NOOVERLAP:
255 overlap = 0;
256 break;
257 case SET_MAXIMUM:
258 case SET_MINIMUM:
259 case SET_MEAN:
260 case SET_STANDARDDEVIATION:
261 case SET_RMS:
262 case SET_SIGMA:
263 case SET_SAMPLE:
264 case SET_MEDIAN:
265 if (scanned[iArg].n_items < 2) {
266 fprintf(stderr, "error: invalid -%s syntax\n", option[code]);
267 exit(EXIT_FAILURE);
268 }
269 if (!
scanItemList(&scanFlags, scanned[iArg].list, &scanned[iArg].n_items,
270 SCANITEMLIST_UNKNOWN_VALUE_OK | SCANITEMLIST_REMOVE_USED_ITEMS | SCANITEMLIST_IGNORE_VALUELESS,
271 "toplimit",
SDDS_DOUBLE, &topLimit, 1, TOPLIMIT_GIVEN,
272 "bottomlimit",
SDDS_DOUBLE, &bottomLimit, 1, BOTTOMLIMIT_GIVEN, NULL)) {
273 sprintf(s, "invalid -%s syntax", scanned[iArg].list[0]);
275 }
276 requests = addStatRequests(&request, requests, scanned[iArg].list + 1, scanned[iArg].n_items - 1, code, scanFlags);
277 request[requests - 1].topLimit = topLimit;
278 request[requests - 1].bottomLimit = bottomLimit;
279 break;
280 case SET_SUM:
281 if (scanned[iArg].n_items < 2) {
282 fprintf(stderr, "error: invalid -%s syntax\n", option[code]);
283 exit(EXIT_FAILURE);
284 }
285 power = 1;
286 if (!
scanItemList(&scanFlags, scanned[iArg].list, &scanned[iArg].n_items,
287 SCANITEMLIST_UNKNOWN_VALUE_OK | SCANITEMLIST_REMOVE_USED_ITEMS | SCANITEMLIST_IGNORE_VALUELESS,
289 "toplimit",
SDDS_DOUBLE, &topLimit, 1, TOPLIMIT_GIVEN,
290 "bottomlimit",
SDDS_DOUBLE, &bottomLimit, 1, BOTTOMLIMIT_GIVEN, NULL))
292 requests = addStatRequests(&request, requests, scanned[iArg].list + 1, scanned[iArg].n_items - 1, code, scanFlags);
293 request[requests - 1].sumPower = power;
294 request[requests - 1].topLimit = topLimit;
295 request[requests - 1].bottomLimit = bottomLimit;
296 break;
297 case SET_SLOPE:
298 if (scanned[iArg].n_items < 2) {
299 fprintf(stderr, "error: invalid -%s syntax\n", option[code]);
300 exit(EXIT_FAILURE);
301 }
302 if (!
scanItemList(&scanFlags, scanned[iArg].list, &scanned[iArg].n_items,
303 SCANITEMLIST_UNKNOWN_VALUE_OK | SCANITEMLIST_REMOVE_USED_ITEMS | SCANITEMLIST_IGNORE_VALUELESS,
304 "independent",
SDDS_STRING, &independent, 1, INDEPENDENT_GIVEN,
305 NULL) ||
306 !(scanFlags & INDEPENDENT_GIVEN))
308 requests = addStatRequests(&request, requests, scanned[iArg].list + 1, scanned[iArg].n_items - 1, code, scanFlags);
309 request[requests - 1].independentColumn = independent;
310 request[requests - 1].topLimit = request[requests - 1].bottomLimit = 0;
311 break;
312 case SET_PIPE:
313 if (!
processPipeOption(scanned[iArg].list + 1, scanned[iArg].n_items - 1, &pipeFlags))
315 break;
316 case SET_WINDOW:
317 windowWidth = -1;
318 windowColumn = NULL;
319 scanned[iArg].n_items -= 1;
320 if (!
scanItemList(&scanFlags, scanned[iArg].list + 1, &scanned[iArg].n_items, 0,
323 !windowColumn ||
324 !strlen(windowColumn) ||
325 windowWidth <= 0)
326 SDDS_Bomb(
"invalid -window syntax/values");
327 break;
328 case SET_PARTIALOK:
329 partialOk = 1;
330 break;
331 default:
332 fprintf(stderr, "error: unknown option '%s' given\n", scanned[iArg].list[0]);
333 exit(EXIT_FAILURE);
334 break;
335 }
336 } else {
337
338 if (!input)
339 input = scanned[iArg].list[0];
340 else if (!output)
341 output = scanned[iArg].list[0];
342 else
344 }
345 }
346
347 if (pointsToStat < 0 && !windowColumn) {
348 pointsToStat = 10;
349 }
351
352 if (!requests)
354
357
358 if (!(stat = compileStatDefinitions(&inData, request, requests)))
360 stats = requests;
361#ifdef DEBUG
362 fprintf(stderr, "%ld stats\n", stats);
363 for (iStat = 0; iStat < stats; iStat++) {
364 for (iColumn = 0; iColumn < stat[iStat].sourceColumns; iColumn++) {
365 fprintf(stderr, "iStat=%ld iColumn=%ld source=%s result=%s\n", iStat, iColumn, stat[iStat].sourceColumn[iColumn], stat[iStat].resultColumn[iColumn]);
366 if (stat[iStat].flags & BOTTOMLIMIT_GIVEN)
367 fprintf(stderr, " bottom = %e\n", stat[iStat].bottomLimit);
368 if (stat[iStat].flags & TOPLIMIT_GIVEN)
369 fprintf(stderr, " top = %e\n", stat[iStat].topLimit);
370 }
371 }
372#endif
373
374 if (!setupOutputFile(&outData, output, &inData, stat, stats, columnMajorOrder))
376
377 if (windowColumn) {
381 SDDS_Bomb(
"Window column is not numeric");
382 }
383
384 outputData = NULL;
385 outputRowsMax = 0;
386 pointsToStat0 = pointsToStat;
389 pointsToStat = pointsToStat0;
390 if (pointsToStat == 0)
391 pointsToStat = rows;
394 }
395 if (!windowColumn) {
396 if (rows < pointsToStat) {
397 if (partialOk)
398 pointsToStat = rows;
399 else
400 continue;
401 }
402 if (overlap) {
403 outputRows = rows - pointsToStat + 1;
404 } else
405 outputRows = rows / pointsToStat;
406 } else
407 outputRows = rows;
408
413 if (outputRows > outputRowsMax &&
414 !(outputData =
SDDS_Realloc(outputData,
sizeof(*outputData) * (outputRowsMax = outputRows))))
416 for (iStat = 0; iStat < stats; iStat++) {
417 for (iColumn = 0; iColumn < stat[iStat].sourceColumns; iColumn++) {
418 double *indepData;
419 lastRegion = 0;
420 windowRef = 0;
421 indepData = NULL;
424 if (stat[iStat].independentColumn &&
427 for (outputRow = startRow = 0; outputRow < outputRows; outputRow++, startRow += (overlap ? 1 : pointsToStat)) {
428 if (windowColumn) {
429 short windowFound = 0;
430 if (overlap) {
431 windowRef += 1;
432 lastRegion = 0;
433 }
434 for (pointsToStat = 1; pointsToStat < outputRows - startRow; pointsToStat++) {
435 region = (windowData[startRow + pointsToStat] - windowData[windowRef]) / windowWidth;
436 if (region != lastRegion) {
437 lastRegion = region;
438 windowFound = 1;
439 break;
440 }
441 }
442 if (!windowFound && pointsToStat < 2)
443 break;
444 if (startRow + pointsToStat > rows) {
445 pointsToStat = rows - startRow - 1;
446 if (pointsToStat <= 0)
447 break;
448 }
449#ifdef DEBUG
450 fprintf(stderr, "row=%" PRId64 " pointsToStat=%" PRId64 " delta=%.9lf (%.9lf -> %.9lf)\n", startRow, pointsToStat, windowData[startRow + pointsToStat - 1] - windowData[startRow], windowData[startRow], windowData[startRow + pointsToStat - 1]);
451#endif
452 }
453 inputDataOffset = inputData + startRow;
454 switch (stat[iStat].optionCode) {
455 case SET_MAXIMUM:
456 result = -DBL_MAX;
457 for (rowOffset = 0; rowOffset < pointsToStat; rowOffset++) {
458 if ((stat[iStat].flags & TOPLIMIT_GIVEN && inputDataOffset[rowOffset] > stat[iStat].topLimit) ||
459 (stat[iStat].flags & BOTTOMLIMIT_GIVEN && inputDataOffset[rowOffset] < stat[iStat].bottomLimit))
460 continue;
461 if (inputDataOffset[rowOffset] > result)
462 result = inputDataOffset[rowOffset];
463 }
464 break;
465 case SET_MINIMUM:
466 result = DBL_MAX;
467 for (rowOffset = 0; rowOffset < pointsToStat; rowOffset++) {
468 if ((stat[iStat].flags & TOPLIMIT_GIVEN && inputDataOffset[rowOffset] > stat[iStat].topLimit) ||
469 (stat[iStat].flags & BOTTOMLIMIT_GIVEN && inputDataOffset[rowOffset] < stat[iStat].bottomLimit))
470 continue;
471 if (inputDataOffset[rowOffset] < result)
472 result = inputDataOffset[rowOffset];
473 }
474 break;
475 case SET_MEAN:
476 result = 0;
477 count = 0;
478 for (rowOffset = 0; rowOffset < pointsToStat; rowOffset++) {
479 if ((stat[iStat].flags & TOPLIMIT_GIVEN && inputDataOffset[rowOffset] > stat[iStat].topLimit) ||
480 (stat[iStat].flags & BOTTOMLIMIT_GIVEN && inputDataOffset[rowOffset] < stat[iStat].bottomLimit))
481 continue;
482 result += inputDataOffset[rowOffset];
483 count++;
484 }
485 if (count)
486 result /= count;
487 else
488 result = DBL_MAX;
489 break;
490 case SET_MEDIAN:
491 result = 0;
492 count = 0;
493 for (rowOffset = 0; rowOffset < pointsToStat; rowOffset++) {
494 if ((stat[iStat].flags & TOPLIMIT_GIVEN && inputDataOffset[rowOffset] > stat[iStat].topLimit) ||
495 (stat[iStat].flags & BOTTOMLIMIT_GIVEN && inputDataOffset[rowOffset] < stat[iStat].bottomLimit))
496 continue;
497 newData =
SDDS_Realloc(newData,
sizeof(*newData) * (count + 1));
498 newData[count] = inputDataOffset[rowOffset];
499 count++;
500 }
501 if (count) {
503 result = DBL_MAX;
504 free(newData);
505 newData = NULL;
506 } else
507 result = DBL_MAX;
508 break;
509 case SET_STANDARDDEVIATION:
510 case SET_SIGMA:
511 sum1 = sum2 = count = 0;
512 for (rowOffset = 0; rowOffset < pointsToStat; rowOffset++) {
513 if ((stat[iStat].flags & TOPLIMIT_GIVEN && inputDataOffset[rowOffset] > stat[iStat].topLimit) ||
514 (stat[iStat].flags & BOTTOMLIMIT_GIVEN && inputDataOffset[rowOffset] < stat[iStat].bottomLimit))
515 continue;
516 sum1 += inputDataOffset[rowOffset];
517 sum2 += inputDataOffset[rowOffset] * inputDataOffset[rowOffset];
518 count++;
519 }
520 if (count > 1) {
521 if ((result = sum2 / count - sqr(sum1 / count)) <= 0)
522 result = 0;
523 else
524 result = sqrt(result * count / (count - 1.0));
525 if (stat[iStat].optionCode == SET_SIGMA)
526 result /= sqrt(count);
527 }
528 break;
529 case SET_RMS:
530 sum2 = count = 0;
531 for (rowOffset = 0; rowOffset < pointsToStat; rowOffset++) {
532 if ((stat[iStat].flags & TOPLIMIT_GIVEN && inputDataOffset[rowOffset] > stat[iStat].topLimit) ||
533 (stat[iStat].flags & BOTTOMLIMIT_GIVEN && inputDataOffset[rowOffset] < stat[iStat].bottomLimit))
534 continue;
535 sum2 += inputDataOffset[rowOffset] * inputDataOffset[rowOffset];
536 count++;
537 }
538 if (count > 0)
539 result = sqrt(sum2 / count);
540 else
541 result = DBL_MAX;
542 break;
543 case SET_SUM:
544 sum1 = count = 0;
545 for (rowOffset = 0; rowOffset < pointsToStat; rowOffset++) {
546 if ((stat[iStat].flags & TOPLIMIT_GIVEN && inputDataOffset[rowOffset] > stat[iStat].topLimit) ||
547 (stat[iStat].flags & BOTTOMLIMIT_GIVEN && inputDataOffset[rowOffset] < stat[iStat].bottomLimit))
548 continue;
549 sum1 +=
ipow(inputDataOffset[rowOffset], stat[iStat].sumPower);
550 count++;
551 }
552 if (count > 0)
553 result = sum1;
554 else
555 result = DBL_MAX;
556 break;
557 case SET_SLOPE:
558 if (!
unweightedLinearFit(indepData + startRow, inputDataOffset, pointsToStat, &slope, &intercept,
559 &variance)) {
560 result = DBL_MAX;
561 } else
562 result = slope;
563 break;
564 case SET_SAMPLE:
565 result = DBL_MAX;
566 for (rowOffset = 0; rowOffset < pointsToStat; rowOffset++) {
567 if ((stat[iStat].flags & TOPLIMIT_GIVEN && inputDataOffset[rowOffset] > stat[iStat].topLimit) ||
568 (stat[iStat].flags & BOTTOMLIMIT_GIVEN && inputDataOffset[rowOffset] < stat[iStat].bottomLimit))
569 continue;
570 result = inputDataOffset[rowOffset];
571 break;
572 }
573 break;
574 default:
575 fprintf(stderr, "Unknown statistics code %ld in sddsrunave\n", stat[iStat].optionCode);
576 exit(EXIT_FAILURE);
577 break;
578 }
579 outputData[outputRow] = result;
580 }
583 free(inputData);
584 }
585 }
586 if (windowColumn)
587 free(windowData);
590 }
593 return EXIT_FAILURE;
594 }
597 return EXIT_FAILURE;
598 }
599 if (outputData)
600 free(outputData);
601 return EXIT_SUCCESS;
602}
int32_t SDDS_CopyParameters(SDDS_DATASET *SDDS_target, SDDS_DATASET *SDDS_source)
int32_t SDDS_CopyArrays(SDDS_DATASET *SDDS_target, SDDS_DATASET *SDDS_source)
int32_t SDDS_StartPage(SDDS_DATASET *SDDS_dataset, int64_t expected_n_rows)
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.
int32_t SDDS_WritePage(SDDS_DATASET *SDDS_dataset)
Writes the current data table to the output file.
int32_t SDDS_GetColumnIndex(SDDS_DATASET *SDDS_dataset, char *name)
Retrieves the index of a named column in the SDDS dataset.
void SDDS_PrintErrors(FILE *fp, int32_t mode)
Prints recorded error messages to a specified file stream.
void SDDS_RegisterProgramName(const char *name)
Registers the executable program name for use in error messages.
int32_t SDDS_GetColumnType(SDDS_DATASET *SDDS_dataset, int32_t index)
Retrieves the data type of a column in the SDDS dataset by its index.
#define SDDS_STRING
Identifier for the string data type.
#define SDDS_LONG
Identifier for the signed 32-bit integer data type.
#define SDDS_DOUBLE
Identifier for the double data type.
#define SDDS_NUMERIC_TYPE(type)
Checks if the given type identifier corresponds to any numeric type.
void bomb(char *error, char *usage)
Reports error messages to the terminal and aborts the program.
double ipow(const double x, const int64_t p)
Compute x raised to the power p (x^p).
long unweightedLinearFit(double *xData, double *yData, long nData, double *slope, double *intercept, double *variance)
Performs an unweighted linear fit on the provided data.
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.
int scanargs(SCANNED_ARG **scanned, int argc, char **argv)
long processPipeOption(char **item, long items, unsigned long *flags)
void processFilenames(char *programName, char **input, char **output, unsigned long pipeFlags, long noWarnings, long *tmpOutputUsed)
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.