docs: 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