16 {
17 if (argc < 3) {
18 printf("Usage: 1.This 2.ndata 3.InFile\n");
19 exit(EXIT_FAILURE);
20 }
25
26 double max, maxData, dtemp;
27 double rand, rand1, rand2;
28 double CDF[MAX];
for(i = 0;i < MAX; ++i)
CDF[i] = 0;
29
32 double mean, stddev, mode;
33 double binwidth=0.05;
34
36 char str[BUFSIZE];
37 char *filename = str;
38
39 FILE *Infp;
40 strcpy(str, argv[2]);
41 Infp = fopen(filename, "r");
42 if (Infp == NULL) {
43 printf("Can't open a file %s\n", filename);
44 exit(EXIT_FAILURE);
45 }
46 skip_header(Infp);
47 while( fscanf(Infp, "%lld %lld %le", &i, &j, &dtemp) != EOF ) {
49 }
50
51 FILE *Outfp;
52 sprintf(str, "Dist_maximum_raw.txt");
53 Outfp = fopen(filename, "w");
54 if (Outfp == NULL) {
55 printf("Can't open a file %s\n", filename);
56 exit(EXIT_FAILURE);
57 }
58 fprintf(Outfp, "#(n=%lld, seed=%lld) Distribution of maximum of |x| \n\n", n, seed);
59
60 init_genrand(seed);
61 maxData = NAN;
62 mean = 0;
63 stddev = 0;
64 for(i = 0;i < ndata; ++i) {
65
66 rand1 = genrand_real3();
67 for(j = 0;j < MAX; ++j)
if(
CDF[j] > rand1)
break;
68 n = j-1;
69
70 max = NAN;
71 for(j = 0;j < n; ++j) {
72 rand1 = genrand_real3();
73 rand = -log(rand1);
74 max = fmax(max, rand);
75
76 }
77 maxData = fmax(maxData, max);
78 mean += max;
79 stddev += max*max;
80 fprintf(Outfp, "%lld %le\n", i, max);
81 }
82 mean /= (double)ndata;
83 stddev /= (double)ndata;
84 stddev = sqrt( stddev -mean*mean );
85 binwidth *= stddev;
86 subMax = (
Integer_t)(maxData/binwidth) +1;
87 histogram = alloc_ivector(subMax); for(i = 0;i < subMax; ++i) histogram[i] = 0;
88
89 sprintf(str, "Dist_maximum_raw.txt");
90 Infp = fopen(filename, "r");
91 if (Infp == NULL) {
92 printf("Can't open a file %s\n", filename);
93 exit(EXIT_FAILURE);
94 }
95 skip_header(Infp);
96
97 count=0;
98 while ( fscanf(Infp, "%lld %le", &i, &max) != EOF ) {
99 histogram[ (
Integer_t)(max/binwidth) ] += 1;
100 count += 1;
101 }
102 maxCount = 0;
103 for(i = 0;i < subMax; ++i) {
104 if(maxCount < histogram[i]) {
105 mode = binwidth*(i+0.5);
106 maxCount = histogram[i];
107 }
108 }
109
110 sprintf(str, "Dist_maximum.txt");
111 Outfp = fopen(filename, "w");
112 if (Outfp == NULL) {
113 printf("Can't open a file %s\n", filename);
114 exit(EXIT_FAILURE);
115 }
116 fprintf(Outfp, "# Input: %s\n", argv[2]);
117 fprintf(Outfp, "# mean=%le, stddev=%le, mode=%le\n", mean, stddev, mode);
118 fprintf(Outfp, "# 1.Maximum value 2.Probability density\n\n");
119 for(i = 0;i < subMax; ++i) {
120 fprintf(Outfp, "%le %le\n", binwidth*(i+0.5), histogram[i]/(binwidth*count));
121 }
122
123 return 0;
124}
double CDF(double x)
Definition maximum.c:16
MKL_INT Integer_t
Definition mytypes.hpp:359