diff --git a/ALPS3_MIGRATION_PLAN.md b/ALPS3_MIGRATION_PLAN.md
index 5e9e573..2af9773 100644
--- a/ALPS3_MIGRATION_PLAN.md
+++ b/ALPS3_MIGRATION_PLAN.md
@@ -361,29 +361,62 @@ and Boost.Math in `legendre_convert`.
 
 ### 2.3 Code fixes
 
-* Fix B1 through B20 (§1.6), including B7b. B21 (logging convention and verbosity) can go here or into a later step; it changes the CLI regression references on purpose.
-* Moved here from 2.2 (they change the CLI references on purpose): replace
-  `boost::diagnostic_information` by `e.what()` (with B1), add a `SEED`
-  parameter for the bootstrap (B10), and delete the commented-out ublas and
-  LAPACK-bindings code.
-* `legendre_convert` (decided 2026-09-24), in this order:
-  1. Add regression cases for `legendre_convert` first (it has no tests), with
-     references from the current Boost build, like the `kk` cases in 2.2.
-  2. Replace Boost.Random by `<random>`. The `mt19937` engine gives identical
-     numbers; the normal variates differ (different algorithm), which is
-     harmless because the error estimate is seeded from the clock.
-  3. Replace `boost::math::factorial` by a product (only small arguments
-     occur).
-  4. Replace `boost::math::legendre_p` by the standard three-term recurrence.
-     `std::legendre` is not an option: libc++ (AppleClang) does not implement
-     the C++17 special math functions.
-  Keep `boost::math::sph_bessel` (also missing in libc++; our own version
-  would need careful checking at large l and argument) and `program_options`
-  (no standard equivalent). B5, B6 and B7 change reachable behavior, so each gets its own commit with a before/after test.
-* Remove all `using namespace boost::numeric;` lines, `#include <alps/config.hpp>`, dead commented-out ublas and lapack-bindings code, and the unused `alps::cast`.
-* Make `eigen_hdf5.hpp`/`eigen_lapack.hpp` functions `inline`, or move them into `.cpp` files (B8).
-* Use `Eigen::Index` for loop indices (B14). Consider `BDCSVD` in place of `JacobiSVD` (B11); that one is a numerics change and needs checking against the references.
-* Keep the `MaxEntSimulation` public getters stable, because the tests use them.
+Grouped by what each item does to the regression references (agreed
+2026-09-24). One PR per group, in this order, each with a Copilot review.
+Already done elsewhere: B2, B3 (2.1), B7 (`xi_bose_bugfix`), B20 (Pade removed).
+
+**2.3A: no change to results** (references must stay bit-identical)
+
+**Done 2026-09-24.** All Maxent-owned code and utilities build with
+`-Wall -Wextra -Wpedantic -Werror`; the fast and full regression suites remain
+unchanged. A text-input regression covers trailing blank lines.
+
+* Fix all compiler warnings (44 at the start of 2.3): B14 (`std::size_t` vs
+  `Eigen::Index` loops), B12 (member initialization order), unused variables,
+  misleading indentation, B13 (deprecated `params::get_origin_name()`; the
+  default `BASENAME` must not change). Then turn on `MAXENT_WERROR` in the `dev`
+  preset.
+* B8: make the functions in `eigen_hdf5.hpp`/`eigen_lapack.hpp` `inline`.
+* B4: make the size check in `Backcont::max_error` unconditional, avoiding an
+  out-of-bounds access in every build type.
+* B5, B6: remove the unreachable kernel branches.
+* B9: the data-file reading loop (`while (datstream)` reads past EOF); same
+  results for valid files, plus a test with trailing blank lines.
+* Remove the commented-out ublas and LAPACK-bindings code, `<alps/config.hpp>`
+  includes and the unused `alps::cast`.
+
+**2.3B: `legendre_convert`** (decided 2026-09-24), in this order:
+
+1. Add regression cases for `legendre_convert` first (it has no tests), with
+   references from the current Boost build, like the `kk` cases in 2.2.
+2. Replace Boost.Random by `<random>`. The `mt19937` engine gives identical
+   numbers; the normal variates differ (different algorithm), which is
+   harmless because the error estimate is seeded from the clock.
+3. Replace `boost::math::factorial` by a product (only small arguments occur).
+4. Replace `boost::math::legendre_p` by the standard three-term recurrence.
+   `std::legendre` is not an option: libc++ (AppleClang) does not implement the
+   C++17 special math functions.
+
+Keep `boost::math::sph_bessel` (also missing in libc++; our own version would
+need careful checking at large l and argument) and `program_options` (no
+standard equivalent).
+
+**2.3C: command-line behavior** (updates the CLI references on purpose)
+
+* B1: non-zero exit code on errors, together with `boost::diagnostic_information`
+  → `e.what()` (moved here from 2.2).
+* B16: accept `half-lorentzian` as well as `half lorentzian`, and show the real
+  `CUT` default in `--help.grids`.
+* B10: a `SEED` parameter for the bootstrap errors (moved here from 2.2).
+
+**2.3D: numerical changes** (each its own commit and `REFERENCE_CHANGES.md` entry)
+
+* B7b: time-bosonic kernel: use the ω→0 limit only at ω = 0.
+* B15: remove the `float` casts in the log grid.
+* B11 (optional): `BDCSVD` instead of `JacobiSVD`; measure against the
+  references before deciding.
+
+Keep the `MaxEntSimulation` public getters stable, because the tests use them.
 
 ### 2.4 Hygiene and CI
 
@@ -406,6 +439,17 @@ This is the last task of step 2, still on ALPSCore:
 2. **HDF5 seam:** move all archive access into `maxent_io_hdf5.{hpp,cpp}`, with functions such as `write_vector(path, Eigen::VectorXd)`, `read_vector(path)` and `write_params(...)`. Nothing else includes `<alps/hdf5...>`.
 3. Tests use `std::filesystem::temp_directory_path()` in place of `alps::temporary_filename`.
 
+### 2.6 Physics and algorithms (after 2.3–2.5)
+
+Moved out of 2.3 (2026-09-24) because they are physics or algorithm changes
+rather than cleanup; each needs its own decision:
+
+* B17: Legendre bosonic kernel (currently the fermionic formula).
+* B18: default grid when `OMEGA_MIN=0` (change the default, or warn).
+* B19: minimizer divergence (step-size control or a trust region).
+* B21: logging convention and verbosity (issue #46).
+* B22: performance of the core at large NFREQ.
+
 **Step 2 acceptance:**
 
 * Configures with CMake 3.22 through 4.x.
@@ -640,3 +684,4 @@ make -j8 && ctest     # 6/6 executables, 35 cases pass, ~9 s
 | 2026-09-23 | 2.x | Pade removed (D7), on branch `modernize/remove-pade`: `pade/` (13 files: 6 built, 5 unbuilt alternatives, header, CMake), `MAXENT_BUILD_PADE`, and the README section. It required GMP (`mpf_class`, 256-bit default precision) and did not compile (B20). |
 | 2026-09-24 | 2.2 | Step 2.2 done on `modernize/step2.2`: kk regression cases, kk spline, closed-form Legendre kernel (GSL removed), Boost utilities → std. No reference changed; 41/41 tests. ALPSCore#668 merged (GoogleTest 1.18, not installed, CMake 3.16). Found B22 (performance). |
 | 2026-09-24 | sync | Synced with ALPSCore `master` `560ae112` (#669: C++14 required and default, ALPSCore compiled with exactly `ALPS_CXX_STD`, Boost.Math dependency removed; #670: Doxygen removed). Maxent `d529a14` against it (installed as C++14, Boost 1.88): configures without `Boost_DIR`, compiles as C++17, 41/41 tests, regression unchanged (only the known Legendre rounding-level differences from 2.2). Note: build ALPSCore in a fresh build directory; an old one keeps a stale compiler identification that Maxent's configure then rejects. |
+| 2026-09-24 | 2.3A | Warning cleanup and no-output-change fixes complete: Eigen index types, initializer order, unused/deprecated code, inline HDF5/LAPACK helpers, unconditional size validation, unreachable kernel branches, robust three- and five-column text-file extraction, and dead bindings/config code removed. The dev preset enables `MAXENT_WERROR`; 44/44 Release tests pass including full regression, and the LAPACK build is warning-free. |
diff --git a/CMakePresets.json b/CMakePresets.json
index 890f970..8b4d5b7 100644
--- a/CMakePresets.json
+++ b/CMakePresets.json
@@ -20,7 +20,10 @@
       "name": "dev",
       "displayName": "Debug",
       "inherits": "base",
-      "cacheVariables": {"CMAKE_BUILD_TYPE": "Debug"}
+      "cacheVariables": {
+        "CMAKE_BUILD_TYPE": "Debug",
+        "MAXENT_WERROR": "ON"
+      }
     },
     {
       "name": "asan",
diff --git a/legendre_convert/legendre_convert.cpp b/legendre_convert/legendre_convert.cpp
index c4f9396..cedbf56 100644
--- a/legendre_convert/legendre_convert.cpp
+++ b/legendre_convert/legendre_convert.cpp
@@ -324,7 +324,7 @@ double GtBoot(vector_type gl_in, void *arg){
 //************************
 int main(int argc, char**argv){
 	namespace po = boost::program_options;
-	double beta,c1,c2,c3;
+	double beta,c1,c2;
 	int lmax,numConvergence,maxit,maxn;
 	bool backcontinue = false,enforceTailOff = false,enforceErrCheck=true,continueMatsubara=false;
 	std::string input_gtau_filename,gl_filename;
diff --git a/src/eigen_hdf5.hpp b/src/eigen_hdf5.hpp
index ad440c1..cdc0073 100644
--- a/src/eigen_hdf5.hpp
+++ b/src/eigen_hdf5.hpp
@@ -11,7 +11,6 @@
 
 #pragma once
 #include <alps/hdf5/archive.hpp>
-#include <alps/utilities/cast.hpp>
 #include <Eigen/Core>
 #include <algorithm>
 
@@ -32,7 +31,7 @@ namespace alps {
                       if (value.size()) {
                           std::vector<std::size_t> first(get_extent(value[0]));
                           if (true)
-                              for(std::size_t i=1;i<value.size();i++) {
+                              for(Eigen::Index i=1;i<value.size();i++) {
                                   std::vector<std::size_t> size(get_extent(value[i]));
                                   if (
                                          first.size() != size.size()
@@ -46,7 +45,7 @@ namespace alps {
                   }
           };
         }
-        void save(
+        inline void save(
             archive &ar
           , std::string const & path
           , Eigen::VectorXd &value
@@ -54,7 +53,6 @@ namespace alps {
           , std::vector<std::size_t> chunk =std::vector<std::size_t>()
           , std::vector<std::size_t> offset = std::vector<std::size_t>()
        ) {
-            using alps::cast;
             if (ar.is_group(path))
                 ar.delete_group(path);
             if (is_continuous<Eigen::VectorXd>::value && value.size() == 0)
diff --git a/src/eigen_lapack.hpp b/src/eigen_lapack.hpp
index 9d55c64..2755c34 100644
--- a/src/eigen_lapack.hpp
+++ b/src/eigen_lapack.hpp
@@ -21,9 +21,7 @@ extern "C" void dgesvd_( const char* jobu, const char* jobvt,
 ///performs a SVD on input matrix K
 /// returns K =  U (S) V^T where S is a vector of the diagonal values of
 /// matrix Sigma
-void lapack_svd(matrix_type &K, vector_type &S, matrix_type &Vt, matrix_type &U){
-
-   /*boost::numeric::bindings::lapack::gesvd('S', 'S', Kt, S, U_, Vt_);*/
+inline void lapack_svd(matrix_type &K, vector_type &S, matrix_type &Vt, matrix_type &U){
   //use 'S' for thin U,Vt matrices
   char jobu = 'S';
   char jobvt = 'S';
diff --git a/src/maxent.cpp b/src/maxent.cpp
index 3ce77f9..8d4639a 100644
--- a/src/maxent.cpp
+++ b/src/maxent.cpp
@@ -89,7 +89,7 @@ int main(int argc,const char** argv)
 
   std::string basename;
   if(parms.defaulted("BASENAME")){
-    basename = alps::fs::remove_extensions(parms.get_origin_name()) + ".out";
+    basename = alps::fs::remove_extensions(alps::origin_name(parms)) + ".out";
     parms["BASENAME"] = basename;
   }
   else
@@ -151,7 +151,7 @@ int main(int argc,const char** argv)
             ofstream_ spec_file;
             spec_file.open((basename+".varspec.dat").c_str());
             spec_file << "#omega mean_maxspec stdev_maxspec mean_avspec stdev_avspec" <<std::endl;
-            for (std::size_t  i=0; i<omega_grid.size(); ++i)
+            for (Eigen::Index i=0; i<omega_grid.size(); ++i)
               spec_file << omega_grid(i)
                 << " " << mean_max(i) << " " << stdev_max(i)
                 << " " << mean_av(i)  << " " << stdev_av(i) 
@@ -171,4 +171,3 @@ int main(int argc,const char** argv)
         std::cerr << "Caught Exception" << boost::current_exception_diagnostic_information();
     }
 }
-
diff --git a/src/maxent_backcont.cpp b/src/maxent_backcont.cpp
index f9dfd45..5b7ab21 100644
--- a/src/maxent_backcont.cpp
+++ b/src/maxent_backcont.cpp
@@ -39,12 +39,10 @@ vector_type Backcont::backcontinue(const vector_type &A){
 }
 double Backcont::max_error(const vector_type &y1, const vector_type &y2){
   double max_err = 0.0;
-#ifdef NDEBUG
   if(y1.size() != y2.size()){
     throw std::runtime_error("backcont vector size mismatch!");
   }
-#endif
-  for(int i=0;i<y1.size();i++){
+  for(Eigen::Index i=0;i<y1.size();i++){
     double delta = std::abs(y1(i)-y2(i));
     if(delta>max_err){
       max_err = delta;
diff --git a/src/maxent_backcont.hpp b/src/maxent_backcont.hpp
index 8acaafc..d3f5edf 100644
--- a/src/maxent_backcont.hpp
+++ b/src/maxent_backcont.hpp
@@ -11,7 +11,6 @@
 
 #pragma once
 #include "maxent.hpp"
-#include <alps/config.hpp> // needed to set up correct bindings
 
 ///This class contains necessities to analytically continue a 
 //real function to the imaginary axis
diff --git a/src/maxent_helper.cpp b/src/maxent_helper.cpp
index 4d7ca75..78f2624 100644
--- a/src/maxent_helper.cpp
+++ b/src/maxent_helper.cpp
@@ -10,7 +10,6 @@
 *****************************************************************************/
 
 #include "maxent.hpp"
-#include <alps/config.hpp> // needed to set up correct bindings
 #include <cmath>
 #include <ctime>
 #include <Eigen/Eigenvalues>
@@ -36,7 +35,7 @@ MaxEntParameters(p) , def_(nfreq()), text_output(p["TEXT_OUTPUT"])
 /// this is needed for transform_into_singular_space
 /// to work safely
 void MaxEntHelper::checkDefaultModel(const vector_type &D) const{
-    for(int i=0;i<D.size();i++){
+    for(Eigen::Index i=0;i<D.size();i++){
         double Di=D(i);
         if(Di==0 || std::isnan(Di))
           throw std::logic_error("Error: Default model point = 0 at omega="
@@ -49,7 +48,7 @@ void MaxEntHelper::checkDefaultModel(const vector_type &D) const{
 vector_type MaxEntHelper::transform_into_singular_space(vector_type A) const
 {
   double D;
-  for (unsigned int i=0; i<A.size(); ++i) {
+  for (Eigen::Index i=0; i<A.size(); ++i) {
     D=Default(i);
     A[i] /= D;
     A[i] = A[i]==0. ? 0. : log(A[i]);
@@ -61,7 +60,7 @@ vector_type MaxEntHelper::transform_into_singular_space(vector_type A) const
 vector_type MaxEntHelper::transform_into_real_space(vector_type u) const
 {
   u = maxent_prec_prod(Vt().transpose(), u);
-  for (unsigned int i=0; i<u.size(); ++i) {
+  for (Eigen::Index i=0; i<u.size(); ++i) {
     u[i] = exp(u[i]);
     u[i] *= Default(i);
   }
@@ -72,7 +71,7 @@ vector_type MaxEntHelper::transform_into_real_space(vector_type u) const
 vector_type MaxEntHelper::get_spectrum(const vector_type& u) const
 {
   vector_type A = transform_into_real_space(u);
-  for (unsigned int i=0; i<A.size(); ++i) 
+  for (Eigen::Index i=0; i<A.size(); ++i)
     A[i] /= delta_omega(i);
   return A;
 }
@@ -84,8 +83,8 @@ matrix_type MaxEntHelper::left_side(const vector_type& u) const
 {
   vector_type A = transform_into_real_space(u);
   matrix_type M = Vt().transpose();
-  for (unsigned int i=0; i<M.rows(); ++i) 
-    for (unsigned int j=0; j<M.cols(); ++j) 
+  for (Eigen::Index i=0; i<M.rows(); ++i)
+    for (Eigen::Index j=0; j<M.cols(); ++j)
       M(i,j) *= A[i];
   M = maxent_prec_prod(Vt(), M);
   M = maxent_prec_prod(Sigma() ,M);
@@ -111,8 +110,8 @@ double MaxEntHelper::step_length(const vector_type& delta, const vector_type& u)
 {
   vector_type A = transform_into_real_space(u);
   matrix_type L = Vt().transpose();
-  for (unsigned int i=0; i<L.rows(); ++i) 
-    for (unsigned int j=0; j<L.cols(); ++j) 
+  for (Eigen::Index i=0; i<L.rows(); ++i)
+    for (Eigen::Index j=0; j<L.cols(); ++j)
       L(i,j) *= A[i];
   L = maxent_prec_prod(Vt(), L);
   return delta.dot(maxent_prec_prod(L, delta));
@@ -121,11 +120,10 @@ double MaxEntHelper::step_length(const vector_type& delta, const vector_type& u)
 //Bryan's paper section 2.3 (or after eq 22)
 double MaxEntHelper::convergence(const vector_type& u, const double alpha) const 
 {
-  //using namespace boost::numeric::ublas;
   vector_type A = transform_into_real_space(u);
   matrix_type L = Vt().transpose();
-  for (unsigned int i=0; i<L.rows(); ++i) 
-    for (unsigned int j=0; j<L.cols(); ++j) 
+  for (Eigen::Index i=0; i<L.rows(); ++i)
+    for (Eigen::Index j=0; j<L.cols(); ++j)
       L(i,j) *= A[i];
   L = maxent_prec_prod(Vt(), L);
   vector_type alpha_dSdu = -alpha*maxent_prec_prod(L, u);
@@ -140,38 +138,35 @@ double MaxEntHelper::log_prob(const vector_type& u, const double alpha) const
 {
   matrix_type L = maxent_prec_prod_trans(K(), K());
   const vector_type A = transform_into_real_space(u);
-  for (unsigned int i=0; i<L.rows(); ++i)
-    for (unsigned int j=0; j<L.cols(); ++j)
+  for (Eigen::Index i=0; i<L.rows(); ++i)
+    for (Eigen::Index j=0; j<L.cols(); ++j)
       L(i,j) *= sqrt(A[i])*sqrt(A[j]);
-  for (unsigned int i=0; i<L.rows(); ++i)
+  for (Eigen::Index i=0; i<L.rows(); ++i)
     L(i,i) += alpha;
-  //boost::numeric::bindings::lapack::potrf(boost::numeric::bindings::lower(L));
-  
   //LAPACK does potrf in place, while Eigen does not
   Eigen::LLT<matrix_type,Eigen::Lower> lltofL(L);
   L=lltofL.matrixL();
 
   double log_det = 0.;
-  for (unsigned int i=0; i<L.rows(); ++i) 
+  for (Eigen::Index i=0; i<L.rows(); ++i)
     log_det  += log(L(i,i)*L(i,i));
   return 0.5*( (nfreq())*log(alpha) - log_det ) - Q(u, alpha);
 }
 
 double MaxEntHelper::chi_scale_factor(vector_type A, const double chi_sq, const double alpha) const
 {
-  for (unsigned int i=0; i<A.size(); ++i) 
+  for (Eigen::Index i=0; i<A.size(); ++i)
     A[i] *= delta_omega(i);
 
   matrix_type L = maxent_prec_prod_trans(K(), K());
-  for (unsigned int i=0; i<L.rows(); ++i)
-    for (unsigned int j=0; j<L.cols(); ++j)
+  for (Eigen::Index i=0; i<L.rows(); ++i)
+    for (Eigen::Index j=0; j<L.cols(); ++j)
       L(i,j) *= sqrt(A[i])*sqrt(A[j]);
   vector_type lambda(L.rows());
-  //bindings::lapack::syev('N', bindings::upper(L) , lambda, bindings::lapack::optimal_workspace());
   Eigen::SelfAdjointEigenSolver<matrix_type> es(L);
   lambda = es.eigenvalues();
   double Ng = 0.;
-  for (unsigned int i=0; i<lambda.size(); ++i) {
+  for (Eigen::Index i=0; i<lambda.size(); ++i) {
     if (lambda[i]>=0) 
       Ng += lambda[i]/(lambda[i]+alpha);
   }
@@ -186,12 +181,12 @@ double MaxEntHelper::chi2(const vector_type& A) const
   vector_type del_G = maxent_prec_prod(K(), A) - y();
   
   /*std::cout<<"in computation of chi2:"<<std::endl;
-   for(int i=0;i<y().size();++i){
+   for(Eigen::Index i=0;i<y().size();++i){
    std::cout<<i<<" "<<maxent_prec_prod(K(), A)[i]<<" "<<y()[i]<<std::endl;
    }*/
   
   double c = 0;
-  for (unsigned int i=0; i<del_G.size(); ++i) 
+  for (Eigen::Index i=0; i<del_G.size(); ++i)
     c += del_G[i]*del_G[i];
   return c;
 }
@@ -201,7 +196,7 @@ void MaxEntHelper::print_chi2(const vector_type& A, std::ostream &os) const
   vector_type backcont=maxent_prec_prod(K(), A);
   vector_type defaultm=maxent_prec_prod(K(), Default());
   os<<"#first column: index (Matsubara frequency). second column: fitted function. third: input data. fourth: default model."<<std::endl;
-  for(int i=0;i<y().size();++i){
+  for(Eigen::Index i=0;i<y().size();++i){
     os<<i<<" "<<backcont[i]<<" "<<y()[i]<<" "<<defaultm[i]<<std::endl;
   }
   os<<std::endl;
@@ -212,7 +207,7 @@ void MaxEntHelper::print_chi2(const vector_type& A, std::ostream &os) const
 double MaxEntHelper::entropy(const vector_type& A) const 
 {
   double S = 0;
-  for (unsigned int i=0; i<A.size(); ++i) {
+  for (Eigen::Index i=0; i<A.size(); ++i) {
     double lg = A[i]==0. ? 0. : log(A[i]/Default(i));
     S += A[i] - Default(i) - A[i]*lg;
   }
@@ -246,7 +241,6 @@ void MaxEntHelper::backcontinue(ofstream_ &os, const vector_type &A_in,const dou
     vector_type G = bc.backcontinue(A);
     kernel_type k_type = pp->getKernelType();
     
-    double beta = 1/(pp->T());
     bool ph_sym = true;
     if(k_type == frequency_fermionic_kernel ||
        k_type == frequency_bosonic_kernel   ||
@@ -256,26 +250,26 @@ void MaxEntHelper::backcontinue(ofstream_ &os, const vector_type &A_in,const dou
     ext_back = G*norm;
     if(text_output){
       if(ph_sym){
-        for(int n=0; n<G.size();n++){
+        for(Eigen::Index n=0; n<G.size();n++){
           os << pp->inputGrid(n) << " " << G(n)*norm << std::endl;
         }
       }
       else{
-        for(int n=0;n<G.size();n+=2){
+        for(Eigen::Index n=0;n<G.size();n+=2){
           os << pp->inputGrid(n/2) << " " << G(n)*norm << " " << G(n+1)*norm << std::endl;
         }
       }
     }
     //scale y by error then determine 'error' of integral
     vector_type y_scaled = y();
-    for(int i=0;i<y_scaled.size();i++){
+    for(Eigen::Index i=0;i<y_scaled.size();i++){
       y_scaled(i) *= sigma(i);
     }
     double max_err = bc.max_error(G,y_scaled); 
 
     //A is missing delta_omega value, add back in for chi2
     vector_type A_chi = A;
-    for (unsigned int i=0; i<A_chi.size(); ++i) 
+    for (Eigen::Index i=0; i<A_chi.size(); ++i)
       A_chi[i] *= delta_omega(i);
 
     double chi_sq = chi2(A_chi);
@@ -297,10 +291,10 @@ void determineVariance(std::vector<vector_type> &in,vector_type &mean, vector_ty
   }
   mean /= in.size();
   //compute stddev  
-  for(int i=0;i<mean.size();i++){
+  for(Eigen::Index i=0;i<mean.size();i++){
     double stddev =0;
     double mean_i = mean(i);
-    for(int v=0;v<in.size();v++){
+    for(std::size_t v=0;v<in.size();v++){
       double val = in[v](i);
       stddev+= (val-mean_i)*(val-mean_i);
     }
@@ -357,7 +351,7 @@ void MaxEntHelper::generateCovariantErr(const vector_type& A, const double alpha
 
     //setup and transform sqrt(A)
     vector_type A_u = A;
-    for (int i=0;i<A.size();i++)
+    for (Eigen::Index i=0;i<A.size();i++)
       A_u(i) = sqrt(A(i));
     A_u = u*A_u;
 
@@ -380,7 +374,7 @@ void MaxEntHelper::generateCovariantErr(const vector_type& A, const double alpha
 
     //save file
     os << "#omega A A_mean approx_err" <<std::endl;
-    for (std::size_t  i=0; i<A.size(); ++i){
+    for (Eigen::Index i=0; i<A.size(); ++i){
         os << omega_coord(i) << " " << A[i] << " " << mean(i)*mean(i) << " " <<A[i]*2*std_err(i)<< std::endl;
     }
   }
diff --git a/src/maxent_kernel.cpp b/src/maxent_kernel.cpp
index cf0bf32..b5ec5c6 100644
--- a/src/maxent_kernel.cpp
+++ b/src/maxent_kernel.cpp
@@ -87,8 +87,6 @@ K_(ndat_,nfreq_)
         }
       }
     }
-    else if(ktype_== time_fermionic_legendre_kernel || ktype_==time_bosonic_legendre_kernel)
-        setup_legendre_kernel(p,freq,ndat_);
     //for zero temperature, only positive frequency matters
     else if (ktype_ == time_boris_kernel) {
       for (int i=0; i<ndat_; ++i) {
@@ -99,15 +97,6 @@ K_(ndat_,nfreq_)
         }
       }
     }
-    else if(ktype_==time_fermionic_kernel){
-        for (int i=0; i<ndat_; ++i) {
-						double tau=tau_points_[i];
-            for (int j=0; j<nfreq_; ++j) {
-                double omega = freq[j];
-                K_(i,j) =  -1.;
-            }
-        }
-    }
     else if(ktype_==frequency_fermionic_ph_kernel) {
     for (int i=0; i<ndat_; ++i) {
       double omegan = (2*i+1)*M_PI*T_;
@@ -201,15 +190,9 @@ void kernel::set_kernel_type(const std::string &dataspace_name, const std::strin
 
   if(dtype_==time_dataspace){
     if(kernel_name=="fermionic")
-            if(dtype_==legendre_dataspace)
-                ktype_=time_fermionic_legendre_kernel;
-            else
-                ktype_=time_fermionic_kernel;
+      ktype_=time_fermionic_kernel;
     else if(kernel_name=="bosonic")
-            if(dtype_==legendre_dataspace)
-                ktype_=time_bosonic_legendre_kernel;
-            else
-                ktype_=time_bosonic_kernel;
+      ktype_=time_bosonic_kernel;
     else if(kernel_name=="tzero")
       ktype_=time_boris_kernel;
     else throw std::invalid_argument("unknown kernel name. In the time domain it should be fermionic, bosonic, or tzero.");
diff --git a/src/maxent_matrix_def.hpp b/src/maxent_matrix_def.hpp
index d3003a1..77a242e 100644
--- a/src/maxent_matrix_def.hpp
+++ b/src/maxent_matrix_def.hpp
@@ -11,7 +11,6 @@
 
 #pragma once
 
-#include <alps/config.hpp> // needed to set up correct bindings
 #include "maxent_config.hpp"
 #include <Eigen/Core>
 
diff --git a/src/maxent_params.cpp b/src/maxent_params.cpp
index bdfb493..5f61e88 100644
--- a/src/maxent_params.cpp
+++ b/src/maxent_params.cpp
@@ -11,7 +11,6 @@
 
 #include "maxent.hpp"
 #include "eigen_lapack.hpp"
-#include <alps/config.hpp> // needed to set up correct bindings
 #include <Eigen/SVD>
 #include <Eigen/Eigenvalues>
 #include <alps/hdf5/vector.hpp>
@@ -51,9 +50,12 @@ void ContiParameters::read_data_from_text_file(const alps::params& p) {
   std::string dataspace = p["DATASPACE"].as<std::string>();
   to_lower(dataspace);
   if(dataspace == "time" || dataspace == "legendre" || p["PARTICLE_HOLE_SYMMETRY"]==true){
-    while (datstream) {
-      double index, X_i, dX_i;
-      datstream >> index >> X_i >> dX_i;
+    double index, X_i, dX_i;
+    while (true) {
+      datstream >> std::ws;
+      if (datstream.eof()) break;
+      if (!(datstream >> index >> X_i >> dX_i))
+        throw std::runtime_error("malformed data record in: " + fname);
       if (datIn < ndat()) {
         inputGrid_(datIn) = index;
         y_(datIn) = X_i / static_cast<double>(p["NORM"]);
@@ -62,7 +64,6 @@ void ContiParameters::read_data_from_text_file(const alps::params& p) {
       }
       expectedDatIn++;
     }
-    expectedDatIn-=1;
   }
   else{
     if(ndat()%2 != 0){
@@ -71,9 +72,12 @@ void ContiParameters::read_data_from_text_file(const alps::params& p) {
                   << std::endl;
         throw std::runtime_error("Your NDAT is odd!");\
     }
-    while (datstream) {
-      double index, X_i_re, dX_i_re, X_i_im, dX_i_im;
-      datstream >> index >> X_i_re >> dX_i_re >> X_i_im >> dX_i_im;
+    double index, X_i_re, dX_i_re, X_i_im, dX_i_im;
+    while (true) {
+      datstream >> std::ws;
+      if (datstream.eof()) break;
+      if (!(datstream >> index >> X_i_re >> dX_i_re >> X_i_im >> dX_i_im))
+        throw std::runtime_error("malformed complex data record in: " + fname);
       if (datIn < ndat()) {
         inputGrid_(datIn) = index;
         inputGrid_(datIn+1) = index; 
@@ -83,11 +87,8 @@ void ContiParameters::read_data_from_text_file(const alps::params& p) {
         sigma_(datIn+1) = dX_i_im / static_cast<double>(p["NORM"]);
         datIn+=2;
       }
-      expectedDatIn++;
+      expectedDatIn+=2;
     }
-    //fix for matsubara length
-    expectedDatIn*=2;
-    expectedDatIn-=1;
   }
   if(p["COVARIANCE_MATRIX"]!="") {
     std::string fname = p["COVARIANCE_MATRIX"];
@@ -117,13 +118,13 @@ void ContiParameters::read_data_from_hdf5_file(const alps::params& p) {
   std::stringstream path;
   path << "/Data";
   ar >> alps::make_pvp(path.str(), tmp);
-  for (std::size_t i = 0; i < ndat(); i++)
+  for (int i = 0; i < ndat(); i++)
     y_(i) = tmp[i] / static_cast<double>(p["NORM"]);
   path.str("");
   if (p["COVARIANCE_MATRIX"]=="") {
     path << "/Error";
     ar >> alps::make_pvp(path.str(), tmp);
-    for (std::size_t i = 0; i < ndat(); i++)
+    for (int i = 0; i < ndat(); i++)
       sigma_(i) = tmp[i] / static_cast<double>(p["NORM"]);
   } else {
     path << "/Covariance";
@@ -131,8 +132,8 @@ void ContiParameters::read_data_from_hdf5_file(const alps::params& p) {
     tmp.clear();
     tmp.resize(ndat() * ndat());
     ar >> alps::make_pvp(path.str(), tmp);
-    for (std::size_t i = 0; i < ndat(); i++)
-      for (std::size_t j = 0; j < ndat(); j++)
+    for (int i = 0; i < ndat(); i++)
+      for (int j = 0; j < ndat(); j++)
         cov_(i, j) = tmp[i * ndat() + j];
   }
 }
@@ -225,7 +226,6 @@ void ContiParameters::read_covariance_matrix_from_text_file(
 void ContiParameters::decompose_covariance_matrix(const alps::params& p){
 
     vector_type var(ndat());
-    //bindings::lapack::syev('V', bindings::upper(cov_) , var, bindings::lapack::optimal_workspace()); 
     //TODO: check if this truly implements lapack's expected overwrite of cov_
     Eigen::SelfAdjointEigenSolver<matrix_type> es(cov_);
     var=es.eigenvalues();
@@ -292,17 +292,12 @@ void MaxEntParameters::truncate_to_singular_space(const vector_type& S) {
 
 void MaxEntParameters::singular_value_decompose_kernel(bool verbose,
     vector_type& S) {
-  /*boost::numeric::bindings::lapack::gesvd('S', 'S', Kt, S, U_, Vt_);
-  */
-
 #ifdef HAVE_LAPACK
 
   matrix_type Kt = K_; // gesvd destroys K!
   lapack_svd(Kt,S, Vt_,U_);
 
 #else
-  const double threshold = std::sqrt(std::numeric_limits<double>::epsilon())
-      * nfreq();
   Eigen::JacobiSVD<matrix_type> svd(K_,Eigen::ComputeThinU | Eigen::ComputeThinV);
   //svd.setThreshold(threshold);
   S=svd.singularValues();
@@ -318,7 +313,7 @@ void MaxEntParameters::singular_value_decompose_kernel(bool verbose,
     std::cout << "# eps = " << sqrt(std::numeric_limits<double>::epsilon())
         << std::endl << "# prec = " << prec << std::endl;
 
-  for (unsigned int s = 0; s < S.size(); ++s) {
+  for (Eigen::Index s = 0; s < S.size(); ++s) {
     if (verbose)
       std::cout << "# " << s << "\t" << S[s] << "\n";
 
@@ -387,10 +382,12 @@ MaxEntParameters::MaxEntParameters(alps::params& p) :
   k_type = ker.getKernelType();
 
   //scale lhs and rhs according to errors, etc.
-  if (p["COVARIANCE_MATRIX"]!="")
+  if (p["COVARIANCE_MATRIX"]!="") {
     decompose_covariance_matrix(p);
-    
-    check_high_frequency_limit(y(),k_type);
+  }
+
+  // This diagnostic applies to both diagonal errors and covariance input.
+  check_high_frequency_limit(y(),k_type);
 
   //Look around Eq. D.5 in Sebastian's thesis. We have sigma_ = sqrt(eigenvalues of covariance matrix) or, 
   //in case of a diagonal covariance matrix, we have sigma_=SIGMA_X. 
@@ -411,6 +408,3 @@ MaxEntParameters::MaxEntParameters(alps::params& p) :
   compute_minimal_chi2();
 }
 
-
-
-
diff --git a/src/maxent_simulation.cpp b/src/maxent_simulation.cpp
index fe4ead8..5282782 100644
--- a/src/maxent_simulation.cpp
+++ b/src/maxent_simulation.cpp
@@ -10,7 +10,6 @@
 *****************************************************************************/
 
 #include "maxent.hpp"
-#include <alps/config.hpp> // needed to set up correct bindings
 #include <alps/hdf5/vector.hpp>
 #include <cmath>
 #include <Eigen/LU>
@@ -27,8 +26,8 @@ MaxEntSimulation::MaxEntSimulation(alps::params &parms)
 , self(parms["SELF"])
 , make_back(parms["BACKCONTINUE"])
 , gen_err(parms["GENERATE_ERR"])
-, qvec((int)parms["N_ALPHA"])
 , nfreq(parms["NFREQ"].as<int>())
+, qvec((int)parms["N_ALPHA"])
 {
   std::string bn=parms["BASENAME"]; name=bn+'.';
 
@@ -36,7 +35,7 @@ MaxEntSimulation::MaxEntSimulation(alps::params &parms)
   const double alpha_min = parms["ALPHA_MIN"];                                          //Smallest value of \alpha that is tried
   const double alpha_max = parms["ALPHA_MAX"];                                          //Largest  value of \alpha that is tried
   alpha[0] = alpha_max;
-  for (std::size_t a=1; a<alpha.size(); ++a)                                            //These are all the alpa values on a log grid
+  for (Eigen::Index a=1; a<alpha.size(); ++a)                                            //These are all the alpa values on a log grid
     alpha[a] =  alpha[a-1] * std::pow(alpha_min/alpha_max, 1./double(alpha.size()-1));
 }
 ///define parameter defaults
@@ -123,7 +122,7 @@ void MaxEntSimulation::run()
   }
   //this loop is the 'core' of the maxent program: iterate over all alphas, compute the spectra, normalization, and probabilities
   //loop over all alpha values
-  for (std::size_t a=0; a<alpha.size(); ++a) {
+  for (Eigen::Index a=0; a<alpha.size(); ++a) {
     std::cerr << "alpha it: " << a << "\t";
     //fitting procedure for 'u'
     u = levenberg_marquardt(u, alpha[a]);
@@ -133,7 +132,7 @@ void MaxEntSimulation::run()
     std::cerr << "norm: " << transform_into_real_space(u).sum() << "\t";
     if (text_output) {
       spectral_function_file<<"# alpha: "<<alpha[a]<<std::endl;
-      for (std::size_t i=0; i<A.size(); ++i)
+      for (Eigen::Index i=0; i<A.size(); ++i)
         spectral_function_file << omega_coord(i) << " " << A[i] << "\n";
       spectral_function_file << "\n";
     }
@@ -149,7 +148,7 @@ void MaxEntSimulation::run()
     qvec(a)=Q(u,alpha[a]); 
   }
     omegaGrid.resize(nfreq);
-    for(std::size_t i=0;i<nfreq;i++)
+    for(Eigen::Index i=0;i<omegaGrid.size();i++)
 	    omegaGrid(i)=omega_coord(i);
 }
   //everything from here on down is evaluation.
@@ -157,13 +156,13 @@ void MaxEntSimulation::evaluate(){
   if (text_output) {
     ofstream_ chi_squared_file;
     chi_squared_file.open((name+"chi2.dat").c_str());
-    for (std::size_t a=0; a<chi_sq.size(); ++a){
+    for (Eigen::Index a=0; a<chi_sq.size(); ++a){
       chi_squared_file << alpha[a] << " " << chi_sq[a] << std::endl;
     }
   }
   int a_chi = 0;
   double diff = std::abs(chi_sq[0]-ndat());
-  for (std::size_t a=1; a<chi_sq.size(); ++a) {
+  for (Eigen::Index a=1; a<chi_sq.size(); ++a) {
     double diff_new = std::abs(chi_sq[a]-ndat());
     if (diff_new < diff) {
       diff = diff_new;
@@ -176,11 +175,10 @@ void MaxEntSimulation::evaluate(){
   if (text_output){
     ofstream_ chispec_file;
     chispec_file.open((name+"chispec.dat").c_str());
-    for (std::size_t i=0; i<spectra[0].size(); ++i){
+    for (Eigen::Index i=0; i<spectra[0].size(); ++i){
       chispec_file << omega_coord(i) << " " << spectra[a_chi][i]*norm << " " << def[i]*norm << std::endl;
     }
   }
-  //boost::numeric::ublas::vector<double>::const_iterator max_lprob = std::max_element(lprob.begin(), lprob.end());  
   //const int max_a = max_lprob-lprob.begin();
   int max_a,nothing; double max_lprob;
   max_lprob=lprob.maxCoeff(&max_a,&nothing);
@@ -191,14 +189,14 @@ void MaxEntSimulation::evaluate(){
 	ar << alps::make_pvp("/alpha/values",alpha);
 
   vector_type om(spectra[0].size());
-  for (int i=0;i<om.size();i++) om[i] = omega_coord(i);        
+  for (Eigen::Index i=0;i<om.size();i++) om[i] = omega_coord(i);
   ar<<alps::make_pvp("/spectrum/omega",om);
 
   //output 'maximum' spectral function (classical maxent metod)
   if (text_output){
     ofstream_ maxspec_file;
     maxspec_file.open((name+"maxspec.dat").c_str());
-    for (std::size_t i=0; i<spectra[0].size(); ++i)
+    for (Eigen::Index i=0; i<spectra[0].size(); ++i)
       maxspec_file << omega_coord(i) << " " << spectra[max_a][i]*norm << " " << def[i]*norm << std::endl;
   }
 	
@@ -208,37 +206,37 @@ void MaxEntSimulation::evaluate(){
   ar << alps::make_pvp("/spectrum/maximum",maxspec);
   
 	vector_type prob(lprob.size());
-  for (std::size_t a=0; a<prob.size(); ++a) 
+  for (Eigen::Index a=0; a<prob.size(); ++a)
     prob[a] = exp(lprob[a]-max_lprob);
   double probnorm = 0;
-  for (std::size_t a=0; a<prob.size()-1; ++a) 
+  for (Eigen::Index a=0; a<prob.size()-1; ++a)
     probnorm += 0.5*(prob[a]+prob[a+1])*(alpha[a]-alpha[a+1]);
   prob /= probnorm;
   ar << alps::make_pvp("/alpha/probability",prob);
   if (text_output){
     ofstream_ prob_str;
     prob_str.open((name+"prob.dat").c_str());
-    for (std::size_t a=0; a<prob.size(); ++a) {
+    for (Eigen::Index a=0; a<prob.size(); ++a) {
       prob_str << alpha[a] << "\t" << prob[a] << "\n";
     }
   }
   postprobdef = 0;
-  for (std::size_t a=0; a<lprob.size()-1; ++a) 
+  for (Eigen::Index a=0; a<lprob.size()-1; ++a)
     postprobdef += 0.5*(exp(lprob[a])+exp(lprob[a+1]))*(alpha[a]-alpha[a+1]);
   std::cout << "posterior probability of the default model: " << postprobdef << std::endl;
 
   //compute 'average' spectral function (Brian's method)
   avspec.resize(spectra[0].size());
-  for (std::size_t i=0; i<avspec.size(); ++i) {
+  for (Eigen::Index i=0; i<avspec.size(); ++i) {
     avspec[i] = 0.;
-    for (std::size_t a=0; a<prob.size()-1; ++a) 
+    for (Eigen::Index a=0; a<prob.size()-1; ++a)
       avspec[i] += 0.5*(prob[a]*spectra[a][i] +prob[a+1]*spectra[a+1][i])*(alpha[a]-alpha[a+1]);
   }
   //Estimate the variance for the spectrum
   vector_type varspec(spectra[0].size());
-  for (std::size_t i=0; i<varspec.size(); ++i) {
+  for (Eigen::Index i=0; i<varspec.size(); ++i) {
     varspec[i] = 0.;
-    for (std::size_t a=0; a<prob.size()-1; ++a)
+    for (Eigen::Index a=0; a<prob.size()-1; ++a)
       varspec[i] += 0.5*(prob[a]*(spectra[a][i]-avspec[i])*(spectra[a][i]-avspec[i]) + prob[a+1]*(spectra[a+1][i]-avspec[i])*(spectra[a+1][i]-avspec[i]))*(alpha[a]-alpha[a+1]);
   }
   avspec *= norm;
@@ -247,7 +245,7 @@ void MaxEntSimulation::evaluate(){
   if (text_output){
     ofstream_ avspec_file;
     avspec_file.open((name+"avspec.dat").c_str());
-    for (std::size_t  i=0; i<avspec.size(); ++i)
+    for (Eigen::Index i=0; i<avspec.size(); ++i)
       avspec_file << omega_coord(i) << " " << avspec[i] << " " << def[i]*norm << std::endl;
   }
   ar << alps::make_pvp("/spectrum/average",avspec);
@@ -265,13 +263,13 @@ void MaxEntSimulation::evaluate(){
     ofstream_ maxspec_anom_str;maxspec_anom_str.open((name+"maxspec_anom.dat").c_str());
     ofstream_ avspec_anom_str; avspec_anom_str.open((name+"avspec_anom.dat").c_str());
     vector_type spec(avspec.size());
-    for (std::size_t  i=0; i<avspec.size(); ++i){ 
+    for (Eigen::Index i=0; i<avspec.size(); ++i){
       //if(omega_coord(i)>=0.)
       spec[i] = avspec[i]*omega_coord(i)*M_PI;
       avspec_anom_str << omega_coord(i) << " " << spec[i]<<std::endl;
     }
     ar << alps::make_pvp("/spectrum/anomalous/average",spec);
-    for (std::size_t i=0; i<spectra[0].size(); ++i){
+    for (Eigen::Index i=0; i<spectra[0].size(); ++i){
       //if(omega_coord(i)>=0.)
       spec[i] = spectra[max_a][i]*norm*omega_coord(i)*M_PI;
       maxspec_anom_str << omega_coord(i) << " " << spec[i] << std::endl;
@@ -280,24 +278,24 @@ void MaxEntSimulation::evaluate(){
   }
   if(Kernel_type=="bosonic"){ //for the anomalous function: use A(Omega_)=Im chi(Omega_)/(pi Omega_) (as for anomalous)
     vector_type spec(avspec.size());
-    for (std::size_t  i=0; i<avspec.size(); ++i){
+    for (Eigen::Index i=0; i<avspec.size(); ++i){
       spec[i] = avspec[i]*omega_coord(i);
     }
     if (text_output) {
       ofstream_ avspec_anom_str;avspec_anom_str.open((name+"maxspec_bose.dat").c_str());
-      for (std::size_t  i=0; i<avspec.size(); ++i){
+      for (Eigen::Index i=0; i<avspec.size(); ++i){
         //if(omega_coord(i)>=0.)
         avspec_anom_str << omega_coord(i) << " " << spec[i]<<std::endl;
       }
     }
     ar << alps::make_pvp("/spectrum/bosonic/average",spec);
-    for (std::size_t i=0; i<spectra[0].size(); ++i){
+    for (Eigen::Index i=0; i<spectra[0].size(); ++i){
       //if(omega_coord(i)>=0.)
       spec[i] = spectra[max_a][i]*norm*omega_coord(i);
     }
     if (text_output) {
       ofstream_ maxspec_anom_str;maxspec_anom_str.open((name+"avspec_bose.dat").c_str());
-      for (std::size_t i=0; i<spectra[0].size(); ++i){
+      for (Eigen::Index i=0; i<spectra[0].size(); ++i){
         maxspec_anom_str << omega_coord(i) << " " << spec[i] << std::endl;
       }
     }
@@ -316,13 +314,13 @@ void MaxEntSimulation::evaluate(){
     ofstream_ maxspec_self_str;maxspec_self_str.open((name+"maxspec_self.dat").c_str());
     ofstream_ avspec_self_str; avspec_self_str.open((name+"avspec_self.dat").c_str());
     ofstream_ chispec_self_str; chispec_self_str.open((name+"chispec_self.dat").c_str());
-    for (std::size_t  i=0; i<avspec.size(); ++i){ 
+    for (Eigen::Index i=0; i<avspec.size(); ++i){
       avspec_self_str << omega_coord(i) << " " << -avspec[i]*M_PI<< " " << -def[i]*norm*M_PI<<std::endl;
     }
-    for (std::size_t i=0; i<spectra[0].size(); ++i){
+    for (Eigen::Index i=0; i<spectra[0].size(); ++i){
       maxspec_self_str << omega_coord(i) << " " << -spectra[max_a][i]*norm*M_PI<< " " << -def[i]*norm*M_PI << std::endl;
     }
-    for (std::size_t  i=0; i<specchi.size(); ++i){ 
+    for (Eigen::Index i=0; i<specchi.size(); ++i){
       chispec_self_str << omega_coord(i) << " " << -specchi[i]*M_PI<< " " << -def[i]*norm*M_PI<<std::endl;
     }
     //for public facing variables
@@ -374,11 +372,6 @@ vector_type MaxEntSimulation::levenberg_marquardt(vector_type u, const double al
     for (it2=0; it2<max_it; ++it2) {
       //compute change vector delta to u
       delta = iteration(u, alpha, mu);
-      /*std::cout<<"delta is: "<<delta<<std::endl;
-      vector_type z=transform_into_real_space(delta);
-      for(int i=0;i<z.size();++i){
-        std::cout<<omega_coord(i)<<" "<<z(i)<<std::endl;
-      }*/
       //compute Q = 1/2 chi^2 - \alpha S
       Q1 = Q(u+delta, alpha);
       if (step_length(delta, u)<=0.02) {
@@ -407,13 +400,12 @@ vector_type MaxEntSimulation::levenberg_marquardt(vector_type u, const double al
 vector_type MaxEntSimulation::iteration(vector_type u, const double alpha, const double mu) const
 {
   matrix_type M = left_side(u);
-  for (std::size_t i=0; i<M.rows(); ++i) 
+  for (Eigen::Index i=0; i<M.rows(); ++i)
     M(i,i) += alpha + mu;
   vector_type b = right_side(u) + alpha*u;
   matrix_type B(b.size(),1);
-  for (std::size_t i=0; i<M.rows(); ++i) 
+  for (Eigen::Index i=0; i<M.rows(); ++i)
     B(i,0) = -b[i];
-   //bindings::lapack::gesv(M, ipiv, B);
   //NOTE: gesv uses LU decomp, but we can switch to a safe QR routine as well
   matrix_type Bp = M.lu().solve(B);
   //may need a transposeInPlace();
diff --git a/test/backcontTest.cpp b/test/backcontTest.cpp
index 6cfe5a0..cc35149 100644
--- a/test/backcontTest.cpp
+++ b/test/backcontTest.cpp
@@ -210,6 +210,7 @@ alps::params p;
 
   //sensibility check
   EXPECT_EQ(G.size(),ndat);
+  EXPECT_THROW(bc.max_error(G, vector_type::Zero(ndat + 1)), std::runtime_error);
 
   //check that backcont is right to first order
   const double threshold = 0.1;
diff --git a/test/paramsTest.cpp b/test/paramsTest.cpp
index b361971..a7ec85a 100644
--- a/test/paramsTest.cpp
+++ b/test/paramsTest.cpp
@@ -101,6 +101,73 @@ std::string pf=alps::temporary_filename("in_file.dat");
   std::remove(pf.c_str());
 }
 
+TEST(Parameters,DataInFileWithTrailingBlankLines){
+  std::string pf=alps::temporary_filename("in_file.dat");
+  write_minimal_input_file(pf);
+  {
+    std::ofstream tempfile(pf.c_str(), std::ios::app);
+    tempfile << "\n\n";
+  }
+
+  alps::params p;
+  MaxEntSimulation::define_parameters(p);
+  p["BETA"]=2;
+  p["DATA"]=pf;
+  p["NDAT"] = 5;
+
+  ContiParameters c(p);
+  EXPECT_EQ(c.ndat(),5);
+  for(int i=0;i<c.ndat();i++){
+    EXPECT_NEAR(c.y(i),(i+1)*0.1,1e-10);
+    EXPECT_EQ(c.sigma(i),0.5);
+  }
+
+  std::remove(pf.c_str());
+}
+
+TEST(Parameters,RejectsIncompleteFinalDataRecord){
+  std::string pf=alps::temporary_filename("in_file.dat");
+  write_minimal_input_file(pf);
+  {
+    std::ofstream tempfile(pf.c_str(), std::ios::app);
+    tempfile << "5 0.6\n"; // Missing the error-bar column.
+  }
+
+  alps::params p;
+  MaxEntSimulation::define_parameters(p);
+  p["BETA"]=2;
+  p["DATA"]=pf;
+  // The five complete records already satisfy NDAT. The malformed tail must
+  // still be rejected rather than mistaken for clean EOF.
+  p["NDAT"] = 5;
+
+  EXPECT_THROW(ContiParameters c(p), std::runtime_error);
+  std::remove(pf.c_str());
+}
+
+TEST(Parameters,RejectsIncompleteFinalComplexDataRecord){
+  std::string pf=alps::temporary_filename("in_file.dat");
+  {
+    std::ofstream tempfile(pf.c_str());
+    tempfile << "0 0.1 0.01 -0.2 0.02\n"
+             << "1 0.3 0.03 -0.4 0.04\n"
+             << "2 0.5 0.05 -0.6\n"; // Missing the imaginary error bar.
+  }
+
+  alps::params p;
+  MaxEntSimulation::define_parameters(p);
+  p["BETA"]=2;
+  p["DATA"]=pf;
+  p["DATASPACE"]="frequency";
+  p["PARTICLE_HOLE_SYMMETRY"]=false;
+  // The two complete complex records already satisfy NDAT. The malformed tail
+  // must still be rejected rather than mistaken for clean EOF.
+  p["NDAT"] = 4;
+
+  EXPECT_THROW(ContiParameters c(p), std::runtime_error);
+  std::remove(pf.c_str());
+}
+
 TEST(Parameters,MaxentParams){
     //set up parameters
 	alps::params p;
@@ -181,6 +248,7 @@ TEST(Parameters,HighFrequencyCheck){
     p["SHIFT"]=2.0;
     
     std::complex<double> G;
+    testing::internal::CaptureStderr();
     for(int i=0;i<numModels;i++){
         p["DEFAULT_MODEL"] = models[i];
         MaxEntParameters c(p);
@@ -196,6 +264,8 @@ TEST(Parameters,HighFrequencyCheck){
 
         EXPECT_TRUE(std::abs(1+limit)<.1);
     }
+    const std::string diagnostics = testing::internal::GetCapturedStderr();
+    EXPECT_NE(diagnostics.find("The high frequency limit is not 1!"), std::string::npos);
 }
 
 TEST(Parameters,HDF5ContiParams){
@@ -326,9 +396,14 @@ TEST(Parameters,CovarianceDataInFile){
   p["DATA"]=pf;
   p["NDAT"] = 5;
   p["COVARIANCE_MATRIX"]=cov;
+  p["DATASPACE"]="frequency";
+  p["PARTICLE_HOLE_SYMMETRY"]=true;
 
   //MaxEntParameters handles covariance scaling
+  testing::internal::CaptureStderr();
   MaxEntParameters c(p);
+  const std::string diagnostics = testing::internal::GetCapturedStderr();
+  EXPECT_NE(diagnostics.find("The high frequency limit is not 1!"), std::string::npos);
   EXPECT_EQ(c.ndat(),5);
   EXPECT_EQ(c.T(),0.5);
 
@@ -402,4 +477,3 @@ TEST(Parameters,CovarianceHDF5Params){
   }
   std::remove(tf.c_str());
 }
-