aboutsummaryrefslogtreecommitdiff
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 v1.2.3 (git 2.44.0) at 2025-01-18 19:32:25 +0000