Internal function for parameter initialization
        Uses a single trait model for initializing variances and
        a pairwise model to initialize correlations

def _getScalesPairwise(self,verbose=False, initDiagonal=False):
        """
        Internal function for parameter initialization
        Uses a single trait model for initializing variances and
        a pairwise model to initialize correlations
        """
        var = sp.zeros((self.P,2))

        if initDiagonal:
            #1. fit single trait model
            if verbose:
                print('.. fit single-trait model for initialization')
            vc = VarianceDecomposition(self.Y[:,0:1])
            for term_i in range(self.n_randEffs):
                if term_i==self.noisPos:
                    vc.addRandomEffect(is_noise=True)
                else:
                    K = self.vd.getTerm(term_i).getK()
                    vc.addRandomEffect(K=K)
            scales0 = sp.sqrt(0.5)*sp.ones(2)

            for p in range(self.P):
                if verbose: print(('   .. trait %d' % p))
                vc.setY(self.Y[:,p:p+1])
                conv = vc.optimize(scales0=scales0)
                if not conv:
                    print('warning initialization not converged')
                var[p,:] = vc.getVarianceComps()[0,:]

        elif fastlmm_present:
            if verbose:
                print('.. fit single-trait model for initialization (using fastlmm)')
            for p in range(self.P):
                if verbose: print(('   .. trait %d' % p))
                covariates = None
                for term_i in range(self.n_randEffs):
                    if term_i==self.noisPos:
                        pass
                    else:
                        K = self.vd.getTerm(term_i).getK()
                varY = sp.var(self.Y[:,p:p+1])
                lmm = fastLMM(X=covariates, Y=self.Y[:,p:p+1], G=None, K=K)
                opt = lmm.findH2(nGridH2=100)
                h2 = opt['h2']
                var[p,:] = h2 * varY
                var[p,self.noisPos] = (1.0-h2) * varY
                #import ipdb;ipdb.set_trace()
        else:
            if verbose:
                print('.. random initialization of diagonal')
            var = sp.random.randn(var.shape[0],var.shape[1])
            var = var*var + 0.001
        #2. fit pairwise model
        if verbose:
            print('.. fit pairwise model for initialization')
        vc = VarianceDecomposition(self.Y[:,0:2])
        for term_i in range(self.n_randEffs):
            if term_i==self.noisPos:
                vc.addRandomEffect(is_noise=True,trait_covar_type='freeform')
            else:
                K = self.vd.getTerm(term_i).getK()
                vc.addRandomEffect(K=K,trait_covar_type='freeform')
        rho_g = sp.ones((self.P,self.P))
        rho_n = sp.ones((self.P,self.P))
        for p1 in range(self.P):
            for p2 in range(p1):
                if verbose:
                    print(('   .. fit pair (%d,%d)'%(p1,p2)))
                vc.setY(self.Y[:,[p1,p2]])
                scales0 = sp.sqrt(sp.array([var[p1,0],1e-4,var[p2,0],1e-4,var[p1,1],1e-4,var[p2,1],1e-4]))
                conv = vc.optimize(scales0=scales0)
                if not conv:
                    print('warning initialization not converged')
                Cg = vc.getTraitCovar(0)
                Cn = vc.getTraitCovar(1)
                rho_g[p1,p2] = Cg[0,1]/sp.sqrt(Cg.diagonal().prod())
                rho_n[p1,p2] = Cn[0,1]/sp.sqrt(Cn.diagonal().prod())
                rho_g[p2,p1] = rho_g[p1,p2]
                rho_n[p2,p1] = rho_n[p1,p2]
        #3. init
        Cg0 = rho_g*sp.dot(sp.sqrt(var[:,0:1]),sp.sqrt(var[:,0:1].T))
        Cn0 = rho_n*sp.dot(sp.sqrt(var[:,1:2]),sp.sqrt(var[:,1:2].T))
        offset_g = abs(sp.minimum(sp.linalg.eigh(Cg0)[0].min(),0))+1e-4
        offset_n = abs(sp.minimum(sp.linalg.eigh(Cn0)[0].min(),0))+1e-4
        Cg0+=offset_g*sp.eye(self.P)
        Cn0+=offset_n*sp.eye(self.P)
        Lg = sp.linalg.cholesky(Cg0)
        Ln = sp.linalg.cholesky(Cn0)
        Cg_params0 = sp.concatenate([Lg[:,p][:p+1] for p in range(self.P)])
        Cn_params0 = sp.concatenate([Ln[:,p][:p+1] for p in range(self.P)])
        scales0 = sp.concatenate([Cg_params0,1e-2*sp.ones(1),Cn_params0,1e-2*sp.ones(1)])

        return scales0