<jupyter_start><jupyter_text>Part 5: Application using the Discretized Misfit calculationBy now, you should have looked through [Part 1](IntroductionToMetric.ipynb), [Part 2](IntroductionToResidual.ipynb) of the introductory notebook series, and [Part 3](OtherIO_options.ipynb) of the introductory notebook series. These introduced the umami `Metric` and `Residual` classes. You should have also worked through [Part 4](ExampleApplication.ipynb), which provides an example application of umami.  ScopeSimilar to [Part 4](ExampleApplication.ipynb), this application will use umami alongside the [terrainbento](https://terrainbento.readthedocs.io/en/latest/) package. As in the prior example, we will define a "synthetic truth" model evaluation with a specific set of input parameters, and then do a grid search letting some of those parameters vary. In this way we will explore which statistics for model-data comparison do best at identifying the "true" parameters. This application focuses on the least intuitive of the umami calculations: the [`discretized_misfit`](https://umami.readthedocs.io/en/latest/umami.calculations.residual.discretized_misfit.html). If you have comments or questions about the notebooks, the best place to get help is through [GitHub Issues](https://github.com/TerrainBento/umami/issues).To begin, we import necessary modules.<jupyter_code>import warnings
warnings.filterwarnings('ignore')

from io import StringIO
from itertools import product
from tqdm import tqdm

import numpy as np

import pandas as pd

import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline

from plotnine import *

import holoviews as hv
hv.notebook_extension('matplotlib')

from landlab import imshow_grid
from terrainbento import Basic
from umami import Metric, Residual<jupyter_output><empty_output><jupyter_text>We begin by defining an input string that defines the terrainbento model.It evolves topography using stream power and linear diffusion. In our application it has a boundary condition of uniform uplift across the core of the model grid. It thus has the following governing equation:$\frac{\partial \eta}{\partial t} = U - KQ^{1/2}S + D\nabla^2 \eta$where $K$ and $D$ are constants, $Q$ is discharge, $S$ is local slope, $U$ is the uplift rate, and $\eta$ is the topography. See the [model Basic documentation](https://terrainbento.readthedocs.io/en/latest/source/terrainbento.derived_models.model_basic.html) for additional information. In this input file we also indicate that the model will run with timesteps of 500 yr and the model grid will have a shape of (50, 80), with grid cell spacing of 100 m. For this application, the model initial condition is for core nodes to have random noise added.A few places in the input file have curly braces around a name. * Two inputs parameters have curly brackets around them: `{duration}` and `{water_erodibility}`. These inputs will be modified using [`str.format`](https://docs.python.org/3/library/stdtypes.htmlstr.format) to set the "truth" model run and to vary the parameters in a grid search numerical experiment. * We also format the `{name}` of output files in order to prevent Windows file permissions errors.<jupyter_code>spec_string = """
# Create the Clock.
clock:
    start: 0
    step: 500
    stop: {duration}

# Create the Grid
grid: 
    RasterModelGrid: 
        - [100, 120]
        - xy_spacing: 50
        - fields: 
            node: 
                topographic__elevation:
                    random:
                        where: CORE_NODE

                        
# Set up Boundary Handlers
boundary_handlers: 
    NotCoreNodeBaselevelHandler: 
        modify_core_nodes: True
        lowering_rate: -{lowering_rate}

# Parameters that control output.
output_interval: 1e3
save_first_timestep: True
output_prefix: 
    disc_resid.{name}.
fields: 
    - topographic__elevation

# Parameters that control process and rates.
water_erodibility: {water_erodibility}
m_sp: 0.5
n_sp: 1.0
regolith_transport_parameter: 0.1
"""<jupyter_output><empty_output><jupyter_text>Next we instantiate the "truth" model and run it.<jupyter_code>truth_duration = 3e4
truth_water_erodibility = 0.0005

lowering_rate = 100 / truth_duration

truth_params = StringIO(
    spec_string.format(duration=truth_duration,
                       water_erodibility=truth_water_erodibility,
                       lowering_rate=lowering_rate,
                       name="truth"))
np.random.seed(42)
truth = Basic.from_file(truth_params)
truth.run()<jupyter_output><empty_output><jupyter_text>The [holoviews](https://holoviews.org) package provides capabilities to visualize the model run.<jupyter_code>ds = truth.to_xarray_dataset(time_unit='years', space_unit='meters')
hvds_topo = hv.Dataset(ds.topographic__elevation)
topo = hvds_topo.to(hv.Image, ['x', 'y'],
                    label='Truth').options(interpolation='bilinear',
                                           cmap='viridis',
                                           colorbar=True)
topo<jupyter_output><empty_output><jupyter_text>As the center nodes are uplifted, a series of ridges and drainage basins form. This model has not yet reached "topographic steady state", in which $\frac{\partial \eta}{\partial t}>\epsilon$ (and $\epsilon$ is small) everywhere. Before moving on, we close the xarray dataset and remove the output netcdf files.<jupyter_code>ds.close()
truth.remove_output_netcdfs()<jupyter_output><empty_output><jupyter_text>Step 2: Define the basis for model-data comparisonThe discretized_misfit takes the following parameters:    Parameters    ----------    model_grid : Landlab model grid    data_grid : Landlab model grid    name : str    misfit_field : str    field_1 : str    field_2 : str    field_1_percentile_edges : list    field_2_percentile_edges : list    This calculation first classifies each grid cell in the landscape into categories based on `field_1`, `field_2` and the percentile edges for each (using the data grid). This results in a set of categories, which may or may not be contiguous in space. For each category, the sum of squared residuals is calculated based on the misfit_field.Since this calculation returns one value for each category, rather than one value in total, a `name` must be provided. This is a string that will be formatted with the values for `{field_1_level}` and `{field_2_level}`. The output is an ordered dictionary with name as the keys, and the sum of squares misfit as the values.The following is the input file (as string) needed to specify a `discretized_misfit` in which the domain is discretized based on `channel__chi_index` (three percentile levels defined by `[0, 30, 60, 100]`), and `topographic__elevation` (two percentile levels defined by `[0, 50, 100]`). Within each of the six category domains, a misfit is made based on the field `topographic__elevation`. Below we will show a plot indicating where each category is located, but first we will specify the numerical experiment.<jupyter_code>residual_string = """
dm:
    _func: discretized_misfit
    name: chi_{field_1_level}.z_{field_2_level}
    misfit_field: topographic__elevation
    field_1: channel__chi_index
    field_2: topographic__elevation
    field_1_percentile_edges:
        - 0
        - 30
        - 60
        - 100
    field_2_percentile_edges:
        - 0
        - 50
        - 100
"""<jupyter_output><empty_output><jupyter_text>Step 3: Create and run the grid search experimentWe will use a grid search to highlight how the misfit values in the `discretized_residual` calculated by umami vary across parameter space. We consider values for `duration` between $10^{3}$ and $10^{5}$ and values for $K$ (`water_erodibility`) between $10^{-4}$ and $10^{-2}$.With a resolution of 10, we evaluate $10^2=100$ simulations. Feel free to change the resolution value, though note that it will impact the run time of this notebook.<jupyter_code>resolution = 10
durations = np.logspace(3, 5, num=resolution)
water_erodibilitys = np.logspace(-4, -2, num=resolution)<jupyter_output><empty_output><jupyter_text>We evaluate each pair of duration and water erodability and save the model output as a dictionary. With the line     np.random.seed(42)commented out, each evaluation uses a different random seed. Feel free to uncomment this line to see how the results change if the *exact same* random seed is used for each model integration.<jupyter_code>out = {}
for i, (duration,
        water_erodibility) in enumerate(tqdm(product(durations,
                                                water_erodibilitys))):
    lowering_rate = 100 / duration
    test_params = StringIO(
        spec_string.format(duration=duration,
                           water_erodibility=water_erodibility,
                           lowering_rate=lowering_rate,
                           name=i))
    #np.random.seed(42)
    test = Basic.from_file(test_params)
    test.run()

    test.remove_output_netcdfs()

    residual = Residual(test.grid,
                        truth.grid,
                        chi_finder_kwds={"min_drainage_area": 1000})
    residual.add_from_file(StringIO(residual_string))
    residual.calculate()

    values = {name: residual.value(name) for name in residual.names}
    out[(duration, water_erodibility)] = values<jupyter_output><empty_output><jupyter_text>Before looking at the results, let's inspect the residual class. The property `residual.names` has a name for each of the six categories. The temporary strings `{field_1_level}` and `{field_2_level}` have been replaced with actual levels.<jupyter_code>residual.names<jupyter_output><empty_output><jupyter_text>The property `residual.values` has one value for each of the six categories.<jupyter_code>residual.values<jupyter_output><empty_output><jupyter_text>We can plot the category using the `residual.category` property. Here each category gets its own panel. For example, the leftmost column represents the cells with `channel__chi_index` values in the lower 30%. The upper left panel is those that have the lower half of elevation values.<jupyter_code>fig, axes = plt.subplots(2, 3, dpi=150)

for i, name in enumerate(residual.names):
    col, row = np.unravel_index(i, (3,2))
    plt.sca(axes[row, col])
    imshow_grid(truth.grid, residual.category==(i+1), cmap="cividis", allow_colorbar=False)
    plt.title(name)
    plt.xticks([])
    plt.yticks([])
    plt.xlabel(None)
    plt.ylabel(None)

plt.tight_layout()<jupyter_output><empty_output><jupyter_text>Next we compile the output into a pandas dataframe and inspect.<jupyter_code>df = pd.DataFrame.from_dict(out, orient="index")
df.index.names = ["duration", "water_erodibility"]
df.head()<jupyter_output><empty_output><jupyter_text>Similar to the earlier notebook, we melt this dataframe in order to plot it. We also rename the column "value" to "sum_of_squared_residuals", as this is what `discretized_misfit` calculates.<jupyter_code>df_melt = df.reset_index().melt(id_vars=["duration", "water_erodibility"])
df_melt = df_melt.rename(columns={"value": "sum_of_squared_residuals"})
df_melt.head()<jupyter_output><empty_output><jupyter_text>We also use the [string.split](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.Series.str.split.html) pandas function in order to turn our variable name into two columns<jupyter_code>variable_split = df_melt["variable"].str.split(".", n=1, expand=True)
df_melt["chi"] = variable_split[0]
df_melt["elev"] = variable_split[1]
df_melt.head()<jupyter_output><empty_output><jupyter_text>Finally we plot the squared residual as a function:<jupyter_code>p1 = (ggplot(df_melt,
             aes(x="duration", y="water_erodibility",
                 fill="sum_of_squared_residuals")) + geom_tile() +
      geom_point(aes(x=truth_duration, y=truth_water_erodibility)) +
      scale_fill_continuous(limits=[1, 100], trans="log10") +
      facet_grid("elev~chi") + theme_bw() + scale_x_log10() + scale_y_log10() +
      coord_equal())
print(p1)<jupyter_output><empty_output>