import fastlmm.association.lrt as lr
import fastlmm.association.score as score
import fastlmm.association.testCV as testCV
import fastlmm.inference as fastlmm
from fastlmm.pyplink.plink import *
from pysnptools.util.pheno import *
from fastlmm.pyplink.altset_list import *
import subprocess, sys, os.path
import fastlmm.util.stats.chi2mixture as c2
from fastlmm.util.distributable import *
from fastlmm.util.runner import *
import pdb
import fastlmm.util.util as utilx
import fastlmm.util.stats as ss
import time
from Result import *
from PairResult import *
from fastlmm.association.tests import *
import fastlmm.util.genphen as gp
import scipy as sp
from itertools import *
from fastlmm.pyplink.snpreader.Bed import *
import logging
import warnings
from tempfile import TemporaryFile
import fastlmm.util.preprocess as util

class FastLmmSet: # implements IDistributable
    '''
    A class for specifying a FastLmmSet and then running it.
    '''

    # member variables with default values, or those not defined in the input file    
    mpheno=1
    mindist = 2.0
    idist=1
    autoselect = True
    extractSim = None
    nperm = 10 # used to permute the data in order to get empirical distrib of null for param fitting of null; these are global
    npermabs = None #used instead of nperm to compute an at least this many permutations (I say at least because it is # of perms per test)
    cache_from_perm=False #If true, then actually does the permutations with the real, all in run_test() so can cache expensive things
    nlocalperm=None #same as nperm, but local to each test (doesn't do pooling)
    fitlocal=False #if true, fit aUD to each local set
    write_lrtperm = False  #when set to True, writes a third file containing only the p-values resulting from nperm
    nullfitfile = None   #file from which to grab the null stats from. should be generated using write_lrtperm
    calseed = None# was 0 # the int that gets added to the random seed 34343 used to generate permutations for lrt null fitting
    sets=None
    permute=None       # permute the SNPs being tested (or, equivalently, the compliment) using this with 1,2,3,4, etc. for different seeds, used for type I error tables
    forcefullrank = False
    nullfit="qq"       #alternative is "ml"    
    qmax=0.1
    minsetsize = 1
    maxsetsize = 10000000
    datestamp = None #use auto to have it auto-generate a date stamp
    covarimp = 'standardize' #use 'standardize' for mean and variance standardization, with mean imputation
    altset_list2=None#for doing pairs of sets
    nullModel = None
    altModel = None
    genphen=None#MUST LEAVE THIS AS THE DEFAULT#{"varE":1.0,"varG":1.0, "varBack":1.0,"varCov":1.0,"link":'linear',"casefrac":0.5,"once":false}      #generate synthetic phen, using LMM, and real SNPs
       
    scoring = None
    greater_is_better = None
    log = None
    detailed_table = False
    signal_ratio = True

    _synthphenfile=None

    alt_snpreader = None

    def addpostifx_to_outfile(self):
        if self.datestamp is not None:
            if self.datestamp=="auto":
                self.datestamp=utilx.datestamp()
            self.outfile=utilx.appendtofilename(self.outfile,self.datestamp,"_")

    def __init__(self, **entries):
        '''
        outfile         : string of the filename to which to write out results (two files will be generated from this, on a record of the run)
        phenofile       : string of the filename containing phenotype
        alt_snpreader   : A snpreader for the SNPs for alternative kernel. If just a file name is given, the Bed reader (with standardization) is used.
        bedfilealt      : (deprecated) same as alt_snpreader
        filenull        : string of the filename (with .bed or .ped suffix) containing SNPs for null kernel. Should be Bed format if autoselect is true.
        extractSim      : string of a filename containing a list of whitespace-delmited SNPs. Only these SNPs will be used when filenull is read.
                               By default all SNPs in filenull are used.
                               It is an error to specify extractSim and not filenull.
                               Currently assumes that filenull is in bed format and that autoselect is not being used.
        altset_list     : list of the altsets
                               By default this is a file that will be read by SnpAndSetNameCollection
                               but a file of nuc ranges can be given again via 'NucRangeSet(filename)'
        altset_list2    : if this is present, then does an interaction test between genes in altset_list and alt_setlist2
        covarfile       : string of the filename containing covariate
        mpheno=1        : integer representing the index of the testing phenotype (starting at 1)
        ipheno=0        : same as mpheno, but starts at 0 (deprecated)
        mindist=2.0     : SNPs within mindist from the alternative SNPs will be removed from
                          the null kernel computation
        idist=1         : the index of the position index to use with mindist
                            1 : genomic distance
                            2 : base-pair distance
        autoselect = True   : Should autoselect to used?
        test="lrt"      :['sc_davies', 'sc_mom', 'sc_liuskat','sc_daviesskat','sc_all','lrt','lrt_up']
	    nperm = 10       : number of pemutations per test
        npermabs = None : absolutely number of permutations
        calseed =   : the int that gets added to the random seed 34343 used to generate permutations for lrt null fitting.
        write_lrtperm   : If true, write the lrtperm vector (dictated by calseed) to a file.
        forcefullrank   :  If true, the covariance is always computed for lrt (default False)
        minsetsize      : If non-None, only include sets at least this large (inclusive)
        maxsetsize      : If non-None, only include sets no more than this large (inclusive)
        datestamp       : Defaults to None. Use "auto" to have it generate a unique date and time string
                          to append to filenames. Otherwise, it appends whatever you set it to.
        log             : (Defaults to not changing logging level) Level of log messages, e.g. logging.CRITICAL, logging.ERROR, logging.WARNING, logging.INFO
        '''
        
        if hasattr(self,"ipheno"):
            assert not hasattr(self,"mpheno")
            self.mpheno = self.ipheno+1

        self.__dict__.update(entries)
        self._check_entries(entries)           
        self._ran_once = False
        self.addpostifx_to_outfile()

        if self.hasNonNoneAttr("alt_snpreader") and self.hasNonNoneAttr("bedfilealt") : raise Exception("alt_snpreader or bedfilealt may be given, but not both")
        if self.hasNonNoneAttr("bedfilealt"):
            warnings.warn("bedfilealt is deprecated. alt_snpreader instead", DeprecationWarning)
            self.alt_snpreader = self.bedfilealt
            delattr(self,"bedfilealt")
        if isinstance(self.alt_snpreader, str):
            self.alt_snpreader = Bed(self.alt_snpreader)


        # convert deprecated "verbose" into "log"        
        if self.hasNonNoneAttr("log") and self.hasNonNoneAttr("verbose") :
            raise Exception("log or verbose may be given, but not both")
        if self.hasNonNoneAttr("verbose"):
            if self.verbose:
                self.log = logging.INFO
            else:
                self.log = logging.CRITICAL
            delattr(self, "verbose")

        if self.hasNonNoneAttr("log"): #If neither were set
            logger = logging.getLogger()
            logger.setLevel(self.log)


        if not hasattr(self, "test") : raise Exception("FastLmmSet must have 'test' set")
        #allowed_tests=['sc_davies', 'sc_mom', 'sc_liuskat','sc_daviesskat','sc_all','lrt','cv_l1','cv_l2']
        #if self.test not in allowed_tests: raise Exception("test must be one of " + str(allowed_tests) +", but found " + self.test)

    # =================================================================
    # BEGIN Implementation of the IDistributable interface described in fastlmm.util.distributable
    # =================================================================

    def hasNonNoneAttr(self, attr):
        return hasattr(self,attr) and getattr(self,attr) != None


    # required by IDistributable
    @property
    def work_count(self):
        '''
        Tells how many work items there are
        '''
        self.run_once() #load files, etc. -- stuff we only want to do once per task (e.g. on the cluster)
                        #no matter how many times we call 'run_test' (and then 'reduce' which of course
                        #only gets run one time, and calls this too, but doesn't actually do the work).

        if self.altset_list2 is None: #singleton sets
            return len(self.altsetlist_filtbysnps) * (1+self.nperm)
        else: #pairs of sets
            return len(self.altsetlist_filtbysnps)*len(self.__altsetlist2_filtbysnps)* (1+self.nperm)   

    def printPhenToFile(self, nInd, SNPsalt, y):
        self._synthphenfile=utilx.appendtofilename(self.outfile,"synthPhenVarBack")
        np.savetxt(self._synthphenfile,y)
        with open(self._synthphenfile,"w") as fp: 
            for n in sp.arange(0,nInd):
                 fp.write(SNPsalt['iid'][n][0] +'\t'+ SNPsalt['iid'][n][1]+ "\t" + '%f' % y[n] + "\n")

    def work_sequence(self):
        '''
        Enumerates a sequence of work items
        Each work item is a lambda expression (i.e. function pointer) that calls 'run_test', returning a list of Results (often just one)
        '''
        self.run_once() #load files, etc. -- stuff we only want to do once per task (e.g. on the cluster)
                        #no matter how many times we call 'run_test' (and then 'reduce' which of course
                        #only gets run one time, and calls this too, but doesn't actually do the work).
        ttt0=time.time()
        y=None        
        y_back=None
        nInd=len(self.__y)               
        haswrittenphen=False;                            
                
        if self.altset_list2 is None: #singleton sets            
            for iset, altset in enumerate(self.altsetlist_filtbysnps):
                for iperm in xrange(-1, self.nperm):   #note that self.nperm is the 'stop', not the 'count'
                    SNPsalt=altset.read()
                    SNPsalt['snps'] = util.standardize(SNPsalt['snps'])
                    G1 = SNPsalt['snps']/sp.sqrt(SNPsalt['snps'].shape[1])  
                    ichrm =  ",".join(sp.array(sp.unique(SNPsalt['pos'][:,0]),dtype=str)) 
                    minpos= str(sp.min(SNPsalt['pos'][:,2]))
                    maxpos= str(sp.max(SNPsalt['pos'][:,2]))
                    iposrange = minpos + "-" + maxpos                     
                    
                    if self.genphen is None:
                        y=self.__y
                    else: 
                        assert self.permute is None, "Error: using permute with genphen--there should be no need, use genphen['seed'] instead"
                        self.__y=None      

                        #if (self.__G0 is not None) and y_G0 is None:  #cache so only do it once
                        if (self.__SNPs0 is not None) and y_G0 is None:  #cache so only do it once   
                            newseed = self.mainseed ^ self.genphen['seed']
                            from numpy.random import RandomState
                            randomstate = RandomState(newseed)
                            nSnp=self.__SNPs0['data']['snps'].shape[1]
                            #good for low rank, other wise, use the dual, as in gp.genphen
                            raise Exception("Bug below. should be randn, is rand")
                            y_G0=sp.sqrt(self.genphen["varBack"]/nSnp)*self.__SNPs0['data']['snps'].dot(randomstate.rand(nSnp,1))    #TODO: CL This seems to be a bug. Should be randn
                            #y_G0=sp.sqrt(self.genphen["varBack"])*self.__SNPs0['data']['snps'].dot(randomstate.rand(nSnp,1))    #TODO: CL This seems to be a bug. Should be randn
                        elif self.__SNPs0 is None:
                            y_G0=0      

                        #always have the same background signal                                                
                        if self.genphen["varBackNullFileGen"] is not None and y_back is None:                                                      
                            nSnp=self.__SNPs0['data']['snps'].shape[1]
                            y_back=sp.sqrt(self.genphen["varBack"]/nSnp)*self.__varBackNullSnpsGen['snps'].dot(sp.random.randn(self.__varBackNullSnpsGen['snps'].shape[1],1))/sp.sqrt(self.__varBackNullSnpsGen['snps'].shape[1]) 
                            #self.printPhenToFile(nInd, SNPsalt, y_back);
                        elif self.genphen.has_key("varBackNullPhenGen") and self.genphen["varBackNullPhenGen"] is not None and y_back is None:
                            y_back=loadPhen(filename = self.genphen["varBackNullPhenGen"])['vals']                            
                        elif y_back is None:          
                            y_back=0                          
                        
                        if self.genphen["once"] and y is None: # only generate the phenotype once for entire run                            
                            newseed = self.mainseed ^ self.genphen['seed']  
                            assert  self.genphen["varG"]==0, "doesn't make sense to have varG>0 and only one phen--not sure why I put this code there"       
                            if self.genphen["varG"]>0 and nInd<=SNPsalt['snps'].shape[1]:                                                       
                                Kall = self.KfromAltSnps(nInd) #useful for full rank
                                y=gp.genphen(y_G0=y_G0+y_back,G1=None,covDat=self.__X,options=self.genphen,nInd=nInd,
                                             K1=Kall,randseed=newseed)                          
                            else:     
                                Kall=None                                     
                                y=gp.genphen(y_G0=y_G0+y_back,G1=G1,covDat=self.__X,options=self.genphen,nInd=nInd,randseed=newseed)
                        elif (not self.genphen["once"]):  #generate for each set in turn                           
                            eachseed=utilx.combineseeds(iset,self.genphen['seed'])                            
                            newseed = self.mainseed ^ eachseed                              
                            y=gp.genphen(y_G0=y_G0+y_back,G1=G1,covDat=self.__X,options=self.genphen,nInd=nInd,randseed=newseed) 
                                                                
                    assert y is not None, "y is None"                                   
                    yield lambda altset=altset,iset=iset,iperm=iperm,y=y,ichrm=ichrm, iposrange=iposrange : self.run_test(SNPs1=SNPsalt,G1=G1, y=y, altset=altset, iset=iset, iperm=iperm, ichrm=ichrm, iposrange=iposrange)
        else: #pairs of sets
            raise Exception("not implemented, started a long time ago and never finished")
            #for iperm in xrange(-1, self.nperm):   #note that self.nperm is the 'stop', not the 'count'
            #    for iset, altset in enumerate(self.__altsetlist_filtbysnps):
            #         for iset2, altset2 in enumerate(self.__altsetlist2_filtbysnps):
            #            if iset!=iset2:
            #                yield lambda altset=altset,iset=iset,altset2=altset2,iset2=iset2,iperm=iperm : self.run_interactiontest(altset,iset,altset2,iset2,iperm)

        ttt1=time.time()
        logging.info("---------------------------------------------------")
        logging.info("Elapsed time for all tests is %.2f seconds" % (ttt1-ttt0))
        logging.info("---------------------------------------------------")
    
    def check_for_None(self, alteqnull, alteqnullperm):
        none_ind = []
        for i in range(len(alteqnull)):
            if alteqnull[i] is None:
                none_ind.append(i)
        none_ind_perm = []
        for i in range(len(alteqnullperm)):
            if alteqnullperm[i] is None:
                none_ind_perm.append(i)
        if len(none_ind)>0:
            raise Exception("found None entries in alteqnull during reduce, which means some Result objects were missing")
        if len(none_ind_perm)>0:
            raise Exception("found None entries in alteqnullperm during reduce, which means some Result objects were missing")

    # required by IDistributable
    def reduce(self, result_list_sequence):
        '''
        Given a sequence of results from 'run_test', create the output report.

        '''
        logging.info("last modified: %s" % time.ctime(os.path.getmtime(__file__)))
        logging.info("-------------------------------------------------")

        # Create the info file and direct all 'stdout' messages to both stdout and to the info file
        infofile=utilx.appendtofilename(self.outfile,"info")
        outfiletab=self.outfile
        if self.datestamp=="auto": self.outfile=utilx.appendtofilename(self.outfile,"tab")

        try:
            utilx.create_directory_if_necessary(infofile)
        except:
            logging.warn("Exception while creating directory for '{0}'. Assuming that other cluster task is creating it.".format(infofile))

        logging_handler=logging.FileHandler(infofile,"w",delay=False)
        logger = logging.getLogger()
        logger.addHandler(logging_handler)

        logging.info("distributable = " + self.__repr__())

        self.run_once() #load files, etc. -- doesn't actually do the work if it's already been done
                         #note, however, that before reduce is called, work_count calls run_once() anyhow

        result_dict = {}

        lrt = SP.nan*SP.ones(len(self.altsetlist_filtbysnps))
        # whether alt and null models give the same margll
        alteqnull = [None]*len(self.altsetlist_filtbysnps)

        lrtperm = SP.nan*SP.ones(len(self.altsetlist_filtbysnps)*self.nperm)
        alteqnullperm = [None]*len(self.altsetlist_filtbysnps)*self.nperm
        setsizeperm = SP.nan*SP.ones(len(self.altsetlist_filtbysnps)*self.nperm)

        npvals=self.test.npvals
        pv_adj= sp.nan*sp.ones((len(self.altsetlist_filtbysnps)))
              
        # results can come in any order, so we have to use iperm and iset to put them in the right place
        # there is one result instance for each combination of test and permutation, and here we are just gathering them
        # into the arrays from above

        for result_list in result_list_sequence:
            for result in result_list:
                if result.iperm < 0:
                    result_dict[result.iset] = result
                    #iset is the index of the test (irrespective of permutation)
                    lrt[result.iset] = self.test.lrt(result)
                    alteqnull[result.iset] = result.alteqnull #equiv to result["alteqnull"]
                    #pv_adj[result.iset,:] = self.test.pv_adj_from_result(result)                
                    pv_adj[result.iset] = self.test.pv_adj_from_result(result)                
                else:
                    isetiperm = result.iset+result.iperm*len(self.altsetlist_filtbysnps)
                    lrtperm[isetiperm] = self.test.lrt(result)
                    alteqnullperm[isetiperm] = result.alteqnull
                    setsizeperm[isetiperm] = result.setsize
                if result.alteqnull is None and str(self.test)[0:3]=="lrt":
                    raise Exception("self.alteqnull is None")

        #look for None to see if anything didn't get filled in, and if so, where
        if str(self.test)[0:3]=="lrt":
            self.check_for_None(alteqnull, alteqnullperm)
            
        if self.write_lrtperm and self.nperm>0:
            self._saveArray("lrtperm", ["2*(LL(alt)-LL(null))","alteqnull","setsize"], (lrtperm,alteqnullperm,setsizeperm))
                
        pv_adj,ind = self.test.pv_adj_and_ind(self.nperm, pv_adj, self.nullfit, 
                                              lrt, lrtperm, alteqnull, alteqnullperm,self.qmax, self.nullfitfile, self.nlocalperm)

        logging.info("writing the result files : " + self.outfile + "")

        #TODO: move the following to a function/object:
        with open(outfiletab,"w") as fp:
            if type(self.test) is Lrt: #dan: this is ugly, will fix soon                
                self.test.write(fp, ind, result_dict, pv_adj, self.detailed_table, self.signal_ratio)
            else:
                self.test.write(fp, ind, result_dict, pv_adj, self.detailed_table)

        logger.removeHandler(logging_handler)
        logging_handler.close()
        return self.outfile
        
    def getRandSnpSignal(self, nSnp, nInd,genphen,newseed):
        from numpy.random import RandomState
        import fastlmm.util.gensnp as gp
        randomstate = RandomState(newseed) 
        nSnp=genphen["numBackSnps"]
        #randsnps=self.alt_snpreader.read(RandomSnpSet(nSnp,newseed))     #this appears to be VERY slow
        #snps=randsnps['snps']         
        snps=gp.gensnps(nInd,nSnp)                        
        randG = snps/sp.sqrt(nSnp) #pre-process for kernel
        y_randG=sp.sqrt(genphen["varBack"])*randG.dot(randomstate.rand(nSnp,1))            
        return y_randG

    def TESTBEFOREUSINGKfromAltSnps(self, N, SNPsalt=None):        
       
        print "constructing K from all SNPs"
        t0=time.time()

        Kall=sp.zeros([N,N])
        nSnpTotal=self.alt_snpreader.snp_count
        print "reading in " +str(nSnpTotal)+ " SNPs and adding up kernels"
        #altsnps = readBED(self.bedfilealt,standardizeSNPs=True)['snps'] #loads all in to memory        
        blocksize=100
        ct=0
        ts=time.time()
        if self.alt_snpreader.snp_count<N: raise Exception("need to adjust code to handle low rank")          
        for start in range(0,self.alt_snpreader.snp_count,blocksize):                                
            ct+=blocksize
            snpSet = PositionRange(start,blocksize)
            snps = self.alt_snpreader.read(snpSet)['snps']
            import fastlmm.util.preprocess as util
            snps = util.standardize(snps)

            #print "start = {0}".format(start)
            Kall=Kall + snps.dot(snps.T)    
            t1=time.time()
            if ct % 50000==0: 
                print "read %s SNPs in %.2f seconds" % (ct, t1-ts)
                ts=time.time()
            #if ct==2: 
            #    break
            #    Kall=sp.rand((N,N))
            #    Kall=Kall.dot(Kall.T)
        Kall=Kall/sp.sqrt(self.alt_snpreader.snp_count)  
        Kall=Kall + 1e-5*sp.eye(N,N)     
        t1=time.time()
        logging.info("%.2f seconds elapsed" % (t1-t0))        
                
        return Kall

    def _saveArray(self, appendNameFile, header, values):
        outfile = utilx.appendtofilename(self.outfile, appendNameFile)
        try:
            utilx.create_directory_if_necessary(outfile)
        except:
            logging.warn("Exception while creating directory for '{0}'. Assuming that other cluster task is creating it.".format(outfile))

        logging.info("writing to file " + outfile + ".")

        assert type(values) == sp.ndarray or type(values) == tuple, 'values can only be sp.ndarray or a tuple of sp.ndarray.'

        if type(values) == sp.ndarray:
            values = (values,)

        values = map( list, zip(*values) )

        with open(outfile, "w") as fp:
            fp.write('\t'.join(header) + "\n")
            for v in values:
                fp.write("{0}\n".format( '\t'.join([str(x) for x in v]) ))
    
    #!! would be nice of this was optional and if not given the OS was asked
    # required by IDistributable
    @property
    def tempdirectory(self):
        return self.outfile + ".work_directory"


    #!! need comments
    def copyinputs(self, copier):
        copier.input(self.phenofile)
        copier.input(self.alt_snpreader)
        copier.input(self.altset_list)
        if (self.altset_list2 is not None):
            copier.input(self.altset_list2)
        copier.input(self.covarfile)
        # similar code elsewhere
        if (self.filenull is not None):
            root, ext = os.path.splitext(self.filenull)
            # allow other formats?
            if ext.lower() == ".bed":
                copier.input(root + ".bed")
                copier.input(root + ".bim")
                copier.input(root + ".fam")
            else:
                copier.input(root + ".ped")
                copier.input(root + ".map")
        if (self.extractSim is not None):
            copier.input(self.extractSim)
        if (self.nullfitfile is not None):                        
            copier.input(self.nullfitfile)
        if self.genphen is not None and self.genphen.has_key("varBackNullFileGen"):
            copier.input(self.genphen["varBackNullFileGen"]+".fam")
            copier.input(self.genphen["varBackNullFileGen"]+".bim")
            copier.input(self.genphen["varBackNullFileGen"]+".bed")
        if self.genphen is not None and self.genphen.has_key("varBackNullPhenGen") and self.genphen["varBackNullPhenGen"] is not None:
            copier.input(self.genphen["varBackNullPhenGen"])

    #Note that the files created are not automatically copied. Instead,
    # whenever we want another file to be created, a second change must be made here so that it will be copied.
    def copyoutputs(self,copier):
        outfiletab=self.outfile
        copier.output(outfiletab)
        infofilename = utilx.appendtofilename(self.outfile,"info")
        copier.output(infofilename)
        if self.write_lrtperm and self.nperm>0:
            lrtperm_outfile = utilx.appendtofilename(self.outfile,"lrtperm")
            copier.output(lrtperm_outfile)
        if self._synthphenfile is not None:
            copier.output(self._synthphenfile)

    # =================================================================
    # END Implementation the IDistributable interface described in fastlmm.util.distributable
    # =================================================================

    def _check_entries(self, entries):
        assert entries.has_key('nullModel') and entries.has_key('altModel')
        assert entries['nullModel'].has_key('effect') and entries['nullModel'].has_key('link')
        assert entries['altModel'].has_key('effect') and entries['altModel'].has_key('link')

        assert entries['nullModel']['effect'] in set(['fixed', 'mixed'])
        assert entries['altModel']['effect'] in set(['fixed', 'mixed'])

        assert entries['nullModel']['link'] in set(['linear', 'logistic', 'erf'])
        assert entries['altModel']['link'] in set(['linear', 'logistic', 'erf'])

        assert all( [(v in set(['effect', 'link', 'penalty', 'approx'])) for v in entries['nullModel']] ),\
               'The only allow keys for nullModel are effect, link, penalty, and approx.'

        assert all( [(v in set(['effect', 'link', 'penalty', 'approx'])) for v in entries['altModel']] ),\
               'The only allow keys for altModel are effect, link, penalty, and approx.'

       

    def _check_params_after(self):        
        self.nperm = self.test.check_nperm(self.nperm)
        if self.minsetsize is not None and self.minsetsize<1: raise Exception("minsetsize must be 1 or greater")
        if self.nullfitfile is not None: assert (self.nperm is None or self.nperm==0), "if using nullfitfile then nperm should be 0"
        assert ((self.nlocalperm is None) or self.nlocalperm==0) or ((self.nperm is None) or self.nperm==0), "cannot use both nperm and nlocalperm"
        assert self.nullfitfile is None or self.nlocalperm==0, "cannot use nullfitfile with nlocalperm"

        if self.genphen is not None:
            if 'casefrac' not in self.genphen.keys(): self.genphen['casefrac']=None
            if 'fracCausal' not in self.genphen.keys(): self.genphen['fracCausal']=1.0
            if 'numBackSnps' not in self.genphen.keys(): self.genphen['numBackSnps']=0
            if 'varBackNullFileGen' not in self.genphen.keys(): self.genphen['varBackNullFileGen']=None  
    
    def _check_params_before(self):
        pass
        #if self.test== "lrt_up" and self.filenull is None:
        #    self.test = "lrt"
                   
    def pv_from_localperm(self, iperm, result, t0, varcomp_test, y):
        '''
        This is not test specific, but is currently only implemented for LRT. No code would 
        need to be changed here for other tests, only the test-specific files elsewhere.
        '''
        assert iperm==-1, "should not be permuting for global null fit here, only local"
        assert str(self.test)=="lrt", "can only use local permutations with lrt, not with " + str(self.test)            
        assert self.altModel["effect"]=="mixed" and self.altModel["link"]=="linear", "nlocalperm only set up for linear lrt at the moment"
        assert self.nullModel["effect"]=="fixed", "nlocalperm only set up for one-kernel" 
        #assert self.__G0 is None, "only 1-kernel implemented"
        assert self.__SNPs0 is None, "only 1-kernel implemented"
        logging.info("running up to " + str(self.nlocalperm) + " local permutations")
        permstatbetter=[]
        allstat=[]
        alteqnullperm=[]
        yorig=y
        tol=0.0
        from numpy.random import RandomState
        randomstate = RandomState(self.rseed)
        checkpoint=1000
        powerneeded=2
        for pm in xrange(self.nlocalperm):
            if pm==checkpoint:
                logging.info('checkpointing '  + str(pm))                    
                numbetter=sp.sum(permstatbetter)
                pv=numbetter/float(pm)                                
                t1=time.time()
                logging.info("Checkpoint time elapsed=%.2f seconds" % (t1-t0))
                if pv*checkpoint>=10**powerneeded:
                    logging.info("stopping due to checkpoint")
                    break;
                checkpoint=checkpoint*2
            
            permutationIndex = utilx.generatePermutation(y.shape[0],randomstate)
            yperm=y[permutationIndex]           
            Xperm=self.__X[permutationIndex]                
                                                                   
            permresult = varcomp_test.testGupdate(yperm,Xperm,self.test)  
            allstat.append(permresult['stat'])
            alteqnullperm.append(permresult['alteqnull'])
            permstatbetter.append(permresult['stat']-result.test['stat']>tol)  
                       
        numbetter=sp.maximum(0.99,sp.sum(permstatbetter))
        pv=numbetter/float(pm) 
        result.test['pv-local']=pv
        result.test['pv']=result.test['pv']
        
        if self.fitlocal and pv*pm<10**powerneeded:
            import fastlmm.association.tests.Cv as cv
            logging.info("fitting aUD to local with " + str(pm+1) + " permutations")
            pv_adj,mixture,scale,dof=cv.lrtpvals_qqfit(nperm=self.nlocalperm, lrt=sp.array([result.test['stat']]), lrtperm=sp.array(allstat), alteqnull=sp.array([result.test['alteqnull']]), alteqnullperm=sp.array(alteqnullperm), qmax=self.qmax)
            result.test['pv-local-aUD']=pv_adj[0]
            logging.info("mixture (non-zero dof)="+str(mixture))
            logging.info("dof="+str(dof))
            logging.info("scale="+str(scale))
        else:
            result.test['pv-local-aUD']=sp.NaN
        
        logging.info("    used " + str(pm+1) + " permutations to compute p=" + str(pv) + ", p50=" + str(result.test['pv']))
        
    def G_exclude(self, i_exclude):
        if self.__SNPs0.has_key("data"):
            G_exclude = self.__SNPs0["data"]["snps"][:,i_exclude]
        else:
            snp_names = self.__SNPs0["reader"].rs[i_exclude]
            snp_set = SnpAndSetName('G_exclude', snp_names)
            G_exclude = self.__SNPs0["reader"].read(snp_set)['snps']
            G_exclude = util.standardize(G_exclude)
            #normalize
            pass
        G_exclude/=sp.sqrt(self.__SNPs0["num_snps"])
        return G_exclude

    def run_test(self, SNPs1, G1, y, altset, iset, ichrm, iposrange, iperm = -1, varcomp_test=None):
        '''
        This function does the main work of the class, and also reads in the SNPs for the alternative model.
        It is called (via a lambda) inside the loops found in 'generate_sequence'.

        Input:
            altset - a set of snps
            iset - index to altset
            iperm - index to permutation (-1 means no permutation)
            varcomp -if not None, assume that it is the correct one for this test, and does not re-compute anything
        Output:
            a list (often just one) of
                instances of the Result class, varcomp (for caching)
        '''
        t0=time.time()            

        self.run_once() #load files, etc. -- stuff we only want to do once no matter how many times we call 'run_test' and then 'reduce.        

        logging.info("Working on permutation index #{0} (where -1 means no permutation) out of {1} total permutations...".format(iperm, self.nperm))
        logging.info("\taltset {0}, {1} of {2}  ".format(altset, iset, len(self.altsetlist_filtbysnps)))
                      
        result = Result(iperm=iperm,iset=iset,setname=str(altset), ichrm=ichrm, iposrange=iposrange)
         
        if varcomp_test is None or not self.cache_from_perm:            
            if G1.shape[1]==0:
                logging.info( "no SNPS in set " + setname )
                result=None
                return [result]
            if sp.isnan(G1.sum()): raise Exception("found missing values in test SNPs that remain after intersection for " + str(altset))

            result.setsize = SNPs1['snps'].shape[1]
                
            logging.info(" (" + str(result.setsize) + " SNPs)")              

            if self.filenull is not None:
                if self.__SNPs0.has_key('data'):
                    pos0=self.__SNPs0['data']['pos']
                else:
                    pos0=self.__SNPs0['snp_set'].pos
                i_exclude =  utilx.excludeinds(pos0, SNPs1['pos'], mindist = self.mindist,idist = self.idist)
                result.nexclude = i_exclude.sum()
                logging.info("numExcluded=" + str(result.nexclude) + "")
            else:
                i_exclude = None
                result.nexclude=0;


        if self.permute is not None: #permute the data to test for type I error
                permutationIndex = utilx.generatePermutation(SNPs1['snps'].shape[0],self.mainseed ^ self.permute)      
        if iperm >= 0 :  #permute the data to create Null-only P-values for lrt null fitting
            eachseed=utilx.combineseeds(self.calseed,iperm)            
            newseed = self.mainseed ^ eachseed  

            if self.permute is not None :
                #Add a left shift so that iperm=1 and self.permute=2 gives a different newseed than 
                #iperm=2 and self.permute=1
                newseed = (newseed << 1) ^ self.permute
            permutationIndex = utilx.generatePermutation(SNPs1['snps'].shape[0], newseed)
        if varcomp_test is None or not self.cache_from_perm:
            #need to make this caching smarter for when background kernel changes in every test (e.g. interactions)
            if (self.permute is not None) or (iperm >=0):
                G1=G1[permutationIndex]

            #if result.nexclude:
            #    logging.info(" (computing needed null info anew) ")
            #    G0_to_use = self.__SNPs0#['data']['snps'][:,~i_exclude]/SP.sqrt(self.__SNPs0['data']['snps'][:,~i_exclude].shape[1])   
            #else: 
            #    G0_to_use=self.__G0
        
            if (str(self.test)!="lrt_up") and (self.genphen is not None) and (not self.genphen['once']) or ((str(self.test)!="lrt_up") and result.nexclude>0) or self.__varcomp_test is None:
                #need to re-construct the test (can't cache)
                #note that if result.nexclude kicked in, then this would happen in clause above anyhow                               
                varcomp_test = self.varcomp_test_setup(y,self.__SNPs0,i_exclude)
            else: #use cached values                
                #G0_to_use=self.__G0
                varcomp_test=self.__varcomp_test
                #varcomp_test.i_exclude=i_exclude
                
                logging.info(" (using cached null info)")
            
            #create G_exclude                         
            if (str(self.test)=="lrt_up" and result.nexclude>0):
                G_exclude = self.G_exclude(i_exclude)
            else:
                G_exclude = None            
         
            result.test = varcomp_test.testG(G1,self.test, i_exclude=i_exclude,G_exclude=G_exclude)             
            logging.info("p=%.2e",result.test['pv'])                   

            # do the permutations here, rather than as they would normally be done with a seperate call to             
            if self.nlocalperm is not None and self.nlocalperm>0:
                self.pv_from_localperm(iperm, result, t0, varcomp_test, y)
        else: #caching for 1-kernel LRT (save expensive stuff from non-permuted run)
            #currently not used, because need to be able to pass back values to work_sequence for caching
            assert iperm!=-1, "should be a permutation run"
            assert str(self.test)=="lrt", "can only use lrt, here now"
            assert self.altModel["effect"]=="mixed" and self.altModel["link"]=="linear", "caching only set up for linear lrt at the moment"
            assert self.nullModel["effect"]=="fixed", "caching only set up for one-kernel" 
            assert self.__G0 is None, "caching only 1-kernel implemented, otherwise need to permute G0 as well here"

            #permute X, y instead of G1 so can use caching                  
            yperm=y[permutationIndex]                
            Xperm=self.__X[permutationIndex]       
            result = varcomp_test.testGupdate(yperm,Xperm,self.test)  
                                                  
        t1=time.time()
        logging.info("%.2f seconds elapsed in run_test" % (t1-t0))
                
        return [result]

    def TESTBEFOREUSING_run_interactiontest(self,altset,iset,altset2,iset2,iperm = -1):
        '''
        This function does the main work of the class (when interaction tests are running).
        It is called (via a lambda) inside the loops found in 'generate_sequence'

        Input:
            altset - a set of snps
            iset - index to altset
            iperm - index to permutation (-1 means no permutation)
        Output:
            an instance of the Result class
        '''

        t0=time.time()

        self.run_once() #load files, etc. -- stuff we only want to do once no matter how many times we call 'run_test' and then 'reduce.

        logging.info("Working on permutation index #{0} out of {1} total permutations...".format(iperm, self.nperm))
        logging.info("\taltset {0}, {1} of {2}  ".format(altset, iset, len(self.altsetlist_filtbysnps)))
        self.stdout.flush()

        result = PairResult(iperm=iperm,iset=iset,setname=str(altset),iset2=iset2,setname2=str(altset2))

        #read the alternative model SNPs (for all sets) from disk.
        SNPs1a = altset.read()      #first set
        import fastlmm.util.preprocess as util
        SNPs1a['snps'] = util.standardize(SNPs1a['snps'])

        SNPs1b = altset2.read()     #second set
        SNPs1b['snps'] = util.standardize(SNPs1b['snps'])

        SNPs1=dict(SNPs1a)
        SNPs1.update(SNPs1b) #concatenate the SNP dictionaries
  
        G1 = SNPs1['snps']/sp.sqrt(SNPs1['snps'].shape[1])
        if G1.shape[1]==0: raise Exception("no SNPs to test")
        result.setsize = SNPs1['snps'].shape[1]

        logging.info(" (" + str(result.setsize) + " SNPs)")

        if self.permute >= 0 :
            permutationIndex = utilx.generatePermutation(SNPs1['snps'].shape[0],self.rseed ^ self.permute)
            G1=G1[permutationIndex]

        if iperm >= 0 :
            newseed = self.rseed ^ iperm
            if self.permute >= 0 :
                #Add a left shift so that iperm=1 and self.permute=2 gives a different newseed than iperm=2 and self.permute=1
                newseed = (newseed << 1) ^ self.permute
            permutationIndex = utilx.generatePermutation(SNPs1['snps'].shape[0], newseed)
            #permute the data to create Null-only P-values
            G1=G1[permutationIndex]

        if self.filenull is not None:
            i_exclude =  utilx.excludeinds(self.__SNPs0['data']['pos'], SNPs1['pos'], mindist = self.mindist,idist = self.idist)
            result.nexclude = i_exclude.sum()
            logging.info("numExcluded=" + str(result.nexclude))
        else:
            result.nexclude=0;

        null_model_changed=(result.nexclude>0)

        #need to make this caching smarter for when background kernel changes in every test (e.g. interactions)
        if null_model_changed:
            logging.info(" (computing needed null info anew) ")
            G0_to_use = self.__SNPs0['data']['snps'][:,~i_exclude]/SP.sqrt(self.__SNPs0['data']['snps'][:,~i_exclude].shape[1])  #excluded SNPs
            null_model=None
        else: #use cached values
            G0_to_use=self.__G0
            null_model=self.__varcomp_test           
            logging.info(" (using cached null info)")

        [result.pv,result.lik0,result.lik1] = self.test.pv_etc(filenull, G0_to_use, G1, y, x, appendBias, null_model, self.__varcomp_test, forcefullrank)

        # where is result.test used?
        result.test=self.test #probably not the best way to do this, but hacking it for now

        t1=time.time()
        logging.info("%.2f seconds elapsed" % (t1-t0))

        return result
  
    def setSNPs0(self):        
        logging.info("Reading SNPs0")
        root, ext = os.path.splitext(self.filenull)
        if (self.autoselect):
            #raise Exception("this autoselect option has been removed. please use the stand-alone FastLmm-SELECT in this same download")
            assert not self.forcefullrank, "cannot currently force full rank with autoselect=True"
            if self.extractSim is not None : raise Exception("extractSim not supported with autoselect")
            #call C++ version of FaST-LMM select to determine background Covariance matrix of the LMM
            osname = sys.platform
            dir = os.path.split(__file__)[0] #__file__ is the pathname of the file from which the module
            if (osname.find("win") >= 0):    #         was loaded, if it was loaded from a file
                fastlmmpath = os.path.join(dir,"Fastlmm_autoselect", "fastlmmc.exe")
            elif (osname.find("linux") >= 0):
                fastlmmpath = os.path.join(dir,"Fastlmm_autoselect", "fastlmmc")
                import stat
                st = os.stat(fastlmmpath)
                os.chmod(fastlmmpath, st.st_mode | stat.S_IEXEC)
            else:
                logging.info('\n\n unsupported operating system!')
            if not os.path.isfile(fastlmmpath) : raise Exception("Expect file {0}".format(fastlmmpath))
            logging.info("fastlmmC path=" + fastlmmpath)

            if (not os.path.exists('select')):
                subprocess.call("mkdir select", shell=True)

            with TemporaryFile() as stdout_temp:
                if (self.covarfile == None):
                    subprocess.call(fastlmmpath + ' -autoselect select/auto -bfilesim ' + root + ' -pheno ' + self.phenofile + ' -pedOutFile ' + "select/null", shell=True, stdout=stdout_temp, stderr=subprocess.STDOUT)
                else:
                    subprocess.call(fastlmmpath + ' -autoselect select/auto -bfilesim ' + root + ' -pheno ' + self.phenofile + ' -covar ' + self.covarfile + ' -pedOutFile ' + "select/null", shell=True, stdout=stdout_temp, stderr=subprocess.STDOUT)
                stdout_temp.seek(0)
                stdout_string = stdout_temp.read()
            logging.info(stdout_string)

            self.__SNPs0 = {
                "data":readPED("select/null"),
                "filename":"select/null",                
                }
            self.__SNPs0["original_iids"] = self.__SNPs0["data"]["iid"]
            self.__SNPs0["num_snps"] = self.__SNPs0["data"]["snps"].shape[1]
        else:                   
            if ext.lower() == ".bed":
                snp_set = CreateSnpSetReaderForFileNull(self.extractSim)  
                bed = Bed(root)
                snp_set_bed = snp_set.addbed(bed)
                num_snps = len(snp_set_bed)                
                self.__SNPs0={
                    "snp_set":snp_set_bed,
                    "reader":snp_set_bed.bed,
                    "filename":root,
                    #"data":readBED(root, snp_set=snp_set)
                    }
                self.__SNPs0["original_iids"] = self.__SNPs0["reader"].original_iids
                self.__SNPs0["num_snps"] = num_snps
                #import fastlmm.util.preprocess as util
                #self.__SNPs0['data']['snps'] = util.standardize(self.__SNPs0['data']['snps'])
            else:
                if self.extractSim is not None : raise Exception("extractSim not supported with filenull in ped format, please use bed format")
                self.__SNPs0 = {"data":readPED(root),
                                "filename":root}
                self.__SNPs0["original_iids"] = self.__SNPs0["data"]["iid"]
                self.__SNPs0["num_snps"] = self.__SNPs0["data"]["snps"].shape[1]

    def varcomp_test_setup(self,y,SNPs0=None, i_exclude=None): 
        #in case we exclude all SNPs (needed for score, not for lrt, not sure about lrt_up)
        nullModelTmp=self.nullModel.copy()           
        if (i_exclude is not None) and SNPs0['num_snps']==i_exclude.sum():
            excluded_all_snps = True           
            nullModelTmp['effect']='fixed'
        else:
            excluded_all_snps = False        

        if self.filenull is not None and (not excluded_all_snps): #aka two-kernel after exclusion has been done        
            varcomp_test = self.test.construct(y,self.__X, forcefullrank = self.forcefullrank,
                             SNPs0 = SNPs0, i_exclude = i_exclude, nullModel = self.nullModel, altModel = self.altModel,
                             scoring = self.scoring, greater_is_better = self.greater_is_better)
        else: #no background kernel
            varcomp_test = self.test.construct_no_backgound_kernel(y,self.__X,
                             forcefullrank = self.forcefullrank, nullModel = nullModelTmp,
                             altModel = self.altModel, scoring = self.scoring,
                             greater_is_better = self.greater_is_better)        
        return varcomp_test

    #so that it only ever gets changed once, in run_once
    @property
    def rseed(self):
        assert self._seed is not None
        return self._seed

    @property
    def mainseed(self):
        '''used to xor indexes with to get seeds for other purposes'''
        return 34343
       
    def run_once(self):
        '''
        Loading datasets (but not SNPs for alternative), etc. that we want to do only once no mater how many times we call 'run_test' and then perhaps 'reduce'.        '''         
               
        if self.calseed is not None:            
            self._seed=utilx.combineseeds(self.calseed,self.mainseed)                        
            #self._seed = self.calseed ^ self.mainseed
        else:
            self._seed = self.mainseed

        if (self._ran_once):
            return
        #self.__snapshot = self.__repr__() # must be before seeting _ran_once
        self._ran_once = True

        logging.info("Running Once")

        self._check_params_before() #must be called before test object Lrt(), Cv(), Sc() is created       
        #!! make this nicer
        if self.test=="lrt":
            self.test = Lrt()
        elif self.test=="lrt_up":
            self.test = lrt.LRT_up()
        elif self.test=="cv":
            self.test = Cv()
        elif self.test=="sc_davies":
            self.test = Sc("davies")
        elif self.test=="sc_mom":
            self.test = Sc("mom")
        else:
            raise Exception("Don't know how to construct {0}".format(self.test))        
        self._check_params_after() #must be called after test object Lrt(), Cv(), Sc() is created         

        phenhead,phentail=os.path.split(self.phenofile)
        logging.info("Reading Pheno: " +  phentail + "")        
        pheno = loadPhen(filename = self.phenofile)     
        pheno['vals']=pheno['vals'][:,self.mpheno-1]#use -1 so compatible with C++ version
        #pheno['header']=pheno['header'][self.mpheno-1] #header is empty 
        goodind=sp.logical_not(sp.isnan(pheno['vals']))
        pheno['vals']=pheno['vals'][goodind]
        pheno['iid']=pheno['iid'][goodind,:]        
           
        if self.genphen is not None:
            #for synthetic data generation--this is the background signal added to the phenotype
            if self.genphen["varBackNullFileGen"] is not None:
                self.__varBackNullSnpsGen=readBED(self.genphen["varBackNullFileGen"])
                import fastlmm.util.preprocess as util
                self.__varBackNullSnpsGen['snps'] = util.standardize(self.__varBackNullSnpsGen['snps'])
            else:
                self.__varBackNullSnpsGen=None

        if self.filenull is not None:
            self.setSNPs0()
            logging.info("Filter SNPs0")
            #prelim intersection so don't read in all null SNP iids if way more than in pheno; now redundant
            #filter(pheno, self.__SNPs0['data'])
        else:
            if self.extractSim is not None : raise Exception("extractSim should be specified only if filenull is specified.")
            self.__SNPs0=None
            #self.__G0=None

        if self.covarfile==None:
            covar=None
        else:
            covar=loadPhen(self.covarfile)
            if self.covarimp=='standardize':
                covar['vals'],fracmissing=utilx.standardize_col(covar['vals'])
            elif self.covarimp is None:
                pass
            else:
                raise Exception("covarimp=" + self.covarimp + " not implemented")          

        covar,pheno,indarr=self.intersect_data(covar, pheno)
        N=pheno['vals'].shape[0];
        if self.covarfile==None:
            self.__X = sp.ones((N,1))
        else:
            #check for covariates which are constant, as the lmm code crashes on these            
            badind = utilx.indof_constfeatures(covar['vals'], axis=0)
            if len(badind)>0:
                raise Exception("found constant covariates with indexes: %s: please remove, or modify code here to do so" % badind)
            self.__X=sp.hstack((sp.ones((N,1)),covar['vals']))
        self.__y = pheno['vals']

        if not covar is None and sp.isnan(covar['vals'].sum()): raise Exception("found missing values in covariates file that remain after intersection")
        if sp.isnan(self.__y.sum()): raise Exception("found missing values in phenotype file that remain after intersection")
        #if hasattr(self,'__G0') and not self.__G0 is None and sp.isnan(self.__G0.sum())): raise Exception("found missing values in background SNPs that remain after intersection")

        #creating sets from set defn files. Looks at bed file to filter down the SNPs to only those present in the bed file
        self.altset_list,self.altsetlist_filtbysnps = self.create_altsetlist_filtbysnps(self.altset_list,self.alt_snpreader)
        if self.altset_list2 is not None:
                self.altset_list2,self.__altsetlist2_filtbysnps = self.create_altsetlist_filtbysnps(self.altset_list2,self.alt_snpreader)

        if self.npermabs is not None and self.nperm is not None:
            assert self.nperm is None, "cannot use both nperm and npermabs"
        if self.npermabs is not None:
            num_tests = len(self.altsetlist_filtbysnps)
            self.nperm = int(sp.ceil(self.npermabs/num_tests))
            logging.info("using npermabs=%i, found %i tests, so using %i permutations per test" % (self.npermabs, num_tests, self.nperm))            
            self.npermabs = None
        
        logging.info("------------------------------------------------")
        logging.info("Found " + str(self.__X.shape[0]) +  " individuals")
        logging.info("Found " + str(self.__X.shape[1]) +  " covariates (including bias)")
        if self.__SNPs0 is not None:
            logging.info("Found " + str(self.__SNPs0['num_snps']) +  " SNPs for null kernel")
        if (len(self.__y.shape)>1): nPhen=self.__y.shape[1]
        else: nPhen=1
        logging.info("Found " + str(nPhen) +  " phenotypes")
        if self.filenull != None:
            logging.info("Running test " +  str(self.test) + " with 1 background kernel")
        else:
            logging.info("Running test " +  str(self.test)  + " with no background kernel")
        #CL3_28_14:The following never evaluates True, as it should be #if str(self.test)=="lrt":
        #but the body is buggy, so it would crash
        #if self.test=="lrt":
        #    logging.info("nullmodel: " + "link=" + self.nullModel['link'] + ", " 
        #                          + "approx=" + self.nullModel['approx'] + ", "
        #                          + "pen="+ self.nullModel['penalty'])
        #    logging.info("altmodel: " + "link=" + self.altModel['link'] + ", " 
        #                          + "approx=" + self.altModel['approx'] + ", "
        #                          + "pen="+ self.altModel['penalty'])        
        logging.info("Found " + str(self.work_count/(self.nperm+1)) + " tests to run (" + str(self.work_count) + " work items, including " + str(self.nperm) + " permutations)")
        if self.nullfitfile is not None:
            logging.info("nullfitfile=" + self.nullfitfile)
        if self.altset_list2 is not None:
            logging.info("(running paired set tests)")
        else:
            logging.info("(running single set tests)")        
        if self.minsetsize is not None and self.maxsetsize is not None:
            logging.info("min set size=%i, max set size=%i", self.minsetsize, self.maxsetsize) 
        if self.mpheno is not None: logging.info("mpheno=%i", self.mpheno)
        if self.extractSim is not None: logging.info("extractSim=%s",self.extractSim)
        logging.info("------------------------------------------------")

        logging.info("Creating null model... ")
        #cache whatever is needed for each case (sometimes likelihood, sometimes rotated items, sometimes nothing)
        t0=time.time()
        
        #if self.filenull is not None: #aka two-kernel
        #    self.__G0 = self.__SNPs0['data']['snps']/sp.sqrt(self.__SNPs0['data']['snps'].shape[1])        
        #else:
        #    self.__G0=None
        
        #for some models, this caches important and expensive information
        #as appropriate, gets re-computed in run_test()        
        if self.genphen is None: #or self.genphen['once']: THIS CAUSES REAL DATA TO GET USED         
            #self.__varcomp_test=self.varcomp_test_setup(self.__y,self.__G0)            
            self.__varcomp_test=self.varcomp_test_setup(self.__y,self.__SNPs0)
        #elif self.genphen is not None and self.__G0 is not None:
        elif self.genphen is not None and self.__SNPs0 is not None:
            raise Exception("this will run if you comment out this line, but I think will be inefficient because it's not caching null info across each gene-set specific phenotypes even when G0 remains constant across these")
            self.__varcomp_test=None
        else:
            # single kernel synthetic, fine to re-do null model computations each time
            self.__varcomp_test=None

        t1=time.time()
        logging.info("done. %.2f seconds elapsed" % (t1-t0))

    def check_id_order(self, covar, pheno):        
        ids_pheno=pheno['iid']
        ids_SNPs=self.alt_snpreader.original_iids
        if (len(ids_pheno)!=len(ids_SNPs)) or (not sp.all(ids_pheno==ids_SNPs)):
          return False
        if self.covarfile is not None:
            ids_X=covar['iid']
            if (len(ids_pheno)!=len(ids_X)) or (not sp.all(ids_pheno==ids_X)): return False
        if self.filenull is not None:
            ids_nullSNPs=self.__SNPs0['original_iids']
            if not sp.all(ids_pheno==ids_nullSNPs): return False
        if hasattr(self,"__varBackNullSnpsGen") and self.__varBackNullSnpsGen is not None: #for synthetic data generation            
            ids_SNPgen=self.__varBackNullSnpsGen["iid"]
            if not sp.all(ids_pheno==ids_SNPgen): return False
        return True

    def intersect_data(self, covar, pheno):
        if self.check_id_order(covar,pheno):
            indarr = SP.arange(pheno['iid'].shape[0])
            logging.info(str(indarr.shape[0]) + " IDs match up across data sets")
        else:            
            logging.info("IDs do not match up, so intersecting the data over individuals")            
            if self.filenull is not None:
                nullsnpids=self.__SNPs0['original_iids']
            else:
                nullsnpids=None
            if self.covarfile is not None:
                covarids=covar['iid']
            else:
                covarids=None
            if hasattr(self,"__varBackNullSnpsGen") and self.__varBackNullSnpsGen is not None:
                varbacknullsnp=self.__varBackNullSnpsGen['iid']
            else:
                varbacknullsnp=None

            #the order of inputs here is reflected in the reordering indexes below of 0,1,2,3,4
            import warnings
            #!!warnings.warn("This intersect_ids is deprecated. Pysnptools includes newer versions of intersect_ids", DeprecationWarning)
            indarr=utilx.intersect_ids([pheno['iid'],self.alt_snpreader.original_iids,covarids,nullsnpids,varbacknullsnp])            
            assert indarr.shape[0]>0, "no individuals remain after intersection, check that ids match in files"
            #sort the indexes so that SNPs ids in their original order (and 
            #therefore we have to move things around in memory the least amount)
            #[indarr[:,3] holds the SNP order    
            sortind=sp.argsort(indarr[:,3])
            indarr=indarr[sortind]

            pheno['iid']=pheno['iid'][indarr[:,0]]
            pheno['vals']=pheno['vals'][indarr[:,0]]

            self.alt_snpreader.ind_used=indarr[:,1]

            if self.covarfile is not None:
                covar['iid']=covar['iid'][indarr[:,2]]
                covar['vals']=covar['vals'][indarr[:,2]]

            if self.filenull is not None:  
                #copy in memory to make later access faster (presumably?):
                #self.__SNPs0['iid']=sp.copy(self.__SNPs0['iid'][indarr[:,3]])
                #self.__SNPs0['snps']=sp.copy(self.__SNPs0['snps'][indarr[:,3]])
                #self.__SNPs0['data']['iid']=self.__SNPs0['data']['iid'][indarr[:,3]]
                if self.__SNPs0.has_key('data'):
                    self.__SNPs0['data']['iid'] = self.__SNPs0['data']['iid'][indarr[:,3]]
                    self.__SNPs0['data']['snps'] = self.__SNPs0['data']['snps'][indarr[:,3]]
                else:
                    self.__SNPs0['reader'].ind_used = indarr[:,3]
                #self.__SNPs0['data']['snps']=self.__SNPs0['data']['snps'][indarr[:,3]]

            if hasattr(self,"_varBackNullSnpsGen") and self.__varBackNullSnpsGen is not None:
                self.__varBackNullSnpsGen['iid']=self.__varBackNullSnpsGen['iid'][indarr[:,4]]
                self.__varBackNullSnpsGen['snps']=self.__varBackNullSnpsGen['snps'][indarr[:,4]]

            logging.info(str(indarr.shape[0]) + " ids left")
        
        return covar,pheno,indarr


    def create_altsetlist_filtbysnps(self, altset_list, altsetbed):
        '''
        returns: altset_list, altsetlist_filtbysnps
        'altset_list' is the raw set defintion read in from file
        'altsetlist_filtbysnps' is the result of filtering altset_list by SNPsin the specified altsetbed SNP file
        Additionally, only allows sets specified by sets are allowed (if this is not None)
        '''
        if isinstance(altset_list, str):
            altset_list = SnpAndSetNameCollection(altset_list)

        if self.sets is not None:  #filter by those sets specified on the command line, stored in self.sets
            altset_list=Subset(altset_list, self.sets)

        if self.minsetsize is not None or self.maxsetsize is not None:
            altset_list=MinMaxSetSize(altset_list,minsetsize=self.minsetsize,maxsetsize=self.maxsetsize) #would be nice if applying 1,None wouldn't wrap at all (becuase that is the identity)

        altsetlist_filtbysnps = altset_list.addbed(altsetbed)

        if len(altsetlist_filtbysnps)==0: raise Exception("Expect altset_list to contain at least one set")

        return altset_list, altsetlist_filtbysnps

    # overridden to provide a nice name on the cluster.
    def __str__(self):
        if self.outfile == None:
            return self.__class__.__name__
        else:
            return "{0} {1}".format(self.__class__.__name__, self.outfile)

    def __repr__(self):
        import cStringIO
        fp = cStringIO.StringIO()
        fp.write("{0}(\n".format(self.__class__.__name__))
        varlist = []
        for f in dir(self):
            if f.startswith("_"): # remove items that start with '_'
                continue
            if type(self.__class__.__dict__.get(f,None)) is property: # remove @properties
                continue
            if callable(getattr(self, f)): # remove methods
                continue
            varlist.append(f)
        for var in varlist[:-1]: #all but last
            fp.write("\t{0}={1},\n".format(var, getattr(self, var).__repr__()))
        var = varlist[-1] # last
        fp.write("\t{0}={1})\n".format(var, getattr(self, var).__repr__()))
        result = fp.getvalue()
        fp.close()
        return result

def CreateSnpSetReaderForFileNull(snp_set):
    if (snp_set is None):
        return AllSnps()
    if isinstance(snp_set, str):        
        return SnpsFromFile(snp_set)
    return snp_set



if __name__ == "__main__":
    logging.basicConfig(level=logging.INFO)
    logging.info("last modified: %s" % time.ctime(os.path.getmtime(__file__)))
    logging.info("-------------------------------------------------")       
    
    # Parse the command line
    if len(sys.argv) == 2:
        sys.argv.append("Local()")
    if len(sys.argv) != 3:
        logging.info("Usage:python FastLmmSet.py <inputfile> [optional <runner>]")
        sys.exit(1)
    file,infile,runner_string = sys.argv
    if not os.path.exists(infile): raise Exception("Fatal Error: " + infile + " does not exist")

    #Run the job specified in infile    
    from fastlmm.association.FastLmmSet import *
    exec("runner = " + runner_string)
    d = open(infile).read()    
    exec(d)     
          
    tt0=time.time()    
    runner.run(distributable)    
    tt1=time.time()
    logging.info("Final elapsed time for all processing is %.2f seconds" % (tt1-tt0))

