fix runtime validation checks in tools

Replace assert-based checks that validate CLI inputs and file/data shape with exceptions so they remain active in release builds. Keep internal invariant assertions unchanged and preserve minimal diffs without broad formatting edits.

Author: rokzitko

c++/coef.hpp

@@ -3,6 +3,7 @@
 
 #include <vector>
 #include <fstream>
+#include <stdexcept>
 #include "traits.hpp"
 #include "params.hpp"
 #include "numerics.hpp" // read_vector
@@ -41,10 +42,10 @@ class coef_table {
 template <scalar S> auto read_matrix_table(std::istream &F, const bool nr_is_max_index = false) {
   const auto nr = read_one<size_t>(F);
   const auto len = nr_is_max_index ? nr+1 : nr;
-  assert(len >= 1);
+  if (len < 1) throw std::runtime_error("read_matrix_table(): expected at least one matrix.");
   const auto size1 = read_one<size_t>(F);
   const auto size2 = read_one<size_t>(F);
-  assert(size1 == size2);
+  if (size1 != size2) throw std::runtime_error("read_matrix_table(): expected square matrices.");
   using t_matrix = Matrix_traits<S>;
   std::vector<t_matrix> vec(len);
   for (auto j = 0; j < len; j++)

tools/adapt/adapt.hpp

@@ -14,6 +14,7 @@
 #include <map>
 #include <algorithm>
 #include <functional>
+#include <stdexcept>
 
 using namespace std;
 using namespace std::string_literals;
@@ -285,24 +286,26 @@ class Adapt {
    void set_parameters(const int PREC = 16) {
      std::cout << std::setprecision(PREC);
      Lambda = LAMBDA(P.P("Lambda", 2.0));
-     assert(Lambda > 1.0);
+     if (!(Lambda > 1.0)) throw std::invalid_argument("Lambda must be greater than 1.");
      adapt = P.Pbool("adapt", false); // Enable adaptable g(x)? Default is false!!
      hardgap  = P.Pbool("hardgap", false); // Exclude an interval around omega=0 ?
      boundary = P.P("boundary", 0.0);      // The boundary of the exclusion interval.
      bandrescale = P.P("bandrescale", 1.0); // band rescaling parameter
      xmax = P.P("xmax", 30); // Integrate over [1..xmax]
-     assert(xmax > 1);
+     if (!(xmax > 1)) throw std::invalid_argument("xmax must be greater than 1.");
      xfine = P.P("xfine", 5); // Fine stepsize integral [1..xfine]
-     assert(xfine > 1);
+     if (!(xfine > 1)) throw std::invalid_argument("xfine must be greater than 1.");
      output_step = P.P("outputstep", 1.0 / 64.0); // Stepsize for output file
-     assert(output_step <= 1.0);
+     if (!(output_step <= 1.0)) throw std::invalid_argument("outputstep must be less than or equal to 1.");
      dx_fine       = P.P("dx_fine", 1e-5);        // Integration stepsize in [1..xfine]
      dx_fast       = P.P("dx_fast", 1e-4);        // Integration stepsize in [xfine..xmax]
      allowed_error = P.P("allowed_error", 1e-10); // error control for adaptable stepsize
      max_subdiv    = P.Pint("max_subdiv", 10);    // maximum nr of integ. step subdivisions
-     assert(dx_fine * pow(0.5, max_subdiv) > DBL_EPSILON);
+     if (!(dx_fine * pow(0.5, max_subdiv) > DBL_EPSILON)) {
+        throw std::invalid_argument("dx_fine and max_subdiv imply a sub-step below machine precision.");
+      }
      max_abs = P.P("max_abs", 100.0); // Maximal |f(x)|
-     assert(max_abs > 0.0);
+     if (!(max_abs > 0.0)) throw std::invalid_argument("max_abs must be greater than 0.");
      convergence_eps = P.P("secant_eps", 1e-4);
      factor0         = 1.0 + P.P("secant_factor", 1e-7);
      max_iter        = P.Pint("secant_max_iter", 10);

tools/binavg/binavg.hpp

@@ -14,6 +14,7 @@
 #include <map>
 #include <string>
 #include <algorithm>
+#include <stdexcept>
 #include <unistd.h> // getopt
 
 namespace NRG::BinAvg {
@@ -79,7 +80,7 @@ class BinAvg {
      }
      name = std::string(argv[optind]); // Name of spectral density files
      Nz = atoi(argv[optind + 1]); // Number of z-values
-     assert(Nz >= 1);
+     if (!(Nz >= 1)) throw std::invalid_argument("Nz must be greater than or equal to 1.");
      std::cout << "Processing: " << name << std::endl;
      std::cout << "Nz=" << Nz << std::endl;
    }
@@ -101,7 +102,7 @@ class BinAvg {
      f.seekg(0, std::ios::end);
      const auto end_pos   = f.tellg();
      const long len       = end_pos - begin_pos;
-     assert(len % (rows * sizeof(double)) == 0);
+     if (len % (rows * sizeof(double)) != 0) throw std::runtime_error("Input binary file has incomplete row data.");
      const int nr = len / (rows * sizeof(double)); // number of lines
      if (verbose) { std::cout << "len=" << len << " nr=" << nr << " data points" << std::endl; }
      // Allocate the read buffer. The data will be kept in memory for the duration of the calculation!

tools/broaden/broaden.hpp

@@ -142,13 +142,13 @@ auto file_size(std::ifstream &f) {
 
 auto nr_doubles(std::ifstream &f) {
   const auto len = file_size(f);
-  assert(len % sizeof(double) == 0);
+  if (len % sizeof(double) != 0) throw std::runtime_error("Binary input size is not a multiple of sizeof(double).");
   return len/sizeof(double);
 }
 
 auto nr_rows(std::ifstream &f, const int cols) {
   const auto d = nr_doubles(f);
-  assert(d % cols == 0);
+  if (d % cols != 0) throw std::runtime_error("Binary input does not contain a whole number of rows.");
   return d / cols;
 }
 
@@ -313,25 +313,25 @@ class Broaden {
      }
      name = std::string(argv[optind]); // Name of spectral density files
      Nz = atoi(argv[optind + 1]); // Number of z-values
-     assert(Nz >= 1);
+     if (!(Nz >= 1)) throw std::invalid_argument("Nz must be greater than or equal to 1.");
      alpha = atof(argv[optind + 2]); // High-energy broadening parameter
-     assert(alpha > 0.0);
+     if (!(alpha > 0.0)) throw std::invalid_argument("alpha must be greater than 0.");
      T = atof(argv[optind + 3]); // Temperature
-     assert(T > 0.0);
+     if (!(T > 0.0)) throw std::invalid_argument("T must be greater than 0.");
      if (remaining == 5) {
        omega0_ratio = atof(argv[optind + 4]); // omega0/T
-       assert(omega0_ratio > 0.0);
+       if (!(omega0_ratio > 0.0)) throw std::invalid_argument("omega0_ratio must be greater than 0.");
        omega0 = omega0_ratio * T;
      }
-      if (remaining == 4) {
-        omega0_ratio = 1e-9; // Effectively zero
-        omega0       = omega0_ratio * T;
-      }
-      if (finalgaussian && ggamma <= 0.0) throw std::invalid_argument("Final Gaussian width must be greater than 0.");
-      if (finalderfd && dgamma <= 0.0) throw std::invalid_argument("Final derFD width must be greater than 0.");
-      if (!meshpositive && !meshnegative) throw std::invalid_argument("Output mesh cannot be empty.");
-      std::cout << "Processing: " << name << std::endl;
-      if (verbose) { std::cout << "Nz=" << Nz << " alpha=" << alpha << " T=" << T << " omega0_ratio=" << omega0_ratio << std::endl; }
+     if (remaining == 4) {
+       omega0_ratio = 1e-9; // Effectively zero
+       omega0       = omega0_ratio * T;
+     }
+     if (finalgaussian && ggamma <= 0.0) throw std::invalid_argument("Final Gaussian width must be greater than 0.");
+     if (finalderfd && dgamma <= 0.0) throw std::invalid_argument("Final derFD width must be greater than 0.");
+     if (!meshpositive && !meshnegative) throw std::invalid_argument("Output mesh cannot be empty.");
+     std::cout << "Processing: " << name << std::endl;
+     if (verbose) { std::cout << "Nz=" << Nz << " alpha=" << alpha << " T=" << T << " omega0_ratio=" << omega0_ratio << std::endl; }
     }
    void check_buffer_normalisation(const std::vector<double> &buffer, const int col_) {
      const auto cols = 1 + nrcol;

tools/bw/bw.cc

@@ -26,6 +26,7 @@
 #include <map>
 #include <string>
 #include <algorithm>
+#include <stdexcept>
 
 #include <unistd.h>
 #include <getopt.h>
@@ -147,9 +148,9 @@ void cmd_line(int argc, char *argv[]) {
   }
   name = string(argv[optind]);   // Name of spectral density files
   b0   = atof(argv[optind + 1]); // Parameter b0 (asymptotic broadening)
-  assert(b0 > 0.0);
+  if (!(b0 > 0.0)) throw std::invalid_argument("b0 must be greater than 0.");
   Nz = atoi(argv[optind + 2]); // Number of z-values
-  assert(Nz >= 1);
+  if (!(Nz >= 1)) throw std::invalid_argument("Nz must be greater than or equal to 1.");
 
   cout << "Processing: " << name << endl;
   cout << "b0=" << b0 << " Nz=" << Nz << endl;
@@ -185,7 +186,7 @@ void load(int i) {
   f.seekg(0, ios::end);
   const ios::pos_type end_pos = f.tellg();
   const long len              = end_pos - begin_pos;
-  assert(len % (2 * sizeof(double)) == 0);
+  if (len % (2 * sizeof(double)) != 0) throw std::runtime_error("Input binary file has incomplete data pairs.");
   const int nr = len / (2 * sizeof(double)); // number of pairs of double
   if (verbose) { cout << "len=" << len << " nr=" << nr << " data points" << endl; }
 

tools/kk/kk.hpp

@@ -27,6 +27,7 @@
 #include <cmath>
 #include <functional>
 #include <memory>
+#include <stdexcept>
 
 #include <gsl/gsl_errno.h>
 #include <gsl/gsl_spline.h>
@@ -135,7 +136,7 @@ class KK {
       if (im.empty()) throw std::runtime_error("No input data points provided.");
       std::sort(im.begin(), im.end());
       len = im.size();
-      assert(len % 2 == 0);
+      if (len % 2 != 0) throw std::runtime_error("Input grid must contain an even number of points.");
       std::tie (Xmin, Xmax) = x_range(im);
       if (mode == MODE::FILES) std::cout << "Range: [" << Xmin << " ; " << Xmax << "]" << std::endl;
       if (gsl_fcmp(-Xmin, Xmax, 1.e-8) != 0) throw std::runtime_error("Only symmetric intervals are supported!");
@@ -157,7 +158,9 @@ class KK {
       if (mode == MODE::FILES) std::cout << "Sum=" << sum << std::endl;
       const auto nr = Xpts.size()/2;
       for (auto i = nr; i < len; i++)
-        assert(gsl_fcmp(Xpts[len - i - 1], -Xpts[i], 1e-8) == 0);     // Check for the symmetry of the grid
+        if (gsl_fcmp(Xpts[len - i - 1], -Xpts[i], 1e-8) != 0) {
+          throw std::runtime_error("Input grid is not symmetric around zero.");
+        }
       Xpos = DVEC(nr);   // Xpos are positive and increasing!
       std::copy(Xpts.begin() + nr, Xpts.end(), Xpos.begin());
       w.reset(gsl_integration_workspace_alloc(workspace_limit));

tools/mats/mats.cc

@@ -20,6 +20,7 @@
 #include <string>
 #include <algorithm>
 #include <complex>
+#include <stdexcept>
 
 #include <unistd.h>
 #include <getopt.h>
@@ -93,11 +94,11 @@ void cmd_line(int argc, char *argv[]) {
   }
   name = string(argv[optind]);   // Name of spectral density files
   Nz   = atoi(argv[optind + 1]); // Number of z-values
-  assert(Nz >= 1);
+  if (!(Nz >= 1)) throw std::invalid_argument("Nz must be greater than or equal to 1.");
   T = atof(argv[optind + 2]); // Temperature
-  assert(T > 0.0);
+  if (!(T > 0.0)) throw std::invalid_argument("T must be greater than 0.");
   nrmats = atoi(argv[optind + 3]); // Number of Matsubara points
-  assert(nrmats >= 1);
+  if (!(nrmats >= 1)) throw std::invalid_argument("nrmats must be greater than or equal to 1.");
   cout << "Processing: " << name << endl;
   cout << "Nz=" << Nz << " T=" << T << " nrmats=" << nrmats << endl;
 }
@@ -131,7 +132,7 @@ void load(int i) {
   f.seekg(0, ios::end);
   const ios::pos_type end_pos = f.tellg();
   const long len              = end_pos - begin_pos;
-  assert(len % (rows * sizeof(double)) == 0);
+  if (len % (rows * sizeof(double)) != 0) throw std::runtime_error("Input binary file has incomplete row data.");
   const int nr = len / (rows * sizeof(double)); // number of lines
   if (verbose) cout << "len=" << len << " nr=" << nr << " data points" << endl;
   // Allocate the read buffer. The data will be kept in memory for the

tools/nrgchain/nrgchain.cc

@@ -15,6 +15,7 @@
 #include <algorithm>
 #include <ctime>
 #include <limits>
+#include <stdexcept>
 
 #include <gmp.h>
 
@@ -145,27 +146,27 @@ void set_parameters() {
   cout << setprecision(PREC);
 
   Lambda = LAMBDA(P("Lambda", 2.0));
-  assert(Lambda > 1.0);
+  if (!(Lambda > 1.0)) throw std::invalid_argument("Lambda must be greater than 1.");
 
   z = P("z", 1.0);
-  assert(0 < z && z <= 1.0);
+  if (!(0 < z && z <= 1.0)) throw std::invalid_argument("z must satisfy 0 < z <= 1.");
 
   adapt = Pbool("adapt", false); // Enable adaptable g(x)? Default is false!!
 
   bandrescale = P("bandrescale", 1.0);
   rescalexi   = Pbool("rescalexi", false);
 
   xmax = P("xmax", 30); // Interval [1..xmax]
-  assert(xmax >= 1.0);
+  if (!(xmax >= 1.0)) throw std::invalid_argument("xmax must be greater than or equal to 1.");
 
   Nmax = Pint("Nmax", 0); // Maximal site index in the Wilson chain
 
   mMAX = Pint("mMAX", 2 * Nmax); // Maximal index of coefficients (e,f)
 
-  assert(mMAX > 0);
+  if (!(mMAX > 0)) throw std::invalid_argument("mMAX must be greater than 0.");
 
   preccpp = Pint("preccpp", 2000); // Precision for GMP
-  assert(preccpp > 10);
+  if (!(preccpp > 10)) throw std::invalid_argument("preccpp must be greater than 10.");
 
   band = Pstr("band", "adapt"); // Default: load FSOL*.dat