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
def calculate ( cls , calc , formula_reg , data_reg , out_reg , timestep = None , idx = None ) : formula = calc [ 'formula' ] func = formula_reg [ formula ] fargs = formula_reg . args . get ( formula , [ ] ) constants = formula_reg . isconstant . get ( formula ) vargs = [ ] if constants is None else [ a for a in fargs if a not in constants ] args = calc [ 'args' ] 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 ) args = [ kwargs . pop ( a ) for a in fargs if a in kwargs ] returns = calc [ 'returns' ] if constants is None : cov = None else : cov = cls . get_covariance ( datargs , outargs , vargs , data_reg . variance , out_reg . variance ) kwargs [ '__covariance__' ] = cov retval = func ( * args , ** kwargs ) if cov is not None : cov , jac = retval [ - 2 : ] retval = retval [ : - 2 ] scale = np . asarray ( [ 1 / r . m if isinstance ( r , UREG . Quantity ) else 1 / r for r in retval ] ) cov = ( np . swapaxes ( ( cov . T * scale ) , 0 , 1 ) * scale ) . T nret = len ( retval ) for m in xrange ( nret ) : a = returns [ m ] 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 len ( retval ) == 1 : retval = retval [ 0 ] if len ( returns ) > 1 : if idx is None : out_reg . update ( zip ( returns , retval ) ) else : for k , v in zip ( returns , retval ) : out_reg [ k ] [ idx ] = v else : if idx is None : out_reg [ returns [ 0 ] ] = retval else : out_reg [ returns [ 0 ] ] [ idx ] = retval
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