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

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




def build_Hx(obsvect, segment_dir, nsample, ddi, ddf):
    """Builds the observation operator matrix based on the simulated 
    values of each member.

    Args:
        obsvect (Plugin): obsvect Plugin
        segment_dir (str): path to the segment directory
        nsample (int): number of members in the ensemble (besides prior and mean members)
        ddi (datetime): initial date of the segment
        ddf (datetime): end date of the segment

    Returns:
        harray (np.array): matrix representing the observation operator for each member
    """

    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, nsample + 3))
    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 observation variables
            if not tracer.isobs:
                continue

            # Load and crop the datastore to the current inversion window
            ds = tracer.datastore

            # Select the part of the datastore matching the segment
            crop = obsvect.crop_monitor(ds, ddi, ddf)
            mask_seg = ds.index.isin(crop.index)

            # Load temporary monitor
            hloc = pd.DataFrame(columns=range(nsample + 3),
                                index=ds.index)
            for isample in range(nsample + 3):
                monitor_file = \
                    "{}/H_matrix/obsvect_{:04d}/{}/{}/monitor.nc"\
                        .format(segment_dir, isample, comp, trcr)
                data = read_datastore(monitor_file)
                hloc.loc[mask_seg, isample] = data.loc[:, ("maindata", "sim")].values

            harray[tracer.ypointer: tracer.ypointer + tracer.dim, :] = hloc

    return harray