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

Functions

double HistogramIntersection (double Dist1[NBIN], Integer_t Dist1Ndata[NBIN], double Dist2[NBIN])
 
double KLDivergence (double Dist1[NBIN], Integer_t Dist1Ndata[NBIN], double Dist2[NBIN])
 
double JSDivergence (double Dist1[NBIN], Integer_t Dist1Ndata[NBIN], double Dist2[NBIN], Integer_t Dist2Ndata[NBIN])
 
double L1Norm (double Dist1[NBIN], Integer_t Dist1Ndata[NBIN], double Dist2[NBIN])
 
double L2Norm (double Dist1[NBIN], Integer_t Dist1Ndata[NBIN], double Dist2[NBIN])
 
int main (int argc, char **argv)
 

Variables

Integer_t const NBIN = 100
 
Integer_t const NMAX = 20
 
Integer_t const Nitem = 4
 

Function Documentation

◆ HistogramIntersection()

double HistogramIntersection ( double  Dist1[NBIN],
Integer_t  Dist1Ndata[NBIN],
double  Dist2[NBIN] 
)
16 {
17 double count = 0;
18 double res = 0;
19 for(Integer_t bin = 0;bin < NBIN; ++bin) {
20 if( Dist1Ndata[bin] == 0 ) continue;
21 if( Dist1[bin] > Dist2[bin] ) res += Dist2[bin];
22 else res += Dist1[bin];
23 count += Dist1[bin];
24 // if( res ) printf("KLDivergence: Dist1[%lld]=%f, Dist2[%lld]=%f\n", i, Dist1[i], i, Dist2[i]);
25 }
26 return res/count;
27}
NBIN
Integer_t const NBIN
Definition OpErgodicity.cpp:12
Integer_t
MKL_INT Integer_t
Definition mytypes.hpp:359

◆ JSDivergence()

double JSDivergence ( double  Dist1[NBIN],
Integer_t  Dist1Ndata[NBIN],
double  Dist2[NBIN],
Integer_t  Dist2Ndata[NBIN] 
)
41 {
42 double res = KLDivergence(Dist1, Dist1Ndata, Dist2);
43 res += KLDivergence(Dist2, Dist2Ndata, Dist1);
44 return res /= 2.0;
45}
KLDivergence
double KLDivergence(double Dist1[NBIN], Integer_t Dist1Ndata[NBIN], double Dist2[NBIN])
Definition OpErgodicity.cpp:29

◆ KLDivergence()

double KLDivergence ( double  Dist1[NBIN],
Integer_t  Dist1Ndata[NBIN],
double  Dist2[NBIN] 
)
29 {
30 double count = 0;
31 double res = 0;
32 for(Integer_t bin = 0;bin < NBIN; ++bin) {
33 if( Dist1[bin] == 0 || Dist2[bin] == 0 ) continue;
34 res += Dist1[bin]*log(Dist1[bin]/Dist2[bin]);
35 // if( isnan(res) ) printf("KLDivergence: Dist1[%lld]=%f, Dist2[%lld]=%f\n", i, Dist1[i], i, Dist2[i]);
36 count += Dist1[bin];
37 }
38 return res/(double)count;
39}

◆ L1Norm()

double L1Norm ( double  Dist1[NBIN],
Integer_t  Dist1Ndata[NBIN],
double  Dist2[NBIN] 
)
47 {
48 double count = 0;
49 double res = 0;
50 for(Integer_t bin = 0;bin < NBIN; ++bin) {
51 // if( Dist1Ndata[bin] == 0 ) continue;
52 res += fabs(Dist1[bin] -Dist2[bin]);
53 // if( isnan(res) ) printf("L1Norm: Dist1[%lld]=%f, Dist2[%lld]=%f\n", i, Dist1[i], i, Dist2[i]);
54 count += 1;
55 }
56 count = NBIN; // binwidthの逆数(Ndata=0のbinもデータとして数える)
57 // printf("L1Norm: res=%f, count=%f\n", res, count);
58 return res/count;
59}

◆ L2Norm()

double L2Norm ( double  Dist1[NBIN],
Integer_t  Dist1Ndata[NBIN],
double  Dist2[NBIN] 
)
61 {
62 Integer_t count = 0;
63 double res = 0;
64 for(Integer_t bin = 0;bin < NBIN; ++bin) {
65 if( Dist1Ndata[bin] == 0 ) Dist1[bin] = 0;
66 if( isnan(Dist2[bin]) ) { count += 1; continue; }
67 res += (Dist1[bin] -Dist2[bin])*(Dist1[bin] -Dist2[bin]);
68 if( isnan(res) ) printf("L2Norm: Dist1[%lld]=%f, Dist2[%lld]=%f\n", bin, Dist1[bin], bin, Dist2[bin]);
69 // count += 1;
70 // printf("count=%lld, Dist1Ndata[%lld]=%lu\n", count, bin, Dist1Ndata[bin]);
71 }
72 count = NBIN -count; // binwidthの逆数(Ndata=0のbinもデータとして数える)
73 if( isinf(res) || count == 0 ) {res = NAN; count = -1;}
74 if( isnan(res) ) printf("L2Norm: res=%f, count=%lld\n", res, count);
75 return sqrt(res/(double)count);
76}

◆ main()

int main ( int  argc,
char **  argv 
)
78 {
79 if(argc < 2) {
80 printf("Usage: 1.This 2.Directory\n");
81 exit(EXIT_FAILURE);
82 }
83 Integer_t i,j;
84 Integer_t Nsample = 0, counter = 0;
85 Integer_t n, bin, itemp;
86 Integer_t Ntotal[NMAX];
87 Integer_t NfuncF[NBIN];
88 double ratio, expval, stddev, fE, dos;
89
90 struct stat statBuf;
91 struct Distance {
92 double HI;
93 double KLdivergence;
94 double JSdivergence;
95 double L1norm;
96 double L2norm;
97 };
98 struct Distance distance[NMAX][Nitem]; // 分布間の距離の指標
99
100 struct EnsembleAverage {
101 Integer_t N[NBIN];
102 double ExpValue[NBIN];
103 double Stddev[NBIN];
104 double funcF[NBIN];
105 double DOS[NBIN];
106 };
107 struct EnsembleAverage Average[NMAX];
108 for(i = 0;i < NMAX; ++i) for(j = 0;j < NBIN; ++j) {
109 Average[i].N[j] = 0;
110 Average[i].ExpValue[j] = 0;
111 Average[i].Stddev[j] = 0;
112 Average[i].funcF[j] = 0;
113 Average[i].DOS[j] = 0;
114 }
115
116 struct Bin {
117 Integer_t N[NBIN];
118 Integer_t Nnan[NBIN];
119 double ExpValue[NBIN];
120 double Stddev[NBIN];
121 double funcF[NBIN];
122 double DOS[NBIN];
123 };
124 struct Bin BinData[NMAX];
125
126 Integer_t const BUFSIZE = 1000;
127 char dstr[BUFSIZE], str[BUFSIZE], stemp[BUFSIZE], outstr[BUFSIZE];
128 char *filename = str;
129 char *stemp_p;
130 char *ext;
131 strcpy(dstr, argv[1]);
132 strcat(dstr, "/Output_New");
133 DIR *dir;
134 FILE *Input, *Output;
135 struct dirent *dp;
136
137 //********** アンサンブル平均の読み込み ********** //
138 strcpy(str, dstr);
139 strcat(str, "/EnsembleAverage.txt");
140 printf(" Input file: %s\n", filename);
141 Input = fopen(filename, "r");
142 if(Input == NULL) {
143 printf("Can't open a file %s\n", filename);
144 exit(EXIT_FAILURE);
145 }
146 fgets(stemp, BUFSIZE-1, Input);
147 if( !strstr(stemp, "NBIN") ) {
148 fprintf(stderr, "Error: Can't find \"NBIN\" in input file \"%s\".\n", filename);
149 exit(EXIT_FAILURE);
150 }
151 stemp_p = strrchr(stemp, ' ') +1;
152 if( atoi(stemp_p) != NBIN ) {
153 fprintf(stderr, "Error: NBIN(%lld) != %lld.\n", atoi(stemp_p), NBIN);
154 exit(EXIT_FAILURE);
155 }
156 skip_header(Input);
157 while( fscanf(Input, "%lld %lf %le %le %le %le %lld", &n, &ratio, &expval, &stddev, &fE, &dos, &itemp) != EOF ) {
158 if(ratio == 1) bin = NBIN-1; else bin = (Integer_t)(ratio*NBIN);
159 Average[n].N[bin] = itemp;
160 Average[n].ExpValue[bin] = expval;
161 Average[n].Stddev[bin] = stddev;
162 Average[n].funcF[bin] = fE;
163 Average[n].DOS[bin] = dos;
164 }
165 fclose(Input);
166 //********** (END)アンサンブル平均の読み込み ********** //
167
168 //********** Output fileを作成 ********** //
169 strcpy(str, dstr);
170 strcat(str, "/");
171 strcat(str, "/ErgodicityFluctuation.txt");
172 strcpy(stemp, str);
173 printf("Output file: %s\n", filename);
174 Output = fopen(filename, "w");
175 if(Output == NULL) {
176 printf("Can't open a file %s\n", filename);
177 exit(EXIT_FAILURE);
178 }
179 fprintf(Output, "# NBIN = %lld\n", NBIN);
180 fprintf(Output, "# 1.Nsample 2.N 3.4.5.6.7ExpValue 8.9.10.11.12.Stddev 13.14.15.16.17 funcF 18.19.20.21.22.DOS\n\n");
181 strcpy(outstr, filename);
182 //********** (END)Output fileを作成 ********** //
183
184 strcpy(str, dstr); strcat(str, "/ProcessedData");
185 dir = opendir(str);
186 for(dp = readdir(dir);dp != NULL;dp = readdir(dir)) {
187 if( strstr(dp->d_name, "Operators_Hist") == NULL ) continue;
188 ext = strrchr(dp->d_name, '.');
189 if( strcmp(ext, ".txt") != 0 ) continue;
190
191 strcpy(str, dstr); strcat(str, "/ProcessedData");
192 strcat(str, "/");
193 strcat(str, dp->d_name);
194 if( stat(filename, &statBuf) == -1 ){
195 printf("Can't get file information: %s.", filename);
196 continue;
197 }
198 counter += 1;
199 printf("(%4d)Input file: %s, size: %lld\n", counter, filename, (long long int)statBuf.st_size);
200 if( statBuf.st_size==0 ) continue;
201 Input = fopen(filename, "r");
202 if(Input == NULL) {
203 printf("Can't open a file %s\n", filename);
204 continue;
205 }
206 fgets(stemp, BUFSIZE-1, Input);
207 if( !strstr(stemp, "NBIN") ) {
208 fprintf(stderr, "Error: Can't find \"NBIN\" in input file \"%s\".\n", filename);
209 exit(EXIT_FAILURE);
210 }
211 stemp_p = strrchr(stemp, ' ') +1;
212 if( atoi(stemp_p) != NBIN ) {
213 fprintf(stderr, "Error: NBIN(%lld) != %lld.\n", atoi(stemp_p), NBIN);
214 exit(EXIT_FAILURE);
215 }
216 skip_header(Input);
217
218 //********** 変数の初期化 **********//
219 for(i = 0;i < NMAX; ++i) Ntotal[i] = 0;
220 for(i = 0;i < NMAX; ++i) for(j = 0;j < Nitem; ++j) {
221 distance[i][j].HI = 0;
222 distance[i][j].KLdivergence = 0;
223 distance[i][j].JSdivergence = 0;
224 distance[i][j].L1norm = 0;
225 distance[i][j].L2norm = 0;
226 }
227 for(i = 0;i < NMAX; ++i) for(j = 0;j < NBIN; ++j) {
228 BinData[i].N[j] = 0;
229 BinData[i].Nnan[j] = 0;
230 BinData[i].ExpValue[j] = 0;
231 BinData[i].Stddev[j] = 0;
232 BinData[i].funcF[j] = 0;
233 BinData[i].DOS[j] = 0;
234 }
235 //********** (END)変数の初期化 **********//
236
237 Nsample += 1;
238 while( fscanf(Input, "%lld %le %le %le %le %le %lld %lld", &n, &ratio, &expval, &stddev, &fE, &dos, &itemp, &itemp) != EOF ) {
239 Ntotal[n] +=1;
240 if(ratio == 1) bin = NBIN-1; else bin = (Integer_t)(ratio*NBIN);
241 BinData[n].N[bin] +=1;
242
243 BinData[n].ExpValue[bin] += expval;
244 if( isnan(stddev) || isnan(fE) ) BinData[n].Nnan[bin] +=1;
245 else {
246 BinData[n].Stddev[bin] += stddev;
247 BinData[n].funcF[bin] += fE;
248 }
249 BinData[n].DOS[bin] += dos;
250 }
251
252 //********** 関数の間の距離を計算 **********//
253 for(n = 0;n <= NMAX; ++n) {
254 if( Ntotal[n] == 0 ) continue;
255 // ExpValue
256 distance[n][0].HI = HistogramIntersection(BinData[n].ExpValue, BinData[n].N, Average[n].ExpValue);
257 distance[n][0].KLdivergence = KLDivergence(BinData[n].ExpValue, BinData[n].N, Average[n].ExpValue);
258 distance[n][0].JSdivergence = JSDivergence(BinData[n].ExpValue, BinData[n].N, Average[n].ExpValue, Average[n].N);
259 distance[n][0].L1norm = L1Norm(BinData[n].ExpValue, BinData[n].N, Average[n].ExpValue);
260 distance[n][0].L2norm = L2Norm(BinData[n].ExpValue, BinData[n].N, Average[n].ExpValue);
261
262 j = 0;
263 for(bin = 0;bin < NBIN; ++bin) {
264 NfuncF[bin] = BinData[n].N[bin] -BinData[n].Nnan[bin];
265 j += NfuncF[bin];
266 }
267 // printf("%lld, NfuncF[%lld]=%lld\n", n, bin, j);
268 // Stddev
269 distance[n][1].HI = HistogramIntersection(BinData[n].Stddev, NfuncF, Average[n].Stddev);
270 distance[n][1].KLdivergence = KLDivergence(BinData[n].Stddev, NfuncF, Average[n].Stddev);
271 distance[n][1].JSdivergence = JSDivergence(BinData[n].Stddev, NfuncF, Average[n].Stddev, Average[n].N);
272 distance[n][1].L1norm = L1Norm(BinData[n].Stddev, NfuncF, Average[n].Stddev);
273 distance[n][1].L2norm = L2Norm(BinData[n].Stddev, NfuncF, Average[n].Stddev);
274 // funcF
275 distance[n][2].HI = HistogramIntersection(BinData[n].funcF, NfuncF, Average[n].funcF);
276 distance[n][2].KLdivergence = KLDivergence(BinData[n].funcF, NfuncF, Average[n].funcF);
277 distance[n][2].JSdivergence = JSDivergence(BinData[n].funcF, NfuncF, Average[n].funcF, Average[n].N);
278 distance[n][2].L1norm = L1Norm(BinData[n].funcF, NfuncF, Average[n].funcF);
279 distance[n][2].L2norm = L2Norm(BinData[n].funcF, NfuncF, Average[n].funcF);
280 // DOS
281 distance[n][3].HI = HistogramIntersection(BinData[n].DOS, BinData[n].N, Average[n].DOS);
282 distance[n][3].KLdivergence = KLDivergence(BinData[n].DOS, BinData[n].N, Average[n].DOS);
283 distance[n][3].JSdivergence = JSDivergence(BinData[n].DOS, BinData[n].N, Average[n].DOS, Average[n].N);
284 distance[n][3].L1norm = L1Norm(BinData[n].DOS, BinData[n].N, Average[n].DOS);
285 distance[n][3].L2norm = L2Norm(BinData[n].DOS, BinData[n].N, Average[n].DOS);
286
287 fprintf(Output, "%lld %lld ", Nsample, n);
288 for(bin = 0;bin < Nitem; ++bin) {
289 fprintf(Output, "%le ", distance[n][bin].HI);
290 fprintf(Output, "%le ", distance[n][bin].KLdivergence);
291 fprintf(Output, "%le ", distance[n][bin].JSdivergence);
292 fprintf(Output, "%le ", distance[n][bin].L1norm);
293 fprintf(Output, "%le ", distance[n][bin].L2norm);
294 }
295 fprintf(Output, "\n");
296 }
297 fclose(Input);
298 }
299
300 printf("Output file: %s\n", outstr);
301 fclose(Output);
302
303 return 0;
304}
NBIN
Integer_t const NBIN
Definition OpAverage.cpp:20
HistogramIntersection
double HistogramIntersection(double Dist1[NBIN], Integer_t Dist1Ndata[NBIN], double Dist2[NBIN])
Definition OpErgodicity.cpp:16
L2Norm
double L2Norm(double Dist1[NBIN], Integer_t Dist1Ndata[NBIN], double Dist2[NBIN])
Definition OpErgodicity.cpp:61
Nitem
Integer_t const Nitem
Definition OpErgodicity.cpp:14
JSDivergence
double JSDivergence(double Dist1[NBIN], Integer_t Dist1Ndata[NBIN], double Dist2[NBIN], Integer_t Dist2Ndata[NBIN])
Definition OpErgodicity.cpp:41
NMAX
Integer_t const NMAX
Definition OpErgodicity.cpp:13
L1Norm
double L1Norm(double Dist1[NBIN], Integer_t Dist1Ndata[NBIN], double Dist2[NBIN])
Definition OpErgodicity.cpp:47
Bin
Integer_t Bin(double x)
Definition Overlap.cpp:59
Distance
Definition OpHistogram.cpp:24
Distance::L1norm
double L1norm
Definition OpHistogram.cpp:29
Distance::L2norm
double L2norm
Definition OpHistogram.cpp:30

Variable Documentation

◆ NBIN

Integer_t const NBIN = 100

◆ Nitem

Integer_t const Nitem = 4

◆ NMAX

Integer_t const NMAX = 20