Source code for pycif.plugins.modes.analytic.build_H

import pandas as pd
import os
import pickle
import numpy as np
from ....utils.datastores.dump import read_datastore




def build_H(controlvect, obsvect, base_dir):
    """Assemble the full observation operator matrix H from pre-computed base functions.

    Reads each base-function obs-vector from *base_dir* and fills the
    corresponding column of the H matrix.  Dumps the result as a pickle file
    (``H.pickle``) in *base_dir*.

    Args:
        controlvect: control vector plugin providing ``dim`` and component
            metadata.
        obsvect: observation vector plugin providing ``dim`` and per-tracer
            pointer / dimension attributes.
        base_dir (str): directory containing ``obsvect_NNNN`` sub-directories
            (one per control-vector dimension) and where ``H.pickle`` will be
            written.

    Returns:
        np.ndarray: H matrix of shape ``(obs_dim, control_dim)``.

    Raises:
        Exception: if the datavect has no ``components`` attribute (no
            observations configured).
    """
    datavect = obsvect.datavect
    
    # If no definition is specified for the control vector in the Yaml,
    # return empty control vector
    if not hasattr(datavect, "components"):
        raise Exception("Computing H operator with no observations")
    
    # Else, carry on initializing
    harray = np.zeros((obsvect.dim, controlvect.dim))
    jacobian = {}
    components = datavect.components
    for comp in components.attributes:
        component = getattr(components, comp)
        
        # Skip if component does not have parameters
        if not hasattr(component, "parameters"):
            continue
        
        for trcr in component.parameters.attributes:
            tracer = getattr(component.parameters, trcr)
            
            # Skip tracers that are not control variables
            tracer.isobs = (
                    comp in obsvect.default_obstypes
                    or tracer.isobs
            )
            
            if not tracer.isobs:
                continue
            
            ds = tracer.datastore
            hloc = pd.DataFrame(columns=range(controlvect.dim),
                                index=ds.index)
            for idim in range(controlvect.dim):
                monitor_file = \
                    "{}/obsvect_{:04d}/{}/{}/monitor.nc"\
                        .format(base_dir, idim, comp, trcr)
                data = read_datastore(monitor_file)
                hloc.loc[:, idim] = data.loc[:, ("maindata", "sim")].values
                
            jacobian[(comp, trcr)] = hloc
            harray[tracer.ypointer: tracer.ypointer + tracer.dim, :] = \
                hloc

    # Dumping H
    H_file = os.path.join(base_dir, "H.pickle")
    with open(H_file, "wb")as f:
        pickle.dump(harray, f, pickle.HIGHEST_PROTOCOL)
    
    return harray