aboutsummaryrefslogtreecommitdiff
path: root/wqflask
diff options
context:
space:
mode:
Diffstat (limited to 'wqflask')
-rw-r--r--wqflask/wqflask/my_pylmm/pyLMM/kinship.py20
-rw-r--r--wqflask/wqflask/my_pylmm/pyLMM/lmm.py13
2 files changed, 27 insertions, 6 deletions
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/kinship.py b/wqflask/wqflask/my_pylmm/pyLMM/kinship.py
index 28f2042d..bb4f5ace 100644
--- a/wqflask/wqflask/my_pylmm/pyLMM/kinship.py
+++ b/wqflask/wqflask/my_pylmm/pyLMM/kinship.py
@@ -21,10 +21,12 @@
import sys
import os
import numpy as np
+from scipy import linalg
import multiprocessing as mp # Multiprocessing is part of the Python stdlib
import Queue
import time
+
from optmatrix import matrix_initialize, matrixMultT
def kinship_full(G):
@@ -153,5 +155,23 @@ def kinship(G,computeSize=1000,numThreads=None,useBLAS=False,verbose=True):
# np.savetxt(outFile+".kve",Kve)
return K
+def kvakve(K, verbose=True):
+ """
+ Obtain eigendecomposition for K and return Kva,Kve where Kva is cleaned
+ of small values < 1e-6 (notably smaller than zero)
+ """
+ if verbose: sys.stderr.write("Obtaining eigendecomposition for %dx%d matrix\n" % (K.shape[0],K.shape[1]) )
+
+ Kva,Kve = linalg.eigh(K)
+ if verbose:
+ print("self.Kva is: ", Kva.shape, Kva)
+ print("self.Kve is: ", Kve.shape, Kve)
+
+ if sum(Kva < 1e-6):
+ if verbose: sys.stderr.write("Cleaning %d eigen values\n" % (sum(Kva < 1e-6)))
+ Kva[Kva < 1e-6] = 1e-6
+ return Kva,Kve
+
+
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/lmm.py b/wqflask/wqflask/my_pylmm/pyLMM/lmm.py
index 1fa7a895..17fe952a 100644
--- a/wqflask/wqflask/my_pylmm/pyLMM/lmm.py
+++ b/wqflask/wqflask/my_pylmm/pyLMM/lmm.py
@@ -49,7 +49,7 @@ has_gn2=True
from utility.benchmark import Bench
from utility import temp_data
-from kinship import kinship, kinship_full
+from kinship import kinship, kinship_full, kvakve
import genotype
try:
@@ -532,11 +532,12 @@ class LMM:
print("this K is:", K.shape, pf(K))
if len(Kva) == 0 or len(Kve) == 0:
- if self.verbose: sys.stderr.write("Obtaining eigendecomposition for %dx%d matrix\n" % (K.shape[0],K.shape[1]) )
- begin = time.time()
- Kva,Kve = linalg.eigh(K)
- end = time.time()
- if self.verbose: sys.stderr.write("Total time: %0.3f\n" % (end - begin))
+ # if self.verbose: sys.stderr.write("Obtaining eigendecomposition for %dx%d matrix\n" % (K.shape[0],K.shape[1]) )
+ # begin = time.time()
+ # Kva,Kve = linalg.eigh(K)
+ # end = time.time()
+ # if self.verbose: sys.stderr.write("Total time: %0.3f\n" % (end - begin))
+ Kva,Kve = kvakve(K)
self.K = K
self.Kva = Kva