From 4931530a6c32eaf56d1586154bc7ea872c0de2ba Mon Sep 17 00:00:00 2001
From: Pjotr Prins
Date: Sat, 14 Mar 2015 15:51:48 +0300
Subject: Tests and indentation

---
 wqflask/wqflask/my_pylmm/pyLMM/kinship.py | 159 +++++++++++++++---------------
 wqflask/wqflask/my_pylmm/pyLMM/runlmm.py  |  13 ++-
 2 files changed, 89 insertions(+), 83 deletions(-)

diff --git a/wqflask/wqflask/my_pylmm/pyLMM/kinship.py b/wqflask/wqflask/my_pylmm/pyLMM/kinship.py
index 1dca1c57..8caa753b 100644
--- a/wqflask/wqflask/my_pylmm/pyLMM/kinship.py
+++ b/wqflask/wqflask/my_pylmm/pyLMM/kinship.py
@@ -57,91 +57,94 @@ def compute_matrixMult(job,W,q = None):
    return job
 
 def f_init(q):
-    compute_matrixMult.q = q
+   compute_matrixMult.q = q
 
 def kinship_full(G):
-    print G.shape
-    m = G.shape[0] # snps
-    n = G.shape[1] # inds
-    sys.stderr.write(str(m)+" SNPs\n")
-    assert m>n, "n should be larger than m (snps>inds)"
-    m = np.dot(G.T,G)
-    m = m/G.shape[0]
-    return m
+   """
+   Calculate the Kinship matrix using a full dot multiplication
+   """
+   print G.shape
+   m = G.shape[0] # snps
+   n = G.shape[1] # inds
+   sys.stderr.write(str(m)+" SNPs\n")
+   assert m>n, "n should be larger than m (snps>inds)"
+   m = np.dot(G.T,G)
+   m = m/G.shape[0]
+   return m
 
 # Calculate the kinship matrix from G (SNPs as rows!), returns K
 #
 def kinship(G,options):
-    numThreads = None
-    if options.numThreads:
-        numThreads = int(options.numThreads)
-    options.computeSize = 1000
-    matrix_initialize(options.useBLAS)
-    
-    sys.stderr.write(str(G.shape)+"\n")
-    n = G.shape[1] # inds
-    inds = n
-    m = G.shape[0] # snps
-    snps = m
-    sys.stderr.write(str(m)+" SNPs\n")
-    assert m>n, "n should be larger than m (snps>inds)"
-
-    q = mp.Queue()
-    p = mp.Pool(numThreads, f_init, [q])
-    iterations = snps/options.computeSize+1
-    if options.testing:
+   numThreads = None
+   if options.numThreads:
+      numThreads = int(options.numThreads)
+   options.computeSize = 1000
+   matrix_initialize(options.useBLAS)
+   
+   sys.stderr.write(str(G.shape)+"\n")
+   n = G.shape[1] # inds
+   inds = n
+   m = G.shape[0] # snps
+   snps = m
+   sys.stderr.write(str(m)+" SNPs\n")
+   assert m>n, "n should be larger than m (snps>inds)"
+
+   q = mp.Queue()
+   p = mp.Pool(numThreads, f_init, [q])
+   iterations = snps/options.computeSize+1
+   if options.testing:
       iterations = 8
-    # jobs = range(0,8) # range(0,iterations)
-
-    results = []
-
-    K = np.zeros((n,n))  # The Kinship matrix has dimension individuals x individuals
-
-    completed = 0
-    for job in range(iterations):
-       if options.verbose:
-          sys.stderr.write("Processing job %d first %d SNPs\n" % (job, ((job+1)*options.computeSize)))
-       W = compute_W(job,G,n,snps,options.computeSize)
-       if numThreads == 1:
-          compute_matrixMult(job,W,q)
-          j,x = q.get()
-          if options.verbose: sys.stderr.write("Job "+str(j)+" finished\n")
-          K_j = x
-          # print j,K_j[:,0]
-          K = K + K_j
-       else:
-          results.append(p.apply_async(compute_matrixMult, (job,W)))
-          # Do we have a result?
-          try: 
-             j,x = q.get_nowait()
-             if options.verbose: sys.stderr.write("Job "+str(j)+" finished\n")
-             K_j = x
-             # print j,K_j[:,0]
-             K = K + K_j
-             completed += 1
-          except Queue.Empty:
-             pass
-
-    if numThreads == None or numThreads > 1:
-       for job in range(len(results)-completed):
-          j,x = q.get(True,15)
-          if options.verbose: sys.stderr.write("Job "+str(j)+" finished\n")
-          K_j = x
-          # print j,K_j[:,0]
-          K = K + K_j
-
-    K = K / float(snps)
-    outFile = 'runtest.kin'
-    if options.verbose: sys.stderr.write("Saving Kinship file to %s\n" % outFile)
-    np.savetxt(outFile,K)
-
-    if options.saveKvaKve:
-       if options.verbose: sys.stderr.write("Obtaining Eigendecomposition\n")
-       Kva,Kve = linalg.eigh(K)
-       if options.verbose: sys.stderr.write("Saving eigendecomposition to %s.[kva | kve]\n" % outFile)
-       np.savetxt(outFile+".kva",Kva)
-       np.savetxt(outFile+".kve",Kve)
-    return K      
+   # jobs = range(0,8) # range(0,iterations)
+
+   results = []
+
+   K = np.zeros((n,n))  # The Kinship matrix has dimension individuals x individuals
+
+   completed = 0
+   for job in range(iterations):
+      if options.verbose:
+         sys.stderr.write("Processing job %d first %d SNPs\n" % (job, ((job+1)*options.computeSize)))
+      W = compute_W(job,G,n,snps,options.computeSize)
+      if numThreads == 1:
+         compute_matrixMult(job,W,q)
+         j,x = q.get()
+         if options.verbose: sys.stderr.write("Job "+str(j)+" finished\n")
+         K_j = x
+         # print j,K_j[:,0]
+         K = K + K_j
+      else:
+         results.append(p.apply_async(compute_matrixMult, (job,W)))
+         # Do we have a result?
+         try: 
+            j,x = q.get_nowait()
+            if options.verbose: sys.stderr.write("Job "+str(j)+" finished\n")
+            K_j = x
+            # print j,K_j[:,0]
+            K = K + K_j
+            completed += 1
+         except Queue.Empty:
+            pass
+
+   if numThreads == None or numThreads > 1:
+      for job in range(len(results)-completed):
+         j,x = q.get(True,15)
+         if options.verbose: sys.stderr.write("Job "+str(j)+" finished\n")
+         K_j = x
+         # print j,K_j[:,0]
+         K = K + K_j
+
+   K = K / float(snps)
+   outFile = 'runtest.kin'
+   if options.verbose: sys.stderr.write("Saving Kinship file to %s\n" % outFile)
+   np.savetxt(outFile,K)
+
+   if options.saveKvaKve:
+      if options.verbose: sys.stderr.write("Obtaining Eigendecomposition\n")
+      Kva,Kve = linalg.eigh(K)
+      if options.verbose: sys.stderr.write("Saving eigendecomposition to %s.[kva | kve]\n" % outFile)
+      np.savetxt(outFile+".kva",Kva)
+      np.savetxt(outFile+".kve",Kve)
+   return K      
 
 
 
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py b/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py
index 26485136..547ac7f1 100644
--- a/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py
+++ b/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py
@@ -140,7 +140,7 @@ elif cmd == 'kinship':
     if not options.skip_genotype_normalization:
         G = np.apply_along_axis( genotype.normalize, axis=1, arr=G)
 
-    if True:
+    if False:
         K = kinship_full(G)
         print "Genotype",G.shape, "\n", G
         print "first Kinship method",K.shape,"\n",K
@@ -157,11 +157,14 @@ elif cmd == 'kinship':
     k3 = round(K3[0][0],4)
     if options.geno == 'data/small.geno':
         assert k1==0.7939, "k1=%f" % k1
-        assert k2==0.7939, "k1=%f" % k1
-        assert k3==0.7939, "k1=%f" % k1
+        assert k2==0.7939, "k2=%f" % k2
+        assert k3==0.7939, "k3=%f" % k3
     if options.geno == 'data/small_na.geno':
         assert k1==0.7172, "k1=%f" % k1
-        assert k2==0.7172, "k1=%f" % k1
-        assert k3==0.7172, "k1=%f" % k1
+        assert k2==0.7172, "k2=%f" % k2
+        assert k3==0.7172, "k3=%f" % k3
+    if options.geno == 'data/test8000.geno':
+        assert k3==1.4352, "k3=%f" % k3
+
 else:
     print "Doing nothing"
-- 
cgit v1.2.3