def addMultiTraitTerm(self,K=None,covar_type='freeform',is_noise=False,normalize=True,Ks=None,offset=1e-4,rank=1,covar_K0=None):
        """
        add multi trait random effects term.
        The inter-trait covariance is parametrized by covar_type, where parameters are optimized.
        
        Args:
            K:      Individual-individual (Intra-Trait) Covariance Matrix [N, N]
                    (K is normalised in the C++ code such that K.trace()=N)
            covar_type: type of covaraince to use. Default 'freeform'. possible values are 
                            'freeform': free form optimization, 
                            'fixed': use a fixed matrix specified in covar_K0,
                            'diag': optimize a diagonal matrix, 
                            'lowrank': optimize a low rank matrix. The rank of the lowrank part is specified in the variable rank,
                            'lowrank_id': optimize a low rank matrix plus the weight of a constant diagonal matrix. The rank of the lowrank part is specified in the variable rank, 
                            'lowrank_diag': optimize a low rank matrix plus a free diagonal matrix. The rank of the lowrank part is specified in the variable rank, 
                            'block': optimize the weight of a constant P x P block matrix of ones,
                            'block_id': optimize the weight of a constant P x P block matrix of ones plus the weight of a constant diagonal matrix,
                            'block_diag': optimize the weight of a constant P x P block matrix of ones plus a free diagonal matrix,                            
            is_noise:   Boolean indicator specifying if the matrix is homoscedastic noise (weighted identity covariance) (default False)
            normalize:  Boolean indicator specifying if K is normalized such that K.trace()=N.
            Ks:			NxNtest cross covariance for predictions
            offset:		diagonal contribution added to trait-to-trait covariance matrices for regularization
            rank:       rank of a possible lowrank component (default 1)
            covar_K0:   PxP matrix for the (predefined) trait-to-trait covariance matrix if fixed type is used
        """
        assert self.P > 1, 'CVarianceDecomposition:: Incompatible number of traits'
        assert K!=None or is_noise, 'CVarianceDecomposition:: Specify covariance structure'
        assert offset>=0, 'CVarianceDecomposition:: offset must be >=0'

        #TODO: check that covar_K0 is correct if fixed typeCF is used..

        if is_noise:
            assert self.noisPos==None, 'CVarianceDecomposition:: noise term already exists'
            K = SP.eye(self.N)
            self.noisPos = self.n_terms
        else:
            assert K.shape[0]==self.N, 'CVarianceDecomposition:: Incompatible shape'
            assert K.shape[1]==self.N, 'CVarianceDecomposition:: Incompatible shape'

        if Ks!=None:
            assert Ks.shape[0]==self.N, 'CVarianceDecomposition:: Incompatible shape'

        if normalize:
            Norm = 1/K.diagonal().mean()
            K *= Norm
            if Ks!=None: Ks *= Norm

        cov = limix.CSumCF()
        if covar_type=='freeform':
            cov.addCovariance(limix.CFreeFormCF(self.P))
            L = SP.eye(self.P)
            diag = SP.concatenate([L[i,:(i+1)] for i in range(self.P)])
        elif covar_type=='fixed':
            cov.addCovariance(limix.CFixedCF(covar_K0))
            diag = SP.zeros(1)
        elif covar_type=='diag':
            cov.addCovariance(limix.CDiagonalCF(self.P))
            diag = SP.ones(self.P)
        elif covar_type=='lowrank':
            cov.addCovariance(limix.CLowRankCF(self.P,rank))
            diag = SP.zeros(self.P*rank)
        elif covar_type=='lowrank_id':
            cov.addCovariance(limix.CLowRankCF(self.P,rank))
            cov.addCovariance(limix.CFixedCF(SP.eye(self.P)))
            diag = SP.concatenate([SP.zeros(self.P*rank),SP.ones(1)])
        elif covar_type=='lowrank_diag':
            cov.addCovariance(limix.CLowRankCF(self.P,rank))
            cov.addCovariance(limix.CDiagonalCF(self.P))
            diag = SP.concatenate([SP.zeros(self.P*rank),SP.ones(self.P)])
        elif covar_type=='block':
            cov.addCovariance(limix.CFixedCF(SP.ones((self.P,self.P))))
            diag = SP.zeros(1)
        elif covar_type=='block_id':
            cov.addCovariance(limix.CFixedCF(SP.ones((self.P,self.P))))
            cov.addCovariance(limix.CFixedCF(SP.eye(self.P)))
            diag = SP.concatenate([SP.zeros(1),SP.ones(1)])
        elif covar_type=='block_diag':
            cov.addCovariance(limix.CFixedCF(SP.ones((self.P,self.P))))
            cov.addCovariance(limix.CDiagonalCF(self.P))
            diag = SP.concatenate([SP.zeros(1),SP.ones(self.P)])
        else:
            assert True==False, 'CVarianceDecomposition:: covar_type not valid'

        if offset>0:
            _cov = limix.CFixedCF(SP.eye(self.P))
            _cov.setParams(SP.array([SP.sqrt(offset)]))
            _cov.setParamMask(SP.zeros(1))
            cov.addCovariance(_cov)
        self.offset.append(offset)

        self.covar_type.append(covar_type)
        self.diag.append(diag)

        self.vd.addTerm(cov,K)
        if Ks!=None: self.setKstar(self.n_terms,Ks)
        self.n_terms+=1

        self.gp      = None
        self.init    = False
        self.fast    = False
        self.optimum = None

        self.cache['Sigma']   = None
        self.cache['Hessian'] = None
        self.cache['Lparams'] = None
        self.cache['paramsST']= None
add multi trait random effects term.
        The inter-trait covariance is parametrized by covar_type, where parameters are optimized.
        
        Args:
            K:      Individual-individual (Intra-Trait) Covariance Matrix [N, N]
                    (K is normalised in the C++ code such that K.trace()=N)
            covar_type: type of covaraince to use. Default 'freeform'. possible values are 
                            'freeform': free form optimization, 
                            'fixed': use a fixed matrix specified in covar_K0,
                            'diag': optimize a diagonal matrix, 
                            'lowrank': optimize a low rank matrix. The rank of the lowrank part is specified in the variable rank,
                            'lowrank_id': optimize a low rank matrix plus the weight of a constant diagonal matrix. The rank of the lowrank part is specified in the variable rank, 
                            'lowrank_diag': optimize a low rank matrix plus a free diagonal matrix. The rank of the lowrank part is specified in the variable rank, 
                            'block': optimize the weight of a constant P x P block matrix of ones,
                            'block_id': optimize the weight of a constant P x P block matrix of ones plus the weight of a constant diagonal matrix,
                            'block_diag': optimize the weight of a constant P x P block matrix of ones plus a free diagonal matrix,                            
            is_noise:   Boolean indicator specifying if the matrix is homoscedastic noise (weighted identity covariance) (default False)
            normalize:  Boolean indicator specifying if K is normalized such that K.trace()=N.
            Ks:			NxNtest cross covariance for predictions
            offset:		diagonal contribution added to trait-to-trait covariance matrices for regularization
            rank:       rank of a possible lowrank component (default 1)
            covar_K0:   PxP matrix for the (predefined) trait-to-trait covariance matrix if fixed type is used