""""""
from collections import OrderedDict
from copy import deepcopy
import numpy as np
import yaml
from numpy.testing import assert_array_equal
import umami.calculations.metric as metric_calcs
import umami.calculations.residual as residual_calcs
from landlab import RasterModelGrid, create_grid
from umami.metric import Metric
from umami.utils.create_landlab_components import _create_landlab_components
from umami.utils.io import _read_input, _write_output
from umami.utils.validate import _validate_fields, _validate_func
_VALID_FUNCS = {}
_VALID_FUNCS.update(residual_calcs.__dict__)
_VALID_FUNCS.update(metric_calcs.__dict__)
[docs]class Residual(object):
"""Create a ``Residual`` class based on a model and data Landlab grid."""
_required_fields = ["topographic__elevation"]
[docs] def __init__(
self,
model,
data,
flow_accumulator_kwds=None,
chi_finder_kwds=None,
residuals=None,
):
"""
Parameters
----------
model : Landlab model grid
data : Landlab model grid
flow_accumulator_kwds : dict
Parameters to pass to the Landlab ``FlowAccumulator`` to specify
flow direction and accumulation.
chi_finder_kwds : dict
Parameters to pass to the Landlab ``ChiFinder`` to specify optional
arguments. `
residuals : dict
A dictionary of desired residuals to calculate. See examples for
required format.
Examples
--------
>>> import numpy as np
>>> from io import StringIO
>>> from landlab import RasterModelGrid
>>> from umami import Residual
>>> np.random.seed(42)
>>> model = RasterModelGrid((10, 10))
>>> z_model = model.add_zeros("node", "topographic__elevation")
>>> z_model += model.x_of_node + model.y_of_node
>>> data = RasterModelGrid((10, 10))
>>> z_data = data.add_zeros("node", "topographic__elevation")
>>> z_data += data.x_of_node + data.y_of_node
>>> z_data[data.core_nodes] += np.random.random(data.core_nodes.shape)
>>> file_like=StringIO('''
... me:
... _func: aggregate
... method: mean
... field: topographic__elevation
... ep10:
... _func: aggregate
... method: percentile
... field: topographic__elevation
... q: 10
... oid1_mean:
... _func: watershed_aggregation
... field: topographic__elevation
... method: mean
... outlet_id: 1
... sn1:
... _func: count_equal
... field: drainage_area
... value: 1
... ''')
>>> residual = Residual(model, data)
>>> residual.add_from_file(file_like)
>>> residual.names
['me', 'ep10', 'oid1_mean', 'sn1']
>>> residual.calculate()
>>> np.testing.assert_array_almost_equal(
... residual.value("me"),
... -0.467,
... decimal=3)
>>> np.testing.assert_array_almost_equal(
... residual.values,
... np.array([ -0.467, -0.151, 3.313, -18. ]),
... decimal=3)
"""
# assert that the model grids have the same x_of_node and y_of_node.
assert_array_equal(data.x_of_node, model.x_of_node)
assert_array_equal(data.y_of_node, model.y_of_node)
assert_array_equal(data.core_nodes, model.core_nodes)
self._data_grid = data
self._model_grid = model
# verify that apppropriate fields are present.
for field in self._required_fields:
if (field not in model.at_node) or (field not in data.at_node):
msg = "umami: Required field: {field} is missing.".format(
field=field
)
raise ValueError(msg)
# run FlowAccumulator and ChiFinder
self._data_fa, self._data_cf = _create_landlab_components(
self._data_grid,
chi_finder_kwds=chi_finder_kwds,
flow_accumulator_kwds=flow_accumulator_kwds,
)
self._model_fa, self._model_cf = _create_landlab_components(
self._model_grid,
chi_finder_kwds=chi_finder_kwds,
flow_accumulator_kwds=flow_accumulator_kwds,
)
# determine which residuals are desired.
self._residuals = OrderedDict(residuals or {})
self._validate_residuals(self._residuals)
self._distinguish_metric_from_resid()
self._category = None
# set up metric and residual objects
self._data_metric = Metric(
data,
flow_accumulator_kwds=flow_accumulator_kwds,
chi_finder_kwds=chi_finder_kwds,
metrics=self._metrics,
)
self._model_metric = Metric(
model,
flow_accumulator_kwds=flow_accumulator_kwds,
chi_finder_kwds=chi_finder_kwds,
metrics=self._metrics,
)
@property
def names(self):
"""Names of residuals in residual order."""
self._names = []
for key, info in self._residuals.items():
if info["_func"] != "discretized_misfit":
self._names.append(key)
else:
n_f1_levels = np.size(info["field_1_percentile_edges"]) - 1
n_f2_levels = np.size(info["field_2_percentile_edges"]) - 1
label = info["name"]
for f1l in range(n_f1_levels):
for f2l in range(n_f2_levels):
n = label.format(field_1_level=f1l, field_2_level=f2l)
self._names.append(n)
return self._names
@property
def category(self):
"""Category labels from discretized_misfit.
Returns
-------
category: ndarray of size (number of nodes,)
The category labels used with ``discretized_misfit`` or ``None`` if
this calculation is not used.
"""
return self._category
[docs] def value(self, name):
"""Get a specific residual value.
Parameters
----------
name: str
Name of desired residual.
"""
return self._values[name]
@property
def values(self):
"""Residual values in residual order."""
return [self._values[key] for key in self.names]
[docs] def add_from_file(self, file):
"""Add residuals to an ``umami.Residual`` from a file.
Parameters
----------
file_like : file path or StringIO
File will be parsed by ``yaml.safe_load`` and converted to an
``OrderedDict``.
"""
params = _read_input(file)
self.add_from_dict(params)
[docs] def add_from_dict(self, params):
"""Add residuals to an ``umami.Residual`` from a dictionary.
Adding residuals through this method does not overwrite already existing
residuals. New residuals are appended to the existing residual list.
Parameters
----------
params : dict or OrderedDict
Keys are residual names and values are a dictionary describing
the creation of the residual. It will be convereted to an OrderedDict
before residuals are added so as to preserve residual order.
"""
new_residuals = OrderedDict(params)
self._validate_residuals(new_residuals)
for key in new_residuals:
self._residuals[key] = new_residuals[key]
self._distinguish_metric_from_resid()
new_metrics = {}
for name, info in self._metrics.items():
if name not in self._data_metric._metrics:
new_metrics[name] = info
if len(new_metrics) > 0:
self._data_metric.add_from_dict(new_metrics)
self._model_metric.add_from_dict(new_metrics)
[docs] def calculate(self):
"""Calculate residual values.
Calculated residual values are stored in the attribute
``Residual.values``.
"""
self._values = OrderedDict()
self._model_metric.calculate()
self._data_metric.calculate()
for key in self._residuals.keys():
info = deepcopy(self._residuals[key])
_func = info.pop("_func")
if key in self._metrics:
resid = (
self._model_metric._values[key]
- self._data_metric._values[key]
)
else:
function = residual_calcs.__dict__[_func]
resid = function(self._model_grid, self._data_grid, **info)
if _func != "discretized_misfit":
self._values[key] = resid
else:
# if
self._category = resid[0]
for label, value in resid[1].items():
self._values[label] = value
[docs] def write_residuals_to_file(self, path, style, decimals=3):
"""Write residuals to a file.
Parameters
----------
path :
style : str
yaml, dakota
decimals: int
Number of decimals to round output to.
Examples
--------
>>> import numpy as np
>>> from io import StringIO
>>> from landlab import RasterModelGrid
>>> from umami import Residual
>>> np.random.seed(42)
>>> model = RasterModelGrid((10, 10))
>>> z_model = model.add_zeros("node", "topographic__elevation")
>>> z_model += model.x_of_node + model.y_of_node
>>> data = RasterModelGrid((10, 10))
>>> z_data = data.add_zeros("node", "topographic__elevation")
>>> z_model += data.x_of_node + data.y_of_node
>>> z_data[data.core_nodes] += np.random.random(data.core_nodes.shape)
>>> file_like=StringIO('''
... me:
... _func: aggregate
... method: mean
... field: topographic__elevation
... ep10:
... _func: aggregate
... method: percentile
... field: topographic__elevation
... q: 10
... oid1_mean:
... _func: watershed_aggregation
... field: topographic__elevation
... method: mean
... outlet_id: 1
... sn1:
... _func: count_equal
... field: drainage_area
... value: 1
... ''')
>>> residual = Residual(model, data)
>>> residual.add_from_file(file_like)
>>> residual.calculate()
First we ouput in *dakota* style, in which each metric is listed on
its own line with its name as a comment.
>>> out = StringIO()
>>> residual.write_residuals_to_file(out, style="dakota", decimals=3)
>>> file_contents = out.getvalue().splitlines()
>>> for line in file_contents:
... print(line.strip())
17.533 me
9.909 ep10
9.813 oid1_mean
-41 sn1
Next we output in *yaml* style, in which each metric is serialized in
YAML format.
>>> out = StringIO()
>>> residual.write_residuals_to_file(out, style="yaml", decimals=3)
>>> file_contents = out.getvalue().splitlines()
>>> for line in file_contents:
... print(line.strip())
me: 17.533
ep10: 9.909
oid1_mean: 9.813
sn1: -41
"""
if style == "dakota":
stream = "\n".join(
[
str(np.round(val, decimals=decimals)) + " " + str(key)
for key, val in self._values.items()
]
)
if style == "yaml":
stream = "\n".join(
[
str(key) + ": " + str(np.round(val, decimals=decimals))
for key, val in self._values.items()
]
)
_write_output(path, stream)
[docs] @classmethod
def from_dict(cls, params):
"""Create an umami ``Residual`` from a dictionary.
Parameters
----------
params : dict or OrderedDict
This dict must contain a key *grid*, the values of which will be
passed to the `Landlab` function ``create_grid`` to create the
model grid. It will be convereted to an OrderedDict before residuals
are added so as to preserve residual order.
Examples
--------
>>> import numpy as np
>>> from io import StringIO
>>> from umami import Residual
>>> np.random.seed(42)
>>> params = {
... "model": {
... "RasterModelGrid": [
... [10, 10],
... {
... "fields": {
... "node": {
... "topographic__elevation": {
... "plane": [
... {"point": [0, 0, 0]},
... {"normal": [-1, -1, 1]},
... ]
... }
... }
... }
... },
... ]
... },
... "data": {
... "RasterModelGrid": [
... [10, 10],
... {
... "fields": {
... "node": {
... "topographic__elevation": {
... "plane": [
... {"point": [0, 0, 0]},
... {"normal": [-1, -1, 1]},
... ],
... "random" : [
... {"where": "CORE_NODE"},
... {"distribution": "standard_normal"},
... ]
... }
... }
... }
... },
... ]
... },
... "residuals": {
... "me": {
... "_func": "aggregate",
... "method": "mean",
... "field": "topographic__elevation",
... },
... "ep10": {
... "_func": "aggregate",
... "method": "percentile",
... "field": "topographic__elevation",
... "q": 10,
... },
... "oid1_mean": {
... "_func": "watershed_aggregation",
... "field": "topographic__elevation",
... "method": "mean",
... "outlet_id": 1,
... },
... "sn1": {
... "_func": "count_equal",
... "field": "drainage_area",
... "value": 1,
... },
... },
... }
>>> residual = Residual.from_dict(params)
>>> residual.names
['me', 'ep10', 'oid1_mean', 'sn1']
>>> residual.calculate()
>>> np.testing.assert_array_almost_equal(
... residual.value("me"),
... 0.158,
... decimal=3)
>>> np.testing.assert_array_almost_equal(
... residual.values,
... np.array([ 0.158, 0.67 , 4.138, -20. ]),
... decimal=3)
"""
model = create_grid(params.pop("model"))
data = create_grid(params.pop("data"))
return cls(model, data, **params)
[docs] @classmethod
def from_file(cls, file_like):
"""Create an umami ``Residual`` from a file-like object.
Parameters
----------
file_like : file path or StringIO
File will be parsed by ``yaml.safe_load`` and converted to an
``OrderedDict``.
Returns
-------
umami.Residual
Examples
--------
>>> import numpy as np
>>> from io import StringIO
>>> from umami import Residual
>>> np.random.seed(42)
>>> file_like=StringIO('''
... model:
... RasterModelGrid:
... - [10, 10]
... - fields:
... node:
... topographic__elevation:
... plane:
... - point: [0, 0, 0]
... - normal: [-1, -1, 1]
... data:
... RasterModelGrid:
... - [10, 10]
... - fields:
... node:
... topographic__elevation:
... plane:
... - point: [0, 0, 0]
... - normal: [-1, -1, 1]
... random:
... distribution: standard_normal
... where: CORE_NODE
... residuals:
... me:
... _func: aggregate
... method: mean
... field: topographic__elevation
... ep10:
... _func: aggregate
... method: percentile
... field: topographic__elevation
... q: 10
... oid1_mean:
... _func: watershed_aggregation
... field: topographic__elevation
... method: mean
... outlet_id: 1
... sn1:
... _func: count_equal
... field: drainage_area
... value: 1
... ''')
>>> residual = Residual.from_file(file_like)
>>> residual.names
['me', 'ep10', 'oid1_mean', 'sn1']
>>> residual.calculate()
>>> np.testing.assert_array_almost_equal(
... residual.value("me"),
... 0.191,
... decimal=3)
>>> np.testing.assert_array_almost_equal(
... residual.values,
... np.array([ 0.191, 0.426, 3.252, -20. ]),
... decimal=3)
"""
params = _read_input(file_like)
return cls.from_dict(params)
def _distinguish_metric_from_resid(self):
self._metrics = {}
for key, info in self._residuals.items():
func = info["_func"]
if func in metric_calcs.__dict__:
self._metrics[key] = info
def _validate_residuals(self, residuals):
""""""
for key in residuals:
info = residuals[key]
_validate_func(key, info, _VALID_FUNCS)
_validate_fields(self._data_grid, info)
_validate_fields(self._model_grid, info)