OpAverage.cpp File Reference

Classes
class	Histogram

Functions
int	main (int argc, char **argv)

Variables
namespace	filesystem = std::experimental::filesystem
Integer_t const	NBIN = 100

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);
  }
  std::string const d_str(argv[1]);
  double const shellWidth = (double)std::atof(argv[2]);
  std::cout << "Directory: " << d_str << "\n";
 
  Integer_t n, Nmin, Nmax=0, bin, count=0, Nline;
  double relEnergy, expval;
  double Emin, Emax, StddevInShell, dtemp;
  std::vector<Integer_t>   Ndata, NdataInShell;
  std::vector<double> MCAverage;
  std::vector<std::vector<double>> eigenEnergy;
  std::vector<std::vector<double>>    ExpValue;
  std::vector<Histogram> Average;
  Histogram BinData;
 
  std::string line;
 
  //********** Output fileを作成 ********** //
  std::string OutProcDIR(d_str); OutProcDIR += "/ProcessedData/";
  filesystem::create_directory(OutProcDIR);
  std::string   OutProcFName;
  //********** (END)Output fileを作成 ********** //
 
  //********** Input directory ********** //
  std::string InDataDir(d_str); InDataDir += "/RawData";
  for(const std::filesystem::directory_entry &data : std::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;
    }
 
    //********** Read data **********//
    ++count;
    Nmin=INT_MAX; Nmax=INT_MIN;
    std::fill(Ndata.begin(),Ndata.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 **********//
 
    OutProcFName = OutProcDIR + data.path().filename().string();
    std::ofstream OutProcfs(OutProcFName);
    OutProcfs.exceptions(std::ios::badbit | std::ios::failbit);
      if( !OutProcfs.is_open() ) {
        std::cerr << "Can't open a file " << OutProcFName << std::endl;
        std::exit(EX_CANTCREAT);
      }
    OutProcfs << "# NBIN = "       << NBIN       << "\n"
              << "# shellWidth = " << shellWidth << "\n"
              << "# 1.(N), 2.(ratio), 3.(ExpValue), 4.(Stddev), 5.(funcF), 6.(DOS), 7.(Ndata), 8.(Nan data)\n"
              << std::scientific
              << std::endl;
    std::cout << "(count=" << count << ") Output file(Histogram): " << OutProcFName << std::endl;
 
    if(Nmax >= Average.size()) Average.resize(Nmax+1);
    for(n = Nmax;n >= Nmin; --n) {
      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);
      }
 
      BinData.clear();
      for(int i = 0;i < Ndata[n]; ++i) {
        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;
        }
        BinData.add(eigenEnergy[n][i], ExpValue[n][i], StddevInShell, NdataInShell[i]);
      }
 
      for(bin = 0;bin < NBIN; ++bin) {
 
        BinData.DOS[bin] = BinData.Ndata[bin]*NBIN/(double)BinData.Ntot;
        OutProcfs << std::noshowpos
                  << std::setw(2)  << n                      << " "
                  << std::showpos
                  << std::setw(13) << (bin+0.5)/(double)NBIN << " "
                  << std::setw(13) << BinData.ExpValue[bin]  << " "
                  << std::setw(13) << BinData.Stddev[bin]    << " "
                  << std::setw(13) << BinData.funcF[bin]     << " "
                  << std::setw(13) << BinData.DOS[bin]       << " "
                  << std::setw(4)  << BinData.Ndata[bin]     << "\n";
      }
      OutProcfs << std::endl;
 
      if(BinData.Ntot > 0 && n==Nmax) Average.at(n).add(BinData);
    }
    OutProcfs.close();
  }
  //********** (END) Input directory ********** //
 
 
  //********** アンサンブル平均のファイルへの書き込み **********//
  std::string OutDIR(d_str); OutDIR += "/Output";
  filesystem::create_directory(OutDIR);
  std::stringstream buff;
  buff.str("");
  buff.clear(std::stringstream::goodbit);
  buff << "/EnsembleAverage" << ".txt";
  std::string   OutAveFName(OutDIR); OutAveFName += buff.str();
  std::ofstream OutAvefs(OutAveFName);
  OutAvefs.exceptions(std::ios::badbit | std::ios::failbit);
    if( !OutAvefs.is_open() ) {
      std::cerr << "Can't open a file " << OutAveFName << std::endl;
      std::exit(EX_CANTCREAT);
    }
  OutAvefs << "# NBIN = "       << NBIN         << "\n"
           << "# Nsample = "    << Average[Nmax].Ntot << "\n"
           << "# shellWidth = " << shellWidth   << "\n"
           << "# 1.(N), 2.(ratio), 3.(ExpValue), 4.(Stddev), 5.(funcF), 6.(DOS), 7.(Ndata), 8.(Nan data)\n"
           << std::scientific
           << std::endl;
  std::cout << "Output file: " << OutAveFName << std::endl;
 
  Integer_t flag;
  for(n = Nmax;n >= Nmin; --n) {
    flag = 0;
    for(bin = 0;bin < NBIN; ++bin) {
      if( Average[n].Ndata[bin] == 0) continue;
      flag = 1;
      OutAvefs << std::noshowpos
               << std::setw(2)  << n << " "
               << std::showpos
               << std::setw(13) << (bin+0.5)/(double)NBIN   << " "
               << std::setw(13) << Average[n].ExpValue[bin] << " "
               << std::setw(13) << Average[n].Stddev[bin]   << " "
               << std::setw(13) << Average[n].funcF[bin]    << " "
               << std::setw(13) << Average[n].DOS[bin]      << " "
               << std::setw(13) << Average[n].Ndata[bin]    << " "
               << std::setw(4)  << Average[n].Nnan[bin]     << "\n";
    }
    if(flag == 1) OutAvefs << std::endl;
  }
  OutAvefs.close();
  //********** (END)アンサンブル平均のファイルへの書き込み **********//
 
  return 0;
}

Variable Documentation

filesystem

namespace filesystem = std::experimental::filesystem

NBIN

Integer_t const NBIN = 100