pylmm has moved out of the GN2 source tree to https://github.com/genenetwork/pylmm_gn2

18 files changed, 0 insertions, 3227 deletions

diff --git a/wqflask/wqflask/my_pylmm/pyLMM/__init__.py b/wqflask/wqflask/my_pylmm/pyLMM/__init__.py
deleted file mode 100644
index f33c4e74..00000000
--- a/wqflask/wqflask/my_pylmm/pyLMM/__init__.py
+++ /dev/null
@@ -1 +0,0 @@
-PYLMM_VERSION="0.51-gn2"
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/benchmark.py b/wqflask/wqflask/my_pylmm/pyLMM/benchmark.py
deleted file mode 100644
index 6c6c9f88..00000000
--- a/wqflask/wqflask/my_pylmm/pyLMM/benchmark.py
+++ /dev/null
@@ -1,44 +0,0 @@
-from __future__ import print_function, division, absolute_import
-
-import collections
-import inspect
-import time
-
-class Bench(object):
-    entries = collections.OrderedDict()
-    
-    def __init__(self, name=None):
-        self.name = name
-
-    def __enter__(self):
-        if self.name:
-            print("Starting benchmark: %s" % (self.name))
-        else:
-            print("Starting benchmark at: %s [%i]" % (inspect.stack()[1][3], inspect.stack()[1][2]))
-        self.start_time = time.time()
-
-    def __exit__(self, type, value, traceback):
-        if self.name:
-            name = self.name
-        else:
-            name = "That"
-
-        time_taken = time.time() - self.start_time
-        print("  %s took: %f seconds" % (name, (time_taken)))
-        
-        if self.name:
-            Bench.entries[self.name] = Bench.entries.get(self.name, 0) + time_taken
-            
-
-    @classmethod
-    def report(cls):
-        total_time = sum((time_taken for time_taken in cls.entries.itervalues()))
-        print("\nTiming report\n")
-        for name, time_taken in cls.entries.iteritems():
-            percent = int(round((time_taken/total_time) * 100))
-            print("[{}%] {}: {}".format(percent, name, time_taken))
-        print()
-        
-    def reset(cls):
-        """Reset the entries"""
-        cls.entries = collections.OrderedDict()
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/chunks.py b/wqflask/wqflask/my_pylmm/pyLMM/chunks.py
deleted file mode 100644
index 9565fb96..00000000
--- a/wqflask/wqflask/my_pylmm/pyLMM/chunks.py
+++ /dev/null
@@ -1,96 +0,0 @@
-from __future__ import absolute_import, print_function, division
-
-import math
-import time
-
-
-def divide_into_chunks(the_list, number_chunks):
-    """Divides a list into approximately number_chunks smaller lists
-    
-    >>> divide_into_chunks([1, 2, 7, 3, 22, 8, 5, 22, 333], 3)
-    [[1, 2, 7], [3, 22, 8], [5, 22, 333]]
-    >>> divide_into_chunks([1, 2, 7, 3, 22, 8, 5, 22, 333], 4)                                                                                                                                                     
-    [[1, 2, 7], [3, 22, 8], [5, 22, 333]]
-    >>> divide_into_chunks([1, 2, 7, 3, 22, 8, 5, 22, 333], 5)                                                                                                                                                     
-    [[1, 2], [7, 3], [22, 8], [5, 22], [333]]
-    >>>
-    
-    """
-    length = len(the_list)
-
-    if length == 0:
-        return [[]]
-    
-    if length <= number_chunks:
-        number_chunks = length
-
-    chunksize = int(math.ceil(length / number_chunks))
-
-    chunks = []
-    for counter in range(0, length, chunksize):
-        chunks.append(the_list[counter:counter+chunksize])
-
-    return chunks
-
-def _confirm_chunk(original, result):
-    all_chunked = []
-    for chunk in result:
-        all_chunked.extend(chunk)
-    print("length of all chunked:", len(all_chunked))
-    assert original == all_chunked, "You didn't chunk right"
-
-
-def _chunk_test(divide_func):
-    import random
-    random.seed(7)
-
-    number_exact = 0
-    total_amount_off = 0
-
-    for test in range(1, 1001):
-        print("\n\ntest:", test)
-        number_chunks = random.randint(1, 20)
-        number_elements = random.randint(0, 100)
-        the_list = list(range(1, number_elements))
-        result = divide_func(the_list, number_chunks)
-
-        print("Dividing list of length {} into approximately {} chunks - got {} chunks".format(
-            len(the_list), number_chunks, len(result)))
-        print("result:", result)
-
-        _confirm_chunk(the_list, result)
-
-        amount_off = abs(number_chunks - len(result))
-        if amount_off == 0:
-            number_exact += 1
-        else:
-            total_amount_off += amount_off
-
-
-        print("\n{} exact out of {}    [Total amount off: {}]".format(number_exact,
-                                                                      test,
-                                                                      total_amount_off))
-    assert number_exact == 558
-    assert total_amount_off == 1580
-    return number_exact, total_amount_off
-
-
-def _main():
-    info = dict()
-    #funcs = (("sam", sam_divide_into_chunks), ("zach", zach_divide_into_chunks))
-    funcs = (("only one", divide_into_chunks),)
-    for name, func in funcs:
-        start = time.time()
-        number_exact, total_amount_off = _chunk_test(func)
-        took = time.time() - start
-        info[name] = dict(number_exact=number_exact,
-                          total_amount_off=total_amount_off,
-                          took=took)
-
-    print("info is:", info)
-
-if __name__ == '__main__':
-    _main()
-    print("\nConfirming doctests...")
-    import doctest
-    doctest.testmod()
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/convertlmm.py b/wqflask/wqflask/my_pylmm/pyLMM/convertlmm.py
deleted file mode 100644
index 4312fed0..00000000
--- a/wqflask/wqflask/my_pylmm/pyLMM/convertlmm.py
+++ /dev/null
@@ -1,184 +0,0 @@
-# This is a converter for common LMM formats, so as to keep file
-# reader complexity outside the main routines.
-
-# Copyright (C) 2015  Pjotr Prins (pjotr.prins@thebird.nl)
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Affero General Public License as
-# published by the Free Software Foundation, either version 3 of the
-# License, or (at your option) any later version.
-
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU Affero General Public License for more details.
-
-# You should have received a copy of the GNU Affero General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-
-from __future__ import print_function
-from optparse import OptionParser
-import sys
-import os
-import numpy as np
-# from lmm import LMM, run_other
-# import input
-import plink
-
-usage = """
-python convertlmm.py [--plink] [--prefix out_basename] [--kinship kfile] [--pheno pname] [--geno gname]
-
-  Convert files for runlmm.py processing. Writes to stdout by default.
-
-  try --help for more information
-
-Examples:
-
-  python ./my_pylmm/pyLMM/convertlmm.py --plink --pheno data/test_snps.132k.clean.noX.fake.phenos > test.pheno
-
-  python ./my_pylmm/pyLMM/convertlmm.py --plink --pheno data/test_snps.132k.clean.noX.fake.phenos --geno data/test_snps.132k.clean.noX > test.geno
-"""
-
-# if len(args) == 0:
-#     print usage
-#     sys.exit(1)
-
-option_parser = OptionParser(usage=usage)
-option_parser.add_option("--kinship", dest="kinship",
-                  help="Parse a kinship file. This is an nxn plain text file and can be computed with the pylmmKinship program")
-option_parser.add_option("--pheno", dest="pheno",
-                         help="Parse a phenotype file (use with --plink only)")
-option_parser.add_option("--geno", dest="geno",
-                         help="Parse a genotype file (use with --plink only)")
-option_parser.add_option("--plink", dest="plink", action="store_true", default=False,
-                  help="Parse PLINK style")
-# option_parser.add_option("--kinship",action="store_false", dest="kinship", default=True,
-#                   help="Parse a kinship file. This is an nxn plain text file and can be computed with the pylmmKinship program.")
-option_parser.add_option("--prefix", dest="prefix",
-                  help="Output prefix for output file(s)")
-option_parser.add_option("-q", "--quiet",
-                  action="store_false", dest="verbose", default=True,
-                  help="don't print status messages to stdout")
-option_parser.add_option("-v", "--verbose",
-                  action="store_true", dest="verbose", default=False,
-                  help="Print extra info")
-
-(options, args) = option_parser.parse_args()
-
-writer = None
-num_inds = None
-snp_names = []
-ind_names = []
-
-def msg(s):
-    sys.stderr.write("INFO: ")
-    sys.stderr.write(s)
-    sys.stderr.write("\n")
-    
-def wr(s):
-    if writer:
-        writer.write(s)
-    else:
-        sys.stdout.write(s)
-
-def wrln(s):
-    wr(s)
-    wr("\n")
-            
-
-if options.pheno:
-    if not options.plink:
-        raise Exception("Use --plink switch")
-    # Because plink does not track size we need to read the whole thing first
-    msg("Converting pheno "+options.pheno)
-    phenos = []
-    count = 0
-    count_pheno = None
-    for line in open(options.pheno,'r'):
-        count += 1
-        list = line.split()
-        pcount = len(list)-2
-        assert(pcount > 0)
-        if count_pheno == None:
-            count_pheno = pcount
-        assert(count_pheno == pcount)
-        row = [list[0]]+list[2:]
-        phenos.append(row)
-
-    writer = None
-    if options.prefix:
-        writer = open(options.prefix+".pheno","w")
-    wrln("# Phenotype format version 1.0")
-    wrln("# Individuals = "+str(count))
-    wrln("# Phenotypes = "+str(count_pheno))
-    for i in range(count_pheno):
-        wr("\t"+str(i+1))
-        wr("\n")
-    for i in range(count):
-        wr("\t".join(phenos[i]))
-        wr("\n")
-    num_inds = count
-    msg(str(count)+" pheno lines written")
-
-if options.kinship:
-    is_header = True
-    count = 0
-    msg("Converting kinship "+options.kinship)
-    writer = None
-    if options.prefix:
-        writer = open(options.prefix+".kin","w")
-    for line in open(options.kinship,'r'):
-        count += 1
-        if is_header:
-            size = len(line.split())
-            wrln("# Kinship format version 1.0")
-            wrln("# Size="+str(size))
-            for i in range(size):
-                wr("\t"+str(i+1))
-            wr("\n")
-            is_header = False
-        wr(str(count))
-        wr("\t")
-        wr("\t".join(line.split()))
-        wr("\n")
-    num_inds = count
-    msg(str(count)+" kinship lines written")
-    
-if options.geno:
-    msg("Converting geno "+options.geno+'.bed')
-    if not options.plink:
-        raise Exception("Use --plink switch")
-    if not num_inds:
-        raise Exception("Can not figure out the number of individuals, use --pheno or --kinship")    
-    bim_snps = plink.readbim(options.geno+'.bim')
-    num_snps = len(bim_snps)
-    writer = None
-    if options.prefix:
-        writer = open(options.prefix+".geno","w")
-    wrln("# Genotype format version 1.0")
-    wrln("# Individuals = "+str(num_inds))
-    wrln("# SNPs = "+str(num_snps))
-    wrln("# Encoding = HAB")
-    for i in range(num_inds):
-        wr("\t"+str(i+1))
-    wr("\n")
-
-    m = []
-    def out(i,x):
-        # wr(str(i)+"\t")
-        # wr("\t".join(x))
-        # wr("\n")
-        m.append(x)
-        
-    snps = plink.readbed(options.geno+'.bed',num_inds, ('A','H','B','-'), out)
-
-    msg("Write transposed genotype matrix")
-    for g in range(num_snps):
-        wr(bim_snps[g][1]+"\t")
-        for i in range(num_inds):
-            wr(m[g][i])
-        wr("\n")
-            
-    msg(str(count)+" geno lines written (with "+str(snps)+" snps)")
-   
-msg("Converting done")
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/genotype.py b/wqflask/wqflask/my_pylmm/pyLMM/genotype.py
deleted file mode 100644
index 49f32e3a..00000000
--- a/wqflask/wqflask/my_pylmm/pyLMM/genotype.py
+++ /dev/null
@@ -1,51 +0,0 @@
-# Genotype routines
-
-# Copyright (C) 2013  Nicholas A. Furlotte (nick.furlotte@gmail.com)
-# Copyright (C) 2015  Pjotr Prins (pjotr.prins@thebird.nl)
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Affero General Public License as
-# published by the Free Software Foundation, either version 3 of the
-# License, or (at your option) any later version.
-
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU Affero General Public License for more details.
-
-# You should have received a copy of the GNU Affero General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-
-import numpy as np
-from collections import Counter
-import operator
-
-def replace_missing_with_MAF(snp_g):
-    """
-    Replace the missing genotype with the minor allele frequency (MAF)
-    in the snp row. It is rather slow!
-    """
-    cnt = Counter(snp_g)
-    tuples = sorted(cnt.items(), key=operator.itemgetter(1))
-    l2 = [t for t in tuples if not np.isnan(t[0])]
-    maf = l2[0][0]
-    res = np.array([maf if np.isnan(snp) else snp for snp in snp_g])
-    return res
-    
-def normalize(ind_g):
-    """
-    Run for every SNP list (for one individual) and return
-    normalized SNP genotype values with missing data filled in
-    """
-    g = np.copy(ind_g)              # copy to avoid side effects
-    missing = np.isnan(g)
-    values = g[True - missing]
-    mean = values.mean()            # Global mean value
-    stddev = np.sqrt(values.var())  # Global stddev
-    g[missing] = mean               # Plug-in mean values for missing data
-    if stddev == 0:
-        g = g - mean                # Subtract the mean
-    else:
-        g = (g - mean) / stddev     # Normalize the deviation
-    return g
-
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/gn2.py b/wqflask/wqflask/my_pylmm/pyLMM/gn2.py
deleted file mode 100644
index 821195c8..00000000
--- a/wqflask/wqflask/my_pylmm/pyLMM/gn2.py
+++ /dev/null
@@ -1,110 +0,0 @@
-# Standalone specific methods and callback handler
-#
-# Copyright (C) 2015  Pjotr Prins (pjotr.prins@thebird.nl)
-#
-# Set the log level with
-#
-#   logging.basicConfig(level=logging.DEBUG)
-
-from __future__ import absolute_import, print_function, division
-
-import numpy as np
-import sys
-import logging
-
-# logger = logging.getLogger(__name__)
-logger = logging.getLogger('lmm2')
-logging.basicConfig(level=logging.DEBUG)
-np.set_printoptions(precision=3,suppress=True)
-
-progress_location = None
-progress_current  = None
-progress_prev_perc     = None
-
-def progress_default_func(location,count,total):
-    global progress_current
-    value = round(count*100.0/total)
-    progress_current = value
-
-progress_func = progress_default_func
-
-def progress_set_func(func):
-    global progress_func
-    progress_func = func
-
-def progress(location, count, total):
-    global progress_location
-    global progress_prev_perc
-
-    perc = round(count*100.0/total)
-    if perc != progress_prev_perc and (location != progress_location or perc > 98 or perc > progress_prev_perc + 5):
-        progress_func(location, count, total)
-        logger.info("Progress: %s %d%%" % (location,perc))
-        progress_location = location
-        progress_prev_perc = perc
-
-def mprint(msg,data):
-    """
-    Array/matrix print function
-    """
-    m = np.array(data)
-    if m.ndim == 1:
-        print(msg,m.shape,"=\n",m[0:3]," ... ",m[-3:])
-    if m.ndim == 2:
-        print(msg,m.shape,"=\n[",
-              m[0][0:3]," ... ",m[0][-3:],"\n ",
-              m[1][0:3]," ... ",m[1][-3:],"\n  ...\n ",
-              m[-2][0:3]," ... ",m[-2][-3:],"\n ",
-              m[-1][0:3]," ... ",m[-1][-3:],"]")
-
-def fatal(msg):
-    logger.critical(msg)
-    raise Exception(msg)
-
-def callbacks():
-    return dict(
-        write = sys.stdout.write,
-        writeln = print,
-        debug = logger.debug,
-        info = logger.info,
-        warning = logger.warning,
-        error = logger.error,
-        critical = logger.critical,
-        fatal = fatal,
-        progress = progress,
-        mprint = mprint
-    )
-
-def uses(*funcs):
-    """
-    Some sugar
-    """
-    return [callbacks()[func] for func in funcs]
-    
-# ----- Minor test cases:
-
-if __name__ == '__main__':
-    # logging.basicConfig(level=logging.DEBUG)
-    logging.debug("Test %i" % (1))
-    d = callbacks()['debug']
-    d("TEST")
-    wrln = callbacks()['writeln']
-    wrln("Hello %i" % 34)
-    progress = callbacks()['progress']
-    progress("I am half way",50,100)
-    list = [0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15,
-            0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15,
-            0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15,
-            0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15,
-            0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15]
-    mprint("list",list)
-    matrix = [[1,0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15],
-            [2,0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15],
-            [3,0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15],
-            [4,0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15],
-            [5,0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15],
-            [6,0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15]]
-    mprint("matrix",matrix)
-    ix,dx = uses("info","debug")
-    ix("ix")
-    dx("dx")
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/gwas.py b/wqflask/wqflask/my_pylmm/pyLMM/gwas.py
deleted file mode 100644
index 247a8729..00000000
--- a/wqflask/wqflask/my_pylmm/pyLMM/gwas.py
+++ /dev/null
@@ -1,165 +0,0 @@
-# pylmm-based GWAS calculation
-#
-# Copyright (C) 2013  Nicholas A. Furlotte (nick.furlotte@gmail.com)
-# Copyright (C) 2015  Pjotr Prins (pjotr.prins@thebird.nl)
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Affero General Public License as
-# published by the Free Software Foundation, either version 3 of the
-# License, or (at your option) any later version.
-
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU Affero General Public License for more details.
-
-# You should have received a copy of the GNU Affero General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-#!/usr/bin/python
-
-import pdb
-import time
-import sys
-import lmm2
-
-import os
-import numpy as np
-import input
-from optmatrix import matrix_initialize
-from lmm2 import LMM2
-
-import multiprocessing as mp # Multiprocessing is part of the Python stdlib
-import Queue 
-
-# ---- A trick to decide on the environment:
-try:
-    from wqflask.my_pylmm.pyLMM import chunks
-    from gn2 import uses
-except ImportError:
-    has_gn2=False
-    from standalone import uses
-
-progress,mprint,debug,info,fatal = uses('progress','mprint','debug','info','fatal')
-
-
-def formatResult(id,beta,betaSD,ts,ps):
-   return "\t".join([str(x) for x in [id,beta,betaSD,ts,ps]]) + "\n"
-
-def compute_snp(j,n,snp_ids,lmm2,REML,q = None):
-   result = []
-   for snp_id in snp_ids:
-      snp,id = snp_id
-      x = snp.reshape((n,1))  # all the SNPs
-      # if refit:
-      #    L.fit(X=snp,REML=REML)
-      ts,ps,beta,betaVar = lmm2.association(x,REML=REML,returnBeta=True)
-      # result.append(formatResult(id,beta,np.sqrt(betaVar).sum(),ts,ps))
-      result.append( (ts,ps) )
-   if not q:
-      q = compute_snp.q
-   q.put([j,result])
-   return j
-
-def f_init(q):
-   compute_snp.q = q
-
-def gwas(Y,G,K,restricted_max_likelihood=True,refit=False,verbose=True):
-   """
-   GWAS. The G matrix should be n inds (cols) x m snps (rows)
-   """
-   info("In gwas.gwas")
-   matrix_initialize()
-   cpu_num = mp.cpu_count()
-   numThreads = None # for now use all available threads
-   kfile2 = False
-   reml = restricted_max_likelihood
-
-   mprint("G",G)
-   n = G.shape[1] # inds
-   inds = n
-   m = G.shape[0] # snps
-   snps = m
-   info("%s SNPs",snps)
-   assert snps>=inds, "snps should be larger than inds (snps=%d,inds=%d)" % (snps,inds)
-
-   # CREATE LMM object for association
-   # if not kfile2:  L = LMM(Y,K,Kva,Kve,X0,verbose=verbose)
-   # else:  L = LMM_withK2(Y,K,Kva,Kve,X0,verbose=verbose,K2=K2)
-
-   lmm2 = LMM2(Y,K) # ,Kva,Kve,X0,verbose=verbose)
-   if not refit:
-      info("Computing fit for null model")
-      lmm2.fit()  # follow GN model in run_other
-      info("heritability=%0.3f, sigma=%0.3f" % (lmm2.optH,lmm2.optSigma))
-            
-   res = []
-
-   # Set up the pool
-   # mp.set_start_method('spawn')
-   q = mp.Queue()
-   p = mp.Pool(numThreads, f_init, [q])
-   collect = []
-
-   count = 0
-   job = 0
-   jobs_running = 0
-   jobs_completed = 0
-   for snp in G:
-      snp_id = (snp,'SNPID')
-      count += 1
-      if count % 1000 == 0:
-         job += 1
-         debug("Job %d At SNP %d" % (job,count))
-         if numThreads == 1:
-            debug("Running on 1 THREAD")
-            compute_snp(job,n,collect,lmm2,reml,q)
-            collect = []
-            j,lst = q.get()
-            debug("Job "+str(j)+" finished")
-            jobs_completed += 1
-            progress("GWAS2",jobs_completed,snps/1000)
-            res.append((j,lst))
-         else:
-            p.apply_async(compute_snp,(job,n,collect,lmm2,reml))
-            jobs_running += 1
-            collect = []
-            while jobs_running > cpu_num:
-               try:
-                  j,lst = q.get_nowait()
-                  debug("Job "+str(j)+" finished")
-                  jobs_completed += 1
-                  progress("GWAS2",jobs_completed,snps/1000)
-                  res.append((j,lst))
-                  jobs_running -= 1
-               except Queue.Empty:
-                  time.sleep(0.1)
-                  pass
-               if jobs_running > cpu_num*2:
-                  time.sleep(1.0)
-               else:
-                  break
-
-      collect.append(snp_id)
-
-   if numThreads==1 or count<1000 or len(collect)>0:
-      job += 1
-      debug("Collect final batch size %i job %i @%i: " % (len(collect), job, count))
-      compute_snp(job,n,collect,lmm2,reml,q)
-      collect = []
-      j,lst = q.get()
-      res.append((j,lst))
-   debug("count=%i running=%i collect=%i" % (count,jobs_running,len(collect)))
-   for job in range(jobs_running):
-      j,lst = q.get(True,15) # time out
-      debug("Job "+str(j)+" finished")
-      jobs_completed += 1
-      progress("GWAS2",jobs_completed,snps/1000)
-      res.append((j,lst))
-
-   mprint("Before sort",[res1[0] for res1 in res])
-   res = sorted(res,key=lambda x: x[0])
-   mprint("After sort",[res1[0] for res1 in res])
-   info([len(res1[1]) for res1 in res])
-   ts = [item[0] for j,res1 in res for item in res1]
-   ps = [item[1] for j,res1 in res for item in res1]
-   return ts,ps
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/input.py b/wqflask/wqflask/my_pylmm/pyLMM/input.py
deleted file mode 100644
index 7063fedf..00000000
--- a/wqflask/wqflask/my_pylmm/pyLMM/input.py
+++ /dev/null
@@ -1,267 +0,0 @@
-# pylmm is a python-based linear mixed-model solver with applications to GWAS
-
-# Copyright (C) 2013  Nicholas A. Furlotte (nick.furlotte@gmail.com)
-
-#The program is free for academic use. Please contact Nick Furlotte
-#<nick.furlotte@gmail.com> if you are interested in using the software for
-#commercial purposes.
-
-#The software must not be modified and distributed without prior
-#permission of the author.
-
-#THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
-#"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
-#LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
-#A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
-#CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
-#EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
-#PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
-#PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
-#LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
-#NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
-#SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
-
-import os
-import sys
-import numpy as np
-import struct
-import pdb
-
-class plink:
-    def __init__(self,fbase,kFile=None,phenoFile=None,type='b',normGenotype=True,readKFile=False):
-        self.fbase = fbase
-        self.type = type
-        self.indivs = self.getIndivs(self.fbase,type)
-        self.kFile = kFile
-        self.phenos = None
-        self.normGenotype = normGenotype
-        self.phenoFile = phenoFile
-        # Originally I was using the fastLMM style that has indiv IDs embedded.
-        # NOW I want to use this module to just read SNPs so I'm allowing 
-        # the programmer to turn off the kinship reading.
-        self.readKFile = readKFile
-
-        if self.kFile:
-            self.K = self.readKinship(self.kFile)
-        elif os.path.isfile("%s.kin" % fbase): 
-            self.kFile = "%s.kin" %fbase
-            if self.readKFile:
-                self.K = self.readKinship(self.kFile)
-        else: 
-            self.kFile = None
-            self.K = None
-
-        self.getPhenos(self.phenoFile)
-
-        self.fhandle = None
-        self.snpFileHandle = None
-
-    def __del__(self): 
-        if self.fhandle: self.fhandle.close()
-        if self.snpFileHandle: self.snpFileHandle.close()
-
-    def getSNPIterator(self):
-        if not self.type == 'b': 
-            sys.stderr.write("Have only implemented this for binary plink files (bed)\n")
-            return
-
-        # get the number of snps
-        file = self.fbase + '.bim'
-        i = 0
-        f = open(file,'r')
-        for line in f: i += 1
-        f.close()
-        self.numSNPs = i
-        self.have_read = 0
-        self.snpFileHandle = open(file,'r')
-
-        self.BytestoRead = self.N / 4 + (self.N % 4 and 1 or 0)
-        self._formatStr = 'c'*self.BytestoRead
-
-        file = self.fbase + '.bed'
-        self.fhandle = open(file,'rb')
-     
-        magicNumber = self.fhandle.read(2)
-        order = self.fhandle.read(1)
-        if not order == '\x01': 
-            sys.stderr.write("This is not in SNP major order - you did not handle this case\n")
-            raise StopIteration
-
-        return self
-
-    def __iter__(self):
-        return self.getSNPIterator()
-
-    def next(self):
-        if self.have_read == self.numSNPs:
-            raise StopIteration
-        X = self.fhandle.read(self.BytestoRead)
-        XX = [bin(ord(x)) for x in struct.unpack(self._formatStr,X)]
-        self.have_read += 1
-        return self.formatBinaryGenotypes(XX,self.normGenotype),self.snpFileHandle.readline().strip().split()[1]
-    
-    def formatBinaryGenotypes(self,X,norm=True):
-        D = { \
-              '00': 0.0, \
-              '10': 0.5, \
-              '11': 1.0, \
-              '01': np.nan \
-           }
-  
-        D_tped = { \
-              '00': '1 1', \
-              '10': '1 2', \
-              '11': '2 2', \
-              '01': '0 0' \
-           }
-  
-        #D = D_tped
-              
-        G = []
-        for x in X:
-            if not len(x) == 10:
-                xx = x[2:]
-                x = '0b' + '0'*(8 - len(xx)) + xx
-            a,b,c,d = (x[8:],x[6:8],x[4:6],x[2:4]) 
-            L = [D[y] for y in [a,b,c,d]]
-            G += L
-        # only take the leading values because whatever is left should be null
-        G = G[:self.N]
-        G = np.array(G)
-        if norm:
-            G = self.normalizeGenotype(G)
-        return G
-    
-    def normalizeGenotype(self,G):
-        # print "Before",G
-        # print G.shape
-        print "call input.normalizeGenotype"
-        raise "This should not be used"
-        x = True - np.isnan(G)
-        m = G[x].mean()
-        s = np.sqrt(G[x].var())
-        G[np.isnan(G)] = m
-        if s == 0: G = G - m
-        else: G = (G - m) / s
-        # print "After",G
-        return G
-    
-    def getPhenos(self,phenoFile=None):
-        if not phenoFile:
-            self.phenoFile = phenoFile = self.fbase+".phenos"
-        if not os.path.isfile(phenoFile): 
-            sys.stderr.write("Could not find phenotype file: %s\n" % (phenoFile))
-            return
-        f = open(phenoFile,'r')
-        keys = []
-        P = []
-        for line in f:
-            v = line.strip().split()
-            keys.append((v[0],v[1]))
-            P.append([(x == 'NA' or x == '-9') and np.nan or float(x) for x in v[2:]])
-        f.close()
-        P = np.array(P)
-     
-        # reorder to match self.indivs
-        D = {}
-        L = []
-        for i in range(len(keys)):
-            D[keys[i]] = i
-        for i in range(len(self.indivs)):
-            if not D.has_key(self.indivs[i]):
-                continue 
-            L.append(D[self.indivs[i]])
-        P = P[L,:]
-     
-        self.phenos = P
-        return P
-    
-    def getIndivs(self,base,type='b'):
-        if type == 't':
-            famFile = "%s.tfam" % base
-        else:
-            famFile = "%s.fam" % base
-        keys = []
-        i = 0
-        f = open(famFile,'r')
-        for line in f:
-           v = line.strip().split()
-           famId = v[0]
-           indivId = v[1]
-           k = (famId.strip(),indivId.strip())
-           keys.append(k)
-           i += 1
-        f.close()
-     
-        self.N = len(keys)
-        sys.stderr.write("Read %d individuals from %s\n" % (self.N, famFile))
-     
-        return keys
-    
-    def readKinship(self,kFile):
-        # Assume the fastLMM style
-        # This will read in the kinship matrix and then reorder it
-        # according to self.indivs - additionally throwing out individuals 
-        # that are not in both sets
-        if self.indivs == None or len(self.indivs) == 0:
-            sys.stderr.write("Did not read any individuals so can't load kinship\n")
-            return 
-     
-        sys.stderr.write("Reading kinship matrix from %s\n" % (kFile) )
-     
-        f = open(kFile,'r')
-        # read indivs 
-        v = f.readline().strip().split("\t")[1:]
-        keys = [tuple(y.split()) for y in v]
-        D = {}
-        for i in range(len(keys)): D[keys[i]] = i
-     
-        # read matrix
-        K = []
-        for line in f:
-            K.append([float(x) for x in line.strip().split("\t")[1:]])
-        f.close()
-        K  = np.array(K)
-     
-        # reorder to match self.indivs
-        L = []
-        KK = []
-        X = []
-        for i in range(len(self.indivs)):
-            if not D.has_key(self.indivs[i]):
-                X.append(self.indivs[i])
-            else: 
-                KK.append(self.indivs[i])
-                L.append(D[self.indivs[i]])
-        K = K[L,:][:,L]
-        self.indivs = KK
-        self.indivs_removed = X
-        if len(self.indivs_removed):
-            sys.stderr.write("Removed %d individuals that did not appear in Kinship\n" % (len(self.indivs_removed)))
-        return K 
-     
-    def getCovariates(self,covFile=None):
-        if not os.path.isfile(covFile): 
-            sys.stderr.write("Could not find covariate file: %s\n" % (phenoFile))
-            return
-        f = open(covFile,'r')
-        keys = []
-        P = []
-        for line in f:
-            v = line.strip().split()
-            keys.append((v[0],v[1]))
-            P.append([x == 'NA' and np.nan or float(x) for x in v[2:]])
-        f.close()
-        P = np.array(P)
-     
-        # reorder to match self.indivs
-        D = {}
-        L = []
-        for i in range(len(keys)):
-            D[keys[i]] = i
-        for i in range(len(self.indivs)):
-            if not D.has_key(self.indivs[i]): continue 
-            L.append(D[self.indivs[i]])
-        P = P[L,:]
-     
-        return P
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/kinship.py b/wqflask/wqflask/my_pylmm/pyLMM/kinship.py
deleted file mode 100644
index 1c157fd8..00000000
--- a/wqflask/wqflask/my_pylmm/pyLMM/kinship.py
+++ /dev/null
@@ -1,168 +0,0 @@
-# pylmm kinship calculation
-#
-# Copyright (C) 2013  Nicholas A. Furlotte (nick.furlotte@gmail.com)
-# Copyright (C) 2015  Pjotr Prins (pjotr.prins@thebird.nl)
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Affero General Public License as
-# published by the Free Software Foundation, either version 3 of the
-# License, or (at your option) any later version.
-
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU Affero General Public License for more details.
-
-# You should have received a copy of the GNU Affero General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-
-# env PYTHONPATH=$pylmm_lib_path:./lib python $pylmm_lib_path/runlmm.py --pheno test.pheno --geno test9000.geno kinship --test
-
-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
-
-# ---- A trick to decide on the environment:
-try:
-    from wqflask.my_pylmm.pyLMM import chunks
-    from gn2 import uses, progress_set_func
-except ImportError:
-    has_gn2=False
-    import standalone as handlers
-    from standalone import uses, progress_set_func
-
-progress,debug,info,mprint = uses('progress','debug','info','mprint')
-
-def kinship_full(G):
-   """
-   Calculate the Kinship matrix using a full dot multiplication
-   """
-   # mprint("kinship_full G",G)
-   m = G.shape[0] # snps
-   n = G.shape[1] # inds
-   info("%d SNPs",m)
-   assert m>n, "n should be larger than m (%d snps > %d inds)" % (m,n)
-   # m = np.dot(G.T,G)
-   m = matrixMultT(G.T)
-   m = m/G.shape[0]
-   # mprint("kinship_full K",m)
-   return m
-
-def compute_W(job,G,n,snps,compute_size):
-   """
-   Read 1000 SNPs at a time into matrix and return the result
-   """
-   m = compute_size
-   W = np.ones((n,m)) * np.nan # W matrix has dimensions individuals x SNPs (initially all NaNs)
-   for j in range(0,compute_size):
-      pos = job*m + j # real position
-      if pos >= snps:
-         W = W[:,range(0,j)]
-         break
-      snp = G[job*compute_size+j]
-      if snp.var() == 0:
-         continue
-      W[:,j] = snp  # set row to list of SNPs
-   return W
-
-def compute_matrixMult(job,W,q = None):
-   """
-   Compute Kinship(W)*j
-
-   For every set of SNPs matrixMult is used to multiply matrices T(W)*W
-   """
-   res = matrixMultT(W)
-   if not q: q=compute_matrixMult.q
-   q.put([job,res])
-   return job
-
-def f_init(q):
-   compute_matrixMult.q = q
-
-# Calculate the kinship matrix from G (SNPs as rows!), returns K
-#
-def kinship(G,computeSize=1000,numThreads=None,useBLAS=False):
-
-   matrix_initialize(useBLAS)
-
-   mprint("G",G)
-   n = G.shape[1] # inds
-   inds = n
-   m = G.shape[0] # snps
-   snps = m
-   info("%i SNPs" % (m))
-   assert snps>=inds, "snps should be larger than inds (%i snps, %i inds)" % (snps,inds)
-
-   q = mp.Queue()
-   p = mp.Pool(numThreads, f_init, [q])
-   cpu_num = mp.cpu_count()
-   info("CPU cores: %i" % cpu_num)
-   iterations = snps/computeSize+1
-
-   results = []
-   K = np.zeros((n,n))  # The Kinship matrix has dimension individuals x individuals
-
-   completed = 0
-   for job in range(iterations):
-      info("Processing job %d first %d SNPs" % (job, ((job+1)*computeSize)))
-      W = compute_W(job,G,n,snps,computeSize)
-      if numThreads == 1:
-         # Single-core
-         compute_matrixMult(job,W,q)
-         j,x = q.get()
-         debug("Job "+str(j)+" finished")
-         progress("kinship",j,iterations)
-         K_j = x
-         K = K + K_j
-      else:
-         # Multi-core
-         results.append(p.apply_async(compute_matrixMult, (job,W)))
-         # Do we have a result?
-         while (len(results)-completed>cpu_num*2):
-            time.sleep(0.1)
-            try: 
-               j,x = q.get_nowait()
-               debug("Job "+str(j)+" finished")
-               K_j = x
-               K = K + K_j
-               completed += 1
-               progress("kinship",completed,iterations)
-            except Queue.Empty:
-               pass
-         
-   if numThreads == None or numThreads > 1:
-      for job in range(len(results)-completed):
-         j,x = q.get(True,15)
-         debug("Job "+str(j)+" finished")
-         K_j = x
-         K = K + K_j
-         completed += 1
-         progress("kinship",completed,iterations)
-
-   K = K / float(snps)
-   return K      
-
-def kvakve(K):
-   """
-   Obtain eigendecomposition for K and return Kva,Kve where Kva is cleaned
-   of small values < 1e-6 (notably smaller than zero)
-   """
-   info("Obtaining eigendecomposition for %dx%d matrix" % (K.shape[0],K.shape[1]) )
-   Kva,Kve = linalg.eigh(K)
-   mprint("Kva",Kva)
-   mprint("Kve",Kve)
-
-   if sum(Kva < 0):
-      info("Cleaning %d eigen values (Kva<0)" % (sum(Kva < 0)))
-      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
deleted file mode 100644
index 2a0c7fdc..00000000
--- a/wqflask/wqflask/my_pylmm/pyLMM/lmm.py
+++ /dev/null
@@ -1,995 +0,0 @@
-# pylmm is a python-based linear mixed-model solver with applications to GWAS
-
-# Copyright (C) 2013  Nicholas A. Furlotte (nick.furlotte@gmail.com)
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Affero General Public License as
-# published by the Free Software Foundation, either version 3 of the
-# License, or (at your option) any later version.
-
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU Affero General Public License for more details.
-
-# You should have received a copy of the GNU Affero General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-
-from __future__ import absolute_import, print_function, division
-
-import sys
-import time
-import uuid
-
-import numpy as np
-from scipy import linalg
-from scipy import optimize
-from scipy import stats
-# import pdb
-
-# import gzip
-# import zlib
-import datetime
-# import cPickle as pickle
-from pprint import pformat as pf
-
-# Add local dir to PYTHONPATH
-import os
-cwd = os.path.dirname(__file__)
-if sys.path[0] != cwd:
-    sys.path.insert(1,cwd)
-
-# pylmm imports
-from kinship import kinship, kinship_full, kvakve
-import genotype
-import phenotype
-import gwas
-from benchmark import Bench
-
-# The following imports are for exchanging data with the webserver
-import simplejson as json
-from redis import Redis
-Redis = Redis()
-import temp_data
-
-has_gn2=None
-
-# sys.stderr.write("INFO: pylmm system path is "+":".join(sys.path)+"\n")
-sys.stderr.write("INFO: pylmm file is "+__file__+"\n")
-
-# ---- A trick to decide on the environment:
-try:
-    sys.stderr.write("INFO: lmm try loading module\n")
-    import utility.formatting # this is never used, just to check the environment
-    sys.stderr.write("INFO: This is a genenetwork2 environment\n")
-    from gn2 import uses, progress_set_func
-    has_gn2=True
-except ImportError:
-    # Failed to load gn2
-    has_gn2=False
-    import standalone as handlers
-    from standalone import uses, progress_set_func
-    sys.stderr.write("WARNING: LMM standalone version missing the Genenetwork2 environment\n")
-
-progress,mprint,debug,info,fatal = uses('progress','mprint','debug','info','fatal')
-
-#np.seterr('raise')
-
-#def run_human(pheno_vector,
-#            covariate_matrix,
-#            plink_input_file,
-#            kinship_matrix,
-#            refit=False,
-#            loading_progress=None):
-
-def run_human(pheno_vector,
-            covariate_matrix,
-            plink_input_file,
-            kinship_matrix,
-            refit=False):
-
-    v = np.isnan(pheno_vector)
-    keep = True - v
-    keep = keep.reshape((len(keep),))
-
-    identifier = str(uuid.uuid4())
-    
-    #print("pheno_vector: ", pf(pheno_vector))
-    #print("kinship_matrix: ", pf(kinship_matrix))
-    #print("kinship_matrix.shape: ", pf(kinship_matrix.shape))
-
-    #lmm_vars = pickle.dumps(dict(
-    #    pheno_vector = pheno_vector,
-    #    covariate_matrix = covariate_matrix,
-    #    kinship_matrix = kinship_matrix
-    #))
-    #Redis.hset(identifier, "lmm_vars", lmm_vars)
-    #Redis.expire(identifier, 60*60)
-
-    if v.sum():
-        pheno_vector = pheno_vector[keep]
-        print("pheno_vector shape is now: ", pf(pheno_vector.shape))
-        covariate_matrix = covariate_matrix[keep,:]
-        print("kinship_matrix shape is: ", pf(kinship_matrix.shape))
-        print("keep is: ", pf(keep.shape))
-        kinship_matrix = kinship_matrix[keep,:][:,keep]
-
-    print("kinship_matrix:", pf(kinship_matrix))
-
-    n = kinship_matrix.shape[0]
-    print("n is:", n)
-    lmm_ob = LMM(pheno_vector,
-                kinship_matrix,
-                covariate_matrix)
-    lmm_ob.fit()
-
-
-    # Buffers for pvalues and t-stats
-    p_values = []
-    t_stats = []
-
-    #print("input_file: ", plink_input_file)
-
-    with Bench("Opening and loading pickle file"):
-        with gzip.open(plink_input_file, "rb") as input_file:
-            data = pickle.load(input_file)
-            
-    plink_input = data['plink_input']
-
-    #plink_input.getSNPIterator()
-    with Bench("Calculating numSNPs"):
-        total_snps = data['numSNPs']
-
-    with Bench("snp iterator loop"):
-        count = 0
-
-        with Bench("Create list of inputs"):
-            inputs = list(plink_input)
-
-        with Bench("Divide into chunks"):
-            results = chunks.divide_into_chunks(inputs, 64)
-
-        result_store = []
-
-        key = "plink_inputs"
-        
-        # Todo: Delete below line when done testing
-        Redis.delete(key)
-        
-        timestamp = datetime.datetime.utcnow().isoformat()
-
-        # Pickle chunks of input SNPs (from Plink interator) and compress them
-        #print("Starting adding loop")
-        for part, result in enumerate(results):
-            #data = pickle.dumps(result, pickle.HIGHEST_PROTOCOL)
-            holder = pickle.dumps(dict(
-                identifier = identifier,
-                part = part,
-                timestamp = timestamp,
-                result = result
-            ), pickle.HIGHEST_PROTOCOL)
-            
-            #print("Adding:", part)
-            Redis.rpush(key, zlib.compress(holder))
-        #print("End adding loop")
-        #print("***** Added to {} queue *****".format(key))
-        for snp, this_id in plink_input:
-            #with Bench("part before association"):
-            #if count > 1000:
-            #    break
-            count += 1
-            progress("human",count,total_snps)
-
-            #with Bench("actual association"):
-            ps, ts = human_association(snp,
-                                       n,
-                                       keep,
-                                       lmm_ob,
-                                       pheno_vector,
-                                       covariate_matrix,
-                                       kinship_matrix,
-                                       refit)
-
-            #with Bench("after association"):
-            p_values.append(ps)
-            t_stats.append(ts)
-        
-    return p_values, t_stats
-
-
-#class HumanAssociation(object):
-#    def __init__(self):
-#        
-
-def human_association(snp,
-                      n,
-                      keep,
-                      lmm_ob,
-                      pheno_vector,
-                      covariate_matrix,
-                      kinship_matrix,
-                      refit):
-
-    x = snp[keep].reshape((n,1))
-    #x[[1,50,100,200,3000],:] = np.nan
-    v = np.isnan(x).reshape((-1,))
-
-    # Check SNPs for missing values
-    if v.sum():
-        keeps = True - v
-        xs = x[keeps,:]
-        # If no variation at this snp or all genotypes missing 
-        if keeps.sum() <= 1 or xs.var() <= 1e-6:
-            return np.nan, np.nan
-            #p_values.append(np.nan)
-            #t_stats.append(np.nan)
-            #continue
-
-        # Its ok to center the genotype -  I used options.normalizeGenotype to
-        # force the removal of missing genotypes as opposed to replacing them with MAF.
-
-        #if not options.normalizeGenotype:
-        #    xs = (xs - xs.mean()) / np.sqrt(xs.var())
-
-        filtered_pheno = pheno_vector[keeps]
-        filtered_covariate_matrix = covariate_matrix[keeps,:]
-        
-        print("kinship_matrix shape is: ", pf(kinship_matrix.shape))
-        print("keeps is: ", pf(keeps.shape))
-        filtered_kinship_matrix = kinship_matrix[keeps,:][:,keeps]
-        filtered_lmm_ob = lmm.LMM(filtered_pheno,filtered_kinship_matrix,X0=filtered_covariate_matrix)
-        if refit:
-            filtered_lmm_ob.fit(X=xs)
-        else:
-            #try:
-            filtered_lmm_ob.fit()
-            #except: pdb.set_trace()
-        ts,ps,beta,betaVar = Ls.association(xs,returnBeta=True)
-    else:
-        if x.var() == 0:
-            return np.nan, np.nan
-            #p_values.append(np.nan)
-            #t_stats.append(np.nan)
-            #continue
-        if refit:
-            lmm_ob.fit(X=x)
-        ts, ps, beta, betaVar = lmm_ob.association(x)
-    return ps, ts
-
-
-#def run(pheno_vector,
-#        genotype_matrix,
-#        restricted_max_likelihood=True,
-#        refit=False,
-#        temp_data=None):
-    
-def run_other_old(pheno_vector,
-        genotype_matrix,
-        restricted_max_likelihood=True,
-        refit=False):
-    
-    """Takes the phenotype vector and genotype matrix and returns a set of p-values and t-statistics
-    
-    restricted_max_likelihood -- whether to use restricted max likelihood; True or False
-    refit -- whether to refit the variance component for each marker
-    
-    """
-    
-    print("Running the original LMM engine in run_other (old)")
-    print("REML=",restricted_max_likelihood," REFIT=",refit)
-    with Bench("Calculate Kinship"):
-        kinship_matrix,genotype_matrix = calculate_kinship_new(genotype_matrix)
-    
-    print("kinship_matrix: ", pf(kinship_matrix))
-    print("kinship_matrix.shape: ", pf(kinship_matrix.shape))
-    
-    # with Bench("Create LMM object"):
-    #     lmm_ob = LMM(pheno_vector, kinship_matrix)
-    
-    # with Bench("LMM_ob fitting"):
-    #     lmm_ob.fit()
-
-    print("run_other_old genotype_matrix: ", genotype_matrix.shape)
-    print(genotype_matrix)
-
-    with Bench("Doing GWAS"):
-        t_stats, p_values = GWAS(pheno_vector,
-                                      genotype_matrix.T,
-                                      kinship_matrix,
-                                      restricted_max_likelihood=True,
-                                      refit=False)
-    Bench().report()
-    return p_values, t_stats
-
-def run_other_new(n,m,pheno_vector,
-        geno,
-        restricted_max_likelihood=True,
-        refit=False):
-    
-    """Takes the phenotype vector and genotype matrix and returns a set of p-values and t-statistics
-    
-    restricted_max_likelihood -- whether to use restricted max likelihood; True or False
-    refit -- whether to refit the variance component for each marker
-    
-    """
-    
-    print("Running the new LMM2 engine in run_other_new")
-    print("REML=",restricted_max_likelihood," REFIT=",refit)
-
-    # Adjust phenotypes
-    n,Y,keep = phenotype.remove_missing_new(n,pheno_vector)
-
-    # if options.maf_normalization:
-    #     G = np.apply_along_axis( genotype.replace_missing_with_MAF, axis=0, arr=g )
-    #     print "MAF replacements: \n",G
-    # if not options.skip_genotype_normalization:
-    # G = np.apply_along_axis( genotype.normalize, axis=1, arr=G)
-
-    geno = geno[:,keep]
-    with Bench("Calculate Kinship"):
-        K,G = calculate_kinship_new(geno)
-    
-    print("kinship_matrix: ", pf(K))
-    print("kinship_matrix.shape: ", pf(K.shape))
-
-    # with Bench("Create LMM object"):
-    #     lmm_ob = lmm2.LMM2(Y,K)
-    # with Bench("LMM_ob fitting"):
-    #     lmm_ob.fit()
-
-    print("run_other_new genotype_matrix: ", G.shape)
-    print(G)
-
-    with Bench("Doing GWAS"):
-        t_stats, p_values = gwas.gwas(Y,
-                                      G,
-                                      K,
-                                      restricted_max_likelihood=True,
-                                      refit=False,verbose=True)
-    Bench().report()
-    return p_values, t_stats
-
-# def matrixMult(A,B):
-#     return np.dot(A,B)
-
-def matrixMult(A,B):
-
-    # If there is no fblas then we will revert to np.dot()
-
-    try:
-        linalg.fblas
-    except AttributeError:
-        return np.dot(A,B)
-
-    #print("A is:", pf(A.shape))
-    #print("B is:", pf(B.shape))
-
-    # If the matrices are in Fortran order then the computations will be faster
-    # when using dgemm.  Otherwise, the function will copy the matrix and that takes time.
-    if not A.flags['F_CONTIGUOUS']:
-       AA = A.T
-       transA = True
-    else:
-       AA = A
-       transA = False
-    
-    if not B.flags['F_CONTIGUOUS']:
-       BB = B.T
-       transB = True
-    else:
-       BB = B
-       transB = False
-
-    return linalg.fblas.dgemm(alpha=1.,a=AA,b=BB,trans_a=transA,trans_b=transB)
-
-def calculate_kinship_new(genotype_matrix):
-    """ 
-    Call the new kinship calculation where genotype_matrix contains
-    inds (columns) by snps (rows).
-    """
-    assert type(genotype_matrix) is np.ndarray
-    info("call genotype.normalize")
-    G = np.apply_along_axis( genotype.normalize, axis=1, arr=genotype_matrix)
-    mprint("G",genotype_matrix)
-    info("call calculate_kinship_new")
-    return kinship(G),G # G gets transposed, we'll turn this into an iterator (FIXME)
-
-def calculate_kinship_iter(geno):
-    """ 
-    Call the new kinship calculation where genotype_matrix contains
-    inds (columns) by snps (rows).
-    """
-    assert type(genotype_matrix) is iter
-    info("call genotype.normalize")
-    G = np.apply_along_axis( genotype.normalize, axis=0, arr=genotype_matrix)
-    info("call calculate_kinship_new")
-    return kinship(G)
-
-def calculate_kinship_old(genotype_matrix):
-    """
-    genotype_matrix is an n x m matrix encoding SNP minor alleles.
-    
-    This function takes a matrix oF SNPs, imputes missing values with the maf,
-    normalizes the resulting vectors and returns the RRM matrix.
-    
-    """
-    info("call calculate_kinship_old")
-    fatal("THE FUNCTION calculate_kinship_old IS OBSOLETE, use calculate_kinship_new instead - see Genotype Normalization Problem on Pjotr's blog")
-    n = genotype_matrix.shape[0]
-    m = genotype_matrix.shape[1]
-    info("genotype 2D matrix n (inds) is: %d" % (n))
-    info("genotype 2D matrix m (snps) is: %d" % (m))
-    assert m>n, "n should be larger than m (snps>inds)"
-    keep = []
-    mprint("G (before old normalize)",genotype_matrix)
-    for counter in range(m):
-        #print("type of genotype_matrix[:,counter]:", pf(genotype_matrix[:,counter]))
-        #Checks if any values in column are not numbers
-        not_number = np.isnan(genotype_matrix[:,counter])
-        
-        #Gets vector of values for column (no values in vector if not all values in col are numbers)
-        marker_values = genotype_matrix[True - not_number, counter]
-        #print("marker_values is:", pf(marker_values))
-        
-        #Gets mean of values in vector
-        values_mean = marker_values.mean()
-
-        genotype_matrix[not_number,counter] = values_mean
-        vr = genotype_matrix[:,counter].var()
-        if vr == 0:
-            continue
-        keep.append(counter)
-        genotype_matrix[:,counter] = (genotype_matrix[:,counter] - values_mean) / np.sqrt(vr)
-        progress('kinship_old normalize genotype',counter,m)
-        
-    genotype_matrix = genotype_matrix[:,keep]
-    mprint("G (after old normalize)",genotype_matrix.T)
-    kinship_matrix = np.dot(genotype_matrix, genotype_matrix.T) * 1.0/float(m)
-    return kinship_matrix,genotype_matrix
-    # return kinship_full(genotype_matrix.T),genotype_matrix
-
-def GWAS(pheno_vector,
-         genotype_matrix,
-         kinship_matrix,
-         kinship_eigen_vals=None,
-         kinship_eigen_vectors=None,
-         covariate_matrix=None,
-         restricted_max_likelihood=True,
-         refit=False,
-         temp_data=None):
-    """
-    Performs a basic GWAS scan using the LMM.  This function
-    uses the LMM module to assess association at each SNP and
-    does some simple cleanup, such as removing missing individuals
-    per SNP and re-computing the eigen-decomp
-
-    pheno_vector - n x 1 phenotype vector
-    genotype_matrix - n x m SNP matrix
-    kinship_matrix - n x n kinship matrix
-    kinship_eigen_vals, kinship_eigen_vectors = linalg.eigh(K) - or the eigen vectors and values for K
-    covariate_matrix - n x q covariate matrix
-    restricted_max_likelihood - use restricted maximum likelihood
-    refit - refit the variance component for each SNP
-    
-    """
-    if kinship_eigen_vals is None:
-        kinship_eigen_vals = []
-    if kinship_eigen_vectors is None:
-        kinship_eigen_vectors = []
-    
-    n = genotype_matrix.shape[0]
-    m = genotype_matrix.shape[1]
-
-    if covariate_matrix == None:
-        covariate_matrix = np.ones((n,1))
-
-    # Remove missing values in pheno_vector and adjust associated parameters
-    v = np.isnan(pheno_vector)
-    if v.sum():
-        keep = True - v
-        print(pheno_vector.shape,pheno_vector)
-        print(keep.shape,keep)
-        pheno_vector = pheno_vector[keep]
-        #genotype_matrix = genotype_matrix[keep,:]
-        #covariate_matrix = covariate_matrix[keep,:]
-        #kinship_matrix = kinship_matrix[keep,:][:,keep]
-        kinship_eigen_vals = []
-        kinship_eigen_vectors = []
-
-    lmm_ob = LMM(pheno_vector,
-                 kinship_matrix,
-                 kinship_eigen_vals,
-                 kinship_eigen_vectors,
-                 covariate_matrix)
-    if not refit:
-        lmm_ob.fit()
-
-    p_values = []
-    t_statistics = []
-    
-    n = genotype_matrix.shape[0]
-    m = genotype_matrix.shape[1]
-    
-    for counter in range(m):
-        x = genotype_matrix[:,counter].reshape((n, 1))
-        v = np.isnan(x).reshape((-1,))
-        if v.sum():
-            keep = True - v
-            xs = x[keep,:]
-            if xs.var() == 0:
-                p_values.append(0)
-                t_statistics.append(np.nan)
-                continue
-            
-            print(genotype_matrix.shape,pheno_vector.shape,keep.shape)
-
-            pheno_vector = pheno_vector[keep]
-            covariate_matrix = covariate_matrix[keep,:]
-            kinship_matrix = kinship_matrix[keep,:][:,keep]
-            lmm_ob_2 = LMM(pheno_vector,
-                           kinship_matrix,
-                           X0=covariate_matrix)
-            if refit:
-                lmm_ob_2.fit(X=xs)
-            else:
-                lmm_ob_2.fit()
-            ts, ps, beta, betaVar = lmm_ob_2.association(xs, REML=restricted_max_likelihood)
-        else:
-            if x.var() == 0:
-                p_values.append(0)
-                t_statistics.append(np.nan)
-                continue
-
-            if refit:
-                lmm_ob.fit(X=x)
-            ts, ps, beta, betaVar = lmm_ob.association(x, REML=restricted_max_likelihood)
-            
-        progress("gwas_old",counter,m)
-        
-        p_values.append(ps)
-        t_statistics.append(ts)
-
-    return t_statistics, p_values
-
-
-class LMM:
-
-    """
-          This is a simple version of EMMA/fastLMM.  
-          The main purpose of this module is to take a phenotype vector (Y), a set of covariates (X) and a kinship matrix (K)
-          and to optimize this model by finding the maximum-likelihood estimates for the model parameters.
-          There are three model parameters: heritability (h), covariate coefficients (beta) and the total
-          phenotypic variance (sigma).
-          Heritability as defined here is the proportion of the total variance (sigma) that is attributed to 
-          the kinship matrix.
- 
-          For simplicity, we assume that everything being input is a numpy array.
-          If this is not the case, the module may throw an error as conversion from list to numpy array
-          is not done consistently.
- 
-    """
-    def __init__(self,Y,K,Kva=[],Kve=[],X0=None,verbose=True):
- 
-       """
-       The constructor takes a phenotype vector or array of size n.
-       It takes a kinship matrix of size n x n.  Kva and Kve can be computed as Kva,Kve = linalg.eigh(K) and cached.
-       If they are not provided, the constructor will calculate them.
-       X0 is an optional covariate matrix of size n x q, where there are q covariates.
-       When this parameter is not provided, the constructor will set X0 to an n x 1 matrix of all ones to represent a mean effect.
-       """
-
-       if X0 is None: X0 = np.ones(len(Y)).reshape(len(Y),1)
-       self.verbose = verbose
- 
-       #x = Y != -9
-       x = True - np.isnan(Y)
-       #pdb.set_trace()
-       if not x.sum() == len(Y):
-          print("Removing %d missing values from Y\n" % ((True - x).sum()))
-          if self.verbose: sys.stderr.write("Removing %d missing values from Y\n" % ((True - x).sum()))
-          Y = Y[x]
-          print("x: ", len(x))
-          print("K: ", K.shape)
-          #K = K[x,:][:,x]
-          X0 = X0[x,:]
-          Kva = []
-          Kve = []
-       self.nonmissing = x
- 
-       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)
-          Kva,Kve = kvakve(K)
-          end = time.time()
-          if self.verbose: sys.stderr.write("Total time: %0.3f\n" % (end - begin))
-          print("sum(Kva),sum(Kve)=",sum(Kva),sum(Kve))
-
-       self.K = K
-       self.Kva = Kva
-       self.Kve = Kve
-       print("self.Kva is: ", self.Kva.shape, pf(self.Kva))
-       print("self.Kve is: ", self.Kve.shape, pf(self.Kve))
-       self.Y = Y
-       self.X0 = X0
-       self.N = self.K.shape[0]
-
-       # ----> Below moved to kinship.kvakve(K)
-       # if sum(self.Kva < 1e-6):
-       #    if self.verbose: sys.stderr.write("Cleaning %d eigen values\n" % (sum(self.Kva < 0)))
-       #    self.Kva[self.Kva < 1e-6] = 1e-6
- 
-       self.transform()
-
-    def transform(self):
- 
-         """
-            Computes a transformation on the phenotype vector and the covariate matrix.
-            The transformation is obtained by left multiplying each parameter by the transpose of the 
-            eigenvector matrix of K (the kinship).
-         """
-     
-         self.Yt = matrixMult(self.Kve.T, self.Y)
-         self.X0t = matrixMult(self.Kve.T, self.X0)
-         self.X0t_stack = np.hstack([self.X0t, np.ones((self.N,1))])
-         self.q = self.X0t.shape[1]
-
-    def getMLSoln(self,h,X):
- 
-       """
-          Obtains the maximum-likelihood estimates for the covariate coefficients (beta),
-          the total variance of the trait (sigma) and also passes intermediates that can 
-          be utilized in other functions. The input parameter h is a value between 0 and 1 and represents
-          the heritability or the proportion of the total variance attributed to genetics.  The X is the 
-          covariate matrix.
-       """
-
-       S = 1.0/(h*self.Kva + (1.0 - h))
-       Xt = X.T*S
-       XX = matrixMult(Xt,X)
-       XX_i = linalg.inv(XX)
-       beta =  matrixMult(matrixMult(XX_i,Xt),self.Yt)
-       Yt = self.Yt - matrixMult(X,beta)
-       Q = np.dot(Yt.T*S,Yt)
-       sigma = Q * 1.0 / (float(self.N) - float(X.shape[1]))
-       return beta,sigma,Q,XX_i,XX
- 
-    def LL_brent(self,h,X=None,REML=False): 
-       #brent will not be bounded by the specified bracket.
-       # I return a large number if we encounter h < 0 to avoid errors in LL computation during the search.
-       if h < 0: return 1e6
-       return -self.LL(h,X,stack=False,REML=REML)[0]
-         
-    def LL(self,h,X=None,stack=True,REML=False):
- 
-        """
-           Computes the log-likelihood for a given heritability (h).  If X==None, then the 
-           default X0t will be used.  If X is set and stack=True, then X0t will be matrix concatenated with
-           the input X.  If stack is false, then X is used in place of X0t in the LL calculation.
-           REML is computed by adding additional terms to the standard LL and can be computed by setting REML=True.
-        """
- 
-        if X is None:
-            X = self.X0t
-        elif stack: 
-            self.X0t_stack[:,(self.q)] = matrixMult(self.Kve.T,X)[:,0]
-            X = self.X0t_stack
- 
-        n = float(self.N)
-        q = float(X.shape[1])
-        beta,sigma,Q,XX_i,XX = self.getMLSoln(h,X)
-        LL = n*np.log(2*np.pi) + np.log(h*self.Kva + (1.0-h)).sum() + n + n*np.log(1.0/n * Q)
-        LL = -0.5 * LL
- 
-        if REML:
-            LL_REML_part = q*np.log(2.0*np.pi*sigma) + np.log(linalg.det(matrixMult(X.T,X))) - np.log(linalg.det(XX))
-            LL = LL + 0.5*LL_REML_part
-
-        return LL,beta,sigma,XX_i
-
-    def getMax(self,H, X=None,REML=False):
- 
-        """
-           Helper functions for .fit(...).  
-           This function takes a set of LLs computed over a grid and finds possible regions 
-           containing a maximum.  Within these regions, a Brent search is performed to find the 
-           optimum.
-  
-        """
-        n = len(self.LLs)
-        HOpt = []
-        for i in range(1,n-2):
-            if self.LLs[i-1] < self.LLs[i] and self.LLs[i] > self.LLs[i+1]: 
-                HOpt.append(optimize.brent(self.LL_brent,args=(X,REML),brack=(H[i-1],H[i+1])))
-                if np.isnan(HOpt[-1][0]):
-                    HOpt[-1][0] = [self.LLs[i-1]]
-
-        if len(HOpt) > 1: 
-            if self.verbose:
-                sys.stderr.write("NOTE: Found multiple optima.  Returning first...\n")
-            return HOpt[0]
-        elif len(HOpt) == 1:
-            return HOpt[0]
-        elif self.LLs[0] > self.LLs[n-1]:
-            return H[0]
-        else:
-            return H[n-1]
-
-    def fit(self,X=None,ngrids=100,REML=True):
- 
-        """
-            Finds the maximum-likelihood solution for the heritability (h) given the current parameters.
-            X can be passed and will transformed and concatenated to X0t.  Otherwise, X0t is used as 
-            the covariate matrix.
-   
-            This function calculates the LLs over a grid and then uses .getMax(...) to find the optimum.
-            Given this optimum, the function computes the LL and associated ML solutions.
-        """
-     
-        if X == None:
-            X = self.X0t
-        else: 
-            #X = np.hstack([self.X0t,matrixMult(self.Kve.T, X)])
-            self.X0t_stack[:,(self.q)] = matrixMult(self.Kve.T,X)[:,0]
-            X = self.X0t_stack
-  
-        H = np.array(range(ngrids)) / float(ngrids)
-        L = np.array([self.LL(h,X,stack=False,REML=REML)[0] for h in H])
-        self.LLs = L
-  
-        hmax = self.getMax(H,X,REML)
-        L,beta,sigma,betaSTDERR = self.LL(hmax,X,stack=False,REML=REML)
-        
-        self.H = H
-        self.optH = hmax
-        self.optLL = L
-        self.optBeta = beta
-        self.optSigma = sigma
-  
-        return hmax,beta,sigma,L
-
-    def association(self,X, h = None, stack=True,REML=True, returnBeta=True):
- 
-        """
-            Calculates association statitics for the SNPs encoded in the vector X of size n.
-            If h == None, the optimal h stored in optH is used.
-  
-        """
-        if stack: 
-           #X = np.hstack([self.X0t,matrixMult(self.Kve.T, X)])
-           self.X0t_stack[:,(self.q)] = matrixMult(self.Kve.T,X)[:,0]
-           X = self.X0t_stack
-           
-        if h == None:
-            h = self.optH
-  
-        L,beta,sigma,betaVAR = self.LL(h,X,stack=False,REML=REML)
-        q  = len(beta)
-        ts,ps = self.tstat(beta[q-1],betaVAR[q-1,q-1],sigma,q)
-        
-        if returnBeta:
-            return ts,ps,beta[q-1].sum(),betaVAR[q-1,q-1].sum()*sigma
-        return ts,ps
-
-    def tstat(self,beta,var,sigma,q): 
- 
-          """
-             Calculates a t-statistic and associated p-value given the estimate of beta and its standard error.
-             This is actually an F-test, but when only one hypothesis is being performed, it reduces to a t-test.
-          """
- 
-          ts = beta / np.sqrt(var * sigma)        
-          ps = 2.0*(1.0 - stats.t.cdf(np.abs(ts), self.N-q))
-          if not len(ts) == 1 or not len(ps) == 1:
-              print("ts=",ts)
-              print("ps=",ps)
-              raise Exception("Something bad happened :(")
-          return ts.sum(),ps.sum()
-
-    def plotFit(self,color='b-',title=''):
- 
-       """
-          Simple function to visualize the likelihood space.  It takes the LLs 
-          calcualted over a grid and normalizes them by subtracting off the mean and exponentiating.
-          The resulting "probabilities" are normalized to one and plotted against heritability.
-          This can be seen as an approximation to the posterior distribuiton of heritability.
- 
-          For diagnostic purposes this lets you see if there is one distinct maximum or multiple 
-          and what the variance of the parameter looks like.
-       """
-       import matplotlib.pyplot as pl
- 
-       mx = self.LLs.max()
-       p = np.exp(self.LLs - mx)
-       p = p/p.sum()
- 
-       pl.plot(self.H,p,color)
-       pl.xlabel("Heritability")
-       pl.ylabel("Probability of data")
-       pl.title(title)
-
-def run_gwas(species,n,m,k,y,geno,cov=None,reml=True,refit=False,inputfn=None,new_code=True):
-    """
-    Invoke pylmm using genotype as a matrix or as a (SNP) iterator.
-    """
-    info("run_gwas")
-    print('pheno', y)
-    
-    if species == "human" :
-        print('kinship', k )
-        ps, ts = run_human(pheno_vector = y,
-                           covariate_matrix = cov,
-                           plink_input_file = inputfn,
-                           kinship_matrix = k,
-                           refit = refit)
-    else:
-        print('geno', geno.shape, geno)
-
-        if new_code:
-            ps, ts = run_other_new(n,m,pheno_vector = y,
-                                   geno = geno,
-                                   restricted_max_likelihood = reml,
-                                   refit = refit)
-        else:
-            ps, ts = run_other_old(pheno_vector = y,
-                               genotype_matrix = geno,
-                               restricted_max_likelihood = reml,
-                               refit = refit)
-    return ps,ts
-
-def gwas_with_redis(key,species,new_code=True):
-    """
-    Invoke pylmm using Redis as a container. new_code runs the new
-    version. All the Redis code goes here!
-    """
-    json_params = Redis.get(key)
-    
-    params = json.loads(json_params)
-    
-    tempdata = temp_data.TempData(params['temp_uuid'])
-    def update_tempdata(loc,i,total):
-        """
-        This is the single method that updates Redis for percentage complete!
-        """
-        tempdata.store("percent_complete",round(i*100.0/total))
-        debug("Updating REDIS percent_complete=%d" % (round(i*100.0/total)))
-    progress_set_func(update_tempdata)
-
-    def narray(t):
-        info("Type is "+t)
-        v = params.get(t)
-        if v is not None:
-            # Note input values can be array of string or float
-            v1 = [x if x != 'NA' else 'nan' for x in v]
-            v = np.array(v1).astype(np.float)
-        return v
-
-    def marray(t):
-        info("Type is "+t)
-        v = params.get(t)
-        if v is not None:
-            m = []
-            for r in v:
-              # Note input values can be array of string or float
-              r1 = [x if x != 'NA' else 'nan' for x in r]
-              m.append(np.array(r1).astype(np.float))
-            return np.array(m)
-        return np.array(v)
-
-    def marrayT(t):
-        m = marray(t)
-        if m is not None:
-            return m.T
-        return m
-    
-    # We are transposing before we enter run_gwas - this should happen on the webserver
-    # side (or when reading data from file)
-    k = marray('kinship_matrix')
-    g = marrayT('genotype_matrix')
-    mprint("geno",g)
-    y = narray('pheno_vector')
-    n = len(y)
-    m = params.get('num_genotypes')
-    if m is None:
-      m = g.shape[0]
-    info("m=%d,n=%d" % (m,n))
-    ps,ts = run_gwas(species,n,m,k,y,g,narray('covariate_matrix'),params['restricted_max_likelihood'],params['refit'],params.get('input_file_name'),new_code)
-        
-    results_key = "pylmm:results:" + params['temp_uuid']
-
-    # fatal(results_key)
-    json_results = json.dumps(dict(p_values = ps,
-                                   t_stats = ts))
-    
-    #Pushing json_results into a list where it is the only item because blpop needs a list
-    Redis.rpush(results_key, json_results)
-    Redis.expire(results_key, 60*60)
-    return ps, ts
-
-def gn2_load_redis(key,species,kinship,pheno,geno,new_code=True):
-    """
-    This function emulates current GN2 behaviour by pre-loading Redis (note the input
-    genotype is transposed to emulate GN2 (FIXME!)
-    """
-    info("Loading Redis from parsed data")
-    if kinship == None:
-        k = None
-    else:
-        k = kinship.tolist()
-    params = dict(pheno_vector = pheno.tolist(),
-                  genotype_matrix = geno.T.tolist(),
-                  num_genotypes = geno.shape[0],
-                  kinship_matrix = k,
-                  covariate_matrix = None,
-                  input_file_name = None,
-                  restricted_max_likelihood = True,
-                  refit = False,
-                  temp_uuid = "testrun_temp_uuid",
-                        
-                  # meta data
-                  timestamp = datetime.datetime.now().isoformat())
-            
-    json_params = json.dumps(params)
-    Redis.set(key, json_params)
-    Redis.expire(key, 60*60)
-
-    return gwas_with_redis(key,species,new_code)
-
-def gn2_load_redis_iter(key,species,kinship,pheno,geno_iterator):
-    """
-    This function emulates GN2 behaviour by pre-loading Redis with
-    a SNP iterator, for this it sets a key for every genotype (SNP)
-    """
-    print("Loading Redis using a SNP iterator")
-    for i,genotypes in enumerate(geno_iterator):
-        gkey = key+'_geno_'+str(i)
-        Redis.set(gkey, genotypes)
-        Redis.expire(gkey, 60*60)
-    
-    if kinship == None:
-        k = None
-    else:
-        k = kinship.tolist()
-    params = dict(pheno_vector = pheno.tolist(),
-                  genotype_matrix = "iterator",
-                  num_genotypes = i,
-                  kinship_matrix = k,
-                  covariate_matrix = None,
-                  input_file_name = None,
-                  restricted_max_likelihood = True,
-                  refit = False,
-                  temp_uuid = "testrun_temp_uuid",
-                        
-                  # meta data
-                  timestamp = datetime.datetime.now().isoformat(),
-    )
-            
-    json_params = json.dumps(params)
-    Redis.set(key, json_params)
-    Redis.expire(key, 60*60)
-    return gwas_with_redis(key,species)
-
-# This is the main function used by Genenetwork2 (with environment)
-#
-# Note that this calling route will become OBSOLETE (we should use runlmm.py
-# instead)
-def gn2_main():
-    import argparse
-    parser = argparse.ArgumentParser(description='Run pyLMM')
-    parser.add_argument('-k', '--key')
-    parser.add_argument('-s', '--species')
-    
-    opts = parser.parse_args()
-    
-    key = opts.key
-    species = opts.species
-
-    gwas_with_redis(key,species)
-
-
-if __name__ == '__main__':
-    print("WARNING: Calling pylmm from lmm.py will become OBSOLETE, use runlmm.py instead!")
-    gn2_main()
-
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/lmm2.py b/wqflask/wqflask/my_pylmm/pyLMM/lmm2.py
deleted file mode 100644
index d871d8d2..00000000
--- a/wqflask/wqflask/my_pylmm/pyLMM/lmm2.py
+++ /dev/null
@@ -1,433 +0,0 @@
-# pylmm is a python-based linear mixed-model solver with applications to GWAS
-
-# Copyright (C) 2013,2014  Nicholas A. Furlotte (nick.furlotte@gmail.com)
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Affero General Public License as
-# published by the Free Software Foundation, either version 3 of the
-# License, or (at your option) any later version.
-
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU Affero General Public License for more details.
-
-# You should have received a copy of the GNU Affero General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-
-import sys
-import time
-import numpy as np
-from scipy.linalg import eigh, inv, det
-import scipy.stats as stats # t-tests
-from scipy import optimize
-from optmatrix import matrixMult
-import kinship
-
-sys.stderr.write("INFO: pylmm (lmm2) system path is "+":".join(sys.path)+"\n")
-sys.stderr.write("INFO: pylmm (lmm2) file is "+__file__+"\n")
-
-# ---- A trick to decide on the environment:
-try:
-    sys.stderr.write("INFO: lmm2 try loading module\n")
-    import utility.formatting # this is never used, just to check the environment
-    sys.stderr.write("INFO: This is a genenetwork2 environment (lmm2)\n")
-    from gn2 import uses, progress_set_func
-except ImportError:
-    # Failed to load gn2
-    has_gn2=False
-    import standalone as handlers
-    from standalone import uses, progress_set_func
-    sys.stderr.write("WARNING: LMM2 standalone version missing the Genenetwork2 environment\n")
-
-progress,mprint,debug,info,fatal = uses('progress','mprint','debug','info','fatal')
-
-
-def calculateKinship(W,center=False):
-      """
-	 W is an n x m matrix encoding SNP minor alleles.
-
-	 This function takes a matrix oF SNPs, imputes missing values with the maf,
-	 normalizes the resulting vectors and returns the RRM matrix.
-      """
-      n = W.shape[0]
-      m = W.shape[1]
-      keep = []
-      for i in range(m):
-	 mn = W[True - np.isnan(W[:,i]),i].mean()
-	 W[np.isnan(W[:,i]),i] = mn
-	 vr = W[:,i].var()
-	 if vr == 0: continue
-
-	 keep.append(i)
-	 W[:,i] = (W[:,i] - mn) / np.sqrt(vr)
-
-      W = W[:,keep]
-      K = matrixMult(W,W.T) * 1.0/float(m)
-      if center:
-	 P = np.diag(np.repeat(1,n)) - 1/float(n) * np.ones((n,n))
-	 S = np.trace(matrixMult(matrixMult(P,K),P))
-	 K_n = (n - 1)*K / S
-	 return K_n
-      return K
-
-def GWAS(Y, X, K, Kva=[], Kve=[], X0=None, REML=True, refit=False):
-      """
-
-        Performs a basic GWAS scan using the LMM.  This function
-        uses the LMM module to assess association at each SNP and 
-        does some simple cleanup, such as removing missing individuals 
-        per SNP and re-computing the eigen-decomp
-
-	Y - n x 1 phenotype vector 
-        X - n x m SNP matrix (genotype matrix)
-	K - n x n kinship matrix
-        Kva,Kve = linalg.eigh(K) - or the eigen vectors and values for K
-        X0 - n x q covariate matrix
-	REML - use restricted maximum likelihood 
-        refit - refit the variance component for each SNP
-
-      """
-      n = X.shape[0]
-      m = X.shape[1]
-      prins("Initialize GWAS")
-      print("genotype matrix n is:", n)
-      print("genotype matrix m is:", m)
-
-      if X0 is None: 
-         X0 = np.ones((n,1))
-      
-      # Remove missing values in Y and adjust associated parameters
-      v = np.isnan(Y)
-      if v.sum():
-	 keep = True - v
-	 keep = keep.reshape((-1,))
-	 Y = Y[keep]
-	 X = X[keep,:]
-	 X0 = X0[keep,:]
-	 K = K[keep,:][:,keep]
-	 Kva = []
-	 Kve = []
-
-      if len(Y) == 0: 
-         return np.ones(m)*np.nan,np.ones(m)*np.nan
-
-      L = LMM(Y,K,Kva,Kve,X0)
-      if not refit: L.fit()
-
-      PS = []
-      TS = []
-
-      n = X.shape[0]
-      m = X.shape[1]
-
-      for i in range(m):
-	 x = X[:,i].reshape((n,1))
-	 v = np.isnan(x).reshape((-1,))
-	 if v.sum():
-	    keep = True - v
-	    xs = x[keep,:]
-	    if xs.var() == 0: 
-	       PS.append(np.nan) 
-	       TS.append(np.nan) 
-	       continue
-
-	    Ys = Y[keep]
-	    X0s = X0[keep,:]
-	    Ks = K[keep,:][:,keep]
-	    Ls = LMM(Ys,Ks,X0=X0s)
-	    if refit: 
-               Ls.fit(X=xs)
-	    else: 
-               Ls.fit()
-	    ts,ps = Ls.association(xs,REML=REML)
-	 else: 
-	    if x.var() == 0: 
-	       PS.append(np.nan) 
-	       TS.append(np.nan) 
-	       continue
-
-	    if refit: 
-               L.fit(X=x)
-	    ts,ps = L.association(x,REML=REML)
-	    
-	 PS.append(ps)
-	 TS.append(ts)
-
-      return TS,PS
-
-class LMM2:
-
-   """This is a simple version of EMMA/fastLMM.
-
-   The main purpose of this module is to take a phenotype vector (Y),
-   a set of covariates (X) and a kinship matrix (K) and to optimize
-   this model by finding the maximum-likelihood estimates for the
-   model parameters.  There are three model parameters: heritability
-   (h), covariate coefficients (beta) and the total phenotypic
-   variance (sigma).  Heritability as defined here is the proportion
-   of the total variance (sigma) that is attributed to the kinship
-   matrix.
-
-   For simplicity, we assume that everything being input is a numpy
-   array.  If this is not the case, the module may throw an error as
-   conversion from list to numpy array is not done consistently.
-
-   """
-   def __init__(self,Y,K,Kva=[],Kve=[],X0=None,verbose=False):
-
-      """The constructor takes a phenotype vector or array Y of size n. It
-      takes a kinship matrix K of size n x n.  Kva and Kve can be
-      computed as Kva,Kve = linalg.eigh(K) and cached.  If they are
-      not provided, the constructor will calculate them.  X0 is an
-      optional covariate matrix of size n x q, where there are q
-      covariates.  When this parameter is not provided, the
-      constructor will set X0 to an n x 1 matrix of all ones to
-      represent a mean effect.
-      """
-
-      if X0 is None: 
-	 X0 = np.ones(len(Y)).reshape(len(Y),1)
-      self.verbose = verbose
-
-      x = True - np.isnan(Y)
-      x = x.reshape(-1,)
-      if not x.sum() == len(Y):
-	 if self.verbose: sys.stderr.write("Removing %d missing values from Y\n" % ((True - x).sum()))
-	 Y = Y[x]
-	 K = K[x,:][:,x]
-	 X0 = X0[x,:]
-	 Kva = []
-	 Kve = []
-      self.nonmissing = x
-
-      print("this K is:", K.shape, 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)
-          Kva,Kve = kinship.kvakve(K)
-          end = time.time()
-          if self.verbose: sys.stderr.write("Total time: %0.3f\n" % (end - begin))
-          print("sum(Kva),sum(Kve)=",sum(Kva),sum(Kve))
-
-      self.K = K
-      self.Kva = Kva
-      self.Kve = Kve
-      self.N = self.K.shape[0]
-      self.Y = Y.reshape((self.N,1))
-      self.X0 = X0
-
-      if sum(self.Kva < 1e-6):
-         if self.verbose: sys.stderr.write("Cleaning %d eigen values\n" % (sum(self.Kva < 0)))
-         self.Kva[self.Kva < 1e-6] = 1e-6
-
-      self.transform()
-
-   def transform(self):
-
-      """
-	 Computes a transformation on the phenotype vector and the covariate matrix.
-	 The transformation is obtained by left multiplying each parameter by the transpose of the 
-	 eigenvector matrix of K (the kinship).
-      """
-
-      self.Yt = matrixMult(self.Kve.T, self.Y)
-      self.X0t = matrixMult(self.Kve.T, self.X0)
-      self.X0t_stack = np.hstack([self.X0t, np.ones((self.N,1))])
-      self.q = self.X0t.shape[1]
-
-   def getMLSoln(self,h,X):
-
-      """
-	 Obtains the maximum-likelihood estimates for the covariate coefficients (beta),
-	 the total variance of the trait (sigma) and also passes intermediates that can 
-	 be utilized in other functions. The input parameter h is a value between 0 and 1 and represents
-	 the heritability or the proportion of the total variance attributed to genetics.  The X is the 
-	 covariate matrix.
-      """
-   
-      S = 1.0/(h*self.Kva + (1.0 - h))
-      Xt = X.T*S
-      XX = matrixMult(Xt,X)
-      XX_i = inv(XX)
-      beta =  matrixMult(matrixMult(XX_i,Xt),self.Yt)
-      Yt = self.Yt - matrixMult(X,beta)
-      Q = np.dot(Yt.T*S,Yt)
-      sigma = Q * 1.0 / (float(self.N) - float(X.shape[1]))
-      return beta,sigma,Q,XX_i,XX
-
-   def LL_brent(self,h,X=None,REML=False): 
-      #brent will not be bounded by the specified bracket.
-      # I return a large number if we encounter h < 0 to avoid errors in LL computation during the search.
-      if h < 0: return 1e6
-      return -self.LL(h,X,stack=False,REML=REML)[0]
-	 
-   def LL(self,h,X=None,stack=True,REML=False):
-
-      """
-	 Computes the log-likelihood for a given heritability (h).  If X==None, then the 
-	 default X0t will be used.  If X is set and stack=True, then X0t will be matrix concatenated with
-	 the input X.  If stack is false, then X is used in place of X0t in the LL calculation.
-	 REML is computed by adding additional terms to the standard LL and can be computed by setting REML=True.
-      """
-
-      if X is None: X = self.X0t
-      elif stack: 
-	 self.X0t_stack[:,(self.q)] = matrixMult(self.Kve.T,X)[:,0]
-	 X = self.X0t_stack
-
-      n = float(self.N)
-      q = float(X.shape[1])
-      beta,sigma,Q,XX_i,XX = self.getMLSoln(h,X)
-      LL = n*np.log(2*np.pi) + np.log(h*self.Kva + (1.0-h)).sum() + n + n*np.log(1.0/n * Q)
-      LL = -0.5 * LL
-
-      if REML:
-	 LL_REML_part = q*np.log(2.0*np.pi*sigma) + np.log(det(matrixMult(X.T,X))) - np.log(det(XX))
-	 LL = LL + 0.5*LL_REML_part
-
-
-      LL = LL.sum()
-      return LL,beta,sigma,XX_i
-
-   def getMax(self,H, X=None,REML=False):
-
-      """
-	 Helper functions for .fit(...).  
-	 This function takes a set of LLs computed over a grid and finds possible regions 
-	 containing a maximum.  Within these regions, a Brent search is performed to find the 
-	 optimum.
-
-      """
-      n = len(self.LLs)
-      HOpt = []
-      for i in range(1,n-2):
-          if self.LLs[i-1] < self.LLs[i] and self.LLs[i] > self.LLs[i+1]:
-	    HOpt.append(optimize.brent(self.LL_brent,args=(X,REML),brack=(H[i-1],H[i+1])))
-	    if np.isnan(HOpt[-1]): HOpt[-1] = H[i-1]
-	    #if np.isnan(HOpt[-1]): HOpt[-1] = self.LLs[i-1]
-	    #if np.isnan(HOpt[-1][0]): HOpt[-1][0] = [self.LLs[i-1]]
-
-      if len(HOpt) > 1: 
-	 if self.verbose: sys.stderr.write("NOTE: Found multiple optima.  Returning first...\n")
-	 return HOpt[0]
-      elif len(HOpt) == 1: return HOpt[0]
-      elif self.LLs[0] > self.LLs[n-1]: return H[0]
-      else: return H[n-1]
-
-
-   def fit(self,X=None,ngrids=100,REML=True):
-
-      """
-	 Finds the maximum-likelihood solution for the heritability (h) given the current parameters.
-	 X can be passed and will transformed and concatenated to X0t.  Otherwise, X0t is used as 
-	 the covariate matrix.
-
-	 This function calculates the LLs over a grid and then uses .getMax(...) to find the optimum.
-	 Given this optimum, the function computes the LL and associated ML solutions.
-      """
-      
-      if X is None: X = self.X0t
-      else: 
-	 #X = np.hstack([self.X0t,matrixMult(self.Kve.T, X)])
-	 self.X0t_stack[:,(self.q)] = matrixMult(self.Kve.T,X)[:,0]
-	 X = self.X0t_stack
-
-      H = np.array(range(ngrids)) / float(ngrids)
-      L = np.array([self.LL(h,X,stack=False,REML=REML)[0] for h in H])
-      self.LLs = L
-
-      hmax = self.getMax(H,X,REML)
-      L,beta,sigma,betaSTDERR = self.LL(hmax,X,stack=False,REML=REML)
-      
-      self.H = H
-      self.optH = hmax.sum()
-      self.optLL = L
-      self.optBeta = beta
-      self.optSigma = sigma.sum()
-
-      return hmax,beta,sigma,L
-
-   def association(self,X,h=None,stack=True,REML=True,returnBeta=False):
-      """
-	Calculates association statitics for the SNPs encoded in the vector X of size n.
-	If h is None, the optimal h stored in optH is used.
-
-      """
-      if False:
-         print "X=",X
-         print "h=",h
-         print "q=",self.q
-         print "self.Kve=",self.Kve
-         print "X0t_stack=",self.X0t_stack.shape,self.X0t_stack
-      
-      if stack:
-	 # X = np.hstack([self.X0t,matrixMult(self.Kve.T, X)])
-         m = matrixMult(self.Kve.T,X)
-         # print "m=",m
-         m = m[:,0]
-         self.X0t_stack[:,(self.q)] = m
-	 X = self.X0t_stack
-	 
-      if h is None: h = self.optH
-
-      L,beta,sigma,betaVAR = self.LL(h,X,stack=False,REML=REML)
-      q  = len(beta)
-      ts,ps = self.tstat(beta[q-1],betaVAR[q-1,q-1],sigma,q)
-      
-      if returnBeta: return ts,ps,beta[q-1].sum(),betaVAR[q-1,q-1].sum()*sigma
-      return ts,ps
-
-   def tstat(self,beta,var,sigma,q,log=False): 
-
-	 """
-	    Calculates a t-statistic and associated p-value given the estimate of beta and its standard error.
-	    This is actually an F-test, but when only one hypothesis is being performed, it reduces to a t-test.
-	 """
-
-	 ts = beta / np.sqrt(var * sigma)	 
-	 #ps = 2.0*(1.0 - stats.t.cdf(np.abs(ts), self.N-q))
-	 # sf == survival function - this is more accurate -- could also use logsf if the precision is not good enough
-	 if log:
-	    ps = 2.0 + (stats.t.logsf(np.abs(ts), self.N-q))
-	 else:
-	    ps = 2.0*(stats.t.sf(np.abs(ts), self.N-q))
-	 if not len(ts) == 1 or not len(ps) == 1: 
-	    raise Exception("Something bad happened :(")
-	 return ts.sum(),ps.sum()
-
-   def plotFit(self,color='b-',title=''):
-
-      """
-	 Simple function to visualize the likelihood space.  It takes the LLs 
-	 calcualted over a grid and normalizes them by subtracting off the mean and exponentiating.
-	 The resulting "probabilities" are normalized to one and plotted against heritability.
-	 This can be seen as an approximation to the posterior distribuiton of heritability.
-
-	 For diagnostic purposes this lets you see if there is one distinct maximum or multiple 
-	 and what the variance of the parameter looks like.
-      """
-      import matplotlib.pyplot as pl
-
-      mx = self.LLs.max()
-      p = np.exp(self.LLs - mx)
-      p = p/p.sum()
-
-      pl.plot(self.H,p,color)
-      pl.xlabel("Heritability")
-      pl.ylabel("Probability of data")
-      pl.title(title)
-   
-   def meanAndVar(self):
-
-      mx = self.LLs.max()
-      p = np.exp(self.LLs - mx)
-      p = p/p.sum()
-
-      mn = (self.H * p).sum()
-      vx = ((self.H - mn)**2 * p).sum()
-
-      return mn,vx
-
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/optmatrix.py b/wqflask/wqflask/my_pylmm/pyLMM/optmatrix.py
deleted file mode 100644
index 5c71db6a..00000000
--- a/wqflask/wqflask/my_pylmm/pyLMM/optmatrix.py
+++ /dev/null
@@ -1,55 +0,0 @@
-import sys
-import time
-import numpy as np
-from numpy.distutils.system_info import get_info;
-from scipy import linalg
-from scipy import optimize
-from scipy import stats
-
-useNumpy = None
-hasBLAS = None
-
-def matrix_initialize(useBLAS=True): 
-    global useNumpy  # module based variable
-    if useBLAS and useNumpy == None:
-        print get_info('blas_opt')
-        try:
-            linalg.fblas
-            sys.stderr.write("INFO: using linalg.fblas\n")
-            useNumpy = False
-            hasBLAS  = True
-        except AttributeError:
-            sys.stderr.write("WARNING: linalg.fblas not found, using numpy.dot instead!\n")
-            useNumpy = True
-    else:
-        sys.stderr.write("INFO: using numpy.dot\n")
-        useNumpy=True
-        
-def matrixMult(A,B):
-   global useNumpy  # module based variable
-
-   if useNumpy:
-       return np.dot(A,B)
-
-   # If the matrices are in Fortran order then the computations will be faster
-   # when using dgemm.  Otherwise, the function will copy the matrix and that takes time.
-   if not A.flags['F_CONTIGUOUS']:
-      AA = A.T
-      transA = True
-   else:
-      AA = A
-      transA = False
-
-   if not B.flags['F_CONTIGUOUS']:
-      BB = B.T
-      transB = True
-   else:
-      BB = B
-      transB = False
-
-   return linalg.fblas.dgemm(alpha=1.,a=AA,b=BB,trans_a=transA,trans_b=transB)
-
-def matrixMultT(M):
-    # res = np.dot(W,W.T)
-    # return linalg.fblas.dgemm(alpha=1.,a=M.T,b=M.T,trans_a=True,trans_b=False)
-    return matrixMult(M,M.T)
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/phenotype.py b/wqflask/wqflask/my_pylmm/pyLMM/phenotype.py
deleted file mode 100644
index 7b652515..00000000
--- a/wqflask/wqflask/my_pylmm/pyLMM/phenotype.py
+++ /dev/null
@@ -1,65 +0,0 @@
-# Phenotype routines
-
-# Copyright (C) 2013  Nicholas A. Furlotte (nick.furlotte@gmail.com)
-# Copyright (C) 2015  Pjotr Prins (pjotr.prins@thebird.nl)
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Affero General Public License as
-# published by the Free Software Foundation, either version 3 of the
-# License, or (at your option) any later version.
-
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU Affero General Public License for more details.
-
-# You should have received a copy of the GNU Affero General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-
-import sys
-import numpy as np
-
-# ---- A trick to decide on the environment:
-try:
-    from wqflask.my_pylmm.pyLMM import chunks
-    from gn2 import uses, progress_set_func
-except ImportError:
-    has_gn2=False
-    import standalone as handlers
-    from standalone import uses, progress_set_func
-
-progress,debug,info,mprint = uses('progress','debug','info','mprint')
-
-def remove_missing(n,y,g):
-    """
-    Remove missing data from matrices, make sure the genotype data has
-    individuals as rows
-    """
-    assert(y is not None)
-    assert y.shape[0] == g.shape[0],"y (n) %d, g (n,m) %s" % (y.shape[0],g.shape)
-
-    y1 = y
-    g1 = g
-    v = np.isnan(y)
-    keep = True - v
-    if v.sum():
-        info("runlmm.py: Cleaning the phenotype vector and genotype matrix by removing %d individuals...\n" % (v.sum()))
-        y1 = y[keep]
-        g1 = g[keep,:]
-        n = y1.shape[0]
-    return n,y1,g1,keep
-
-def remove_missing_new(n,y):
-    """
-    Remove missing data. Returns new n,y,keep
-    """
-    assert(y is not None)
-    y1 = y
-    v = np.isnan(y)
-    keep = True - v
-    if v.sum():
-        info("runlmm.py: Cleaning the phenotype vector by removing %d individuals" % (v.sum()))
-        y1 = y[keep]
-        n = y1.shape[0]
-    return n,y1,keep
-
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/plink.py b/wqflask/wqflask/my_pylmm/pyLMM/plink.py
deleted file mode 100644
index 7bd2df91..00000000
--- a/wqflask/wqflask/my_pylmm/pyLMM/plink.py
+++ /dev/null
@@ -1,107 +0,0 @@
-# Plink module
-#
-# Copyright (C) 2015  Pjotr Prins (pjotr.prins@thebird.nl)
-# Some of the BED file parsing came from pylmm:
-# Copyright (C) 2013  Nicholas A. Furlotte (nick.furlotte@gmail.com)
-
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Affero General Public License as
-# published by the Free Software Foundation, either version 3 of the
-# License, or (at your option) any later version.
-
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU Affero General Public License for more details.
-
-# You should have received a copy of the GNU Affero General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-
-# According to the PLINK information
-
-# Parse a textual BIM file and return the contents as a list of tuples
-# 
-# Extended variant information file accompanying a .bed binary genotype table.
-# 
-# A text file with no header line, and one line per variant with the following six fields:
-# 
-#     Chromosome code (either an integer, or 'X'/'Y'/'XY'/'MT'; '0' indicates unknown) or name
-#     Variant identifier
-#     Position in morgans or centimorgans (safe to use dummy value of '0')
-#     Base-pair coordinate (normally 1-based, but 0 ok; limited to 231-2)
-#     Allele 1 (corresponding to clear bits in .bed; usually minor)
-#     Allele 2 (corresponding to set bits in .bed; usually major)
-# 
-# Allele codes can contain more than one character. Variants with negative bp coordinates are ignored by PLINK. Example
-#
-#     1 mm37-1-3125499  0 3125499 1 2
-#     1 mm37-1-3125701  0 3125701 1 2
-#     1 mm37-1-3187481  0 3187481 1 2
-
-import struct
-# import numpy as np
-
-def readbim(fn):
-    res = []
-    for line in open(fn):
-        list = line.split()
-        if len([True for e in list if e == 'nan']) == 0:
-          res.append( (list[0],list[1],int(list[2]),int(list[3]),int(list[4]),int(list[5])) )
-        else:
-          res.append( (list[0],list[1],list[2],float('nan'),float('nan'),float('nan')) )
-    return res
-
-# .bed (PLINK binary biallelic genotype table)
-# 
-# Primary representation of genotype calls at biallelic variants. Must
-# be accompanied by .bim and .fam files. Basically contains num SNP
-# blocks containing IND (compressed 4 IND into a byte)
-#
-# Since it is a biallelic format it supports for every individual
-# whether the first allele is homozygous (b00), the second allele is
-# homozygous (b11), it is heterozygous (b10) or that it is missing
-# (b01).
-
-# http://pngu.mgh.harvard.edu/~purcell/plink2/formats.html#bed
-# http://pngu.mgh.harvard.edu/~purcell/plink2/formats.html#fam
-# http://pngu.mgh.harvard.edu/~purcell/plink2/formats.html#bim
-
-def readbed(fn,inds,encoding,func=None):
-
-    # For every SNP block fetch the individual genotypes using values
-    # 0.0 and 1.0 for homozygous and 0.5 for heterozygous alleles
-    def fetchGenotypes(X):
-        # D = { \
-        #      '00': 0.0, \
-        #      '10': 0.5, \
-        #      '11': 1.0, \
-        #      '01': float('nan') \
-        #   }
-
-        Didx = { '00': 0, '10': 1, '11': 2, '01': 3 }
-        G = []
-        for x in X:
-            if not len(x) == 10:
-                xx = x[2:]
-                x = '0b' + '0'*(8 - len(xx)) + xx
-            a,b,c,d = (x[8:],x[6:8],x[4:6],x[2:4]) 
-            L = [encoding[Didx[y]] for y in [a,b,c,d]]
-            G += L
-        G = G[:inds]
-        # G = np.array(G)
-        return G
-
-    bytes = inds / 4 + (inds % 4 and 1 or 0)
-    format = 'c'*bytes
-    count = 0
-    with open(fn,'rb') as f:
-        magic = f.read(3)
-        assert( ":".join("{:02x}".format(ord(c)) for c in magic) == "6c:1b:01")
-        while True:
-            count += 1
-            X = f.read(bytes)
-            if not X:
-                return(count-1)
-            XX = [bin(ord(x)) for x in struct.unpack(format,X)]
-            xs = fetchGenotypes(XX)
-            func(count,xs)
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py b/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py
deleted file mode 100644
index 6b241cd6..00000000
--- a/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py
+++ /dev/null
@@ -1,229 +0,0 @@
-# This is the LMM runner that calls the possible methods using command line
-# switches. It acts as a multiplexer where all the invocation complexity
-# is kept outside the main LMM routines.
-#
-# Copyright (C) 2015  Pjotr Prins (pjotr.prins@thebird.nl)
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Affero General Public License as
-# published by the Free Software Foundation, either version 3 of the
-# License, or (at your option) any later version.
-
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU Affero General Public License for more details.
-
-# You should have received a copy of the GNU Affero General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-
-from optparse import OptionParser
-import sys
-import tsvreader
-import numpy as np
-
-# Add local dir to PYTHONPATH
-import os
-cwd = os.path.dirname(__file__)
-if sys.path[0] != cwd:
-    sys.path.insert(1,cwd)
-
-# pylmm modules
-from lmm import gn2_load_redis, gn2_load_redis_iter, calculate_kinship_new, run_gwas
-from kinship import kinship, kinship_full
-import genotype
-import phenotype
-from standalone import uses
-
-progress,mprint,debug,info,fatal = uses('progress','mprint','debug','info','fatal')
-
-usage = """
-python runlmm.py [options] command
-
-  runlmm.py processing multiplexer reads standardised input formats
-  and calls the different routines (writes to stdout)
-
-  Current commands are:
-
-    parse        : only parse input files
-    redis        : use Redis to call into GN2
-    kinship      : calculate (new) kinship matrix
-
-  try --help for more information
-"""
-
-
-parser = OptionParser(usage=usage)
-# parser.add_option("-f", "--file", dest="input file",
-#                   help="In", metavar="FILE")
-parser.add_option("--kinship",dest="kinship",
-                  help="Kinship file format 1.0")
-parser.add_option("--pheno",dest="pheno",
-                  help="Phenotype file format 1.0")
-parser.add_option("--geno",dest="geno",
-                  help="Genotype file format 1.0")
-parser.add_option("--maf-normalization",
-                  action="store_true", dest="maf_normalization", default=False,
-                  help="Apply MAF genotype normalization")
-parser.add_option("--genotype-normalization",
-                  action="store_true", dest="genotype_normalization", default=False,
-                  help="Force genotype normalization")
-parser.add_option("--remove-missing-phenotypes",
-                  action="store_true", dest="remove_missing_phenotypes", default=False,
-                  help="Remove missing phenotypes")
-parser.add_option("-q", "--quiet",
-                  action="store_false", dest="verbose", default=True,
-                  help="don't print status messages to stdout")
-parser.add_option("--blas", action="store_true", default=False, dest="useBLAS", help="Use BLAS instead of numpy matrix multiplication")
-parser.add_option("-t", "--threads",
-                  type="int", dest="numThreads", 
-                  help="Threads to use")
-parser.add_option("--saveKvaKve",
-                  action="store_true", dest="saveKvaKve", default=False,
-                  help="Testing mode")
-parser.add_option("--test",
-                  action="store_true", dest="testing", default=False,
-                  help="Testing mode")
-parser.add_option("--test-kinship",
-                  action="store_true", dest="test_kinship", default=False,
-                  help="Testing mode for Kinship calculation")
-
-(options, args) = parser.parse_args()
-
-if len(args) != 1:
-    print usage
-    sys.exit(1)
-
-cmd = args[0]
-print "Command: ",cmd
-
-k = None
-y = None
-g = None
-
-if options.kinship:
-    k = tsvreader.kinship(options.kinship)
-    print k.shape
-
-if options.pheno:
-    y = tsvreader.pheno(options.pheno)
-    print y.shape
-
-if options.geno and cmd != 'iterator':
-    g = tsvreader.geno(options.geno)
-    print g.shape
-
-def check_results(ps,ts):
-    print np.array(ps)
-    print len(ps),sum(ps)
-    p1 = round(ps[0],4)
-    p2 = round(ps[-1],4)
-    if options.geno == 'data/small.geno':
-        info("Validating results for "+options.geno)
-        assert p1==0.7387, "p1=%f" % p1
-        assert p2==0.7387, "p2=%f" % p2
-    if options.geno == 'data/small_na.geno':
-        info("Validating results for "+options.geno)
-        assert p1==0.062, "p1=%f" % p1
-        assert p2==0.062, "p2=%f" % p2
-    if options.geno == 'data/test8000.geno':
-        info("Validating results for "+options.geno)
-        assert round(sum(ps)) == 4070
-        assert len(ps) == 8000
-    info("Run completed")
-
-if y is not None:
-    n = y.shape[0]
-
-if cmd == 'run':
-    if options.remove_missing_phenotypes:
-        raise Exception('Can not use --remove-missing-phenotypes with LMM2')
-    n = len(y)
-    m = g.shape[1]
-    ps, ts = run_gwas('other',n,m,k,y,g)  # <--- pass in geno by SNP
-    check_results(ps,ts)
-elif cmd == 'iterator':
-    if options.remove_missing_phenotypes:
-        raise Exception('Can not use --remove-missing-phenotypes with LMM2')
-    geno_iterator =  tsvreader.geno_iter(options.geno)
-    ps, ts = gn2_load_redis_iter('testrun_iter','other',k,y,geno_iterator)
-    check_results(ps,ts)
-elif cmd == 'redis_new':
-    # The main difference between redis_new and redis is that missing
-    # phenotypes are handled by the first
-    if options.remove_missing_phenotypes:
-        raise Exception('Can not use --remove-missing-phenotypes with LMM2')
-    Y = y
-    G = g
-    print "Original G",G.shape, "\n", G
-    # gt = G.T
-    # G = None
-    ps, ts = gn2_load_redis('testrun','other',k,Y,G,new_code=True)
-    check_results(ps,ts)
-elif cmd == 'redis':
-    # Emulating the redis setup of GN2
-    G = g
-    print "Original G",G.shape, "\n", G
-    if y is not None and options.remove_missing_phenotypes:
-        gnt = np.array(g).T
-        n,Y,g,keep = phenotype.remove_missing(n,y,gnt)
-        G = g.T
-        print "Removed missing phenotypes",G.shape, "\n", G
-    else:
-        Y = y
-    if options.maf_normalization:
-        G = np.apply_along_axis( genotype.replace_missing_with_MAF, axis=0, arr=g )
-        print "MAF replacements: \n",G
-    if options.genotype_normalization:
-        G = np.apply_along_axis( genotype.normalize, axis=1, arr=G)
-    g = None
-    gnt = None
-
-    # gt = G.T
-    # G = None
-    ps, ts = gn2_load_redis('testrun','other',k,Y,G, new_code=False)
-    check_results(ps,ts)
-elif cmd == 'kinship':
-    G = g
-    print "Original G",G.shape, "\n", G
-    if y != None and options.remove_missing_phenotypes:
-        gnt = np.array(g).T
-        n,Y,g,keep = phenotype.remove_missing(n,y,g.T)
-        G = g.T
-        print "Removed missing phenotypes",G.shape, "\n", G
-    if options.maf_normalization:
-        G = np.apply_along_axis( genotype.replace_missing_with_MAF, axis=0, arr=g )
-        print "MAF replacements: \n",G
-    if options.genotype_normalization:
-        G = np.apply_along_axis( genotype.normalize, axis=1, arr=G)
-    g = None
-    gnt = None
-
-    if options.test_kinship:
-        K = kinship_full(np.copy(G))
-        print "Genotype",G.shape, "\n", G
-        print "first Kinship method",K.shape,"\n",K
-        k1 = round(K[0][0],4)
-        K2,G = calculate_kinship_new(np.copy(G))
-        print "Genotype",G.shape, "\n", G
-        print "GN2 Kinship method",K2.shape,"\n",K2
-        k2 = round(K2[0][0],4)
-    
-    print "Genotype",G.shape, "\n", G
-    K3 = kinship(G)
-    print "third Kinship method",K3.shape,"\n",K3
-    sys.stderr.write(options.geno+"\n")
-    k3 = round(K3[0][0],4)
-    if options.geno == 'data/small.geno':
-        assert k1==0.8333, "k1=%f" % k1
-        assert k2==0.9375, "k2=%f" % k2
-        assert k3==0.9375, "k3=%f" % k3
-    if options.geno == 'data/small_na.geno':
-        assert k1==0.8333, "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:
-    fatal("Doing nothing")
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/standalone.py b/wqflask/wqflask/my_pylmm/pyLMM/standalone.py
deleted file mode 100644
index 40b2021d..00000000
--- a/wqflask/wqflask/my_pylmm/pyLMM/standalone.py
+++ /dev/null
@@ -1,110 +0,0 @@
-# Standalone specific methods and callback handler
-#
-# Copyright (C) 2015  Pjotr Prins (pjotr.prins@thebird.nl)
-#
-# Set the log level with
-#
-#   logging.basicConfig(level=logging.DEBUG)
-
-from __future__ import absolute_import, print_function, division
-
-import numpy as np
-import sys
-import logging
-
-# logger = logging.getLogger(__name__)
-logger = logging.getLogger('lmm2')
-logging.basicConfig(level=logging.DEBUG)
-np.set_printoptions(precision=3,suppress=True)
-
-progress_location = None 
-progress_current  = None
-progress_prev_perc     = None
-
-def progress_default_func(location,count,total):
-    global progress_current
-    value = round(count*100.0/total)
-    progress_current = value
-    
-progress_func = progress_default_func
-
-def progress_set_func(func):
-    global progress_func
-    progress_func = func
-    
-def progress(location, count, total):
-    global progress_location
-    global progress_prev_perc
-    
-    perc = round(count*100.0/total)
-    if perc != progress_prev_perc and (location != progress_location or perc > 98 or perc > progress_prev_perc + 5):
-        progress_func(location, count, total)
-        logger.info("Progress: %s %d%%" % (location,perc))
-        progress_location = location
-        progress_prev_perc = perc
-
-def mprint(msg,data):
-    """
-    Array/matrix print function
-    """
-    m = np.array(data)
-    if m.ndim == 1:
-        print(msg,m.shape,"=\n",m[0:3]," ... ",m[-3:])
-    if m.ndim == 2:
-        print(msg,m.shape,"=\n[",
-              m[0][0:3]," ... ",m[0][-3:],"\n ",
-              m[1][0:3]," ... ",m[1][-3:],"\n  ...\n ",
-              m[-2][0:3]," ... ",m[-2][-3:],"\n ",
-              m[-1][0:3]," ... ",m[-1][-3:],"]")
-
-def fatal(msg):
-    logger.critical(msg)
-    raise Exception(msg)
-
-def callbacks():
-    return dict(
-        write = sys.stdout.write,
-        writeln = print,
-        debug = logger.debug,
-        info = logger.info,
-        warning = logger.warning,
-        error = logger.error,
-        critical = logger.critical,
-        fatal = fatal,
-        progress = progress,
-        mprint = mprint
-    )
-
-def uses(*funcs):
-    """
-    Some sugar
-    """
-    return [callbacks()[func] for func in funcs]
-    
-# ----- Minor test cases:
-
-if __name__ == '__main__':
-    # logging.basicConfig(level=logging.DEBUG)
-    logging.debug("Test %i" % (1))
-    d = callbacks()['debug']
-    d("TEST")
-    wrln = callbacks()['writeln']
-    wrln("Hello %i" % 34)
-    progress = callbacks()['progress']
-    progress("I am half way",50,100)
-    list = [0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15,
-            0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15,
-            0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15,
-            0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15,
-            0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15]
-    mprint("list",list)
-    matrix = [[1,0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15],
-            [2,0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15],
-            [3,0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15],
-            [4,0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15],
-            [5,0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15],
-            [6,0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15]]
-    mprint("matrix",matrix)
-    ix,dx = uses("info","debug")
-    ix("ix")
-    dx("dx")
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/temp_data.py b/wqflask/wqflask/my_pylmm/pyLMM/temp_data.py
deleted file mode 100644
index 004d45c6..00000000
--- a/wqflask/wqflask/my_pylmm/pyLMM/temp_data.py
+++ /dev/null
@@ -1,25 +0,0 @@
-from __future__ import print_function, division, absolute_import
-from redis import Redis
-
-import simplejson as json
-
-class TempData(object):
-    
-    def __init__(self, temp_uuid):
-        self.temp_uuid = temp_uuid
-        self.redis = Redis()
-        self.key = "tempdata:{}".format(self.temp_uuid)
-        
-    def store(self, field, value):
-        self.redis.hset(self.key, field, value)
-        self.redis.expire(self.key, 60*15)  # Expire in 15 minutes
-        
-    def get_all(self):
-        return self.redis.hgetall(self.key)
-
-
-if __name__ == "__main__":
-    redis = Redis()
-    for key in redis.keys():
-        for field in redis.hkeys(key):
-            print("{}.{}={}".format(key, field, redis.hget(key, field)))
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/tsvreader.py b/wqflask/wqflask/my_pylmm/pyLMM/tsvreader.py
deleted file mode 100644
index 66b34ee2..00000000
--- a/wqflask/wqflask/my_pylmm/pyLMM/tsvreader.py
+++ /dev/null
@@ -1,122 +0,0 @@
-# Standard file readers
-#
-# Copyright (C) 2015  Pjotr Prins (pjotr.prins@thebird.nl)
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Affero General Public License as
-# published by the Free Software Foundation, either version 3 of the
-# License, or (at your option) any later version.
-
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU Affero General Public License for more details.
-
-# You should have received a copy of the GNU Affero General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-
-import sys
-import os
-import numpy as np
-import csv
-
-def kinship(fn):
-    K1 = []
-    print fn
-    with open(fn,'r') as tsvin:
-        assert(tsvin.readline().strip() == "# Kinship format version 1.0")
-        tsvin.readline()
-        tsvin.readline()
-        tsv = csv.reader(tsvin, delimiter='\t')
-        for row in tsv:
-            ns = np.genfromtxt(row[1:])
-            K1.append(ns) # <--- slow
-    K = np.array(K1)
-    return K
-
-def pheno(fn):
-    Y1 = []
-    print fn
-    with open(fn,'r') as tsvin:
-        assert(tsvin.readline().strip() == "# Phenotype format version 1.0")
-        tsvin.readline()
-        tsvin.readline()
-        tsvin.readline()
-        tsv = csv.reader(tsvin, delimiter='\t')
-        for row in tsv:
-            ns = np.genfromtxt(row[1:])
-            Y1.append(ns) # <--- slow
-    Y = np.array(Y1)
-    return Y
-
-def geno(fn):
-    G1 = []
-    hab_mapper = {'A':0,'H':1,'B':2,'-':3}
-    pylmm_mapper = [ 0.0, 0.5, 1.0, float('nan') ]
-
-    print fn
-    with open(fn,'r') as tsvin:
-        line = tsvin.readline().strip()
-        assert line == "# Genotype format version 1.0", line
-        tsvin.readline()
-        tsvin.readline()
-        tsvin.readline()
-        tsvin.readline()
-        tsv = csv.reader(tsvin, delimiter='\t')
-        for row in tsv:
-            # print(row)
-            id = row[0]
-            gs = list(row[1])
-            # print id,gs
-            gs2 = [pylmm_mapper[hab_mapper[g]] for g in gs]
-            # print id,gs2
-            # ns = np.genfromtxt(row[1:])
-            G1.append(gs2) # <--- slow
-    G = np.array(G1)
-    return G
-
-def geno(fn):
-    G1 = []
-    for id,values in geno_iter(fn):
-        G1.append(values) # <--- slow
-    G = np.array(G1)
-    return G
-
-def geno_callback(fn,func):
-    hab_mapper = {'A':0,'H':1,'B':2,'-':3}
-    pylmm_mapper = [ 0.0, 0.5, 1.0, float('nan') ]
-
-    print fn
-    with open(fn,'r') as tsvin:
-        assert(tsvin.readline().strip() == "# Genotype format version 1.0")
-        tsvin.readline()
-        tsvin.readline()
-        tsvin.readline()
-        tsvin.readline()
-        tsv = csv.reader(tsvin, delimiter='\t')
-        for row in tsv:
-            id = row[0]
-            gs = list(row[1])
-            gs2 = [pylmm_mapper[hab_mapper[g]] for g in gs]
-            func(id,gs2) 
-
-def geno_iter(fn):
-    """
-    Yield a tuple of snpid and values
-    """
-    hab_mapper = {'A':0,'H':1,'B':2,'-':3}
-    pylmm_mapper = [ 0.0, 0.5, 1.0, float('nan') ]
-
-    print fn
-    with open(fn,'r') as tsvin:
-        assert(tsvin.readline().strip() == "# Genotype format version 1.0")
-        tsvin.readline()
-        tsvin.readline()
-        tsvin.readline()
-        tsvin.readline()
-        tsv = csv.reader(tsvin, delimiter='\t')
-        for row in tsv:
-            id = row[0]
-            gs = list(row[1])
-            gs2 = [pylmm_mapper[hab_mapper[g]] for g in gs]
-            yield (id,gs2)