docs: def calculate(cls, calc, formula_reg, data_reg, out_reg,
                    timestep=None, idx=None):
          """
          Execute calculation
  
          :param calc: calculation, with formula, args and return keys
          :type calc: dict
          :param formula_reg: Registry of formulas.
          :type formula_reg: :class:`~simkit.core.FormulaRegistry`
          :param data_reg: Data registry.
          :type data_reg: :class:`~simkit.core.data_sources.DataRegistry`
          :param out_reg: Outputs registry.
          :type out_reg: :class:`~simkit.core.outputs.OutputRegistry`
          :param timestep: simulation interval length [time], default is ``None``
          :param idx: interval index, default is ``None``
          :type idx: int
          """
          # get the formula-key from each static calc
          formula = calc['formula']  # name of formula in calculation
          func = formula_reg[formula]  # formula function object
          fargs = formula_reg.args.get(formula, [])  # formula arguments
          constants = formula_reg.isconstant.get(formula)  # constant args
          # formula arguments that are not constant
          vargs = [] if constants is None else [a for a in fargs if a not in constants]
          args = calc['args']  # calculation arguments
          # separate data and output arguments
          datargs, outargs = args.get('data', {}), args.get('outputs', {})
          data = index_registry(datargs, data_reg, timestep, idx)
          outputs = index_registry(outargs, out_reg, timestep, idx)
          kwargs = dict(data, **outputs)  # combined data and output args
          args = [kwargs.pop(a) for a in fargs if a in kwargs]
          returns = calc['returns']  # return arguments
          # if constants is None then the covariance should also be None
          # TODO: except other values, eg: "all" to indicate no covariance
          if constants is None:
              cov = None  # do not propagate uncertainty
          else:
              # get covariance matrix
              cov = cls.get_covariance(datargs, outargs, vargs,
                                       data_reg.variance, out_reg.variance)
              # update kwargs with covariance if it exists
              kwargs['__covariance__'] = cov
          retval = func(*args, **kwargs)  # calculate function
          # update output registry with covariance and jacobian
          if cov is not None:
              # split uncertainty and jacobian from return values
              cov, jac = retval[-2:]
              retval = retval[:-2]
              # scale covariance
              scale = np.asarray(
                  [1 / r.m if isinstance(r, UREG.Quantity) else 1 / r
                   for r in retval]
              )  # use magnitudes if quantities
              cov = (np.swapaxes((cov.T * scale), 0, 1) * scale).T
              nret = len(retval)  # number of return output
              for m in xrange(nret):
                  a = returns[m]  # name in output registry
                  out_reg.variance[a] = {}
                  out_reg.uncertainty[a] = {}
                  out_reg.jacobian[a] = {}
                  for n in xrange(nret):
                      b = returns[n]
                      out_reg.variance[a][b] = cov[:, m, n]
                      if a == b:
                          unc = np.sqrt(cov[:, m, n]) * 100 * UREG.percent
                          out_reg.uncertainty[a][b] = unc
                  for n in xrange(len(vargs)):
                      b = vargs[n]
                      try:
                          b = datargs[b]
                      except (KeyError, TypeError):
                          b = outargs[b]
                      out_reg.jacobian[a][b] = jac[:, m, n]
                  LOGGER.debug('%s cov:\n%r', a, out_reg.variance[a])
                  LOGGER.debug('%s jac:\n%r', a, out_reg.jacobian[a])
                  LOGGER.debug('%s unc:\n%r', a, out_reg.uncertainty[a])
          # if there's only one return value, squeeze out extra dimensions
          if len(retval) == 1:
              retval = retval[0]
          # put return values into output registry
          if len(returns) > 1:
              # more than one return, zip them up
              if idx is None:
                  out_reg.update(zip(returns, retval))
              else:
                  for k, v in zip(returns, retval):
                      out_reg[k][idx] = v
          else:
              # only one return, get it by index at 0
              if idx is None:
                  out_reg[returns[0]] = retval
              else:
                  out_reg[returns[0]][idx] = retval