StatMech
Loading...
Searching...
No Matches
Functions | Variables
Extreme_Compare.cpp File Reference

Functions

double residual (double a, Integer_t const NExt, double const *ExVal, double const *ExCDF, Integer_t const Ndata, double const *DatVal, double const *DatCDF)
 
int main (int argc, char **argv)
 

Variables

double DM
 

Function Documentation

◆ main()

int main ( int  argc,
char **  argv 
)
30 {
31 if(argc != 3) {
32 printf("Usage: 1.This 2.DataDir 3.ExtremeDir\n");
33 exit(EX_USAGE);
34 }
35 double const epsilon = 1.0e-8;
36 Integer_t Nmin=100, Nmax=0, n, DimSell, Ndata, total, itemp;
37 Integer_t NExtreme, NData;
38 double value, cdf, norm, mean, stddev, dtemp;
39 double amin, amid, amax;
40 double ResMin, ResMid, ResMax;
41
42 const char* DataDir=argv[1];
43 const char* ExtremeDir=argv[2];
44 double *ExtremeVal, *ExtremeCDF, ExtremeMean;
45 double *DataVal, *DataCDF, DataMean;
46
47 Integer_t const BUFSIZE = 1000;
48 char str[BUFSIZE], stemp[BUFSIZE];
49 char *filename = str;
50
51 FILE *StatFp, *DShellFp, *HistFp, *ExtremeFp;
52 strcpy(str, DataDir);
53 sprintf(stemp, "/Stat_MaxFluc.txt");
54 strcat(str, stemp);
55 StatFp = fopen(filename, "r");
56 if( InputFileCheck(stderr, StatFp, filename) != 0 ) exit(EX_NOINPUT);
57 if( skip_header(StatFp)!=0 ) exit(EX_NOINPUT);
58
59 while( fscanf(StatFp, "%lld %le %le %lld %le", &n, &dtemp, &dtemp, &itemp, &dtemp) != EOF ) {
60 if(n > Nmax) Nmax = n; if(n < Nmin) Nmin = n;
61 }
62 // printf("Nmin=%lld, Nmax=%lld\n", Nmin, Nmax);
63 fclose(StatFp);
64
65 for(n = Nmin;n <= Nmax; ++n) {
66 strcpy(str, DataDir);
67 sprintf(stemp, "/Hist_DimShell_n%lld.txt", n);
68 strcat(str, stemp);
69 DShellFp = fopen(filename, "r");
70 if( InputFileCheck(stderr, DShellFp, filename) != 0 ) continue;
71 if( skip_header(DShellFp)!=0 ) continue;
72 total = 0; mean = 0;
73 while( fscanf(DShellFp, "%lld %lld %le", &DimSell, &Ndata, &dtemp) != EOF ) {
74 total += Ndata;
75 mean += (double)DimSell*Ndata;
76 }
77 mean /= (double)total;
78 // printf("Mean of DimShell = %lld\n", (Integer_t)mean);
79 fclose(DShellFp);
80
81
82 strcpy(str, ExtremeDir);
83 sprintf(stemp, "/Dist_maximum_%lld.txt", (Integer_t)mean);
84 strcat(str, stemp);
85 ExtremeFp = fopen(filename, "r");
86 if( InputFileCheck(stderr, ExtremeFp, filename) != 0 ) continue;
87 if( skip_header(ExtremeFp)!=0 ) continue;
88 NExtreme = 0;
89 while( fscanf(ExtremeFp, "%le %le %le", &dtemp, &dtemp, &dtemp) != EOF ) NExtreme += 1;
90 rewind(ExtremeFp); if( skip_header(ExtremeFp)!=0 ) continue;
91
92 ExtremeVal = alloc_dvector(NExtreme);
93 ExtremeCDF = alloc_dvector(NExtreme);
94 total = 0; norm = 0; ExtremeMean = 0;
95 while( fscanf(ExtremeFp, "%le %le %le", &value, &dtemp, &cdf) != EOF ) {
96 if(total >= NExtreme) exit(EXIT_FAILURE);
97 ExtremeVal[total] = value;
98 ExtremeCDF[total] = 1-cdf;
99 ExtremeMean += ExtremeVal[total]*dtemp;
100 norm += dtemp;
101 total += 1;
102 }
103 ExtremeMean /= norm;
104 // printf("ExtremeMean = %le, NExtreme = %lld\n", ExtremeMean, NExtreme);
105 fclose(ExtremeFp);
106
107
108 strcpy(str, DataDir);
109 sprintf(stemp, "/Hist_MaxFluc_n%lld.txt", n);
110 strcat(str, stemp);
111 HistFp = fopen(filename, "r");
112 if( InputFileCheck(stderr, HistFp, filename) != 0 ) continue;
113 if( skip_header(HistFp)!=0 ) continue;
114 NData = 0;
115 while( fscanf(HistFp, "%le %le %le", &value, &dtemp, &cdf) != EOF ) NData += 1;
116 rewind(HistFp); if( skip_header(HistFp)!=0 ) continue;
117
118 DataVal = alloc_dvector(NData);
119 DataCDF = alloc_dvector(NData);
120 total = 0; norm = 0; DataMean = 0; stddev = 0;
121 while( fscanf(HistFp, "%le %le %le", &value, &dtemp, &cdf) != EOF ) {
122 if(total >= NData) exit(EXIT_FAILURE);
123 total += 1;
124 DataVal[total] = value;
125 DataCDF[total] = 1-cdf;
126 DataMean += DataVal[total]*dtemp;
127 // stddev += DataVal[total]*DataVal[total]*dtemp;
128 norm += dtemp;
129 // printf("DataVal[%lld] = %le, dtemp = %le\n", total, DataVal[total], dtemp);
130 }
131 DataMean /= norm;
132 // stddev /= norm; stddev=sqrt( stddev -DataMean*DataMean );
133 DM = DataMean +2*stddev;
134 // printf("DataMean = %le, NData = %lld\n", DataMean, NData);
135 fclose(HistFp);
136
137
138 amid = DataMean/ExtremeMean;
139 amin = 0.95*amid; amax = 1.05*amid;
140 total = 0;
141 // printf("DataMean/ExtremeMean = %le, amin = %le, amax = %le\n", amid, amin, amax);
142 while(1==1) {
143 amid = (amax+amin)/2.0;
144 ResMin = residual(amin, NExtreme, ExtremeVal, ExtremeCDF, NData, DataVal, DataCDF);
145 ResMid = residual(amid, NExtreme, ExtremeVal, ExtremeCDF, NData, DataVal, DataCDF);
146 ResMax = residual(amax, NExtreme, ExtremeVal, ExtremeCDF, NData, DataVal, DataCDF);
147
148 // printf("n=%lld, total=%lld, amin=%le, amid=%le, amax=%le, ResMin=%le, ResMid=%le, ResMax=%le\r", n, total, amin, amid, amax, ResMin, ResMid, ResMax);
149 if(ResMin>ResMid && ResMid>ResMax) amax *= 1.1;
150 else if(ResMin<ResMid && ResMid<ResMax) amin *= 0.9;
151 else if(ResMin<ResMid && ResMid>ResMax) { amid = DataMean/ExtremeMean; amin += 0.1*amid; amax += 0.1*amid; }
152 else break;
153 }
154 // printf("\n");
155 while( fabs(amax-amin) > epsilon ) {
156 amid = (amax+amin)/2.0;
157 ResMin = residual(amin, NExtreme, ExtremeVal, ExtremeCDF, NData, DataVal, DataCDF);
158 ResMid = residual(amid, NExtreme, ExtremeVal, ExtremeCDF, NData, DataVal, DataCDF);
159 ResMax = residual(amax, NExtreme, ExtremeVal, ExtremeCDF, NData, DataVal, DataCDF);
160
161 dtemp = (amax+amid)/2.0;
162 if(ResMid < residual(dtemp, NExtreme, ExtremeVal, ExtremeCDF, NData, DataVal, DataCDF)) amax = dtemp;
163 else amin = amid;
164 // dtemp = (amin+amid)/2.0;
165 // if(ResMid < residual(dtemp, NExtreme, ExtremeVal, ExtremeCDF, NData, DataVal, DataCDF)) amin = dtemp;
166 // else amax = amid;
167 total += 1;
168 // printf("n=%lld, total=%lld, amin=%le, amid=%le, amax=%le\r", n, total, amin, amid, amax);
169 }
170 printf("%lld %le %lld\n", n, amid, total);
171
172 // free(ExtremeVal);
173 // free(ExtremeCDF);
174 // free(DataVal);
175 // free(DataCDF);
176 }
177
178 return 0;
179}
DM
double DM
Definition Extreme_Compare.cpp:11
residual
double residual(double a, Integer_t const NExt, double const *ExVal, double const *ExCDF, Integer_t const Ndata, double const *DatVal, double const *DatCDF)
Definition Extreme_Compare.cpp:13
Integer_t
MKL_INT Integer_t
Definition mytypes.hpp:359

◆ residual()

double residual ( double  a,
Integer_t const  NExt,
double const *  ExVal,
double const *  ExCDF,
Integer_t const  Ndata,
double const *  DatVal,
double const *  DatCDF 
)
13 {
14 Integer_t j;
15 double p;
16 double res=0, approx, count=0;
17 for(Integer_t i = 0;i < Ndata; ++i) {
18 // if(DatVal[i] > 1.5*DM) break;
19 for(j = 1;j<NExt && DatVal[i]>a*ExVal[j]; ++j);
20 p = (DatVal[i] -a*ExVal[j-1]) / (a*ExVal[j] -a*ExVal[j-1]);
21 if(p<0 || p>1) continue;
22 approx = (1-p)*ExCDF[j-1] +p*ExCDF[j];
23 res += (log(DatCDF[i]/approx))*(log(DatCDF[i]/approx));
24 // res += (DatCDF[i] -approx)*(DatCDF[i] -approx);
25 count +=1 ;
26 }
27 return res / (double)count;
28}

Variable Documentation

◆ DM

double DM