Theory IV: Cycles in the graph

Not all model collections are acyclic: models can, directly or indirectly through other models, depend on each other.

We’ll re-use the same Model1, Model2, Model3 as before, and add a cycle:

Show code cell source

Hide code cell source

from math import sin, exp
import inspect
import networkx as nx



def Model1(a,b):
    
    # Complicated *pure* internal calculation
    x = a**2
    y = sin(b)
    z = exp(x * y)
    output = z / (1 + x)

    #--------- Outputs
    return output

    
def Model2(model1_output, c):
    
    # something else complicated and *pure*
    iteration_var = 0.531
    for i in range(1, 4):
        iteration_var  = c * iteration_var + model1_output

    #--------- Outputs
    return iteration_var


def Model3(model2_output, e, c, d):
    
    # something else also complicated and *pure*
    tolerance = 1e-6
    error = float('inf')
    x = model2_output

    while error > tolerance:
        old_x = x

        # Some nonlinear update
        x = (old_x - e) / (1 + abs(old_x - e))
        x = x + c * sin(x) + d
        error = abs(x - old_x)

    #--------- Outputs
    return x


def interface(f):
    return {f.__name__ : {"inputs" : list(inspect.signature(f).parameters.keys()), "outputs" : [f.__name__.lower() + "_output"]}}


def external_inputs(models : set):
    total_inputs  = set()
    model_outputs = set()

    for model in models:
        interface_i = interface(model)[model.__name__]
        total_inputs.update(interface_i["inputs"])
        model_outputs.update(interface_i["outputs"])

    return total_inputs - model_outputs

def Model4(model5_output, a, model1_output):
    return model5_output + a + model1_output

def Model5(model1_output, b, model6_output):
    return model1_output + b * model6_output

def Model6(d, model4_output):
    return d + model4_output

These form a chain: Model4 requires Model5, which requires Model6, which again requires Model4.

There is now no order at all in which the models can be run one-by-one: the topological sort will simply error:

Show code cell source

Hide code cell source

def models_to_graph(models):
    graph = nx.DiGraph()

    # Map each output to the model that produces it

    models = [interface(i) for i in models]
    
    output_to_model = {}

    for model_dict in models:
        for model_name, config in model_dict.items():
            graph.add_node(model_name)

            for output in config.get("outputs", []):
                output_to_model[output] = model_name

    # Create model -> model edges based on inputs
    for model_dict in models:
        for model_name, config in model_dict.items():

            for input_name in config.get("inputs", []):
                producer = output_to_model.get(input_name)

                if producer is not None:
                    # Several variables can join the same pair of models, so the edge
                    # keeps all of them rather than only the last one written.
                    if graph.has_edge(producer, model_name):
                        graph[producer][model_name]["inputs"].append(input_name)
                    else:
                        graph.add_edge(producer, model_name, inputs=[input_name])

    return graph

cycle_graph = models_to_graph([Model1, Model2, Model3, Model4, Model5, Model6])

try:
    list(nx.topological_sort(cycle_graph))
except nx.NetworkXUnfeasible as e:
    print("No valid ordering:", e)

No valid ordering: Graph contains a cycle or graph changed during iteration

Such cycles can be resolved with a fixed point iteration. In order to do this, we do need to know what exactly the cycle is. In this case, it’s straightforward, but it’s not always the case.

A standard graph algorithm that finds such cycles is the strongly connected component algorithm.

Instead of running model-by-model, this forces us to run ‘block’ by block. A block is a model or a set of models that are run together: in the case of a cycle, it’s because they have to be solved together.

If you have an acyclic graph, each model is itself just a block, but not here:

def blocks_in_order(nx_digraph):
    """The models grouped into blocks, and the blocks put in a valid run order."""
    blocks   = list(nx.strongly_connected_components(nx_digraph))
    block_of = {model: i for i, block in enumerate(blocks) for model in block}

    # Order each block:
    block_graph = nx.DiGraph()
    block_graph.add_nodes_from(range(len(blocks)))
    block_graph.add_edges_from((block_of[producer], block_of[consumer])
                               for producer, consumer in nx_digraph.edges
                               if block_of[producer] != block_of[consumer])

    return [sorted(blocks[i]) for i in nx.topological_sort(block_graph)]


blocks_in_order(cycle_graph)

[['Model1'], ['Model2'], ['Model4', 'Model5', 'Model6'], ['Model3']]

We see that we can run Model1 and Model2 first, then have to solve the block cycle, and then run Model3 (which could also be run before the cycle).

Solving a block

Running a single model is unchanged from the previous notebook:

def run_model(model, data_structure):

    # Get I/O
    interface_i = interface(model)[model.__name__]

    # Prepare inputs:
    inputs = {input_name: data_structure[input_name] for input_name in interface_i["inputs"]}

    # Return updated datastructure:
    return {**data_structure, interface_i["outputs"][0]: model(**inputs)}

A block of several models cannot be run like that: whichever model we start with is missing one of its inputs. The way out is to assume a value for one of the variables inside the cycle [1] and make that variable the iteration variable.

Given that assumption, in a block, that model output will no longer be needed and thus it can now be ordered:

def block_order(block, nx_digraph, assumed):
    """The order to run a block's models in, once `assumed` is taken as given."""
    block_graph = nx_digraph.subgraph(block).copy()

    # Taking `assumed` as given means the edges carrying it are no longer dependencies.
    block_graph.remove_edges_from([(producer, consumer)
                                   for producer, consumer, edge in block_graph.edges(data=True)
                                   if assumed in edge["inputs"]])

    return list(nx.topological_sort(block_graph))


cycle = ['Model4', 'Model5', 'Model6']

for assumed in ('model4_output', 'model5_output', 'model6_output'):
    print(f"assuming {assumed:15s} ->", block_order(cycle, cycle_graph, assumed))

assuming model4_output   -> ['Model6', 'Model5', 'Model4']
assuming model5_output   -> ['Model4', 'Model6', 'Model5']
assuming model6_output   -> ['Model5', 'Model4', 'Model6']

Each choice gives a different (valid) order. Now the block can be run: assume a value, run the models in that order, and see what comes back out. Convergence is achieved when the resulting value is the same as the assumed value:

def solve_block(block, data_structure, models, nx_digraph, assumed,
                start = 1.0, tolerance = 1e-10, max_iterations = 100):
    by_name = {model.__name__: model for model in models}
    order   = block_order(block, nx_digraph, assumed)

    # We have to start somewhere: nothing has produced `assumed` yet.
    ds = {**data_structure, assumed: start}

    for iteration in range(max_iterations):
        old_value = ds[assumed]

        for model_name in order:
            ds = run_model(by_name[model_name], ds)

        # Converged when running the block reproduces the value we assumed.
        if abs(ds[assumed] - old_value) < tolerance:
            print(f"{block}: converged in {iteration + 1} iterations, assuming '{assumed}'")
            return ds

    print(f"{block}: did NOT converge in {max_iterations} iterations, assuming '{assumed}'")
    return ds

And that is enough to run any collection of models, cyclic or not:

def TotalModel(data_structure, models : set, assumed):
    ds = data_structure.copy()

    if not set(external_inputs(models)) == set(ds.keys()):
        raise ValueError(f"The provided data structure keys {set(ds.keys())} " +
                         f"do not match the required external inputs {set(external_inputs(models))} for the given models.")

    nx_digraph = models_to_graph(models)
    by_name    = {model.__name__: model for model in models}

    for block in blocks_in_order(nx_digraph):

        if len(block) == 1:
            # An ordinary model: everything it reads is available already.
            ds = run_model(by_name[block[0]], ds)
        else:
            # A cycle: it has to be solved. (note we start with a value of 1)
            ds = solve_block(block, ds, models, nx_digraph, assumed)

    return ds


all_models = {Model1, Model2, Model3, Model4, Model5, Model6}
ds_start   = {"a": 1.0, "b": 0.5, "c": 2.0, "d": 0.1, "e": 0.2}

TotalModel(ds_start, all_models, assumed = "model5_output")

['Model4', 'Model5', 'Model6']: converged in 35 iterations, assuming 'model5_output'

{'a': 1.0,
 'b': 0.5,
 'c': 2.0,
 'd': 0.1,
 'e': 0.2,
 'model1_output': 0.8075731482210419,
 'model2_output': 9.901012037547295,
 'model5_output': 3.5227194445897045,
 'model4_output': 5.330292592737326,
 'model6_output': 5.430292592737326,
 'model3_output': 1.9113901980206618}

This works, and the nice thing is it works for any model collection like this! You just have to assume a model output and then the fixed point iteration will handle it.

However, there were two big hidden assumptions:

1. We chose which variable is the iteration variable. All three would have worked, so why this particular one?

TotalModel(ds_start, all_models, assumed = "model4_output")

['Model4', 'Model5', 'Model6']: converged in 36 iterations, assuming 'model4_output'

{'a': 1.0,
 'b': 0.5,
 'c': 2.0,
 'd': 0.1,
 'e': 0.2,
 'model1_output': 0.8075731482210419,
 'model2_output': 9.901012037547295,
 'model4_output': 5.330292592821154,
 'model6_output': 5.430292592758139,
 'model5_output': 3.5227194446001113,
 'model3_output': 1.9113901980206618}

TotalModel(ds_start, all_models, assumed = "model6_output")

['Model4', 'Model5', 'Model6']: converged in 36 iterations, assuming 'model6_output'

{'a': 1.0,
 'b': 0.5,
 'c': 2.0,
 'd': 0.1,
 'e': 0.2,
 'model1_output': 0.8075731482210419,
 'model2_output': 9.901012037547295,
 'model6_output': 5.430292592819699,
 'model5_output': 3.5227194445986565,
 'model4_output': 5.330292592819699,
 'model3_output': 1.9113901980206618}

2. The algorithm: why should we use a fixed point iteration always? In some cases this might not converge, and we don’t know whether the failure is algorithmic or whether the problem really has no answer: [2]

TotalModel({**ds_start, "b": 1.1}, all_models, assumed = "model5_output")

['Model4', 'Model5', 'Model6']: did NOT converge in 100 iterations, assuming 'model5_output'

{'a': 1.0,
 'b': 1.1,
 'c': 2.0,
 'd': 0.1,
 'e': 0.2,
 'model1_output': 1.2190357525781628,
 'model2_output': 12.78125026804714,
 'model5_output': 533268.563728798,
 'model4_output': 484788.39517549577,
 'model6_output': 484788.49517549574,
 'model3_output': 1.911390209140226}

These choices cannot be made just on the basis of the graph (i.e. algorithmically); there’s always some ambiguity.

And there are more choices still: whether two cycles that happen to sit next to each other should be one block or two, or whether the cycle should not be solved here at all, but handed to an optimizer as a consistency constraint.

Summary

A cyclic graph has no valid model ordering, but it does have a valid block ordering: the strongly connected components, ordered by sorting a graph of the blocks themselves. A block of one model is run; a block of several is a cycle and has to be solved. The graph-based representation enables us to look at the minimal cycles and automatically assign a solver (in this notebook a fixed point iteration) to it.

Solving one requires assuming a variable inside it, which cuts the cycle and makes the block orderable again. Unlike the ordering itself, that choice (and the choice of solver, and of what is solved together) cannot be read off the graph.

In other words, we need to, as users, take a look at the graph and make an informed decision.

In the next notebook, we’ll look at how to represent the graph in an easier-to-understand way, and look at other possible cycles.

---

[1]

In later notebooks we will formalize this notion in the sense of MDAO architectures; for now this demonstrates the mechanism as well.

[2]

Here it does have one. The three models in the cycle are affine, so the block has exactly one solution (\(\text{model5\_output} = -37.70\)). The sweep gain is \(b\) whichever variable is assumed, so at \(b = 1.1\) Gauss–Seidel diverges from all three choices; Newton, which is exact in one step on an affine residual, finds it immediately.