ABACUS/src/EXECS/Histogram_RAW_File.cc

/**********************************************************

This software is part of J.-S. Caux's ABACUS library.

Copyright (c) J.-S. Caux.

-----------------------------------------------------------

File:  Histogram_RAW_File.cc

Purpose:  Produce a histogram of matrix elements from a RAW file.

***********************************************************/

#include "ABACUS.h"

using namespace std;
using namespace ABACUS;


int main(int argc, char* argv[])
{
  if (argc != 7) {
    cout << "Argument needed: rawfile, iKmin, iKmax, omega, domega, logscale_int" << endl;
    ABACUSerror("");
  }

  int n = 1;
  const char* rawfilename = argv[n++];
  int iKmin = atoi(argv[n++]);
  int iKmax = atoi(argv[n++]);
  DP om = atof(argv[n++]);
  DP dom = atof(argv[n++]);
  int logscale_int = atoi(argv[n++]); // see below for meaning

  ifstream RAW_infile;
  RAW_infile.open(rawfilename);
  if (RAW_infile.fail()) {
    cout << rawfilename << endl;
    ABACUSerror("Could not open RAW_infile... ");
  }

  stringstream outfilename;
  string outfilename_string;
  outfilename << rawfilename << "_hist";
  if (iKmin != iKmax) outfilename << "_iKmin_" << iKmin << "_iKmax_" << iKmax;
  outfilename << "_om_" << om << "_dom_" << dom << "_li_" << logscale_int;
  outfilename_string = outfilename.str();
  const char* outfilename_c_str = outfilename_string.c_str();

  ofstream outfile;
  outfile.open(outfilename_c_str);
  outfile.precision(16);



  // Use the binning in current ABACUS version's threads,
  // with logscale unit = 2^{1/64} \simeq 1.01
  int npts = 6400/logscale_int;
  DP logscale_used = logscale_int * (1.0/64) * log(2.0);

  DP omega_min = om - 0.5* dom;
  DP omega_max = om + 0.5* dom;

  DP omega;
  int iK;
  DP FF;
  DP dev;
  string label;

  Vect<int> nFF (0, npts);

  int nread = 0;
  int nfound = 0;
  int il = 0;

  while (RAW_infile.peek() != EOF) {
    RAW_infile >> omega >> iK >> FF >> dev >> label;

    nread++;

    if (iK >= iKmin && iK <= iKmax && omega > omega_min && omega < omega_max) {
      nfound++;
      il = int(-log(FF*FF)/logscale_used);
      if (il < 0) il=0;
      if (il > npts) continue;
      nFF[il] += 1;
    }
  }
  RAW_infile.close();

  if (nfound == 0) {
    cout << "Didn't find any contributing state." << endl;
    return(1);
  }

  for (int i = 0; i < npts; ++i) {
    if (i > 0) outfile << endl;
    outfile << i << "\t" << exp(-i * logscale_used) << "\t" << nFF[i] << "\t" << DP(nFF[i])/nfound;
  }
  outfile.close();


  // Now that we know the total number of entries, produce animated data:
  // the total number of entries is sliced into 100 pieces
  // and binned in the same way as before

  int stepsize = nfound/100;

  outfilename << "_sliced";
  outfilename_string = outfilename.str();
  const char* sliced_c_str = outfilename_string.c_str();

  ofstream sliced_file;
  sliced_file.open(sliced_c_str);
  sliced_file.precision(16);

  RAW_infile.open(rawfilename);
  if (RAW_infile.fail()) {
    cout << rawfilename << endl;
    ABACUSerror("Could not open RAW_infile... ");
  }

  nread = 0;
  nfound = 0;
  int nsteps = 0;
  for (int i = 0; i < npts; ++i) nFF[i] = 0;

  while (RAW_infile.peek() != EOF) {
    RAW_infile >> omega >> iK >> FF >> dev >> label;

    nread++;

    if (iK >= iKmin && iK <= iKmax && omega > omega_min && omega < omega_max) {
      nfound++;
      il = int(-log(FF*FF)/logscale_used);
      if (il < 0) il=0;
      if (il > npts) continue;
      nFF[il] += 1;
    }

    if (nfound == stepsize) {// || RAW_infile.peek() == EOF) {
      if (nsteps++ > 0) sliced_file << endl;
      sliced_file << DP(nFF[0])/stepsize;
      for (int i = 1; i < npts; ++i) sliced_file << "\t" << DP(nFF[i])/stepsize;
      for (int i = 0; i < npts; ++i) nFF[i] = 0;
      nfound = 0;
    }
  }

  RAW_infile.close();

  sliced_file.close();

  return(0);
}