SrednickiAnsatz.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
	)

                               {
  if(argc < 5) {
    std::cerr << "Usage: 1.This 2.Directory 3.Nmax 4.Nmin 5.dE (Op1).SampleMin (Op2).SampleMax" << std::endl;
    std::exit(EX_USAGE);
  }
  constexpr int NDIV = 100;
     int const Nmax = (argc>=3 ? std::stoi(argv[2]) :16);
     int const Nmin = (argc>=4 ? std::stoi(argv[3]) :16);
  double const   dE = (argc>=5 ? std::stod(argv[4]) :0.02);
     int const SampleMin = (argc>=6 ? std::stod(argv[5]) :-1);
     int const SampleMax = (argc>=7 ? std::stod(argv[6]) :-1);
  if(Nmax < Nmin) {
    std::cerr << "Error: Nmax(" << Nmax << ") < Nmin(" << Nmin << ")" << std::endl;
    std::exit(EX_USAGE);
  }
 
 
  std::string const baseDir(argv[1]);
  std::cout << "# Directory: " << baseDir << "\n";
  std::string InDataDir(baseDir); InDataDir += "/RawData";
  std::vector<std::string> dirList;
  {
    std::string str;
    int dataNo;
    int const idShift=(SampleMin>0? SampleMin: 0);
    if(SampleMin>0 && SampleMax>0) dirList.resize(SampleMax-idShift+1);
    for(const filesystem::directory_entry &entry : filesystem::directory_iterator(InDataDir)) {
      if ( entry.is_directory() ) {
        str = entry.path().string();
        dataNo = std::atoi(str.substr(str.find("_No")+3).c_str());
        if((SampleMin>0 || SampleMax>0) && (dataNo<SampleMin || SampleMax<dataNo)) continue;
        debug_print("# dir=" << str << ", dataNo=" << dataNo);
        if( (SampleMin<=0 || SampleMax<=0)  && dirList.size()<=dataNo) dirList.resize(dataNo+1);
        debug_print("dataNo-idShift=" << dataNo-idShift << ", dataNo=" << dataNo << ", idShift=" << idShift);
        dirList.at(dataNo-idShift) = str;
      }
    }
  }
 
 
  #pragma omp parallel for schedule(dynamic, 10)
  for(size_t sample = 0;sample < dirList.size(); ++sample) {
  // for(size_t sample = 0;sample < 1; ++sample) {
    int n, bin, binMax, binMin;
    Integer_t dim;
    std::string command;
    double OpRange, dtemp;
    std::vector<double> eigenEnergy, eigenExpValue, MCAverage, binStddev(NDIV);
    std::vector<Integer_t> NdataInShell;
    std::vector<std::vector<Integer_t>> binElem(NDIV);
    constexpr auto Bin = [&](double x) { return (x==1 ? NDIV-1: (int)(x*NDIV)); };
    std::string InExpValueFname, InEnergyFname, InMatElemsFname, str;
    std::ifstream InExpValueFs, InEnergyFs, InMatElemsFs;
    std::string OutFname;
    std::ofstream OutFs;
    std::vector<std::string> candidateFileNames;
 
 
    // No = std::stoi(entry.path().filename().string<char>().substr(9));
    std::string OutDir(std::regex_replace(dirList[sample], std::regex("/RawData"), "/ProcessedData"));
    filesystem::create_directories(OutDir);
 
    bool flag = false;
    for(n = Nmin;n <= Nmax; ++n) {
      dim = 0;
      std::stringstream buff(""); buff << "_N" << n << ".txt";
      { //-------------------- Load EigenExpValues --------------------//
        std::stringstream ss;
        candidateFileNames.resize(3);
        candidateFileNames[0] = dirList[sample] + "/EigenExpValueZ" + buff.str();
        candidateFileNames[1] = dirList[sample] + "/EigenExpValueC" + buff.str();
        candidateFileNames[2] = dirList[sample] + "/EigenExpValue"  + buff.str();
        if( InExpValueFs.is_open() ) { InExpValueFs.close(); InExpValueFs.clear(); }
        for(auto name : candidateFileNames) {
          InExpValueFname = name;
          InExpValueFs.open(InExpValueFname);
          if( InExpValueFs.is_open() ) break;
        }
        if( !InExpValueFs.is_open() ) continue;
        flag = true;
        InExpValueFs.exceptions(std::ios::badbit | std::ios::failbit);
        while( std::getline(InExpValueFs, str) ) if(str.find("OpRange") != std::string::npos) break;
        OpRange = std::stod(str.substr(str.find("OpRange=")+8));
 
          eigenEnergy.resize(0);
        eigenExpValue.resize(0);
        try{
          skip_header(InExpValueFs);
          for(getline(InExpValueFs, str); !InExpValueFs.eof(); getline(InExpValueFs, str)) {
            ss.str(str); ss.clear(std::stringstream::goodbit);
            ss >> dtemp; eigenEnergy.push_back(dtemp);
            ss >> dtemp;
            ss >> dtemp; eigenExpValue.push_back(dtemp/OpRange);
          }
        } catch(std::exception& e) {
          if( !InExpValueFs.eof() )
          std::cerr << "Error(" << InExpValueFs.eof() << InExpValueFs.fail() << InExpValueFs.bad() << "):" << e.what() << std::endl;
        }
        InExpValueFs.close(); InExpValueFs.clear();
        if(eigenEnergy.size() != eigenExpValue.size()) {
          std::cerr << "Error: eigenEnergy.size() != eigenExpValue.size() for a file\n"
                    << "\t" << InExpValueFname << std::endl;
          continue;
        }
        dim = eigenEnergy.size();
      }
      // if(dim == 0) { //-------------------- Load eigen-energies --------------------//
      //   bool xzFlag;
      //   std::string::size_type pos;
      //   double gE, EnergyRange;
      //   candidateFileNames.resize(2);
      //   candidateFileNames[0] = dirList[sample] + "/EigenEnergiesZ" + buff.str();
      //   candidateFileNames[1] = dirList[sample] + "/EigenEnergiesC" + buff.str();
      //   if( InEnergyFs.is_open() ) { InEnergyFs.close(); InEnergyFs.clear(); }
      //   for(auto name : candidateFileNames) {
      //     InEnergyFname = name;
      //     InEnergyFs.open(InEnergyFname);
      //     if( InEnergyFs.is_open() ) break;
      //   }
      //   if( !InEnergyFs.is_open() ) continue;
      //   while(std::getline(InEnergyFs, str)) if( (pos = str.find("dim = ")) != std::string::npos) break;
      //   dim = std::stoi(str.substr(pos+6));
      //   eigenEnergy.resize(dim);
      //   for(int index = 0;index<dim && std::getline(InEnergyFs, str); ++index) eigenEnergy.at(index) = std::stod(str);
      //   debug_print("\t" << "  InEnergyFname=" << InEnergyFname);
      //   //-------------------- Load eigen matrix elements --------------------//
      //   candidateFileNames.resize(4);
      //   candidateFileNames[0] = dirList[sample] + "/EigenMatrixElementsZ" + buff.str();
      //   candidateFileNames[1] = dirList[sample] + "/EigenMatrixElementsC" + buff.str();
      //   candidateFileNames[2] = dirList[sample] + "/EigenMatrixElementsZ" + buff.str() + ".xz";
      //   candidateFileNames[3] = dirList[sample] + "/EigenMatrixElementsC" + buff.str() + ".xz";
      //   if( InMatElemsFs.is_open() ) { InMatElemsFs.close(); InMatElemsFs.clear(); }
      //   for(auto name : candidateFileNames) {
      //     InMatElemsFname = name;
      //     InMatElemsFs.open(InMatElemsFname);
      //     if( InMatElemsFs.is_open() ) break;
      //   }
      //   if( !InMatElemsFs.is_open() ) continue;
      //   if( (pos = InMatElemsFname.find(".xz")) != std::string::npos ) {
      //     xzFlag = true;
      //     command = "unxz -vk " + InMatElemsFname;
      //     if(system(command.c_str()) == 1) {
      //       std::cerr << "Error: Command \"" << command << "\" failed." << std::endl;
      //       std::exit(-1);
      //     };
      //
      //     InMatElemsFs.close(); InMatElemsFs.clear();
      //     InMatElemsFname = InMatElemsFname.substr(0,pos);
      //     InMatElemsFs.open(InMatElemsFname); checkIsFileOpen(InMatElemsFs, InMatElemsFname);
      //   } else { xzFlag = false; }
      //   while(std::getline(InMatElemsFs, str, '\n')) if(str.find("OpRange") != std::string::npos) break;
      //   OpRange = std::stod(str.substr(str.find("OpRange=")+8));
      //   if( std::get<0>( readDiagonal(InMatElemsFs, eigenExpValue) ) != dim ) {
      //     std::cerr << "# Error: dimension does not match between the two input files.\n"
      //               << "         " << InEnergyFname << "\n"
      //               << "         " << InMatElemsFname << std::endl;
      //     continue;
      //   };
      //   debug_print("\t" << "InMatElemsFname=" << InMatElemsFname);
      //
      //   InMatElemsFs.close(); InMatElemsFs.clear();
      //   if(xzFlag) {
      //     command = "rm -f " + InMatElemsFname;
      //     if(system(command.c_str()) != 0) {
      //       std::cerr << "Error: Command \"" << command << "\" failed." << std::endl;
      //       std::exit(-1);
      //     };
      //   }
      //
      //   gE = *std::min_element(eigenEnergy.begin(), eigenEnergy.end());
      //   EnergyRange = *std::max_element(eigenEnergy.begin(), eigenEnergy.end()) - gE;
      //   for(size_t i = 0;i < eigenEnergy.size(); ++i) {
      //     eigenEnergy[i] = (eigenEnergy[i] - gE)/EnergyRange;
      //     eigenExpValue[i] /= OpRange;
      //   }
      // }
    #ifndef NDEBUG
      std::cout << std::scientific;
      std::cout << "# OpRange = " << OpRange << std::endl;
      print(eigenEnergy, eigenEnergy.size());
      print(eigenExpValue, eigenExpValue.size());
    #endif
      debug_print("########## Some process ##########");
 
      // Calculate Microcanonical averages with energy corresponding to each eigenstate.
      MCAverage.resize(dim); NdataInShell.resize(dim);
      std::fill(   MCAverage.begin(),    MCAverage.end(), 0.0);
      std::fill(NdataInShell.begin(), NdataInShell.end(), 0.0);
      for(Integer_t i = 0;i < dim; ++i) {
        MCAverage[i] = MicroCanonicalAverage(NdataInShell[i], eigenEnergy[i], dE, dim, eigenEnergy, eigenExpValue, i);
      }
 
      for(size_t bin = 0;bin < NDIV; ++bin) binElem[bin].resize(0);
      for(size_t i = 0;i < dim; ++i) {
        if( NdataInShell[i] <= 1 ) continue;
        bin = Bin(eigenEnergy[i]);
        for(binMin = bin;binMin>=0   && std::abs(eigenEnergy[i] - (binMin+0.5)/NDIV) < dE; --binMin) {}; ++binMin;
        for(binMax = bin;binMax<NDIV && std::abs(eigenEnergy[i] - (binMax+0.5)/NDIV) < dE; ++binMax) {}; --binMax;
        if(binMin < 0) binMin = 0;
        debug_print("binMin=" << binMin << ", binMax=" << binMax << ", Energy=" << eigenEnergy[i]);
        for(bin = binMin;bin <= binMax; ++bin) binElem.at(bin).push_back(i);
      }
      std::fill(binStddev.begin(), binStddev.end(), 0.0);
 
      #pragma omp parallel for
      for(size_t bin = 0;bin < NDIV; ++bin) {
        for(auto const& i : binElem[bin]) {
          double dtemp = eigenExpValue[i] - MCAverage[i];
          binStddev[bin] += dtemp*dtemp;
        }
        if(binElem[bin].size() == 0) binStddev[bin] = NAN;
        else binStddev[bin] = sqrt(binStddev[bin]/binElem[bin].size());
      }
      debug_print("########## (END)Some process ##########\n");
 
      {
        std::stringstream buff(""); buff << "_N" << n << ".txt";
        OutFname = OutDir + "/SrednickiAnsatz_Diagonal" + buff.str();
      }
      if( OutFs.is_open() ) { OutFs.close(); OutFs.clear(); }
      OutFs.open(OutFname); checkIsFileOpen(OutFs, OutFname);
      OutFs << "# NDIV= " << NDIV << "\n"
            << "# shellWidth= " << dE << "\n"
            << "# 1.(N), 2.(Normalized energy), 3.(ExpValue), 4.(Stddev), 5.(funcF), 6.(DOS), 7.(Ndata)" << "\n\n";
      OutFs << std::scientific << std::showpos;
      for(size_t bin = 0;bin < NDIV; ++bin) {
        OutFs << n              << " "
              << (bin+0.5)/NDIV << " "
              << "nan"          << " "
              << binStddev[bin] << " "
              << binStddev[bin]*sqrt((double)binElem[bin].size()) << " "
              << "nan"          << " "
              << binElem[bin].size()
              << std::endl;
      }
      OutFs.close();
    }
    if(flag) std::cout << "# " << dirList[sample] << std::endl;
  }
 
  return EXIT_SUCCESS;
}

Variable Documentation

filesystem

namespace filesystem = std::experimental::filesystem