From b42a02d02b3d9384b1da55bd091f0f89c808b626 Mon Sep 17 00:00:00 2001 From: Pjotr Prins Date: Fri, 20 Oct 2017 08:10:39 +0000 Subject: Travis-ci: - Disabled gcc-6 since we develop with later tools anyway - Turned the release test into a simple integration test - Adding MacOSX on Travis-ci Tests: Adding tests for Plink w. LOCO Safety: Introduce strtok_safe to get rid of segfaults --- .travis.yml | 58 ++++++++++-------- Makefile | 23 ++++--- src/bslmmdap.cpp | 19 +++--- src/debug.cpp | 11 +++- src/debug.h | 15 ++++- src/io.cpp | 160 ++++++++++++++++++++++++------------------------- src/lm.cpp | 12 ++-- src/lmm.cpp | 20 +++---- src/mvlmm.cpp | 20 +++---- src/param.cpp | 8 +-- src/prdt.cpp | 2 +- src/vc.cpp | 35 ++++++----- test/dev_test_suite.sh | 43 +++++++++++++ 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 &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 &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 &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 &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 &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 &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 &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) { 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> &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 &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 &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 &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 &indicator_idv, vector 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 &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 &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 &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 mapID2ID; @@ -1193,11 +1184,11 @@ void ReadFile_kin(const string &file_kin, vector &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 &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 &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 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 &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 &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 &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 &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 &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 &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 &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 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 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 &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 &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 &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 &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 . */ +#include +#include #include #include #include #include -#include -#include #include #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 &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 &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 &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 &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 &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 -- cgit v1.2.3