about summary refs log tree commit diff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--.travis.yml58
-rw-r--r--Makefile23
-rw-r--r--src/bslmmdap.cpp19
-rw-r--r--src/debug.cpp11
-rw-r--r--src/debug.h15
-rw-r--r--src/io.cpp160
-rw-r--r--src/lm.cpp12
-rw-r--r--src/lmm.cpp20
-rw-r--r--src/mvlmm.cpp20
-rw-r--r--src/param.cpp8
-rw-r--r--src/prdt.cpp2
-rw-r--r--src/vc.cpp35
-rwxr-xr-xtest/dev_test_suite.sh43
13 files changed, 251 insertions, 175 deletions
diff --git a/.travis.yml b/.travis.yml
index bf1cf11..3607992 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -1,46 +1,52 @@
 language: C++
-compiler: gcc
 matrix:
+  # OSX testing is under development
+  # allow_failures:
+  #   - os: osx
   include:
     - os: linux
+      compiler: gcc
       addons:
         apt:
           sources:
             - ubuntu-toolchain-r-test
           packages:
+            # Our dev environment is a more recent GNU C++ and GSL2
             - g++-4.9
+            - libopenblas-dev 
+            - zlib1g-dev
+            - libeigen3-dev
+            - libgsl0-dev
+            - liblapack-dev
+            # - gfortran-dev for static
       env:
-         - MATRIX_EVAL="CC=gcc-4.9 && CXX=g++-4.9"
-    - os: linux
-      addons:
-        apt:
-          sources:
-            - ubuntu-toolchain-r-test
-          packages:
-            - g++-6
+         - MATRIX_EVAL="CC=gcc-4.9 && CXX=g++-4.9 && EIGEN_INCLUDE_PATH=/usr/include/eigen3"
+    - os: osx
+      compiler: clang
       env:
-        - MATRIX_EVAL="CC=gcc-6 && CXX=g++-6"
+         - MATRIX_EVAL="EIGEN_INCLUDE_PATH=/usr/local/include/eigen3"
+#     - os: linux
+#       addons:
+#         apt:
+#           sources:
+#             - ubuntu-toolchain-r-test
+#           packages:
+#             - g++-6
+#       env:
+#         - MATRIX_EVAL="CC=gcc-6 && CXX=g++-6"
 before_install:
-  - sudo apt-get -qq update
-  - sudo apt-get install -y libopenblas-dev zlib1g-dev
-  - sudo apt-get install -y libeigen3-dev
-  - sudo apt-get install -y libgsl0-dev
-  - sudo apt-get install -y liblapack-dev
-  # for the static release version we need the following
-  # - sudo apt-get install -y gfortran-dev
-  - dpkg -l
-  - eval "${MATRIX_EVAL}"
-  - $CXX --version
+  - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then brew update && brew install gsl openblas zlib eigen lapack ; fi
 script:
+  - echo $MATRIX_EVAL
   - eval "${MATRIX_EVAL}"
   - $CXX --version
   # build and test debug version
-  - make CXX=$CXX WITH_OPENBLAS=1 OPENBLAS_LEGACY=1 -j 4
-  - time make CXX=$CXX WITH_OPENBLAS=1 check
-  - make clean
-  # build and test release version
-  - make CXX=$CXX DEBUG= FORCE_DYNAMIC=1 WITH_OPENBLAS=1 OPENBLAS_LEGACY=1 -j 4
-  - time make CXX=$CXX DEBUG= WITH_OPENBLAS=1 check
+  - make CXX=$CXX EIGEN_INCLUDE_PATH=$EIGEN_INCLUDE_PATH WITH_OPENBLAS=1 OPENBLAS_LEGACY=1 -j 4
+  - time make CXX=$CXX WITH_OPENBLAS=1 EIGEN_INCLUDE_PATH=$EIGEN_INCLUDE_PATH WITH_OPENBLAS=1 OPENBLAS_LEGACY=1  check
+  # - make clean
+  # build and test release version (integration test mostly)
+  # - make CXX=$CXX EIGEN_INCLUDE_PATH=$EIGEN_INCLUDE_PATH DEBUG= FORCE_DYNAMIC=1 WITH_OPENBLAS=1 OPENBLAS_LEGACY=1 -j 4
+  # - time make CXX=$CXX DEBUG= WITH_OPENBLAS=1 fast-check
   # build static release (fast-check only)
   # - make clean
   # - make CXX=$CXX TRAVIS_CI=1 -j 4 fast-check
diff --git a/Makefile b/Makefile
index dc9a0d1..c57567b 100644
--- a/Makefile
+++ b/Makefile
@@ -42,15 +42,15 @@
 SYS                    = LNX # LNX|MAC (Linux is the default)
 # Leave blank after "=" to disable; put "= 1" to enable
 DIST_NAME              = gemma-0.97.3
-DEBUG                  = 1   # DEBUG mode, set DEBUG=0 for a release
+DEBUG                  = 1                # DEBUG mode, set DEBUG=0 for a release
 SHOW_COMPILER_WARNINGS =
 WITH_LAPACK            = 1
-WITH_OPENBLAS          =     # Defaults to LAPACK - OPENBLAS may be faster
-OPENBLAS_LEGACY        =     # Using older OpenBlas
-FORCE_STATIC           =     # Static linking of libraries
-GCC_FLAGS              = -O3 # extra flags -Wl,--allow-multiple-definition
-TRAVIS_CI              =     # used by TRAVIS for testing
-EIGEN_INCLUDE_PATH=/usr/include/eigen3
+WITH_OPENBLAS          =                  # Defaults to LAPACK - OPENBLAS may be faster
+OPENBLAS_LEGACY        =                  # Using older OpenBlas
+FORCE_STATIC           =                  # Static linking of libraries
+GCC_FLAGS              = -O3 -std=gnu++11 # extra flags -Wl,--allow-multiple-definition
+TRAVIS_CI              =                  # used by TRAVIS for testing
+EIGEN_INCLUDE_PATH     = /usr/include/eigen3
 
 # --------------------------------------------------------------------
 # Edit below this line with caution
@@ -67,6 +67,11 @@ else
   CPP = g++
 endif
 
+ifeq ($(CPP), clang++)
+  # macOS Homebrew settings (as used on Travis-CI)
+  GCC_FLAGS=-O3 -std=c++11 -stdlib=libc++ -isystem//usr/local/opt/openblas/include -isystem//usr/local/include/eigen3 -Wl,-L/usr/local/opt/openblas/lib
+endif
+
 ifdef WITH_OPENBLAS
   OPENBLAS=1
   # WITH_LAPACK =  # OPENBLAS usually includes LAPACK
@@ -77,10 +82,10 @@ ifdef WITH_OPENBLAS
 endif
 
 ifdef DEBUG
-  CPPFLAGS += -g $(GCC_FLAGS) -std=gnu++11 -isystem/$(EIGEN_INCLUDE_PATH) -Icontrib/catch-1.9.7 -Isrc
+  CPPFLAGS += -g $(GCC_FLAGS) -isystem/$(EIGEN_INCLUDE_PATH) -Icontrib/catch-1.9.7 -Isrc
 else
   # release mode
-  CPPFLAGS += -DNDEBUG $(GCC_FLAGS) -std=gnu++11 -isystem/$(EIGEN_INCLUDE_PATH) -Icontrib/catch-1.9.7 -Isrc
+  CPPFLAGS += -DNDEBUG $(GCC_FLAGS) -isystem/$(EIGEN_INCLUDE_PATH) -Icontrib/catch-1.9.7 -Isrc
 endif
 
 ifdef SHOW_COMPILER_WARNINGS
diff --git a/src/bslmmdap.cpp b/src/bslmmdap.cpp
index 7aac1d4..6f9aba7 100644
--- a/src/bslmmdap.cpp
+++ b/src/bslmmdap.cpp
@@ -116,16 +116,16 @@ void ReadFile_hyb(const string &file_hyp, vector<double> &vec_sa2,
   getline(infile, line);
 
   while (!safeGetline(infile, line).eof()) {
-    ch_ptr = strtok((char *)line.c_str(), " , \t");
-    ch_ptr = strtok(NULL, " , \t");
+    ch_ptr = strtok_safe((char *)line.c_str(), " , \t");
+    ch_ptr = strtok_safe(NULL, " , \t");
 
-    ch_ptr = strtok(NULL, " , \t");
+    ch_ptr = strtok_safe(NULL, " , \t");
     vec_sa2.push_back(atof(ch_ptr));
 
-    ch_ptr = strtok(NULL, " , \t");
+    ch_ptr = strtok_safe(NULL, " , \t");
     vec_sb2.push_back(atof(ch_ptr));
 
-    ch_ptr = strtok(NULL, " , \t");
+    ch_ptr = strtok_safe(NULL, " , \t");
     vec_wab.push_back(atof(ch_ptr));
   }
 
@@ -160,11 +160,11 @@ void ReadFile_bf(const string &file_bf, vector<string> &vec_rs,
   while (!safeGetline(infile, line).eof()) {
     flag_block = 0;
 
-    ch_ptr = strtok((char *)line.c_str(), " , \t");
+    ch_ptr = strtok_safe((char *)line.c_str(), " , \t");
     rs = ch_ptr;
     vec_rs.push_back(rs);
 
-    ch_ptr = strtok(NULL, " , \t");
+    ch_ptr = strtok_safe(NULL, " , \t");
     if (t == 0) {
       block = ch_ptr;
     } else {
@@ -223,7 +223,7 @@ void ReadFile_cat(const string &file_cat, const vector<string> &vec_rs,
 
   // Read header.
   HEADER header;
-  !safeGetline(infile, line).eof();
+  safeGetline(infile, line).eof();
   ReadHeader_io(line, header);
 
   // Use the header to determine the number of categories.
@@ -238,7 +238,7 @@ void ReadFile_cat(const string &file_cat, const vector<string> &vec_rs,
 
   // Read the following lines to record mapRS2cat.
   while (!safeGetline(infile, line).eof()) {
-    ch_ptr = strtok((char *)line.c_str(), " , \t");
+    ch_ptr = strtok_safe((char *)line.c_str(), " , \t");
 
     if (header.rs_col == 0) {
       rs = chr + ":" + pos;
@@ -248,6 +248,7 @@ void ReadFile_cat(const string &file_cat, const vector<string> &vec_rs,
     catd.clear();
 
     for (size_t i = 0; i < header.coln; i++) {
+      enforce(ch_ptr);
       if (header.rs_col != 0 && header.rs_col == i + 1) {
         rs = ch_ptr;
       } else if (header.chr_col != 0 && header.chr_col == i + 1) {
diff --git a/src/debug.cpp b/src/debug.cpp
index 6bd834e..82d2be0 100644
--- a/src/debug.cpp
+++ b/src/debug.cpp
@@ -72,6 +72,13 @@ gsl_vector *gsl_vector_safe_alloc(size_t n) {
   return v;
 }
 
+char *do_strtok_safe(char *tokenize, const char *delimiters, const char *__pretty_function, const char *__file, int __line) {
+  auto token = strtok(tokenize,delimiters);
+  if (token == NULL && (is_debug_mode() || is_strict_mode()))
+    fail_at_msg(__file,__line,string("strtok failed in ") + __pretty_function);
+  return token;
+}
+
 // Helper function called by macro validate_K(K, check). K is validated
 // unless -no-check option is used.
 void do_validate_K(const gsl_matrix *K, const char *__file, int __line) {
@@ -88,13 +95,13 @@ void do_validate_K(const gsl_matrix *K, const char *__file, int __line) {
     if (!isMatrixIllConditioned(eigenvalues))
       warning_at_msg(__file,__line,"K is ill conditioned!");
     if (!isMatrixSymmetric(K))
-      fail_at_msg(is_strict_mode(),__file,__line,"K is not symmetric!" );
+      warnfail_at_msg(is_strict_mode(),__file,__line,"K is not symmetric!" );
     const bool negative_values = has_negative_values_but_one(eigenvalues);
     if (negative_values) {
       warning_at_msg(__file,__line,"K has more than one negative eigenvalues!");
     }
     if (count_small>1 && negative_values && !isMatrixPositiveDefinite(K))
-      fail_at_msg(is_strict_mode(),__file,__line,"K is not positive definite!");
+      warnfail_at_msg(is_strict_mode(),__file,__line,"K is not positive definite!");
     gsl_vector_free(eigenvalues);
   }
 }
diff --git a/src/debug.h b/src/debug.h
index c127630..b3ec17b 100644
--- a/src/debug.h
+++ b/src/debug.h
@@ -28,6 +28,9 @@ bool is_legacy_mode();
 gsl_matrix *gsl_matrix_safe_alloc(size_t rows,size_t cols);
 gsl_vector *gsl_vector_safe_alloc(size_t n);
 
+char *do_strtok_safe(char *tokenize, const char *delimiters, const char *__pretty_function, const char *__file, int __line);
+#define strtok_safe(string,delimiters) do_strtok_safe(string,delimiters,__PRETTY_FUNCTION__,__FILE__,__LINE__)
+
 // Validation routines
 void do_validate_K(const gsl_matrix *K, const char *__file, int __line);
 
@@ -36,7 +39,7 @@ void do_validate_K(const gsl_matrix *K, const char *__file, int __line);
 
 #define warning_at_msg(__file,__line,msg) cerr << "**** WARNING: " << msg << " in " << __file << " at line " << __line << endl;
 
-inline void fail_at_msg(bool strict, const char *__file, int __line, const char *msg) {
+inline void warnfail_at_msg(bool strict, const char *__file, int __line, const char *msg) {
   if (strict)
     std::cerr << "**** STRICT FAIL: ";
   else
@@ -46,6 +49,11 @@ inline void fail_at_msg(bool strict, const char *__file, int __line, const char
     exit(1);
 }
 
+inline void fail_at_msg(const char *__file, int __line, std::string msg) {
+  std::cerr << msg << " in " << __file << " at line " << __line << std::endl;
+  exit(1);
+}
+
 # ifndef __ASSERT_VOID_CAST
 # define __ASSERT_VOID_CAST (void)
 # endif
@@ -55,6 +63,11 @@ inline void fail_msg(const char *msg) {
   exit(5);
 }
 
+inline void fail_msg(std::string msg) {
+  std::cerr << "**** FAILED: " << msg << std::endl;
+  exit(5);
+}
+
 #if defined NDEBUG
 
 #define warning_msg(msg) cerr << "**** WARNING: " << msg << endl;
diff --git a/src/io.cpp b/src/io.cpp
index 8abdeec..35a59ee 100644
--- a/src/io.cpp
+++ b/src/io.cpp
@@ -152,7 +152,7 @@ bool ReadFile_snps_header(const string &file_snps, set<string> &setSnps) {
 
   // Read header.
   HEADER header;
-  !safeGetline(infile, line).eof();
+  safeGetline(infile, line).eof();
   ReadHeader_io(line, header);
 
   if (header.rs_col == 0 && (header.chr_col == 0 || header.pos_col == 0)) {
@@ -208,7 +208,7 @@ bool ReadFile_log(const string &file_log, double &pheno_mean) {
   size_t flag = 0;
 
   while (getline(infile, line)) {
-    ch_ptr = strtok((char *)line.c_str(), " , \t");
+    ch_ptr = strtok_safe((char *)line.c_str(), " , \t");
     ch_ptr = strtok(NULL, " , \t");
 
     if (ch_ptr != NULL && strcmp(ch_ptr, "estimated") == 0) {
@@ -216,7 +216,7 @@ bool ReadFile_log(const string &file_log, double &pheno_mean) {
       if (ch_ptr != NULL && strcmp(ch_ptr, "mean") == 0) {
         ch_ptr = strtok(NULL, " , \t");
         if (ch_ptr != NULL && strcmp(ch_ptr, "=") == 0) {
-          ch_ptr = strtok(NULL, " , \t");
+          ch_ptr = strtok_safe(NULL, " , \t");
           pheno_mean = atof(ch_ptr);
           flag = 1;
         }
@@ -314,7 +314,7 @@ bool ReadFile_column(const string &file_pheno, vector<int> &indicator_idv,
   string id;
   double p;
   while (!safeGetline(infile, line).eof()) {
-    ch_ptr = strtok((char *)line.c_str(), " , \t");
+    ch_ptr = strtok_safe((char *)line.c_str(), " , \t");
     for (int i = 0; i < (p_column - 1); ++i) {
       ch_ptr = strtok(NULL, " , \t");
     }
@@ -486,17 +486,17 @@ bool ReadFile_bim(const string &file_bim, vector<SNPINFO> &snpInfo) {
   string minor;
 
   while (getline(infile, line)) {
-    ch_ptr = strtok((char *)line.c_str(), " \t");
+    ch_ptr = strtok_safe((char *)line.c_str(), " \t");
     chr = ch_ptr;
-    ch_ptr = strtok(NULL, " \t");
+    ch_ptr = strtok_safe(NULL, " \t");
     rs = ch_ptr;
-    ch_ptr = strtok(NULL, " \t");
+    ch_ptr = strtok_safe(NULL, " \t");
     cM = atof(ch_ptr);
-    ch_ptr = strtok(NULL, " \t");
+    ch_ptr = strtok_safe(NULL, " \t");
     b_pos = atol(ch_ptr);
-    ch_ptr = strtok(NULL, " \t");
+    ch_ptr = strtok_safe(NULL, " \t");
     minor = ch_ptr;
-    ch_ptr = strtok(NULL, " \t");
+    ch_ptr = strtok_safe(NULL, " \t");
     major = ch_ptr;
 
     SNPINFO sInfo = {chr, rs, cM, b_pos, minor, major, 0, -9, -9, 0, 0, 0};
@@ -542,12 +542,12 @@ bool ReadFile_fam(const string &file_fam, vector<vector<int>> &indicator_pheno,
   }
 
   while (!safeGetline(infile, line).eof()) {
-    ch_ptr = strtok((char *)line.c_str(), " \t");
-    ch_ptr = strtok(NULL, " \t");
+    ch_ptr = strtok_safe((char *)line.c_str(), " \t");
+    ch_ptr = strtok_safe(NULL, " \t");
     id = ch_ptr;
-    ch_ptr = strtok(NULL, " \t");
-    ch_ptr = strtok(NULL, " \t");
-    ch_ptr = strtok(NULL, " \t");
+    ch_ptr = strtok_safe(NULL, " \t");
+    ch_ptr = strtok_safe(NULL, " \t");
+    ch_ptr = strtok_safe(NULL, " \t");
     ch_ptr = strtok(NULL, " \t");
 
     size_t i = 0;
@@ -649,11 +649,11 @@ bool ReadFile_geno(const string &file_geno, const set<string> &setSnps,
   file_pos = 0;
   auto count_warnings = 0;
   while (!safeGetline(infile, line).eof()) {
-    ch_ptr = strtok((char *)line.c_str(), " , \t");
+    ch_ptr = strtok_safe((char *)line.c_str(), " , \t");
     rs = ch_ptr;
-    ch_ptr = strtok(NULL, " , \t");
+    ch_ptr = strtok_safe(NULL, " , \t");
     minor = ch_ptr;
-    ch_ptr = strtok(NULL, " , \t");
+    ch_ptr = strtok_safe(NULL, " , \t");
     major = ch_ptr;
 
     if (setSnps.size() != 0 && setSnps.count(rs) == 0) {
@@ -694,7 +694,7 @@ bool ReadFile_geno(const string &file_geno, const set<string> &setSnps,
     c_idv = 0;
     gsl_vector_set_zero(genotype_miss);
     for (int i = 0; i < ni_total; ++i) {
-      ch_ptr = strtok(NULL, " , \t");
+      ch_ptr = strtok_safe(NULL, " , \t");
       if (indicator_idv[i] == 0)
         continue;
 
@@ -1004,13 +1004,13 @@ bool Bimbam_ReadOneSNP(const size_t inc, const vector<int> &indicator_idv,
   bool flag = false;
 
   for (size_t i = 0; i < inc; i++) {
-    !safeGetline(infile, line).eof();
+    safeGetline(infile, line).eof();
   }
 
   if (!safeGetline(infile, line).eof()) {
-    ch_ptr = strtok((char *)line.c_str(), " , \t");
-    ch_ptr = strtok(NULL, " , \t");
-    ch_ptr = strtok(NULL, " , \t");
+    ch_ptr = strtok_safe((char *)line.c_str(), " , \t");
+    ch_ptr = strtok_safe(NULL, " , \t");
+    ch_ptr = strtok_safe(NULL, " , \t");
 
     geno_mean = 0.0;
     double d;
@@ -1018,7 +1018,7 @@ bool Bimbam_ReadOneSNP(const size_t inc, const vector<int> &indicator_idv,
     vector<size_t> geno_miss;
 
     for (size_t i = 0; i < ni_total; ++i) {
-      ch_ptr = strtok(NULL, " , \t");
+      ch_ptr = strtok_safe(NULL, " , \t");
       if (indicator_idv[i] == 0) {
         continue;
       }
@@ -1134,9 +1134,7 @@ void ReadFile_kin(const string &file_kin, vector<int> &indicator_idv,
     size_t i_test = 0, i_total = 0, j_test = 0, j_total = 0;
     while (getline(infile, line)) {
       if (i_total == ni_total) {
-        cout << "error! number of rows in the kinship "
-             << "file is larger than the number of phentypes." << endl;
-        error = true;
+        fail_msg("number of rows in the kinship file is larger than the number of phentypes");
       }
 
       if (indicator_idv[i_total] == 0) {
@@ -1149,10 +1147,7 @@ void ReadFile_kin(const string &file_kin, vector<int> &indicator_idv,
       ch_ptr = strtok((char *)line.c_str(), " , \t");
       while (ch_ptr != NULL) {
         if (j_total == ni_total) {
-          cout << "error! number of columns in the "
-               << "kinship file is larger than the number"
-               << " of phenotypes for row = " << i_total << endl;
-          error = true;
+          fail_msg(string("number of columns in the kinship file is larger than the number of individuals for row = ")+to_string(i_total));
         }
 
         d = atof(ch_ptr);
@@ -1165,18 +1160,14 @@ void ReadFile_kin(const string &file_kin, vector<int> &indicator_idv,
         ch_ptr = strtok(NULL, " , \t");
       }
       if (j_total != ni_total) {
-        cout << "error! number of columns in the kinship "
-             << "file do not match the number of phentypes for "
-             << "row = " << i_total << endl;
-        error = true;
+        string msg = "number of columns in the kinship file does not match the number of individuals for row = " + to_string( i_total );
+        fail_msg(msg);
       }
       i_total++;
       i_test++;
     }
     if (i_total != ni_total) {
-      cout << "error! number of rows in the kinship file do "
-           << "not match the number of phenotypes." << endl;
-      error = true;
+      fail_msg("number of rows in the kinship file does not match the number of individuals.");
     }
   } else {
     map<size_t, size_t> mapID2ID;
@@ -1193,11 +1184,11 @@ void ReadFile_kin(const string &file_kin, vector<int> &indicator_idv,
     size_t n_id1, n_id2;
 
     while (getline(infile, line)) {
-      ch_ptr = strtok((char *)line.c_str(), " , \t");
+      ch_ptr = strtok_safe((char *)line.c_str(), " , \t");
       id1 = ch_ptr;
-      ch_ptr = strtok(NULL, " , \t");
+      ch_ptr = strtok_safe(NULL, " , \t");
       id2 = ch_ptr;
-      ch_ptr = strtok(NULL, " , \t");
+      ch_ptr = strtok_safe(NULL, " , \t");
       d = atof(ch_ptr);
       if (mapID2num.count(id1) == 0 || mapID2num.count(id2) == 0) {
         continue;
@@ -1212,9 +1203,10 @@ void ReadFile_kin(const string &file_kin, vector<int> &indicator_idv,
 
       Cov_d = gsl_matrix_get(G, n_id1, n_id2);
       if (Cov_d != 0 && Cov_d != d) {
-        cout << "error! redundant and unequal terms in the "
+        cerr << "error! redundant and unequal terms in the "
              << "kinship file, for id1 = " << id1 << " and id2 = " << id2
              << endl;
+        fail_msg("");
       } else {
         gsl_matrix_set(G, n_id1, n_id2, d);
         gsl_matrix_set(G, n_id2, n_id1, d);
@@ -1253,7 +1245,6 @@ void ReadFile_mk(const string &file_mk, vector<int> &indicator_idv,
 
   infile.close();
   infile.clear();
-  return;
 }
 
 void ReadFile_eigenU(const string &file_ku, bool &error, gsl_matrix *U) {
@@ -1329,7 +1320,7 @@ void ReadFile_eigenD(const string &file_kd, bool &error, gsl_vector *eval) {
       error = true;
     }
 
-    ch_ptr = strtok((char *)line.c_str(), " , \t");
+    ch_ptr = strtok_safe((char *)line.c_str(), " , \t");
     d = atof(ch_ptr);
 
     ch_ptr = strtok(NULL, " , \t");
@@ -1380,7 +1371,7 @@ bool BimbamKin(const string file_geno, const set<string> ksnps,
   size_t ns_test = 0;
   for (size_t t = 0; t < indicator_snp.size(); ++t) {
     string line;
-    !safeGetline(infile, line).eof();
+    safeGetline(infile, line).eof();
     if (t % display_pace == 0 || t == (indicator_snp.size() - 1)) {
       ProgressBar("Reading SNPs", t, indicator_snp.size() - 1);
     }
@@ -1662,21 +1653,21 @@ bool ReadFile_geno(const string file_geno, vector<int> &indicator_idv,
   int c_idv = 0, c_snp = 0;
 
   for (int i = 0; i < ns_total; ++i) {
-    !safeGetline(infile, line).eof();
+    safeGetline(infile, line).eof();
     if (indicator_snp[i] == 0) {
       continue;
     }
 
-    ch_ptr = strtok((char *)line.c_str(), " , \t");
-    ch_ptr = strtok(NULL, " , \t");
-    ch_ptr = strtok(NULL, " , \t");
+    ch_ptr = strtok_safe((char *)line.c_str(), " , \t");
+    ch_ptr = strtok_safe(NULL, " , \t");
+    ch_ptr = strtok_safe(NULL, " , \t");
 
     c_idv = 0;
     geno_mean = 0;
     n_miss = 0;
     gsl_vector_set_zero(genotype_miss);
     for (int j = 0; j < ni_total; ++j) {
-      ch_ptr = strtok(NULL, " , \t");
+      ch_ptr = strtok_safe(NULL, " , \t");
       if (indicator_idv[j] == 0) {
         continue;
       }
@@ -1771,21 +1762,21 @@ bool ReadFile_geno(const string &file_geno, vector<int> &indicator_idv,
   size_t c_idv = 0, c_snp = 0;
 
   for (size_t i = 0; i < ns_total; ++i) {
-    !safeGetline(infile, line).eof();
+    safeGetline(infile, line).eof();
     if (indicator_snp[i] == 0) {
       continue;
     }
 
-    ch_ptr = strtok((char *)line.c_str(), " , \t");
-    ch_ptr = strtok(NULL, " , \t");
-    ch_ptr = strtok(NULL, " , \t");
+    ch_ptr = strtok_safe((char *)line.c_str(), " , \t");
+    ch_ptr = strtok_safe(NULL, " , \t");
+    ch_ptr = strtok_safe(NULL, " , \t");
 
     c_idv = 0;
     geno_mean = 0;
     n_miss = 0;
     gsl_vector_set_zero(genotype_miss);
     for (uint j = 0; j < ni_total; ++j) {
-      ch_ptr = strtok(NULL, " , \t");
+      ch_ptr = strtok_safe(NULL, " , \t");
       if (indicator_idv[j] == 0) {
         continue;
       }
@@ -2135,7 +2126,7 @@ bool ReadFile_est(const string &file_est, const vector<size_t> &est_column,
   size_t n = *max_element(est_column.begin(), est_column.end());
 
   while (getline(infile, line)) {
-    ch_ptr = strtok((char *)line.c_str(), " \t");
+    ch_ptr = strtok_safe((char *)line.c_str(), " \t");
 
     alpha = 0.0;
     beta = 0.0;
@@ -2154,7 +2145,7 @@ bool ReadFile_est(const string &file_est, const vector<size_t> &est_column,
         gamma = atof(ch_ptr);
       }
       if (i < n) {
-        ch_ptr = strtok(NULL, " \t");
+        ch_ptr = strtok_safe(NULL, " \t");
       }
     }
 
@@ -2212,7 +2203,7 @@ bool ReadFile_gene(const string &file_gene, vector<double> &vec_read,
   getline(infile, line);
 
   while (getline(infile, line)) {
-    ch_ptr = strtok((char *)line.c_str(), " , \t");
+    ch_ptr = strtok_safe((char *)line.c_str(), " , \t");
     rs = ch_ptr;
 
     ch_ptr = strtok(NULL, " , \t");
@@ -2536,7 +2527,7 @@ bool ReadFile_cat(const string &file_cat, map<string, size_t> &mapRS2cat,
 
   // Read header.
   HEADER header;
-  !safeGetline(infile, line).eof();
+  safeGetline(infile, line).eof();
   ReadHeader_io(line, header);
 
   // Use the header to count the number of categories.
@@ -2562,10 +2553,11 @@ bool ReadFile_cat(const string &file_cat, map<string, size_t> &mapRS2cat,
 
   // Read the following lines to record mapRS2cat.
   while (!safeGetline(infile, line).eof()) {
-    ch_ptr = strtok((char *)line.c_str(), " , \t");
+    ch_ptr = strtok_safe((char *)line.c_str(), " , \t");
 
     i_cat = 0;
     for (size_t i = 0; i < header.coln; i++) {
+      enforce(ch_ptr);
       if (header.rs_col != 0 && header.rs_col == i + 1) {
         rs = ch_ptr;
       } else if (header.chr_col != 0 && header.chr_col == i + 1) {
@@ -2686,16 +2678,16 @@ bool BimbamKinUncentered(const string &file_geno, const set<string> ksnps,
 
   size_t ns_test = 0;
   for (size_t t = 0; t < indicator_snp.size(); ++t) {
-    !safeGetline(infile, line).eof();
+    safeGetline(infile, line).eof();
     if (t % display_pace == 0 || t == (indicator_snp.size() - 1)) {
       ProgressBar("Reading SNPs", t, indicator_snp.size() - 1);
     }
     if (indicator_snp[t] == 0)
       continue;
 
-    ch_ptr = strtok((char *)line.c_str(), " , \t");
-    ch_ptr = strtok(NULL, " , \t");
-    ch_ptr = strtok(NULL, " , \t");
+    ch_ptr = strtok_safe((char *)line.c_str(), " , \t");
+    ch_ptr = strtok_safe(NULL, " , \t");
+    ch_ptr = strtok_safe(NULL, " , \t");
 
     rs = snpInfo[t].rs_number; // This line is new.
 
@@ -2709,7 +2701,7 @@ bool BimbamKinUncentered(const string &file_geno, const set<string> ksnps,
       if (indicator_idv[i] == 0) {
         continue;
       }
-      ch_ptr = strtok(NULL, " , \t");
+      ch_ptr = strtok_safe(NULL, " , \t");
       if (strcmp(ch_ptr, "NA") == 0) {
         gsl_vector_set(geno_miss, i, 0);
         n_miss++;
@@ -3151,9 +3143,9 @@ bool ReadFile_wsnp(const string &file_wsnp, map<string, double> &mapRS2weight) {
   double weight;
 
   while (!safeGetline(infile, line).eof()) {
-    ch_ptr = strtok((char *)line.c_str(), " , \t");
+    ch_ptr = strtok_safe((char *)line.c_str(), " , \t");
     rs = ch_ptr;
-    ch_ptr = strtok(NULL, " , \t");
+    ch_ptr = strtok_safe(NULL, " , \t");
     weight = atof(ch_ptr);
     mapRS2weight[rs] = weight;
   }
@@ -3182,17 +3174,18 @@ bool ReadFile_wsnp(const string &file_wcat, const size_t n_vc,
 
   // Read header.
   HEADER header;
-  !safeGetline(infile, line).eof();
+  safeGetline(infile, line).eof();
   ReadHeader_io(line, header);
 
   while (!safeGetline(infile, line).eof()) {
     if (isBlankLine(line)) {
       continue;
     }
-    ch_ptr = strtok((char *)line.c_str(), " , \t");
+    ch_ptr = strtok_safe((char *)line.c_str(), " , \t");
 
     size_t t = 0;
     for (size_t i = 0; i < header.coln; i++) {
+      enforce(ch_ptr);
       if (header.rs_col != 0 && header.rs_col == i + 1) {
         rs = ch_ptr;
       } else if (header.chr_col != 0 && header.chr_col == i + 1) {
@@ -3274,7 +3267,7 @@ void ReadFile_beta(const string &file_beta,
 
   // Read header.
   HEADER header;
-  !safeGetline(infile, line).eof();
+  safeGetline(infile, line).eof();
   ReadHeader_io(line, header);
 
   if (header.n_col == 0) {
@@ -3296,7 +3289,7 @@ void ReadFile_beta(const string &file_beta,
     if (isBlankLine(line)) {
       continue;
     }
-    ch_ptr = strtok((char *)line.c_str(), " , \t");
+    ch_ptr = strtok_safe((char *)line.c_str(), " , \t");
 
     z = 0;
     beta = 0;
@@ -3311,6 +3304,7 @@ void ReadFile_beta(const string &file_beta,
     af = 0;
     var_x = 0;
     for (size_t i = 0; i < header.coln; i++) {
+      enforce(ch_ptr);
       if (header.rs_col != 0 && header.rs_col == i + 1) {
         rs = ch_ptr;
       }
@@ -3456,7 +3450,7 @@ void ReadFile_beta(const string &file_beta, const map<string, double> &mapRS2wA,
 
   // Read header.
   HEADER header;
-  !safeGetline(infile, line).eof();
+  safeGetline(infile, line).eof();
   ReadHeader_io(line, header);
 
   if (header.n_col == 0) {
@@ -3477,7 +3471,7 @@ void ReadFile_beta(const string &file_beta, const map<string, double> &mapRS2wA,
     if (isBlankLine(line)) {
       continue;
     }
-    ch_ptr = strtok((char *)line.c_str(), " , \t");
+    ch_ptr = strtok_safe((char *)line.c_str(), " , \t");
 
     z = 0;
     beta = 0;
@@ -3492,6 +3486,7 @@ void ReadFile_beta(const string &file_beta, const map<string, double> &mapRS2wA,
     af = 0;
     var_x = 0;
     for (size_t i = 0; i < header.coln; i++) {
+      enforce(ch_ptr);
       if (header.rs_col != 0 && header.rs_col == i + 1) {
         rs = ch_ptr;
       }
@@ -3762,8 +3757,8 @@ void ReadFile_vector(const string &file_vec, gsl_vector *vec) {
   char *ch_ptr;
 
   for (size_t i = 0; i < vec->size; i++) {
-    !safeGetline(infile, line).eof();
-    ch_ptr = strtok((char *)line.c_str(), " , \t");
+    safeGetline(infile, line).eof();
+    ch_ptr = strtok_safe((char *)line.c_str(), " , \t");
     gsl_vector_set(vec, i, atof(ch_ptr));
   }
 
@@ -3785,9 +3780,10 @@ void ReadFile_matrix(const string &file_mat, gsl_matrix *mat) {
   char *ch_ptr;
 
   for (size_t i = 0; i < mat->size1; i++) {
-    !safeGetline(infile, line).eof();
-    ch_ptr = strtok((char *)line.c_str(), " , \t");
+    safeGetline(infile, line).eof();
+    ch_ptr = strtok_safe((char *)line.c_str(), " , \t");
     for (size_t j = 0; j < mat->size2; j++) {
+      enforce(ch_ptr);
       gsl_matrix_set(mat, i, j, atof(ch_ptr));
       ch_ptr = strtok(NULL, " , \t");
     }
@@ -3812,18 +3808,20 @@ void ReadFile_matrix(const string &file_mat, gsl_matrix *mat1,
   char *ch_ptr;
 
   for (size_t i = 0; i < mat1->size1; i++) {
-    !safeGetline(infile, line).eof();
-    ch_ptr = strtok((char *)line.c_str(), " , \t");
+    safeGetline(infile, line).eof();
+    ch_ptr = strtok_safe((char *)line.c_str(), " , \t");
     for (size_t j = 0; j < mat1->size2; j++) {
+      enforce(ch_ptr);
       gsl_matrix_set(mat1, i, j, atof(ch_ptr));
       ch_ptr = strtok(NULL, " , \t");
     }
   }
 
   for (size_t i = 0; i < mat2->size1; i++) {
-    !safeGetline(infile, line).eof();
-    ch_ptr = strtok((char *)line.c_str(), " , \t");
+    safeGetline(infile, line).eof();
+    ch_ptr = strtok_safe((char *)line.c_str(), " , \t");
     for (size_t j = 0; j < mat2->size2; j++) {
+      enforce(ch_ptr);
       gsl_matrix_set(mat2, i, j, atof(ch_ptr));
       ch_ptr = strtok(NULL, " , \t");
     }
diff --git a/src/lm.cpp b/src/lm.cpp
index ea628cf..b94a426 100644
--- a/src/lm.cpp
+++ b/src/lm.cpp
@@ -333,12 +333,12 @@ void LM::AnalyzeGene(const gsl_matrix *W, const gsl_vector *x) {
     if (t % d_pace == 0 || t == ng_total - 1) {
       ProgressBar("Performing Analysis", t, ng_total - 1);
     }
-    ch_ptr = strtok((char *)line.c_str(), " , \t");
+    ch_ptr = strtok_safe((char *)line.c_str(), " , \t");
     rs = ch_ptr;
 
     c_phen = 0;
     for (size_t i = 0; i < indicator_idv.size(); ++i) {
-      ch_ptr = strtok(NULL, " , \t");
+      ch_ptr = strtok_safe(NULL, " , \t");
       if (indicator_idv[i] == 0) {
         continue;
       }
@@ -427,16 +427,16 @@ void LM::AnalyzeBimbam(const gsl_matrix *W, const gsl_vector *y) {
       continue;
     }
 
-    ch_ptr = strtok((char *)line.c_str(), " , \t");
-    ch_ptr = strtok(NULL, " , \t");
-    ch_ptr = strtok(NULL, " , \t");
+    ch_ptr = strtok_safe((char *)line.c_str(), " , \t");
+    ch_ptr = strtok_safe(NULL, " , \t");
+    ch_ptr = strtok_safe(NULL, " , \t");
 
     x_mean = 0.0;
     c_phen = 0;
     n_miss = 0;
     gsl_vector_set_zero(x_miss);
     for (size_t i = 0; i < ni_total; ++i) {
-      ch_ptr = strtok(NULL, " , \t");
+      ch_ptr = strtok_safe(NULL, " , \t");
       if (indicator_idv[i] == 0) {
         continue;
       }
diff --git a/src/lmm.cpp b/src/lmm.cpp
index e80cd76..ae8b747 100644
--- a/src/lmm.cpp
+++ b/src/lmm.cpp
@@ -1197,16 +1197,16 @@ void LMM::AnalyzeGene(const gsl_matrix *U, const gsl_vector *eval,
   getline(infile, line);
 
   for (size_t t = 0; t < ng_total; t++) {
-    !safeGetline(infile, line).eof();
+    safeGetline(infile, line).eof();
     if (t % d_pace == 0 || t == ng_total - 1) {
       ProgressBar("Performing Analysis", t, ng_total - 1);
     }
-    ch_ptr = strtok((char *)line.c_str(), " , \t");
+    ch_ptr = strtok_safe((char *)line.c_str(), " , \t");
     rs = ch_ptr;
 
     c_phen = 0;
     for (size_t i = 0; i < indicator_idv.size(); ++i) {
-      ch_ptr = strtok(NULL, " , \t");
+      ch_ptr = strtok_safe(NULL, " , \t");
       if (indicator_idv[i] == 0) {
         continue;
       }
@@ -1472,8 +1472,8 @@ void LMM::AnalyzeBimbam(const gsl_matrix *U, const gsl_vector *eval,
     enforce_msg(ch_ptr, "Parsing BIMBAM genofile"); // ch_ptr should not be NULL
 
     auto snp = string(ch_ptr);
-    ch_ptr = strtok(NULL, " , \t"); // skip column
-    ch_ptr = strtok(NULL, " , \t"); // skip column
+    ch_ptr = strtok_safe(NULL, " , \t"); // skip column
+    ch_ptr = strtok_safe(NULL, " , \t"); // skip column
 
     gs.assign (ni_total,nan("")); // wipe values
 
@@ -1937,7 +1937,7 @@ void LMM::AnalyzeBimbamGXE(const gsl_matrix *U, const gsl_vector *eval,
 
   // Start reading genotypes and analyze.
   for (size_t t = 0; t < indicator_snp.size(); ++t) {
-    !safeGetline(infile, line).eof();
+    safeGetline(infile, line).eof();
     if (t % d_pace == 0 || t == (ns_total - 1)) {
       ProgressBar("Reading SNPs", t, ns_total - 1);
     }
@@ -1945,16 +1945,16 @@ void LMM::AnalyzeBimbamGXE(const gsl_matrix *U, const gsl_vector *eval,
       continue;
     }
 
-    ch_ptr = strtok((char *)line.c_str(), " , \t");
-    ch_ptr = strtok(NULL, " , \t");
-    ch_ptr = strtok(NULL, " , \t");
+    ch_ptr = strtok_safe((char *)line.c_str(), " , \t");
+    ch_ptr = strtok_safe(NULL, " , \t");
+    ch_ptr = strtok_safe(NULL, " , \t");
 
     x_mean = 0.0;
     c_phen = 0;
     n_miss = 0;
     gsl_vector_set_zero(x_miss);
     for (size_t i = 0; i < ni_total; ++i) {
-      ch_ptr = strtok(NULL, " , \t");
+      ch_ptr = strtok_safe(NULL, " , \t");
       if (indicator_idv[i] == 0) {
         continue;
       }
diff --git a/src/mvlmm.cpp b/src/mvlmm.cpp
index 95ef14b..bdcbe5b 100644
--- a/src/mvlmm.cpp
+++ b/src/mvlmm.cpp
@@ -3189,7 +3189,7 @@ void MVLMM::AnalyzeBimbam(const gsl_matrix *U, const gsl_vector *eval,
     t_last++;
   }
   for (size_t t = 0; t < indicator_snp.size(); ++t) {
-    !safeGetline(infile, line).eof();
+    safeGetline(infile, line).eof();
     if (t % d_pace == 0 || t == (ns_total - 1)) {
       ProgressBar("Reading SNPs", t, ns_total - 1);
     }
@@ -3197,16 +3197,16 @@ void MVLMM::AnalyzeBimbam(const gsl_matrix *U, const gsl_vector *eval,
       continue;
     }
 
-    ch_ptr = strtok((char *)line.c_str(), " , \t");
-    ch_ptr = strtok(NULL, " , \t");
-    ch_ptr = strtok(NULL, " , \t");
+    ch_ptr = strtok_safe((char *)line.c_str(), " , \t");
+    ch_ptr = strtok_safe(NULL, " , \t");
+    ch_ptr = strtok_safe(NULL, " , \t");
 
     x_mean = 0.0;
     c_phen = 0;
     n_miss = 0;
     gsl_vector_set_zero(x_miss);
     for (size_t i = 0; i < ni_total; ++i) {
-      ch_ptr = strtok(NULL, " , \t");
+      ch_ptr = strtok_safe(NULL, " , \t");
       if (indicator_idv[i] == 0) {
         continue;
       }
@@ -4166,7 +4166,7 @@ void MVLMM::AnalyzeBimbamGXE(const gsl_matrix *U, const gsl_vector *eval,
 
   // Start reading genotypes and analyze.
   for (size_t t = 0; t < indicator_snp.size(); ++t) {
-    !safeGetline(infile, line).eof();
+    safeGetline(infile, line).eof();
     if (t % d_pace == 0 || t == (ns_total - 1)) {
       ProgressBar("Reading SNPs", t, ns_total - 1);
     }
@@ -4174,16 +4174,16 @@ void MVLMM::AnalyzeBimbamGXE(const gsl_matrix *U, const gsl_vector *eval,
       continue;
     }
 
-    ch_ptr = strtok((char *)line.c_str(), " , \t");
-    ch_ptr = strtok(NULL, " , \t");
-    ch_ptr = strtok(NULL, " , \t");
+    ch_ptr = strtok_safe((char *)line.c_str(), " , \t");
+    ch_ptr = strtok_safe(NULL, " , \t");
+    ch_ptr = strtok_safe(NULL, " , \t");
 
     x_mean = 0.0;
     c_phen = 0;
     n_miss = 0;
     gsl_vector_set_zero(x_miss);
     for (size_t i = 0; i < ni_total; ++i) {
-      ch_ptr = strtok(NULL, " , \t");
+      ch_ptr = strtok_safe(NULL, " , \t");
       if (indicator_idv[i] == 0) {
         continue;
       }
diff --git a/src/param.cpp b/src/param.cpp
index 8aa7a64..1a27a53 100644
--- a/src/param.cpp
+++ b/src/param.cpp
@@ -16,12 +16,12 @@
     along with this program. If not, see <http://www.gnu.org/licenses/>.
 */
 
+#include <iostream>
+#include <string>
 #include <algorithm>
 #include <cmath>
 #include <cstring>
 #include <fstream>
-#include <iostream>
-#include <string>
 #include <sys/stat.h>
 
 #include "gsl/gsl_blas.h"
@@ -899,7 +899,7 @@ void PARAM::CheckParam(void) {
   if (!loco.empty()) {
     enforce_msg((a_mode >= 1 && a_mode <= 4) || a_mode == 21 || a_mode == 22,
                 "LOCO only works with LMM and K");
-    enforce_msg(file_bfile.empty(), "LOCO does not work with PLink (yet)");
+    // enforce_msg(file_bfile.empty(), "LOCO does not work with PLink (yet)");
     enforce_msg(file_gxe.empty(), "LOCO does not support GXE (yet)");
     enforce_msg(!file_anno.empty(),
                 "LOCO requires annotation file (-a switch)");
@@ -1317,7 +1317,7 @@ void PARAM::CalcKin(gsl_matrix *matrix_kin) {
 
   if (!file_bfile.empty()) {
     file_str = file_bfile + ".bed";
-    enforce_msg(loco.empty(), "FIXME: LOCO nyi");
+    // enforce_msg(loco.empty(), "FIXME: LOCO nyi");
     if (PlinkKin(file_str, indicator_snp, a_mode - 20, d_pace, matrix_kin) ==
         false) {
       error = true;
diff --git a/src/prdt.cpp b/src/prdt.cpp
index 9dc84bc..fc0abe8 100644
--- a/src/prdt.cpp
+++ b/src/prdt.cpp
@@ -227,7 +227,7 @@ void PRDT::AnalyzeBimbam(gsl_vector *y_prdt) {
 
   // Start reading genotypes and analyze.
   for (size_t t = 0; t < ns_total; ++t) {
-    !safeGetline(infile, line).eof();
+    safeGetline(infile, line).eof();
     if (t % d_pace == 0 || t == (ns_total - 1)) {
       ProgressBar("Reading SNPs  ", t, ns_total - 1);
     }
diff --git a/src/vc.cpp b/src/vc.cpp
index f8cc2b5..f4cd650 100644
--- a/src/vc.cpp
+++ b/src/vc.cpp
@@ -663,7 +663,7 @@ void ReadFile_cor(const string &file_cor, const set<string> &setSnps,
   HEADER header;
 
   // Header.
-  !safeGetline(infile, line).eof();
+  safeGetline(infile, line).eof();
   ReadHeader_vc(line, header);
 
   if (header.n_col == 0) {
@@ -678,7 +678,7 @@ void ReadFile_cor(const string &file_cor, const set<string> &setSnps,
   while (!safeGetline(infile, line).eof()) {
 
     // do not read cor values this time; upto col_n-1.
-    ch_ptr = strtok((char *)line.c_str(), " , \t");
+    ch_ptr = strtok_safe((char *)line.c_str(), " , \t");
 
     n_total = 0;
     n_mis = 0;
@@ -688,6 +688,7 @@ void ReadFile_cor(const string &file_cor, const set<string> &setSnps,
     d_cm = 0;
     d_pos = 0;
     for (size_t i = 0; i < header.coln - 1; i++) {
+      enforce(ch_ptr);
       if (header.rs_col != 0 && header.rs_col == i + 1) {
         rs = ch_ptr;
       }
@@ -822,7 +823,7 @@ void ReadFile_beta(const bool flag_priorscale, const string &file_beta,
 
   // Read header.
   HEADER header;
-  !safeGetline(infile, line).eof();
+  safeGetline(infile, line).eof();
   ReadHeader_vc(line, header);
 
   if (header.n_col == 0) {
@@ -844,7 +845,7 @@ void ReadFile_beta(const bool flag_priorscale, const string &file_beta,
   }
 
   while (!safeGetline(infile, line).eof()) {
-    ch_ptr = strtok((char *)line.c_str(), " , \t");
+    ch_ptr = strtok_safe((char *)line.c_str(), " , \t");
 
     z = 0;
     beta = 0;
@@ -857,6 +858,7 @@ void ReadFile_beta(const bool flag_priorscale, const string &file_beta,
     af = 0;
     var_x = 0;
     for (size_t i = 0; i < header.coln; i++) {
+      enforce(ch_ptr);
       if (header.rs_col != 0 && header.rs_col == i + 1) {
         rs = ch_ptr;
       }
@@ -1055,7 +1057,7 @@ void ReadFile_cor(const string &file_cor, const vector<string> &vec_rs,
   // Header.
   HEADER header;
 
-  !safeGetline(infile, line).eof();
+  safeGetline(infile, line).eof();
   ReadHeader_vc(line, header);
 
   while (!safeGetline(infile, line).eof()) {
@@ -1063,8 +1065,9 @@ void ReadFile_cor(const string &file_cor, const vector<string> &vec_rs,
     // Do not read cor values this time; upto col_n-1.
     d_pos1 = 0;
     d_cm1 = 0;
-    ch_ptr = strtok((char *)line.c_str(), " , \t");
+    ch_ptr = strtok_safe((char *)line.c_str(), " , \t");
     for (size_t i = 0; i < header.coln - 1; i++) {
+      enforce(ch_ptr);
       if (header.rs_col != 0 && header.rs_col == i + 1) {
         rs = ch_ptr;
       }
@@ -2238,7 +2241,7 @@ bool BimbamXwz(const string &file_geno, const int display_pace,
   gsl_vector_mul(wz, w);
 
   for (size_t t = 0; t < indicator_snp.size(); ++t) {
-    !safeGetline(infile, line).eof();
+    safeGetline(infile, line).eof();
     if (t % display_pace == 0 || t == (indicator_snp.size() - 1)) {
       ProgressBar("Reading SNPs  ", t, indicator_snp.size() - 1);
     }
@@ -2246,9 +2249,9 @@ bool BimbamXwz(const string &file_geno, const int display_pace,
       continue;
     }
 
-    ch_ptr = strtok((char *)line.c_str(), " , \t");
-    ch_ptr = strtok(NULL, " , \t");
-    ch_ptr = strtok(NULL, " , \t");
+    ch_ptr = strtok_safe((char *)line.c_str(), " , \t");
+    ch_ptr = strtok_safe(NULL, " , \t");
+    ch_ptr = strtok_safe(NULL, " , \t");
 
     geno_mean = 0.0;
     n_miss = 0;
@@ -2260,7 +2263,7 @@ bool BimbamXwz(const string &file_geno, const int display_pace,
       if (indicator_idv[i] == 0) {
         continue;
       }
-      ch_ptr = strtok(NULL, " , \t");
+      ch_ptr = strtok_safe(NULL, " , \t");
       if (strcmp(ch_ptr, "NA") == 0) {
         gsl_vector_set(geno_miss, i, 0);
         n_miss++;
@@ -2491,7 +2494,7 @@ bool BimbamXtXwz(const string &file_geno, const int display_pace,
   gsl_vector *geno_miss = gsl_vector_alloc(ni_test);
 
   for (size_t t = 0; t < indicator_snp.size(); ++t) {
-    !safeGetline(infile, line).eof();
+    safeGetline(infile, line).eof();
     if (t % display_pace == 0 || t == (indicator_snp.size() - 1)) {
       ProgressBar("Reading SNPs  ", t, indicator_snp.size() - 1);
     }
@@ -2499,9 +2502,9 @@ bool BimbamXtXwz(const string &file_geno, const int display_pace,
       continue;
     }
 
-    ch_ptr = strtok((char *)line.c_str(), " , \t");
-    ch_ptr = strtok(NULL, " , \t");
-    ch_ptr = strtok(NULL, " , \t");
+    ch_ptr = strtok_safe((char *)line.c_str(), " , \t");
+    ch_ptr = strtok_safe(NULL, " , \t");
+    ch_ptr = strtok_safe(NULL, " , \t");
 
     geno_mean = 0.0;
     n_miss = 0;
@@ -2513,7 +2516,7 @@ bool BimbamXtXwz(const string &file_geno, const int display_pace,
       if (indicator_idv[i] == 0) {
         continue;
       }
-      ch_ptr = strtok(NULL, " , \t");
+      ch_ptr = strtok_safe(NULL, " , \t");
       if (strcmp(ch_ptr, "NA") == 0) {
         gsl_vector_set(geno_miss, i, 0);
         n_miss++;
diff --git a/test/dev_test_suite.sh b/test/dev_test_suite.sh
index 77604be..284c9aa 100755
--- a/test/dev_test_suite.sh
+++ b/test/dev_test_suite.sh
@@ -84,6 +84,49 @@ testUnivariateLinearMixedModelLOCO1() {
     assertEquals "15465346.22" `perl -nle 'foreach $x (split(/\s+/,$_)) { $sum += sprintf("%.2f",(substr($x,,0,6))) } END { printf "%.2f",$sum }' $outfn`
 }
 
+testPlinkCenteredRelatednessMatrixKLOCO1() {
+    return 0
+    outn=mouse_hs1940_Plink_LOCO1
+    rm -f output/$outn.*
+    $gemma -bfile ../example/mouse_hs1940 \
+           -a ../example/mouse_hs1940.anno.txt \
+           -snps ../example/mouse_hs1940_snps.txt \
+           -nind 400 \
+           -loco 1 \
+           -gk \
+           -debug \
+           -o $outn
+    assertEquals 0 $?
+    grep "total computation time" < output/$outn.log.txt
+    outfn=output/$outn.cXX.txt
+    assertEquals 0 $?
+    assertEquals "400" `wc -l < $outfn`
+    assertEquals "0.312" `head -c 5 $outfn`
+    assertEquals "71.03" `perl -nle 'foreach $x (split(/\s+/,$_)) { $sum += sprintf("%.2f",(substr($x,,0,6))) } END { printf "%.2f",$sum }' $outfn`
+}
+
+
+testPlinkUnivariateLinearMixedModelLOCO1() {
+    return 0
+    outn=mouse_hs1940_CD8_Plink_LOCO1_lmm
+    rm -f output/$outn.*
+    $gemma -bfile ../example/mouse_hs1940 \
+	   -n 1 \
+	   -loco 1 \
+           -k ./output/mouse_hs1940_Plink_LOCO1.cXX.txt \
+           -a ../example/mouse_hs1940.anno.txt \
+	   -snps ../example/mouse_hs1940_snps.txt -lmm \
+	   -nind 400 \
+	   -debug \
+           -o $outn
+    assertEquals 0 $?
+    grep "total computation time" < output/$outn.log.txt
+    assertEquals 0 $?
+    outfn=output/$outn.assoc.txt
+    assertEquals "68" `wc -l < $outfn`
+    assertEquals "15465346.22" `perl -nle 'foreach $x (split(/\s+/,$_)) { $sum += sprintf("%.2f",(substr($x,,0,6))) } END { printf "%.2f",$sum }' $outfn`
+}
+
 shunit2=`which shunit2`
 
 if [ -x "$shunit2" ]; then
generated by cgit-pink 1.4.1 (git 2.36.1) at 2025-10-14 05:03:08 +0000