trandcounts.c - numtools - perform numerical operations on vectors and matrices in unix pipes
(HTM) git clone git://src.adamsgaard.dk/numtools
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
(DIR) LICENSE
---
trandcounts.c (2893B)
---
1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <err.h>
4 #include <limits.h>
5 #include <string.h>
6 #include <math.h>
7 #include <sys/time.h>
8 #include <unistd.h>
9
10 #include "util.h"
11
12 static void
13 usage(void)
14 {
15 errx(1, "usage: randcounts [-d delimstr] [-n] [-N num] [-p prec]"
16 "[-r repeats] [-R] [-s seed] weight1 [weight2 ...]");
17 }
18
19 int
20 main(int argc, char *argv[])
21 {
22 int i, ch, nbins, prec = 17, finalnl = 1, s = 0;
23 long j, seed, *counts = NULL, n = 1, r, repeats = 1, ratio = 0;
24 double val = 0.0, weightsum = 0.0, *weights = NULL;
25 char *delimstr = "\t";
26 const char *errstr;
27 struct timeval t1;
28
29 if (pledge("stdio", NULL) == -1)
30 err(2, "pledge");
31
32 while ((ch = getopt(argc, argv, "d:nN:r:Rp:s:")) != -1) {
33 switch (ch) {
34 case 'd':
35 delimstr = optarg;
36 break;
37 case 'n':
38 finalnl = 0;
39 break;
40 case 'N':
41 n = strtonum(optarg, 0, LONG_MAX, &errstr);
42 if (errstr != NULL)
43 errx(1, "bad num value, %s: %s", errstr, optarg);
44 break;
45 case 'r':
46 repeats = strtonum(optarg, 0, LONG_MAX, &errstr);
47 if (errstr != NULL)
48 errx(1, "bad repeats value, %s: %s", errstr, optarg);
49 break;
50 case 'R':
51 ratio = 1;
52 break;
53 case 'p':
54 prec = strtonum(optarg, 0, INT_MAX, &errstr);
55 if (errstr != NULL)
56 errx(1, "bad precision value, %s: %s", errstr, optarg);
57 break;
58 case 's':
59 seed = strtonum(optarg, 0, INT_MAX, &errstr);
60 if (errstr != NULL)
61 errx(1, "bad seed value, %s: %s", errstr, optarg);
62 break;
63 default:
64 usage();
65 }
66 }
67 argc -= optind;
68 argv += optind;
69 if (argc < 1)
70 usage();
71
72 nbins = argc;
73 if (!(weights = calloc(nbins, sizeof(double))) ||
74 !(counts = calloc(nbins, sizeof(long))))
75 err(1, "calloc");
76
77 for (i = 0; i < nbins; i++) {
78 if (!sscanf(argv[i], "%lf", &weights[i]))
79 errx(1, "bad weight value: %s", argv[i]);
80 if (weights[i] <= 0.0)
81 errx(1, "weight %d is not positive (%g)", i, weights[i]);
82 if (weights[i] > 1.0)
83 errx(1, "weight %d is greater than 1 (%g)", i, weights[i]);
84 }
85
86 for (i = 0; i < nbins; i++)
87 weightsum += weights[i];
88 if (fabs(weightsum - 1.0) > 1e-3)
89 errx(1, "weights do not sum to 1 (%g)", weightsum);
90
91 if (s)
92 #ifdef __OpenBSD__
93 srand48_deterministic(seed);
94 #else
95 srand48(seed);
96 else {
97 gettimeofday(&t1, NULL);
98 srand48(t1.tv_sec * t1.tv_usec); /* Caveat: identical seed for each ms */
99 }
100 #endif
101
102 for (r = 0; r < repeats; r++) {
103 for (i = 0; i < nbins; i++)
104 counts[i] = 0;
105 for (j = 0; j < n; j++) {
106 val = drand48();
107 weightsum = 0.0;
108 for (i = 0; i < nbins; i++) {
109 weightsum += weights[i];
110 if (val <= weightsum) {
111 counts[i]++;
112 break;
113 }
114 }
115 }
116 for (i = 0; i < nbins; i++) {
117 if (ratio)
118 printf("%.*g", prec, (double)counts[i] / n);
119 else
120 printf("%ld", counts[i]);
121 if (i < nbins - 1)
122 fputs(delimstr, stdout);
123 else
124 if (finalnl)
125 putchar('\n');
126 }
127 }
128
129 free(counts);
130 free(weights);
131
132 return 0;
133 }