code: def getCovariance(self,normalize=True,i0=None,i1=None,pos0=None,pos1=None,chrom=None,center=True,unit=True,pos_cum0=None,pos_cum1=None,blocksize=None,X=None,**kw_args):
          
          if X is not None:
              K=X.dot(X.T)
              Nsnp=X.shape[1]
          else:
              if (i0 is None) and (i1 is None) and ((pos0 is not None) & (pos1 is not None) & (chrom is not None)) or ((pos_cum0 is not None) & (pos_cum1 is not None)):
                  i0,i1=self.getGenoIndex(pos0=pos0,pos1=pos1,chrom=chrom,pos_cum0=pos_cum0,pos_cum1=pose_cum1)
  
              [N,M]=self.genoM.shape
              if blocksize is None:
                  blocksize=M
              if i0 is None:
                  i0=0
              if i1 is None:
                  i1=M
              nread = i0
              K=None
              Nsnp=i1-i0
              while nread<i1:
                  thisblock=min(blocksize,i1-nread)
                  X=self.getGenotypes(i0=nread,i1=(nread+thisblock),center=center,unit=unit,**kw_args)    
                  if K is None:
                      K=X.dot(X.T)
                  else:
                      K+=X.dot(X.T)
                  nread+=thisblock
          if normalize:
              K/=(K.diagonal().mean())
          else:#divide by number of SNPs in K
              K/=Nsnp
          return K
docs: calculate the empirical genotype covariance in a region