SDDS ToolKit Programs and Libraries for C and Python
Loading...
Searching...
No Matches
sddshist.c
Go to the documentation of this file.
1/**
2 * @file sddshist.c
3 * @brief Command-line tool for generating histograms from SDDS-formatted data.
4 *
5 * @details
6 * This program is a command-line tool designed to generate histograms from data contained in SDDS (Self Describing Data Set) files.
7 * It supports various configurations like normalization, filtering, multi-threading for performance, and the ability to define
8 * histograms based on specific regions in the data.
9 *
10 * @section Usage
11 * ```
12 * sddshist [<inputfile>] [<outputfile>]
13 * [-pipe=[input][,output]]
14 * -dataColumn=<column-name>
15 * [{
16 * -bins=<number> |
17 * -sizeOfBins=<value> |
18 * -regions=filename=<filename>,position=<columnName>,name=<columnName>
19 * }]
20 * [-lowerLimit=<value>]
21 * [-upperLimit=<value>]
22 * [-expand=<factor>]
23 * [-filter=<column-name>,<lower-limit>,<upper-limit>]
24 * [-weightColumn=<column-name>]
25 * [-sides[=<points>]]
26 * [-normalize[={sum|area|peak}]]
27 * [-cdf[=only]]
28 * [-threads=<number>]
29 * [-statistics]
30 * [-verbose]
31 * [-majorOrder=row|column]
32 * ```
33 *
34 * @section Options
35 * | Required | Description |
36 * |----------|-----------------------------------------------------------------------------|
37 * | `-dataColumn=<column-name>` | Specifies the column to be histogrammed. |
38 *
39 * | Optional | Description |
40 * |----------|-----------------------------------------------------------------------------|
41 * | `-pipe` | Use pipe for input and/or output. |
42 * | `-bins` | Set the number of bins for the histogram. |
43 * | `-sizeOfBins` | Set the size of each bin. |
44 * | `-regions` | Define region-based histogramming with file and column information. |
45 * | `-lowerLimit` | Set the lower limit for the histogram. |
46 * | `-upperLimit` | Set the upper limit for the histogram. |
47 * | `-expand` | Expand the histogram range by a given factor. |
48 * | `-filter` | Filter data points based on values in a specified column. |
49 * | `-weightColumn` | Specify a column to weight the histogram. |
50 * | `-sides` | Extend the histogram to zero level. |
51 * | `-normalize` | Normalize the histogram (`sum`, `area`, or `peak`). |
52 * | `-cdf` | Include the Cumulative Distribution Function (CDF) in the output. |
53 * | `-threads` | Set the number of threads for processing. |
54 * | `-statistics` | Include statistical details in the output. |
55 * | `-verbose` | Print additional processing details. |
56 * | `-majorOrder` | Set data major order (row or column). |
57 *
58 * @subsection Incompatibilities
59 * - `-bins`, `-sizeOfBins`, and `-regions` are mutually exclusive.
60 * - `-normalize` and `-cdf=only` cannot be used together.
61 *
62 * @copyright
63 * - (c) 2002 The University of Chicago, as Operator of Argonne National Laboratory.
64 * - (c) 2002 The Regents of the University of California, as Operator of Los Alamos National Laboratory.
65 *
66 * @license
67 * This file is distributed under the terms of the Software License Agreement
68 * found in the file LICENSE included with this distribution.
69 *
70 * @authors
71 * M. Borland, C. Saunders, R. Soliday, H. Shang
72 */
73
74#include "mdb.h"
75#include "scan.h"
76#include "SDDS.h"
77
78/* Enumeration for option types */
79enum option_type {
80 SET_BINS,
81 SET_LOWERLIMIT,
82 SET_UPPERLIMIT,
83 SET_DATACOLUMN,
84 SET_FILTER,
85 SET_BINSIZE,
86 SET_WEIGHTCOLUMN,
87 SET_NORMALIZE,
88 SET_STATISTICS,
89 SET_SIDES,
90 SET_VERBOSE,
91 SET_PIPE,
92 SET_CDF,
93 SET_EXPAND,
94 SET_MAJOR_ORDER,
95 SET_REGION_FILE,
96 SET_THREADS,
97 N_OPTIONS
98};
99
100char *option[N_OPTIONS] = {
101 "bins", "lowerlimit", "upperlimit", "datacolumn",
102 "filter", "sizeofbins", "weightcolumn",
103 "normalize", "statistics", "sides", "verbose", "pipe", "cdf", "expand",
104 "majorOrder", "regions", "threads"
105};
106
107char *USAGE =
108 "Usage: sddshist [<inputfile>] [<outputfile>]\n"
109 " [-pipe=[input][,output]]\n"
110 " -dataColumn=<column-name>\n"
111 " [{\n"
112 " -bins=<number> |\n"
113 " -sizeOfBins=<value> |\n"
114 " -regions=filename=<filename>,position=<columnName>,name=<columnName>\n"
115 " }]\n"
116 " [-lowerLimit=<value>]\n"
117 " [-upperLimit=<value>]\n"
118 " [-expand=<factor>]\n"
119 " [-filter=<column-name>,<lower-limit>,<upper-limit>]\n"
120 " [-weightColumn=<column-name>]\n"
121 " [-sides[=<points>]]\n"
122 " [-normalize[={sum|area|peak}]]\n"
123 " [-cdf[=only]]\n"
124 " [-threads=<number>]\n"
125 " [-statistics]\n"
126 " [-verbose]\n"
127 " [-majorOrder=row|column]\n"
128 "Options:\n"
129 " -pipe=[input][,output] Use pipe for input and/or output.\n"
130 " -dataColumn=<column-name> Specify the column to histogram.\n"
131 " -bins=<number> Set the number of bins for the histogram.\n"
132 " -sizeOfBins=<value> Set the size of each bin.\n"
133 " -regions=filename=<filename>,position=<columnName>,name=<columnName>\n"
134 " Define region-based histogramming.\n"
135 " -lowerLimit=<value> Set the lower limit of the histogram.\n"
136 " -upperLimit=<value> Set the upper limit of the histogram.\n"
137 " -expand=<factor> Expand the range of the histogram by the given factor.\n"
138 " -filter=<column-name>,<lower>,<upper> Filter data points based on column values.\n"
139 " -weightColumn=<column-name> Weight the histogram with the specified column.\n"
140 " -sides[=<points>] Add sides to the histogram down to zero level.\n"
141 " -normalize[={sum|area|peak}] Normalize the histogram.\n"
142 " -cdf[=only] Include the CDF in the output. Use 'only' to exclude the histogram.\n"
143 " -threads=<number> Specify the number of threads to use.\n"
144 " -statistics Include statistical information in the output.\n"
145 " -verbose Enable informational printouts during processing.\n"
146 " -majorOrder=row|column Set the major order of data.\n\n"
147 "Program by Michael Borland. (" __DATE__ " " __TIME__ ", SVN revision: " SVN_VERSION ")\n";
148
149#define NORMALIZE_PEAK 0
150#define NORMALIZE_AREA 1
151#define NORMALIZE_SUM 2
152#define NORMALIZE_NO 3
153#define N_NORMALIZE_OPTIONS 4
154char *normalize_option[N_NORMALIZE_OPTIONS] = {
155 "peak", "area", "sum", "no"
156};
157
158static int64_t filter(double *x, double *y, double *filterData, int64_t npts, double lower_filter, double upper_filter);
159static long setupOutputFile(SDDS_DATASET *outTable, char *outputfile, SDDS_DATASET *inTable, char *inputfile,
160 char *dataColumn, char *weightColumn, char *filterColumn, double lowerFilter,
161 double upperFilter, SDDS_DATASET *regionTable, char *regionNameColumn, long doStats,
162 int64_t bins, double binSize, long normalizeMode, short columnMajorOrder);
163
164/* Column and parameter indices for output file */
165static long iIndep, iFreq, iBins, iBinSize, iLoFilter, iUpFilter, iMean, iRMS, iStDev, iPoints, iCdf;
166static short cdfOnly, freOnly;
167
168int64_t readRegionFile(SDDS_DATASET *SDDSin, char *filename, char *positionColumn, char *nameColumn, double **regionPosition, char ***regionName);
169void classifyByRegion(double *data, double *weight, int64_t points, double *histogram, double *regionPosition, int64_t bins);
170static long make_histogram_threaded(double *hist, long n_bins, double lo, double hi, double *data,
171 int64_t n_pts, long new_start, int threads);
172static long make_histogram_weighted_threaded(double *hist, long n_bins, double lo, double hi, double *data,
173 int64_t n_pts, long new_start, double *weight, int threads);
174
175int main(int argc, char **argv) {
176 /* Flags to keep track of what is set in command line */
177 long binsGiven, lowerLimitGiven, upperLimitGiven;
178 SDDS_DATASET inTable, outTable;
179 double *data; /* Pointer to the array to histogram */
180 double *filterData; /* Pointer to the filter data */
181 double *weightData; /* Pointer to the weight data */
182 double *hist, *hist1; /* To store the histogram */
183 double *CDF, *CDF1; /* To store the CDF */
184 double sum; /* Total of histogram */
185 double *indep; /* Values of bin centers */
186 double lowerLimit, upperLimit; /* Lower and upper limits in histogram */
187 double givenLowerLimit, givenUpperLimit; /* Given lower and upper limits */
188 double range, binSize;
189 int64_t bins; /* Number of bins in the histogram */
190 long doStats; /* Include statistics in output file */
191 double mean, rms, standDev, mad;
192 char *filterColumn, *dataColumn, *weightColumn;
193 double lowerFilter = 0, upperFilter = 0; /* Filter range */
194 int64_t points; /* Number of data points after filtering */
195 SCANNED_ARG *scanned; /* Scanned argument structure */
196 char *inputfile, *outputfile; /* Input and output filenames */
197 double dx; /* Spacing of bins in histogram */
198 int64_t i; /* Loop variable */
199 long pointsBinned; /* Number of points in histogram */
200 long normalizeMode, doSides, verbose, readCode;
201 int64_t rows;
202 unsigned long pipeFlags, majorOrderFlag, regionFlags = 0;
203 char *cdf;
204 double expansionFactor = 0;
205 short columnMajorOrder = -1;
206 char *regionFilename = NULL, *regionPositionColumn = NULL, *regionNameColumn = NULL;
207 double *regionPosition = NULL;
208 int64_t nRegions = 0;
209 char **regionName = NULL;
210 SDDS_DATASET SDDSregion;
211 int threads = 1;
212
214 argc = scanargs(&scanned, argc, argv);
215 if (argc < 3) {
216 fprintf(stderr, "%s\n", USAGE);
217 exit(EXIT_FAILURE);
218 }
219
220 binsGiven = lowerLimitGiven = upperLimitGiven = 0;
221 binSize = doSides = 0;
222 inputfile = outputfile = NULL;
223 dataColumn = filterColumn = weightColumn = NULL;
224 doStats = verbose = 0;
225 normalizeMode = NORMALIZE_NO;
226 pipeFlags = 0;
227 dx = 0;
228 cdfOnly = 0;
229 freOnly = 1;
230
231 for (i = 1; i < argc; i++) {
232 if (scanned[i].arg_type == OPTION) {
233 switch (match_string(scanned[i].list[0], option, N_OPTIONS, 0)) {
234 case SET_MAJOR_ORDER:
235 majorOrderFlag = 0;
236 scanned[i].n_items--;
237 if (scanned[i].n_items > 0 && (!scanItemList(&majorOrderFlag, scanned[i].list + 1, &scanned[i].n_items, 0, "row", -1, NULL, 0, SDDS_ROW_MAJOR_ORDER, "column", -1, NULL, 0, SDDS_COLUMN_MAJOR_ORDER, NULL)))
238 SDDS_Bomb("invalid -majorOrder syntax/values");
239 if (majorOrderFlag & SDDS_COLUMN_MAJOR_ORDER)
240 columnMajorOrder = 1;
241 else if (majorOrderFlag & SDDS_ROW_MAJOR_ORDER)
242 columnMajorOrder = 0;
243 break;
244 case SET_BINS: /* Set number of bins */
245 if (binsGiven)
246 SDDS_Bomb("-bins specified more than once");
247 binsGiven = 1;
248 if (sscanf(scanned[i].list[1], "%" SCNd64, &bins) != 1 || bins <= 0)
249 SDDS_Bomb("invalid value for bins");
250 break;
251 case SET_LOWERLIMIT:
252 if (lowerLimitGiven)
253 SDDS_Bomb("-lowerLimit specified more than once");
254 lowerLimitGiven = 1;
255 if (sscanf(scanned[i].list[1], "%lf", &givenLowerLimit) != 1)
256 SDDS_Bomb("invalid value for lowerLimit");
257 break;
258 case SET_UPPERLIMIT:
259 if (upperLimitGiven)
260 SDDS_Bomb("-upperLimit specified more than once");
261 upperLimitGiven = 1;
262 if (sscanf(scanned[i].list[1], "%lf", &givenUpperLimit) != 1)
263 SDDS_Bomb("invalid value for upperLimit");
264 break;
265 case SET_EXPAND:
266 expansionFactor = 0;
267 if (sscanf(scanned[i].list[1], "%lf", &expansionFactor) != 1 || expansionFactor <= 0)
268 SDDS_Bomb("invalid value for expand");
269 break;
270 case SET_DATACOLUMN:
271 if (dataColumn)
272 SDDS_Bomb("-dataColumn specified more than once");
273 if (scanned[i].n_items != 2)
274 SDDS_Bomb("invalid -dataColumn syntax---supply name");
275 dataColumn = scanned[i].list[1];
276 break;
277 case SET_FILTER:
278 if (filterColumn)
279 SDDS_Bomb("multiple filter specifications not allowed");
280 if (scanned[i].n_items != 4 || sscanf(scanned[i].list[2], "%lf", &lowerFilter) != 1 ||
281 sscanf(scanned[i].list[3], "%lf", &upperFilter) != 1 || lowerFilter > upperFilter)
282 SDDS_Bomb("invalid -filter syntax/values");
283 filterColumn = scanned[i].list[1];
284 break;
285 case SET_WEIGHTCOLUMN:
286 if (weightColumn)
287 SDDS_Bomb("multiple weighting columns not allowed");
288 if (scanned[i].n_items != 2)
289 SDDS_Bomb("-weightColumn requires a column name");
290 weightColumn = scanned[i].list[1];
291 break;
292 case SET_NORMALIZE:
293 if (scanned[i].n_items == 1)
294 normalizeMode = NORMALIZE_SUM;
295 else if (scanned[i].n_items != 2 || (normalizeMode = match_string(scanned[i].list[1], normalize_option, N_NORMALIZE_OPTIONS, 0)) < 0)
296 SDDS_Bomb("invalid -normalize syntax");
297 break;
298 case SET_STATISTICS:
299 doStats = 1;
300 break;
301 case SET_SIDES:
302 if (scanned[i].n_items == 1)
303 doSides = 1;
304 else if (scanned[i].n_items > 2 || (sscanf(scanned[i].list[1], "%ld", &doSides) != 1 || doSides <= 0))
305 SDDS_Bomb("invalid -sides syntax");
306 break;
307 case SET_VERBOSE:
308 verbose = 1;
309 break;
310 case SET_BINSIZE:
311 if (sscanf(scanned[i].list[1], "%le", &binSize) != 1 || binSize <= 0)
312 SDDS_Bomb("invalid value for bin size");
313 break;
314 case SET_PIPE:
315 if (!processPipeOption(scanned[i].list + 1, scanned[i].n_items - 1, &pipeFlags))
316 SDDS_Bomb("invalid -pipe syntax");
317 break;
318 case SET_CDF:
319 if (scanned[i].n_items == 1)
320 cdfOnly = 0;
321 else {
322 if (scanned[i].n_items != 2)
323 SDDS_Bomb("invalid -cdf syntax");
324 cdf = scanned[i].list[1];
325 if (strcmp(cdf, "only") != 0)
326 SDDS_Bomb("invalid -cdf value, it should be -cdf or -cdf=only");
327 cdfOnly = 1;
328 }
329 freOnly = 0;
330 break;
331 case SET_REGION_FILE:
332 if (scanned[i].n_items != 4)
333 SDDS_Bomb("invalid -regionFile syntax");
334 regionFlags = 0;
335 scanned[i].n_items -= 1;
336 if (!scanItemList(&regionFlags, scanned[i].list + 1, &scanned[i].n_items, 0,
337 "filename", SDDS_STRING, &regionFilename, 1, 1,
338 "position", SDDS_STRING, &regionPositionColumn, 1, 2,
339 "name", SDDS_STRING, &regionNameColumn, 1, 4, NULL) ||
340 regionFlags != (1 + 2 + 4) || !regionFilename || !regionPositionColumn || !regionNameColumn)
341 SDDS_Bomb("invalid -regionFile syntax");
342 break;
343 case SET_THREADS:
344 if (scanned[i].n_items != 2 ||
345 !sscanf(scanned[i].list[1], "%d", &threads) || threads < 1)
346 SDDS_Bomb("invalid -threads syntax");
347 break;
348 default:
349 fprintf(stderr, "Error: option %s not recognized\n", scanned[i].list[0]);
350 exit(EXIT_FAILURE);
351 break;
352 }
353 } else {
354 /* Argument is filename */
355 if (!inputfile)
356 inputfile = scanned[i].list[0];
357 else if (!outputfile)
358 outputfile = scanned[i].list[0];
359 else
360 SDDS_Bomb("too many filenames seen");
361 }
362 }
363
364 processFilenames("sddshist", &inputfile, &outputfile, pipeFlags, 0, NULL);
365
366 if (binSize && binsGiven && regionFlags)
367 SDDS_Bomb("Provide only one of -bins, -sizeOfBins, or -regions");
368 if (!binsGiven)
369 bins = 20;
370 if (!dataColumn)
371 SDDS_Bomb("-dataColumn must be specified");
372
373 if (regionFlags) {
374 if (!(nRegions = readRegionFile(&SDDSregion, regionFilename, regionPositionColumn, regionNameColumn, &regionPosition, &regionName)))
375 SDDS_Bomb("Problem with region file. Check existence and type of columns");
376 doSides = 0;
377 bins = nRegions + 1;
378 }
379
380 hist = tmalloc(sizeof(*hist) * (bins + 2 * doSides));
381 CDF = CDF1 = tmalloc(sizeof(*hist) * (bins + 2 * doSides));
382 indep = tmalloc(sizeof(*indep) * (bins + 2 * doSides));
383 pointsBinned = 0;
384
385 if (!SDDS_InitializeInput(&inTable, inputfile) ||
386 SDDS_GetColumnIndex(&inTable, dataColumn) < 0 ||
387 (weightColumn && SDDS_GetColumnIndex(&inTable, weightColumn) < 0) ||
388 (filterColumn && SDDS_GetColumnIndex(&inTable, filterColumn) < 0))
389 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
390 if (!setupOutputFile(&outTable, outputfile, &inTable, inputfile, dataColumn, weightColumn, filterColumn, lowerFilter, upperFilter, &SDDSregion, regionNameColumn, doStats, bins, binSize, normalizeMode, columnMajorOrder))
391 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
392
393 data = weightData = filterData = NULL;
394 while ((readCode = SDDS_ReadPage(&inTable)) > 0) {
395 if ((rows = SDDS_CountRowsOfInterest(&inTable)) < 0)
396 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
397 if (rows && (!(data = SDDS_GetColumnInDoubles(&inTable, dataColumn)) ||
398 (weightColumn && !(weightData = SDDS_GetColumnInDoubles(&inTable, weightColumn))) ||
399 (filterColumn && !(filterData = SDDS_GetColumnInDoubles(&inTable, filterColumn)))))
400 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
401
402 if (rows && filterColumn)
403 points = filter(data, weightData, filterData, rows, lowerFilter, upperFilter);
404 else
405 points = rows;
406
407 pointsBinned = 0;
408 if (points) {
409 if (doStats) {
410 if (!weightColumn)
411 computeMomentsThreaded(&mean, &rms, &standDev, &mad, data, points, threads);
412 else
413 computeWeightedMomentsThreaded(&mean, &rms, &standDev, &mad, data, weightData, points, threads);
414 }
415
416 if (regionFlags) {
417 classifyByRegion(data, weightData, points, hist, regionPosition, bins);
418 hist1 = hist;
419 } else {
420 if (!lowerLimitGiven) {
421 lowerLimit = (points > 0) ? data[0] : 0;
422 for (i = 0; i < points; i++)
423 if (lowerLimit > data[i])
424 lowerLimit = data[i];
425 } else {
426 lowerLimit = givenLowerLimit;
427 }
428 if (!upperLimitGiven) {
429 upperLimit = (points > 0) ? data[0] : 0;
430 for (i = 0; i < points; i++)
431 if (upperLimit < data[i])
432 upperLimit = data[i];
433 } else {
434 upperLimit = givenUpperLimit;
435 }
436
437 range = upperLimit - lowerLimit;
438 if (!lowerLimitGiven)
439 lowerLimit -= range * 1e-7;
440 if (!upperLimitGiven)
441 upperLimit += range * 1e-7;
442 if (upperLimit == lowerLimit) {
443 if (binSize) {
444 upperLimit += binSize / 2;
445 lowerLimit -= binSize / 2;
446 } else if (fabs(upperLimit) < sqrt(DBL_MIN)) {
447 upperLimit = sqrt(DBL_MIN);
448 lowerLimit = -sqrt(DBL_MIN);
449 } else {
450 upperLimit += upperLimit * (1 + 2 * DBL_EPSILON);
451 lowerLimit -= upperLimit * (1 - 2 * DBL_EPSILON);
452 }
453 }
454 if (expansionFactor > 0) {
455 double center = (upperLimit + lowerLimit) / 2;
456 range = expansionFactor * (upperLimit - lowerLimit);
457 lowerLimit = center - range / 2;
458 upperLimit = center + range / 2;
459 }
460 dx = (upperLimit - lowerLimit) / bins;
461
462 if (binSize) {
463 double middle;
464 range = ((range / binSize) + 1) * binSize;
465 middle = (lowerLimit + upperLimit) / 2;
466 lowerLimit = middle - range / 2;
467 upperLimit = middle + range / 2;
468 dx = binSize;
469 bins = range / binSize + 0.5;
470 if (bins < 1 && !doSides)
471 bins = 2 * doSides;
472 indep = trealloc(indep, sizeof(*indep) * (bins + 2 * doSides));
473 hist = trealloc(hist, sizeof(*hist) * (bins + 2 * doSides));
474 CDF = trealloc(CDF, sizeof(*hist) * (bins + 2 * doSides));
475 }
476
477 for (i = -doSides; i < bins + doSides; i++)
478 indep[i + doSides] = (i + 0.5) * dx + lowerLimit;
479 hist1 = hist + doSides;
480 CDF1 = CDF + doSides;
481 if (doSides) {
482 hist[0] = hist[bins + doSides] = 0;
483 }
484
485 if (!weightColumn)
486 pointsBinned = make_histogram_threaded(hist1, bins, lowerLimit, upperLimit, data, points, 1, threads);
487 else
488 pointsBinned = make_histogram_weighted_threaded(hist1, bins, lowerLimit, upperLimit, data, points, 1, weightData, threads);
489 }
490
491 sum = 0;
492 for (i = 0; i < bins + doSides; i++) {
493 sum += hist1[i];
494 }
495 CDF1[0] = hist1[0] / sum;
496 for (i = 1; i < bins + doSides; i++) {
497 CDF1[i] = CDF1[i - 1] + hist1[i] / sum;
498 }
499
500 if (verbose)
501 fprintf(stderr, "%ld points of %" PRId64 " from page %ld histogrammed in %" PRId64 " bins\n", pointsBinned, rows, readCode, bins);
502 if (!cdfOnly) {
503 if (normalizeMode != NORMALIZE_NO) {
504 double norm = 0;
505 switch (normalizeMode) {
506 case NORMALIZE_PEAK:
507 norm = max_in_array(hist1, bins);
508 break;
509 case NORMALIZE_AREA:
510 case NORMALIZE_SUM:
511 for (i = 0; i < bins; i++)
512 norm += hist1[i];
513 if (normalizeMode == NORMALIZE_AREA)
514 norm *= dx;
515 break;
516 default:
517 SDDS_Bomb("invalid normalize mode--consult programmer.");
518 break;
519 }
520 if (norm)
521 for (i = 0; i < bins; i++)
522 hist1[i] /= norm;
523 }
524 }
525 }
526
527 if (regionFlags) {
528 if (!SDDS_StartPage(&outTable, bins) ||
529 !SDDS_CopyParameters(&outTable, &inTable) ||
530 !SDDS_SetParameters(&outTable, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE, iBins, bins, iBinSize, dx, iPoints, pointsBinned, -1))
531 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
532 if (points) {
533 if (!SDDS_SetColumn(&outTable, SDDS_SET_BY_INDEX, regionPosition, bins, iIndep) ||
534 !SDDS_SetColumn(&outTable, SDDS_SET_BY_NAME, regionName, bins, regionNameColumn))
535 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
536 if (!freOnly && !SDDS_SetColumn(&outTable, SDDS_SET_BY_INDEX, CDF, bins, iCdf))
537 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
538 if (!cdfOnly && !SDDS_SetColumn(&outTable, SDDS_SET_BY_INDEX, hist, bins, iFreq))
539 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
540 }
541 } else {
542 if (!SDDS_StartPage(&outTable, bins + 2 * doSides) ||
543 !SDDS_CopyParameters(&outTable, &inTable) ||
544 (points && (!SDDS_SetColumn(&outTable, SDDS_SET_BY_INDEX, indep, bins + 2 * doSides, iIndep))) ||
545 !SDDS_SetParameters(&outTable, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE, iBins, bins, iBinSize, dx, iPoints, pointsBinned, -1))
546 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
547 if (!freOnly) {
548 if (points && !SDDS_SetColumn(&outTable, SDDS_SET_BY_INDEX, CDF, bins + 2 * doSides, iCdf))
549 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
550 }
551 if (!cdfOnly) {
552 if (points && !SDDS_SetColumn(&outTable, SDDS_SET_BY_INDEX, hist, bins + 2 * doSides, iFreq))
553 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
554 }
555 }
556
557 if (filterColumn && points &&
558 !SDDS_SetParameters(&outTable, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE, iLoFilter, lowerFilter, iUpFilter, upperFilter, -1))
559 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
560 if (doStats && points &&
561 !SDDS_SetParameters(&outTable, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE, iMean, mean, iRMS, rms, iStDev, standDev, -1))
562 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
563
564 if (!SDDS_WritePage(&outTable))
565 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
566 if (data)
567 free(data);
568 if (weightData)
569 free(weightData);
570 if (filterData)
571 free(filterData);
572 data = weightData = filterData = NULL;
573 }
574
575 if (!SDDS_Terminate(&inTable)) {
576 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
577 exit(EXIT_FAILURE);
578 }
579 if (!SDDS_Terminate(&outTable)) {
580 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
581 exit(EXIT_FAILURE);
582 }
583 return EXIT_SUCCESS;
584}
585
586/* routine: filter()
587 * purpose: filter set of points {(x, y)} with a specified
588 * window in one of the variables. Works in place.
589 */
590static int64_t filter(double *x, double *y, double *filterData, int64_t npts, double lower_filter, double upper_filter) {
591 int64_t i, j;
592 static char *keep = NULL;
593 static int64_t maxPoints = 0;
594
595 if (maxPoints < npts)
596 keep = trealloc(keep, sizeof(*keep) * (maxPoints = npts));
597
598 for (i = 0; i < npts; i++) {
599 if (filterData[i] < lower_filter || filterData[i] > upper_filter)
600 keep[i] = 0;
601 else
602 keep[i] = 1;
603 }
604
605 for (i = j = 0; i < npts; i++)
606 if (keep[i]) {
607 if (i != j) {
608 if (x)
609 x[j] = x[i];
610 if (y)
611 y[j] = y[i];
612 filterData[j] = filterData[i];
613 }
614 j++;
615 }
616
617 return j;
618}
619
620static long setupOutputFile(SDDS_DATASET *outTable, char *outputfile, SDDS_DATASET *inTable, char *inputfile,
621 char *dataColumn, char *weightColumn, char *filterColumn, double lowerFilter,
622 double upperFilter, SDDS_DATASET *regionTable, char *regionNameColumn, long doStats,
623 int64_t bins, double binSize, long normalizeMode, short columnMajorOrder) {
624 char *symbol, *units, *dataUnits, *outputFormat;
625 int32_t outputType;
626 char s[1024];
627
628 if (!SDDS_InitializeOutput(outTable, SDDS_BINARY, 0, NULL, "sddshist output", outputfile))
629 return 0;
630 if (columnMajorOrder != -1)
631 outTable->layout.data_mode.column_major = columnMajorOrder;
632 else
633 outTable->layout.data_mode.column_major = inTable->layout.data_mode.column_major;
634 if (!SDDS_GetColumnInformation(inTable, "units", &dataUnits, SDDS_GET_BY_NAME, dataColumn))
635 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
636
637 /* Define output columns */
638 outputType = SDDS_DOUBLE;
639 outputFormat = NULL;
640
641 if (!SDDS_TransferColumnDefinition(outTable, inTable, dataColumn, NULL) ||
642 !SDDS_ChangeColumnInformation(outTable, "type", &outputType, SDDS_BY_NAME, dataColumn) ||
643 !SDDS_ChangeColumnInformation(outTable, "format_string", &outputFormat, SDDS_BY_NAME, dataColumn) ||
644 (iIndep = SDDS_GetColumnIndex(outTable, dataColumn)) < 0)
645 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
646 if (regionNameColumn && !SDDS_TransferColumnDefinition(outTable, regionTable, regionNameColumn, NULL))
647 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
648 if (!cdfOnly) {
649 switch (normalizeMode) {
650 case NORMALIZE_PEAK:
651 symbol = "RelativeFrequency";
652 units = NULL;
653 break;
654 case NORMALIZE_AREA:
655 symbol = "NormalizedFrequency";
656 if (dataUnits && !SDDS_StringIsBlank(dataUnits)) {
657 units = tmalloc(sizeof(*units) * (strlen(dataUnits) + 5));
658 if (strchr(dataUnits, ' '))
659 sprintf(units, "1/(%s)", dataUnits);
660 else
661 sprintf(units, "1/%s", dataUnits);
662 } else
663 units = NULL;
664 break;
665 case NORMALIZE_SUM:
666 symbol = "FractionalFrequency";
667 units = NULL;
668 break;
669 default:
670 if (weightColumn) {
671 char *weightUnits = NULL;
672 if (weightColumn && !SDDS_GetColumnInformation(inTable, "units", &weightUnits, SDDS_GET_BY_NAME, weightColumn))
673 return 0;
674 symbol = "WeightedNumberOfOccurrences";
675 units = weightUnits;
676 } else {
677 symbol = "NumberOfOccurrences";
678 units = NULL;
679 }
680 break;
681 }
682
683 if ((iFreq = SDDS_DefineColumn(outTable, "frequency", symbol, units, NULL, NULL, SDDS_DOUBLE, 0)) < 0)
684 return 0;
685 free(units);
686 units = NULL;
687 }
688 if (!freOnly) {
689 sprintf(s, "%sCdf", dataColumn);
690 if ((iCdf = SDDS_DefineColumn(outTable, s, NULL, NULL, NULL, NULL, SDDS_DOUBLE, 0)) < 0)
691 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
692 }
693
694 /* Define output parameters */
695 if (SDDS_DefineParameter(outTable, "sddshistInput", NULL, NULL, NULL, NULL, SDDS_STRING,
696 inputfile) < 0 ||
697 (weightColumn && SDDS_DefineParameter(outTable, "sddshistWeight", NULL, NULL, NULL, NULL,
698 SDDS_STRING, weightColumn) < 0) ||
699 (iBins = SDDS_DefineParameter(outTable, "sddshistBins", NULL, NULL, NULL, NULL, SDDS_LONG,
700 NULL)) < 0 ||
701 (iBinSize = SDDS_DefineParameter(outTable, "sddshistBinSize", NULL, NULL, NULL, NULL, SDDS_DOUBLE, NULL)) < 0 ||
702 (iPoints = SDDS_DefineParameter(outTable, "sddshistBinned", NULL, NULL, NULL, NULL, SDDS_LONG, NULL)) < 0)
703 return 0;
704 if (filterColumn) {
705 char *filterUnits;
706 if (!SDDS_GetColumnInformation(inTable, "units", &filterUnits, SDDS_GET_BY_NAME, filterColumn) ||
707 SDDS_DefineParameter(outTable, "sddshistFilter", NULL, NULL, NULL, NULL, SDDS_STRING,
708 filterColumn) < 0 ||
709 (iLoFilter = SDDS_DefineParameter(outTable, "sddshistLowerFilter", NULL, filterUnits, NULL, NULL, SDDS_DOUBLE, NULL)) < 0 ||
710 (iUpFilter = SDDS_DefineParameter(outTable, "sddshistUpperFilter", NULL, filterUnits, NULL, NULL, SDDS_DOUBLE, NULL)) < 0)
711 return 0;
712 if (filterUnits)
713 free(filterUnits);
714 }
715 if (doStats) {
716 char *buffer;
717 buffer = tmalloc(sizeof(*buffer) * (strlen(dataColumn) + 20));
718 sprintf(buffer, "%sMean", dataColumn);
719 if ((iMean = SDDS_DefineParameter(outTable, buffer, NULL, dataUnits, NULL, NULL, SDDS_DOUBLE, NULL)) < 0)
720 return 0;
721 sprintf(buffer, "%sRms", dataColumn);
722 if ((iRMS = SDDS_DefineParameter(outTable, buffer, NULL, dataUnits, NULL, NULL, SDDS_DOUBLE, NULL)) < 0)
723 return 0;
724 sprintf(buffer, "%sStDev", dataColumn);
725 if ((iStDev = SDDS_DefineParameter(outTable, buffer, NULL, dataUnits, NULL, NULL, SDDS_DOUBLE, NULL)) < 0)
726 return 0;
727 free(buffer);
728 }
729 if (SDDS_DefineParameter(outTable, "sddshistNormMode", NULL, NULL, NULL, NULL, SDDS_STRING, normalize_option[normalizeMode]) < 0 ||
730 !SDDS_TransferAllParameterDefinitions(outTable, inTable, SDDS_TRANSFER_KEEPOLD) ||
731 !SDDS_WriteLayout(outTable))
732 return 0;
733 return 1;
734}
735
736int64_t readRegionFile(SDDS_DATASET *SDDSin, char *filename, char *positionColumn, char *nameColumn, double **regionPosition, char ***regionName) {
737 int64_t i, rows = 0;
738 if (!SDDS_InitializeInput(SDDSin, filename) ||
739 SDDS_ReadPage(SDDSin) != 1 || (rows = SDDS_RowCount(SDDSin)) < 1 ||
740 !(*regionPosition = SDDS_GetColumnInDoubles(SDDSin, positionColumn)) ||
741 !(*regionName = (char **)SDDS_GetColumn(SDDSin, nameColumn)))
742 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
743 for (i = 1; i < rows; i++)
744 if ((*regionPosition)[i] <= (*regionPosition)[i - 1]) {
745 fprintf(stderr, "sddshist: Error in region position data: row %" PRId64 " is %21.15e while row %" PRId64 " is %21.15e\n",
746 i - 1, (*regionPosition)[i - 1], i, (*regionPosition)[i]);
747 exit(EXIT_FAILURE);
748 }
749 *regionPosition = SDDS_Realloc(*regionPosition, sizeof(**regionPosition) * (rows + 1));
750 (*regionPosition)[rows] = DBL_MAX;
751 *regionName = SDDS_Realloc(*regionName, sizeof(**regionName) * (rows + 1));
752 cp_str(&(*regionName)[rows], "Beyond");
753 return rows;
754}
755
756void classifyByRegion(double *data, double *weight, int64_t points, double *histogram, double *regionPosition, int64_t bins) {
757 int64_t iData, iBin;
758
759 for (iBin = 0; iBin < bins; iBin++)
760 histogram[iBin] = 0;
761
762 for (iData = 0; iData < points; iData++) {
763 /* Note that bins is 1 greater than the number of positions */
764 for (iBin = 0; iBin < bins - 1; iBin++) {
765 if (data[iData] < regionPosition[iBin])
766 break;
767 }
768 if (weight) {
769 histogram[iBin] += weight[iData];
770 } else {
771 histogram[iBin] += 1;
772 }
773 }
774}
775
776static long make_histogram_threaded(double *hist, long n_bins, double lo, double hi, double *data,
777 int64_t n_pts, long new_start, int threads) {
778 double bin_size, *partial;
779 long *counts, count = 0;
780 int activeThreads, thread;
781 int64_t i;
782
783 if (threads <= 1 || n_pts <= 0 || n_bins <= 0)
784 return make_histogram(hist, n_bins, lo, hi, data, n_pts, new_start);
785
786 activeThreads = threads;
787 if (activeThreads > n_pts)
788 activeThreads = (int)n_pts;
789 if (activeThreads <= 1)
790 return make_histogram(hist, n_bins, lo, hi, data, n_pts, new_start);
791
792 if (new_start)
793 for (i = 0; i < n_bins; i++)
794 hist[i] = 0;
795 bin_size = (hi - lo) / n_bins;
796 partial = calloc((size_t)activeThreads * n_bins, sizeof(*partial));
797 counts = calloc(activeThreads, sizeof(*counts));
798 if (!partial || !counts)
799 SDDS_Bomb("memory allocation failure");
800
801#pragma omp parallel for if (activeThreads > 1) num_threads(activeThreads)
802 for (thread = 0; thread < activeThreads; thread++) {
803 int64_t start = thread * (n_pts / activeThreads);
804 int64_t end = (thread == activeThreads - 1) ? n_pts : (thread + 1) * (n_pts / activeThreads);
805 double *local = partial + (size_t)thread * n_bins;
806 long localCount = 0;
807 for (int64_t point = start; point < end; point++) {
808 double dbin = (data[point] - lo) / bin_size;
809 long bin = dbin;
810 if (dbin < 0)
811 continue;
812 if (bin < 0 || bin >= n_bins)
813 continue;
814 local[bin] += 1;
815 localCount++;
816 }
817 counts[thread] = localCount;
818 }
819
820 for (thread = 0; thread < activeThreads; thread++) {
821 double *local = partial + (size_t)thread * n_bins;
822 count += counts[thread];
823 for (i = 0; i < n_bins; i++)
824 hist[i] += local[i];
825 }
826 free(partial);
827 free(counts);
828 return count;
829}
830
831static long make_histogram_weighted_threaded(double *hist, long n_bins, double lo, double hi, double *data,
832 int64_t n_pts, long new_start, double *weight, int threads) {
833 double bin_size, *partial;
834 long *counts, count = 0;
835 int activeThreads, thread;
836 int64_t i;
837
838 if (threads <= 1 || n_pts <= 0 || n_bins <= 0)
839 return make_histogram_weighted(hist, n_bins, lo, hi, data, n_pts, new_start, weight);
840
841 activeThreads = threads;
842 if (activeThreads > n_pts)
843 activeThreads = (int)n_pts;
844 if (activeThreads <= 1)
845 return make_histogram_weighted(hist, n_bins, lo, hi, data, n_pts, new_start, weight);
846
847 if (new_start)
848 for (i = 0; i < n_bins; i++)
849 hist[i] = 0;
850 bin_size = (hi - lo) / n_bins;
851 partial = calloc((size_t)activeThreads * n_bins, sizeof(*partial));
852 counts = calloc(activeThreads, sizeof(*counts));
853 if (!partial || !counts)
854 SDDS_Bomb("memory allocation failure");
855
856#pragma omp parallel for if (activeThreads > 1) num_threads(activeThreads)
857 for (thread = 0; thread < activeThreads; thread++) {
858 int64_t start = thread * (n_pts / activeThreads);
859 int64_t end = (thread == activeThreads - 1) ? n_pts : (thread + 1) * (n_pts / activeThreads);
860 double *local = partial + (size_t)thread * n_bins;
861 long localCount = 0;
862 for (int64_t point = start; point < end; point++) {
863 double dbin = (data[point] - lo) / bin_size;
864 long bin = dbin;
865 if (dbin < 0)
866 continue;
867 if (bin < 0 || bin >= n_bins)
868 continue;
869 local[bin] += weight[point];
870 localCount++;
871 }
872 counts[thread] = localCount;
873 }
874
875 for (thread = 0; thread < activeThreads; thread++) {
876 double *local = partial + (size_t)thread * n_bins;
877 count += counts[thread];
878 for (i = 0; i < n_bins; i++)
879 hist[i] += local[i];
880 }
881 free(partial);
882 free(counts);
883 return count;
884}
SDDS (Self Describing Data Set) Data Types Definitions and Function Prototypes.
int32_t SDDS_CopyParameters(SDDS_DATASET *SDDS_target, SDDS_DATASET *SDDS_source)
Definition SDDS_copy.c:286
int32_t SDDS_StartPage(SDDS_DATASET *SDDS_dataset, int64_t expected_n_rows)
int32_t SDDS_SetParameters(SDDS_DATASET *SDDS_dataset, int32_t mode,...)
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.
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".
int64_t SDDS_CountRowsOfInterest(SDDS_DATASET *SDDS_dataset)
Counts the number of rows marked as "of interest" in the current data table.
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...
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
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
int32_t SDDS_InitializeInput(SDDS_DATASET *SDDS_dataset, char *filename)
Definition SDDS_input.c:50
int32_t SDDS_Terminate(SDDS_DATASET *SDDS_dataset)
int32_t SDDS_ReadPage(SDDS_DATASET *SDDS_dataset)
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.
int32_t SDDS_WritePage(SDDS_DATASET *SDDS_dataset)
Writes the current data table to the output file.
int32_t SDDS_DefineColumn(SDDS_DATASET *SDDS_dataset, const char *name, const char *symbol, const char *units, const char *description, const char *format_string, int32_t type, int32_t field_length)
Defines a data column within the SDDS dataset.
int32_t SDDS_WriteLayout(SDDS_DATASET *SDDS_dataset)
Writes the SDDS layout header to the output file.
int32_t SDDS_DefineParameter(SDDS_DATASET *SDDS_dataset, const char *name, const char *symbol, const char *units, const char *description, const char *format_string, int32_t type, char *fixed_value)
Defines a data parameter with a fixed string value.
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.
int32_t SDDS_TransferAllParameterDefinitions(SDDS_DATASET *SDDS_target, SDDS_DATASET *SDDS_source, uint32_t mode)
Transfers all parameter definitions from a source dataset to a target dataset.
int32_t SDDS_GetColumnIndex(SDDS_DATASET *SDDS_dataset, char *name)
Retrieves the index of a named column in the SDDS dataset.
void SDDS_PrintErrors(FILE *fp, int32_t mode)
Prints recorded error messages to a specified file stream.
Definition SDDS_utils.c:474
void SDDS_RegisterProgramName(const char *name)
Registers the executable program name for use in error messages.
Definition SDDS_utils.c:318
int32_t SDDS_StringIsBlank(char *s)
Checks if a string is blank (contains only whitespace characters).
void SDDS_Bomb(char *message)
Terminates the program after printing an error message and recorded errors.
Definition SDDS_utils.c:380
void * SDDS_Realloc(void *old_ptr, size_t new_size)
Reallocates memory to a new size.
Definition SDDS_utils.c:743
#define SDDS_STRING
Identifier for the string data type.
Definition SDDStypes.h:85
#define SDDS_LONG
Identifier for the signed 32-bit integer data type.
Definition SDDStypes.h:61
#define SDDS_DOUBLE
Identifier for the double data type.
Definition SDDStypes.h:37
void * trealloc(void *old_ptr, uint64_t size_of_block)
Reallocates a memory block to a new size.
Definition array.c:190
void * tmalloc(uint64_t size_of_block)
Allocates a memory block of the specified size with zero initialization.
Definition array.c:65
char * cp_str(char **s, char *t)
Copies a string, allocating memory for storage.
Definition cp_str.c:28
double max_in_array(double *array, long n)
Finds the maximum value in an array of doubles.
Definition findMinMax.c:318
long make_histogram_weighted(double *hist, long n_bins, double lo, double hi, double *data, long n_pts, long new_start, double *weight)
Compiles a weighted histogram from data points.
long make_histogram(double *hist, long n_bins, double lo, double hi, double *data, int64_t n_pts, long new_start)
Compiles a histogram from data points.
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 computeWeightedMomentsThreaded(double *mean, double *rms, double *standDev, double *meanAbsoluteDev, double *x, double *w, long n, long numThreads)
Computes weighted statistical moments of an array using multiple threads.
Definition moments.c:240
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...
Definition moments.c:127
int scanargs(SCANNED_ARG **scanned, int argc, char **argv)
Definition scanargs.c:36
long processPipeOption(char **item, long items, unsigned long *flags)
Definition scanargs.c:357
void processFilenames(char *programName, char **input, char **output, unsigned long pipeFlags, long noWarnings, long *tmpOutputUsed)
Definition scanargs.c:391
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.