StatMech
Loading...
Searching...
No Matches
maximum.c File Reference

Functions

int main (int argc, char **argv)
 

Function Documentation

◆ main()

int main ( int  argc,
char **  argv 
)
16 {
17 if (argc < 3) {
18 printf("Usage: 1.This 2.n 3.ndata\n");
19 exit(EXIT_FAILURE);
20 }
21 Integer_t i, j, count;
22 Integer_t const n = atoi(argv[1]);
23 Integer_t const ndata = atoi(argv[2]);
24 Integer_t seed=0;
25 double max, maxData;
26 double rand, rand1, rand2;
27 double x1 = 2;
28 double thereshold = x1/(1.0+x1);
29 double coeff = exp(x1)/(1+x1);
30
31 Integer_t subMax;
32 Integer_t *histogram;
33 double mean, stddev;
34 double binwidth=0.05;
35
36 char str[100];
37 char *filename = str;
38
39 FILE *Infp;
40 FILE *Outfp;
41 sprintf(str, "Dist_maximum_%lld_raw.txt", n);
42 Outfp = fopen(filename, "w");
43 if (Outfp == NULL) {
44 printf("Can't open a file %s\n", filename);
45 exit(EXIT_FAILURE);
46 }
47 fprintf(Outfp, "#(n=%lld, seed=%lld) Distribution of maximum of |x| \n\n", n, seed);
48
49 init_genrand(seed);
50 maxData = NAN;
51 mean = 0;
52 stddev = 0;
53 for(i = 0;i < ndata; ++i) {
54 max = NAN;
55 for(j = 0;j < n; ++j) {
56 rand1 = genrand_real3();
57 if(rand1 < thereshold) rand = (1.0+x1)*rand1;
58 else rand = x1-log((1-rand1)*(1.0+x1));
59 max = fmax(max, rand);
60 }
61 maxData = fmax(maxData, max);
62 mean += max;
63 stddev += max*max;
64 fprintf(Outfp, "%lld %le\n", i, max);
65 }
66 mean /= (double)ndata;
67 stddev /= (double)ndata;
68 stddev = sqrt( stddev -mean*mean );
69 binwidth *= stddev;
70 subMax = (Integer_t)(maxData/binwidth) +1;
71 histogram = alloc_ivector(subMax); for(i = 0;i < subMax; ++i) histogram[i] = 0;
72
73 sprintf(str, "Dist_maximum_%lld_raw.txt", n);
74 Infp = fopen(filename, "r");
75 if (Infp == NULL) {
76 printf("Can't open a file %s\n", filename);
77 exit(EXIT_FAILURE);
78 }
79 skip_header(Infp);
80
81 count=0;
82 while ( fscanf(Infp, "%lld %le", &i, &max) != EOF ) {
83 histogram[ (Integer_t)(max/binwidth) ] += 1;
84 count += 1;
85 }
86
87 sprintf(str, "Dist_maximum_%lld.txt", n);
88 Outfp = fopen(filename, "w");
89 if (Outfp == NULL) {
90 printf("Can't open a file %s\n", filename);
91 exit(EXIT_FAILURE);
92 }
93 fprintf(Outfp, "# mean=%le, stddev=%le\n\n", mean, stddev);
94 fprintf(Outfp, "# 1.Maximum value 2.Probability density\n\n");
95 for(i = 0;i < subMax; ++i) {
96 fprintf(Outfp, "%le %le\n", binwidth*(i+0.5), histogram[i]/(binwidth*count));
97 }
98
99 return 0;
100}
MKL_INT Integer_t
Definition mytypes.hpp:359