Add histogram-making functions and plotting scripts

Makefile

@@ -123,6 +123,7 @@ EXECS = $(BINDIR)LiebLin_DSF $(BINDIR)LiebLin_Data_Daemon $(BINDIR)LiebLin_RAW_F
 lib$(VERSION).a : $(Objects_ALL)
 	ar -cru lib$(VERSION).a $(Objects_ALL)
 	mv lib$(VERSION).a $(BASEDIR)lib/
+	$(COMPILE) $(EXECSDIR)Histogram_RAW_File.cc -o $(BINDIR)Histogram_RAW_File -l$(VERSION)
 	$(COMPILE) $(EXECSDIR)LiebLin_DSF.cc -o $(BINDIR)LiebLin_DSF -l$(VERSION)
 	$(COMPILE) $(EXECSDIR)LiebLin_Data_Daemon.cc -o $(BINDIR)LiebLin_Data_Daemon -l$(VERSION)
 	$(COMPILE) $(EXECSDIR)LiebLin_Data_Daemon_Nscaling.cc -o $(BINDIR)LiebLin_Data_Daemon_Nscaling -l$(VERSION)

scripts/DSF_iK.py

@@ -0,0 +1,20 @@
+import argparse
+import plotly.graph_objects as go
+import numpy
+
+parser = argparse.ArgumentParser()
+parser.add_argument('omegafile', help='Omega filename')
+parser.add_argument('dsffile', help='DSF filename')
+args = parser.parse_args()
+
+omega = numpy.loadtxt(args.omegafile)
+dsf = numpy.loadtxt(args.dsffile)
+
+x = [o for o in omega]
+y = [d for d in dsf]
+
+fig = go.Figure(data=go.Scatter(x=x, y=y))
+fig.update_layout(title=args.dsffile.rpartition('/')[2])
+fig.show()
+
+

scripts/Histograms.py

@@ -0,0 +1,21 @@
+import argparse
+import plotly.graph_objects as go
+import numpy
+
+parser = argparse.ArgumentParser()
+parser.add_argument('filenames', nargs='*')
+args = parser.parse_args()
+
+print(args)
+
+fig = go.Figure()
+
+for datafilename in args.filenames:
+	data = numpy.loadtxt(datafilename, delimiter="\t", usecols=[0,1,2,3])
+	x = [line[0] for line in data]
+	y = [line[3] for line in data]
+	fig.add_trace(go.Bar(x=x, y=y))
+
+fig.show()
+
+

scripts/Histograms_Animated.py

@@ -0,0 +1,38 @@
+import argparse
+import plotly.graph_objects as go
+import numpy
+
+parser = argparse.ArgumentParser()
+parser.add_argument('filename', help='Enter the sliced data filename')
+args = parser.parse_args()
+
+data = numpy.loadtxt(args.filename, delimiter="\t")
+# The data is in row format
+x = [i for i in range(len(data[0]))]
+
+fig = go.Figure()
+
+for line in data:
+	fig.add_trace(go.Bar(x=x, y=[d for d in line]))
+
+steps = []
+for i in range(len(fig.data)):
+	step = dict(
+		method="update",
+		args=[{'visible': [False] * len(fig.data)},
+			{'title': "Slider switched to step: " + str(i)}],)
+	step["args"][0]["visible"][i] = True
+	steps.append(step)
+
+sliders = [dict(
+	active=0,
+	currentvalue={"prefix": "bin box "},
+	pad={"t": 50},
+	steps=steps)]
+
+fig.update_yaxes(range=[0,0.2])
+fig.update_layout(
+	sliders=sliders
+)
+fig.show()
+

src/EXECS/Histogram_RAW_File.cc

@@ -0,0 +1,156 @@
+/**********************************************************
+
+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);
+}