SDDS ToolKit Programs and Libraries for C and Python
Loading...
Searching...
No Matches
sddsmultihist.c
Go to the documentation of this file.
1/**
2 * @file sddsmultihist.c
3 * @brief SDDS-format multi-column histogramming program.
4 *
5 * @details
6 * This program generates histograms for specified columns of an SDDS-formatted input file, offering customization
7 * options for bin sizing, normalization, boundaries, and more.
8 *
9 * @section Usage
10 * ```
11 * sddsmultihist [<inputfile>] [<outputfile>]
12 * [-pipe=[input][,output]]
13 * -columns=<name>[,...]
14 * -abscissa=<name>[,...]
15 * [-exclude=<name>[,...]]
16 * [-bins=<integer>]
17 * [-sizeOfBins=<value>]
18 * [-autobins=target=<number>[,minimum=<integer>][,maximum=<integer>]]
19 * [-boundaryData=<filename>,<column>]
20 * [-sides[=close|against]]
21 * [-expand=<fraction>]
22 * [-lowerLimit=<value>[,...]]
23 * [-upperLimit=<value>[,...]]
24 * [-separate]
25 * [-cdf=[only]]
26 * [-weightColumn=<name>]
27 * [-majorOrder=row|column]
28 * [-normalize={sum|peak|no}]
29 * [-threads=<number>]
30 * ```
31 *
32 * @section Options
33 * | Required | Description |
34 * |---------------------------------------|---------------------------------------------------------------------------------------|
35 * | `-columns` | Names of columns from the input to be histogrammed. Supports wildcards. |
36 * | `-abscissa` | Names of abscissas in the output file. |
37 *
38 * | Optional | Description |
39 * |---------------------------------------|---------------------------------------------------------------------------------------|
40 * | `-pipe` | Use standard SDDS Toolkit pipe for input/output. |
41 * | `-exclude` | Column names to exclude from histogramming. |
42 * | `-bins` | Number of bins for the histogram. |
43 * | `-sizeOfBins` | Size of each bin for the histogram. |
44 * | `-autobins` | Auto-determine bin count based on target samples per bin. |
45 * | `-boundaryData` | Specifies irregular bin boundaries from a file. |
46 * | `-sides` | Adds zero-height bins at the histogram's ends. |
47 * | `-expand` | Expands histogram range by the given fraction. |
48 * | `-lowerLimit` | Sets lower limits for histograms. |
49 * | `-upperLimit` | Sets upper limits for histograms. |
50 * | `-separate` | Creates separate abscissas for each histogram. |
51 * | `-cdf` | Includes CDF in the output; "only" includes only CDF. |
52 * | `-weightColumn` | Specifies a column to weight the histogram. |
53 * | `-majorOrder` | Sets the output data order to row-major or column-major. |
54 * | `-normalize` | Normalizes the histogram based on sum, peak, or none. |
55 * | `-threads` | Number of threads for regular histogram binning. |
56 *
57 * @subsection Incompatibilities
58 * -boundaryData is incompatible with:
59 * - -separate
60 * - -abscissa
61 * - Only one of the following may be specified:
62 * - -bins
63 * - -sizeOfBins
64 * - -autobins
65 * - For -separate:
66 * - The number of -lowerLimit or -upperLimit values must match the number of columns.
67 *
68 * @copyright
69 * - (c) 2002 The University of Chicago, as Operator of Argonne National Laboratory.
70 * - (c) 2002 The Regents of the University of California, as Operator of Los Alamos National Laboratory.
71 *
72 * @license
73 * This file is distributed under the terms of the Software License Agreement
74 * found in the file LICENSE included with this distribution.
75 *
76 * @authors
77 * M. Borland, L. Emery, R. Soliday, H. Shang
78 */
79
80#include "mdb.h"
81#include "SDDS.h"
82#include "scan.h"
83#include "SDDSutils.h"
84#include <ctype.h>
85
86/* Enumeration for option types */
87enum option_type {
88 SET_COLUMNS,
89 SET_PIPE,
90 SET_EXCLUDE,
91 SET_ABSCISSA,
92 SET_BINS,
93 SET_SIZEOFBINS,
94 SET_LOWERLIMIT,
95 SET_UPPERLIMIT,
96 SET_SIDES,
97 SET_SEPARATE,
98 SET_EXPAND,
99 SET_CDF,
100 SET_AUTOBINS,
101 SET_MAJOR_ORDER,
102 SET_BOUNDARYDATA,
103 SET_WEIGHT,
104 SET_NORMALIZE,
105 SET_THREADS,
106 N_OPTIONS
107};
108
109char *option[N_OPTIONS] = {
110 "columns",
111 "pipe",
112 "exclude",
113 "abscissa",
114 "bins",
115 "sizeofbins",
116 "lowerlimit",
117 "upperlimit",
118 "sides",
119 "separate",
120 "expand",
121 "cdf",
122 "autobins",
123 "majorOrder",
124 "boundaryData",
125 "weightColumn",
126 "normalize",
127 "threads",
128};
129
130static char *USAGE =
131 "Usage: sddsmultihist [<inputfile>] [<outputfile>]\n"
132 " [-pipe=[input][,output]]\n"
133 " -columns=<name>[,...]\n"
134 " -abscissa=<name>[,...]\n"
135 " [-exclude=<name>[,...]]\n"
136 " [-bins=<integer>]\n"
137 " [-sizeOfBins=<value>]\n"
138 " [-autobins=target=<number>[,minimum=<integer>][,maximum=<integer>]]\n"
139 " [-boundaryData=<filename>,<column>]\n"
140 " [-sides[=close|against]]\n"
141 " [-expand=<fraction>]\n"
142 " [-lowerLimit=<value>[,...]]\n"
143 " [-upperLimit=<value>[,...]]\n"
144 " [-separate]\n"
145 " [-cdf=[only]]\n"
146 " [-weightColumn=<name>]\n"
147 " [-majorOrder=row|column]\n"
148 " [-normalize={sum|peak|no}]\n"
149 " [-threads=<number>]\n"
150 "Options:\n"
151 " -pipe=[input][,output] The standard SDDS Toolkit pipe option.\n"
152 " -columns=<name>[,...] Specifies the names of columns from the input to be histogrammed.\n"
153 " Names may contain wildcards.\n"
154 " -abscissa=<name>[,...] Specifies the names of the abscissas in the output file.\n"
155 " When using column names as abscissa names,\n"
156 " the -abscissa option is not required (use -separate).\n"
157 " At least one abscissa name must be supplied if -separate is not used.\n"
158 " -exclude=<name>[,...] (Optional) Specifies column names to exclude from histogramming.\n"
159 " -bins=<integer> Sets the number of bins for the histogram.\n"
160 " -sizeOfBins=<value> Sets the size of each bin for the histogram.\n"
161 " -autobins=target=<number>[,minimum=<integer>][,maximum=<integer>]\n"
162 " Automatically determines the number of bins based on the target number of samples per bin.\n"
163 " Optionally specify minimum and maximum number of bins.\n"
164 " -boundaryData=<filename>,<column> Specifies irregular bin boundaries from a file.\n"
165 " Incompatible with -separate and -abscissa.\n"
166 " -sides[=close|against] Adds zero-height bins at the ends of the histogram.\n"
167 " 'close' centers the first and last bins.\n"
168 " 'against' aligns the first and last bins with the data range.\n"
169 " -expand=<fraction> Expands the range of the histogram by the given fraction.\n"
170 " -lowerLimit=<value>[,...] Sets lower limits for the histograms.\n"
171 " -upperLimit=<value>[,...] Sets upper limits for the histograms.\n"
172 " -separate Creates separate abscissas for each histogram in the output file.\n"
173 " -cdf=[only] Includes the Cumulative Distribution Function (CDF) in the output.\n"
174 " 'only' includes only the CDF, excluding the histogram.\n"
175 " -weightColumn=<name> Specifies a column to weight the histogram.\n"
176 " -majorOrder=row|column Sets the output file's data order to row-major or column-major.\n"
177 " -normalize={sum|peak|no} Normalizes the histogram.\n"
178 " 'sum' normalizes so that the sum of all bins equals 1.\n"
179 " 'peak' normalizes so that the peak bin equals 1.\n"
180 " 'no' applies no normalization.\n"
181 " -threads=<number> Number of threads for regular histogram binning.\n"
182 "\n"
183 "Program by Michael Borland. (" __DATE__ " " __TIME__ ", SVN revision: " SVN_VERSION ")\n";
184
185#define DO_SIDES 1
186#define DO_CLOSE_SIDES 2
187#define DO_AGAINST_SIDES 3
188
189#define NORMALIZE_PEAK 0
190#define NORMALIZE_SUM 1
191#define NORMALIZE_NO 2
192#define N_NORMALIZE_OPTIONS 3
193char *normalize_option[N_NORMALIZE_OPTIONS] = {
194 "peak", "sum", "no"};
195
196void SetUpOutput(SDDS_DATASET *SDDSout, SDDS_DATASET *SDDSin, long **abscissaIndex, long **cdfIndex,
197 long **histogramIndex, char *output, char **columnName, long columnNames,
198 char **abscissaName, long abscissaNames, char *boundaryColumn, char *boundaryColumnUnits,
199 short columnMajorOrder, short normMode);
200double *ReadBoundaryData(char *file, char *column, int64_t *n, char **units);
201void MakeBoundaryHistogram(double *histogram, double *cdf, double *boundaryValue, int64_t nBoundaryValues,
202 double *data, double *weight, int64_t nData);
203void NormalizeHistogram(double *hist, int64_t bins, short mode);
204static long make_histogram_threaded(double *hist, long n_bins, double lo, double hi, double *data,
205 int64_t n_pts, long new_start, int threads);
206static long make_histogram_weighted_threaded(double *hist, long n_bins, double lo, double hi, double *data,
207 int64_t n_pts, long new_start, double *weight, int threads);
208
209static short cdfOnly, freOnly;
210
211int main(int argc, char **argv) {
212 int iArg;
213 char *boundaryFile, *boundaryColumn, *boundaryColumnUnits;
214 double *boundaryValue;
215 int64_t nBoundaryValues;
216 char **abscissaName, **columnName, **excludeName;
217 long columnNames, excludeNames, abscissaNames, column, offset;
218 long givenLowerLimits, givenUpperLimits;
219 char *input, *output;
220 long readCode, binsGiven;
221 int64_t i, rows, bins, writeBins;
222 long lowerLimitGiven, upperLimitGiven, doSides, doSeparate;
223 double autoBinsTarget;
224 long autoBinsMinimum, autoBinsMaximum;
225 char *weightColumn;
226 double *weightData;
227 short normMode = NORMALIZE_NO;
228 unsigned long pipeFlags;
229 SCANNED_ARG *scanned;
230 SDDS_DATASET SDDSin, SDDSout;
231 double binSize, *lowerLimit = NULL, *upperLimit = NULL, *givenLowerLimit, *givenUpperLimit, *dx = NULL, range, middle;
232 double **inputData, *abscissa, *histogram, *minValue, *maxValue, transferLimit, *cdf, sum;
233 long *abscissaIndex, *histogramIndex, *cdfIndex;
234 double expandRange, maxRange;
235 char *CDFONLY;
236 unsigned long dummyFlags, majorOrderFlag;
237 short columnMajorOrder = -1;
238 int threads = 1;
239
240 SDDS_RegisterProgramName(argv[0]);
241 argc = scanargs(&scanned, argc, argv);
242 if (argc < 3 || argc > (3 + N_OPTIONS)) {
243 fprintf(stderr, "%s", USAGE);
244 exit(EXIT_FAILURE);
245 }
246
247 minValue = maxValue = NULL;
248 output = input = NULL;
249 pipeFlags = 0;
250 abscissaName = NULL;
251 boundaryFile = boundaryColumn = boundaryColumnUnits = NULL;
252 boundaryValue = NULL;
253 offset = 0;
254 columnName = excludeName = NULL;
255 columnNames = excludeNames = abscissaNames = 0;
256 givenLowerLimits = givenUpperLimits = 0;
257 givenLowerLimit = givenUpperLimit = NULL;
258 bins = binsGiven = binSize = doSides = doSeparate = 0;
259 lowerLimitGiven = upperLimitGiven = 0;
260 expandRange = 0;
261 cdfOnly = 0;
262 freOnly = 1;
263 autoBinsTarget = 0;
264 autoBinsMinimum = autoBinsMaximum = 0;
265 weightColumn = NULL;
266 weightData = NULL;
267
268 for (iArg = 1; iArg < argc; iArg++) {
269 if (scanned[iArg].arg_type == OPTION) {
270 /* process options here */
271 switch (match_string(scanned[iArg].list[0], option, N_OPTIONS, 0)) {
272 case SET_MAJOR_ORDER:
273 majorOrderFlag = 0;
274 scanned[iArg].n_items--;
275 if (scanned[iArg].n_items > 0 && (!scanItemList(&majorOrderFlag, scanned[iArg].list + 1, &scanned[iArg].n_items, 0, "row", -1, NULL, 0, SDDS_ROW_MAJOR_ORDER, "column", -1, NULL, 0, SDDS_COLUMN_MAJOR_ORDER, NULL)))
276 SDDS_Bomb("invalid -majorOrder syntax/values");
277 if (majorOrderFlag & SDDS_COLUMN_MAJOR_ORDER)
278 columnMajorOrder = 1;
279 else if (majorOrderFlag & SDDS_ROW_MAJOR_ORDER)
280 columnMajorOrder = 0;
281 break;
282 case SET_PIPE:
283 if (!processPipeOption(scanned[iArg].list + 1, scanned[iArg].n_items - 1, &pipeFlags))
284 SDDS_Bomb("invalid -pipe syntax");
285 break;
286 case SET_COLUMNS:
287 if (columnName)
288 SDDS_Bomb("only one -columns option may be given");
289 if (scanned[iArg].n_items < 2)
290 SDDS_Bomb("invalid -columns syntax");
291 if (!(columnName = SDDS_Realloc(columnName, sizeof(*columnName) * (columnNames + scanned[iArg].n_items - 1))))
292 SDDS_Bomb("memory allocation failure");
293 for (i = 1; i < scanned[iArg].n_items; i++)
294 columnName[columnNames + i - 1] = scanned[iArg].list[i];
295 columnNames += scanned[iArg].n_items - 1;
296 break;
297 case SET_ABSCISSA:
298 if (abscissaName)
299 SDDS_Bomb("only one -abscissa option may be given");
300 if (scanned[iArg].n_items >= 2) {
301 if (!(abscissaName = SDDS_Realloc(abscissaName, sizeof(*abscissaName) * (abscissaNames + scanned[iArg].n_items - 1))))
302 SDDS_Bomb("memory allocation failure");
303 for (i = 1; i < scanned[iArg].n_items; i++)
304 abscissaName[abscissaNames + i - 1] = scanned[iArg].list[i];
305 abscissaNames += scanned[iArg].n_items - 1;
306 }
307 break;
308 case SET_BINS:
309 if (binsGiven)
310 SDDS_Bomb("-bins specified more than once");
311 binsGiven = 1;
312 if (sscanf(scanned[iArg].list[1], "%" SCNd64, &bins) != 1 || bins <= 0)
313 SDDS_Bomb("invalid value for bins---give a positive value");
314 break;
315 case SET_SIZEOFBINS:
316 if (sscanf(scanned[iArg].list[1], "%le", &binSize) != 1 || binSize <= 0)
317 SDDS_Bomb("invalid value for bin size---give a positive value");
318 break;
319 case SET_AUTOBINS:
320 if (scanned[iArg].n_items < 2)
321 SDDS_Bomb("incorrect -autoBins syntax");
322 scanned[iArg].n_items -= 1;
323 autoBinsTarget = autoBinsMinimum = autoBinsMaximum = 0;
324 if (!scanItemList(&dummyFlags, scanned[iArg].list + 1, &scanned[iArg].n_items, 0,
325 "target", SDDS_DOUBLE, &autoBinsTarget, 1, 0,
326 "minimum", SDDS_LONG, &autoBinsMinimum, 1, 0,
327 "maximum", SDDS_LONG, &autoBinsMaximum, 1, 0, NULL) ||
328 autoBinsTarget <= 0 || autoBinsMinimum < 0 || autoBinsMaximum < 0)
329 SDDS_Bomb("incorrect -autoBins syntax or values");
330 break;
331 case SET_EXCLUDE:
332 if (excludeName)
333 SDDS_Bomb("only one -exclude option may be given");
334 if (scanned[iArg].n_items < 2)
335 SDDS_Bomb("invalid -exclude syntax");
336 if (!(excludeName = SDDS_Realloc(excludeName, sizeof(*excludeName) * (excludeNames + scanned[iArg].n_items - 1))))
337 SDDS_Bomb("memory allocation failure");
338 for (i = 1; i < scanned[iArg].n_items; i++)
339 excludeName[excludeNames + i - 1] = scanned[iArg].list[i];
340 excludeNames += scanned[iArg].n_items - 1;
341 break;
342 case SET_LOWERLIMIT:
343 if (lowerLimitGiven)
344 SDDS_Bomb("-lowerLimit specified more than once");
345 lowerLimitGiven = 1;
346 if (!(givenLowerLimit = SDDS_Realloc(givenLowerLimit, sizeof(*givenLowerLimit) * (givenLowerLimits + scanned[iArg].n_items - 1))))
347 SDDS_Bomb("SET_LOWERLIMIT: memory allocation failure");
348 for (i = 1; i < scanned[iArg].n_items; i++) {
349 if (sscanf(scanned[iArg].list[i], "%lf", &transferLimit) != 1)
350 SDDS_Bomb("invalid value for -lowerLimit");
351 givenLowerLimit[givenLowerLimits + i - 1] = transferLimit;
352 }
353 givenLowerLimits += scanned[iArg].n_items - 1;
354 break;
355 case SET_UPPERLIMIT:
356 if (upperLimitGiven)
357 SDDS_Bomb("-upperLimit specified more than once");
358 upperLimitGiven = 1;
359 if (!(givenUpperLimit = SDDS_Realloc(givenUpperLimit, sizeof(*givenUpperLimit) * (givenUpperLimits + scanned[iArg].n_items - 1))))
360 SDDS_Bomb("SET_UPPERLIMIT: memory allocation failure");
361 for (i = 1; i < scanned[iArg].n_items; i++) {
362 if (sscanf(scanned[iArg].list[i], "%lf", &transferLimit) != 1)
363 SDDS_Bomb("invalid value for -upperLimit");
364 givenUpperLimit[givenUpperLimits + i - 1] = transferLimit;
365 }
366 givenUpperLimits += scanned[iArg].n_items - 1;
367 break;
368 case SET_SIDES:
369 doSides = DO_SIDES;
370 if (scanned[iArg].n_items == 2) {
371 static char *sideOpt[2] = {"close", "against"};
372 switch (match_string(scanned[iArg].list[1], sideOpt, 2, 0)) {
373 case 0:
374 doSides = DO_CLOSE_SIDES;
375 break;
376 case 1:
377 doSides = DO_AGAINST_SIDES;
378 break;
379 default:
380 SDDS_Bomb("invalid value for -sides");
381 break;
382 }
383 }
384 break;
385 case SET_SEPARATE:
386 doSeparate = 1;
387 break;
388 case SET_EXPAND:
389 if (scanned[iArg].n_items != 2 ||
390 sscanf(scanned[iArg].list[1], "%lf", &expandRange) != 1 || expandRange <= 0)
391 SDDS_Bomb("invalid -expand syntax");
392 break;
393 case SET_CDF:
394 if (scanned[iArg].n_items == 1)
395 cdfOnly = 0;
396 else {
397 if (scanned[iArg].n_items > 2)
398 SDDS_Bomb("invalid -cdf syntax");
399 CDFONLY = scanned[iArg].list[1];
400 if (strcmp(CDFONLY, "only") != 0)
401 SDDS_Bomb("invalid -cdf value, it should be -cdf or -cdf=only");
402 cdfOnly = 1;
403 }
404 freOnly = 0;
405 break;
406 case SET_BOUNDARYDATA:
407 if (scanned[iArg].n_items != 3 ||
408 !(boundaryFile = scanned[iArg].list[1]) ||
409 !strlen(boundaryFile) ||
410 !(boundaryColumn = scanned[iArg].list[2]) ||
411 !strlen(boundaryColumn))
412 SDDS_Bomb("invalid -boundaryData syntax or values");
413 break;
414 case SET_WEIGHT:
415 if (scanned[iArg].n_items != 2 ||
416 !(weightColumn = scanned[iArg].list[1]) ||
417 !strlen(weightColumn))
418 SDDS_Bomb("invalid -weightColumn syntax or values");
419 break;
420 case SET_NORMALIZE:
421 if (scanned[iArg].n_items == 1)
422 normMode = NORMALIZE_SUM;
423 else if (scanned[iArg].n_items != 2 ||
424 (normMode = match_string(scanned[iArg].list[1], normalize_option, N_NORMALIZE_OPTIONS, 0)) < 0)
425 SDDS_Bomb("invalid -normalize syntax");
426 break;
427 case SET_THREADS:
428 if (scanned[iArg].n_items != 2 ||
429 sscanf(scanned[iArg].list[1], "%d", &threads) != 1 || threads < 1)
430 SDDS_Bomb("invalid -threads syntax");
431 break;
432 default:
433 fprintf(stderr, "Error: unknown or ambiguous option: %s\n", scanned[iArg].list[0]);
434 fprintf(stderr, "%s", USAGE);
435 exit(EXIT_FAILURE);
436 break;
437 }
438 } else {
439 if (!input)
440 input = scanned[iArg].list[0];
441 else if (!output)
442 output = scanned[iArg].list[0];
443 else
444 SDDS_Bomb("too many filenames seen");
445 }
446 }
447
448 if (boundaryColumn && (abscissaNames > 0 || doSeparate))
449 SDDS_Bomb("-boundaryData option is incompatible with -abscissa and -separate options");
450
451 if (columnNames <= 0)
452 SDDS_Bomb("Supply the names of columns to histogram with -columns");
453
454 processFilenames("sddsmultihist", &input, &output, pipeFlags, 0, NULL);
455
456 if (!SDDS_InitializeInput(&SDDSin, input)) {
457 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
458 exit(EXIT_FAILURE);
459 }
460
461 if ((columnNames = expandColumnPairNames(&SDDSin, &columnName, NULL, columnNames, excludeName, excludeNames, FIND_NUMERIC_TYPE, 0)) <= 0) {
462 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
463 SDDS_Bomb("No quantities selected to histogram.");
464 }
465
466 if (doSeparate) {
467 if (abscissaNames <= 0) {
468 if (!(abscissaName = SDDS_Realloc(abscissaName, sizeof(*abscissaName) * (abscissaNames + columnNames))))
469 SDDS_Bomb("memory allocation failure");
470 abscissaNames = 0;
471 for (i = 0; i < columnNames; i++) {
472 abscissaName[abscissaNames + i] = columnName[i];
473 }
474 abscissaNames += columnNames;
475 }
476 if (columnNames > 1) {
477 if (abscissaNames > 0) {
478 if (columnNames != abscissaNames)
479 SDDS_Bomb("the number of abscissa names must match the number of columns");
480 }
481 if (givenLowerLimits)
482 if (columnNames != givenLowerLimits)
483 SDDS_Bomb("the number of lower limits must match the number of columns");
484 if (givenUpperLimits)
485 if (columnNames != givenUpperLimits)
486 SDDS_Bomb("the number of upper limits must match the number of columns");
487 } else {
488 /* Handle single column with multiple option names */
489 if (abscissaNames > 0) {
490 if (columnNames != abscissaNames)
491 abscissaNames = 1;
492 }
493 if (givenLowerLimits)
494 if (columnNames != givenLowerLimits)
495 givenLowerLimits = 1;
496 if (givenUpperLimits)
497 if (columnNames != givenUpperLimits)
498 givenUpperLimits = 1;
499 }
500 } else if (boundaryFile) {
501 if (!(boundaryValue = ReadBoundaryData(boundaryFile, boundaryColumn, &nBoundaryValues, &boundaryColumnUnits))) {
502 SDDS_Bomb("Problem reading boundary data");
503 }
504 } else {
505 if (abscissaNames <= 0)
506 SDDS_Bomb("Supply the name of the abscissa with -abscissaName");
507 abscissaNames = 1;
508 }
509
510 SetUpOutput(&SDDSout, &SDDSin, &abscissaIndex, &cdfIndex, &histogramIndex, output, columnName, columnNames, abscissaName, abscissaNames, boundaryColumn, boundaryColumnUnits, columnMajorOrder, normMode);
511
512 if (!(inputData = (double **)malloc(columnNames * sizeof(*inputData))) ||
513 !(minValue = (double *)malloc(columnNames * sizeof(*minValue))) ||
514 !(maxValue = (double *)malloc(columnNames * sizeof(*maxValue))))
515 SDDS_Bomb("memory allocation failure");
516
517 if (((binSize ? 1 : 0) + (binsGiven ? 1 : 0) + (autoBinsTarget ? 1 : 0)) > 1)
518 SDDS_Bomb("Specify only one of -binSize, -bins, or -autoBins");
519 if (!binSize && !binsGiven && !autoBinsTarget) {
520 binsGiven = 1;
521 bins = 20;
522 }
523
524 abscissa = histogram = cdf = NULL;
525
526 while ((readCode = SDDS_ReadPage(&SDDSin)) > 0) {
527 if ((rows = SDDS_CountRowsOfInterest(&SDDSin)) < 0)
528 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
529
530 if (rows) {
531 if (weightColumn && !(weightData = SDDS_GetColumnInDoubles(&SDDSin, weightColumn)))
532 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
533
534 for (column = 0; column < columnNames; column++) {
535 if (!(inputData[column] = SDDS_GetColumnInDoubles(&SDDSin, columnName[column])))
536 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
537 }
538
539 if (!boundaryColumn) {
540 if (!(lowerLimit = (double *)malloc(sizeof(*lowerLimit) * columnNames)))
541 SDDS_Bomb("memory allocation failure");
542 if (!(upperLimit = (double *)malloc(sizeof(*upperLimit) * columnNames)))
543 SDDS_Bomb("memory allocation failure");
544 if (!(dx = (double *)malloc(sizeof(*dx) * columnNames)))
545 SDDS_Bomb("memory allocation failure");
546
547 if (doSeparate) {
548 for (column = 0; column < columnNames; column++) {
549 find_min_max(&minValue[column], &maxValue[column], inputData[column], rows);
550 lowerLimit[column] = givenLowerLimits ? givenLowerLimit[column] : minValue[column];
551 upperLimit[column] = givenUpperLimits ? givenUpperLimit[column] : maxValue[column];
552 }
553 } else {
554 for (column = 0; column < columnNames; column++)
555 find_min_max(&minValue[column], &maxValue[column], inputData[column], rows);
556 lowerLimit[0] = givenLowerLimits ? givenLowerLimit[0] : min_in_array(minValue, columnNames);
557 upperLimit[0] = givenUpperLimits ? givenUpperLimit[0] : max_in_array(maxValue, columnNames);
558 for (column = 1; column < columnNames; column++) {
559 lowerLimit[column] = lowerLimit[0];
560 upperLimit[column] = upperLimit[0];
561 }
562 }
563
564 range = (1 + expandRange) * (upperLimit[0] - lowerLimit[0]);
565 if (autoBinsTarget) {
566 bins = (int64_t)(rows / autoBinsTarget);
567 if (autoBinsMinimum) {
568 if (bins < autoBinsMinimum)
569 bins = autoBinsMinimum;
570 } else if (bins < 5) {
571 bins = 5;
572 }
573 if (autoBinsMaximum) {
574 if (bins > autoBinsMaximum)
575 bins = autoBinsMaximum;
576 } else if (bins > rows)
577 bins = rows;
578 }
579
580 if (binSize) {
581 maxRange = range;
582 for (column = 0; column < columnNames; column++) {
583 range = (1 + expandRange) * (upperLimit[column] - lowerLimit[column]);
584 range = ((range / binSize) + 1) * binSize;
585 if (range > maxRange)
586 maxRange = range;
587 middle = (lowerLimit[column] + upperLimit[column]) / 2;
588 lowerLimit[column] = middle - range / 2;
589 upperLimit[column] = middle + range / 2;
590 }
591 range = maxRange;
592 if (!doSeparate) {
593 lowerLimit[0] = min_in_array(lowerLimit, columnNames);
594 upperLimit[0] = max_in_array(upperLimit, columnNames);
595 dx[0] = binSize;
596 }
597 bins = maxRange / binSize + 0.5;
598 if (bins < 1 && !doSides)
599 bins = 2;
600 if (doSeparate)
601 for (column = 0; column < columnNames; column++) {
602 range = upperLimit[column] - lowerLimit[column];
603 upperLimit[column] += (maxRange - range) / 2;
604 lowerLimit[column] -= (maxRange - range) / 2;
605 dx[column] = binSize;
606 }
607 } else {
608 if (doSeparate) {
609 for (column = 0; column < columnNames; column++) {
610 range = (1 + expandRange) * (upperLimit[column] - lowerLimit[column]);
611 middle = (upperLimit[column] + lowerLimit[column]) / 2;
612 upperLimit[column] = middle + range / 2;
613 lowerLimit[column] = middle - range / 2;
614 if (upperLimit[column] == lowerLimit[column]) {
615 if (fabs(upperLimit[column]) < sqrt(DBL_MIN)) {
616 upperLimit[column] = sqrt(DBL_MIN);
617 lowerLimit[column] = -sqrt(DBL_MIN);
618 } else {
619 lowerLimit[column] = upperLimit[column] * (1 - 10000 * DBL_EPSILON);
620 upperLimit[column] = upperLimit[column] * (1 + 10000 * DBL_EPSILON);
621 }
622 }
623 dx[column] = (upperLimit[column] - lowerLimit[column]) / bins;
624 }
625 } else {
626 range = (1 + expandRange) * (upperLimit[0] - lowerLimit[0]);
627 middle = (upperLimit[0] + lowerLimit[0]) / 2;
628 upperLimit[0] = middle + range / 2;
629 lowerLimit[0] = middle - range / 2;
630 if (upperLimit[0] == lowerLimit[0]) {
631 if (fabs(upperLimit[0]) < sqrt(DBL_MIN)) {
632 upperLimit[0] = sqrt(DBL_MIN);
633 lowerLimit[0] = -sqrt(DBL_MIN);
634 } else {
635 lowerLimit[0] = upperLimit[0] * (1 - 10000 * DBL_EPSILON);
636 upperLimit[0] = upperLimit[0] * (1 + 10000 * DBL_EPSILON);
637 }
638 }
639 dx[0] = (upperLimit[0] - lowerLimit[0]) / bins;
640 }
641 }
642 if (!binsGiven || !abscissa) {
643 if (!(abscissa = SDDS_Realloc(abscissa, sizeof(*abscissa) * (bins + 2))) ||
644 !(cdf = SDDS_Realloc(cdf, sizeof(*cdf) * (bins + 2))) ||
645 !(histogram = SDDS_Realloc(histogram, sizeof(*histogram) * (bins + 2))))
646 SDDS_Bomb("memory allocation failure");
647 }
648 writeBins = bins + (doSides ? 2 : 0);
649 offset = doSides ? 0 : 1;
650 } else {
651 bins = writeBins = nBoundaryValues;
652 abscissa = NULL;
653 if (!(cdf = SDDS_Realloc(cdf, sizeof(*cdf) * bins)) ||
654 !(histogram = SDDS_Realloc(histogram, sizeof(*histogram) * bins)))
655 SDDS_Bomb("memory allocation failure");
656 }
657
658 if (!SDDS_StartPage(&SDDSout, writeBins) || !SDDS_CopyParameters(&SDDSout, &SDDSin))
659 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
660
661 if (boundaryColumn) {
662 for (column = 0; column < columnNames; column++) {
663 MakeBoundaryHistogram(histogram, cdf, boundaryValue, nBoundaryValues, inputData[column], weightData, rows);
664 NormalizeHistogram(histogram, nBoundaryValues, normMode);
665 if (!cdfOnly && !SDDS_SetColumn(&SDDSout, SDDS_SET_BY_INDEX, histogram, writeBins, histogramIndex[column]))
666 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
667 if (!freOnly && !SDDS_SetColumn(&SDDSout, SDDS_SET_BY_INDEX, cdf, writeBins, cdfIndex[column]))
668 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
669 free(inputData[column]);
670 }
671 if (!SDDS_SetColumn(&SDDSout, SDDS_SET_BY_NAME, boundaryValue, writeBins, boundaryColumn))
672 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
673 } else if (!doSeparate) {
674 for (i = -1; i < bins + 1; i++)
675 abscissa[i + 1] = (i + 0.5) * dx[0] + lowerLimit[0];
676 switch (doSides) {
677 case DO_CLOSE_SIDES:
678 abscissa[0] = abscissa[1] - dx[0] / 2;
679 abscissa[bins + 1] = abscissa[bins] + dx[0] / 2;
680 break;
681 case DO_AGAINST_SIDES:
682 abscissa[0] = abscissa[1];
683 abscissa[bins + 1] = abscissa[bins];
684 break;
685 }
686 if (!SDDS_SetColumn(&SDDSout, SDDS_SET_BY_INDEX, abscissa + offset, writeBins, abscissaIndex[0]))
687 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
688 for (column = 0; column < columnNames; column++) {
689 histogram[0] = histogram[bins + 1] = 0;
690 if (!weightColumn)
691 make_histogram_threaded(histogram + 1, bins, lowerLimit[0], upperLimit[0], inputData[column], rows, 1, threads);
692 else
693 make_histogram_weighted_threaded(histogram + 1, bins, lowerLimit[0], upperLimit[0], inputData[column], rows, 1, weightData, threads);
694 NormalizeHistogram(histogram, bins, normMode);
695 sum = 0;
696 for (i = 0; i <= bins + 1; i++)
697 sum += histogram[i];
698 for (i = 0; i <= bins + 1; i++) {
699 if (!i)
700 cdf[i] = histogram[i] / sum;
701 else
702 cdf[i] = cdf[i - 1] + histogram[i] / sum;
703 }
704 if (!cdfOnly) {
705 if (!SDDS_SetColumn(&SDDSout, SDDS_SET_BY_INDEX, histogram + offset, writeBins, histogramIndex[column]))
706 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
707 }
708 if (!freOnly) {
709 if (!SDDS_SetColumn(&SDDSout, SDDS_SET_BY_INDEX, cdf + offset, writeBins, cdfIndex[column]))
710 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
711 }
712 free(inputData[column]);
713 }
714 } else {
715 for (column = 0; column < columnNames; column++) {
716 for (i = -1; i < bins + 1; i++)
717 abscissa[i + 1] = (i + 0.5) * dx[column] + lowerLimit[column];
718 switch (doSides) {
719 case DO_CLOSE_SIDES:
720 abscissa[0] = abscissa[1] - dx[column] / 2;
721 abscissa[bins + 1] = abscissa[bins] + dx[column] / 2;
722 break;
723 case DO_AGAINST_SIDES:
724 abscissa[0] = abscissa[1];
725 abscissa[bins + 1] = abscissa[bins];
726 break;
727 }
728 if (!SDDS_SetColumn(&SDDSout, SDDS_SET_BY_INDEX, abscissa + offset, writeBins, abscissaIndex[column]))
729 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
730 histogram[0] = histogram[bins + 1] = 0;
731 if (!weightColumn)
732 make_histogram_threaded(histogram + 1, bins, lowerLimit[column], upperLimit[column], inputData[column], rows, 1, threads);
733 else
734 make_histogram_weighted_threaded(histogram + 1, bins, lowerLimit[column], upperLimit[column], inputData[column], rows, 1, weightData, threads);
735 NormalizeHistogram(histogram, bins + 2, normMode);
736 sum = 0;
737 for (i = 0; i <= bins + 1; i++)
738 sum += histogram[i];
739 for (i = 0; i <= bins + 1; i++) {
740 if (!i)
741 cdf[i] = histogram[i] / sum;
742 else
743 cdf[i] = cdf[i - 1] + histogram[i] / sum;
744 }
745 if (!cdfOnly) {
746 if (!SDDS_SetColumn(&SDDSout, SDDS_SET_BY_INDEX, histogram + offset, writeBins, histogramIndex[column]))
747 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
748 }
749 if (!freOnly) {
750 if (!SDDS_SetColumn(&SDDSout, SDDS_SET_BY_INDEX, cdf + offset, writeBins, cdfIndex[column]))
751 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
752 }
753 free(inputData[column]);
754 }
755 }
756 } else {
757 if (!SDDS_StartPage(&SDDSout, rows ? bins : 0))
758 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
759 }
760 if (weightData) {
761 free(weightData);
762 weightData = NULL;
763 }
764 if (!SDDS_WritePage(&SDDSout))
765 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
766 }
767
768 if (!SDDS_Terminate(&SDDSin)) {
769 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
770 exit(EXIT_FAILURE);
771 }
772 if (!SDDS_Terminate(&SDDSout)) {
773 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
774 exit(EXIT_FAILURE);
775 }
776
777 free(histogram);
778 free(cdf);
779 free(abscissa);
780 if (minValue)
781 free(minValue);
782 if (maxValue)
783 free(maxValue);
784 if (lowerLimit)
785 free(lowerLimit);
786 if (upperLimit)
787 free(upperLimit);
788 if (dx)
789 free(dx);
790
791 return EXIT_SUCCESS;
792}
793
794void SetUpOutput(SDDS_DATASET *SDDSout, SDDS_DATASET *SDDSin, long **abscissaIndex, long **cdfIndex,
795 long **histogramIndex, char *output, char **columnName, long columnNames,
796 char **abscissaName, long abscissaNames, char *boundaryColumn, char *boundaryUnits,
797 short columnMajorOrder, short normMode) {
798 long column;
799 char s[SDDS_MAXLINE];
800 int32_t outputType = SDDS_DOUBLE;
801 char *blankString = "";
802
803 if (!SDDS_InitializeOutput(SDDSout, SDDS_BINARY, 0, NULL, "sddsmultihist output", output))
804 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
805
806 if (columnMajorOrder != -1)
807 SDDSout->layout.data_mode.column_major = columnMajorOrder;
808 else
809 SDDSout->layout.data_mode.column_major = SDDSin->layout.data_mode.column_major;
810
811 if (!(*cdfIndex = (long *)malloc(sizeof(**cdfIndex) * columnNames)))
812 SDDS_Bomb("memory allocation failure");
813 if (!cdfOnly) {
814 if (!(*histogramIndex = (long *)malloc(sizeof(**histogramIndex) * columnNames)))
815 SDDS_Bomb("memory allocation failure");
816 }
817 if (!boundaryColumn) {
818 if (!(*abscissaIndex = (long *)malloc(sizeof(**abscissaIndex) * columnNames)))
819 SDDS_Bomb("memory allocation failure");
820
821 for (column = 0; column < abscissaNames; column++) {
822 if (!SDDS_TransferColumnDefinition(SDDSout, SDDSin, columnName[column], abscissaName[column]) ||
823 !SDDS_ChangeColumnInformation(SDDSout, "type", &outputType, SDDS_BY_NAME, abscissaName[column]) ||
824 ((*abscissaIndex)[column] = SDDS_GetColumnIndex(SDDSout, abscissaName[column])) < 0)
825 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
826 if (strcmp(columnName[column], abscissaName[column]) != 0 &&
827 (!SDDS_ChangeColumnInformation(SDDSout, "description", blankString, SDDS_BY_NAME, abscissaName[column]) ||
828 !SDDS_ChangeColumnInformation(SDDSout, "symbol", blankString, SDDS_BY_NAME, abscissaName[column])))
829 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
830 }
831 } else {
832 if (!SDDS_DefineSimpleColumn(SDDSout, boundaryColumn, boundaryUnits, SDDS_DOUBLE))
833 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
834 }
835
836 for (column = 0; column < columnNames; column++) {
837 if (!freOnly) {
838 sprintf(s, "%sCdf", columnName[column]);
839 if (((*cdfIndex)[column] = SDDS_DefineColumn(SDDSout, s, NULL, NULL, NULL, NULL, SDDS_DOUBLE, 0)) < 0)
840 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
841 }
842 if (!cdfOnly) {
843 switch (normMode) {
844 case NORMALIZE_PEAK:
845 sprintf(s, "%sRelativeFrequency", columnName[column]);
846 break;
847 case NORMALIZE_SUM:
848 sprintf(s, "%sFractionalFrequency", columnName[column]);
849 break;
850 default:
851 sprintf(s, "%sFrequency", columnName[column]);
852 break;
853 }
854 if (((*histogramIndex)[column] = SDDS_DefineColumn(SDDSout, s, NULL, NULL, NULL, NULL, SDDS_DOUBLE, 0)) < 0)
855 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
856 }
857 }
858 if (!SDDS_TransferAllParameterDefinitions(SDDSout, SDDSin, SDDS_TRANSFER_KEEPOLD) ||
859 !SDDS_WriteLayout(SDDSout))
860 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
861}
862
863double *ReadBoundaryData(char *file, char *column, int64_t *n, char **units) {
864 SDDS_DATASET SDDSin;
865 double *data;
866 int64_t j;
867
868 if (!SDDS_InitializeInput(&SDDSin, file)) {
869 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
870 }
871
872 if (SDDS_GetColumnInformation(&SDDSin, "units", units, SDDS_GET_BY_NAME, column) != SDDS_STRING)
873 return NULL;
874
875 if (SDDS_ReadPage(&SDDSin) <= 0) {
876 SDDS_SetError("No pages in boundary data file");
877 return NULL;
878 }
879
880 *n = SDDS_RowCount(&SDDSin);
881 if (*n <= 0)
882 return NULL;
883 if (!(data = SDDS_GetColumnInDoubles(&SDDSin, column)))
884 return NULL;
885 for (j = 1; j < (*n); j++) {
886 if (data[j] <= data[j - 1]) {
887 memmove(data + j, data + j + 1, sizeof(*data) * ((*n) - 1 - j));
888 *n -= 1;
889 j -= 1;
890 }
891 }
892 return data;
893}
894
895void MakeBoundaryHistogram(double *histogram, double *cdf, double *boundaryValue, int64_t nBoundaryValues,
896 double *data, double *weight, int64_t nData) {
897 int64_t i, j;
898 for (i = 0; i < nBoundaryValues; i++)
899 histogram[i] = cdf[i] = 0;
900 for (i = 0; i < nData; i++) {
901 j = binaryArraySearch(boundaryValue, sizeof(double), nBoundaryValues, &data[i], double_cmpasc, 1) + 1;
902 if (j < nBoundaryValues && j >= 0)
903 histogram[j] += weight ? fabs(weight[i]) : 1;
904 }
905 cdf[0] = histogram[0];
906 for (j = 1; j < nBoundaryValues; j++) {
907 cdf[j] = cdf[j - 1] + histogram[j];
908 }
909 if (cdf[nBoundaryValues - 1] > 0)
910 for (j = 0; j < nBoundaryValues; j++)
911 cdf[j] /= cdf[nBoundaryValues - 1];
912}
913
914void NormalizeHistogram(double *hist, int64_t bins, short mode) {
915 int64_t i;
916 double value = 0;
917 switch (mode) {
918 case NORMALIZE_SUM:
919 for (i = 0; i < bins; i++)
920 value += hist[i];
921 break;
922 case NORMALIZE_PEAK:
923 for (i = 0; i < bins; i++)
924 if (hist[i] > value)
925 value = hist[i];
926 break;
927 default:
928 return;
929 break;
930 }
931 if (value)
932 for (i = 0; i < bins; i++)
933 hist[i] /= value;
934}
935
936static long make_histogram_threaded(double *hist, long n_bins, double lo, double hi, double *data,
937 int64_t n_pts, long new_start, int threads) {
938 double bin_size, *partial;
939 long *counts, count = 0;
940 int activeThreads, thread;
941 int64_t i;
942
943 if (threads <= 1 || n_pts <= 0 || n_bins <= 0)
944 return make_histogram(hist, n_bins, lo, hi, data, n_pts, new_start);
945
946 activeThreads = threads;
947 if (activeThreads > n_pts)
948 activeThreads = (int)n_pts;
949 if (activeThreads <= 1)
950 return make_histogram(hist, n_bins, lo, hi, data, n_pts, new_start);
951
952 if (new_start)
953 for (i = 0; i < n_bins; i++)
954 hist[i] = 0;
955 bin_size = (hi - lo) / n_bins;
956 partial = calloc((size_t)activeThreads * n_bins, sizeof(*partial));
957 counts = calloc(activeThreads, sizeof(*counts));
958 if (!partial || !counts)
959 SDDS_Bomb("memory allocation failure");
960
961#pragma omp parallel for if (activeThreads > 1) num_threads(activeThreads)
962 for (thread = 0; thread < activeThreads; thread++) {
963 int64_t start = thread * (n_pts / activeThreads);
964 int64_t end = (thread == activeThreads - 1) ? n_pts : (thread + 1) * (n_pts / activeThreads);
965 double *local = partial + (size_t)thread * n_bins;
966 long localCount = 0;
967 for (int64_t point = start; point < end; point++) {
968 double dbin = (data[point] - lo) / bin_size;
969 long bin = dbin;
970 if (dbin < 0)
971 continue;
972 if (bin < 0 || bin >= n_bins)
973 continue;
974 local[bin] += 1;
975 localCount++;
976 }
977 counts[thread] = localCount;
978 }
979
980 for (thread = 0; thread < activeThreads; thread++) {
981 double *local = partial + (size_t)thread * n_bins;
982 count += counts[thread];
983 for (i = 0; i < n_bins; i++)
984 hist[i] += local[i];
985 }
986 free(partial);
987 free(counts);
988 return count;
989}
990
991static long make_histogram_weighted_threaded(double *hist, long n_bins, double lo, double hi, double *data,
992 int64_t n_pts, long new_start, double *weight, int threads) {
993 double bin_size, *partial;
994 long *counts, count = 0;
995 int activeThreads, thread;
996 int64_t i;
997
998 if (threads <= 1 || n_pts <= 0 || n_bins <= 0)
999 return make_histogram_weighted(hist, n_bins, lo, hi, data, n_pts, new_start, weight);
1000
1001 activeThreads = threads;
1002 if (activeThreads > n_pts)
1003 activeThreads = (int)n_pts;
1004 if (activeThreads <= 1)
1005 return make_histogram_weighted(hist, n_bins, lo, hi, data, n_pts, new_start, weight);
1006
1007 if (new_start)
1008 for (i = 0; i < n_bins; i++)
1009 hist[i] = 0;
1010 bin_size = (hi - lo) / n_bins;
1011 partial = calloc((size_t)activeThreads * n_bins, sizeof(*partial));
1012 counts = calloc(activeThreads, sizeof(*counts));
1013 if (!partial || !counts)
1014 SDDS_Bomb("memory allocation failure");
1015
1016#pragma omp parallel for if (activeThreads > 1) num_threads(activeThreads)
1017 for (thread = 0; thread < activeThreads; thread++) {
1018 int64_t start = thread * (n_pts / activeThreads);
1019 int64_t end = (thread == activeThreads - 1) ? n_pts : (thread + 1) * (n_pts / activeThreads);
1020 double *local = partial + (size_t)thread * n_bins;
1021 long localCount = 0;
1022 for (int64_t point = start; point < end; point++) {
1023 double dbin = (data[point] - lo) / bin_size;
1024 long bin = dbin;
1025 if (dbin < 0)
1026 continue;
1027 if (bin < 0 || bin >= n_bins)
1028 continue;
1029 local[bin] += weight[point];
1030 localCount++;
1031 }
1032 counts[thread] = localCount;
1033 }
1034
1035 for (thread = 0; thread < activeThreads; thread++) {
1036 double *local = partial + (size_t)thread * n_bins;
1037 count += counts[thread];
1038 for (i = 0; i < n_bins; i++)
1039 hist[i] += local[i];
1040 }
1041 free(partial);
1042 free(counts);
1043 return count;
1044}
SDDS.h
SDDS (Self Describing Data Set) Data Types Definitions and Function Prototypes.
SDDS_CopyParameters
int32_t SDDS_CopyParameters(SDDS_DATASET *SDDS_target, SDDS_DATASET *SDDS_source)
Definition SDDS_copy.c:286
SDDS_StartPage
int32_t SDDS_StartPage(SDDS_DATASET *SDDS_dataset, int64_t expected_n_rows)
Definition SDDS_dataprep.c:70
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_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_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_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_DefineSimpleColumn
int32_t SDDS_DefineSimpleColumn(SDDS_DATASET *SDDS_dataset, const char *name, const char *unit, int32_t type)
Defines a simple data column within the SDDS dataset.
Definition SDDS_output.c:1853
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_DefineColumn
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.
Definition SDDS_output.c:1716
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_TransferAllParameterDefinitions
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.
Definition SDDS_transfer.c:343
SDDS_SetError
void SDDS_SetError(char *error_text)
Records an error message in the SDDS error stack.
Definition SDDS_utils.c:421
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_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_Bomb
void SDDS_Bomb(char *message)
Terminates the program after printing an error message and recorded errors.
Definition SDDS_utils.c:380
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
SDDS_DOUBLE
#define SDDS_DOUBLE
Identifier for the double data type.
Definition SDDStypes.h:37
SDDSutils.h
Utility functions for SDDS dataset manipulation and string array operations.
binaryArraySearch
long binaryArraySearch(void *array, size_t elemSize, long members, void *key, int(*compare)(const void *c1, const void *c2), long bracket)
Searches for a key in a sorted array of data values using binary search.
Definition binsert.c:156
find_min_max
int find_min_max(double *min, double *max, double *list, int64_t n)
Finds the minimum and maximum values in a list of doubles.
Definition findMinMax.c:33
max_in_array
double max_in_array(double *array, long n)
Finds the maximum value in an array of doubles.
Definition findMinMax.c:318
min_in_array
double min_in_array(double *array, long n)
Finds the minimum value in an array of doubles.
Definition findMinMax.c:336
make_histogram_weighted
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.
Definition makeHistogram.c:74
make_histogram
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.
Definition makeHistogram.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
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
double_cmpasc
int double_cmpasc(const void *a, const void *b)
Compare two doubles in ascending order.
Definition sortfunctions.c:28
SDDS_DATASET
Definition SDDS.h:308
Generated by doxygen 1.12.0