SDDS ToolKit Programs and Libraries for C and Python
Loading...
Searching...
No Matches
Functions
sddshist.c File Reference

Detailed Description

Command-line tool for generating histograms from SDDS-formatted data.

This program is a command-line tool designed to generate histograms from data contained in SDDS (Self Describing Data Set) files. It supports various configurations like normalization, filtering, multi-threading for performance, and the ability to define histograms based on specific regions in the data.

Usage

sddshist [<inputfile>] [<outputfile>]
[-pipe=[input][,output]]
-dataColumn=<column-name>
[{
-bins=<number> |
-sizeOfBins=<value> |
-regions=filename=<filename>,position=<columnName>,name=<columnName>
}]
[-lowerLimit=<value>]
[-upperLimit=<value>]
[-expand=<factor>]
[-filter=<column-name>,<lower-limit>,<upper-limit>]
[-weightColumn=<column-name>]
[-sides[=<points>]]
[-normalize[={sum|area|peak}]]
[-cdf[=only]]
[-threads=<number>]
[-statistics]
[-verbose]
[-majorOrder=row|column]

Options

Required Description
-dataColumn=<column-name> Specifies the column to be histogrammed.
Optional Description
-pipe Use pipe for input and/or output.
-bins Set the number of bins for the histogram.
-sizeOfBins Set the size of each bin.
-regions Define region-based histogramming with file and column information.
-lowerLimit Set the lower limit for the histogram.
-upperLimit Set the upper limit for the histogram.
-expand Expand the histogram range by a given factor.
-filter Filter data points based on values in a specified column.
-weightColumn Specify a column to weight the histogram.
-sides Extend the histogram to zero level.
-normalize Normalize the histogram (sum, area, or peak).
-cdf Include the Cumulative Distribution Function (CDF) in the output.
-threads Set the number of threads for processing.
-statistics Include statistical details in the output.
-verbose Print additional processing details.
-majorOrder Set data major order (row or column).

Incompatibilities

  • -bins, -sizeOfBins, and -regions are mutually exclusive.
  • -normalize and -cdf=only cannot be used together.
Copyright
  • (c) 2002 The University of Chicago, as Operator of Argonne National Laboratory.
  • (c) 2002 The Regents of the University of California, as Operator of Los Alamos National Laboratory.
License
This file is distributed under the terms of the Software License Agreement found in the file LICENSE included with this distribution.
Authors
M. Borland, C. Saunders, R. Soliday, H. Shang

Definition in file sddshist.c.

#include "mdb.h"
#include "scan.h"
#include "SDDS.h"

Go to the source code of this file.

Functions

static int64_t filter (double *x, double *y, double *filterData, int64_t npts, double lower_filter, double upper_filter)
 
static long setupOutputFile (SDDS_DATASET *outTable, char *outputfile, SDDS_DATASET *inTable, char *inputfile, char *dataColumn, char *weightColumn, char *filterColumn, double lowerFilter, double upperFilter, SDDS_DATASET *regionTable, char *regionNameColumn, long doStats, int64_t bins, double binSize, long normalizeMode, short columnMajorOrder)
 
int64_t readRegionFile (SDDS_DATASET *SDDSin, char *filename, char *positionColumn, char *nameColumn, double **regionPosition, char ***regionName)
 
void classifyByRegion (double *data, double *weight, int64_t points, double *histogram, double *regionPosition, int64_t bins)
 
static long make_histogram_threaded (double *hist, long n_bins, double lo, double hi, double *data, int64_t n_pts, long new_start, int threads)
 
static long make_histogram_weighted_threaded (double *hist, long n_bins, double lo, double hi, double *data, int64_t n_pts, long new_start, double *weight, int threads)
 
int main (int argc, char **argv)
 

Function Documentation

◆ classifyByRegion()

void classifyByRegion ( double * data,
double * weight,
int64_t points,
double * histogram,
double * regionPosition,
int64_t bins )

Definition at line 756 of file sddshist.c.

756 {
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}

◆ filter()

static int64_t filter ( double * x,
double * y,
double * filterData,
int64_t npts,
double lower_filter,
double upper_filter )
static

Definition at line 590 of file sddshist.c.

590 {
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}
trealloc
void * trealloc(void *old_ptr, uint64_t size_of_block)
Reallocates a memory block to a new size.
Definition array.c:190

◆ main()

int main ( int argc,
char ** argv )

Definition at line 175 of file sddshist.c.

175 {
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
213 SDDS_RegisterProgramName(argv[0]);
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}
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_SetParameters
int32_t SDDS_SetParameters(SDDS_DATASET *SDDS_dataset, int32_t mode,...)
Definition SDDS_dataprep.c:375
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_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_WritePage
int32_t SDDS_WritePage(SDDS_DATASET *SDDS_dataset)
Writes the current data table to the output file.
Definition SDDS_output.c:1229
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_STRING
#define SDDS_STRING
Identifier for the string data type.
Definition SDDStypes.h:85
tmalloc
void * tmalloc(uint64_t size_of_block)
Allocates a memory block of the specified size with zero initialization.
Definition array.c:65
max_in_array
double max_in_array(double *array, long n)
Finds the maximum value in an array of doubles.
Definition findMinMax.c:318
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
computeWeightedMomentsThreaded
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
computeMomentsThreaded
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
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

◆ make_histogram_threaded()

static long make_histogram_threaded ( double * hist,
long n_bins,
double lo,
double hi,
double * data,
int64_t n_pts,
long new_start,
int threads )
static

Definition at line 776 of file sddshist.c.

777 {
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}
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

◆ make_histogram_weighted_threaded()

static long make_histogram_weighted_threaded ( double * hist,
long n_bins,
double lo,
double hi,
double * data,
int64_t n_pts,
long new_start,
double * weight,
int threads )
static

Definition at line 831 of file sddshist.c.

832 {
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}
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

◆ readRegionFile()

int64_t readRegionFile ( SDDS_DATASET * SDDSin,
char * filename,
char * positionColumn,
char * nameColumn,
double ** regionPosition,
char *** regionName )

Definition at line 736 of file sddshist.c.

736 {
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}
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_Realloc
void * SDDS_Realloc(void *old_ptr, size_t new_size)
Reallocates memory to a new size.
Definition SDDS_utils.c:743
cp_str
char * cp_str(char **s, char *t)
Copies a string, allocating memory for storage.
Definition cp_str.c:28

◆ setupOutputFile()

static long setupOutputFile ( SDDS_DATASET * outTable,
char * outputfile,
SDDS_DATASET * inTable,
char * inputfile,
char * dataColumn,
char * weightColumn,
char * filterColumn,
double lowerFilter,
double upperFilter,
SDDS_DATASET * regionTable,
char * regionNameColumn,
long doStats,
int64_t bins,
double binSize,
long normalizeMode,
short columnMajorOrder )
static

Definition at line 620 of file sddshist.c.

623 {
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}
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_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_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_DefineParameter
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.
Definition SDDS_output.c:1473
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_StringIsBlank
int32_t SDDS_StringIsBlank(char *s)
Checks if a string is blank (contains only whitespace characters).
Definition SDDS_utils.c:2470
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
Generated by doxygen 1.12.0