def genInitSigmaFactor(self):
     """ depending on the algorithm settings, we start out with in identity matrix, or perturb it """
     if self.perturbedInitSigma:
         res = mat(eye(self.xdim)*self.initSigmaCoeff+randn(self.xdim, self.xdim)*self.initSigmaRandCoeff)            
     else:
         res = mat(eye(self.xdim)*self.initSigmaCoeff)
     return res   

    def GetAugMat(self, s, sym=False):
        """Return the augmented element transfer matrix for the
        AVS1_kp element."""
        N = self.maxsize
        if sym:
            myparams = self.symparams
            matout = eye(N + 1, dtype="f")
            matout = matout.astype("S30")
        else:
            matout = eye(N + 1, dtype="D")
            myparams = self.params
        myrow = 1  # hard coding for now#(self.params['axis']-1)*4+1#axis should be 1, 2, or 3
        Gact = self.Gact_func(s, self.params)
        Gth = self.kp
        k_spring = self.params["k_spring"]
        c_spring = self.params["c_spring"]
        H = self.params["H"]
        term1 = 1.0 / ((1.0 + Gact * Gth * H) * (k_spring + c_spring * s))
        term2 = Gact * Gth / (1.0 + Gact * Gth * H)
        # term1 = 1.0/(k_spring + c_spring*s + Gact*Gth*k_spring + Gact*Gth*c_spring*s)
        # term2 = Gact*Gth/(1.0 + Gact*Gth)

        myrow = 1  # hard coding for now#(self.params['axis']-1)*4+1#axis should be 1, 2, or 3
        matout[myrow, 2] = term1
        matout[myrow, N] = term2
        return matout

def compute_diagonal_loading(mat, svd, target_cond=SUFFICIENT_CONDITION,
                             overwrite_mat=False):
    """tries to condition :mat: by imposing a spherical constraint on the
    covariance ellipsoid (adding alpha*eye)

    solves: cond(mat + alpha*I) = target_cond for alpha

    Note: this is a noop if the condition is already >= target_cond!

    :type mat: ndarray
    :param mat: input matrix
    :type svd: tuple
    :param svd: return tuple of svd(:mat:) - consistency will not be checked!
    :type target_cond: float
    :param target_cond: condition number to archive after loading
    :type overwrite_mat: bool
    :param overwrite_mat: if True, operate inplace and overwrite :mat:
    :returns: ndarray - matrix like :mat: conditioned s.t. cond = target_cond
    """

    sv = svd[1]
    if target_cond == 1.0:
        return sp.eye(mat.shape[0], mat.shape[1])
    if target_cond > compute_matrix_cond(sv):
        return mat
    if overwrite_mat is True:
        rval = mat
    else:
        rval = mat.copy()
    alpha = (sv[0] - target_cond * sv[-1]) / (target_cond - 1)
    return rval + alpha * sp.eye(rval.shape[0], rval.shape[1])

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 __iter__(self):
        dim = self.wrt.shape[0]
        I = scipy.eye(dim)

        # Square root of covariance matrix.
        A = scipy.eye(dim)
        center = self.wrt.copy()
        n_evals = 0
        best_wrt = None
        best_x = float("-inf")
        for i, (args, kwargs) in enumerate(self.args):
            # Draw samples, evaluate and update best solution if a better one
            # was found.
            samples = scipy.random.standard_normal((self.batch_size, dim))
            samples = scipy.dot(samples, A) + center
            fitnesses = [self.f(samples[j], *args, **kwargs) for j in range(samples.shape[0])]
            fitnesses = scipy.array(fitnesses).flatten()

            if fitnesses.max() > best_x:
                best_loss = fitnesses.max()
                self.wrt[:] = samples[fitnesses.argmax()]

            # Update center and variances.
            utilities = self.compute_utilities(fitnesses)
            center += scipy.dot(scipy.dot(utilities, samples), A)
            # TODO: vectorize this
            cov_gradient = sum([u * (scipy.outer(s, s) - I) for (s, u) in zip(samples, utilities)])
            update = scipy.linalg.expm2(A * cov_gradient * self.step_rate * 0.5)
            A[:] = scipy.dot(A, update)

            yield dict(loss=-best_x, n_iter=i)

    def testFill(self):
        """test buffer filling"""

        self.rb.fill(sp.eye(4))
        self.assertEqual(len(self.rb), 6)
        for item in self.rb:
            assert_equal(item, sp.eye(4))

 def GetMat(self,s,sym=False):
     """Return the element transfer matrix for the
     TorsionalSpringDamper element.  If sym=True, 's' must be a
     symbolic string and a matrix of strings will be returned.
     Otherwise, 's' is a numeric value (probably complex) and the
     matrix returned will be complex."""
     N=self.maxsize
     if sym:
         myparams=self.symparams
     else:
         myparams=self.params
     k=myparams['k']
     c=myparams['c']
     springterm=1/(k[0]+c[0]*s)
     if sym:
         maxlen=len(springterm)+10
         matout=eye(N,dtype='f')
         matout=matout.astype('S%d'%maxlen)
     else:
         matout=eye(N,dtype='D')
     matout[1,2]=springterm
     if max(shape(k))>1 and self.maxsize>=8:
         matout[5,6]=1/(k[1]+c[1]*s)
     if max(shape(k))>2 and self.maxsize>=12:
         matout[9,10]=1/(k[2]+c[2]*s)
     return matout

    def _maximum_likelihood(self, X):
        n_samples, n_features = X.shape if X.ndim > 1 else (1, X.shape[0])
        n_components = self.n_components

        # Predict mean
        mu = X.mean(axis=0)

        # Predict covariance
        cov = sp.cov(X, rowvar=0)
        eigvals, eigvecs = self._eig_decomposition(cov)
        sigma2 = ((sp.sum(cov.diagonal()) - sp.sum(eigvals.sum())) /
                  (n_features - n_components))  # FIXME: M < D?

        weight = sp.dot(eigvecs, sp.diag(sp.sqrt(eigvals - sigma2)))
        M = sp.dot(weight.T, weight) + sigma2 * sp.eye(n_components)
        inv_M = spla.inv(M)

        self.eigvals = eigvals
        self.eigvecs = eigvecs
        self.predict_mean = mu
        self.predict_cov = sp.dot(weight, weight.T) + sigma2 * sp.eye(n_features)
        self.latent_mean = sp.transpose(sp.dot(inv_M, sp.dot(weight.T, X.T - mu[:, sp.newaxis])))
        self.latent_cov = sigma2 * inv_M
        self.sigma2 = sigma2    # FIXME!
        self.weight = weight
        self.inv_M = inv_M

        return self.latent_mean

def find_homog_trans(points_a, points_b, err_threshold=0, rot_0=None):
    """Finds a homogeneous transformation matrix that, when applied to 
    the points in points_a, minimizes the squared Euclidean distance 
    between the transformed points and the corresponding points in 
    points_b. Both points_a and points_b are (n, 3) arrays.
    """
    #Align the centroids of the two point clouds
    cent_a = sp.average(points_a, axis=0)
    cent_b = sp.average(points_b, axis=0)
    points_a = points_a - cent_a
    points_b = points_b - cent_b
    
    #Define the error as a function of a rotation vector in R^3
    rot_cost = lambda rot: (sp.dot(vec_to_rot(rot), points_a.T).T
                    - points_b).flatten()**2
    
    #Run the optimization
    if rot_0 == None:
        rot_0 = sp.zeros(3)
    rot = opt.leastsq(rot_cost, rot_0)[0]
    
    #Compute the final homogeneous transformation matrix
    homog_1 = sp.eye(4)
    homog_1[0:3, 3] = -cent_a
    homog_2 = sp.eye(4)
    homog_2[0:3,0:3] = vec_to_rot(rot)
    homog_3 = sp.eye(4)
    homog_3[0:3,3] = cent_b
    homog = sp.dot(homog_3, sp.dot(homog_2, homog_1))
    return homog, rot

def sample_moments( X, k ):
    """Get the sample moments from data"""
    N, d = X.shape

    # Partition X into two halves to independently estimate M2 and M3
    X1, X2 = X[:N/2], X[N/2:]

    # Get the moments  
    M1 = X1.mean(0)
    M1_ = X2.mean(0)
    M2 = Pairs( X1, X1 ) 
    M3 = lambda theta: TriplesP( X2, X2, X2, theta )
    #M3 = Triples( X2, X2, X2 )

    # TODO: Ah, not computing sigma2! 
    # Estimate \sigma^2 = k-th eigenvalue of  M2 - mu mu^T
    sigma2 = svdvals( M2 - outer( M1, M1 ) )[k-1]
    assert( sc.isreal( sigma2 ) and sigma2 > 0 )
    # P (M_2) is the best kth rank apprximation to M2 - sigma^2 I
    P = approxk( M2 - sigma2 * eye( d ), k )

    B = matrix_tensorify( eye(d), M1_ )
    T = lambda theta: M3(theta) - sigma2 * ( M1_.dot(theta) * eye( d ) + outer( M1_, theta ) + outer( theta, M1_ ) )
    #T = M3 - sigma2 * ( B + B.swapaxes(2, 1) + B.swapaxes(2, 0) )

    return P, T    

def initialize_batch(X_bar0, P_bar0, x_bar0):
    """
    Generate t=0 values for a new iteration from an initial state, covariance
    and a-priori estimate.
    """
    # Get initial state and STM and initialize integrator
    X_bar0_list = X_bar0.T.tolist()[0]
    stm0 = sp.matrix(sp.eye(18))
    stm0_list = sp.eye(18).reshape(1,324).tolist()[0]

    eom = ode(Udot).set_integrator('dop853', atol=1.0E-10, rtol=1.0E-9)
    eom.set_initial_value(X_bar0_list + stm0_list, 0)

    # Accumulate measurement at t=0
    obs0 = OBS[0]
    stn0 = obs0[0]
    comp0, Htilda0 = Htilda_matrix(X_bar0_list, 0, stn0)
    resid0 = [ obs0[1] - float(comp0[0]),
               obs0[2] - float(comp0[1]) ]
    y0 = sp.matrix([resid0]).T
    H0 = Htilda0 * stm0

    L0 = P_bar0.I + H0.T * W * H0
    N0 = P_bar0.I * x_bar0 + H0.T * W * y0

    return [stm0, comp0, resid0, Htilda0, H0, L0, N0, eom]

def xNES(f, x0, maxEvals=1e6, verbose=False, targetFitness= -1e-10):
    """ Exponential NES (xNES), as described in 
    Glasmachers, Schaul, Sun, Wierstra and Schmidhuber (GECCO'10).
    Maximizes a function f. 
    Returns (best solution found, corresponding fitness).
    """
    dim = len(x0)  
    I = eye(dim)
    learningRate = 0.6 * (3 + log(dim)) / dim / sqrt(dim)
    batchSize = 4 + int(floor(3 * log(dim)))    
    center = x0.copy()
    A = eye(dim)  # sqrt of the covariance matrix
    numEvals = 0
    bestFound = None
    bestFitness = -Inf
    while numEvals + batchSize <= maxEvals and bestFitness < targetFitness:
        # produce and evaluate samples
        samples = [randn(dim) for _ in range(batchSize)]
        fitnesses = [f(dot(A, s) + center) for s in samples]
        if max(fitnesses) > bestFitness:
            bestFitness = max(fitnesses)
            bestFound = samples[argmax(fitnesses)]
        numEvals += batchSize 
        if verbose: print "Step", numEvals / batchSize, ":", max(fitnesses), "best:", bestFitness
        #print A
        # update center and variances
        utilities = computeUtilities(fitnesses)
        center += dot(A, dot(utilities, samples))
        covGradient = sum([u * (outer(s, s) - I) for (s, u) in zip(samples, utilities)])
        A = dot(A, expm2(0.5 * learningRate * covGradient))                      

    return bestFound, bestFitness

    def testSnrFuncs(self):
        """test for signal to noise ratio functions"""

        # trivial
        data_triv = sp.ones((3, 10))
        snr_triv_test = sp.ones(3)
        assert_equal(
            snr_peak(data_triv, 1.0),
            snr_triv_test)
        assert_equal(
            snr_power(data_triv, 1.0),
            snr_triv_test)
        assert_equal(
            snr_maha(data_triv, sp.eye(data_triv.shape[1])),
            snr_triv_test)

        # application
        data = sp.array([
            sp.sin(sp.linspace(0.0, 2 * sp.pi, 100)),
            sp.sin(sp.linspace(0.0, 2 * sp.pi, 100)) * 2,
            sp.sin(sp.linspace(0.0, 2 * sp.pi, 100)) * 5,
        ])
        assert_equal(
            snr_peak(data, 1.0),
            sp.absolute(data).max(axis=1))
        assert_equal(
            snr_power(data, 1.0),
            sp.sqrt((data * data).sum(axis=1) / data.shape[1]))
        assert_almost_equal(
            snr_maha(data, sp.eye(data.shape[1])),
            sp.sqrt((data * data).sum(axis=1) / data.shape[1]))

 def __init__(self, basef, 
              translate=True, 
              rotate=False, 
              conditioning=None, 
              asymmetry=None,
              oscillate=False, 
              penalize=None,
              ):
     FunctionEnvironment.__init__(self, basef.xdim, basef.xopt)
     self.desiredValue = basef.desiredValue            
     self.toBeMinimized = basef.toBeMinimized
     
     if translate:            
         self.xopt = (rand(self.xdim) - 0.5) * 9.8
         
     self._diags = eye(self.xdim)            
     self._R = eye(self.xdim)            
     self._Q = eye(self.xdim)            
     
     if conditioning is not None:
         self._diags = generateDiags(conditioning, self.xdim)
     if rotate:
         self._R = orth(rand(basef.xdim, basef.xdim))        
         if conditioning:
             self._Q = orth(rand(basef.xdim, basef.xdim))
                     
     tmp = lambda x: dot(self._Q, dot(self._diags, dot(self._R, x-self.xopt)))
     if asymmetry is not None:
         tmp2 = tmp
         tmp = lambda x: asymmetrify(tmp2(x), asymmetry)
     if oscillate:
         tmp3 = tmp
         tmp = lambda x: oscillatify(tmp3(x))
     
     self.f = lambda x: basef.f(tmp(x))

def sqrtm3(X):
    M = sp.copy(X)
    m, fb, fe = block_structure(M)
    n = M.shape[0]
    for i in range(0,m):
        M[fb[i]:fe[i],fb[i]:fe[i]] = twobytworoot(M[fb[i]:fe[i],fb[i]:fe[i]])
        #print M

    for j in range(1,m):
        for i in range(0,m-j):
            #print M[fb[i]:fe[i],fb[JJ]:fe[JJ]]
            JJ = i+j
            Tnoto = M[fb[i]:fe[i],fb[JJ]:fe[JJ]] #dopo togliere il copy
            #print "Tnot: "
            #print Tnoto
            for k in range(i+1,JJ):
                Tnoto -= (M[fb[i]:fe[i],fb[k]:fe[k]]).dot(M[fb[k]:fe[k],fb[JJ]:fe[JJ]])
                #print M[fb[i]:fe[i],fb[k]:fe[k]]
                #print M[fb[k]:fe[k],fb[JJ]:fe[JJ]]

            if((M[fb[i]:fe[i],fb[JJ]:fe[JJ]]).shape==(1,1)):
                #print "forma 1"
                #print M[fb[i]:fe[i],fb[JJ]:fe[JJ]]           #  Uij
                #print M[fb[i]:fe[i],fb[i]:fe[i]]               #  Uii
                #print M[fb[JJ]:fe[JJ],fb[JJ]:fe[JJ]]       #  Ujj
                M[fb[i]:fe[i],fb[JJ]:fe[JJ]] = Tnoto/(M[fb[i]:fe[i],fb[i]:fe[i]] + M[fb[JJ]:fe[JJ],fb[JJ]:fe[JJ]])

            else:
                Uii = M[fb[i]:fe[i],fb[i]:fe[i]]
                Ujj = M[fb[JJ]:fe[JJ],fb[JJ]:fe[JJ]]
                shapeUii = Uii.shape[0]
                shapeUjj = Ujj.shape[0]
                """
                print "------------"
                print Tnoto
                print Tnoto.shape
                print sp.kron(sp.eye(shapeUjj),Uii)
                print sp.kron(Ujj.T,sp.eye(shapeUii))
                print Tnoto
                """
                #M[fb[i]:fe[i],fb[JJ]:fe[JJ]] = sp.linalg.solve_sylvester(Uii, Ujj, Tnoto)

                """
                x, scale, info = dtrsyl(Uii, Ujj, Tnoto

                if (scale==1.0):
                     = x

                else:
                    M[fb[i]:fe[i],fb[JJ]:fe[JJ]] = x*scale
                    print "scale!=0"
                """
                Tnoto = Tnoto.reshape((shapeUii*shapeUjj),1,order="F")
                M[fb[i]:fe[i],fb[JJ]:fe[JJ]] = \
                linalg.solve(sp.kron(sp.eye(shapeUjj),Uii) +
                sp.kron(Ujj.T,sp.eye(shapeUii)),
                Tnoto).reshape(shapeUii,shapeUjj,order="F")


    return M

    def __init__(self, evaluator, evaluable, **parameters):
        BlackBoxOptimizer.__init__(self, evaluator, evaluable, **parameters)
        self.alphas = ones(self.numberOfCenters)/self.numberOfCenters
        self.mus = []
        self.sigmas = []

        self.tau = 1.
        if self.rangemins == None:
            self.rangemins = -ones(self.xdim)
        if self.rangemaxs == None:
            self.rangemaxs = ones(self.xdim)
        if self.initCovariances == None:
            self.initCovariances = eye(self.xdim)
            
        if self.elitist and self.numberOfCenters == 1 and not self.noisyEvaluator:
            # in the elitist case seperate evaluations are not necessary. 
            # CHECKME: maybe in the noisy case?
            self.evalMus = False
            
        assert not(self.useCauchy and self.numberOfCenters > 1)
            
        for dummy in range(self.numberOfCenters):
            self.mus.append(rand(self.xdim) * (self.rangemaxs-self.rangemins) + self.rangemins)
            self.sigmas.append(dot(eye(self.xdim), self.initCovariances))
        self.reset()

    def LMLdebug(self):
        """
        LML function for debug
        """
        assert self.N*self.P<5000, 'gp2kronSum:: N*P>=5000'

        y = SP.reshape(self.Y,(self.N*self.P), order='F') 
        V = SP.kron(SP.eye(self.P),self.F)

        XX = SP.dot(self.Xr,self.Xr.T)
        K  = SP.kron(self.Cr.K(),XX)
        K += SP.kron(self.Cn.K()+self.offset*SP.eye(self.P),SP.eye(self.N))

        # inverse of K
        cholK = LA.cholesky(K)
        Ki = LA.cho_solve((cholK,False),SP.eye(self.N*self.P))

        # Areml and inverse
        Areml = SP.dot(V.T,SP.dot(Ki,V))
        cholAreml = LA.cholesky(Areml)
        Areml_i = LA.cho_solve((cholAreml,False),SP.eye(self.K*self.P))

        # effect sizes and z
        b = SP.dot(Areml_i,SP.dot(V.T,SP.dot(Ki,y)))
        z = y-SP.dot(V,b)
        Kiz = SP.dot(Ki,z)

        # lml
        lml  = y.shape[0]*SP.log(2*SP.pi)
        lml += 2*SP.log(SP.diag(cholK)).sum()
        lml += 2*SP.log(SP.diag(cholAreml)).sum()
        lml += SP.dot(z,Kiz)
        lml *= 0.5

        return lml

    def __init__(self, evaluator, evaluable, **parameters):
        BlackBoxOptimizer.__init__(self, evaluator, evaluable, **parameters)
        
        self.numParams = self.xdim + self.xdim * (self.xdim+1) / 2
                
        if self.momentum != None:
            self.momentumVector = zeros(self.numParams)
        if self.learningRateSigma == None:
            self.learningRateSigma = self.learningRate
        
        if self.rangemins == None:
            self.rangemins = -ones(self.xdim)
        if self.rangemaxs == None:
            self.rangemaxs = ones(self.xdim)
        if self.initCovariances == None:
            if self.diagonalOnly:
                self.initCovariances = ones(self.xdim)
            else:
                self.initCovariances = eye(self.xdim)

        self.x = rand(self.xdim) * (self.rangemaxs-self.rangemins) + self.rangemins
        self.sigma = dot(eye(self.xdim), self.initCovariances)
        self.factorSigma = cholesky(self.sigma)
        
        self.reset()

    def test_lowrank_iso(self):
        theta = SP.array(SP.random.randn(2)**2)
        theta_hat = SP.exp(2*theta)

        _K = theta_hat[0]*SP.dot(self.Xtrain,self.Xtrain.T) + theta_hat[1]*SP.eye(self.n_train)
        _Kcross = theta_hat[0]*SP.dot(self.Xtrain,self.Xtest.T)
        _Kgrad_theta = []
        _Kgrad_theta.append(2*theta_hat[0]*SP.dot(self.Xtrain,self.Xtrain.T) )
        _Kgrad_theta.append(2*theta_hat[1]*SP.eye(self.n_train))

        cov = lowrank.LowRankCF(self.n_dimensions)
        cov.X = self.Xtrain
        cov.Xcross = self.Xtest
        
        K = cov.K(theta)
        Kcross = cov.Kcross(theta)

        assert SP.allclose(K,_K), 'ouch, covariance matrix is wrong'
        assert SP.allclose(Kcross,_Kcross), 'ouch, cross covariance matrix is wrong'
        assert SP.allclose(_Kgrad_theta[0],cov.Kgrad_theta(theta,0))
        assert SP.allclose(_Kgrad_theta[1],cov.Kgrad_theta(theta,1))

        # gradient with respect to latent factors
        for i in range(self.n_dimensions):
            for j in range(self.n_train):
                Xgrad = SP.zeros(self.Xtrain.shape)
                Xgrad[j,i] = 1
                _Kgrad_x =  theta_hat[0]*(SP.dot(Xgrad,self.Xtrain.T) + SP.dot(self.Xtrain,Xgrad.T))
                Kgrad_x = cov.Kgrad_x(theta,i,j)
                assert SP.allclose(Kgrad_x,_Kgrad_x), 'ouch, gradient with respect to x is wrong for entry [%d,%d]'%(i,j)

    def _LMLgrad_covar_debug(self,covar):

        assert self.N*self.P<2000, 'gp2kronSum:: N*P>=2000'

        y  = SP.reshape(self.Y,(self.N*self.P), order='F') 

        K  = SP.kron(self.Cg.K(),self.XX)
        K += SP.kron(self.Cn.K()+self.offset*SP.eye(self.P),SP.eye(self.N))

        cholK = LA.cholesky(K).T
        Ki  = LA.cho_solve((cholK,True),SP.eye(y.shape[0]))
        Kiy   = LA.cho_solve((cholK,True),y)

        if covar=='Cr':     n_params = self.Cr.getNumberParams()
        elif covar=='Cg':   n_params = self.Cg.getNumberParams()
        elif covar=='Cn':   n_params = self.Cn.getNumberParams()

        RV = SP.zeros(n_params)

        for i in range(n_params):
            #0. calc grad_i
            if covar=='Cg':
                C   = self.Cg.Kgrad_param(i)
                Kgrad  = SP.kron(C,self.XX)
            elif covar=='Cn':
                C   = self.Cn.Kgrad_param(i)
                Kgrad  = SP.kron(C,SP.eye(self.N))

            #1. der of log det
            RV[i]  = 0.5*(Ki*Kgrad).sum()
            
            #2. der of quad form
            RV[i] -= 0.5*(Kiy*SP.dot(Kgrad,Kiy)).sum()

        return RV