OpHistogram.cpp File Reference

Classes
class	Distance

Functions
double	NormalDist (double x)
int	main (int argc, char **argv)

Variables
namespace	filesystem = std::experimental::filesystem

Function Documentation

main()

int main	(	int	argc,
		char **	argv
	)

                               {
  if(argc < 3) {
    printf("Usage: 1.This 2.Directory 3.Shell width\n");
    exit(EXIT_FAILURE);
  }
  double const cutoff   = 0;
  double const binratio = 0.05;
  std::string const d_str(argv[1]);
  double const shellWidth = (double)std::atof(argv[2]);
  std::cout << "Directory: " << d_str << "\n";
 
  Integer_t n, ndata, Nmin, Nmax, bin, binMax, count=0, Nline;
  // Integer_t ndata, itemp;
  double relEnergy, expval, binwidth, ave, aveOld, var;
  double Emin, Emax, StddevInShell, minFluc, maxFluc, Fluc, dtemp;
  std::vector<Integer_t>   Ndata, Nnan, NdataInShell;
  std::vector<double> MCAverage, NormalizedFluc, histogram, ProbDensity, gauss;
  std::vector<std::vector<double>> eigenEnergy, ExpValue;
  Distance dist;
  std::string line;
 
 
  //********** Output fileを作成 ********** //
  std::string OutDIR(d_str); OutDIR += "/Output";
  filesystem::create_directory(OutDIR);
  std::string OutDevFName(OutDIR); OutDevFName += "/EF_deviation.txt";
  std::ofstream OutDevfs(OutDevFName);
  OutDevfs.exceptions(std::ios::badbit | std::ios::failbit);
    if( !OutDevfs.is_open() ) {
      std::cerr << "Can't open a file " << OutDevFName << std::endl;
      std::exit(EX_CANTCREAT);
    }
  OutDevfs << "# shellWidth = " << shellWidth << ", binratio = " << binratio << "\n"
           << "# 1.(Sample No.) 2.(System Size) 3.(Histogram Intersection) 4.(KLdivergence) 5.(JSdivergence) 6.(L1norm) 7.(L2norm)\n"
           << std::endl;
  std::cout << "Output file: " << OutDevFName << std::endl;
 
  OutDIR = d_str; OutDIR += "/ProcessedData";
  filesystem::create_directory(OutDIR);
  //********** (END)Output fileを作成 ********** //
 
  std::string InDataDir(d_str); InDataDir += "/RawData";
  for(const filesystem::directory_entry &data : filesystem::directory_iterator(InDataDir)) {
    if( data.path().string().find("Sample_No") == std::string::npos ) continue;
    if( data.path().extension() != ".txt" ) continue;
    std::ifstream Infs(data.path());
    Infs.exceptions(std::ios::badbit | std::ios::failbit);
      if( !Infs.is_open() ) {
        std::cerr << "Can't open a file \"" << data.path() << "\"" << std::endl;
        continue;
      }
    try { skip_header(Infs); }
    catch(std::exception& e) {
      std::cerr << "Warning(skip_header): " << e.what() << std::endl;
      std::cerr << "File: \"" << data.path() << "\" is empty." << std::endl;
      continue;
    }
 
    //********** Create output file for Histogram of eigenstate fluctuation **********//
    std::string OutHistFName(OutDIR);
    OutHistFName += "/EFHistNormalized";
    OutHistFName += data.path().string().substr( data.path().string().find("_No") );
    std::ofstream OutHistfs(OutHistFName);
    OutHistfs.exceptions(std::ios::badbit | std::ios::failbit);
      if( !OutHistfs.is_open() ) {
        std::cerr << "Can't open a file " << OutHistFName << std::endl;
        std::exit(EX_CANTCREAT);
      }
    OutHistfs << "# shellWidth = " << shellWidth << ", binratio = " << binratio << "\n"
              << "# 1.(System Size) 2.(Normalized fluctuation) 3.(Probability density) 4.(data# in bin) 5.(Gauss)\n"
              << std::scientific
              << std::endl;
    std::cout << "OutEFHist: " << OutHistFName << std::endl;
    //********** (END) Create output file for Histogram of eigenstate fluctuation **********//
 
    //********** Read data **********//
    ++count;
    Nmin=INT_MAX; Nmax=INT_MIN;
    std::fill(Ndata.begin(),Ndata.end(),0);
    std::fill( Nnan.begin(), Nnan.end(),0);
    for(n = 0;n < eigenEnergy.size(); ++n) eigenEnergy.at(n).resize(0);
    for(n = 0;n <    ExpValue.size(); ++n)    ExpValue.at(n).resize(0);
    try{
      Nline = 0;
      for(std::getline(Infs, line); !Infs.eof(); std::getline(Infs, line)) {
        ++Nline;
        std::cout << "(count=" << count << ", line=" << Nline << "): " << line << ": " << Infs.eof() << "\r";
        std::istringstream inStream(line);
        inStream >> n >> dtemp >> relEnergy >> expval;
        if( isinf(expval) || isnan(expval) ) std::cerr << expval << std::endl;
        if(n < Nmin) Nmin = n;
        if(n > Nmax) Nmax = n;
        if(n >=       Ndata.size())       Ndata.resize(n+1);
        if(n >= eigenEnergy.size()) eigenEnergy.resize(n+1);
        if(n >=    ExpValue.size())    ExpValue.resize(n+1);
        Ndata.at(n) += 1;
        eigenEnergy.at(n).push_back(relEnergy);
           ExpValue.at(n).push_back(expval);
      }
    } catch(std::exception& e) {
      if( !Infs.eof() ) {
        std::cout << std::endl;
        std::cerr << "Error(" << Infs.eof() << Infs.fail() << Infs.bad() << "): " << e.what() << std::endl;
      }
    }
    Infs.close();
    //********** (END) Read data **********//
 
    Nnan.resize(Nmax+1);   std::fill(          Nnan.begin(),           Nnan.end(), 0);
    for(n = Nmax;n >= Nmin; --n){
      if(Ndata[n] == 0) continue;
      if(Ndata[n] != eigenEnergy[n].size()) {
        std::cerr << "Error: Ndata["     << n << "](" << Ndata[n]
                  << ") != eigenEnergy[" << n << "].size (" <<  eigenEnergy[n].size() << ")"
                  << std::endl;
        std::exit(EX_DATAERR);
      }
      if(Ndata[n] != ExpValue[n].size()) {
        std::cerr << "Error: Ndata["  << n << "](" << Ndata[n]
                  << ") != ExpValue[" << n << "].size (" <<  ExpValue[n].size() << ")"
                  << std::endl;
        std::exit(EX_DATAERR);
      }
 
      MCAverage.resize(Ndata[n]); NdataInShell.resize(Ndata[n]);
      std::fill(   MCAverage.begin(),    MCAverage.end(), 0.0);
      std::fill(NdataInShell.begin(), NdataInShell.end(), 0.0);
      for(Integer_t i = 0;i < Ndata[n]; ++i) {
        MCAverage[i] = MicroCanonicalAverage(NdataInShell[i], eigenEnergy[n][i], shellWidth, Ndata[n], eigenEnergy[n], ExpValue[n], i);
      }
 
      NormalizedFluc.resize(0);
      minFluc = NAN; maxFluc = NAN;
      for(int i = 0;i<Ndata[n] && eigenEnergy[n][i]<=1-cutoff; ++i) {
        if(eigenEnergy[n][i] < cutoff) continue;
        Emin = eigenEnergy[n][i] -shellWidth;
        Emax = eigenEnergy[n][i] +shellWidth;
        if((bin=MeasureOfETH(StddevInShell, dtemp, dtemp, Emin, Emax, Ndata[n], eigenEnergy[n], ExpValue[n], MCAverage)) != NdataInShell[i]) {
          std::cerr << "Warning:(n,Nmax,i)=(" << n << "," << Nmax << "," << i << ") MeasureOfETH(" << bin << ") != NdataInShell[" << NdataInShell[i] << "]" << std::endl;
        }
        if( isinf(StddevInShell) ) StddevInShell = NAN;
        if( isnan(StddevInShell) ) { Nnan[n] += 1; continue; }
        dtemp = ( ExpValue[n][i] -MCAverage[i] ) / StddevInShell;
        if(bin >= 1 && !isnan(dtemp)) {
          NormalizedFluc.push_back(dtemp);
        }
        minFluc = std::fmin(minFluc, NormalizedFluc.back());
        maxFluc = std::fmax(maxFluc, NormalizedFluc.back());
      }
      if(minFluc > maxFluc) continue;
      if(NormalizedFluc.size() == 0) continue;
 
      ndata = 0; ave = 0.0; var = 0.0;
      for(int i = 0;i < NormalizedFluc.size(); ++i) {
        if( isnan(NormalizedFluc[i]) ) continue;
        ++ndata;
        dtemp = 1/(double)ndata;
        aveOld = ave;
        ave = (1-dtemp)*aveOld + dtemp*NormalizedFluc[i];
        var += (NormalizedFluc[i]-ave)*(NormalizedFluc[i]-aveOld);
      }
      var /= (double)(ndata-1);
      if(ndata == 0) std::cerr << "ndata == 0, NormalizedFluc.size(" << NormalizedFluc.size() << ")" << std::endl;
 
      // ave = std::reduce(NormalizedFluc.begin(), NormalizedFluc.end(), 0.0) / NormalizedFluc.size();
      // var = std::accumulate(NormalizedFluc.begin(), NormalizedFluc.end(), 0.0, [ave](double sum, const auto& e){
      //           const auto temp = e - ave;
      //           return sum + temp * temp;
      //         }) / NormalizedFluc.size();
 
      //********** Create histogram ********** //
      if( var == 0 ) continue;
      binwidth = binratio *sqrt(var);
      binMax   = (Integer_t)((maxFluc-minFluc)/binwidth) +1;
      if(binMax <= 0) {
        std::cerr << "      # Warning(n=" << n << "): binMax(" << binMax << ") is not positive." << std::endl;
        continue;
       }
        histogram.resize(binMax); std::fill(  histogram.begin(),   histogram.end(), 0);
      ProbDensity.resize(binMax); std::fill(ProbDensity.begin(), ProbDensity.end(), 0);
            gauss.resize(binMax); std::fill(      gauss.begin(),       gauss.end(), 0);
 
      ndata = 0;
      for(Integer_t i = 0;i < NormalizedFluc.size(); ++i) {
        if( isnan(NormalizedFluc[i]) ) continue;
        bin = (Integer_t)( (NormalizedFluc[i] -minFluc)/binwidth );
        if(bin >= binMax || bin<0) {
          std::cerr << "      # Warning: bin(" << bin << ") > binMax(" << binMax << "): NormalizedFluc[" << i << "]=" << NormalizedFluc[i] << std::endl;
          std::exit(EX_DATAERR);
        }
        histogram[bin] += 1;
        ndata += 1;
      }
      // if(ndata != Ndata[n]) std::cerr << "# Error: ndata(" << ndata << ") != Ndata[" << n << "](" << Ndata[n] << ")" << std::endl;
      for(bin = 0;bin < binMax; ++bin) {
        ProbDensity[bin] = histogram[bin]/((double)Ndata[n]*binwidth);
        if( isinf(ProbDensity[bin]) ) ProbDensity[bin] = NAN;
      }
      //********** (END) Create histogram ********** //
 
      //********** Create discrete Gauss distribution ********** //
      dtemp = 0;
      for(bin = 0;bin < binMax; ++bin) {
        Fluc = (bin+0.5)*binwidth +minFluc;
        gauss[bin] = NormalDist(Fluc);
        dtemp += gauss[bin];
      }
      for(bin = 0;bin < binMax; ++bin) gauss[bin] /= dtemp*binwidth;
      //********** (END) Create discrete Gauss distribution ********** //
 
      //********** Write histogram to a file ********** //
      for(bin = 0;bin < binMax; ++bin) {
        Fluc  = (bin+0.5)*binwidth +minFluc;
        OutHistfs << std::noshowpos
                  << std::setw(2)  << n     << " "
                  << std::showpos
                  << std::setw(13) << Fluc << " "
                  << std::setw(13) << ProbDensity[bin] << " "
                  << std::noshowpos
                  << std::setw(2)  << histogram[bin] << " "
                  << std::showpos
                  << std::setw(13) << gauss[bin] << "\n";
      }
      OutHistfs << std::endl;
      //********** (END) Write histogram to a file ********** //
 
      //********** Write distance to a file ********** //
      dist.calculate(ProbDensity, gauss, minFluc, binwidth);
      OutDevfs << count << " "
               << n     << " "
               << dist.HistIntersection << " "
               << dist.KLDivergence     << " "
               << dist.JSDivergence     << " "
               << dist.L1norm           << " "
               << dist.L2norm           << "\n";
      //********** (END) Write distance to a file ********** //
    }
    OutDevfs << std::endl;
  }
 
  return 0;
}

NormalDist()

double NormalDist

(

double

x

)

                           {
  return 1.0/sqrt(2*M_PI)*exp(-x*x/2.0);
}

Variable Documentation

filesystem

namespace filesystem = std::experimental::filesystem