aboutsummaryrefslogtreecommitdiff
path: root/test
diff options
context:
space:
mode:
Diffstat (limited to 'test')
-rwxr-xr-xtest/dev_test_suite.sh79
-rw-r--r--test/src/unittests-math.cpp117
-rwxr-xr-xtest/test_suite.sh111
3 files changed, 275 insertions, 32 deletions
diff --git a/test/dev_test_suite.sh b/test/dev_test_suite.sh
index 0fc4423..0d3d8a0 100755
--- a/test/dev_test_suite.sh
+++ b/test/dev_test_suite.sh
@@ -1,29 +1,31 @@
#!/usr/bin/env bash
gemma=../bin/gemma
+# gemmaopts="-debug -strict"
+gemmaopts="-debug"
# Related to https://github.com/genetics-statistics/GEMMA/issues/78
testBXDStandardRelatednessMatrixKSingularError() {
outn=BXDerr
rm -f output/$outn.*
- $gemma -g ../example/BXD_geno.txt.gz \
+ $gemma $gemmaopts \
+ -g ../example/BXD_geno.txt.gz \
-p ../example/BXD_pheno.txt \
-c ../example/BXD_covariates.txt \
-a ../example/BXD_snps.txt \
-gk \
- -debug -o $outn
+ -o $outn
assertEquals 22 $? # should show singular error
}
testBXDStandardRelatednessMatrixK() {
outn=BXD
rm -f output/$outn.*
- $gemma -g ../example/BXD_geno.txt.gz \
+ $gemma $gemmaopts -g ../example/BXD_geno.txt.gz \
-p ../example/BXD_pheno.txt \
-c ../example/BXD_covariates2.txt \
-a ../example/BXD_snps.txt \
-gk \
- -debug \
-o $outn
assertEquals 0 $?
outfn=output/$outn.cXX.txt
@@ -31,28 +33,43 @@ testBXDStandardRelatednessMatrixK() {
assertEquals "-116.11" `perl -nle 'foreach $x (split(/\s+/,$_)) { $sum += sprintf("%.2f",(substr($x,,0,6))) } END { printf "%.2f",$sum }' $outfn`
}
+testBXDLMLikelihoodRatio() {
+ outn=BXD_LM_LR
+ $gemma $gemmaopts -g ../example/BXD_geno.txt.gz \
+ -p ../example/BXD_pheno.txt \
+ -c ../example/BXD_covariates2.txt \
+ -a ../example/BXD_snps.txt \
+ -k ./output/BXD.cXX.txt \
+ -lm 4 -maf 0.1 \
+ -o $outn
+ assertEquals 0 $?
+
+ outfn=output/$outn.assoc.txt
+ assertEquals "95134" `wc -w < $outfn`
+ assertEquals "3089042886.28" `perl -nle 'foreach $x (split(/\s+/,$_)) { $sum += sprintf("%.2f",(substr($x,,0,6))) } END { printf "%.2f",$sum }' $outfn`
+}
+
testBXDLMMLikelihoodRatio() {
outn=BXD_LMM_LR
- $gemma -g ../example/BXD_geno.txt.gz \
+ $gemma $gemmaopts -g ../example/BXD_geno.txt.gz \
-p ../example/BXD_pheno.txt \
-c ../example/BXD_covariates2.txt \
-a ../example/BXD_snps.txt \
-k ./output/BXD.cXX.txt \
-lmm 2 -maf 0.1 \
- -debug \
-o $outn
assertEquals 0 $?
outfn=output/$outn.assoc.txt
- assertEquals "80498" `wc -w < $outfn`
+ assertEquals "73180" `wc -w < $outfn`
assertEquals "3088458212.93" `perl -nle 'foreach $x (split(/\s+/,$_)) { $sum += sprintf("%.2f",(substr($x,,0,6))) } END { printf "%.2f",$sum }' $outfn`
}
testCenteredRelatednessMatrixKLOCO1() {
outn=mouse_hs1940_LOCO1
rm -f output/$outn.*
- $gemma -g ../example/mouse_hs1940.geno.txt.gz -p ../example/mouse_hs1940.pheno.txt \
- -a ../example/mouse_hs1940.anno.txt -snps ../example/mouse_hs1940_snps.txt -nind 400 -loco 1 -gk -debug -o $outn
+ $gemma $gemmaopts -g ../example/mouse_hs1940.geno.txt.gz -p ../example/mouse_hs1940.pheno.txt \
+ -a ../example/mouse_hs1940.anno.txt -snps ../example/mouse_hs1940_snps.txt -nind 400 -loco 1 -gk -o $outn
assertEquals 0 $?
grep "total computation time" < output/$outn.log.txt
outfn=output/$outn.cXX.txt
@@ -65,7 +82,7 @@ testCenteredRelatednessMatrixKLOCO1() {
testUnivariateLinearMixedModelLOCO1() {
outn=mouse_hs1940_CD8_LOCO1_lmm
rm -f output/$outn.*
- $gemma -g ../example/mouse_hs1940.geno.txt.gz \
+ $gemma $gemmaopts -g ../example/mouse_hs1940.geno.txt.gz \
-p ../example/mouse_hs1940.pheno.txt \
-n 1 \
-loco 1 \
@@ -73,7 +90,47 @@ testUnivariateLinearMixedModelLOCO1() {
-k ./output/mouse_hs1940_LOCO1.cXX.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`
+}
+
+testPlinkCenteredRelatednessMatrixKLOCO1() {
+ return 0
+ outn=mouse_hs1940_Plink_LOCO1
+ rm -f output/$outn.*
+ $gemma $gemmaopts -bfile ../example/mouse_hs1940 \
+ -a ../example/mouse_hs1940.anno.txt \
+ -snps ../example/mouse_hs1940_snps.txt \
+ -nind 400 \
+ -loco 1 \
+ -gk \
+ -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 $gemmaopts -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 \
-o $outn
assertEquals 0 $?
grep "total computation time" < output/$outn.log.txt
diff --git a/test/src/unittests-math.cpp b/test/src/unittests-math.cpp
index ac4c180..757c2dc 100644
--- a/test/src/unittests-math.cpp
+++ b/test/src/unittests-math.cpp
@@ -1,14 +1,23 @@
#include <catch.hpp>
#include <iostream>
#include "gsl/gsl_matrix.h"
-#include "mathfunc.h"
+#include <cblas.h>
+
#include <algorithm>
#include <limits>
#include <numeric>
+#include "debug.h"
+#include "mathfunc.h"
+#include "fastblas.h"
+#include "fastopenblas.h"
+
using namespace std;
TEST_CASE( "Math functions", "[math]" ) {
+ debug_set_debug_mode(true);
+ debug_set_no_check_mode(false);
+ debug_set_strict_mode(true);
double data[] = { 2,-1, 0,
-1, 2,-1,
0,-1, 2};
@@ -51,3 +60,109 @@ TEST_CASE( "Math functions", "[math]" ) {
REQUIRE (std::isnan(v3[2]));
REQUIRE(has_nan(v3));
}
+
+TEST_CASE("cblas_dgemm", "[math]") {
+ double *A, *B, *C;
+ int m, n, k, i, j;
+ double alpha, beta;
+
+ printf ("\n This example computes real matrix C=alpha*A*B+beta*C using \n"
+ " Intel(R) MKL function dgemm, where A, B, and C are matrices and \n"
+ " alpha and beta are double precision scalars\n\n");
+
+ m = 2000, k = 200, n = 1000;
+ printf (" Initializing data for matrix multiplication C=A*B for matrix \n"
+ " A(%ix%i) and matrix B(%ix%i)\n\n", m, k, k, n);
+ alpha = 1.0; beta = 0.0;
+
+ printf (" Allocating memory for matrices aligned on 64-byte boundary for better \n"
+ " performance \n\n");
+ A = (double *)malloc( m*k*sizeof( double ));
+ B = (double *)malloc( k*n*sizeof( double ));
+ C = (double *)malloc( m*n*sizeof( double ));
+
+ printf (" Intializing matrix data \n\n");
+ for (i = 0; i < (m*k); i++) {
+ A[i] = (double)(i+1);
+ }
+
+ for (i = 0; i < (k*n); i++) {
+ B[i] = (double)(-i-1);
+ }
+
+ for (i = 0; i < (m*n); i++) {
+ C[i] = 0.0;
+ }
+
+ printf (" Computing matrix product using Intel(R) MKL dgemm function via CBLAS interface \n\n");
+ assert(m==2000);
+ assert(k==200);
+ assert(n==1000);
+ //cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
+ // m, n, k, alpha, A, k, B, n, beta, C, n);
+ fast_cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
+ m, n, k, alpha, A, k, B, n, beta, C, n);
+
+ REQUIRE(trunc(C[0]) == -2666620100.0 );
+ REQUIRE(trunc(C[1]) == -2666640200.0 );
+ REQUIRE(trunc(C[2003]) == -10627000400.0 );
+
+}
+
+TEST_CASE("fast_dgemm", "[math]") {
+ double *A, *B, *C;
+ int m, n, k, i, j;
+ double alpha, beta;
+
+ printf ("\n This example computes real matrix C=alpha*A*B+beta*C using \n"
+ " Intel(R) MKL function dgemm, where A, B, and C are matrices and \n"
+ " alpha and beta are double precision scalars\n\n");
+
+ m = 2000, k = 200, n = 1000;
+ printf (" Initializing data for matrix multiplication C=A*B for matrix \n"
+ " A(%ix%i) and matrix B(%ix%i)\n\n", m, k, k, n);
+ alpha = 1.0; beta = 0.0;
+
+ printf (" Allocating memory for matrices aligned on 64-byte boundary for better \n"
+ " performance \n\n");
+ A = (double *)malloc( m*k*sizeof( double ));
+ B = (double *)malloc( k*n*sizeof( double ));
+ C = (double *)malloc( m*n*sizeof( double ));
+
+ printf (" Intializing matrix data \n\n");
+ for (i = 0; i < (m*k); i++) {
+ A[i] = (double)(i+1);
+ }
+
+ for (i = 0; i < (k*n); i++) {
+ B[i] = (double)(-i-1);
+ }
+
+ for (i = 0; i < (m*n); i++) {
+ C[i] = 0.0;
+ }
+
+ printf (" Computing matrix product using Intel(R) MKL dgemm function via CBLAS interface \n\n");
+ // cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
+ // m, n, k, alpha, A, k, B, n, beta, C, n);
+ // eigenlib_dgemm(const char *TransA, const char *TransB, const double alpha,
+ // const gsl_matrix *A, const gsl_matrix *B, const double beta,
+ // gsl_matrix *C) {
+ gsl_matrix *AM = gsl_matrix_safe_alloc(m,k); // rows x cols
+ gsl_matrix *BM = gsl_matrix_safe_alloc(k,n);
+ gsl_matrix *CM = gsl_matrix_calloc(m,n);
+
+ fast_copy(AM,A);
+ fast_copy(BM,B);
+ fast_copy(CM,C);
+ fast_dgemm("N","N",alpha,AM,BM,beta,CM);
+ printf ("\n Computations completed.\n\n");
+ A = AM->data;
+ B = BM->data;
+ C = CM->data;
+
+ REQUIRE(trunc(C[0]) == -2666620100.0 );
+ REQUIRE(trunc(C[1]) == -2666640200.0 );
+ REQUIRE(trunc(C[2003]) == -10627000400.0 );
+
+}
diff --git a/test/test_suite.sh b/test/test_suite.sh
index 350fc27..7af33aa 100755
--- a/test/test_suite.sh
+++ b/test/test_suite.sh
@@ -1,13 +1,89 @@
#!/usr/bin/env bash
gemma=../bin/gemma
+gemmaopts="-debug"
+
+testBslmm1() {
+ outn=mouse_hs1940_CD8_bslmm
+ $gemma $gemmaopts -g ../example/mouse_hs1940.geno.txt.gz \
+ -p ../example/mouse_hs1940.pheno.txt \
+ -n 2 -a ../example/mouse_hs1940.anno.txt \
+ -bslmm \
+ -o $outn -w 1000 -s 10000 -seed 1
+ assertEquals 0 $?
+ outfn1=output/$outn.hyp.txt
+ outfn2=output/$outn.param.txt
+ # assertEquals "45181" `perl -nle 'foreach $x (split(/\s+/,$_)) { $sum += sprintf("%.0f",(substr($x,,0,6))) } END { printf "%.0f",$sum }' $outfn1`
+ # assertEquals "4043967139.42" `perl -nle 'foreach $x (split(/\s+/,$_)) { $sum += sprintf("%.2f",(substr($x,,0,6))) } END { printf "%.2f",$sum }' $outfn2`
+}
+
+testBslmm2() {
+ outn=mouse_hs1940_CD8_train
+ $gemma $gemmaopts -g ../example/mouse_hs1940.geno.txt.gz \
+ -p ../example/mouse_hs1940.pheno.txt \
+ -n 2 \
+ -a ../example/mouse_hs1940.anno.txt \
+ -gk 1 -o $outn
+ assertEquals 0 $?
+ outfn=output/$outn.cXX.txt
+ assertEquals "579.66" `perl -nle 'foreach $x (split(/\s+/,$_)) { $sum += sprintf("%.2f",(substr($x,,0,6))) } END { printf "%.2f",$sum }' $outfn`
+}
+
+testBslmm3() {
+ ## Fit a binary trait using a linear model
+ outn=mouse_hs1940_CD8_bslmm_cc1
+ $gemma $gemmaopts \
+ -g ../example/mouse_hs1940.geno.txt.gz \
+ -p ../example/mouse_hs1940.pheno.txt \
+ -n 4 \
+ -a ../example/mouse_hs1940.anno.txt \
+ -bslmm \
+ -o $outn \
+ -w 1000 -s 10000 -seed 1
+ assertEquals 0 $?
+ outfn=output/$outn.hyp.txt
+ # assertEquals "291" `perl -nle 'foreach $x (split(/\s+/,$_)) { $sum += sprintf("%.0f",(substr($x,,0,6))) } END { printf "%.0f",100*$sum }' $outfn`
+}
+
+testBslmm4() {
+ outn=mouse_hs1940_CD8_prdt_k
+ $gemma $gemmaopts -g ../example/mouse_hs1940.geno.txt.gz \
+ -p ../example/mouse_hs1940.pheno.txt \
+ -n 2 \
+ -epm ./output/mouse_hs1940_CD8_bslmm.param.txt \
+ -emu ./output/mouse_hs1940_CD8_bslmm.log.txt \
+ -ebv ./output/mouse_hs1940_CD8_bslmm.bv.txt \
+ -k ./output/mouse_hs1940_CD8_train.cXX.txt \
+ -predict \
+ -o $outn
+ assertEquals 0 $?
+ outfn=output/$outn.prdt.txt
+ # assertEquals "-60.33" `perl -nle 'foreach $x (split(/\s+/,$_)) { $sum += sprintf("%.2f",(substr($x,,0,6))) } END { printf "%.2f",$sum }' $outfn`
+}
+
+testBslmm5() {
+ ## Now, do prediction in the test set for the binary traits
+ ## If the traits were fitted using the linear model, then:
+ outn=mouse_hs1940_CD8_prdt_cc1
+ $gemma $gemmaopts \
+ -g ../example/mouse_hs1940.geno.txt.gz \
+ -p ../example/mouse_hs1940.pheno.txt \
+ -n 4 \
+ -epm ./output/mouse_hs1940_CD8_bslmm_cc1.param.txt \
+ -emu ./output/mouse_hs1940_CD8_bslmm_cc1.log.txt \
+ -predict \
+ -o $outn
+ assertEquals 0 $?
+ outfn=output/$outn.prdt.txt
+ assertEquals "550.67" `perl -nle 'foreach $x (split(/\s+/,$_)) { $sum += sprintf("%.2f",(substr($x,,0,6))) } END { printf "%.2f",$sum }' $outfn`
+}
testCenteredRelatednessMatrixKFullLOCO1() {
outn=mouse_hs1940_full_LOCO1
- $gemma -g ../example/mouse_hs1940.geno.txt.gz \
+ $gemma $gemmaopts -g ../example/mouse_hs1940.geno.txt.gz \
-p ../example/mouse_hs1940.pheno.txt \
-a ../example/mouse_hs1940.anno.txt \
- -loco 1 -gk -debug -o $outn
+ -loco 1 -gk -o $outn
assertEquals 0 $?
outfn=output/$outn.cXX.txt
assertEquals "1940" `wc -l < $outfn`
@@ -16,14 +92,13 @@ testCenteredRelatednessMatrixKFullLOCO1() {
testUnivariateLinearMixedModelFullLOCO1() {
outn=mouse_hs1940_CD8_full_LOCO1_lmm
- $gemma -g ../example/mouse_hs1940.geno.txt.gz \
+ $gemma $gemmaopts -g ../example/mouse_hs1940.geno.txt.gz \
-p ../example/mouse_hs1940.pheno.txt \
-n 1 \
-loco 1 \
-a ../example/mouse_hs1940.anno.txt \
-k ./output/mouse_hs1940_full_LOCO1.cXX.txt \
-lmm \
- -debug \
-o $outn
assertEquals 0 $?
grep "total computation time" < output/$outn.log.txt
@@ -34,9 +109,9 @@ testUnivariateLinearMixedModelFullLOCO1() {
}
testCenteredRelatednessMatrixK() {
- $gemma -g ../example/mouse_hs1940.geno.txt.gz \
+ $gemma $gemmaopts -g ../example/mouse_hs1940.geno.txt.gz \
-p ../example/mouse_hs1940.pheno.txt \
- -gk -o mouse_hs1940 -debug
+ -gk -o mouse_hs1940
assertEquals 0 $?
outfn=output/mouse_hs1940.cXX.txt
assertEquals "1940" `wc -l < $outfn`
@@ -46,14 +121,13 @@ testCenteredRelatednessMatrixK() {
}
testUnivariateLinearMixedModel() {
- $gemma -g ../example/mouse_hs1940.geno.txt.gz \
+ $gemma $gemmaopts -g ../example/mouse_hs1940.geno.txt.gz \
-p ../example/mouse_hs1940.pheno.txt \
-n 1 \
-a ../example/mouse_hs1940.anno.txt \
-k ./output/mouse_hs1940.cXX.txt \
-lmm \
- -o mouse_hs1940_CD8_lmm \
- -debug
+ -o mouse_hs1940_CD8_lmm
assertEquals 0 $?
grep "total computation time" < output/mouse_hs1940_CD8_lmm.log.txt
assertEquals 0 $?
@@ -62,14 +136,13 @@ testUnivariateLinearMixedModel() {
assertEquals "4038540440.86" `perl -nle 'foreach $x (split(/\s+/,$_)) { $sum += sprintf("%.2f",(substr($x,,0,6))) } END { printf "%.2f",$sum }' $outfn`
}
-testMultivariateLinearMixedModel() {
- $gemma -g ../example/mouse_hs1940.geno.txt.gz \
+testLinearMixedModelPhenotypes() {
+ $gemma $gemmaopts -g ../example/mouse_hs1940.geno.txt.gz \
-p ../example/mouse_hs1940.pheno.txt \
-n 1 6 \
-a ../example/mouse_hs1940.anno.txt \
-k ./output/mouse_hs1940.cXX.txt \
- -lmm -o mouse_hs1940_CD8MCH_lmm \
- -debug
+ -lmm -o mouse_hs1940_CD8MCH_lmm
assertEquals 0 $?
outfn=output/mouse_hs1940_CD8MCH_lmm.assoc.txt
@@ -82,9 +155,8 @@ testPlinkStandardRelatednessMatrixK() {
datadir=../example
outfn=output/$testname.sXX.txt
rm -f $outfn
- $gemma -bfile $datadir/HLC \
- -gk 2 -o $testname \
- -debug
+ $gemma $gemmaopts -bfile $datadir/HLC \
+ -gk 2 -o $testname
assertEquals 0 $?
assertEquals "427" `wc -l < $outfn`
assertEquals "-358.07" `perl -nle 'foreach $x (split(/\s+/,$_)) { $sum += sprintf("%.2f",(substr($x,,0,6))) } END { printf "%.2f",$sum }' $outfn`
@@ -92,15 +164,14 @@ testPlinkStandardRelatednessMatrixK() {
# Test for https://github.com/genetics-statistics/GEMMA/issues/58
# fixed GSLv2 NaN's that appeared with covariates.
-testPlinkMultivariateLinearMixedModel() {
- testname=testPlinkMultivariateLinearMixedModel
+testPlinkLinearMixedModelCovariates() {
+ testname=testPlinkLinearMixedModelCovariates
datadir=../example
- $gemma -bfile $datadir/HLC \
+ $gemma $gemmaopts -bfile $datadir/HLC \
-k output/testPlinkStandardRelatednessMatrixK.sXX.txt \
-lmm 1 \
-maf 0.1 \
-c $datadir/HLC_covariates.txt \
- -debug \
-o $testname
assertEquals 0 $?
outfn=output/$testname.assoc.txt
generated by cgit v1.2.3 (git 2.47.1) at 2025-07-02 04:42:00 +0000