docs: def get_covariance(datargs, outargs, vargs, datvar, outvar):
          """
          Get covariance matrix.
  
          :param datargs: data arguments
          :param outargs: output arguments
          :param vargs: variable arguments
          :param datvar: variance of data arguments
          :param outvar: variance of output arguments
          :return: covariance
          """
          # number of formula arguments that are not constant
          argn = len(vargs)
          # number of observations must be the same for all vargs
          nobs = 1
          c = []
          # FIXME: can just loop ver varg, don't need indices I think, do we?
          for m in xrange(argn):
              a = vargs[m]  # get the variable formula arg in vargs at idx=m
              try:
                  a = datargs[a]  # get the calculation data arg
              except (KeyError, TypeError):
                  a = outargs[a]  # get the calculation output arg
                  if not isinstance(a, basestring):
                      # calculation arg might be sequence (arg, idx, [unit])
                      a = a[0]  # if a is a sequence, get just the arg from a[0]
                      LOGGER.debug('using output variance key: %r', a)
                  avar = outvar[a]  # get variance from output registry
              else:
                  if not isinstance(a, basestring):
                      # calculation arg might be sequence (arg, idx, [unit])
                      a = a[0]  # if a is a sequence, get just the arg from a[0]
                      LOGGER.debug('using data variance key: %r', a)
                  avar = datvar[a]  # get variance from data registry
              d = []
              # avar is a dictionary with the variance of "a" vs all other vargs
              for n in xrange(argn):
                  # FIXME: just get all of the calculation args one time
                  b = vargs[n]  # get the variable formula arg in vargs at idx=n
                  try:
                      b = datargs[b]  # get the calculation data arg
                  except (KeyError, TypeError):
                      b = outargs[b]  # get variance from output registry
                  if not isinstance(b, basestring):
                      # calculation arg might be sequence (arg, idx, [unit])
                      b = b[0]  # if a is a sequence, get just the arg from b[0]
                      LOGGER.debug('using variance key: %r', b)
                  d.append(avar.get(b, 0.0))  # add covariance to sequence
                  # figure out number of observations from longest covariance
                  # only works if nobs is in one of the covariance
                  # fails if nobs > 1, but covariance shape doesn't have nobs!!!
                  # eg: if variance for data is uniform for all observations!!!
                  try:
                      nobs = max(nobs, len(d[-1]))
                  except (TypeError, ValueError):
                      LOGGER.debug('c of %s vs %s = %g', a, b, d[-1])
              LOGGER.debug('d:\n%r', d)
              c.append(d)
          # covariance matrix is initially zeros
          cov = np.zeros((nobs, argn, argn))
          # loop over arguments in both directions, fill in covariance
          for m in xrange(argn):
              d = c.pop()
              LOGGER.debug('pop row %d:\n%r', argn-1-m, d)
              for n in xrange(argn):
                  LOGGER.debug('pop col %d:\n%r', argn - 1 - n, d[-1])
                  cov[:, argn-1-m, argn-1-n] = d.pop()
          if nobs == 1:
              cov = cov.squeeze()  # squeeze out any extra dimensions
          LOGGER.debug('covariance:\n%r', cov)
          return cov