Allow re-linking variables at Model creation

[benbovy]

Links (references) to foreign variables are defined at process class creation, but sometimes it might be needed to change these references depending on which processes are included in a model.

As a use case, let's consider the fastscape model implemented in xarray-topo. FlowRouterD8 process needs topographic elevation values to compute flow routing. FlowRouterD8.elevation currently refers to Topography.elevation for which block uplift is applied independently at each time step. However, it might make sense (especially with large time steps) to force computing flow routing based on topography that has been block uplifted at the current time step. A nice, modular solution to this problem would be to add a new process that somehow "steals" the FlowRouterD8.elevation reference and yield a value that is based on block uplifted topography. "stealing" a reference defined in a external process is not yet supported in xarray-simlab.

Taking the example above and with the new design discussed in #19 (and #20) in mind, a possible solution would look like:

@xsimlab.process
def UpliftBeforeFlow(object):
    elevation = xsimlab.link(Topography, 'elevation')
    block_uplift = xsimlab.link(BlockUplift, 'uplift')
    elevation_flow = xsimlab.link(FlowRouterD8, 'elevation',
                                  yield_value=True, relink=True)

    def run_step(self, dt):
        self.elevation_flow = self.elevation + self.block_uplift

Here relink=True specifies that FlowRouterD8.elevation has to be re-linked to UpliftBeforeFlow.elevation_flow (thus discarding its original link to Topography.elevation).

relink=True is allowed only if yield_value=True and if the referenced external variable (i.e., FlowRouterD8.elevation in this example) is itself declared as a link with yield_value=False. An error should be raised when creating the process class if it these conditions are not filled.

Unfortunately, this solution is not straightforward to implement internally. A constraint is that we must leave FlowRouterD8 untouched, which looks like:

def FlowRouterD8(object):
    def __init__(self):
        # ...
        # later will be instance of Topography OR UpliftBeforeFlow
        self._xsimlab_elevation_obj = None

      @property
      def elevation(self):
         return self._xsimlab_elevation_obj.elevation

The problem is that elevation is already declared in UpliftBeforeFlow and is different from elevation_flow:

def UpliftBeforeFlow(object):
    def __init__(self):
        # ...
        # later will be instance of Topography
        self._xsimlab_elevation_obj = None
        # later will be instance of FlowRouterD8
        self._xsimlab_elevation_flow_obj = None

    @property
    def elevation(self):
        return self._xsimlab_elevation_obj.elevation

[benbovy]

If subclassing process classes is supported, then the case explained above can be easily addressed:

@xsimlab.process
class FlowRouterAfterUplift(FlowRouterD8):
    original_elevation = xsimlab.link(Topography, 'elevation')
    uplift = xsimlab.link(BlockUplift, 'uplift')
    elevation = xsimlab.variable(('y', 'x'), yield_value=True)

    def run_step(self, dt):
        self.elevation = self.original_elevation + self.uplift
        super(FlowRouterAfterUplift, self).run_step(dt)

alt_model = fastscape_model.update_processes(
    {'flow_routing': FlowRouterAfterUplift}
)

Additionally, to avoid having to update processes that depend on flow routing, we can chain linked variables, i.e., create a link that itself refers to a link. However, this does not really solves the problem. In the example below, elevation is still linked to the base class used for flow routing.

@xsimlab.process
class StreamPower(object):
    # ...
    elevation = xsimlab.link(FlowRouterD8, 'elevation')

A possible solution would be to accept a string for the 1st argument of xsimlab.link that would correspond to a key of the processes dictionary passed to the Model constructor, e.g.,

@xsimlab.process
class StreamPower(object):
    # ...
    elevation = xsimlab.link('flow_routing', 'elevation')

[benbovy]

Another workaround would be to keep in the example above the link to FlowRouterD8, 'elevation' untouched, and allow at model creation to look for a subclass when a class is not found, i.e., instead of raising an Error because FlowRouterD8 is not found in the given process dict, it will look if there a subclass (here FlowRouterAfterUplift) in the process dict and use it for the link.

Maybe this would impliy that a process class and one of its subclass cannot be both added to a model (which probably doesn't make sense anyway).