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

import copy
import os
import glob
import re
import pandas as pd
import numpy as np
import scipy
import pickle
from logging import info
from ....utils.datastores.dump import read_datastore
from .basefunctions import compute_base_functions, check_base_functions
from .build_H import build_H




def execute(self, **kwargs):
    """Performs an analytical inversions, including computation of base
    functions"""
    
    # Working directory
    workdir = self.workdir

    # Control vector
    controlvect = self.controlvect
    xb_ref = copy.deepcopy(controlvect.xb)
    
    # Observation operator
    obsoper = self.obsoperator
    
    # Observation vector
    obsvect = self.obsvect

    # Simulation window
    datei = self.datei
    datef = self.datef

    # Some verbose
    towrite = """
        Running an analytical inversion with the following modules:
            Observation operator: {}
            Model: {}
        """.format(
        obsoper.plugin.name, obsoper.model.plugin.name
    )
    info(towrite)
    
    # Check if H already computed
    base_dir = "{}/base_functions/H_matrix/".format(workdir)
    H_file = os.path.join(base_dir, "H.pickle")
    try:
        with open(H_file, "rb") as f:
            harray = pickle.load(f)
    
    except IOError:
        # Check base functions and re-compute if necessary
        missing_base_functions = \
            check_base_functions(base_dir, controlvect.dim)
        compute_base_functions(self, controlvect,
                               missing_base_functions,
                               self.dryrun, self.sequential)
        
        # Re-construct the H operator
        harray = build_H(controlvect, obsvect, base_dir)
    
    # Do analytical inversion
    # Xa = Xb + K(y-Hx)
    # K = BH*(R + HBH*)-1
    # Pa = B - KHB
    controlvect.xb[:] = xb_ref[:]
    dy = obsvect.yobs - np.dot(harray, controlvect.xb)
    if ~self.resp_func_only:
        bfull = controlvect.build_b(controlvect)
        rfull = obsvect.build_r(obsvect)

    # Crop matrices with is_obsvect
    mask_obsvect = obsvect.obsvect_mask
    dy = dy[mask_obsvect]
    harray = harray[mask_obsvect]
    list_ind = np.argwhere(mask_obsvect).flatten()
    if self.resp_func_only: 
      return
      
    
    rfull = rfull[np.ix_(list_ind, list_ind)]

    # Compute the inversion
    srm = scipy.linalg.inv(rfull + harray.dot(bfull.dot(harray.T)))
    kmatrix = bfull.dot(harray.T.dot(srm))

    xa = xb_ref + kmatrix.dot(dy)
    pa = bfull - kmatrix.dot(harray.dot(bfull))
    
    # Saving posterior vector and uncertainties
    controlvect.x = xa
    controlvect.pa = pa
    
    dir_control = "{}/controlvect/".format(workdir)
    file_xa = "{}/controlvect_post.pickle".format(dir_control)
    controlvect.dump(file_xa, to_netcdf=True, dir_netcdf=dir_control)
    
    # Saving prior simulations
    obsvect.ysim = np.dot(harray, controlvect.xb)
    dir_dump = "{}/obsvect_prior/".format(workdir)
    obsvect.dump(dir_dump)
    
    # Saving posterior simulations
    obsvect.ysim = np.dot(harray, xa)
    dir_dump = "{}/obsvect_posterior/".format(workdir)
    obsvect.dump(dir_dump)
    
    # Return optimized control vector and observation vector
    return controlvect, obsvect