2 files changed, 0 insertions, 218 deletions
diff --git a/topics/lmms/GSL_ERROR_function_value_is_not_finite_in_brent.c-issue210.gmi b/topics/lmms/GSL_ERROR_function_value_is_not_finite_in_brent.c-issue210.gmi
deleted file mode 100644
index 23ad3e5..0000000
--- a/topics/lmms/GSL_ERROR_function_value_is_not_finite_in_brent.c-issue210.gmi
+++ /dev/null
@@ -1,93 +0,0 @@
-# GEMMA GSL ERROR: function value is not finite in brent.c
-
-In github issue 210 this comes up. @hrwang was kind enough to send me a dataset that brings up this error. Now we stand a chance of doing something about it! As in
-
-```
-napoli:~/iwrk/opensource/code/genetics/gemma/tmp/issue210$ ../../bin/gemma -bfile 20 -k BS_kinship.sXX.txt -lmm 1 -c Covariate_PC3.txt -n 1 -o outname
-GEMMA 0.99.0-pre1 (2021-08-11) by Xiang Zhou, Pjotr Prins and team (C) 2012-2021
-Reading Files ...
-## number of total individuals = 727
-## number of analyzed individuals = 57
-## number of covariates = 4
-## number of phenotypes = 1
-## number of total SNPs/var        =   215451
-## number of analyzed SNPs         =    72731
-Start Eigen-Decomposition...
-pve estimate =0.956403
-se(pve) =0.255427
-GSL ERROR: function value is not finite in newton.c at line 88 errno 9
-```
-
-Next I ran with the debug switch and it stops with the slightly more informative
-
-```
-../../bin/gemma -bfile 20 -k BS_kinship.sXX.txt -lmm 1 -c Covariate_PC3.txt -n 1 -o outname -debug
-GEMMA 0.99.0-pre1 (2021-08-11) by Xiang Zhou, Pjotr Prins and team (C) 2012-2021
-GSL random generator type: mt19937; seed = 25479 (option -1); first value = 191757793
-Reading Files ...
-## number of total individuals = 727
-## number of analyzed individuals = 57
-## number of covariates = 4
-## number of phenotypes = 1
-## number of total SNPs/var        =   215451
-## number of analyzed SNPs         =    72731
-Start Eigen-Decomposition...
-pve estimate =0.956403
-se(pve) =0.255427
-ERROR: Enforce failed for Trying to take the log of 0.000000 in src/mathfunc.cpp at line 118 in safe_log
-```
-
-Next step is rebuilding gemma with debug info `make debug` and rerun in the debugger. A backtrace of
-
-```
-gdb --args ../../bin/gemma -bfile 20 -k BS_kinship.sXX.txt -lmm 1 -c Covariate_PC3.txt -n 1 -o outname -debug
-```
-
-shows
-
-```
-#2  0x000000000047dcb0 in safe_log (d=d@entry=0) at src/mathfunc.cpp:118
-#3  0x00000000004707ea in LogRL_f (l=l@entry=1.0000000000000001e-05, params=params@entry=0x7fffffffcc20)
-    at src/lmm.cpp:851
-#4  0x0000000000471b05 in CalcLambda (func_name=func_name@entry=82 'R', params=...,
-    l_min=l_min@entry=1.0000000000000001e-05, l_max=100000, n_region=10,
-    lambda=@0x7fffffffccc8: 3.2054654756518093, logf=@0x7fffffffcca0: 0.83040722111380116)
-    at src/lmm.cpp:1970
-#5  0x0000000000475cac in LMM::AnalyzePlink (this=this@entry=0x7fffffffd1a0, U=U@entry=0x4fcc80,
-    eval=eval@entry=0x52dff0, UtW=UtW@entry=0x4fccc0, Uty=Uty@entry=0x7fffffffd140, W=W@entry=0x4fe140,
-    y=<optimized out>, gwasnps=...) at src/lmm.cpp:1860
-#6  0x00000000004359f5 in GEMMA::BatchRun (this=this@entry=0x7fffffffe720, cPar=...)
-    at src/gemma.cpp:2798
-#7  0x000000000047bed0 in main (argc=14, argv=0x7fffffffe888) at src/main.cpp:86
-```
-
-In this case a SNP for some reason balks. Another run, after disabling the log check renders
-
-```
-#0  0x00007ffff5e1deca in raise ()
-   from /gnu/store/fa6wj5bxkj5ll1d7292a70knmyl7a0cr-glibc-2.31/lib/libpthread.so.0
-#1  0x00000000004245ba in gemma_gsl_error_handler (reason=<optimized out>, file=<optimized out>,
-    line=88, gsl_errno=<optimized out>) at src/gemma.cpp:75
-#2  0x00007ffff7e59f6d in newton_iterate ()
-   from /gnu/store/qsr7pw4cbka4qhxm3h0bkwwl7h08qfmw-profile/lib/libgsl.so.25
-#3  0x0000000000471cba in CalcLambda (func_name=func_name@entry=82 'R', params=...,
-    l_min=l_min@entry=1.0000000000000001e-05, l_max=100000, n_region=<optimized out>,
-    lambda=@0x7fffffffccc8: 1.0000000000000001e-05, logf=@0x7fffffffcca0: 507.88508012472289)
-    at src/lmm.cpp:2040
-```
-
-points out that in line 2040 `gsl_root_fdfsolver_iterate` is not happy executing a Newton root solver. What is interesting is that, in this case, the gemma_gsl_error_handler is called and stops processing. We now have three ways to stop GEMMA from the same error!
-
-The original version bails out inside GSL:
-
-```
-GSL ERROR: function value is not finite in newton.c at line 88 errno 9
-
-#0  0x00007ffff5e1deca in raise ()
-   from /gnu/store/fa6wj5bxkj5ll1d7292a70knmyl7a0cr-glibc-2.31/lib/libpthread.so.0
-#1  0x00007ffff7e59f6d in newton_iterate ()
-   from /gnu/store/qsr7pw4cbka4qhxm3h0bkwwl7h08qfmw-profile/lib/libgsl.so.25
-#2  0x000000000047a9c5 in CalcLambda(char, FUNC_PARAM&, double, double, unsigned long, double&, double&) [clone .constprop.671] ()
-```
-
-@hrwang had a good suggestion to make sure that the output for the SNP is made invalid and to continue processing. Disabling the error handler let GEMMA continue. I added a check for the status flag to make sure output was set to invalid (nan). In the test set 178 out of 72732 SNPs failed.
diff --git a/topics/lmms/issue243.gmi b/topics/lmms/issue243.gmi
deleted file mode 100644
index 255de8b..0000000
--- a/topics/lmms/issue243.gmi
+++ /dev/null
@@ -1,125 +0,0 @@
-# Multivariate GEMMA hangs (issue 243)
-
-The following command does finish
-
-```
-time ../../bin/gemma -bfile multivariate_2traits -k output/multivariate.cXX.txt -maf 0.0000001 -lmm 4 -n 1 2 -o gemma.polygenic.result
-real    133m29.533s
-user    133m31.746s
-sys     0m5.064s
-```
-
-But is slow. A print statement shows close to a second is spent on every SNP:
-
-```
-0.83 9994,10000
-0.82 9995,10000
-0.81 9996,10000
-0.80 9997,10000
-0.80 9998,10000
-0.80 9999,10000
-```
-
-Running the debugger on and hitting Ctrl-C repeatedly stops in the following places
-
-```
-#0  0x00000000004807f9 in MphCalcLogL (eval=eval@entry=0x5233c0, xHiy=xHiy@entry=0x59e5d0,
-    D_l=D_l@entry=0x5ac840, UltVehiY=UltVehiY@entry=0x53e500, Qi=Qi@entry=0x59e2e0) at src/mvlmm.cpp:579
-#1  0x000000000048104b in MphEM (func_name=func_name@entry=82 'R', max_iter=1000, max_prec=0.001,
-    eval=eval@entry=0x5233c0, X=X@entry=0x540930, Y=Y@entry=0x5408d0, U_hat=U_hat@entry=0x53e380,
-    E_hat=0x53e3e0, OmegaU=0x53e440, OmegaE=0x53e4a0, UltVehiY=0x53e500, UltVehiBX=0x53e560,
-    UltVehiU=0x53e5c0, UltVehiE=0x53e620, V_g=0x540990, V_e=0x540a20, B=0x540ab0) at src/mvlmm.cpp:662
-#2  0x000000000048f924 in MVLMM::AnalyzePlink (this=this@entry=0x7fffffffcf20, U=U@entry=0x523380,
-    eval=eval@entry=0x5233c0, UtW=UtW@entry=0x5233f0, UtY=UtY@entry=0x523430) at src/mvlmm.cpp:3803
-#3  0x0000000000435c2e in GEMMA::BatchRun (this=this@entry=0x7fffffffe4a0, cPar=...)
-    at src/gemma.cpp:2830
-
-#0  0x00007ffff7df8690 in gsl_matrix_sub ()
-   from /gnu/store/pmmbmrbizz668plrrvfyvqh5ymvjy5p4-profile/lib/libgsl.so.25
-#1  0x000000000047fc91 in UpdateU (OmegaE=OmegaE@entry=0x53e4a0, UltVehiY=UltVehiY@entry=0x53e500,
-    UltVehiBX=UltVehiBX@entry=0x53e560, UltVehiU=UltVehiU@entry=0x53e5c0) at src/mvlmm.cpp:387
-#2  0x0000000000480e76 in MphEM (func_name=func_name@entry=76 'L', max_iter=1000, max_prec=0.001,
-    eval=eval@entry=0x5233c0, X=X@entry=0x540930, Y=Y@entry=0x5408d0, U_hat=U_hat@entry=0x53e380,
-    E_hat=0x53e3e0, OmegaU=0x53e440, OmegaE=0x53e4a0, UltVehiY=0x53e500, UltVehiBX=0x53e560,
-    UltVehiU=0x53e5c0, UltVehiE=0x53e620, V_g=0x540990, V_e=0x540a20, B=0x540ab0) at src/mvlmm.cpp:686
-#3  0x000000000048faed in MVLMM::AnalyzePlink (this=this@entry=0x7fffffffcf20, U=U@entry=0x523380,
-    eval=eval@entry=0x5233c0, UtW=UtW@entry=0x5233f0, UtY=UtY@entry=0x523430) at src/mvlmm.cpp:3778
-#4  0x0000000000435c2e in GEMMA::BatchRun (this=this@entry=0x7fffffffe4a0, cPar=...)
-
-at src/gemma.cpp:2830
-
-#0  0x00007ffff6300718 in dgemm_nt ()
-   from /gnu/store/pmmbmrbizz668plrrvfyvqh5ymvjy5p4-profile/lib/libopenblas.so.0
-#1  0x00007ffff622e6f3 in cblas_dgemm ()
-   from /gnu/store/pmmbmrbizz668plrrvfyvqh5ymvjy5p4-profile/lib/libopenblas.so.0
-#2  0x00007ffff7d6af44 in gsl_blas_dgemm ()
-   from /gnu/store/pmmbmrbizz668plrrvfyvqh5ymvjy5p4-profile/lib/libgsl.so.25
-#3  0x00000000004806a4 in CalcSigma (func_name=func_name@entry=82 'R', eval=eval@entry=0x5233c0,
-    D_l=D_l@entry=0x59e5d0, X=X@entry=0x540930, OmegaU=OmegaU@entry=0x53e440,
-    OmegaE=OmegaE@entry=0x53e4a0, UltVeh=0x5abd30, Qi=0x5abe50, Sigma_uu=0x540fb0, Sigma_ee=0x5ac7e0)
-    at src/mvlmm.cpp:542
-#4  0x0000000000480f62 in MphEM (func_name=func_name@entry=82 'R', max_iter=1000, max_prec=0.001,
-    eval=eval@entry=0x5233c0, X=X@entry=0x540930, Y=Y@entry=0x5408d0, U_hat=U_hat@entry=0x53e380,
-    E_hat=0x53e3e0, OmegaU=0x53e440, OmegaE=0x53e4a0, UltVehiY=0x53e500, UltVehiBX=0x53e560,
-    UltVehiU=0x53e5c0, UltVehiE=0x53e620, V_g=0x540990, V_e=0x540a20, B=0x540ab0) at src/mvlmm.cpp:704
-#5  0x000000000048f924 in MVLMM::AnalyzePlink (this=this@entry=0x7fffffffcf20, U=U@entry=0x523380,
-    eval=eval@entry=0x5233c0, UtW=UtW@entry=0x5233f0, UtY=UtY@entry=0x523430) at src/mvlmm.cpp:3803
-#6  0x0000000000435c2e in GEMMA::BatchRun (this=this@entry=0x7fffffffe4a0, cPar=...)
-    at src/gemma.cpp:2830
-
-#0  0x00007ffff5e16f00 in __pthread_mutex_unlock_usercnt ()
-   from /gnu/store/fa6wj5bxkj5ll1d7292a70knmyl7a0cr-glibc-2.31/lib/libpthread.so.0
-#1  0x00007ffff622e6fb in cblas_dgemm ()
-   from /gnu/store/pmmbmrbizz668plrrvfyvqh5ymvjy5p4-profile/lib/libopenblas.so.0
-#2  0x00007ffff7d6af44 in gsl_blas_dgemm ()
-   from /gnu/store/pmmbmrbizz668plrrvfyvqh5ymvjy5p4-profile/lib/libgsl.so.25
-#3  0x000000000048067b in CalcSigma (func_name=func_name@entry=82 'R', eval=eval@entry=0x5233c0,
-    D_l=D_l@entry=0x5ac970, X=X@entry=0x540930, OmegaU=OmegaU@entry=0x53e440,
-    OmegaE=OmegaE@entry=0x53e4a0, UltVeh=0x5abdb0, Qi=0x59eb10, Sigma_uu=0x5ac3e0, Sigma_ee=0x5ab290)
-    at src/mvlmm.cpp:541
-
-#0  0x00007ffff5e16f19 in __pthread_mutex_unlock_usercnt ()
-   from /gnu/store/fa6wj5bxkj5ll1d7292a70knmyl7a0cr-glibc-2.31/lib/libpthread.so.0
-#1  0x00007ffff643d667 in blas_memory_alloc ()
-   from /gnu/store/pmmbmrbizz668plrrvfyvqh5ymvjy5p4-profile/lib/libopenblas.so.0
-#2  0x00007ffff622e648 in cblas_dgemm ()
-   from /gnu/store/pmmbmrbizz668plrrvfyvqh5ymvjy5p4-profile/lib/libopenblas.so.0
-#3  0x00007ffff7d6af44 in gsl_blas_dgemm ()
-   from /gnu/store/pmmbmrbizz668plrrvfyvqh5ymvjy5p4-profile/lib/libgsl.so.25
-#4  0x0000000000480652 in CalcSigma (func_name=func_name@entry=82 'R', eval=eval@entry=0x5233c0,
-    D_l=D_l@entry=0x5ab1e0, X=X@entry=0x540930, OmegaU=OmegaU@entry=0x53e440,
-    OmegaE=OmegaE@entry=0x53e4a0, UltVeh=0x5a0990, Qi=0x59e570, Sigma_uu=0x59e2e0, Sigma_ee=0x5abe30)
-    at src/mvlmm.cpp:539
-
-#0  0x000000000047f8e7 in gsl_matrix_get (j=389, i=0, m=0x53e500)
-    at /gnu/store/pmmbmrbizz668plrrvfyvqh5ymvjy5p4-profile/include/gsl/gsl_matrix_double.h:286
-#1  CalcXHiY (eval=eval@entry=0x5233c0, D_l=D_l@entry=0x59f980, X=X@entry=0x540930,
-    UltVehiY=UltVehiY@entry=0x53e500, xHiy=xHiy@entry=0x5aaf40) at src/mvlmm.cpp:349
-#2  0x0000000000481033 in MphEM (func_name=func_name@entry=82 'R', max_iter=1000, max_prec=0.001,
-    eval=eval@entry=0x5233c0, X=X@entry=0x540930, Y=Y@entry=0x5408d0, U_hat=U_hat@entry=0x53e380,
-    E_hat=0x53e3e0, OmegaU=0x53e440, OmegaE=0x53e4a0, UltVehiY=0x53e500, UltVehiBX=0x53e560,
-    UltVehiU=0x53e5c0, UltVehiE=0x53e620, V_g=0x540990, V_e=0x540a20, B=0x540ab0) at src/mvlmm.cpp:658
-#3  0x000000000048f9da in MVLMM::AnalyzePlink (this=this@entry=0x7fffffffcf20, U=U@entry=0x523380,
-    eval=eval@entry=0x5233c0, UtW=UtW@entry=0x5233f0, UtY=UtY@entry=0x523430) at src/mvlmm.cpp:3805
-#4  0x0000000000435c2e in GEMMA::BatchRun (this=this@entry=0x7fffffffe4a0, cPar=...)
-    at src/gemma.cpp:2830
-
-#0  0x00007ffff62ffd0b in dgemm_nn ()
-   from /gnu/store/pmmbmrbizz668plrrvfyvqh5ymvjy5p4-profile/lib/libopenblas.so.0
-#1  0x00007ffff622e6f3 in cblas_dgemm ()
-   from /gnu/store/pmmbmrbizz668plrrvfyvqh5ymvjy5p4-profile/lib/libopenblas.so.0
-#2  0x00007ffff7d6af44 in gsl_blas_dgemm ()
-   from /gnu/store/pmmbmrbizz668plrrvfyvqh5ymvjy5p4-profile/lib/libgsl.so.25
-#3  0x000000000048067b in CalcSigma (func_name=func_name@entry=82 'R', eval=eval@entry=0x5233c0,
-    D_l=D_l@entry=0x5a13f0, X=X@entry=0x540930, OmegaU=OmegaU@entry=0x53e440,
-    OmegaE=OmegaE@entry=0x53e4a0, UltVeh=0x5a4410, Qi=0x5a8320, Sigma_uu=0x59fec0, Sigma_ee=0x5acca0)
-    at src/mvlmm.cpp:541
-#4  0x0000000000480f62 in MphEM (func_name=func_name@entry=82 'R', max_iter=1000, max_prec=0.001,
-    eval=eval@entry=0x5233c0, X=X@entry=0x540930, Y=Y@entry=0x5408d0, U_hat=U_hat@entry=0x53e380,
-    E_hat=0x53e3e0, OmegaU=0x53e440, OmegaE=0x53e4a0, UltVehiY=0x53e500, UltVehiBX=0x53e560,
-    UltVehiU=0x53e5c0, UltVehiE=0x53e620, V_g=0x540990, V_e=0x540a20, B=0x540ab0) at src/mvlmm.cpp:704
-#5  0x000000000048fa1d in MVLMM::AnalyzePlink (this=this@entry=0x7fffffffcf20, U=U@entry=0x523380,
-    eval=eval@entry=0x5233c0, UtW=UtW@entry=0x5233f0, UtY=UtY@entry=0x523430) at src/mvlmm.cpp:3806
-```
-
-A profiler should help pinpoint where time is spent. The issue submitter admits that this is a contrived edge-case so I am not going to work on it until I start optimizing mvlmm.