def getPL(self,r,RSSStd):
        """ Get Power Level from a given distance

        Parameters
        ----------

        r : range
        RSSStd : range standard deviation

        Examples
        --------

        >>> M = PLSmodel(f=0.3,rssnp=2.64,d0=1,sigrss=3,method='mode')
        >>> PL =  M.getPL(16,1)

        """

        if self.method =='OneSlope':
            PL=self.OneSlope(r)

        elif self.method == 'mode' or self.method == 'median' or self.method == 'mean':
            PLmean          = self.getPLmean(r)
            try:
                shPLmean        = np.shape(PLmean)
                Xrand           = RSSStd*sp.randn(shPLmean[0])
            except:
                Xrand           = RSSStd*sp.randn()
            PL        = PLmean+Xrand

        else :
            raise NameError('Pathloss method name')

        return(PL)

 def _generate_null_shift_errors(self):
     """
     Returns null shift (constant bias) entries in a 6 entry array:
        array([gx_n, gy_n, gz_n, ax_n, ay_n, az_n])
     
     where the subscript 'n' stands for 'null shift'.
     
     Note
     -----
     The constant bias is constant during the run, but varies from run-to-run.  
     
     If, for experimental purposes, you'ld like the SAME constant bias generated
     every run, then this function should be modified to use the standard deviation
     value (not multiplied by random number).
     """
     # If the null-shift value has been generated once, then that value should be used.
     # Otherwise, it will be generated and saved for the next call.        
     try:
         return self._constant_bias
     except AttributeError:          
         # Original code used a non-random null-shift.  This could be acceptable if
         # using a single vehicle, but unrealistic if used for a entire community.  So now the null-shift is random.
         accel_null = self._sqd['sigma_n_f'] * sp.randn(3)
         gyro_null  = self._sqd['sigma_n_g'] * sp.randn(3)
         
         self._constant_bias =  np.hstack((gyro_null, accel_null)) 
         return self._constant_bias

def fitPairwiseModel(Y,XX=None,S_XX=None,U_XX=None,verbose=False):
    N,P = Y.shape
    """ initilizes parameters """
    RV = fitSingleTraitModel(Y,XX=XX,S_XX=S_XX,U_XX=U_XX,verbose=verbose)
    Cg = covariance.freeform(2)
    Cn = covariance.freeform(2)
    gp = gp2kronSum(mean(Y[:,0:2]),Cg,Cn,XX=XX,S_XX=S_XX,U_XX=U_XX)
    conv2 = SP.ones((P,P),dtype=bool)
    rho_g = SP.ones((P,P))
    rho_n = SP.ones((P,P))
    for p1 in range(P):
        for p2 in range(p1):
            if verbose:
                print '.. fitting correlation (%d,%d)'%(p1,p2)
            gp.setY(Y[:,[p1,p2]])
            Cg_params0 = SP.array([SP.sqrt(RV['varST'][p1,0]),1e-6*SP.randn(),SP.sqrt(RV['varST'][p2,0])])
            Cn_params0 = SP.array([SP.sqrt(RV['varST'][p1,1]),1e-6*SP.randn(),SP.sqrt(RV['varST'][p2,1])])
            params0 = {'Cg':Cg_params0,'Cn':Cn_params0}
            conv2[p1,p2],info = OPT.opt_hyper(gp,params0,factr=1e3)
            rho_g[p1,p2] = Cg.K()[0,1]/SP.sqrt(Cg.K().diagonal().prod())
            rho_n[p1,p2] = Cn.K()[0,1]/SP.sqrt(Cn.K().diagonal().prod())
            conv2[p2,p1] = conv2[p1,p2]; rho_g[p2,p1] = rho_g[p1,p2]; rho_n[p2,p1] = rho_n[p1,p2]
    RV['Cg0'] = rho_g*SP.dot(SP.sqrt(RV['varST'][:,0:1]),SP.sqrt(RV['varST'][:,0:1].T))
    RV['Cn0'] = rho_n*SP.dot(SP.sqrt(RV['varST'][:,1:2]),SP.sqrt(RV['varST'][:,1:2].T))
    RV['conv2'] = conv2
    #3. regularizes covariance matrices
    offset_g = abs(SP.minimum(LA.eigh(RV['Cg0'])[0].min(),0))+1e-4
    offset_n = abs(SP.minimum(LA.eigh(RV['Cn0'])[0].min(),0))+1e-4
    RV['Cg0_reg'] = RV['Cg0']+offset_g*SP.eye(P)
    RV['Cn0_reg'] = RV['Cn0']+offset_n*SP.eye(P)
    RV['params0_Cg']=LA.cholesky(RV['Cg0_reg'])[SP.tril_indices(P)]
    RV['params0_Cn']=LA.cholesky(RV['Cn0_reg'])[SP.tril_indices(P)]
    return RV

    def setUp(self):
        np.random.seed(1)

        # generate data
        N = 400
        s_x = 0.05
        s_y = 0.1
        X = (sp.linspace(0, 2, N) + s_x * sp.randn(N))[:, sp.newaxis]
        Y = sp.sin(X) + s_y * sp.randn(N, 1)
        Y -= Y.mean(0)
        Y /= Y.std(0)

        Xstar = sp.linspace(0, 2, 1000)[:, sp.newaxis]

        # define mean term
        F = 1.0 * (sp.rand(N, 2) < 0.2)
        mean = lin_mean(Y, F)

        # define covariance matrices
        covar1 = SQExpCov(X, Xstar=Xstar)
        covar2 = FixedCov(sp.eye(N))
        covar = SumCov(covar1, covar2)

        # define gp
        self._gp = GP(covar=covar, mean=mean)

def fit_krr_dskl_rks_result(X,Y,Xtest,Ytest,its=100,eta=.1,C=.001,nPredSamples=30,nExpandSamples=10, kernel=(GaussianKernel,(1.))):
    # random gaussian for rks
    Zrks = sp.randn(len(Y),X.shape[0]) / (kernel[1]**2)
    Wrks = sp.randn(len(Y))
    for it in range(1,its+1):
        Wrks = step_dskl_rks_krr(X,Y,Wrks,Zrks,eta/it,C,nPredSamples,nExpandSamples)
    return predict_krr_rks(Xtest,Wrks,Zrks) 

def test_converting_to_factors():

    test_data = DataFrame(
        {
            'colA': Series(randn(1, 5000).flatten() > 0),
            'colB': Series(100 * randn(1, 5000).flatten()),
            'colC': Series(100 + randn(1, 5000).flatten()),
            'colD': Series(randn(1, 5000).flatten() > 0),
        },
    )

    test_data['colA'] = test_data['colA'].map(str)
    test_data['colD'] = test_data['colD'].map(str)

    factor_cols = [('colA', 'True'),
                   ('colD', 'True')]

    rpy_test_df = com.convert_to_r_dataframe(test_data)

    rpy_out_df = Rtools.convert_columns_to_factors(rpy_test_df, factor_cols)
    test_cols = [('colA', 'factor'),
                 ('colB', 'numeric'),
                 ('colC', 'numeric'),
                 ('colD', 'factor')]

    for col, typ in test_cols:
        if typ == 'factor':
            yield eq_, rpy_out_df.rx2(col).nlevels, 2
        elif typ == 'numeric':
            yield ok_, (not hasattr(rpy_out_df.rx2(col), 'nlevels'))

 def test_far_apart_clusters_estimate_all(self):
     cluster1 = sp.randn(40,1000)
     cluster2 = sp.randn(40,1000) * 2
     cluster2[0,:] += 10
     clusterList1 = [cluster1[:,i]
                     for i in xrange(sp.size(cluster1,1))]
     clusterList2 = [cluster2[:,i]
                     for i in xrange(sp.size(cluster2,1))]
     total, pair = qa.overlap_fp_fn(
         {1: clusterList1, 2: clusterList2})
     self.assertLess(total[1][0], 1e-4)
     self.assertLess(total[1][1], 1e-4)
     self.assertLess(total[2][0], 1e-4)
     self.assertLess(total[2][1], 1e-4)
     self.assertLess(pair[1][2][0], 1e-4)
     self.assertLess(pair[1][2][1], 1e-4)
     self.assertLess(pair[2][1][0], 1e-4)
     self.assertLess(pair[2][1][1], 1e-4)
     self.assertGreater(total[1][0], 0.0)
     self.assertGreater(total[1][1], 0.0)
     self.assertGreater(total[2][0], 0.0)
     self.assertGreater(total[2][1], 0.0)
     self.assertGreater(pair[1][2][0], 0.0)
     self.assertGreater(pair[1][2][1], 0.0)
     self.assertGreater(pair[2][1][0], 0.0)
     self.assertGreater(pair[2][1][1], 0.0)

 def simulate(self,standardize=True):
     self._update_cache()
     RV = SP.zeros((self.N,self.P))
     # region
     Z = SP.randn(self.S,self.P)
     Sc,Uc = LA.eigh(self.Cr.K())
     Sc[Sc<1e-9] = 0
     USh_c = Uc*Sc[SP.newaxis,:]**0.5 
     RV += SP.dot(SP.dot(self.Xr,Z),USh_c.T)
     # background
     Z = SP.randn(self.N,self.P)
     USh_r = self.cache['Lr'].T*self.cache['Srstar'][SP.newaxis,:]**0.5
     Sc,Uc = LA.eigh(self.Cg.K())
     Sc[Sc<1e-9] = 0
     USh_c = Uc*Sc[SP.newaxis,:]**0.5
     RV += SP.dot(SP.dot(USh_r,Z),USh_c.T)
     # noise
     Z = SP.randn(self.N,self.P)
     Sc,Uc = LA.eigh(self.Cn.K())
     Sc[Sc<1e-9] = 0
     USh_c = Uc*Sc[SP.newaxis,:]**0.5 
     RV += SP.dot(Z,USh_c.T)
     # standardize
     if standardize:
         RV-=RV.mean(0)
         RV/=RV.std(0) 
     return RV

    def _genBgTerm_fromXX(self,vTot,vCommon,XX,a=None,c=None):
        """
        generate background term from SNPs

        Args:
            vTot: variance of Yc+Yi
            vCommon: variance of Yc
            XX: kinship matrix
            a: common scales, it can be set for debugging purposes
            c: indipendent scales, it can be set for debugging purposes
        """
        vSpecific = vTot-vCommon

        SP.random.seed(0)
        if c==None: c = SP.randn(self.P)
        XX += 1e-3 * SP.eye(XX.shape[0])
        L = LA.cholesky(XX,lower=True)

        # common effect
        R = self.genWeights(self.N,self.P)
        A = self.genTraitEffect()
        if a is not None: A[0,:] = a
        Yc = SP.dot(L,SP.dot(R,A))
        Yc*= SP.sqrt(vCommon)/SP.sqrt(Yc.var(0).mean())

        # specific effect
        R = SP.randn(self.N,self.P)
        Yi = SP.dot(L,SP.dot(R,SP.diag(c)))
        Yi*= SP.sqrt(vSpecific)/SP.sqrt(Yi.var(0).mean())

        return Yc, Yi

 def _sim_from(self, set_covar='block', seed=None, qq=False):
     ##1. region term
     if set_covar=='block':
         Cr = self.block['Cr']
         Cg = self.block['Cg']
         Cn = self.block['Cn']
     if set_covar=='rank1':
         Cr = self.lr['Cr']
         Cg = self.lr['Cg']
         Cn = self.lr['Cn']
     Lc = msqrt(Cr)
     U, Sh, V = nla.svd(self.Xr, full_matrices=0)
     Lr = sp.zeros((self.Y.shape[0], self.Y.shape[0]))
     Lr[:, :Sh.shape[0]] = U * Sh[sp.newaxis, :]
     Z = sp.randn(*self.Y.shape)
     Yr = sp.dot(Lr, sp.dot(Z, Lc.T))
     ##2. bg term
     Lc = msqrt(Cg)
     Lr = self.XXh
     Z = sp.randn(*self.Y.shape)
     Yg = sp.dot(Lr, sp.dot(Z, Lc.T))
     # noise terms
     Lc = msqrt(Cn)
     Z = sp.randn(*self.Y.shape)
     Yn = sp.dot(Z, Lc.T)
     # normalize
     Y = Yr + Yg + Yn
     if qq:
         Y = gaussianize(Y)
         Y-= Y.mean(0)
         Y/= Y.std(0)
     return Y

def test_mixed_model():

    test_data = DataFrame(
        {
            'colA': Series(randn(1, 5000).flatten() > 0),
            'colB': Series(100 * randn(1, 5000).flatten()),
            'colC': Series(100 + randn(1, 5000).flatten()),
            'colD': Series(randn(1, 5000).flatten() > 0),
            },
        )

    test_data['colA'] = test_data['colA'].map(str)
    test_data['colD'] = test_data['colD'].map(str)

    factor_cols = [('colA', 'True'),
                   ('colD', 'True')]

    rpy_test_df = com.convert_to_r_dataframe(test_data)
    rpy_test_df = Rtools.convert_columns_to_factors(rpy_test_df, factor_cols)

    base_formula = Formula('colC ~ as.factor(colA) + colB')
    rand_formula = Formula('~1|colD')

    results = Rtools.R_linear_mixed_effects_model(rpy_test_df, base_formula, rand_formula)

    print results['tTable']
    ok_(('tTable' in results), 'Did not have the tTable in the results')
    ok_(('as.factor(colA)False' in results['tTable'].index), 'Did not have the factor in the tTable')
    ok_(('colB' in results['tTable'].index), 'Did not have the variable in the tTable')

def gendat(TWO_KERNEL,nInd,nSnp,nCovar,minMaf=0.05,maxMaf=0.4,minSigE2=0.5,maxSigE2=1,minSigG2=0.5,maxSigG2=1):
    '''
    Generate synthetic SNPs and phenotype.
    SNPs are iid, and there is no population structure.
    Phenotype y is generated from a LMM with SNPs in a PS kernel.

    Returns:
        covDat
        y
        psSnps
    '''

    if TWO_KERNEL:
        psSnps=gensnps(nInd,nSnp,minMaf,maxMaf)
        psK=psSnps.dot(psSnps.T)
        psK+=1e-5*sp.eye(nInd)
        psKchol=la.cholesky(psK)
    else:
        psSnps=None

    covDat=sp.random.uniform(0,1,(nInd,nCovar))
    covWeights=sp.random.uniform(-0.5,0.5,(nCovar,1))

    sigE2=sp.random.uniform(low=minSigE2,high=maxSigE2)
    sigG2=sp.random.uniform(low=minSigG2,high=maxSigG2)

    ##generate the phenotype using the background kernel and covariates
    if TWO_KERNEL:
        y_pop=sp.sqrt(sigG2)*psKchol.dot(sp.randn(nInd,1))
    else:
        y_pop=0
    y_noise=sp.randn(nInd,1)*sp.sqrt(sigE2)
    y=(covDat.dot(covWeights) + y_pop + y_noise).squeeze()
    return covDat, y, psSnps

def writeBackRepsAddNoise(wc1,wc2,y1,y2,geneName,n_reps):
    for i in range(n_reps):
        c1 = [geneName]
        c2 = [geneName]
        c1.extend(y1+SP.randn(y1.shape[0])*.1)
        c2.extend(y2+SP.randn(y2.shape[0])*.1)
        wc1.writerow(c1)
        wc2.writerow(c2)

def awgn(sig,snrdb,sigpower=0):
    """Additive white gaussian noise.  Assumes signal power is 0 dBW"""
    if sp.iscomplexobj(sig):
        noise = (sp.randn(*sig.shape) + 1j*sp.randn(*sig.shape))/math.sqrt(2)
    else:
        noise = sp.randn(*sig.shape)
    noisev = 10**((sigpower - snrdb)/20)
    return sig + noise*noisev

 def simulate_pheno(self):
     Yc = sp.dot(self.mean.F[0], sp.dot(self.mean.B[0], self.mean.A[0].T))
     Z = sp.randn(self.covar.G.shape[1], self.covar.Cr.X.shape[1])
     Yr = sp.dot(self.covar.G, sp.dot(Z, self.covar.Cr.X.T))
     _S, _U = LA.eigh(self.covar.Cn.K()); _S[_S<0] = 0
     Cn_h = _U*_S**0.5
     Yn = sp.dot(sp.randn(*self.mean.Y.shape), Cn_h.T)
     RV = Yc+Yr+Yn
     return RV

 def sim_bivariate(self, N):
     R = np.empty(N)
     X = np.empty(N)
     X[0] = 1
     for t in range(N-1):
         R[t] = sqrt(X[t]) * randn(1)
         X[t+1] = self.a0 + self.b * X[t] + self.a1 * R[t]**2
     R[N-1] = sqrt(X[N-1]) * randn(1)
     return R, X

 def ts(self, n):
     R = np.empty(n)
     X = np.empty(n)
     X[0] = 1
     for t in range(n-1):
         R[t] = np.sqrt(X[t]) * randn(1)
         X[t+1] = self.a0 + self.b * X[t] + self.a1 * R[t]**2
     R[n-1] = np.sqrt(X[n-1]) * randn(1)
     return R

	def _initParams_random(self):
		""" 
		initialize the gp parameters randomly
		"""
		# gp hyper params
		params = limix.CGPHyperParams()
		if self.interaction:	params['covar'] = SP.concatenate([SP.randn(self.N*self.k+1),SP.ones(1),SP.randn(1)])
		else:					params['covar'] = SP.randn(self.N*self.k+1)
		params['lik'] = SP.randn(1)
		return params

 def _additionalInit(self):
     phi_size = self.num_actions * self.num_features
     if self.randomInit:
         self._A = randn(phi_size, phi_size) / 100.
         self._b = randn(phi_size) / 100.
     else:
         self._A = zeros((phi_size, phi_size))
         self._b = zeros(phi_size)          
     self._untouched = ones(phi_size, dtype=bool)
     self._count = 0