From d922a0dae1c283752172e3ba840579044d6ee996 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jean-S=C3=A9bastien?= <git@jscaux.org>
Date: Mon, 6 Dec 2021 21:18:17 +0100
Subject: [PATCH] Add basic benchmarks facilities

---
 scripts/plot_benchmarks.py      |  30 +++++++
 src/EXECS/Benchmark_RAW_File.cc | 142 ++++++++++++++++++++++++++++++++
 2 files changed, 172 insertions(+)
 create mode 100644 scripts/plot_benchmarks.py
 create mode 100644 src/EXECS/Benchmark_RAW_File.cc

diff --git a/scripts/plot_benchmarks.py b/scripts/plot_benchmarks.py
new file mode 100644
index 0000000..f8334a1
--- /dev/null
+++ b/scripts/plot_benchmarks.py
@@ -0,0 +1,30 @@
+#!/usr/bin/python
+
+"""
+Plots benchmark files produced from Benchmark_RAW_File.
+
+Usage: python plot_benchmarks.py [benchmark file name].
+
+Example: python plot_benchmarks.py LiebLin_Rho_Rho_c_4_L_64_N_64_64_0__iK_32.raw_bmk
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+import sys
+
+filename = str(sys.argv[1])
+
+index, srcont = np.loadtxt(filename, usecols=(0,1), unpack=True)
+index_srt, srcont_srt = np.loadtxt(filename + '_srt', usecols=(0, 1), unpack=True)
+
+plt.semilogy(index, srcont, '.', markersize=1, label='Realized order')
+plt.semilogy(index, srcont_srt, '.', markersize=1, label='Ideal order') # careful: index, not index_srt 
+
+plt.xlabel('Order')
+plt.ylabel('ln(ME)')
+plt.title('Comparison of realized and ideal order\n' + filename)
+plt.legend()
+
+plt.savefig(filename.replace('.', '_') + '_plot.png')
+
+plt.show()
diff --git a/src/EXECS/Benchmark_RAW_File.cc b/src/EXECS/Benchmark_RAW_File.cc
new file mode 100644
index 0000000..b071b52
--- /dev/null
+++ b/src/EXECS/Benchmark_RAW_File.cc
@@ -0,0 +1,142 @@
+/**********************************************************
+
+This software is part of J.-S. Caux's ABACUS library.
+
+Copyright (c) J.-S. Caux.
+
+-----------------------------------------------------------
+
+File:  Benchmark_RAW_File.cc
+
+Purpose:  produces benchmarking files.
+    Two files are produced, each with line output: label sr_cont order,
+    - one in order of original RAW file
+    - one in order of decreasing sr_cont
+
+***********************************************************/
+
+#include "ABACUS.h"
+
+using namespace std;
+using namespace ABACUS;
+
+
+namespace ABACUS {
+
+  void Benchmark_RAW_File (const char ff_file[], char whichDSF)
+  {
+    // Check size of raw file:
+    struct stat statbuf;
+
+    stat (ff_file, &statbuf);
+    int filesize = statbuf.st_size;
+
+    // Determine the number of entries approximately
+    int entry_size = 2* sizeof(float) + 2*sizeof(int);
+
+    //const int MAXDATA = 50000000;
+    const int MAXDATA = filesize/entry_size + 10;
+
+    DP* omega = new DP[MAXDATA];
+    int* iK = new int[MAXDATA];
+    DP* ff = new DP[MAXDATA];
+    DP* ff_im = new DP[MAXDATA];
+    DP* dev = new DP[MAXDATA];
+    string* label = new string[MAXDATA];
+
+    ifstream infile;
+    infile.open(ff_file);
+
+    if (infile.fail()) ABACUSerror("The input file was not opened successfully in Benchmark_RAW_File. ");
+
+    int Ndata = 0;
+    while (((infile.peek()) != EOF) && (Ndata < MAXDATA)) {
+
+      infile >> omega[Ndata];
+      infile >> iK[Ndata];
+      if (whichDSF != 'Z') infile >> ff[Ndata];
+      if (whichDSF == 'q') infile >> ff_im[Ndata];  // imaginary part of overlap, quench case
+      infile >> dev[Ndata];
+      infile >> label[Ndata];
+
+      Ndata++;
+    }
+    infile.close();
+
+    int* index = new int[Ndata];
+    // Sum rule contribution
+    DP* sr_cont = new DP[MAXDATA];
+
+    for (int i = 0; i < Ndata; ++i) index[i] = i;
+
+    if (whichDSF == 'd') {
+      for (int i = 0; i < Ndata; ++i) sr_cont[i] = fabs(omega[i]) * ff[i] * ff[i];
+    }
+
+    stringstream outfilename;
+    string outfilename_string;
+    outfilename << ff_file << "_bmk";
+    outfilename_string = outfilename.str();
+    const char* outfilename_c_str = outfilename_string.c_str();
+
+    ofstream outfile1;
+    outfile1.open(outfilename_c_str);
+    outfile1.precision(16);
+
+    for (int i = 0; i < Ndata; i++) {
+
+      if (i > 0) outfile1 << endl;
+      outfile1 << index[i] << "\t" << sr_cont[i] << "\t" << label[i];
+    }
+    outfile1.close();
+
+    // Now the ordered one
+    QuickSort(sr_cont, index, 0, Ndata - 1);
+
+    outfilename << "_srt";
+    outfilename_string = outfilename.str();
+    const char* outfilename2_c_str = outfilename_string.c_str();
+
+    ofstream outfile2;
+    outfile2.open(outfilename2_c_str);
+    outfile2.precision(16);
+
+    for (int i = 0; i < Ndata; i++) {
+
+      if (i > 0) outfile2 << endl;
+      outfile2 << index[Ndata - 1 - i] << "\t" << sr_cont[Ndata - i - 1] << "\t" << label[index[Ndata - 1 - i] ];
+    }
+    outfile2.close();
+
+    delete[] omega;
+    delete[] iK;
+    delete[] ff;
+    delete[] ff_im;
+    delete[] dev;
+    delete[] label;
+
+    delete[] index;
+    delete[] sr_cont;
+
+    return;
+  }
+
+} // namespace ABACUS
+
+
+
+
+int main(int argc, char* argv[])
+{
+  if (argc != 3) {
+    cout << "Arguments needed: rawfile, whichDSF." << endl;
+    ABACUSerror("");
+  }
+
+  const char* rawfilename = argv[1];
+  char whichDSF = *argv[2];
+
+  Benchmark_RAW_File (rawfilename, whichDSF);
+
+  return(0);
+}