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

Detailed Description

Generates sampled distributions based on input SDDS files or direct specifications.

This program generates samples from specified distributions (Gaussian, Uniform, Poisson) or based on cumulative distribution functions (CDF/DF) provided in input SDDS files. It supports various options to customize the sampling process, including the number of samples, random seed, verbosity, and output ordering.

Usage

sddssampledist [<inputfile>] [<outputfile>]
[-pipe=[in][,out]]
-columns=independentVariable=<name>,{cdf=<CDFName> | df=<DFName>}[,output=<name>][,units=<string>][,factor=<value>][,offset=<value>][,datafile=<filename>][,haltonRadix=<primeNumber>[,haltonOffset=<integer>][,randomize[,group=<groupID>]]]
[-columns=...]
-samples=<integer>
[-seed=<integer>]
[-verbose]
[-gaussian=columnName=<columnName>[,meanValue=<value>|@<parameter_name>][,sigmaValue=<value>|@<parameter_name>][,units=<string>]]
[-uniform=columnName=<columnName>[,minimumValue=<value>|@<parameter_name>][,maximumValue=<value>|@<parameter_name>][,units=<string>]]
[-poisson=columnName=<columnName>[,meanValue=<value>|@<parameter_name>][,units=<string>]]
[-optimalHalton]
[-majorOrder=row|column]
[-threads=<number>]

Options

Required Description
-columns Defines independent variable and distribution function with customization options.
-samples Specifies the number of samples to generate.
Optional Description
-pipe Use standard input and/or output streams.
-seed Specifies the seed for the random number generator.
-verbose Enables verbose output.
-gaussian Samples from a Gaussian distribution.
-uniform Samples from a Uniform distribution.
-poisson Samples from a Poisson distribution.
-majorOrder Specifies the output file order as row-major or column-major.
-threads Number of datafile-backed distributions to read concurrently.

Incompatibilities

  • -columns
    • Requires at least one independentVariable and exactly one of cdf or df.
    • If haltonRadix is used, the value must be a prime number.
    • Randomization requires valid group identifiers if group=<groupID> is specified.
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, R. Soliday, H. Shang

Definition in file sddssampledist.c.

#include "mdb.h"
#include "scan.h"
#include "SDDS.h"
#include <time.h>

Go to the source code of this file.

Functions

long CreatePoissonDistributionTable (double **x, double **pos_CDF, double mean)
 
int main (int argc, char **argv)
 

Function Documentation

◆ CreatePoissonDistributionTable()

long CreatePoissonDistributionTable ( double ** x,
double ** pos_CDF,
double mean )

Definition at line 818 of file sddssampledist.c.

818 {
819 long i, npoints = 20, count = 0;
820 double *pos = NULL;
821 /* SDDS_DATASET pos_out; */
822
823 *x = *pos_CDF = NULL;
824 if (!(*x = malloc(sizeof(**x) * npoints)) ||
825 !(pos = malloc(sizeof(*pos) * npoints)) ||
826 !(*pos_CDF = malloc(sizeof(**pos_CDF) * npoints)))
827 SDDS_Bomb("memory allocation failure.");
828 i = count = 0;
829 while (1) {
830 if (count + 2 >= npoints) {
831 npoints += 20;
832 *x = SDDS_Realloc(*x, sizeof(**x) * npoints);
833 *pos_CDF = SDDS_Realloc(*pos_CDF, sizeof(**pos_CDF) * npoints);
834 pos = SDDS_Realloc(pos, sizeof(*pos) * npoints);
835 }
836
837 (*x)[count] = i;
838 if (!i) {
839 pos[i] = exp(-mean);
840 (*pos_CDF)[count] = pos[i];
841 count++;
842 } else {
843 pos[i] = pos[i - 1] * mean / i;
844 (*pos_CDF)[count] = (*pos_CDF)[count - 1];
845 (*pos_CDF)[count + 1] = (*pos_CDF)[count - 1] + pos[i];
846 (*x)[count + 1] = i;
847 if (1.0 - (*pos_CDF)[count + 1] <= 1.0e-15)
848 break;
849 count += 2;
850 }
851 i++;
852 }
853 /* fprintf(stderr,"lamda=%f\n", mean);
854 if (!SDDS_InitializeOutput(&pos_out, SDDS_BINARY, 0, NULL, NULL, "pos_dist.sdds"))
855 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors|SDDS_EXIT_PrintErrors);
856 if (!SDDS_DefineSimpleColumn(&pos_out, "Count", NULL, SDDS_DOUBLE) ||
857 !SDDS_DefineSimpleColumn(&pos_out, "P", NULL, SDDS_DOUBLE))
858 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors|SDDS_EXIT_PrintErrors);
859 if (!SDDS_SaveLayout(&pos_out) || !SDDS_WriteLayout(&pos_out))
860 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors|SDDS_EXIT_PrintErrors);
861 if (!SDDS_StartPage(&pos_out, count) ||
862 !SDDS_SetColumnFromDoubles(&pos_out, SDDS_SET_BY_NAME, *x, count, "Count") ||
863 !SDDS_SetColumnFromDoubles(&pos_out, SDDS_SET_BY_NAME, *pos_CDF, count, "P"))
864 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors|SDDS_EXIT_PrintErrors);
865 if (!SDDS_WritePage(&pos_out) || !SDDS_Terminate(&pos_out))
866 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors|SDDS_EXIT_PrintErrors);
867 */
868 free(pos);
869 return count;
870}
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

◆ main()

int main ( int argc,
char ** argv )

Definition at line 191 of file sddssampledist.c.

191 {
192 int iArg;
193 char *input, *output, *meanPar, *sigmaPar, *maxPar, *minPar;
194 long i, mainInputOpened, haltonID = 0, requireInput = 0;
195 unsigned long pipeFlags, majorOrderFlag;
196 SCANNED_ARG *scanned;
197 SDDS_DATASET SDDSin, SDDSout, *SDDSptr;
198 long randomNumberSeed = 0;
199 SEQ_REQUEST *seqRequest;
200 long samples, seqRequests, randomizationGroups = 0;
201 int64_t j, values;
202 double *sample, *IVValue, *CDFValue;
203 char msgBuffer[1000];
204 RANDOMIZED_ORDER *randomizationData = NULL;
205 long verbose, optimalHalton = 0;
206 short columnMajorOrder = -1;
207 int threads = 1;
208 int *dataFilePageStatus = NULL;
209
211 argc = scanargs(&scanned, argc, argv);
212 if (argc < 2) {
213 fprintf(stderr, "%s%s%s\n", USAGE1, USAGE2, USAGE3);
214 return EXIT_FAILURE;
215 }
216 seqRequest = NULL;
217 seqRequests = 0;
218 output = input = NULL;
219 pipeFlags = 0;
220 samples = values = 0;
221 sample = IVValue = CDFValue = NULL;
222 verbose = 0;
223 maxPar = minPar = meanPar = sigmaPar = NULL;
224
225 for (iArg = 1; iArg < argc; iArg++) {
226 if (scanned[iArg].arg_type == OPTION) {
227 /* process options here */
228 switch (match_string(scanned[iArg].list[0], option, CLO_OPTIONS, 0)) {
229 case CLO_MAJOR_ORDER:
230 majorOrderFlag = 0;
231 scanned[iArg].n_items--;
232 if (scanned[iArg].n_items > 0 &&
233 (!scanItemList(&majorOrderFlag, scanned[iArg].list + 1, &scanned[iArg].n_items, 0,
234 "row", -1, NULL, 0, SDDS_ROW_MAJOR_ORDER,
235 "column", -1, NULL, 0, SDDS_COLUMN_MAJOR_ORDER, NULL)))
236 SDDS_Bomb("invalid -majorOrder syntax/values");
237 if (majorOrderFlag & SDDS_COLUMN_MAJOR_ORDER)
238 columnMajorOrder = 1;
239 else if (majorOrderFlag & SDDS_ROW_MAJOR_ORDER)
240 columnMajorOrder = 0;
241 break;
242 case CLO_COLUMNS:
243 if (scanned[iArg].n_items < 3)
244 SDDS_Bomb("invalid -columns syntax");
245 if (!(seqRequest = SDDS_Realloc(seqRequest, sizeof(*seqRequest) * (seqRequests + 1))))
246 SDDS_Bomb("memory allocation failure");
247 scanned[iArg].n_items -= 1;
248 memset(seqRequest + seqRequests, 0, sizeof(*seqRequest));
249 /* remove following pointer initialization because memset already initializes them */
250 seqRequest[seqRequests].randomizationGroup = -1;
251 seqRequest[seqRequests].factor = 1;
252 seqRequest[seqRequests].offset = 0;
253 if (!scanItemList(&seqRequest[seqRequests].flags,
254 scanned[iArg].list + 1, &scanned[iArg].n_items, 0,
255 "datafile", SDDS_STRING, &seqRequest[seqRequests].dataFileName, 1, SEQ_DATAFILE,
256 "independentvariable", SDDS_STRING, &seqRequest[seqRequests].indepName, 1, SEQ_INDEPNAME,
257 "cdf", SDDS_STRING, &seqRequest[seqRequests].CDFName, 1, SEQ_CDFNAME,
258 "df", SDDS_STRING, &seqRequest[seqRequests].DFName, 1, SEQ_DFNAME,
259 "output", SDDS_STRING, &seqRequest[seqRequests].outputName, 1, SEQ_OUTPUTNAME,
260 "units", SDDS_STRING, &seqRequest[seqRequests].units, 1, SEQ_UNITSGIVEN,
261 "haltonradix", SDDS_LONG, &seqRequest[seqRequests].haltonRadix, 1, SEQ_HALTONRADIX,
262 "haltonoffset", SDDS_LONG, &seqRequest[seqRequests].haltonOffset, 1, SEQ_HALTONOFFSET,
263 "randomize", -1, NULL, 0, SEQ_RANDOMIZE,
264 "group", SDDS_LONG, &seqRequest[seqRequests].randomizationGroup, 1, SEQ_RANDOMGROUP,
265 "factor", SDDS_DOUBLE, &seqRequest[seqRequests].factor, 1, 0,
266 "offset", SDDS_DOUBLE, &seqRequest[seqRequests].offset, 1, 0, NULL) ||
267 bitsSet(seqRequest[seqRequests].flags & (SEQ_INDEPNAME + SEQ_CDFNAME + SEQ_DFNAME)) != 2)
268 SDDS_Bomb("invalid -columns syntax");
269 if (seqRequest[seqRequests].flags & SEQ_RANDOMGROUP && seqRequest[seqRequests].randomizationGroup <= 0)
270 SDDS_Bomb("use a positive integer for the randomization group ID");
271 if (seqRequest[seqRequests].flags & SEQ_CDFNAME && seqRequest[seqRequests].flags & SEQ_DFNAME)
272 SDDS_Bomb("give df or cdf for -columns, not both");
273 if (seqRequest[seqRequests].flags & SEQ_HALTONRADIX && !is_prime(seqRequest[seqRequests].haltonRadix))
274 SDDS_Bomb("halton radix must be a prime number");
275 seqRequests++;
276 scanned[iArg].n_items += 1;
277 break;
278 case CLO_GAUSSIAN:
279 if (scanned[iArg].n_items < 2)
280 SDDS_Bomb("invalid -gaussian syntax");
281 if (!(seqRequest = SDDS_Realloc(seqRequest, sizeof(*seqRequest) * (seqRequests + 1))))
282 SDDS_Bomb("memory allocation failure");
283 memset(seqRequest + seqRequests, 0, sizeof(*seqRequest));
284 scanned[iArg].n_items -= 1;
285 seqRequest[seqRequests].randomizationGroup = -1;
286 seqRequest[seqRequests].mean = 0;
287 seqRequest[seqRequests].sigma = 1;
288 if (!scanItemList(&seqRequest[seqRequests].flags,
289 scanned[iArg].list + 1, &scanned[iArg].n_items, 0,
290 "columnName", SDDS_STRING, &seqRequest[seqRequests].outputName, 1, SEQ_OUTPUTNAME,
291 "meanValue", SDDS_STRING, &meanPar, 1, 0,
292 "sigmaValue", SDDS_STRING, &sigmaPar, 1, 0,
293 "units", SDDS_STRING, &seqRequest[seqRequests].units, 1, SEQ_UNITSGIVEN, NULL))
294 SDDS_Bomb("invalid -gaussian syntax");
295 seqRequest[seqRequests].flags |= SEQ_DIRECT_GAUSSIAN;
296 if (!(seqRequest[seqRequests].flags & SEQ_OUTPUTNAME) || !(seqRequest[seqRequests].outputName))
297 SDDS_Bomb("columnName is not provided for gaussian distribution/");
298 if (meanPar) {
299 if (wild_match(meanPar, "@*"))
300 SDDS_CopyString(&seqRequest[seqRequests].meanPar, meanPar + 1);
301 else if (!get_double(&seqRequest[seqRequests].mean, meanPar))
302 SDDS_Bomb("Invalid value given for mean value of -gaussian distribution.");
303 free(meanPar);
304 meanPar = NULL;
305 }
306 if (sigmaPar) {
307 if (wild_match(sigmaPar, "@*"))
308 SDDS_CopyString(&seqRequest[seqRequests].sigmaPar, sigmaPar + 1);
309 else if (!get_double(&seqRequest[seqRequests].sigma, sigmaPar))
310 SDDS_Bomb("Invalid value given for sigma value of -gaussian distribution.");
311 free(sigmaPar);
312 sigmaPar = NULL;
313 }
314 seqRequests++;
315 scanned[iArg].n_items += 1;
316 break;
317 case CLO_UNIFORM:
318 if (scanned[iArg].n_items < 2)
319 SDDS_Bomb("invalid -uniform syntax");
320 if (!(seqRequest = SDDS_Realloc(seqRequest, sizeof(*seqRequest) * (seqRequests + 1))))
321 SDDS_Bomb("memory allocation failure");
322 memset(seqRequest + seqRequests, 0, sizeof(*seqRequest));
323 scanned[iArg].n_items -= 1;
324 memset(seqRequest + seqRequests, 0, sizeof(*seqRequest));
325 seqRequest[seqRequests].randomizationGroup = -1;
326 seqRequest[seqRequests].min = 0;
327 seqRequest[seqRequests].max = 1;
328 if (!scanItemList(&seqRequest[seqRequests].flags,
329 scanned[iArg].list + 1, &scanned[iArg].n_items, 0,
330 "columnName", SDDS_STRING, &seqRequest[seqRequests].outputName, 1, SEQ_OUTPUTNAME,
331 "minimumValue", SDDS_STRING, &minPar, 1, 0,
332 "maximumValue", SDDS_STRING, &maxPar, 1, 0,
333 "units", SDDS_STRING, &seqRequest[seqRequests].units, 1, SEQ_UNITSGIVEN, NULL))
334 SDDS_Bomb("invalid -uniform syntax");
335 seqRequest[seqRequests].flags |= SEQ_DIRECT_UNIFORM;
336 if (!(seqRequest[seqRequests].flags & SEQ_OUTPUTNAME) ||
337 !(seqRequest[seqRequests].outputName))
338 SDDS_Bomb("columnName is not provided for uniform distribution/");
339 if (minPar) {
340 if (wild_match(minPar, "@*"))
341 SDDS_CopyString(&seqRequest[seqRequests].minPar, minPar + 1);
342 else if (!get_double(&seqRequest[seqRequests].min, minPar))
343 SDDS_Bomb("Invalid value given for minimum value of -uniform distribution.");
344 free(minPar);
345 minPar = NULL;
346 }
347 if (maxPar) {
348 if (wild_match(maxPar, "@*"))
349 SDDS_CopyString(&seqRequest[seqRequests].maxPar, maxPar + 1);
350 else if (!get_double(&seqRequest[seqRequests].max, maxPar))
351 SDDS_Bomb("Invalid value given for maximum value of -uniform distribution.");
352 free(maxPar);
353 maxPar = NULL;
354 }
355 seqRequests++;
356 scanned[iArg].n_items += 1;
357 break;
358 case CLO_POISSON:
359 if (scanned[iArg].n_items < 2)
360 SDDS_Bomb("invalid -poisson syntax");
361 if (!(seqRequest = SDDS_Realloc(seqRequest, sizeof(*seqRequest) * (seqRequests + 1))))
362 SDDS_Bomb("memory allocation failure");
363 memset(seqRequest + seqRequests, 0, sizeof(*seqRequest));
364 scanned[iArg].n_items -= 1;
365 memset(seqRequest + seqRequests, 0, sizeof(*seqRequest));
366 seqRequest[seqRequests].randomizationGroup = -1;
367 seqRequest[seqRequests].mean = 1;
368 if (!scanItemList(&seqRequest[seqRequests].flags,
369 scanned[iArg].list + 1, &scanned[iArg].n_items, 0,
370 "columnName", SDDS_STRING, &seqRequest[seqRequests].outputName, 1, SEQ_OUTPUTNAME,
371 "meanValue", SDDS_STRING, &meanPar, 1, 0,
372 "units", SDDS_STRING, &seqRequest[seqRequests].units, 1, SEQ_UNITSGIVEN, NULL))
373 SDDS_Bomb("invalid -poisson syntax");
374 seqRequest[seqRequests].flags |= SEQ_DIRECT_POISSON;
375 if (!(seqRequest[seqRequests].flags & SEQ_OUTPUTNAME) || !(seqRequest[seqRequests].outputName))
376 SDDS_Bomb("columnName is not provided for poisson distribution/");
377 if (meanPar) {
378 if (wild_match(meanPar, "@*"))
379 SDDS_CopyString(&seqRequest[seqRequests].meanPar, meanPar + 1);
380 else if (!get_double(&seqRequest[seqRequests].mean, meanPar))
381 SDDS_Bomb("Invalid value given for mean value of -poisson distribution.");
382 free(meanPar);
383 meanPar = NULL;
384 }
385 seqRequests++;
386 scanned[iArg].n_items += 1;
387 break;
388 case CLO_SAMPLES:
389 if (scanned[iArg].n_items != 2 ||
390 sscanf(scanned[iArg].list[1], "%ld", &samples) != 1 ||
391 samples <= 0)
392 SDDS_Bomb("invalid -samples syntax");
393 break;
394 case CLO_SEED:
395 if (scanned[iArg].n_items != 2 ||
396 sscanf(scanned[iArg].list[1], "%ld", &randomNumberSeed) != 1)
397 SDDS_Bomb("invalid -seed syntax");
398 break;
399 case CLO_PIPE:
400 if (!processPipeOption(scanned[iArg].list + 1, scanned[iArg].n_items - 1, &pipeFlags))
401 SDDS_Bomb("invalid -pipe syntax");
402 break;
403 case CLO_VERBOSE:
404 verbose = 1;
405 break;
406 case CLO_OPTIMAL_HALTON:
407 optimalHalton = 1;
408 break;
409 case CLO_THREADS:
410 if (scanned[iArg].n_items != 2 ||
411 sscanf(scanned[iArg].list[1], "%d", &threads) != 1 ||
412 threads < 1)
413 SDDS_Bomb("invalid -threads syntax");
414 break;
415 default:
416 fprintf(stderr, "error: unknown/ambiguous option: %s\n", scanned[iArg].list[0]);
417 exit(EXIT_FAILURE);
418 break;
419 }
420 } else {
421 if (!input)
422 input = scanned[iArg].list[0];
423 else if (!output)
424 output = scanned[iArg].list[0];
425 else
426 SDDS_Bomb("too many filenames seen");
427 }
428 }
429
430 if (!seqRequests)
431 SDDS_Bomb("give one or more -columns options");
432 if (samples < 1)
433 SDDS_Bomb("-samples option not given");
434
435 for (i = 0; i < seqRequests; i++) {
436 if (!(seqRequest[i].flags & (SEQ_DATAFILE | SEQ_DIRECT_GAUSSIAN | SEQ_DIRECT_UNIFORM | SEQ_DIRECT_POISSON)))
437 break;
438 }
439 if (i == seqRequests) {
440 /* all columns options have either their own input files or else use
441 * one of the "direct" distributions. Hence, we don't expect an input
442 * file.
443 */
444 if (!input)
445 pipeFlags |= USE_STDIN; /* not really, but fakes out processFilenames */
446 if (input && !output) {
447 output = input;
448 input = NULL;
449 pipeFlags |= USE_STDIN;
450 if (fexists(output)) {
451 sprintf(msgBuffer, "%s exists already (sddssampledist)", output);
452 SDDS_Bomb(msgBuffer);
453 }
454 }
455 }
456
457 processFilenames("sddssampledist", &input, &output, pipeFlags, 0, NULL);
458
459 if (!SDDS_InitializeOutput(&SDDSout, SDDS_BINARY, 0, NULL, NULL, output))
460 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
461
462 if (verbose)
463 fprintf(stderr, "Initialized output file %s\n", output);
464
465 /* open and check input files */
466 for (i = mainInputOpened = 0; i < seqRequests; i++) {
467 if (seqRequest[i].flags & SEQ_DIRECT_GAUSSIAN) {
468 if (seqRequest[i].meanPar || seqRequest[i].sigmaPar) {
469 if (!mainInputOpened) {
470 if (!SDDS_InitializeInput(&SDDSin, input) ||
471 !SDDS_TransferAllParameterDefinitions(&SDDSout, &SDDSin, 0))
472 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
473 mainInputOpened = 1;
474 }
475 requireInput = 1;
476 SDDSptr = &SDDSin;
477 if ((seqRequest[i].meanPar &&
478 SDDS_CheckParameter(SDDSptr, seqRequest[i].meanPar, NULL, SDDS_ANY_NUMERIC_TYPE, stderr) != SDDS_CHECK_OK) ||
479 (seqRequest[i].sigmaPar &&
480 SDDS_CheckParameter(SDDSptr, seqRequest[i].sigmaPar, NULL, SDDS_ANY_NUMERIC_TYPE, stderr) != SDDS_CHECK_OK)) {
481 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
482 exit(EXIT_FAILURE);
483 }
484 }
485 if (!SDDS_DefineSimpleColumn(&SDDSout, seqRequest[i].outputName, NULL, SDDS_DOUBLE))
486 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
487 } else if (seqRequest[i].flags & SEQ_DIRECT_UNIFORM) {
488 if (seqRequest[i].minPar || seqRequest[i].maxPar) {
489 if (!mainInputOpened) {
490 if (!SDDS_InitializeInput(&SDDSin, input) ||
491 !SDDS_TransferAllParameterDefinitions(&SDDSout, &SDDSin, 0))
492 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
493 mainInputOpened = 1;
494 }
495 requireInput = 1;
496 SDDSptr = &SDDSin;
497 if ((seqRequest[i].minPar &&
498 SDDS_CheckParameter(SDDSptr, seqRequest[i].minPar, NULL, SDDS_ANY_NUMERIC_TYPE, stderr) != SDDS_CHECK_OK) ||
499 (seqRequest[i].maxPar &&
500 SDDS_CheckParameter(SDDSptr, seqRequest[i].maxPar, NULL, SDDS_ANY_NUMERIC_TYPE, stderr) != SDDS_CHECK_OK)) {
501 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
502 exit(EXIT_FAILURE);
503 }
504 }
505 if (!SDDS_DefineSimpleColumn(&SDDSout, seqRequest[i].outputName, NULL, SDDS_DOUBLE))
506 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
507 } else if (seqRequest[i].flags & SEQ_DIRECT_POISSON) {
508 if (seqRequest[i].meanPar) {
509 if (!mainInputOpened) {
510 if (!SDDS_InitializeInput(&SDDSin, input) ||
511 !SDDS_TransferAllParameterDefinitions(&SDDSout, &SDDSin, 0))
512 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
513 mainInputOpened = 1;
514 }
515 requireInput = 1;
516 SDDSptr = &SDDSin;
517 if (SDDS_CheckParameter(SDDSptr, seqRequest[i].meanPar, NULL, SDDS_ANY_NUMERIC_TYPE, stderr) != SDDS_CHECK_OK) {
518 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
519 exit(EXIT_FAILURE);
520 }
521 }
522 if (!SDDS_DefineSimpleColumn(&SDDSout, seqRequest[i].outputName, NULL, SDDS_LONG))
523 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
524 } else {
525 if (seqRequest[i].flags & SEQ_RANDOMIZE) {
526 long newGroupID = 0;
527 /* define randomization groups */
528 if (seqRequest[i].flags & SEQ_RANDOMGROUP) {
529 newGroupID = seqRequest[i].randomizationGroup;
530 for (j = 0; j < randomizationGroups; j++)
531 if (randomizationData[j].group == newGroupID) {
532 newGroupID = 0;
533 break;
534 }
535 } else {
536 seqRequest[i].randomizationGroup = newGroupID = -(i + 1);
537 }
538 if (newGroupID != 0) {
539 if (!(randomizationData = SDDS_Realloc(randomizationData, sizeof(*randomizationData) * (randomizationGroups + 1))))
540 SDDS_Bomb("memory allocation failure");
541 randomizationData[randomizationGroups].group = newGroupID;
542 randomizationData[randomizationGroups].order = NULL;
543 randomizationGroups++;
544 }
545 }
546 if (seqRequest[i].flags & SEQ_DATAFILE) {
547 if (!SDDS_InitializeInput(&seqRequest[i].SDDSin, seqRequest[i].dataFileName))
548 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
549 SDDSptr = &seqRequest[i].SDDSin;
550 } else {
551 if (!mainInputOpened) {
552 if (!SDDS_InitializeInput(&SDDSin, input) ||
553 !SDDS_TransferAllParameterDefinitions(&SDDSout, &SDDSin, 0))
554 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
555 mainInputOpened = 1;
556 }
557 requireInput = 1;
558 SDDSptr = &SDDSin;
559 }
560 if (SDDS_CheckColumn(SDDSptr, seqRequest[i].indepName, NULL, SDDS_ANY_NUMERIC_TYPE, stderr) != SDDS_CHECK_OK ||
561 ((seqRequest[i].flags & SEQ_CDFNAME) &&
562 SDDS_CheckColumn(SDDSptr, seqRequest[i].CDFName, NULL, SDDS_ANY_NUMERIC_TYPE, stderr) != SDDS_CHECK_OK) ||
563 ((seqRequest[i].flags & SEQ_DFNAME) &&
564 SDDS_CheckColumn(SDDSptr, seqRequest[i].DFName, NULL, SDDS_ANY_NUMERIC_TYPE, stderr) != SDDS_CHECK_OK) ||
565 !SDDS_TransferColumnDefinition(&SDDSout, SDDSptr, seqRequest[i].indepName, seqRequest[i].outputName)) {
566 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
567 exit(EXIT_FAILURE);
568 }
569 }
570
571 if (seqRequest[i].flags & SEQ_UNITSGIVEN &&
572 !SDDS_ChangeColumnInformation(&SDDSout, "units", seqRequest[i].units, SDDS_SET_BY_NAME, seqRequest[i].outputName))
573 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
574 }
575
576 if (verbose)
577 fprintf(stderr, "Initialized input files\n");
578
579 if (columnMajorOrder != -1)
580 SDDSout.layout.data_mode.column_major = columnMajorOrder;
581 else
582 SDDSout.layout.data_mode.column_major = 0;
583
584 if (!SDDS_WriteLayout(&SDDSout))
585 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
586
587 if (randomNumberSeed <= 0) {
588 randomNumberSeed = (long)time((time_t *)NULL);
589 randomNumberSeed = 2 * (randomNumberSeed / 2) + 1;
590#if defined(_WIN32) || defined(__APPLE__)
591 random_1(-labs(randomNumberSeed));
592#else
593 random_1(-fabs(randomNumberSeed));
594#endif
595 } else
596 random_1(-randomNumberSeed);
597
598 if (!(sample = calloc(samples, sizeof(*sample))))
599 SDDS_Bomb("memory allocation failure");
600 if (!(dataFilePageStatus = calloc(seqRequests, sizeof(*dataFilePageStatus))))
601 SDDS_Bomb("memory allocation failure");
602 while (1) {
603 if (verbose)
604 fprintf(stderr, "Beginning page loop\n");
605 if (input && SDDS_ReadPage(&SDDSin) <= 0)
606 break;
607 for (i = 0; i < seqRequests; i++)
608 dataFilePageStatus[i] = 1;
609#pragma omp parallel for if (threads > 1 && seqRequests > 1) num_threads(threads) schedule(dynamic)
610 for (i = 0; i < seqRequests; i++) {
611 if (seqRequest[i].flags & SEQ_DATAFILE)
612 dataFilePageStatus[i] = SDDS_ReadPage(&seqRequest[i].SDDSin);
613 }
614 for (i = 0; i < seqRequests; i++)
615 if ((seqRequest[i].flags & SEQ_DATAFILE) && dataFilePageStatus[i] <= 0)
616 break;
617 if (i != seqRequests)
618 break;
619 if (!SDDS_StartPage(&SDDSout, samples) || (input && !SDDS_CopyParameters(&SDDSout, &SDDSin)))
620 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
621 if (verbose)
622 fprintf(stderr, "Defining randomization tables\n");
623 /* define randomization tables */
624 for (i = 0; i < randomizationGroups; i++) {
625 if (!(randomizationData[i].order = SDDS_Malloc(sizeof(*randomizationData[i].order) * samples)))
626 SDDS_Bomb("memory allocation failure");
627 for (j = 0; j < samples; j++)
628 randomizationData[i].order[j] = j;
629 randomizeOrder((char *)randomizationData[i].order, sizeof(*randomizationData[i].order), samples, 0, random_1);
630 }
631 if (verbose)
632 fprintf(stderr, "Beginning loop over sequence requests\n");
633 for (i = 0; i < seqRequests; i++) {
634 if (verbose)
635 fprintf(stderr, "Processing sequence request %ld\n", i);
636 if (seqRequest[i].flags & SEQ_DIRECT_GAUSSIAN) {
637 if ((seqRequest[i].meanPar &&
638 !SDDS_GetParameterAsDouble(&SDDSin, seqRequest[i].meanPar, &seqRequest[i].mean)) ||
639 (seqRequest[i].sigmaPar &&
640 !SDDS_GetParameterAsDouble(&SDDSin, seqRequest[i].sigmaPar, &seqRequest[i].sigma)))
641 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
642 for (j = 0; j < samples; j++)
643 sample[j] = gauss_rn_lim(seqRequest[i].mean, seqRequest[i].sigma, -1, random_1);
644 } else if (seqRequest[i].flags & SEQ_DIRECT_UNIFORM) {
645 if ((seqRequest[i].minPar &&
646 !SDDS_GetParameterAsDouble(&SDDSin, seqRequest[i].minPar, &seqRequest[i].min)) ||
647 (seqRequest[i].maxPar &&
648 !SDDS_GetParameterAsDouble(&SDDSin, seqRequest[i].maxPar, &seqRequest[i].max)))
649 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
650 for (j = 0; j < samples; j++)
651 sample[j] = seqRequest[i].min + (seqRequest[i].max - seqRequest[i].min) * random_1(1);
652 } else if (seqRequest[i].flags & SEQ_DIRECT_POISSON) {
653 double *pos_x, *pos_cdf, CDF;
654 long pos_points, code;
655 pos_x = pos_cdf = NULL;
656 if ((seqRequest[i].meanPar &&
657 !SDDS_GetParameterAsDouble(&SDDSin, seqRequest[i].meanPar, &seqRequest[i].mean)))
658 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
659 pos_points = CreatePoissonDistributionTable(&pos_x, &pos_cdf, seqRequest[i].mean);
660
661 for (j = 0; j < samples; j++) {
662 CDF = random_1(1);
663 sample[j] = (int)(interp(pos_x, pos_cdf, pos_points, CDF, 0, 1, &code));
664 /* fprintf(stderr, "%ld, cdf=%f, sample=%f\n", j, CDF, sample[j]); */
665 }
666 free(pos_x);
667 free(pos_cdf);
668 } else {
669 if (input && !(seqRequest[i].flags & SEQ_DATAFILE))
670 SDDSptr = &SDDSin;
671 else
672 SDDSptr = &seqRequest[i].SDDSin;
673 if ((values = SDDS_CountRowsOfInterest(SDDSptr))) {
674 if (!(IVValue = SDDS_GetColumnInDoubles(SDDSptr, seqRequest[i].indepName)) ||
675 (seqRequest[i].flags & SEQ_CDFNAME &&
676 !(CDFValue = SDDS_GetColumnInDoubles(SDDSptr, seqRequest[i].CDFName))) ||
677 (seqRequest[i].flags & SEQ_DFNAME &&
678 !(CDFValue = SDDS_GetColumnInDoubles(SDDSptr, seqRequest[i].DFName))))
679 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
680 } else {
681 sprintf(msgBuffer, "empty page for file %s\n",
682 (seqRequest[i].flags & SEQ_DATAFILE) ? seqRequest[i].dataFileName : input);
683 SDDS_Bomb(msgBuffer);
684 }
685 if (verbose)
686 fprintf(stderr, "Checking and converting CDF/DF values\n");
687 /* check/convert CDF/DF values */
688 for (j = 1; j < values; j++) {
689 if (IVValue[j - 1] > IVValue[j]) {
690 sprintf(msgBuffer, "random variate values not monotonically increasing for %s",
691 (seqRequest[i].flags & SEQ_DATAFILE) ? seqRequest[i].dataFileName : input);
692 SDDS_Bomb(msgBuffer);
693 }
694 if (seqRequest[i].flags & SEQ_DFNAME)
695 /* convert DF to CDF */
696 CDFValue[j] += CDFValue[j - 1];
697 if (CDFValue[j] < CDFValue[j - 1]) {
698 sprintf(msgBuffer, "CDF values decreasing for %s",
699 (seqRequest[i].flags & SEQ_DATAFILE) ? seqRequest[i].dataFileName : input);
700 SDDS_Bomb(msgBuffer);
701 }
702 }
703 if (verbose)
704 fprintf(stderr, "Normalizing CDF\n");
705 /* normalize the CDF */
706 if (CDFValue[values - 1] <= 0) {
707 sprintf(msgBuffer, "CDF not valid for %s\n", seqRequest[i].dataFileName);
708 SDDS_Bomb(msgBuffer);
709 }
710 for (j = 0; j < values; j++)
711 CDFValue[j] /= CDFValue[values - 1];
712 if (seqRequest[i].flags & SEQ_HALTONRADIX) {
713 if (verbose)
714 fprintf(stderr, "Starting halton sequence, offset=%" PRId32 "\n", seqRequest[i].haltonOffset);
715 if (!optimalHalton)
716 haltonID = startHaltonSequence(&seqRequest[i].haltonRadix, 0.5);
717 else
718 haltonID = startModHaltonSequence(&seqRequest[i].haltonRadix, 0);
719 while (seqRequest[i].haltonOffset-- > 0) {
720 if (!optimalHalton)
721 nextHaltonSequencePoint(haltonID);
722 else
724 }
725 }
726 if (verbose)
727 fprintf(stderr, "Generating samples\n");
728 for (j = 0; j < samples; j++) {
729 double CDF;
730 long code;
731 while (1) {
732 if (seqRequest[i].flags & SEQ_HALTONRADIX) {
733 if (!optimalHalton)
734 CDF = nextHaltonSequencePoint(haltonID);
735 else
736 CDF = nextModHaltonSequencePoint(haltonID);
737 } else
738 CDF = random_1(1);
739 if (CDF <= CDFValue[values - 1] && CDF >= CDFValue[0])
740 break;
741 }
742 sample[j] = seqRequest[i].factor * interp(IVValue, CDFValue, values, CDF, 0, 1, &code) + seqRequest[i].offset;
743 }
744 if (seqRequest[i].flags & SEQ_RANDOMIZE) {
745 long k, l;
746 double *sample1;
747 if (verbose)
748 fprintf(stderr, "Randomizing order of values\n");
749 if (!(sample1 = malloc(sizeof(*sample1) * samples)))
750 SDDS_Bomb("memory allocation failure");
751 for (l = 0; l < randomizationGroups; l++)
752 if (randomizationData[l].group == seqRequest[i].randomizationGroup)
753 break;
754 if (l == randomizationGroups)
755 SDDS_Bomb("problem with construction of randomization groups!");
756 for (k = 0; k < samples; k++)
757 sample1[k] = sample[randomizationData[l].order[k]];
758 free(sample);
759 sample = sample1;
760 }
761 free(IVValue);
762 free(CDFValue);
763 }
764 if (verbose)
765 fprintf(stderr, "Setting SDDS column values\n");
766 if (!SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_NAME, sample, samples,
767 seqRequest[i].outputName ? seqRequest[i].outputName : seqRequest[i].indepName))
768 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
769 }
770 if (verbose)
771 fprintf(stderr, "Writing data page\n");
772 if (!SDDS_WritePage(&SDDSout))
773 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
774 if (!requireInput)
775 break;
776 }
777 if (verbose)
778 fprintf(stderr, "Exited read loop\n");
779 free(sample);
780 free(dataFilePageStatus);
781 if ((input && !SDDS_Terminate(&SDDSin)) || !SDDS_Terminate(&SDDSout))
782 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
783 for (i = 0; i < seqRequests; i++) {
784 if (seqRequest[i].dataFileName)
785 free(seqRequest[i].dataFileName);
786 if (seqRequest[i].indepName)
787 free(seqRequest[i].indepName);
788 if (seqRequest[i].outputName)
789 free(seqRequest[i].outputName);
790 if (seqRequest[i].DFName)
791 free(seqRequest[i].DFName);
792 if (seqRequest[i].CDFName)
793 free(seqRequest[i].CDFName);
794 if (seqRequest[i].units)
795 free(seqRequest[i].units);
796 if (seqRequest[i].meanPar)
797 free(seqRequest[i].meanPar);
798 if (seqRequest[i].sigmaPar)
799 free(seqRequest[i].sigmaPar);
800 if (seqRequest[i].minPar)
801 free(seqRequest[i].minPar);
802 if (seqRequest[i].maxPar)
803 free(seqRequest[i].maxPar);
804 if (seqRequest[i].flags & SEQ_DATAFILE && !SDDS_Terminate(&(seqRequest[i].SDDSin)))
805 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
806 }
807 free(seqRequest);
808 for (i = 0; i < randomizationGroups; i++)
809 free(randomizationData[i].order);
810 if (randomizationData)
811 free(randomizationData);
812
813 free_scanargs(&scanned, argc);
814
815 return EXIT_SUCCESS;
816}
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_SetColumnFromDoubles(SDDS_DATASET *SDDS_dataset, int32_t mode, double *data, int64_t rows,...)
Sets the values for a single data column using double-precision floating-point numbers.
double * SDDS_GetParameterAsDouble(SDDS_DATASET *SDDS_dataset, char *parameter_name, double *memory)
Retrieves the value of a specified parameter as a double from the current data table of an SDDS datas...
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_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_DefineSimpleColumn(SDDS_DATASET *SDDS_dataset, const char *name, const char *unit, int32_t type)
Defines a simple data column within the SDDS dataset.
int32_t SDDS_WritePage(SDDS_DATASET *SDDS_dataset)
Writes the current data table to the output file.
int32_t SDDS_WriteLayout(SDDS_DATASET *SDDS_dataset)
Writes the SDDS layout header to the output file.
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_CheckColumn(SDDS_DATASET *SDDS_dataset, char *name, char *units, int32_t type, FILE *fp_message)
Checks if a column exists in the SDDS dataset with the specified name, units, and type.
void SDDS_PrintErrors(FILE *fp, int32_t mode)
Prints recorded error messages to a specified file stream.
Definition SDDS_utils.c:474
void * SDDS_Malloc(size_t size)
Allocates memory of a specified size.
Definition SDDS_utils.c:705
void SDDS_RegisterProgramName(const char *name)
Registers the executable program name for use in error messages.
Definition SDDS_utils.c:318
int32_t SDDS_CheckParameter(SDDS_DATASET *SDDS_dataset, char *name, char *units, int32_t type, FILE *fp_message)
Checks if a parameter exists in the SDDS dataset with the specified name, units, and type.
int32_t SDDS_CopyString(char **target, const char *source)
Copies a source string to a target string with memory allocation.
Definition SDDS_utils.c:922
#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_ANY_NUMERIC_TYPE
Special identifier used by SDDS_Check*() routines to accept any numeric type.
Definition SDDStypes.h:157
#define SDDS_DOUBLE
Identifier for the double data type.
Definition SDDStypes.h:37
long bitsSet(unsigned long data)
Counts the number of set bits (1s) in the given data.
Definition binary.c:52
int get_double(double *dptr, char *s)
Parses a double value from the given string.
Definition data_scan.c:40
double random_1(long iseed)
Generate a uniform random double in [0,1] using a custom seed initialization.
Definition drand.c:230
long randomizeOrder(char *ptr, long size, long length, long iseed, double(*urandom)(long iseed1))
Randomize the order of an array of elements.
Definition drand.c:520
double gauss_rn_lim(double mean, double sigma, double limit_in_sigmas, double(*urandom)(long iseed))
Generate a Gaussian-distributed random number with specified mean, sigma, and optional cutoff.
Definition drand.c:433
int64_t is_prime(int64_t number)
Determine if a number is prime.
Definition factorize.c:26
long fexists(const char *filename)
Checks if a file exists.
Definition fexists.c:27
int32_t startModHaltonSequence(int32_t *radix, double tiny)
Start a modified Halton sequence.
Definition halton.c:509
double nextModHaltonSequencePoint(long ID)
Retrieve the next point from the modified Halton sequence.
Definition halton.c:625
int32_t startHaltonSequence(int32_t *radix, double value)
Initialize and start a new Halton sequence.
Definition halton.c:38
double nextHaltonSequencePoint(long ID)
Get the next point in a Halton sequence.
Definition halton.c:107
double interp(double *f, double *x, long n, double xo, long warnings, long order, long *returnCode)
Performs simple linear interpolation of data.
Definition interp.c:34
long match_string(char *string, char **option, long n_options, long mode)
Matches a given string against an array of option strings based on specified modes.
int scanargs(SCANNED_ARG **scanned, int argc, char **argv)
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
void free_scanargs(SCANNED_ARG **scanned, int argc)
Definition scanargs.c:588
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.
int wild_match(char *string, char *template)
Determine whether one string is a wildcard match for another.
Definition wild_match.c:49