218 {
220 long stats;
222 long requests;
223 SCANNED_ARG *scanned;
225 long i_arg, code, power, iStat, pages, nowarnings = 0;
226 int64_t i, rows, firstRows;
227 char *input, *output;
228 double *inputData, *otherData, indepData, *weight;
229 unsigned long pipeFlags, majorOrderFlag;
230 double decilePoint[2] = {10.0, 90.0}, decileResult[2];
231 double percentilePoint, percentileResult;
232 double percentile;
233 short columnMajorOrder = -1;
234
236 argc =
scanargs(&scanned, argc, argv);
237 if (argc < 2) {
238 bomb(
"too few arguments", USAGE);
239 exit(EXIT_FAILURE);
240 }
241 weight = NULL;
242 input = output = NULL;
243 stat = NULL;
244 request = NULL;
245 stats = requests = pipeFlags = 0;
246 rows = firstRows = i = 0;
247
248 for (i_arg = 1; i_arg < argc; i_arg++) {
249 if (scanned[i_arg].arg_type == OPTION) {
250
251 switch (code =
match_string(scanned[i_arg].list[0], option, N_OPTIONS, 0)) {
252 case SET_COPY:
253 case SET_MINIMA:
254 case SET_MAXIMA:
255 case SET_LARGEST:
256 case SET_SIGNEDLARGEST:
257 case SET_MEANS:
258 case SET_SDS:
259 case SET_SIGMAS:
260 case SET_RMSS:
261 case SET_MEDIAN:
262 case SET_DRANGE:
263 case SET_EXMM_MEAN:
264 if (scanned[i_arg].n_items < 2) {
265 fprintf(stderr, "error: invalid -%s syntax\n", option[code]);
266 exit(EXIT_FAILURE);
267 }
268 requests = addStatRequests(&request, requests, scanned[i_arg].list + 1, scanned[i_arg].n_items - 1, code, 0, 0, NULL, 0, NULL);
269 break;
270 case SET_WMEANS:
271 case SET_WSDS:
272 case SET_WRMSS:
273 case SET_WSIGMAS:
274 if (scanned[i_arg].n_items < 3) {
275 fprintf(stderr, "error: invalid -%s syntax\n", option[code]);
276 exit(EXIT_FAILURE);
277 }
278
279 requests = addStatRequests(&request, requests, scanned[i_arg].list + 1, scanned[i_arg].n_items - 1, code, 0, 0, NULL, 1, NULL);
280 break;
281 case SET_SUMS:
282 if (scanned[i_arg].n_items < 3) {
283 fprintf(stderr, "error: invalid -%s syntax\n", option[code]);
284 exit(EXIT_FAILURE);
285 }
286 if (sscanf(scanned[i_arg].list[1], "%ld", &power) != 1 || power < 1) {
287 fprintf(stderr, "error: invalid -%s syntax--bad power in field %s\n", option[code], scanned[i_arg].list[2]);
288 exit(EXIT_FAILURE);
289 }
290 requests = addStatRequests(&request, requests, scanned[i_arg].list + 2, scanned[i_arg].n_items - 2, code, 0, power, NULL, 0, NULL);
291 break;
292 case SET_PERCENTILE:
293 if (scanned[i_arg].n_items < 3) {
294 fprintf(stderr, "error: invalid -%s syntax\n", option[code]);
295 exit(EXIT_FAILURE);
296 }
297 if (sscanf(scanned[i_arg].list[1], "%lf", &percentile) != 1 || percentile < 0 || percentile > 100) {
298 fprintf(stderr, "error: invalid -%s syntax--bad percentage in field %s\n", option[code], scanned[i_arg].list[1]);
299 exit(EXIT_FAILURE);
300 }
301 requests = addStatRequests(&request, requests, scanned[i_arg].list + 2, scanned[i_arg].n_items - 2, code, percentile, 0, NULL, 0, scanned[i_arg].list[1]);
302 break;
303 case SET_SLOPE:
304 case SET_INTERCEPT:
305 case SET_PMINIMA:
306 case SET_PMAXIMA:
307 case SET_CMINIMA:
308 case SET_CMAXIMA:
309 if (scanned[i_arg].n_items < 3) {
310 fprintf(stderr, "error: invalid -%s syntax\n", option[code]);
311 exit(EXIT_FAILURE);
312 }
313 requests = addStatRequests(&request, requests, scanned[i_arg].list + 2, scanned[i_arg].n_items - 2, code, 0, 0, scanned[i_arg].list[1], 0, NULL);
314 break;
315 case SET_PIPE:
316 if (!
processPipeOption(scanned[i_arg].list + 1, scanned[i_arg].n_items - 1, &pipeFlags))
318 break;
319 case SET_NOWARNINGS:
320 nowarnings = 1;
321 break;
322 case SET_MAJOR_ORDER:
323 majorOrderFlag = 0;
324 scanned[i_arg].n_items--;
325 if (scanned[i_arg].n_items > 0 && (!
scanItemList(&majorOrderFlag, scanned[i_arg].list + 1, &scanned[i_arg].n_items, 0,
"row", -1, NULL, 0, SDDS_ROW_MAJOR_ORDER,
"column", -1, NULL, 0, SDDS_COLUMN_MAJOR_ORDER, NULL)))
326 SDDS_Bomb(
"invalid -majorOrder syntax/values");
327 if (majorOrderFlag & SDDS_COLUMN_MAJOR_ORDER)
328 columnMajorOrder = 1;
329 else if (majorOrderFlag & SDDS_ROW_MAJOR_ORDER)
330 columnMajorOrder = 0;
331 break;
332 default:
333 fprintf(stderr, "error: unknown option '%s' given\n", scanned[i_arg].list[0]);
334 exit(EXIT_FAILURE);
335 break;
336 }
337 } else {
338
339 if (!input)
340 input = scanned[i_arg].list[0];
341 else if (!output)
342 output = scanned[i_arg].list[0];
343 else
345 }
346 }
347
349
350 if (!requests)
352
355
356 pages = 0;
358 pages++;
360 SDDS_Bomb(
"empty data page in input file");
361 if (code == 1) {
362 firstRows = rows;
363 if (!(stat = compileStatDefinitions(&inTable, &stats, request, requests))) {
365 exit(EXIT_FAILURE);
366 }
367 if (!setupOutputFile(&outTable, output, &inTable, stat, stats, rows, columnMajorOrder)) {
370 else
371 fprintf(stderr, "Error setting up output file.\n");
372 exit(EXIT_FAILURE);
373 }
374 } else if (firstRows != rows)
375 SDDS_Bomb(
"inconsistent number of rows in input file");
376 for (iStat = 0; iStat < stats; iStat++) {
377 if (stat[iStat].optionCode == SET_COPY) {
378 if (code == 1 && !(stat[iStat].copy =
SDDS_GetColumn(&inTable, stat[iStat].sourceColumn)))
380 continue;
381 }
382 stat[iStat].copy = NULL;
385 switch (stat[iStat].optionCode) {
386 case SET_MINIMA:
387 if (code == 1)
388 for (i = 0; i < rows; i++)
389 stat[iStat].value1[i] = inputData[i];
390 else
391 for (i = 0; i < rows; i++)
392 if (stat[iStat].value1[i] > inputData[i])
393 stat[iStat].value1[i] = inputData[i];
394 break;
395 case SET_MAXIMA:
396 if (code == 1)
397 for (i = 0; i < rows; i++)
398 stat[iStat].value1[i] = inputData[i];
399 else
400 for (i = 0; i < rows; i++)
401 if (stat[iStat].value1[i] < inputData[i])
402 stat[iStat].value1[i] = inputData[i];
403 break;
404 case SET_CMINIMA:
405
408 if (code == 1)
409 for (i = 0; i < rows; i++) {
410 stat[iStat].value2[i] = inputData[i];
411 stat[iStat].value1[i] = otherData[i];
412 }
413 else
414 for (i = 0; i < rows; i++)
415 if (stat[iStat].value2[i] > inputData[i]) {
416 stat[iStat].value2[i] = inputData[i];
417 stat[iStat].value1[i] = otherData[i];
418 }
419 free(otherData);
420 break;
421 case SET_CMAXIMA:
422
425 if (code == 1)
426 for (i = 0; i < rows; i++) {
427 stat[iStat].value2[i] = inputData[i];
428 stat[iStat].value1[i] = otherData[i];
429 }
430 else
431 for (i = 0; i < rows; i++)
432 if (stat[iStat].value2[i] < inputData[i]) {
433 stat[iStat].value2[i] = inputData[i];
434 stat[iStat].value1[i] = otherData[i];
435 }
436 free(otherData);
437 break;
438 case SET_LARGEST:
439 if (code == 1)
440 for (i = 0; i < rows; i++)
441 stat[iStat].value1[i] = fabs(inputData[i]);
442 else
443 for (i = 0; i < rows; i++)
444 if (stat[iStat].value1[i] < fabs(inputData[i]))
445 stat[iStat].value1[i] = fabs(inputData[i]);
446 break;
447 case SET_SIGNEDLARGEST:
448 if (code == 1)
449 for (i = 0; i < rows; i++)
450 stat[iStat].value1[i] = inputData[i];
451 else
452 for (i = 0; i < rows; i++)
453 if (fabs(stat[iStat].value1[i]) < fabs(inputData[i]))
454 stat[iStat].value1[i] = inputData[i];
455 break;
456 case SET_MEANS:
457 if (code == 1)
458 for (i = 0; i < rows; i++)
459 stat[iStat].value1[i] = inputData[i];
460 else
461 for (i = 0; i < rows; i++)
462 stat[iStat].value1[i] += inputData[i];
463 break;
464 case SET_WMEANS:
467 for (i = 0; i < rows; i++) {
468 stat[iStat].sumWeight[i] += weight[i];
469 stat[iStat].value1[i] += inputData[i] * weight[i];
470 }
471 free(weight);
472 break;
473 case SET_SDS:
474 case SET_SIGMAS:
475 if (code == 1)
476 for (i = 0; i < rows; i++) {
477 stat[iStat].value1[i] = inputData[i];
478 stat[iStat].value2[i] = inputData[i] * inputData[i];
479 }
480 else
481 for (i = 0; i < rows; i++) {
482 stat[iStat].value1[i] += inputData[i];
483 stat[iStat].value2[i] += inputData[i] * inputData[i];
484 }
485 break;
486 case SET_WSDS:
487 case SET_WSIGMAS:
490 for (i = 0; i < rows; i++) {
491 stat[iStat].sumWeight[i] += weight[i];
492 stat[iStat].value1[i] += inputData[i] * weight[i];
493 stat[iStat].value2[i] += inputData[i] * inputData[i] * weight[i];
494 }
495 free(weight);
496 break;
497 case SET_RMSS:
498 if (code == 1)
499 for (i = 0; i < rows; i++)
500 stat[iStat].value1[i] = inputData[i] * inputData[i];
501 else
502 for (i = 0; i < rows; i++)
503 stat[iStat].value1[i] += inputData[i] * inputData[i];
504 break;
505 case SET_WRMSS:
508 for (i = 0; i < rows; i++) {
509 stat[iStat].sumWeight[i] += weight[i];
510 stat[iStat].value1[i] += inputData[i] * inputData[i] * weight[i];
511 }
512 free(weight);
513 break;
514 case SET_SUMS:
515 if (code == 1)
516 for (i = 0; i < rows; i++)
517 stat[iStat].value1[i] =
ipow(inputData[i], stat[iStat].sumPower);
518 else
519 for (i = 0; i < rows; i++)
520 stat[iStat].value1[i] +=
ipow(inputData[i], stat[iStat].sumPower);
521 break;
522 case SET_PMINIMA:
524 fprintf(stderr, "error: unable to get value of parameter %s\n", stat[iStat].functionOf);
526 }
527 if (code == 1)
528 for (i = 0; i < rows; i++) {
529 stat[iStat].value2[i] = inputData[i];
530 stat[iStat].value1[i] = indepData;
531 }
532 else
533 for (i = 0; i < rows; i++) {
534 if (stat[iStat].value2[i] > inputData[i]) {
535 stat[iStat].value2[i] = inputData[i];
536 stat[iStat].value1[i] = indepData;
537 }
538 }
539 break;
540 case SET_PMAXIMA:
542 fprintf(stderr, "error: unable to get value of parameter %s\n", stat[iStat].functionOf);
544 }
545 if (code == 1)
546 for (i = 0; i < rows; i++) {
547 stat[iStat].value2[i] = inputData[i];
548 stat[iStat].value1[i] = indepData;
549 }
550 else
551 for (i = 0; i < rows; i++) {
552 if (stat[iStat].value2[i] < inputData[i]) {
553 stat[iStat].value2[i] = inputData[i];
554 stat[iStat].value1[i] = indepData;
555 }
556 }
557 break;
558 case SET_SLOPE:
559 case SET_INTERCEPT:
561 fprintf(stderr, "error: unable to get value of parameter %s\n", stat[iStat].functionOf);
563 }
564
565
566
567
568
569
570 if (code == 1)
571 for (i = 0; i < rows; i++) {
572 stat[iStat].value1[i] = indepData;
573 stat[iStat].value2[i] = indepData * indepData;
574 stat[iStat].value3[i] = inputData[i];
575 stat[iStat].value4[i] = indepData * inputData[i];
576 }
577 else
578 for (i = 0; i < rows; i++) {
579 stat[iStat].value1[i] += indepData;
580 stat[iStat].value2[i] += indepData * indepData;
581 stat[iStat].value3[i] += inputData[i];
582 stat[iStat].value4[i] += indepData * inputData[i];
583 }
584 break;
585 case SET_MEDIAN:
586 case SET_DRANGE:
587 case SET_PERCENTILE:
588 case SET_EXMM_MEAN:
589 if (code == 1)
590 for (i = 0; i < rows; i++) {
591 stat[iStat].array[i] =
tmalloc(
sizeof(*stat[iStat].array[i]));
592 stat[iStat].array[i][pages - 1] = inputData[i];
593 }
594 else {
595 for (i = 0; i < rows; i++) {
596 stat[iStat].array[i] =
SDDS_Realloc(stat[iStat].array[i],
sizeof(*stat[iStat].array[i]) * pages);
597 stat[iStat].array[i][pages - 1] = inputData[i];
598 }
599 }
600 break;
601 default:
602 SDDS_Bomb(
"invalid statistic code (accumulation loop)");
603 break;
604 }
605 free(inputData);
606 }
607 }
608 if (pages == 0)
610 for (iStat = 0; iStat < stats; iStat++) {
611 switch (stat[iStat].optionCode) {
612 case SET_COPY:
613 case SET_MINIMA:
614 case SET_MAXIMA:
615 case SET_PMINIMA:
616 case SET_PMAXIMA:
617 case SET_CMINIMA:
618 case SET_CMAXIMA:
619 case SET_LARGEST:
620 case SET_SIGNEDLARGEST:
621 case SET_SUMS:
622 break;
623 case SET_MEANS:
624 for (i = 0; i < rows; i++)
625 stat[iStat].value1[i] /= pages;
626 break;
627 case SET_WMEANS:
628 for (i = 0; i < rows; i++)
629 if (stat[iStat].sumWeight[i])
630 stat[iStat].value1[i] /= stat[iStat].sumWeight[i];
631 else {
632 if (!nowarnings)
633 fprintf(stderr, "warning: the total weight for the %" PRId64 "th row of %s is zero.\n", i + 1, stat[iStat].sourceColumn);
634 stat[iStat].value1[i] = DBL_MAX;
635 }
636 break;
637 case SET_SDS:
638 if (pages < 2)
639 stat[iStat].value1[i] = DBL_MAX;
640 else
641 for (i = 0; i < rows; i++) {
642 double tmp1;
643 if ((tmp1 = stat[iStat].value2[i] / pages - sqr(stat[iStat].value1[i] / pages)) <= 0)
644 stat[iStat].value1[i] = 0;
645 else
646 stat[iStat].value1[i] = sqrt(tmp1 * pages / (pages - 1.0));
647 }
648 break;
649 case SET_WSDS:
650 if (pages < 2)
651 stat[iStat].value1[i] = DBL_MAX;
652 else
653 for (i = 0; i < rows; i++) {
654 double tmp1;
655 if (stat[iStat].sumWeight[i]) {
656 if ((tmp1 = stat[iStat].value2[i] / stat[iStat].sumWeight[i] - sqr(stat[iStat].value1[i] / stat[iStat].sumWeight[i])) <= 0)
657 stat[iStat].value1[i] = 0;
658 else
659 stat[iStat].value1[i] = sqrt(tmp1 * pages / (pages - 1.0));
660 } else {
661 if (!nowarnings)
662 fprintf(stderr, "Warning, the total weight for the %" PRId64 "th row of %s is zero.\n", i + 1, stat[iStat].sourceColumn);
663 stat[iStat].value1[i] = DBL_MAX;
664 }
665 }
666 break;
667 case SET_SIGMAS:
668 if (pages < 2)
669 stat[iStat].value1[i] = DBL_MAX;
670 else
671 for (i = 0; i < rows; i++) {
672 double tmp1;
673 if ((tmp1 = stat[iStat].value2[i] / pages - sqr(stat[iStat].value1[i] / pages)) <= 0)
674 stat[iStat].value1[i] = 0;
675 else
676 stat[iStat].value1[i] = sqrt(tmp1 / (pages - 1.0));
677 }
678 break;
679 case SET_WSIGMAS:
680 if (pages < 2)
681 stat[iStat].value1[i] = DBL_MAX;
682 else
683 for (i = 0; i < rows; i++) {
684 double tmp1;
685 if (stat[iStat].sumWeight[i]) {
686 if ((tmp1 = stat[iStat].value2[i] / stat[iStat].sumWeight[i] - sqr(stat[iStat].value1[i] / stat[iStat].sumWeight[i])) <= 0)
687 stat[iStat].value1[i] = 0;
688 else
689 stat[iStat].value1[i] = sqrt(tmp1 / (pages - 1.0));
690 } else {
691 if (!nowarnings)
692 fprintf(stderr, "Warning, the total weight for the %" PRId64 "th row of %s is zero.\n", i + 1, stat[iStat].sourceColumn);
693 stat[iStat].value1[i] = DBL_MAX;
694 }
695 }
696 break;
697 case SET_RMSS:
698 for (i = 0; i < rows; i++)
699 stat[iStat].value1[i] = sqrt(stat[iStat].value1[i] / pages);
700 break;
701 case SET_WRMSS:
702 for (i = 0; i < rows; i++) {
703 if (stat[iStat].sumWeight[i])
704 stat[iStat].value1[i] = sqrt(stat[iStat].value1[i] / stat[iStat].sumWeight[i]);
705 else {
706 if (!nowarnings)
707 fprintf(stderr, "Warning, the total weight for the %" PRId64 "th row of %s is zero.\n", i + 1, stat[iStat].sourceColumn);
708 stat[iStat].value1[i] = DBL_MAX;
709 }
710 }
711 break;
712 case SET_SLOPE:
713 for (i = 0; i < rows; i++) {
714 double D;
715 D = pages * stat[iStat].value2[i] - stat[iStat].value1[i] * stat[iStat].value1[i];
716 stat[iStat].value1[i] = (pages * stat[iStat].value4[i] - stat[iStat].value1[i] * stat[iStat].value3[i]) / D;
717 }
718 break;
719 case SET_INTERCEPT:
720 for (i = 0; i < rows; i++) {
721 double D;
722 D = pages * stat[iStat].value2[i] - stat[iStat].value1[i] * stat[iStat].value1[i];
723 stat[iStat].value1[i] = (stat[iStat].value2[i] * stat[iStat].value3[i] - stat[iStat].value1[i] * stat[iStat].value4[i]) / D;
724 }
725 break;
726 case SET_MEDIAN:
727 for (i = 0; i < rows; i++) {
728 compute_median(&stat[iStat].value1[i], stat[iStat].array[i], pages);
729 }
730 break;
731 case SET_DRANGE:
732 for (i = 0; i < rows; i++) {
734 stat[iStat].value1[i] = 0;
735 else
736 stat[iStat].value1[i] = decileResult[1] - decileResult[0];
737 }
738 break;
739 case SET_PERCENTILE:
740 percentilePoint = stat[iStat].percentile;
741 for (i = 0; i < rows; i++) {
742 if (!
compute_percentiles(&percentileResult, &percentilePoint, 1, stat[iStat].array[i], pages))
743 stat[iStat].value1[i] = 0;
744 else
745 stat[iStat].value1[i] = percentileResult;
746 }
747 break;
748 case SET_EXMM_MEAN:
749 for (i = 0; i < rows; i++)
750 if (!compute_mean_exclude_min_max(&(stat[iStat].value1[i]), stat[iStat].array[i], pages))
751 stat[iStat].value1[i] = 0;
752 break;
753 default:
754 SDDS_Bomb(
"invalid statistic code (final loop)");
755 break;
756 }
757 if (stat[iStat].optionCode == SET_COPY) {
758 if (!
SDDS_SetColumn(&outTable, SDDS_SET_BY_NAME, stat[iStat].copy, rows, stat[iStat].resultColumn)) {
759 fprintf(stderr, "error setting column values for column %s\n", stat[iStat].resultColumn);
761 }
762 }
else if (!
SDDS_SetColumnFromDoubles(&outTable, SDDS_SET_BY_NAME, stat[iStat].value1, rows, stat[iStat].resultColumn)) {
763 fprintf(stderr, "error setting column values for column %s\n", stat[iStat].resultColumn);
765 }
766 if (stat[iStat].value1)
767 free(stat[iStat].value1);
768 if (stat[iStat].value2)
769 free(stat[iStat].value2);
770 if (stat[iStat].value3)
771 free(stat[iStat].value3);
772 if (stat[iStat].value4)
773 free(stat[iStat].value4);
774 if (stat[iStat].copy)
775 free(stat[iStat].copy);
776 if (stat[iStat].array) {
777 for (i = 0; i < rows; i++) {
778 free(stat[iStat].array[i]);
779 }
780 free(stat[iStat].array);
781 }
782 if (stat[iStat].sumWeight)
783 free(stat[iStat].sumWeight);
784 free(stat[iStat].sourceColumn);
785 free(stat[iStat].resultColumn);
786 stat[iStat].value1 = stat[iStat].value2 = stat[iStat].value3 = stat[iStat].value4 = NULL;
787 stat[iStat].copy = NULL;
788 stat[iStat].array = NULL;
789 }
790 free(stat);
793 return EXIT_SUCCESS;
794}
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_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_RegisterProgramName(const char *name)
Registers the executable program name for use in error messages.
int32_t SDDS_NumberOfErrors()
Retrieves the number of errors recorded by SDDS library routines.
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.
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.