SDDS ToolKit Programs and Libraries for C and Python
Loading...
Searching...
No Matches
sddsenvelope.c
Go to the documentation of this file.
1/**
2 * @file sddsenvelope.c
3 * @brief Combine data from SDDS pages to create a new file with computed statistics.
4 *
5 * @details
6 * This program processes SDDS (Self Describing Data Set) files, performing statistical computations
7 * such as maximum, minimum, mean, and others across pages of data. The resulting statistics are stored
8 * in an output SDDS file for further analysis.
9 *
10 * @section Usage
11 * ```
12 * sddsenvelope [<inputfile>] [<outputfile>]
13 * [-pipe=[input][,output]]
14 * [-nowarnings]
15 * [-maximum=<column-names>]
16 * [-minimum=<column-names>]
17 * [-cmaximum=<indep-column>,<column-names>]
18 * [-cminimum=<indep-column>,<column-names>]
19 * [-pmaximum=<indep-parameter>,<column-names>]
20 * [-pminimum=<indep-parameter>,<column-names>]
21 * [-largest=<column-names>]
22 * [-signedLargest=<column-names>]
23 * [-mean=<column-names>]
24 * [-sum=<power>,<column-names>]
25 * [-median=<column-names>]
26 * [-decilerange=<column-names>]
27 * [-percentile=<percentage>,<column-names>]
28 * [-standarddeviation=<column-names>]
29 * [-rms=<column-names>]
30 * [-sigma=<column-names>]
31 * [-slope=<indep-parameter>,<column-names>]
32 * [-intercept=<indep-parameter>,<column-names>]
33 * [-wmean=<weightColumn>,<columnNames>]
34 * [-wstandarddeviation=<weightColumn>,<columnNames>]
35 * [-wrms=<weightColumn>,<columnNames>]
36 * [-wsigma=<weightColumn>,<columnNames>]
37 * [-majorOrder=row|column]
38 * [-threads=<number>]
39 * ```
40 *
41 * @section Options
42 * | Option | Description |
43 * |------------------------------------|--------------------------------------------------------------------------|
44 * | `-pipe` | Use pipe for input/output. |
45 * | `-nowarnings` | Suppress warnings. |
46 * | `-maximum` | Compute maximum values for specified columns. |
47 * | `-minimum` | Compute minimum values for specified columns. |
48 * | `-cmaximum` | Compute conditional maximum based on an independent column. |
49 * | `-cminimum` | Compute conditional minimum based on an independent column. |
50 * | `-pmaximum` | Parameter-based maximum. |
51 * | `-pminimum` | Parameter-based minimum. |
52 * | `-largest` | Compute largest absolute values. |
53 * | `-signedLargest` | Compute largest signed values. |
54 * | `-mean` | Compute mean values. |
55 * | `-sum` | Compute sum raised to a specified power. |
56 * | `-median` | Compute median values. |
57 * | `-decilerange` | Compute decile range. |
58 * | `-percentile` | Compute specified percentile. |
59 * | `-standarddeviation` | Compute standard deviation. |
60 * | `-rms` | Compute root mean square (RMS) values. |
61 * | `-sigma` | Compute sigma values. |
62 * | `-slope` | Compute slope for linear fit. |
63 * | `-intercept` | Compute intercept for linear fit. |
64 * | `-wmean` | Compute weighted mean. |
65 * | `-wstandarddeviation` | Compute weighted standard deviation. |
66 * | `-wrms` | Compute weighted RMS. |
67 * | `-wsigma` | Compute weighted sigma. |
68 * | `-majorOrder` | Set data ordering for output file. |
69 * | `-threads` | Number of threads for row-wise statistic calculations. |
70 *
71 * @subsection Incompatibilities
72 * - Only one of the following options may be specified:
73 * - `-mean`
74 * - `-median`
75 * - `-percentile`
76 * - `-standarddeviation`
77 *
78 * @subsection SR Specific Requirements
79 * - `-percentile` requires a percentage value between 0 and 100.
80 *
81 * @authors
82 * M. Borland, C. Saunders, R. Soliday, H. Shang
83 *
84 * @license
85 * This file is distributed under the terms of the Software License Agreement
86 * found in the file LICENSE included with this distribution.
87 */
88
89#include "mdb.h"
90#include "scan.h"
91#include "SDDS.h"
92#include <ctype.h>
93
94/* Enumeration for option types */
95enum option_type {
96 SET_COPY,
97 SET_MAXIMA,
98 SET_MINIMA,
99 SET_MEANS,
100 SET_SDS,
101 SET_RMSS,
102 SET_SUMS,
103 SET_SLOPE,
104 SET_INTERCEPT,
105 SET_PIPE,
106 SET_SIGMAS,
107 SET_MEDIAN,
108 SET_DRANGE,
109 SET_WMEANS,
110 SET_WSDS,
111 SET_WRMSS,
112 SET_WSIGMAS,
113 SET_NOWARNINGS,
114 SET_LARGEST,
115 SET_PERCENTILE,
116 SET_SIGNEDLARGEST,
117 SET_PMAXIMA,
118 SET_PMINIMA,
119 SET_MAJOR_ORDER,
120 SET_EXMM_MEAN,
121 SET_CMAXIMA,
122 SET_CMINIMA,
123 SET_THREADS,
124 N_OPTIONS
125};
126
127char *option[N_OPTIONS] = {
128 "copy",
129 "maximum",
130 "minimum",
131 "mean",
132 "standarddeviations",
133 "rms",
134 "sum",
135 "slope",
136 "intercept",
137 "pipe",
138 "sigmas",
139 "median",
140 "decilerange",
141 "wmean",
142 "wstandarddeviations",
143 "wrms",
144 "wsigma",
145 "nowarnings",
146 "largest",
147 "percentile",
148 "signedlargest",
149 "pmaximum",
150 "pminimum",
151 "majorOrder",
152 "exmmMean",
153 "cmaximum",
154 "cminimum",
155 "threads",
156};
157
158char *optionSuffix[N_OPTIONS] = {
159 "",
160 "Max",
161 "Min",
162 "Mean",
163 "StDev",
164 "Rms",
165 "Sum",
166 "Slope",
167 "Intercept",
168 "",
169 "Sigma",
170 "Median",
171 "DRange",
172 "WMean",
173 "WStDev",
174 "WRms",
175 "WSigma",
176 "",
177 "Largest",
178 "Percentile",
179 "SignedLargest",
180 "PMaximum",
181 "PMinimum",
182 "",
183 "ExmmMean",
184 "CMaximum",
185 "CMinimum",
186 "",
187};
188
189/* this structure stores a command-line request for statistics computation */
190/* columnName may contain wildcards */
191typedef struct
192{
193 char *columnName;
194 char *weightColumnName;
195 long optionCode, sumPower;
196 double percentile;
197 char *percentileString;
198 char *functionOf;
199} STAT_REQUEST;
200
201/* this structure stores data necessary for accessing/creating SDDS columns and
202 * for computing a statistic
203 */
204typedef struct
205{
206 char *sourceColumn, *weightColumn, *resultColumn, *functionOf;
207 long optionCode, resultIndex, sumPower;
208 double percentile;
209 char *percentileString;
210 /* these store intermediate values during processing */
211 void *copy;
212 double *value1, *value2, *value3, *value4;
213 double **array;
214 double *sumWeight;
215} STAT_DEFINITION;
216
217long addStatRequests(STAT_REQUEST **statRequest, long requests, char **item, long items, long code, double percentile, long power, char *functionOf, long weighted, char *percentileString);
218/*weighted=0, no weighted column; else, weighted statistic, the weight factor is given by
219 weightedColumn*/
220STAT_DEFINITION *compileStatDefinitions(SDDS_DATASET *inTable, long *stats, STAT_REQUEST *request, long requests);
221long setupOutputFile(SDDS_DATASET *outTable, char *output, SDDS_DATASET *inTable, STAT_DEFINITION *stat, long stats, int64_t rows, short columnMajorOrder);
222int compute_mean_exclude_min_max(double *value, double *data, long n);
223
224static char *USAGE = "sddsenvelope [<input>] [<output>] [options]\n"
225 " [-pipe=[input][,output]]\n"
226 " [-nowarnings]\n"
227 " [-maximum=<column-names>]\n"
228 " [-minimum=<column-names>]\n"
229 " [-cmaximum=<indep-column>,<column-names>]\n"
230 " [-cminimum=<indep-column>,<column-names>]\n"
231 " [-pmaximum=<indep-parameter>,<column-names>]\n"
232 " [-pminimum=<indep-parameter>,<column-names>]\n"
233 " [-largest=<column-names>]\n"
234 " [-signedLargest=<column-names>]\n"
235 " [-mean=<column-names>]\n"
236 " [-sum=<power>,<column-names>]\n"
237 " [-median=<column-names>]\n"
238 " [-decilerange=<column-names>]\n"
239 " [-percentile=<percentage>,<column-names>]\n"
240 " [-standarddeviation=<column-names>]\n"
241 " [-rms=<column-names>]\n"
242 " [-sigma=<column-names>]\n"
243 " [-slope=<indep-parameter>,<column-names>]\n"
244 " [-intercept=<indep-parameter>,<column-names>]\n"
245 " [-wmean=<weightColumn>,<columnNames>]\n"
246 " [-wstandarddeviation=<weightColumn>,<columnNames>]\n"
247 " [-wrms=<weightColumn>,<columnNames>]\n"
248 " [-wsigma=<weightColumn>,<columnNames>]\n"
249 " [-majorOrder=row|column]\n"
250 " [-threads=<number>]\n"
251 "Options:\n"
252 " -copy=<column-names> Copy specified columns.\n"
253 " -pipe=[input][,output] Use pipe for input/output.\n"
254 " -nowarnings Suppress warnings.\n"
255 " -maximum=<column-names> Compute maximum values.\n"
256 " -minimum=<column-names> Compute minimum values.\n"
257 " -cmaximum=<indep-column>,<column-names> Conditional maximum based on an independent column.\n"
258 " -cminimum=<indep-column>,<column-names> Conditional minimum based on an independent column.\n"
259 " -pmaximum=<indep-parameter>,<column-names> Parameter-based maximum.\n"
260 " -pminimum=<indep-parameter>,<column-names> Parameter-based minimum.\n"
261 " -largest=<column-names> Compute the largest absolute values.\n"
262 " -signedLargest=<column-names> Compute the largest signed values.\n"
263 " -mean=<column-names> Compute mean values.\n"
264 " -sum=<power>,<column-names> Compute sum with power.\n"
265 " -median=<column-names> Compute median values.\n"
266 " -decilerange=<column-names> Compute decile range.\n"
267 " -percentile=<percentage>,<column-names> Compute specified percentile.\n"
268 " -standarddeviation=<column-names> Compute standard deviations.\n"
269 " -rms=<column-names> Compute RMS values.\n"
270 " -sigma=<column-names> Compute sigma values.\n"
271 " -slope=<indep-parameter>,<column-names> Compute slope for linear fit.\n"
272 " -intercept=<indep-parameter>,<column-names> Compute intercept for linear fit.\n"
273 " -wmean=<weightColumn>,<columnNames> Compute weighted mean.\n"
274 " -wstandarddeviation=<weightColumn>,<columnNames> Compute weighted standard deviation.\n"
275 " -wrms=<weightColumn>,<columnNames> Compute weighted RMS.\n"
276 " -wsigma=<weightColumn>,<columnNames> Compute weighted sigma.\n"
277 " -majorOrder=row|column Set major order.\n\n"
278 " -threads=<number> Number of threads for row-wise statistics. The default is 1.\n"
279 "Processes pages from <input> to produce <output> with\n"
280 "one page containing the specified quantities across pages\n"
281 "for each row of the specified columns.\n"
282 "Program by Michael Borland. (" __DATE__ " " __TIME__ ", SVN revision: " SVN_VERSION ")";
283
284int main(int argc, char **argv) {
285 STAT_DEFINITION *stat;
286 long stats;
287 STAT_REQUEST *request;
288 long requests;
289 SCANNED_ARG *scanned; /* structure for scanned arguments */
290 SDDS_DATASET inTable, outTable;
291 long i_arg, code, power, iStat, pages, nowarnings = 0;
292 int64_t i, rows, firstRows;
293 char *input, *output;
294 double *inputData, *otherData, indepData, *weight;
295 unsigned long pipeFlags, majorOrderFlag;
296 double decilePoint[2] = {10.0, 90.0};
297 double percentilePoint;
298 double percentile;
299 short columnMajorOrder = -1;
300 int threads = 1;
301
302 SDDS_RegisterProgramName(argv[0]);
303 argc = scanargs(&scanned, argc, argv);
304 if (argc < 2) {
305 bomb("too few arguments", USAGE);
306 exit(EXIT_FAILURE);
307 }
308 weight = NULL;
309 input = output = NULL;
310 stat = NULL;
311 request = NULL;
312 stats = requests = pipeFlags = 0;
313 rows = firstRows = i = 0;
314
315 for (i_arg = 1; i_arg < argc; i_arg++) {
316 if (scanned[i_arg].arg_type == OPTION) {
317 /* process options here */
318 switch (code = match_string(scanned[i_arg].list[0], option, N_OPTIONS, 0)) {
319 case SET_COPY:
320 case SET_MINIMA:
321 case SET_MAXIMA:
322 case SET_LARGEST:
323 case SET_SIGNEDLARGEST:
324 case SET_MEANS:
325 case SET_SDS:
326 case SET_SIGMAS:
327 case SET_RMSS:
328 case SET_MEDIAN:
329 case SET_DRANGE:
330 case SET_EXMM_MEAN:
331 if (scanned[i_arg].n_items < 2) {
332 fprintf(stderr, "error: invalid -%s syntax\n", option[code]);
333 exit(EXIT_FAILURE);
334 }
335 requests = addStatRequests(&request, requests, scanned[i_arg].list + 1, scanned[i_arg].n_items - 1, code, 0, 0, NULL, 0, NULL);
336 break;
337 case SET_WMEANS:
338 case SET_WSDS:
339 case SET_WRMSS:
340 case SET_WSIGMAS:
341 if (scanned[i_arg].n_items < 3) {
342 fprintf(stderr, "error: invalid -%s syntax\n", option[code]);
343 exit(EXIT_FAILURE);
344 }
345 /*note here, items=scanned[i_arg].n_items-2, because the weightedColumn should be excluded */
346 requests = addStatRequests(&request, requests, scanned[i_arg].list + 1, scanned[i_arg].n_items - 1, code, 0, 0, NULL, 1, NULL);
347 break;
348 case SET_SUMS:
349 if (scanned[i_arg].n_items < 3) {
350 fprintf(stderr, "error: invalid -%s syntax\n", option[code]);
351 exit(EXIT_FAILURE);
352 }
353 if (sscanf(scanned[i_arg].list[1], "%ld", &power) != 1 || power < 1) {
354 fprintf(stderr, "error: invalid -%s syntax--bad power in field %s\n", option[code], scanned[i_arg].list[2]);
355 exit(EXIT_FAILURE);
356 }
357 requests = addStatRequests(&request, requests, scanned[i_arg].list + 2, scanned[i_arg].n_items - 2, code, 0, power, NULL, 0, NULL);
358 break;
359 case SET_PERCENTILE:
360 if (scanned[i_arg].n_items < 3) {
361 fprintf(stderr, "error: invalid -%s syntax\n", option[code]);
362 exit(EXIT_FAILURE);
363 }
364 if (sscanf(scanned[i_arg].list[1], "%lf", &percentile) != 1 || percentile < 0 || percentile > 100) {
365 fprintf(stderr, "error: invalid -%s syntax--bad percentage in field %s\n", option[code], scanned[i_arg].list[1]);
366 exit(EXIT_FAILURE);
367 }
368 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]);
369 break;
370 case SET_SLOPE:
371 case SET_INTERCEPT:
372 case SET_PMINIMA:
373 case SET_PMAXIMA:
374 case SET_CMINIMA:
375 case SET_CMAXIMA:
376 if (scanned[i_arg].n_items < 3) {
377 fprintf(stderr, "error: invalid -%s syntax\n", option[code]);
378 exit(EXIT_FAILURE);
379 }
380 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);
381 break;
382 case SET_PIPE:
383 if (!processPipeOption(scanned[i_arg].list + 1, scanned[i_arg].n_items - 1, &pipeFlags))
384 SDDS_Bomb("invalid -pipe syntax");
385 break;
386 case SET_NOWARNINGS:
387 nowarnings = 1;
388 break;
389 case SET_MAJOR_ORDER:
390 majorOrderFlag = 0;
391 scanned[i_arg].n_items--;
392 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)))
393 SDDS_Bomb("invalid -majorOrder syntax/values");
394 if (majorOrderFlag & SDDS_COLUMN_MAJOR_ORDER)
395 columnMajorOrder = 1;
396 else if (majorOrderFlag & SDDS_ROW_MAJOR_ORDER)
397 columnMajorOrder = 0;
398 break;
399 case SET_THREADS:
400 if (scanned[i_arg].n_items != 2 ||
401 sscanf(scanned[i_arg].list[1], "%d", &threads) != 1 ||
402 threads < 1)
403 SDDS_Bomb("invalid -threads syntax");
404 break;
405 default:
406 fprintf(stderr, "error: unknown option '%s' given\n", scanned[i_arg].list[0]);
407 exit(EXIT_FAILURE);
408 break;
409 }
410 } else {
411 /* argument is filename */
412 if (!input)
413 input = scanned[i_arg].list[0];
414 else if (!output)
415 output = scanned[i_arg].list[0];
416 else
417 SDDS_Bomb("too many filenames seen");
418 }
419 }
420
421 processFilenames("sddsenvelope", &input, &output, pipeFlags, 0, NULL);
422
423 if (!requests)
424 SDDS_Bomb("no statistics requested");
425
426 if (!SDDS_InitializeInput(&inTable, input))
427 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
428
429 pages = 0;
430 while ((code = SDDS_ReadPage(&inTable)) > 0) {
431 pages++;
432 if (!(rows = SDDS_CountRowsOfInterest(&inTable)))
433 SDDS_Bomb("empty data page in input file");
434 if (code == 1) {
435 firstRows = rows;
436 if (!(stat = compileStatDefinitions(&inTable, &stats, request, requests))) {
437 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
438 exit(EXIT_FAILURE);
439 }
440 if (!setupOutputFile(&outTable, output, &inTable, stat, stats, rows, columnMajorOrder)) {
441 if (SDDS_NumberOfErrors())
442 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
443 else
444 fprintf(stderr, "Error setting up output file.\n");
445 exit(EXIT_FAILURE);
446 }
447 } else if (firstRows != rows)
448 SDDS_Bomb("inconsistent number of rows in input file");
449 for (iStat = 0; iStat < stats; iStat++) {
450 if (stat[iStat].optionCode == SET_COPY) {
451 if (code == 1 && !(stat[iStat].copy = SDDS_GetColumn(&inTable, stat[iStat].sourceColumn)))
452 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
453 continue;
454 }
455 stat[iStat].copy = NULL;
456 if (!(inputData = SDDS_GetColumnInDoubles(&inTable, stat[iStat].sourceColumn)))
457 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
458 switch (stat[iStat].optionCode) {
459 case SET_MINIMA:
460 if (code == 1)
461#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
462 for (i = 0; i < rows; i++)
463 stat[iStat].value1[i] = inputData[i];
464 else
465#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
466 for (i = 0; i < rows; i++)
467 if (stat[iStat].value1[i] > inputData[i])
468 stat[iStat].value1[i] = inputData[i];
469 break;
470 case SET_MAXIMA:
471 if (code == 1)
472#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
473 for (i = 0; i < rows; i++)
474 stat[iStat].value1[i] = inputData[i];
475 else
476#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
477 for (i = 0; i < rows; i++)
478 if (stat[iStat].value1[i] < inputData[i])
479 stat[iStat].value1[i] = inputData[i];
480 break;
481 case SET_CMINIMA:
482 /* Value from another column corresponding to minimum value in main column */
483 if (!(otherData = SDDS_GetColumnInDoubles(&inTable, stat[iStat].functionOf)))
484 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
485 if (code == 1)
486#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
487 for (i = 0; i < rows; i++) {
488 stat[iStat].value2[i] = inputData[i];
489 stat[iStat].value1[i] = otherData[i];
490 }
491 else
492#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
493 for (i = 0; i < rows; i++)
494 if (stat[iStat].value2[i] > inputData[i]) {
495 stat[iStat].value2[i] = inputData[i];
496 stat[iStat].value1[i] = otherData[i];
497 }
498 free(otherData);
499 break;
500 case SET_CMAXIMA:
501 /* Value from another column corresponding to maximum value in main column */
502 if (!(otherData = SDDS_GetColumnInDoubles(&inTable, stat[iStat].functionOf)))
503 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
504 if (code == 1)
505#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
506 for (i = 0; i < rows; i++) {
507 stat[iStat].value2[i] = inputData[i];
508 stat[iStat].value1[i] = otherData[i];
509 }
510 else
511#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
512 for (i = 0; i < rows; i++)
513 if (stat[iStat].value2[i] < inputData[i]) {
514 stat[iStat].value2[i] = inputData[i];
515 stat[iStat].value1[i] = otherData[i];
516 }
517 free(otherData);
518 break;
519 case SET_LARGEST:
520 if (code == 1)
521#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
522 for (i = 0; i < rows; i++)
523 stat[iStat].value1[i] = fabs(inputData[i]);
524 else
525#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
526 for (i = 0; i < rows; i++)
527 if (stat[iStat].value1[i] < fabs(inputData[i]))
528 stat[iStat].value1[i] = fabs(inputData[i]);
529 break;
530 case SET_SIGNEDLARGEST:
531 if (code == 1)
532#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
533 for (i = 0; i < rows; i++)
534 stat[iStat].value1[i] = inputData[i];
535 else
536#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
537 for (i = 0; i < rows; i++)
538 if (fabs(stat[iStat].value1[i]) < fabs(inputData[i]))
539 stat[iStat].value1[i] = inputData[i];
540 break;
541 case SET_MEANS:
542 if (code == 1)
543#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
544 for (i = 0; i < rows; i++)
545 stat[iStat].value1[i] = inputData[i];
546 else
547#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
548 for (i = 0; i < rows; i++)
549 stat[iStat].value1[i] += inputData[i];
550 break;
551 case SET_WMEANS:
552 if (!(weight = SDDS_GetColumnInDoubles(&inTable, stat[iStat].weightColumn)))
553 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
554#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
555 for (i = 0; i < rows; i++) {
556 stat[iStat].sumWeight[i] += weight[i];
557 stat[iStat].value1[i] += inputData[i] * weight[i];
558 }
559 free(weight);
560 break;
561 case SET_SDS:
562 case SET_SIGMAS:
563 if (code == 1)
564#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
565 for (i = 0; i < rows; i++) {
566 stat[iStat].value1[i] = inputData[i];
567 stat[iStat].value2[i] = inputData[i] * inputData[i];
568 }
569 else
570#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
571 for (i = 0; i < rows; i++) {
572 stat[iStat].value1[i] += inputData[i];
573 stat[iStat].value2[i] += inputData[i] * inputData[i];
574 }
575 break;
576 case SET_WSDS:
577 case SET_WSIGMAS:
578 if (!(weight = SDDS_GetColumnInDoubles(&inTable, stat[iStat].weightColumn)))
579 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
580#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
581 for (i = 0; i < rows; i++) {
582 stat[iStat].sumWeight[i] += weight[i];
583 stat[iStat].value1[i] += inputData[i] * weight[i];
584 stat[iStat].value2[i] += inputData[i] * inputData[i] * weight[i];
585 }
586 free(weight);
587 break;
588 case SET_RMSS:
589 if (code == 1)
590#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
591 for (i = 0; i < rows; i++)
592 stat[iStat].value1[i] = inputData[i] * inputData[i];
593 else
594#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
595 for (i = 0; i < rows; i++)
596 stat[iStat].value1[i] += inputData[i] * inputData[i];
597 break;
598 case SET_WRMSS:
599 if (!(weight = SDDS_GetColumnInDoubles(&inTable, stat[iStat].weightColumn)))
600 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
601#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
602 for (i = 0; i < rows; i++) {
603 stat[iStat].sumWeight[i] += weight[i];
604 stat[iStat].value1[i] += inputData[i] * inputData[i] * weight[i];
605 }
606 free(weight);
607 break;
608 case SET_SUMS:
609 if (code == 1)
610#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
611 for (i = 0; i < rows; i++)
612 stat[iStat].value1[i] = ipow(inputData[i], stat[iStat].sumPower);
613 else
614#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
615 for (i = 0; i < rows; i++)
616 stat[iStat].value1[i] += ipow(inputData[i], stat[iStat].sumPower);
617 break;
618 case SET_PMINIMA:
619 if (!SDDS_GetParameterAsDouble(&inTable, stat[iStat].functionOf, &indepData)) {
620 fprintf(stderr, "error: unable to get value of parameter %s\n", stat[iStat].functionOf);
621 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
622 }
623 if (code == 1)
624#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
625 for (i = 0; i < rows; i++) {
626 stat[iStat].value2[i] = inputData[i];
627 stat[iStat].value1[i] = indepData;
628 }
629 else
630#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
631 for (i = 0; i < rows; i++) {
632 if (stat[iStat].value2[i] > inputData[i]) {
633 stat[iStat].value2[i] = inputData[i];
634 stat[iStat].value1[i] = indepData;
635 }
636 }
637 break;
638 case SET_PMAXIMA:
639 if (!SDDS_GetParameterAsDouble(&inTable, stat[iStat].functionOf, &indepData)) {
640 fprintf(stderr, "error: unable to get value of parameter %s\n", stat[iStat].functionOf);
641 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
642 }
643 if (code == 1)
644#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
645 for (i = 0; i < rows; i++) {
646 stat[iStat].value2[i] = inputData[i];
647 stat[iStat].value1[i] = indepData;
648 }
649 else
650#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
651 for (i = 0; i < rows; i++) {
652 if (stat[iStat].value2[i] < inputData[i]) {
653 stat[iStat].value2[i] = inputData[i];
654 stat[iStat].value1[i] = indepData;
655 }
656 }
657 break;
658 case SET_SLOPE:
659 case SET_INTERCEPT:
660 if (!SDDS_GetParameterAsDouble(&inTable, stat[iStat].functionOf, &indepData)) {
661 fprintf(stderr, "error: unable to get value of parameter %s\n", stat[iStat].functionOf);
662 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
663 }
664 /* linear fit:
665 y = a + bx
666 a = (S x^2 Sy - S x S xy)/D
667 b = (N S xy - Sx Sy)/D
668 D = N S x^2 - (S x)^2
669 */
670 if (code == 1)
671#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
672 for (i = 0; i < rows; i++) {
673 stat[iStat].value1[i] = indepData; /* Sum x */
674 stat[iStat].value2[i] = indepData * indepData; /* Sum x^2 */
675 stat[iStat].value3[i] = inputData[i]; /* Sum y */
676 stat[iStat].value4[i] = indepData * inputData[i]; /* Sum xy */
677 }
678 else
679#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
680 for (i = 0; i < rows; i++) {
681 stat[iStat].value1[i] += indepData;
682 stat[iStat].value2[i] += indepData * indepData;
683 stat[iStat].value3[i] += inputData[i];
684 stat[iStat].value4[i] += indepData * inputData[i];
685 }
686 break;
687 case SET_MEDIAN:
688 case SET_DRANGE:
689 case SET_PERCENTILE:
690 case SET_EXMM_MEAN:
691 if (code == 1)
692#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
693 for (i = 0; i < rows; i++) {
694 stat[iStat].array[i] = tmalloc(sizeof(*stat[iStat].array[i]));
695 stat[iStat].array[i][pages - 1] = inputData[i];
696 }
697 else {
698#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
699 for (i = 0; i < rows; i++) {
700 stat[iStat].array[i] = SDDS_Realloc(stat[iStat].array[i], sizeof(*stat[iStat].array[i]) * pages);
701 stat[iStat].array[i][pages - 1] = inputData[i];
702 }
703 }
704 break;
705 default:
706 SDDS_Bomb("invalid statistic code (accumulation loop)");
707 break;
708 }
709 free(inputData);
710 }
711 }
712 if (pages == 0)
713 SDDS_Bomb("no pages in input");
714 for (iStat = 0; iStat < stats; iStat++) {
715 switch (stat[iStat].optionCode) {
716 case SET_COPY:
717 case SET_MINIMA:
718 case SET_MAXIMA:
719 case SET_PMINIMA:
720 case SET_PMAXIMA:
721 case SET_CMINIMA:
722 case SET_CMAXIMA:
723 case SET_LARGEST:
724 case SET_SIGNEDLARGEST:
725 case SET_SUMS:
726 break;
727 case SET_MEANS:
728#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
729 for (i = 0; i < rows; i++)
730 stat[iStat].value1[i] /= pages;
731 break;
732 case SET_WMEANS:
733 for (i = 0; i < rows; i++)
734 if (stat[iStat].sumWeight[i])
735 stat[iStat].value1[i] /= stat[iStat].sumWeight[i];
736 else {
737 if (!nowarnings)
738 fprintf(stderr, "warning: the total weight for the %" PRId64 "th row of %s is zero.\n", i + 1, stat[iStat].sourceColumn);
739 stat[iStat].value1[i] = DBL_MAX;
740 }
741 break;
742 case SET_SDS:
743 if (pages < 2) {
744#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
745 for (i = 0; i < rows; i++)
746 stat[iStat].value1[i] = DBL_MAX;
747 } else {
748#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
749 for (i = 0; i < rows; i++) {
750 double tmp1;
751 if ((tmp1 = stat[iStat].value2[i] / pages - sqr(stat[iStat].value1[i] / pages)) <= 0)
752 stat[iStat].value1[i] = 0;
753 else
754 stat[iStat].value1[i] = sqrt(tmp1 * pages / (pages - 1.0));
755 }
756 }
757 break;
758 case SET_WSDS:
759 if (pages < 2) {
760#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
761 for (i = 0; i < rows; i++)
762 stat[iStat].value1[i] = DBL_MAX;
763 } else {
764 for (i = 0; i < rows; i++) {
765 double tmp1;
766 if (stat[iStat].sumWeight[i]) {
767 if ((tmp1 = stat[iStat].value2[i] / stat[iStat].sumWeight[i] - sqr(stat[iStat].value1[i] / stat[iStat].sumWeight[i])) <= 0)
768 stat[iStat].value1[i] = 0;
769 else
770 stat[iStat].value1[i] = sqrt(tmp1 * pages / (pages - 1.0));
771 } else {
772 if (!nowarnings)
773 fprintf(stderr, "Warning, the total weight for the %" PRId64 "th row of %s is zero.\n", i + 1, stat[iStat].sourceColumn);
774 stat[iStat].value1[i] = DBL_MAX;
775 }
776 }
777 }
778 break;
779 case SET_SIGMAS:
780 if (pages < 2) {
781#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
782 for (i = 0; i < rows; i++)
783 stat[iStat].value1[i] = DBL_MAX;
784 } else {
785#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
786 for (i = 0; i < rows; i++) {
787 double tmp1;
788 if ((tmp1 = stat[iStat].value2[i] / pages - sqr(stat[iStat].value1[i] / pages)) <= 0)
789 stat[iStat].value1[i] = 0;
790 else
791 stat[iStat].value1[i] = sqrt(tmp1 / (pages - 1.0));
792 }
793 }
794 break;
795 case SET_WSIGMAS:
796 if (pages < 2) {
797#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
798 for (i = 0; i < rows; i++)
799 stat[iStat].value1[i] = DBL_MAX;
800 } else {
801 for (i = 0; i < rows; i++) {
802 double tmp1;
803 if (stat[iStat].sumWeight[i]) {
804 if ((tmp1 = stat[iStat].value2[i] / stat[iStat].sumWeight[i] - sqr(stat[iStat].value1[i] / stat[iStat].sumWeight[i])) <= 0)
805 stat[iStat].value1[i] = 0;
806 else
807 stat[iStat].value1[i] = sqrt(tmp1 / (pages - 1.0));
808 } else {
809 if (!nowarnings)
810 fprintf(stderr, "Warning, the total weight for the %" PRId64 "th row of %s is zero.\n", i + 1, stat[iStat].sourceColumn);
811 stat[iStat].value1[i] = DBL_MAX;
812 }
813 }
814 }
815 break;
816 case SET_RMSS:
817#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
818 for (i = 0; i < rows; i++)
819 stat[iStat].value1[i] = sqrt(stat[iStat].value1[i] / pages);
820 break;
821 case SET_WRMSS:
822 for (i = 0; i < rows; i++) {
823 if (stat[iStat].sumWeight[i])
824 stat[iStat].value1[i] = sqrt(stat[iStat].value1[i] / stat[iStat].sumWeight[i]);
825 else {
826 if (!nowarnings)
827 fprintf(stderr, "Warning, the total weight for the %" PRId64 "th row of %s is zero.\n", i + 1, stat[iStat].sourceColumn);
828 stat[iStat].value1[i] = DBL_MAX;
829 }
830 }
831 break;
832 case SET_SLOPE:
833#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
834 for (i = 0; i < rows; i++) {
835 double D;
836 D = pages * stat[iStat].value2[i] - stat[iStat].value1[i] * stat[iStat].value1[i];
837 stat[iStat].value1[i] = (pages * stat[iStat].value4[i] - stat[iStat].value1[i] * stat[iStat].value3[i]) / D;
838 }
839 break;
840 case SET_INTERCEPT:
841#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
842 for (i = 0; i < rows; i++) {
843 double D;
844 D = pages * stat[iStat].value2[i] - stat[iStat].value1[i] * stat[iStat].value1[i];
845 stat[iStat].value1[i] = (stat[iStat].value2[i] * stat[iStat].value3[i] - stat[iStat].value1[i] * stat[iStat].value4[i]) / D;
846 }
847 break;
848 case SET_MEDIAN:
849#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
850 for (i = 0; i < rows; i++) {
851 compute_median(&stat[iStat].value1[i], stat[iStat].array[i], pages);
852 }
853 break;
854 case SET_DRANGE:
855#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
856 for (i = 0; i < rows; i++) {
857 double decileResult[2];
858 if (!compute_percentiles(decileResult, decilePoint, 2, stat[iStat].array[i], pages))
859 stat[iStat].value1[i] = 0;
860 else
861 stat[iStat].value1[i] = decileResult[1] - decileResult[0];
862 }
863 break;
864 case SET_PERCENTILE:
865 percentilePoint = stat[iStat].percentile;
866#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
867 for (i = 0; i < rows; i++) {
868 double percentileResult;
869 if (!compute_percentiles(&percentileResult, &percentilePoint, 1, stat[iStat].array[i], pages))
870 stat[iStat].value1[i] = 0;
871 else
872 stat[iStat].value1[i] = percentileResult;
873 }
874 break;
875 case SET_EXMM_MEAN:
876#pragma omp parallel for if (threads > 1 && rows > 1) num_threads(threads)
877 for (i = 0; i < rows; i++)
878 if (!compute_mean_exclude_min_max(&(stat[iStat].value1[i]), stat[iStat].array[i], pages))
879 stat[iStat].value1[i] = 0;
880 break;
881 default:
882 SDDS_Bomb("invalid statistic code (final loop)");
883 break;
884 }
885 if (stat[iStat].optionCode == SET_COPY) {
886 if (!SDDS_SetColumn(&outTable, SDDS_SET_BY_NAME, stat[iStat].copy, rows, stat[iStat].resultColumn)) {
887 fprintf(stderr, "error setting column values for column %s\n", stat[iStat].resultColumn);
888 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
889 }
890 } else if (!SDDS_SetColumnFromDoubles(&outTable, SDDS_SET_BY_NAME, stat[iStat].value1, rows, stat[iStat].resultColumn)) {
891 fprintf(stderr, "error setting column values for column %s\n", stat[iStat].resultColumn);
892 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
893 }
894 if (stat[iStat].value1)
895 free(stat[iStat].value1);
896 if (stat[iStat].value2)
897 free(stat[iStat].value2);
898 if (stat[iStat].value3)
899 free(stat[iStat].value3);
900 if (stat[iStat].value4)
901 free(stat[iStat].value4);
902 if (stat[iStat].copy)
903 free(stat[iStat].copy);
904 if (stat[iStat].array) {
905 for (i = 0; i < rows; i++) {
906 free(stat[iStat].array[i]);
907 }
908 free(stat[iStat].array);
909 }
910 if (stat[iStat].sumWeight)
911 free(stat[iStat].sumWeight);
912 free(stat[iStat].sourceColumn);
913 free(stat[iStat].resultColumn);
914 stat[iStat].value1 = stat[iStat].value2 = stat[iStat].value3 = stat[iStat].value4 = NULL;
915 stat[iStat].copy = NULL;
916 stat[iStat].array = NULL;
917 }
918 free(stat);
919 if (!SDDS_WritePage(&outTable) || !SDDS_Terminate(&inTable) || !SDDS_Terminate(&outTable))
920 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
921 return EXIT_SUCCESS;
922}
923
924long addStatRequests(STAT_REQUEST **statRequest, long requests, char **item, long items, long code, double percentile, long power, char *functionOf, long weighted, char *percentileString) {
925 long i;
926 /*weighted factor should be either 0 or 1 */
927 if (weighted != 0 && weighted != 1)
928 SDDS_Bomb("addStatRequests: weighted parameter should be either 0 or 1");
929 if (code == SET_PERCENTILE && (!percentileString || !strlen(percentileString))) {
930 fprintf(stderr, "Percentile specification is incorrect: percentile=%e, percentileString=%s\n", percentile, percentileString ? percentileString : "NULL");
931 exit(EXIT_FAILURE);
932 }
933 *statRequest = SDDS_Realloc(*statRequest, sizeof(**statRequest) * (requests + items - weighted));
934 for (i = 0; i < items - weighted; i++) {
935 if (weighted)
936 (*statRequest)[requests + i].weightColumnName = item[0];
937 else
938 (*statRequest)[requests + i].weightColumnName = NULL;
939 (*statRequest)[i + requests].columnName = item[i + weighted];
940 (*statRequest)[i + requests].optionCode = code;
941 (*statRequest)[i + requests].sumPower = power;
942 (*statRequest)[i + requests].percentile = percentile;
943 (*statRequest)[i + requests].percentileString = percentileString;
944 (*statRequest)[i + requests].functionOf = functionOf;
945 }
946
947 return items + requests - weighted;
948}
949
950STAT_DEFINITION *compileStatDefinitions(SDDS_DATASET *inTable, long *stats, STAT_REQUEST *request, long requests) {
951 STAT_DEFINITION *stat;
952 long iReq, iStat, iName;
953 int32_t columnNames;
954 char s[SDDS_MAXLINE];
955 char **columnName;
956
957 *stats = iStat = 0;
958 stat = tmalloc(sizeof(*stat) * requests);
959 for (iReq = 0; iReq < requests; iReq++) {
960 if (iStat >= *stats)
961 stat = SDDS_Realloc(stat, sizeof(*stat) * (*stats += 10));
962 if (!has_wildcards(request[iReq].columnName)) {
963 if (SDDS_GetColumnIndex(inTable, request[iReq].columnName) < 0) {
964 sprintf(s, "error: column %s not found input file", request[iReq].columnName);
965 SDDS_SetError(s);
966 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
967 }
968 stat[iStat].weightColumn = request[iReq].weightColumnName;
969 stat[iStat].sourceColumn = request[iReq].columnName;
970 stat[iStat].optionCode = request[iReq].optionCode;
971 stat[iStat].percentile = request[iReq].percentile;
972 stat[iStat].percentileString = request[iReq].percentileString;
973 stat[iStat].sumPower = request[iReq].sumPower;
974 stat[iStat].value1 = stat[iStat].value2 = stat[iStat].value3 = stat[iStat].value4 = NULL;
975 stat[iStat].array = NULL;
976 stat[iStat].copy = NULL;
977 stat[iStat].sumWeight = NULL;
978 if ((stat[iStat].functionOf = request[iReq].functionOf)) {
979 if (stat[iStat].optionCode != SET_CMAXIMA && stat[iStat].optionCode != SET_CMINIMA) {
980 if (SDDS_GetParameterIndex(inTable, request[iReq].functionOf) < 0) {
981 sprintf(s, "error: parameter %s not found input file (1)", request[iReq].functionOf);
982 SDDS_SetError(s);
983 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
984 }
985 } else {
986 if (SDDS_GetColumnIndex(inTable, request[iReq].functionOf) < 0) {
987 sprintf(s, "error: column %s not found input file (1)", request[iReq].functionOf);
988 SDDS_SetError(s);
989 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
990 }
991 }
992 }
993 iStat++;
994 } else {
995 SDDS_SetColumnFlags(inTable, 0);
996 if (!SDDS_SetColumnsOfInterest(inTable, SDDS_MATCH_STRING, request[iReq].columnName, SDDS_OR))
997 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
998 if (!(columnName = SDDS_GetColumnNames(inTable, &columnNames))) {
999 sprintf(s, "no columns selected for wildcard sequence %s", request[iReq].columnName);
1000 SDDS_SetError(s);
1001 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1002 }
1003 if (iStat + columnNames > *stats)
1004 stat = SDDS_Realloc(stat, sizeof(*stat) * (*stats = iStat + columnNames + 10));
1005 for (iName = 0; iName < columnNames; iName++) {
1006 stat[iStat + iName].weightColumn = request[iReq].weightColumnName;
1007 stat[iStat + iName].sourceColumn = columnName[iName];
1008 stat[iStat + iName].optionCode = request[iReq].optionCode;
1009 stat[iStat + iName].sumPower = request[iReq].sumPower;
1010 stat[iStat + iName].percentile = request[iReq].percentile;
1011 stat[iStat + iName].percentileString = request[iReq].percentileString;
1012 stat[iStat + iName].value1 = stat[iStat + iName].value2 = stat[iStat + iName].value3 = stat[iStat + iName].value4 = NULL;
1013 stat[iStat + iName].array = NULL;
1014 stat[iStat + iName].copy = NULL;
1015 stat[iStat + iName].sumWeight = NULL;
1016 if ((stat[iStat + iName].functionOf = request[iReq].functionOf) && iName == 0) {
1017 if (stat[iStat + iName].optionCode != SET_CMAXIMA && stat[iStat + iName].optionCode != SET_CMINIMA) {
1018 if (SDDS_GetParameterIndex(inTable, request[iReq].functionOf) < 0) {
1019 sprintf(s, "error: parameter %s not found input file (2)", request[iReq].functionOf);
1020 SDDS_SetError(s);
1021 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1022 }
1023 } else {
1024 if (SDDS_GetColumnIndex(inTable, request[iReq].functionOf) < 0) {
1025 sprintf(s, "error: column %s not found input file (2)", request[iReq].functionOf);
1026 SDDS_SetError(s);
1027 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1028 }
1029 }
1030 }
1031 }
1032 iStat += columnNames;
1033 free(columnName);
1034 }
1035 }
1036
1037 *stats = iStat;
1038 for (iStat = 0; iStat < *stats; iStat++) {
1039 switch (stat[iStat].optionCode) {
1040 case SET_COPY:
1041 strcpy(s, stat[iStat].sourceColumn);
1042 break;
1043 case SET_SUMS:
1044 if (stat[iStat].sumPower == 1)
1045 sprintf(s, "%s%s", stat[iStat].sourceColumn, optionSuffix[stat[iStat].optionCode]);
1046 else
1047 sprintf(s, "%s%ld%s", stat[iStat].sourceColumn, stat[iStat].sumPower, optionSuffix[stat[iStat].optionCode]);
1048 break;
1049 case SET_PERCENTILE:
1050 sprintf(s, "%s%s%s", stat[iStat].sourceColumn, stat[iStat].percentileString, optionSuffix[stat[iStat].optionCode]);
1051 break;
1052 case SET_PMAXIMA:
1053 case SET_PMINIMA:
1054 sprintf(s, "%s%s%s", stat[iStat].functionOf, optionSuffix[stat[iStat].optionCode], stat[iStat].sourceColumn);
1055 break;
1056 case SET_CMAXIMA:
1057 case SET_CMINIMA:
1058 sprintf(s, "%s%s%s", stat[iStat].functionOf, optionSuffix[stat[iStat].optionCode], stat[iStat].sourceColumn);
1059 break;
1060 default:
1061 sprintf(s, "%s%s", stat[iStat].sourceColumn, optionSuffix[stat[iStat].optionCode]);
1062 break;
1063 }
1064 if (!SDDS_CopyString(&stat[iStat].resultColumn, s))
1065 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
1066 }
1067 return stat;
1068}
1069
1070long setupOutputFile(SDDS_DATASET *outTable, char *output, SDDS_DATASET *inTable, STAT_DEFINITION *stat, long stats, int64_t rows, short columnMajorOrder) {
1071 long column;
1072 char s[SDDS_MAXLINE], *symbol, *symbol1, *units1;
1073
1074 if (!SDDS_InitializeOutput(outTable, SDDS_BINARY, 0, NULL, "sddsenvelope output", output))
1075 return 0;
1076 if (columnMajorOrder != -1)
1077 outTable->layout.data_mode.column_major = columnMajorOrder;
1078 else
1079 outTable->layout.data_mode.column_major = inTable->layout.data_mode.column_major;
1080 for (column = 0; column < stats; column++) {
1081 stat[column].value1 = calloc(sizeof(*stat[column].value1), rows);
1082 stat[column].value2 = stat[column].value3 = stat[column].value4 = NULL;
1083 if (stat[column].optionCode == SET_SDS || stat[column].optionCode == SET_SIGMAS ||
1084 stat[column].optionCode == SET_WSDS || stat[column].optionCode == SET_WSIGMAS ||
1085 stat[column].optionCode == SET_PMINIMA || stat[column].optionCode == SET_PMAXIMA ||
1086 stat[column].optionCode == SET_CMINIMA || stat[column].optionCode == SET_CMAXIMA)
1087 stat[column].value2 = calloc(sizeof(*stat[column].value2), rows);
1088 if (stat[column].optionCode == SET_INTERCEPT || stat[column].optionCode == SET_SLOPE) {
1089 stat[column].value2 = malloc(sizeof(*stat[column].value2) * rows);
1090 stat[column].value3 = malloc(sizeof(*stat[column].value3) * rows);
1091 stat[column].value4 = malloc(sizeof(*stat[column].value4) * rows);
1092 }
1093 if (stat[column].optionCode == SET_WSDS || stat[column].optionCode == SET_WSIGMAS ||
1094 stat[column].optionCode == SET_WRMSS || stat[column].optionCode == SET_WMEANS)
1095 stat[column].sumWeight = calloc(sizeof(*stat[column].sumWeight), rows);
1096 if (stat[column].optionCode == SET_MEDIAN || stat[column].optionCode == SET_DRANGE ||
1097 stat[column].optionCode == SET_PERCENTILE || stat[column].optionCode == SET_EXMM_MEAN) {
1098 stat[column].array = tmalloc(sizeof(*stat[column].array) * rows);
1099 }
1100 if (!SDDS_TransferColumnDefinition(outTable, inTable, stat[column].sourceColumn, stat[column].resultColumn)) {
1101 sprintf(s, "Problem transferring definition of column %s to %s\n", stat[column].sourceColumn, stat[column].resultColumn);
1102 SDDS_SetError(s);
1103 return 0;
1104 }
1105 if (SDDS_ChangeColumnInformation(outTable, "description", NULL, SDDS_SET_BY_NAME, stat[column].resultColumn) != SDDS_STRING ||
1106 SDDS_GetColumnInformation(outTable, "symbol", &symbol, SDDS_BY_NAME, stat[column].resultColumn) != SDDS_STRING) {
1107 fprintf(stderr, "Error: problem setting description for column %s\n", stat[column].resultColumn);
1108 return 0;
1109 }
1110 if (stat[column].optionCode > 0) {
1111 if (SDDS_ChangeColumnInformation(outTable, "type", "double", SDDS_PASS_BY_STRING | SDDS_SET_BY_NAME, stat[column].resultColumn) != SDDS_LONG) {
1112 fprintf(stderr, "Error: problem setting type for column %s\n", stat[column].resultColumn);
1113 return 0;
1114 }
1115 }
1116 if (!symbol)
1117 SDDS_CopyString(&symbol, stat[column].sourceColumn);
1118 switch (stat[column].optionCode) {
1119 case SET_COPY:
1120 strcpy(s, symbol);
1121 break;
1122 case SET_SUMS:
1123 if (stat[column].sumPower == 1)
1124 sprintf(s, "%s[%s]", optionSuffix[stat[column].optionCode], symbol);
1125 else
1126 sprintf(s, "%s[%s$a%ld$n]", optionSuffix[stat[column].optionCode], symbol, stat[column].sumPower);
1127 break;
1128 case SET_PERCENTILE:
1129 sprintf(s, "%s[%s,%g]", optionSuffix[stat[column].optionCode], symbol, stat[column].percentile);
1130 break;
1131 case SET_PMINIMA:
1132 case SET_PMAXIMA:
1133 if (SDDS_GetParameterInformation(inTable, "symbol", &symbol1, SDDS_BY_NAME, stat[column].functionOf) != SDDS_STRING ||
1134 !symbol1 ||
1135 !strlen(symbol1))
1136 symbol1 = stat[column].functionOf;
1137 sprintf(s, "%s[%s:%s]", optionSuffix[stat[column].optionCode], symbol, symbol1);
1138 if (SDDS_GetParameterInformation(inTable, "units", &units1, SDDS_BY_NAME, stat[column].functionOf) != SDDS_STRING)
1139 return 0;
1140 if (units1) {
1141 if (SDDS_ChangeColumnInformation(outTable, "units", units1, SDDS_BY_NAME, stat[column].resultColumn) != SDDS_STRING) {
1142 fprintf(stderr, "Error: problem setting units for column %s (1)\n", stat[column].resultColumn);
1143 return 0;
1144 }
1145 } else if (SDDS_ChangeColumnInformation(outTable, "units", "", SDDS_BY_NAME, stat[column].resultColumn) != SDDS_STRING) {
1146 fprintf(stderr, "Error: problem setting units for column %s (2)\n", stat[column].resultColumn);
1147 return 0;
1148 }
1149 break;
1150 case SET_CMINIMA:
1151 case SET_CMAXIMA:
1152 if (SDDS_GetColumnInformation(inTable, "symbol", &symbol1, SDDS_BY_NAME, stat[column].functionOf) != SDDS_STRING ||
1153 !symbol1 ||
1154 !strlen(symbol1))
1155 symbol1 = stat[column].functionOf;
1156 sprintf(s, "%s[%s:%s]", optionSuffix[stat[column].optionCode], symbol, symbol1);
1157 if (SDDS_GetColumnInformation(inTable, "units", &units1, SDDS_BY_NAME, stat[column].resultColumn) != SDDS_STRING)
1158 return 0;
1159 if (units1) {
1160 if (SDDS_ChangeColumnInformation(outTable, "units", units1, SDDS_BY_NAME, stat[column].resultColumn) != SDDS_STRING) {
1161 fprintf(stderr, "Error: problem setting units for column %s\n", stat[column].resultColumn);
1162 return 0;
1163 }
1164 } else if (SDDS_ChangeColumnInformation(outTable, "units", "", SDDS_BY_NAME, stat[column].resultColumn) != SDDS_STRING) {
1165 fprintf(stderr, "Error: problem setting units for column %s\n", stat[column].resultColumn);
1166 return 0;
1167 }
1168 break;
1169 case SET_INTERCEPT:
1170 case SET_SLOPE:
1171 if (SDDS_GetParameterInformation(inTable, "symbol", &symbol1, SDDS_BY_NAME, stat[column].functionOf) != SDDS_STRING ||
1172 !symbol1 ||
1173 !strlen(symbol1))
1174 symbol1 = stat[column].functionOf;
1175 sprintf(s, "%s[%s:%s]", optionSuffix[stat[column].optionCode], symbol, symbol1);
1176 break;
1177 default:
1178 sprintf(s, "%s[%s]", optionSuffix[stat[column].optionCode], symbol);
1179 break;
1180 }
1181 free(symbol);
1182 if (SDDS_ChangeColumnInformation(outTable, "symbol", s, SDDS_BY_NAME, stat[column].resultColumn) != SDDS_STRING) {
1183 fprintf(stderr, "Error: problem setting symbol for column %s\n", stat[column].resultColumn);
1184 return 0;
1185 }
1186 }
1187 if (!SDDS_WriteLayout(outTable) || !SDDS_StartPage(outTable, rows))
1188 return 0;
1189 return 1;
1190}
1191
1192int compute_mean_exclude_min_max(double *value, double *data, long n) {
1193 double sum = 0;
1194 long i, count = 0;
1195 double min, max;
1196 max = -(min = DBL_MAX);
1197 if (n <= 0 || !data || !value)
1198 return 0;
1199 for (i = 0; i < n; i++) {
1200 if (min > data[i])
1201 min = data[i];
1202 if (max < data[i])
1203 max = data[i];
1204 }
1205 for (i = 0; i < n; i++) {
1206 if (data[i] == min || data[i] == max)
1207 continue;
1208 count++;
1209 sum += data[i];
1210 }
1211 if (count == 0)
1212 *value = min;
1213 else
1214 *value = sum / count;
1215 return 1;
1216}
SDDS.h
SDDS (Self Describing Data Set) Data Types Definitions and Function Prototypes.
SDDS_StartPage
int32_t SDDS_StartPage(SDDS_DATASET *SDDS_dataset, int64_t expected_n_rows)
Definition SDDS_dataprep.c:70
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:1648
SDDS_SetColumn
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.
Definition SDDS_dataprep.c:1542
SDDS_GetColumn
void * SDDS_GetColumn(SDDS_DATASET *SDDS_dataset, char *column_name)
Retrieves a copy of the data for a specified column, including only rows marked as "of interest".
Definition SDDS_extract.c:611
SDDS_GetParameterAsDouble
double * SDDS_GetParameterAsDouble(SDDS_DATASET *SDDS_dataset, char *parameter_name, double *memory)
Retrieves the value of a specified parameter as a double from the current data table of an SDDS datas...
Definition SDDS_extract.c:2783
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_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_SetColumnFlags
int32_t SDDS_SetColumnFlags(SDDS_DATASET *SDDS_dataset, int32_t column_flag_value)
Sets the acceptance flags for all columns in the current data table of a data set.
Definition SDDS_extract.c:218
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_GetParameterInformation
int32_t SDDS_GetParameterInformation(SDDS_DATASET *SDDS_dataset, char *field_name, void *memory, int32_t mode,...)
Retrieves information about a specified parameter in the SDDS dataset.
Definition SDDS_info.c:117
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_InitializeInput
int32_t SDDS_InitializeInput(SDDS_DATASET *SDDS_dataset, char *filename)
Definition SDDS_input.c:50
SDDS_Terminate
int32_t SDDS_Terminate(SDDS_DATASET *SDDS_dataset)
Definition SDDS_input.c:1471
SDDS_ReadPage
int32_t SDDS_ReadPage(SDDS_DATASET *SDDS_dataset)
Definition SDDS_input.c:1007
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:602
SDDS_WritePage
int32_t SDDS_WritePage(SDDS_DATASET *SDDS_dataset)
Writes the current data table to the output file.
Definition SDDS_output.c:1229
SDDS_WriteLayout
int32_t SDDS_WriteLayout(SDDS_DATASET *SDDS_dataset)
Writes the SDDS layout header to the output file.
Definition SDDS_output.c:900
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_SetError
void SDDS_SetError(char *error_text)
Records an error message in the SDDS error stack.
Definition SDDS_utils.c:421
SDDS_GetParameterIndex
int32_t SDDS_GetParameterIndex(SDDS_DATASET *SDDS_dataset, char *name)
Retrieves the index of a named parameter in the SDDS dataset.
Definition SDDS_utils.c:1404
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:1375
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:2526
SDDS_PrintErrors
void SDDS_PrintErrors(FILE *fp, int32_t mode)
Prints recorded error messages to a specified file stream.
Definition SDDS_utils.c:474
SDDS_RegisterProgramName
void SDDS_RegisterProgramName(const char *name)
Registers the executable program name for use in error messages.
Definition SDDS_utils.c:318
SDDS_NumberOfErrors
int32_t SDDS_NumberOfErrors()
Retrieves the number of errors recorded by SDDS library routines.
Definition SDDS_utils.c:340
SDDS_Bomb
void SDDS_Bomb(char *message)
Terminates the program after printing an error message and recorded errors.
Definition SDDS_utils.c:380
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:922
SDDS_Realloc
void * SDDS_Realloc(void *old_ptr, size_t new_size)
Reallocates memory to a new size.
Definition SDDS_utils.c:743
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
tmalloc
void * tmalloc(uint64_t size_of_block)
Allocates a memory block of the specified size with zero initialization.
Definition array.c:65
bomb
void bomb(char *error, char *usage)
Reports error messages to the terminal and aborts the program.
Definition bomb.c:26
ipow
double ipow(const double x, const int64_t p)
Compute x raised to the power p (x^p).
Definition ipow.c:33
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
compute_percentiles
long compute_percentiles(double *position, double *percent, long positions, double *x, long n)
Computes multiple percentiles of an array of doubles.
Definition median.c:84
compute_median
long compute_median(double *value, double *x, long n)
Computes the median of an array of doubles.
Definition median.c:29
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:357
processFilenames
void processFilenames(char *programName, char **input, char **output, unsigned long pipeFlags, long noWarnings, long *tmpOutputUsed)
Definition scanargs.c:391
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
STAT_DEFINITION
Definition sddsenvelope.c:205
STAT_REQUEST
Definition sddsenvelope.c:192
has_wildcards
int has_wildcards(char *template)
Check if a template string contains any wildcard characters.
Definition wild_match.c:498
Generated by doxygen 1.12.0