Source code for pycif.plugins.transforms.complex.isotopes.forward

import numpy as np
import pandas as pd
import copy




def forward(
    transf,
    data,
    mapper,
    ddi,
    ddf,
    mode,
    runsubdir,
    workdir,
    **kwargs
):
    """Compute two-isotopologue concentrations from a δ-value and total concentration.

    Given one isotopic signature δ and total concentration :math:`C_{tot}`:

    .. math::

        a = (1 + \\delta/1000) \\cdot R_{std}

        [iso_1] = C_{tot} / (1 + a), \\quad
        [iso_2] = C_{tot} \\cdot a / (1 + a)

    Saves ``signature`` and ``spec_data`` on ``transf.fwd_data`` for use
    by the adjoint.  Applies an optional mass correction (``unit = 'mass'``).

    Args:
        transf (Plugin): isotopes plugin instance (carries ``parameters_in``
            with ``standard``, ``parameters_out`` with ``iso_mass`` /
            ``spec_mass``, ``unit``, and ``component``).
        data: datastore object carrying ``.datastore`` dict.
        mapper (dict): transform mapper.
        ddi (datetime): sub-simulation start date.
        ddf (datetime): sub-simulation end date.
        mode (str): ``'fwd'`` or ``'tl'``.
        runsubdir (str): unused.
        workdir (str): unused.
        **kwargs: unused.

    Returns:
        data: updated datastore object.
    """

    xmod = data.datastore

    inputs = transf.parameters_in.names
    outputs = transf.parameters_out.names
    r_std = transf.parameters_in.standard
    iso_mass = transf.parameters_out.iso_mass
    spec_mass = transf.parameters_out.spec_mass
    unit = transf.unit

    in_types = transf.component
    out_types = 2 * [transf.component[1]]

    keys = ['spec']
    if mode == 'tl':
        keys.append('incr')

    for (trid, trid_out) in zip(mapper["inputs"], mapper["outputs"]):
        xmod[trid_out] = {k: xmod[trid][k] for k in xmod[trid]}
        xmod[trid_out]["spec"] = copy.deepcopy(xmod[trid]["spec"])
        if mode == 'tl':
            xmod[trid_out]["incr"] = copy.deepcopy(xmod[trid]["incr"])

    signature = xmod[(in_types[0], inputs[0])]['spec']
    spec_data = xmod[(in_types[1], inputs[1])]['spec']

    # Save signature and data for later use by adjoint
    transf.fwd_data = {
        'signature': signature.values,
        'spec_data': spec_data.values}

    a_factor = (1 + signature / 1000) * r_std

    xmod[(out_types[0], outputs[0])]['spec'] = \
        spec_data / (1 + a_factor)
    xmod[(out_types[1], outputs[1])]['spec'] = \
        spec_data * a_factor / (1 + a_factor)

    # Applying tangent linear
    if mode == 'tl':
        spec_data_tl = xmod[(in_types[1], inputs[1])]['incr']
        signature_tl = xmod[(in_types[0], inputs[0])]['incr']

        xmod[(out_types[0], outputs[0])]['incr'] = \
            spec_data_tl / (1 + a_factor) \
            - spec_data * r_std * signature_tl \
            / (1000 * (1 + a_factor) ** 2)
        xmod[(out_types[1], outputs[1])]['incr'] = \
            spec_data_tl * a_factor / (1 + a_factor) \
            + spec_data * r_std * signature_tl \
            / (1000 * (1 + a_factor) ** 2)

    # Mass correction if units are in mass
    if unit == 'mass':
        xmod[(out_types[0], outputs[0])]['spec'] *= \
            (iso_mass[0] / spec_mass)
        xmod[(out_types[1], outputs[1])]['spec'] *= \
            (iso_mass[1] / spec_mass)

        if mode == 'tl':
            xmod[(out_types[0], outputs[0])]['incr'] *= \
                (iso_mass[0] / spec_mass)
            xmod[(out_types[1], outputs[1])]['incr'] *= \
                (iso_mass[1] / spec_mass)

    # del xmod[(in_types[0], inputs[0])]
    # del xmod[(in_types[1], inputs[1])]

    data.datastore = xmod
    return data