233 {
235 long stats;
237 long requests;
238 SCANNED_ARG *scanned;
240
241 int32_t power;
242 long i_arg, code, iStat, tmpFileUsed, iColumn, posColIndex;
243 int64_t rows, row, count;
244 long noWarnings, maxSourceColumns;
245 char *input, *output, *positionColumn, **posColumnName;
246 double **inputData, *outputData, value1, value2, topLimit, bottomLimit;
247 unsigned long pipeFlags, scanFlags, majorOrderFlag;
248 char s[100];
249 double *statWorkArray;
250 double quartilePoint[2] = {25.0, 75.0}, quartileResult[2];
251 double decilePoint[2] = {10.0, 90.0}, decileResult[2];
252 double percent;
253 short columnMajorOrder = -1;
254 int threads = 1;
255
257 argc =
scanargs(&scanned, argc, argv);
258 if (argc < 2) {
259 bomb(
"too few arguments", USAGE);
260 }
261
262 posColumnName = NULL;
263 input = output = positionColumn = NULL;
264 stat = NULL;
265 request = NULL;
266 stats = requests = pipeFlags = 0;
267 noWarnings = 0;
268 outputData = NULL;
269 statWorkArray = NULL;
270
271 for (i_arg = 1; i_arg < argc; i_arg++) {
272 scanFlags = 0;
273 if (scanned[i_arg].arg_type == OPTION) {
274
275 switch (code =
match_string(scanned[i_arg].list[0], option, N_OPTIONS, 0)) {
276 case SET_MAXIMUM:
277 case SET_MINIMUM:
278 case SET_MEAN:
279 case SET_MEDIAN:
280 case SET_STANDARDDEVIATION:
281 case SET_RMS:
282 case SET_SIGMA:
283 case SET_MAD:
284 case SET_COUNT:
285 case SET_DRANGE:
286 case SET_QRANGE:
287 case SET_SMALLEST:
288 case SET_LARGEST:
289 case SET_SPREAD_ARG:
290 if (scanned[i_arg].n_items < 3) {
291 fprintf(stderr, "error: invalid -%s syntax\n", option[code]);
292 exit(EXIT_FAILURE);
293 }
294 if (!
scanItemList(&scanFlags, scanned[i_arg].list, &scanned[i_arg].n_items,
295 SCANITEMLIST_UNKNOWN_VALUE_OK | SCANITEMLIST_REMOVE_USED_ITEMS |
296 SCANITEMLIST_IGNORE_VALUELESS,
297 "positionColumn",
SDDS_STRING, &positionColumn, 1, POSITIONCOLUMN_GIVEN,
298 "toplimit",
SDDS_DOUBLE, &topLimit, 1, TOPLIMIT_GIVEN,
299 "bottomlimit",
SDDS_DOUBLE, &bottomLimit, 1, BOTTOMLIMIT_GIVEN, NULL)) {
300 sprintf(s, "invalid -%s syntax", scanned[i_arg].list[0]);
302 }
303 requests = addStatRequests(&request, requests, scanned[i_arg].list + 1, scanned[i_arg].n_items - 1, code, scanFlags);
304 request[requests - 1].topLimit = topLimit;
305 request[requests - 1].bottomLimit = bottomLimit;
306 if (positionColumn) {
307 if (code == SET_MAXIMUM || code == SET_MINIMUM || code == SET_LARGEST || code == SET_SMALLEST)
308 SDDS_CopyString(&request[requests - 1].positionColumn, positionColumn);
309 free(positionColumn);
310 positionColumn = NULL;
311 }
312 break;
313 case SET_PERCENTILE:
314 if (scanned[i_arg].n_items < 3) {
315 fprintf(stderr, "error: invalid -%s syntax\n", option[code]);
316 exit(EXIT_FAILURE);
317 }
318 if (!
scanItemList(&scanFlags, scanned[i_arg].list, &scanned[i_arg].n_items,
319 SCANITEMLIST_UNKNOWN_VALUE_OK | SCANITEMLIST_REMOVE_USED_ITEMS |
320 SCANITEMLIST_IGNORE_VALUELESS,
322 "toplimit",
SDDS_DOUBLE, &topLimit, 1, TOPLIMIT_GIVEN,
323 "bottomlimit",
SDDS_DOUBLE, &bottomLimit, 1, BOTTOMLIMIT_GIVEN, NULL) ||
324 !(scanFlags & PERCENT_GIVEN) || percent <= 0 || percent >= 100)
326 requests = addStatRequests(&request, requests, scanned[i_arg].list + 1, scanned[i_arg].n_items - 1, code, scanFlags);
327 request[requests - 1].percent = percent;
328 request[requests - 1].topLimit = topLimit;
329 request[requests - 1].bottomLimit = bottomLimit;
330 break;
331 case SET_SUM:
332 if (scanned[i_arg].n_items < 3) {
333 fprintf(stderr, "error: invalid -%s syntax\n", option[code]);
334 exit(EXIT_FAILURE);
335 }
336 power = 1;
337 if (!
scanItemList(&scanFlags, scanned[i_arg].list, &scanned[i_arg].n_items,
338 SCANITEMLIST_UNKNOWN_VALUE_OK | SCANITEMLIST_REMOVE_USED_ITEMS | SCANITEMLIST_IGNORE_VALUELESS,
340 "toplimit",
SDDS_DOUBLE, &topLimit, 1, TOPLIMIT_GIVEN,
341 "bottomlimit",
SDDS_DOUBLE, &bottomLimit, 1, BOTTOMLIMIT_GIVEN, NULL))
343 requests = addStatRequests(&request, requests, scanned[i_arg].list + 1, scanned[i_arg].n_items - 1, code, scanFlags);
344 request[requests - 1].sumPower = power;
345 request[requests - 1].topLimit = topLimit;
346 request[requests - 1].bottomLimit = bottomLimit;
347 break;
348 case SET_PIPE:
349 if (!
processPipeOption(scanned[i_arg].list + 1, scanned[i_arg].n_items - 1, &pipeFlags))
351 break;
352 case SET_NOWARNINGS:
353 noWarnings = 1;
354 break;
355 case SET_MAJOR_ORDER:
356 majorOrderFlag = 0;
357 scanned[i_arg].n_items--;
358 if (scanned[i_arg].n_items > 0 &&
359 (!
scanItemList(&majorOrderFlag, scanned[i_arg].list + 1, &scanned[i_arg].n_items, 0,
360 "row", -1, NULL, 0, SDDS_ROW_MAJOR_ORDER,
361 "column", -1, NULL, 0, SDDS_COLUMN_MAJOR_ORDER, NULL)))
362 SDDS_Bomb(
"invalid -majorOrder syntax/values");
363 if (majorOrderFlag & SDDS_COLUMN_MAJOR_ORDER)
364 columnMajorOrder = 1;
365 else if (majorOrderFlag & SDDS_ROW_MAJOR_ORDER)
366 columnMajorOrder = 0;
367 break;
368 case SET_THREADS:
369 if (scanned[i_arg].n_items != 2 ||
370 !sscanf(scanned[i_arg].list[1], "%d", &threads) || threads < 1)
372 break;
373 default:
374 fprintf(stderr, "error: unknown option '%s' given\n", scanned[i_arg].list[0]);
375 exit(EXIT_FAILURE);
376 break;
377 }
378 } else {
379
380 if (!input)
381 input = scanned[i_arg].list[0];
382 else if (!output)
383 output = scanned[i_arg].list[0];
384 else
386 }
387 }
388
389 processFilenames(
"sddsrowstats", &input, &output, pipeFlags, noWarnings, &tmpFileUsed);
390
391 if (!requests)
393
396
397 if (!(stat = compileStatDefinitions(&inData, request, requests, &stats, noWarnings))) {
399 exit(EXIT_FAILURE);
400 }
401 if (stats < 0)
402 SDDS_Bomb(
"No valid statistics requests.");
403 for (iStat = maxSourceColumns = 0; iStat < stats; iStat++) {
404 if (stat[iStat].sourceColumns > maxSourceColumns)
405 maxSourceColumns = stat[iStat].sourceColumns;
406 }
407 if (!(statWorkArray = malloc(sizeof(*statWorkArray) * maxSourceColumns)))
408 SDDS_Bomb(
"allocation failure (statWorkArray)");
409
410 if (!setupOutputFile(&outData, output, &inData, stat, stats, columnMajorOrder)) {
412 exit(EXIT_FAILURE);
413 }
414
415 inputData = NULL;
416
420 exit(EXIT_FAILURE);
421 }
423 if (!(outputData = (double *)malloc(sizeof(*outputData) * rows))) {
425 }
426 if (!(posColumnName = (char **)malloc(sizeof(*posColumnName) * rows))) {
428 }
429 for (iStat = 0; iStat < stats; iStat++) {
430 if (!(inputData = (double **)malloc(sizeof(*inputData) * stat[iStat].sourceColumns))) {
432 }
433 for (iColumn = 0; iColumn < stat[iStat].sourceColumns; iColumn++) {
436 }
437 }
438 for (row = 0; row < rows; row++)
439 outputData[row] = DBL_MAX;
440 switch (stat[iStat].optionCode) {
441 case SET_MINIMUM:
442 for (row = 0; row < rows; row++) {
443 value1 = DBL_MAX;
444 posColIndex = 0;
445 posColumnName[row] = NULL;
446 for (iColumn = 0; iColumn < stat[iStat].sourceColumns; iColumn++) {
447 if (stat[iStat].flags & TOPLIMIT_GIVEN && inputData[iColumn][row] > stat[iStat].topLimit)
448 continue;
449 if (stat[iStat].flags & BOTTOMLIMIT_GIVEN && inputData[iColumn][row] < stat[iStat].bottomLimit)
450 continue;
451 if (inputData[iColumn][row] < value1) {
452 value1 = inputData[iColumn][row];
453 posColIndex = iColumn;
454 }
455 }
456 outputData[row] = value1;
457 if (stat[iStat].positionColumn)
458 posColumnName[row] = stat[iStat].sourceColumn[posColIndex];
459 }
460 break;
461 case SET_MAXIMUM:
462 for (row = 0; row < rows; row++) {
463 posColIndex = 0;
464 value1 = -DBL_MAX;
465 posColumnName[row] = NULL;
466 for (iColumn = 0; iColumn < stat[iStat].sourceColumns; iColumn++) {
467 if (stat[iStat].flags & TOPLIMIT_GIVEN && inputData[iColumn][row] > stat[iStat].topLimit)
468 continue;
469 if (stat[iStat].flags & BOTTOMLIMIT_GIVEN && inputData[iColumn][row] < stat[iStat].bottomLimit)
470 continue;
471 if (inputData[iColumn][row] > value1) {
472 posColIndex = iColumn;
473 value1 = inputData[iColumn][row];
474 }
475 }
476 outputData[row] = value1;
477 if (stat[iStat].positionColumn)
478 posColumnName[row] = stat[iStat].sourceColumn[posColIndex];
479 }
480 break;
481 case SET_MEAN:
482 for (row = 0; row < rows; row++) {
483 value1 = 0;
484 count = 0;
485 for (iColumn = 0; iColumn < stat[iStat].sourceColumns; iColumn++) {
486 if (stat[iStat].flags & TOPLIMIT_GIVEN && inputData[iColumn][row] > stat[iStat].topLimit)
487 continue;
488 if (stat[iStat].flags & BOTTOMLIMIT_GIVEN && inputData[iColumn][row] < stat[iStat].bottomLimit)
489 continue;
490 value1 += inputData[iColumn][row];
491 count++;
492 }
493 if (count)
494 outputData[row] = value1 / count;
495 }
496 break;
497 case SET_MEDIAN:
498 for (row = 0; row < rows; row++) {
499 for (iColumn = count = 0; iColumn < stat[iStat].sourceColumns; iColumn++) {
500 if (stat[iStat].flags & TOPLIMIT_GIVEN && inputData[iColumn][row] > stat[iStat].topLimit)
501 continue;
502 if (stat[iStat].flags & BOTTOMLIMIT_GIVEN && inputData[iColumn][row] < stat[iStat].bottomLimit)
503 continue;
504 statWorkArray[count] = inputData[iColumn][row];
505 count++;
506 }
507 if (count)
509 }
510 break;
511 case SET_STANDARDDEVIATION:
512 for (row = 0; row < rows; row++) {
513 value1 = 0;
514 value2 = 0;
515 count = 0;
516 for (iColumn = 0; iColumn < stat[iStat].sourceColumns; iColumn++) {
517 if (stat[iStat].flags & TOPLIMIT_GIVEN && inputData[iColumn][row] > stat[iStat].topLimit)
518 continue;
519 if (stat[iStat].flags & BOTTOMLIMIT_GIVEN && inputData[iColumn][row] < stat[iStat].bottomLimit)
520 continue;
521 value1 += inputData[iColumn][row];
522 value2 += inputData[iColumn][row] * inputData[iColumn][row];
523 count++;
524 }
525 if (count > 1) {
526 if ((value1 = value2 / count - sqr(value1 / count)) <= 0)
527 outputData[row] = 0;
528 else
529 outputData[row] = sqrt(value1 * count / (count - 1.0));
530 }
531 }
532 break;
533 case SET_SIGMA:
534 for (row = 0; row < rows; row++) {
535 value1 = 0;
536 value2 = 0;
537 count = 0;
538 for (iColumn = 0; iColumn < stat[iStat].sourceColumns; iColumn++) {
539 if (stat[iStat].flags & TOPLIMIT_GIVEN && inputData[iColumn][row] > stat[iStat].topLimit)
540 continue;
541 if (stat[iStat].flags & BOTTOMLIMIT_GIVEN && inputData[iColumn][row] < stat[iStat].bottomLimit)
542 continue;
543 value1 += inputData[iColumn][row];
544 value2 += inputData[iColumn][row] * inputData[iColumn][row];
545 count++;
546 }
547 if (count > 1) {
548 if ((value1 = value2 / count - sqr(value1 / count)) <= 0)
549 outputData[row] = 0;
550 else
551 outputData[row] = sqrt(value1 / (count - 1.0));
552 }
553 }
554 break;
555 case SET_RMS:
556 for (row = 0; row < rows; row++) {
557 value1 = 0;
558 count = 0;
559 for (iColumn = 0; iColumn < stat[iStat].sourceColumns; iColumn++) {
560 if (stat[iStat].flags & TOPLIMIT_GIVEN && inputData[iColumn][row] > stat[iStat].topLimit)
561 continue;
562 if (stat[iStat].flags & BOTTOMLIMIT_GIVEN && inputData[iColumn][row] < stat[iStat].bottomLimit)
563 continue;
564 value1 += sqr(inputData[iColumn][row]);
565 count++;
566 }
567 if (count)
568 outputData[row] = sqrt(value1 / count);
569 }
570 break;
571 case SET_SUM:
572 for (row = 0; row < rows; row++) {
573 value1 = 0;
574 count = 0;
575 for (iColumn = 0; iColumn < stat[iStat].sourceColumns; iColumn++) {
576 if (stat[iStat].flags & TOPLIMIT_GIVEN && inputData[iColumn][row] > stat[iStat].topLimit)
577 continue;
578 if (stat[iStat].flags & BOTTOMLIMIT_GIVEN && inputData[iColumn][row] < stat[iStat].bottomLimit)
579 continue;
580 value1 +=
ipow(inputData[iColumn][row], stat[iStat].sumPower);
581 count++;
582 }
583 if (count)
584 outputData[row] = value1;
585 }
586 break;
587 case SET_COUNT:
588 for (row = 0; row < rows; row++) {
589 count = 0;
590 for (iColumn = 0; iColumn < stat[iStat].sourceColumns; iColumn++) {
591 if (stat[iStat].flags & TOPLIMIT_GIVEN && inputData[iColumn][row] > stat[iStat].topLimit)
592 continue;
593 if (stat[iStat].flags & BOTTOMLIMIT_GIVEN && inputData[iColumn][row] < stat[iStat].bottomLimit)
594 continue;
595 count++;
596 }
597 outputData[row] = count;
598 }
599 break;
600 case SET_MAD:
601 for (row = 0; row < rows; row++) {
602 for (iColumn = count = value1 = 0; iColumn < stat[iStat].sourceColumns; iColumn++) {
603 if (stat[iStat].flags & TOPLIMIT_GIVEN && inputData[iColumn][row] > stat[iStat].topLimit)
604 continue;
605 if (stat[iStat].flags & BOTTOMLIMIT_GIVEN && inputData[iColumn][row] < stat[iStat].bottomLimit)
606 continue;
607 statWorkArray[count] = inputData[iColumn][row];
608 count++;
609 }
610 if (count)
612 }
613 break;
614 case SET_DRANGE:
615 for (row = 0; row < rows; row++) {
616 for (iColumn = count = value1 = 0; iColumn < stat[iStat].sourceColumns; iColumn++) {
617 if (stat[iStat].flags & TOPLIMIT_GIVEN && inputData[iColumn][row] > stat[iStat].topLimit)
618 continue;
619 if (stat[iStat].flags & BOTTOMLIMIT_GIVEN && inputData[iColumn][row] < stat[iStat].bottomLimit)
620 continue;
621 statWorkArray[count] = inputData[iColumn][row];
622 count++;
623 }
625 outputData[row] = decileResult[1] - decileResult[0];
626 }
627 break;
628 case SET_QRANGE:
629 for (row = 0; row < rows; row++) {
630 for (iColumn = count = value1 = 0; iColumn < stat[iStat].sourceColumns; iColumn++) {
631 if (stat[iStat].flags & TOPLIMIT_GIVEN && inputData[iColumn][row] > stat[iStat].topLimit)
632 continue;
633 if (stat[iStat].flags & BOTTOMLIMIT_GIVEN && inputData[iColumn][row] < stat[iStat].bottomLimit)
634 continue;
635 statWorkArray[count] = inputData[iColumn][row];
636 count++;
637 }
639 outputData[row] = quartileResult[1] - quartileResult[0];
640 }
641 break;
642 case SET_SMALLEST:
643 for (row = 0; row < rows; row++) {
644 value1 = DBL_MAX;
645 posColIndex = 0;
646 posColumnName[row] = NULL;
647 for (iColumn = 0; iColumn < stat[iStat].sourceColumns; iColumn++) {
648 if (stat[iStat].flags & TOPLIMIT_GIVEN && inputData[iColumn][row] > stat[iStat].topLimit)
649 continue;
650 if (stat[iStat].flags & BOTTOMLIMIT_GIVEN && inputData[iColumn][row] < stat[iStat].bottomLimit)
651 continue;
652 if ((value2 = fabs(inputData[iColumn][row])) < value1) {
653 posColIndex = iColumn;
654 value1 = value2;
655 }
656 }
657 outputData[row] = value1;
658 if (stat[iStat].positionColumn)
659 posColumnName[row] = stat[iStat].sourceColumn[posColIndex];
660 }
661 break;
662 case SET_LARGEST:
663 for (row = 0; row < rows; row++) {
664 value1 = 0;
665 posColIndex = 0;
666 posColumnName[row] = NULL;
667 for (iColumn = 0; iColumn < stat[iStat].sourceColumns; iColumn++) {
668 if (stat[iStat].flags & TOPLIMIT_GIVEN && inputData[iColumn][row] > stat[iStat].topLimit)
669 continue;
670 if (stat[iStat].flags & BOTTOMLIMIT_GIVEN && inputData[iColumn][row] < stat[iStat].bottomLimit)
671 continue;
672 if ((value2 = fabs(inputData[iColumn][row])) > value1) {
673 posColIndex = iColumn;
674 value1 = value2;
675 }
676 }
677 outputData[row] = value1;
678 if (stat[iStat].positionColumn)
679 posColumnName[row] = stat[iStat].sourceColumn[posColIndex];
680 }
681 break;
682 case SET_SPREAD_ARG:
683 for (row = 0; row < rows; row++) {
684 value1 = DBL_MAX;
685 value2 = -DBL_MAX;
686 for (iColumn = 0; iColumn < stat[iStat].sourceColumns; iColumn++) {
687 if (stat[iStat].flags & TOPLIMIT_GIVEN && inputData[iColumn][row] > stat[iStat].topLimit)
688 continue;
689 if (stat[iStat].flags & BOTTOMLIMIT_GIVEN && inputData[iColumn][row] < stat[iStat].bottomLimit)
690 continue;
691 if (inputData[iColumn][row] < value1)
692 value1 = inputData[iColumn][row];
693 if (inputData[iColumn][row] > value2)
694 value2 = inputData[iColumn][row];
695 }
696 outputData[row] = value2 - value1;
697 }
698 break;
699 case SET_PERCENTILE:
700 for (row = 0; row < rows; row++) {
701 for (iColumn = count = value1 = 0; iColumn < stat[iStat].sourceColumns; iColumn++) {
702 if (stat[iStat].flags & TOPLIMIT_GIVEN && inputData[iColumn][row] > stat[iStat].topLimit)
703 continue;
704 if (stat[iStat].flags & BOTTOMLIMIT_GIVEN && inputData[iColumn][row] < stat[iStat].bottomLimit)
705 continue;
706 statWorkArray[count] = inputData[iColumn][row];
707 count++;
708 }
709 outputData[row] = HUGE_VAL;
710 if (count)
712 }
713 break;
714 default:
715 SDDS_Bomb(
"invalid statistic code (accumulation loop)");
716 break;
717 }
718 if (!
SDDS_SetColumn(&outData, SDDS_SET_BY_INDEX, outputData, rows, stat[iStat].resultIndex))
720
721 if (stat[iStat].positionColumn) {
722 if (!
SDDS_SetColumn(&outData, SDDS_SET_BY_INDEX, posColumnName, rows, stat[iStat].positionColumnIndex))
724 }
725 for (iColumn = 0; iColumn < stat[iStat].sourceColumns; iColumn++)
726 free(inputData[iColumn]);
727 free(inputData);
728 inputData = NULL;
729 }
730 free(outputData);
731 outputData = NULL;
732 free(posColumnName);
733 posColumnName = NULL;
734 }
737 }
739 for (iStat = 0; iStat < stats; iStat++) {
740 if (stat[iStat].positionColumn)
741 free(stat[iStat].positionColumn);
742 for (iColumn = 0; iColumn < stat[iStat].sourceColumns; iColumn++)
743 free(stat[iStat].sourceColumn[iColumn]);
744 free(stat[iStat].sourceColumn);
745 }
746 free(request);
747 free(stat);
748 if (statWorkArray)
749 free(statWorkArray);
750
753 exit(EXIT_FAILURE);
754 }
756 exit(EXIT_FAILURE);
757 return EXIT_SUCCESS;
758}
int32_t SDDS_CopyPage(SDDS_DATASET *SDDS_target, SDDS_DATASET *SDDS_source)
int32_t SDDS_SetColumn(SDDS_DATASET *SDDS_dataset, int32_t mode, void *data, int64_t rows,...)
Sets the values for one data column in the current data table of an SDDS dataset.
int32_t SDDS_WritePage(SDDS_DATASET *SDDS_dataset)
Writes the current data table to the output file.
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.
#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.
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 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.
long computeMomentsThreaded(double *mean, double *rms, double *standDev, double *meanAbsoluteDev, double *x, long n, long numThreads)
Computes the mean, RMS, standard deviation, and mean absolute deviation of an array using multiple th...
long replaceFileAndBackUp(char *file, char *replacement)
Replaces a file with a replacement file and creates a backup of the original.
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)
void free_scanargs(SCANNED_ARG **scanned, int argc)
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.