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

Functions

int main (int argc, char **argv)
 

Variables

namespace filesystem = std::experimental::filesystem
 

Function Documentation

◆ main()

int main ( int  argc,
char **  argv 
)
55 {
56 if(argc < 8) {
57 std::cerr << "Usage: 1.This 2.Directory 3.Nmax 4.Nmin 5.dE 6.Emin 7.Emax 8.minNdata (Op1).SampleMin (Op2).SampleMax" << std::endl;
58 std::exit(EX_USAGE);
59 }
60 int const Nmax = (argc>=3 ? std::stoi(argv[2]) :16);
61 int const Nmin = (argc>=4 ? std::stoi(argv[3]) :16);
62 double const dE = (argc>=5 ? std::stod(argv[4]) :0.02);
63 double const Emin = (argc>=6 ? std::stod(argv[5]) :0.45);
64 double const Emax = (argc>=7 ? std::stod(argv[6]) :0.55);
65 int const minNdata = (argc>=8 ? std::stod(argv[7]) :2);
66 int const SampleMin = (argc>=9 ? std::stod(argv[8]) :-1);
67 int const SampleMax = (argc>=10? std::stod(argv[9]) :-1);
68 if(Nmax < Nmin) {
69 std::cerr << "Error: Nmax(" << Nmax << ") < Nmin(" << Nmin << ")" << std::endl;
70 std::exit(EX_USAGE);
71 }
72 std::cout << "# Input: Nmax=" << Nmax
73 << ", Nmin=" << Nmin
74 << ", dE=" << dE
75 << ", Emin=" << Emin
76 << ", Emax=" << Emax
77 << ", minNdata=" << minNdata
78 << std::endl;
79
80 std::string const baseDir(argv[1]);
81 std::cout << "# Directory: " << baseDir << "\n";
82 std::string InDataDir(baseDir+"/RawData");
83 std::vector<std::string> dirList;
84 {
85 std::string str;
86 int dataNo;
87 int const idShift=(SampleMin>=0? SampleMin: 0);
88 if(SampleMin>=0 && SampleMax>0) dirList.resize(SampleMax-idShift+1);
89 for(const filesystem::directory_entry &entry : filesystem::directory_iterator(InDataDir)) {
90 if ( entry.is_directory() ) {
91 str = entry.path().string();
92 dataNo = std::atoi(str.substr(str.find("_No")+3).c_str());
93 if((SampleMin>=0 || SampleMax>0) && (dataNo<SampleMin || SampleMax<dataNo)) continue;
94 debug_print("# dir=" << str << ", dataNo=" << dataNo);
95 if( (SampleMin<0 || SampleMax<=0) && dirList.size()<=dataNo) dirList.resize(dataNo+1);
96 debug_print("dataNo-idShift=" << dataNo-idShift << ", dataNo=" << dataNo << ", idShift=" << idShift);
97 dirList.at(dataNo-idShift) = str;
98 }
99 }
100 }
101
102 debug_print("# dirList.size() = " << dirList.size() << ", Nmax=" << Nmax);
103
104
105 matrix<double> L2norm(Nmax+1, dirList.size(), NAN);
106 matrix<double> L2norm_err(Nmax+1, dirList.size(), NAN);
107 matrix<double> KLdive(Nmax+1, dirList.size(), NAN);
108 matrix<double> KLdive_err(Nmax+1, dirList.size(), NAN);
109 matrix<Integer_t> dataNum(Nmax+1, dirList.size(), -1);
110 matrix<std::vector<double>> varR(Nmax+1, dirList.size());
111 std::vector<std::vector<double>> varRmin(Nmax+1, std::vector<double>(dirList.size(), NAN));
112 std::vector<std::vector<double>> varRmax(Nmax+1, std::vector<double>(dirList.size(), NAN));
113
114 #pragma omp parallel for schedule(dynamic, 10)
115 for(size_t sample = 0;sample < dirList.size(); ++sample) {
116 // #pragma omp critical
117 // std::cout << "# (for) Directory: " << dirList[sample] << "\n";
118 // for(size_t sample = 0;sample < 1; ++sample) {
119 int n, idMin, idMax;
120 Integer_t dim, Ndata;
121 std::string command;
122 double OpRange, shellMin, shellMax, dtemp;
123 std::vector<double> eigenEnergy, eigenExpValue, MCAverage, Stddev;
124 std::vector<Integer_t> NdataInShell, Histogram;
125 std::string InExpValueFname, InEnergyFname, InMatElemsFname, str;
126 std::ifstream InExpValueFs, InEnergyFs, InMatElemsFs;
127 std::string OutFname;
128 std::ofstream OutFs;
129 std::vector<std::string> candidateFileNames;
130
131
132 // No = std::stoi(entry.path().filename().string<char>().substr(9));
133 std::string OutDir(std::regex_replace(dirList[sample], std::regex("/RawData"), "/ProcessedData"));
134 filesystem::create_directories(OutDir);
135
136 bool flag = false;
137 for(n = Nmax;n >= Nmin; --n) {
138 dim = 0;
139 std::stringstream buff(""); buff << "_N" << n << ".txt";
140 { //-------------------- Load EigenExpValues --------------------//
141 std::stringstream ss;
142 candidateFileNames.resize(3);
143 candidateFileNames[0] = dirList[sample] + "/EigenExpValueZ" + buff.str();
144 candidateFileNames[1] = dirList[sample] + "/EigenExpValueC" + buff.str();
145 candidateFileNames[2] = dirList[sample] + "/EigenExpValue" + buff.str();
146 if( InExpValueFs.is_open() ) { InExpValueFs.close(); InExpValueFs.clear(); }
147 for(auto name : candidateFileNames) {
148 InExpValueFname = name;
149 InExpValueFs.open(InExpValueFname);
150 if( InExpValueFs.is_open() ) break;
151 }
152 if( !InExpValueFs.is_open() ) continue;
153 flag = true;
154 InExpValueFs.exceptions(std::ios::badbit | std::ios::failbit);
155 while( std::getline(InExpValueFs, str) ) if(str.find("OpRange") != std::string::npos) break;
156 OpRange = std::stod(str.substr(str.find("OpRange=")+8));
157
158 eigenEnergy.resize(0);
159 eigenExpValue.resize(0);
160 try{
161 skip_header(InExpValueFs);
162 for(getline(InExpValueFs, str); !InExpValueFs.eof(); getline(InExpValueFs, str)) {
163 ss.str(str); ss.clear(std::stringstream::goodbit);
164 ss >> dtemp; eigenEnergy.push_back(dtemp);
165 ss >> dtemp;
166 ss >> dtemp; eigenExpValue.push_back(dtemp/OpRange);
167 }
168 } catch(std::exception& e) {
169 if( !InExpValueFs.eof() )
170 std::cerr << "Error(" << InExpValueFs.eof() << InExpValueFs.fail() << InExpValueFs.bad() << "):" << e.what() << std::endl;
171 }
172 InExpValueFs.close(); InExpValueFs.clear();
173 if(eigenEnergy.size() != eigenExpValue.size()) {
174 std::cerr << "Error: eigenEnergy.size() != eigenExpValue.size() for a file\n"
175 << "\t" << InExpValueFname << std::endl;
176 continue;
177 }
178 dim = eigenEnergy.size();
179 }
180
181 #ifndef NDEBUG
182 std::cout << std::scientific;
183 std::cout << "# OpRange = " << OpRange << std::endl;
184 // print(eigenEnergy, eigenEnergy.size());
185 // print(eigenExpValue, eigenExpValue.size());
186 #endif
187 debug_print("########## Some process ##########");
188
189 for(idMin=0; idMin<dim && eigenEnergy.at(idMin)<Emin; ++idMin) { };
190 for(idMax=dim-1; idMax>=0 && eigenEnergy.at(idMax)>Emax; --idMax) { };
191 debug_print("# eigenEnergy[" << idMin << "]=" << eigenEnergy[idMin] << ", eigenEnergy[" << idMax << "]=" << eigenEnergy[idMax]);
192
193 MCAverage.resize(dim); NdataInShell.resize(dim);
194 std::fill( MCAverage.begin(), MCAverage.end(), 0.0);
195 std::fill(NdataInShell.begin(), NdataInShell.end(), 0.0);
196 for(Integer_t i = 0;i < dim; ++i) {
197 MCAverage[i] = MicroCanonicalAverage(NdataInShell[i], eigenEnergy[i], dE, dim, eigenEnergy, eigenExpValue, i);
198 }
199
200 debug_print("# Calculating the standard deviation of eigenstate expectation values within a shell.");
201 Stddev.resize(dim);
202 std::fill(Stddev.begin(), Stddev.end(), NAN);
203 // #pragma omp parallel for
204 for(size_t i = idMin;i <= idMax; ++i) {
205 if(NdataInShell[i] < minNdata) { Stddev[i] = NAN; continue; }
206 shellMin = eigenEnergy[i] - dE;
207 shellMax = eigenEnergy[i] + dE;
208 if(NdataInShell[i] != (Ndata = MeasureOfETH(Stddev[i], dtemp, dtemp, shellMin, shellMax, dim, eigenEnergy, eigenExpValue, MCAverage)) ) {
209 std::cerr << "Warning: NdataInShell[" << i << "](" << NdataInShell[i] << ") != Ndata(" << Ndata <<")." << std::endl;
210 // std::exit(EX_DATAERR);
211 }
212 // if( sample == 1070 ) {
213 // std::cerr << "# NdataInShell[" << i << "] = " << NdataInShell[i] << ", eigenEnergy[i]=" << eigenEnergy[i] << ", Stddev[i]=" << Stddev[i] << std::endl;
214 // }
215 }
216
217 debug_print("# Calculating variable R and its average and standard deviation.");
218 varR.at(n,sample).resize(0);
219 std::fill(varR.at(n,sample).begin(), varR.at(n,sample).end(), NAN);
220 double ave = 0, var = 0;
221 Ndata = 0;
222 for(size_t i = idMin;i <= idMax; ++i) {
223 if( isnan(Stddev[i]) ) continue;
224 dtemp = (Stddev[i] == 0 ? 0 : (eigenExpValue[i] - MCAverage[i])/Stddev[i] );
225 ave += dtemp;
226 Ndata += 1;
227 varRmin[n][sample] = std::fmin(dtemp,varRmin[n][sample]);
228 varRmax[n][sample] = std::fmax(dtemp,varRmax[n][sample]);
229 varR.at(n,sample).push_back(dtemp);
230 }
231 // if( Ndata < minNdata ) continue;
232 ave /= Ndata;
233 for(size_t i = 0;i < varR.at(n,sample).size(); ++i) {
234 if( isnan(varR.at(n,sample)[i]) ) continue;
235 var += (varR.at(n,sample)[i] - ave)*(varR.at(n,sample)[i] - ave)/(Ndata-1);
236 }
237 // for(size_t i = idMin;i <= idMax; ++i) varR[i] /= sqrt(var);
238 debug_print("# n=" << n << ", Ndata=" << Ndata << ", ave=" << ave << ", var=" << var);
239 if( !(var > 0) ) continue;
240 double binwidth = 0.1;
241 varRmin[n][sample] = std::floor(varRmin[n][sample]/binwidth)*binwidth;
242 varRmax[n][sample] = std::ceil( varRmax[n][sample]/binwidth)*binwidth;
243 Integer_t Nbin = (Integer_t)((varRmax[n][sample]-varRmin[n][sample])/binwidth + 1);
244 auto Bin = [&](double x) { Integer_t bin = (Integer_t)((x-varRmin[n][sample])/binwidth); return (bin==Nbin ? Nbin-1 : bin); };
245
246 debug_print("# varRmin=" << varRmin[n][sample] << ", varRmax=" << varRmax[n][sample] << ", binwidth=" << binwidth << ", Nbin=" << Nbin << ", Ndata=" << Ndata);
247 debug_print("# Calculate the probability density (Nbin=" << Nbin << ")");
248 Histogram.resize(Nbin);
249 std::fill(Histogram.begin(), Histogram.end(), 0);
250 for(size_t i = 0;i < varR.at(n,sample).size(); ++i) {
251 if( isnan(varR.at(n,sample)[i]) ) continue;
252 Histogram.at( Bin(varR.at(n,sample)[i]) ) += 1;
253 }
254
255 debug_print("# Calculate the L2-norm and the KL-divergence.");
256 Measure(L2norm.at(n,sample), L2norm_err.at(n,sample), KLdive.at(n,sample), KLdive_err.at(n,sample), dataNum.at(n,sample), Histogram, binwidth, varRmin[n][sample]);
257 if( L2norm.at(n,sample) > 1.0E+3 || Ndata == 2 ) {
258 std::cerr << "# (sample=" << sample << ", n=" << n << ") Nbin=" << Nbin << ", Ndata=" << dataNum.at(n,sample) << ", delta=" << L2norm.at(n,sample) << ", binwidth=" << binwidth << ", binwidth*Ndata=" << binwidth*dataNum.at(n,sample) << ", var=" << var << std::endl;
259 }
260 debug_print("########## (END)Some process ##########\n");
261
262 {
263 std::stringstream buff(""); buff << "_minNdata" << minNdata << "_N" << n << ".txt";
264 OutFname = OutDir + "/SrednickiAnsatz_varR" + buff.str();
265 }
266 if( OutFs.is_open() ) { OutFs.close(); OutFs.clear(); }
267 OutFs.open(OutFname); checkIsFileOpen(OutFs, OutFname);
268 OutFs << std::scientific << std::showpos << std::setprecision(6);
269 OutFs << "# shellWidth= " << dE << "\n"
270 << "# L2norm= " << L2norm.at(n,sample) << "\n"
271 << "# L2norm_err= " << L2norm_err.at(n,sample) << "\n"
272 << "# KLdive= " << KLdive.at(n,sample) << "\n"
273 << "# KLdive_err= " << KLdive_err.at(n,sample) << "\n"
274 << "# Ndata= " << dataNum.at(n,sample) << "\n"
275 << "# 1.(varR) 2.(PDF) 3.(Histogram)" << "\n\n";
276 for(size_t bin = 0;bin < Nbin; ++bin) {
277 OutFs << (bin+0.5)*binwidth + varRmin[n][sample] << " "
278 << Histogram[bin]/(binwidth*dataNum.at(n,sample)) << " "
279 << Histogram[bin]
280 << std::endl;
281 }
282 OutFs.close();
283
284
285 {
286 std::stringstream buff(""); buff << "_minNdata" << minNdata << "_N" << n << ".txt";
287 OutFname = OutDir + "/MCAverage" + buff.str();
288 }
289 if( OutFs.is_open() ) { OutFs.close(); OutFs.clear(); }
290 OutFs.open(OutFname); checkIsFileOpen(OutFs, OutFname);
291 OutFs << std::scientific << std::showpos << std::setprecision(6);
292 OutFs << "# shellWidth= " << dE << "\n"
293 << "# 1.(eigenEnergy) 2.(eigenExpValue) 3.(MCAverage) 4.(Stddev)" << "\n\n";
294 for(Integer_t i = 0;i < dim; ++i) {
295 OutFs << eigenEnergy[i] << " "
296 << eigenExpValue[i] << " "
297 << MCAverage[i] << " "
298 << Stddev[i] << std::endl;
299 }
300 OutFs.close();
301 }
302
303 if(flag) {
304 #pragma omp critical
305 std::cout << "# " << dirList[sample] << std::endl;
306 }
307 }
308
309
310 debug_print("# Calculate the statistics.");
311 std::string OutDir(baseDir);
312 {
313 std::stringstream buff("");
314 buff << "/Measure_Emin" << Emin << "_Emax" << Emax << "_dE" << dE;
315 OutDir = baseDir+buff.str();
316 filesystem::create_directories(OutDir);
317 filesystem::create_directories(OutDir+"/Output");
318 }
319 {
320 std::string OutStatFname;
321 std::stringstream buff(""); buff << "_minNdata" << minNdata << ".txt";
322 OutStatFname = OutDir + "/Output/SrednickiAnsatz_varR_L2norm" + buff.str();
323 std::ofstream OutStat(OutStatFname); checkIsFileOpen(OutStat, OutStatFname);
324 OutStat << "# shellWidth= " << dE << ", minNdata=" << minNdata << "\n"
325 << "# 1.(System size) 2.(average L2norm) 3.(L2norm_err) 4.(stddev of L2norm) 5.(NdataL2) 6.(average KLdive) 7.(KLdive_err) 8.(stddev of KLdive) 9.(NdataKL)" << "\n\n";
326 OutStat << std::scientific << std::setprecision(6);
327 std::string OutFname;
328 for(int n = Nmax;n >= Nmin; --n) {
329 int NdataL2=0, NdataKL=0;
330 double meanL2=0, stddevL2=0, meanKL=0, stddevKL=0;
331 double meanL2_err=0, meanKL_err=0;
332 for(size_t sample = 0;sample < dirList.size(); ++sample) {
333 if( !isnan(L2norm.at(n,sample)) ) { meanL2 += L2norm.at(n,sample); NdataL2 += 1; }
334 if( !isnan(KLdive.at(n,sample)) ) { meanKL += KLdive.at(n,sample); NdataKL += 1; }
335 }
336 meanL2 /= NdataL2; meanKL /= NdataKL;
337 for(size_t sample = 0;sample < dirList.size(); ++sample) {
338 if( !isnan(L2norm.at(n,sample)) ) { stddevL2 += (meanL2 - L2norm.at(n,sample))*(meanL2 - L2norm.at(n,sample)); }
339 if( !isnan(KLdive.at(n,sample)) ) { stddevKL += (meanKL - KLdive.at(n,sample))*(meanKL - KLdive.at(n,sample)); }
340 if( !isnan(L2norm_err.at(n,sample)) ) { meanL2_err += L2norm_err.at(n,sample)*L2norm_err.at(n,sample) / NdataL2; }
341 if( !isnan(KLdive_err.at(n,sample)) ) { meanKL_err += KLdive_err.at(n,sample)*KLdive_err.at(n,sample) / NdataKL; }
342 }
343 stddevL2 = sqrt(stddevL2/(NdataL2-1)); stddevKL = sqrt(stddevKL/(NdataKL-1));
344 meanL2_err = sqrt(meanL2_err+stddevL2*stddevL2)/sqrt(NdataL2);
345 meanKL_err = sqrt(meanKL_err+stddevKL*stddevKL)/sqrt(NdataKL);
346 // if( !(Ndata > 0) ) continue;
347 OutStat << n << " "
348 << meanL2 << " " << meanL2_err << " " << stddevL2 << " " << NdataL2 << " "
349 << meanKL << " " << meanKL_err << " " << stddevKL << " " << NdataKL << std::endl;
350
351 {
352 std::stringstream buff(""); buff << "_minNdata" << minNdata << "_N" << n << ".txt";
353 OutFname = OutDir + "/SrednickiAnsatz_varR_L2norm" + buff.str();
354 }
355 std::ofstream OutFs(OutFname); checkIsFileOpen(OutFs, OutFname);
356 OutFs << "# shellWidth= " << dE << ", minNdata=" << minNdata << "\n"
357 << "# 1.(sample) 2.(L2norm) 3.(KLdive) 4.(dataNum)" << "\n\n";
358 OutFs << std::scientific << std::setprecision(6);
359 for(size_t sample = 0;sample < dirList.size(); ++sample) {
360 OutFs << sample << " "
361 << L2norm.at(n,sample) << " "
362 << KLdive.at(n,sample) << " "
363 << dataNum.at(n,sample) << "\n";
364 }
365 }
366 }
367
368 debug_print("# Calculate the ensemble average.");
369 {
370 double binwidth = 0.1;
371 std::string OutFname;
372 for(int n = Nmax;n >= Nmin; --n) {
373 debug_print("# Calculate the ensemble average. (N=" << n << ")");
374 // print(varRmin[n], varRmin[n].size());
375 // print(varRmax[n], varRmax[n].size());
376 double VarRMin = *std::min_element(varRmin[n].begin(), varRmin[n].end());
377 double VarRMax = *std::max_element(varRmax[n].begin(), varRmax[n].end());
378
379 if( isnan(VarRMin) || isnan(VarRMax) ) continue;
380 Integer_t Nbin = (Integer_t)((VarRMax-VarRMin)/binwidth + 1);
381 auto Bin = [&](double x) { Integer_t bin = (Integer_t)((x-VarRMin)/binwidth); return (bin==Nbin ? Nbin-1 : bin); };
382
383 std::cout << "VarRMin = " << VarRMin << ", VarRmax = " << VarRMax << ", Nbin = " << Nbin << std::endl;
384 debug_print("# Calculate the ensemble average.");
385 Integer_t Ndata = 0;
386 std::vector<Integer_t> Histogram(Nbin,0);
387 for(size_t sample = 0;sample < dirList.size(); ++sample) {
388 for(size_t j = 0;j < varR.at(n,sample).size(); ++j) {
389 if( isnan(varR.at(n,sample)[j]) ) continue;
390
391 if( Bin(varR.at(n,sample)[j]) >= Nbin ) {
392 std::cerr << "Error: Bin(varR.at(" << n << "," << sample << ")[" << j << "]) = " << Bin(varR.at(n,sample)[j]) << " >= " << Nbin << std::endl;
393 }
394 Histogram.at( Bin(varR.at(n,sample)[j]) ) += 1;
395 Ndata += 1;
396 }
397 }
398
399 debug_print("# Calculate the L2-norm and the KL-divergence.");
400 double L2=0, KL=0, L2_err=0, KL_err=0;
401 Measure(L2, L2_err, KL, KL_err, Ndata, Histogram, binwidth, VarRMin);
402
403 {
404 std::stringstream buff(""); buff << "_minNdata" << minNdata << "_N" << n << ".txt";
405 OutFname = OutDir + "/SrednickiAnsatz_varR_AverageHist" + buff.str();
406 }
407 std::ofstream OutFs(OutFname); checkIsFileOpen(OutFs, OutFname);
408 OutFs << std::scientific << std::showpos << std::setprecision(6);
409 OutFs << "# shellWidth= " << dE << "\n"
410 << "# L2norm= " << L2 << "\n"
411 << "# L2norm_err= " << L2_err<< "\n"
412 << "# KLdive= " << KL << "\n"
413 << "# KLdive_err= " << KL_err<< "\n"
414 << "# Ndata= " << Ndata << "\n"
415 << "# 1.(varR) 2.(PDF) 3.(Histogram)" << "\n\n";
416 for(size_t bin = 0;bin < Histogram.size(); ++bin) {
417 OutFs << (bin+0.5)*binwidth + VarRMin << " "
418 << Histogram[bin]/(binwidth*Ndata) << " "
419 << Histogram[bin]
420 << std::endl;
421 }
422 OutFs.close();
423 }
424 }
425
426 return EXIT_SUCCESS;
427}
Nbin
Integer_t const Nbin
Definition OpErgodicity_Compare.cpp:13
Bin
Integer_t Bin(double x)
Definition Overlap.cpp:59
filesystem
namespace filesystem
Definition SrednickiAnsatz_varR.cpp:18
Histogram
Definition OpAverage.cpp:21
MicroCanonicalAverage
Calculate the microcanonical averages with respect to a given sorted vector 'eigVal'.
Definition MicroCanonicalAverage.hpp:25
checkIsFileOpen
bool checkIsFileOpen(std::ifstream &file, std::string const &filename)
Definition file_util.hpp:22
debug_print
debug_print("# Determining GPU configuration.")
Integer_t
MKL_INT Integer_t
Definition mytypes.hpp:359
Emin
constexpr double Emin
Definition setVariablesForMCAverage.cpp:4
Emax
constexpr double Emax
Definition setVariablesForMCAverage.cpp:5
dE
double const dE
Definition setVariablesForMCAverage.cpp:2
MeasureOfETH
Integer_t MeasureOfETH(double &Fluc_2, double &Fluc_inf, double &MaxPos, double const Emin, double const Emax, Integer_t const dim, const std::vector< double > &eigenEnergy, const std::vector< double > &EXPvalue, const std::vector< double > &MCAverage)
Definition statmech.cpp:244
size
__constant__ int size
Definition testEigenOnGPU.cu:4
buff
std::stringstream buff("")

Variable Documentation

◆ filesystem

namespace filesystem = std::experimental::filesystem