From 15e9fdb5d488a6073d1268b9c1594a4737ac5333 Mon Sep 17 00:00:00 2001 From: "J.-S. Caux" Date: Fri, 15 May 2020 14:55:21 +0200 Subject: [PATCH] Add histogram-making functions and plotting scripts --- Makefile | 1 + scripts/DSF_iK.py | 20 ++++ scripts/Histograms.py | 21 +++++ scripts/Histograms_Animated.py | 38 ++++++++ src/EXECS/Histogram_RAW_File.cc | 156 ++++++++++++++++++++++++++++++++ 5 files changed, 236 insertions(+) create mode 100755 scripts/DSF_iK.py create mode 100755 scripts/Histograms.py create mode 100755 scripts/Histograms_Animated.py create mode 100644 src/EXECS/Histogram_RAW_File.cc diff --git a/Makefile b/Makefile index 65bafc2..2555e07 100644 --- a/Makefile +++ b/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) diff --git a/scripts/DSF_iK.py b/scripts/DSF_iK.py new file mode 100755 index 0000000..09e1fd0 --- /dev/null +++ b/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() + + diff --git a/scripts/Histograms.py b/scripts/Histograms.py new file mode 100755 index 0000000..0e117b4 --- /dev/null +++ b/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() + + diff --git a/scripts/Histograms_Animated.py b/scripts/Histograms_Animated.py new file mode 100755 index 0000000..c6e6114 --- /dev/null +++ b/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() + diff --git a/src/EXECS/Histogram_RAW_File.cc b/src/EXECS/Histogram_RAW_File.cc new file mode 100644 index 0000000..7b6135e --- /dev/null +++ b/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 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); +}