Running
#######

.. role:: bash(code)
  :language: bash

Restarting a crashed simulation
===============================

If the :bash:`reload_results` input argument is set to ``true``, an interrupted
``response-functions`` mode simulation will not recompute already computed
response function simulation and resume where it was interrupted when restarted.

The :doc:`observation operator plugin</documentation/plugins/obsoperators/standard>`
:bash:`autorestart` input argument can also be used to resume the simulation
that was interrupted.


Exploit results
###############

The following code examples can be used to read and exploit the response functions results:

:math:`\mathbf{H}` matrix
=========================

Reading the :math:`\mathbf{H}` matrix (with the control vector coordinates) from
the observation vector decomposition.

.. note:: 
  The response function simulation need to be run with the input argument
  :bash:`dump_obsvect_decompostion` set to :bash:`true` to the CIF produces
  the necessary outputs.

.. code-block:: python

  import os

  import numpy as np
  import pandas as pd
  import xarray as xr
  from pycif.utils.sparse_array import to_dense_dataset

  def read_decomp(path, region_names=None, use_sparse_dtype=False):
      with xr.open_dataset(path) as ds:
          if 'sparse_sim_tl' in ds:
              ds = to_dense_dataset(ds)

          ds = ds.isel(vresol=0).drop_vars('vresol')  # Dropping vertical levels

          # Replacing hresol with region names
          if region_names is not None and 'region_id' in ds:
              ds['hresol'] = (['hresol'], region_names[ds.region_id.values])

          decomp = ds.sim_tl.to_series().unstack(['date', 'hresol'])

      if use_sparse_dtype:
          decomp = decomp.astype(pd.SparseDtype("float", 0.0))

      return decomp

  def read_h_matrix(workdir, region_names=None):
      h_dict = {}
      dirname = os.path.join(workdir, "base_functions/decomposition")

      for component in os.listdir(dirname):
          comp_dirname = os.path.join(dirname, component)        
            for parameter in os.listdir(comp_dirname):
                print(f"reading H matrix for {component}/{parameter}")
                path = os.path.join(comp_dirname, parameter, "obsvect_decomp.nc")
                h_dict[(component, parameter)] = read_decomp(path, region_names, True)

      df = pd.concat(
          h_dict.values(),
          keys=h_dict.keys(),
          names=['component', 'parameter'],
          axis='columns',
      )

      return df


  with xr.open_dataset("region_mask.nc") as ds_mask:
      region_names = ds_mask.regions_name.values

  h_matrix = read_h_matrix("path/to/workdir", region_names)


Variational inv
###############

* **TODO**: Output
* **TODO**: Rerun with another B, R