Theory I: Imperative programming problems

Imperative programming is one of the most dominant scientific computing paradigms. The program explicitly describes a sequence of operations:

from math import sin, exp
def total_program(inputs):

    start_var = inputs[0]

    for i in range(1, 4):
        start_var = inputs[i] * start_var
    print("Total model output:", start_var)
    return start_var

total_program([1.0, 2.0, 3.0, 4.0])

Total model output: 24.0

24.0

This style of programming is very flexible and places little restrictions on what you are allowed to do inside a function or class.

Let’s make a simple collection of models that take a data structure and update it with new values[1]:

def Model1(data_structure):
    a = data_structure['a']
    b = data_structure['b']

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

    data_structure['model1_output'] = output

def Model2(data_structure):
    model_1_output = data_structure['model1_output']
    c              = data_structure['c']

    data_structure['model2_output'] = 0.531

    # something else complicated
    for i in range(1, 4):
        data_structure['model2_output']  = c * data_structure['model2_output'] + model_1_output



def Model3(data_structure):
    model_2_output = data_structure['model2_output']
    e              = data_structure['e']
    c              = data_structure['c']
    d              = data_structure['d']

    # something else also complicated
    tolerance = 1e-6
    error = float('inf')
    x = model_2_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)

    data_structure['model3_output'] = x

It’s very easy to see how these can be combined to create a total model:

def TotalModel(data_structure):
    Model1(data_structure)
    Model2(data_structure)
    Model3(data_structure)

def base_data_structure():
    # We create a new datastructure with just the inputs.
    return {'a': 1.0, 'b' : 2.0, "c" : 3.0, "d" : 4.0, "e" : 5.0}


TotalModel(cds_1 := base_data_structure())

cds_1

{'a': 1.0,
 'b': 2.0,
 'c': 3.0,
 'd': 4.0,
 'e': 5.0,
 'model1_output': 1.2412888640075004,
 'model2_output': 30.473755232097503,
 'model3_output': 1.0547835473689968}

For such small model collections this is fine; but for larger and more complex models, this type of programming becomes intractable to do efficiently.

Ordering and flexibility

Let’s say we want to add a new model, which now produces a value that was previously an input (e)

def Model4(data_structure):
    model_1_output = data_structure['model1_output']
    f              = data_structure['f']
    data_structure['e'] = model_1_output + f

This requires rewriting a new TotalModel. But where do we put it? We have to know what each model internally does, up to the smallest details; since there’s no guarantee that model 3 doesn’t actually do something with f or model2_output.

In this case, we know what each model does, so we know we can put Model4 after Model1 and before Model3.

def TotalModelWithModel4(data_structure):
    Model1(data_structure)
    Model4(data_structure)
    Model2(data_structure)    
    Model3(data_structure)

try:
    TotalModelWithModel4(cds_2 := base_data_structure())
except KeyError as e:
    print("Key Not Found:", e)

Key Not Found: 'f'

This doesn’t work yet: we forgot to make sure the data structure also knows about what each model internally accesses. There’s no ‘f’ in it.

We’ll have to create a new data structure that has ‘f’ in it.

def new_data_structure():
    return {    'a': 1.0, 'b' : 2.0, "c" : 3.0, "d" : 4.0, "e" : 5.0, "f" : 5.0}


TotalModelWithModel4(cds_3 := new_data_structure())

cds_3

{'a': 1.0,
 'b': 2.0,
 'c': 3.0,
 'd': 4.0,
 'e': 6.2412888640075,
 'f': 5.0,
 'model1_output': 1.2412888640075004,
 'model2_output': 30.473755232097503,
 'model3_output': 0.9207100747728094}

For now this works: but imagine you have a large team, with everyone developing models independently. You’d also have to know what those others models do internally to figure out where your model fits in the order (which involves reading the entire source code of all models).

Furthermore, let’s imagine we want to make some models behave differently depending on some choice of model.

def AlternativeModel4(data_structure, behaviour : int):
    # Behaviour takes different input:
    if behaviour == 0: 
        aa = data_structure['model2_output']
    elif behaviour == 1:
        aa = data_structure['model1_output']        

    c = data_structure['c']
    data_structure['e'] = aa + c


def AlternativeModel2(data_structure, behaviour : int):

    # Behaviour takes different input:
    if behaviour == 0:
        dd = data_structure['f'] 
    elif behaviour == 1:
        dd = data_structure['e']  # model 4 produces this!       

    model_1_output = data_structure['model1_output']
    data_structure['model2_output'] = model_1_output * dd **0.5
    

How do we order this? Let’s try:

def TotalModel(data_structure, behaviour : int):

    Model1(data_structure)
    AlternativeModel2(data_structure , behaviour)
    AlternativeModel4(data_structure, behaviour)
    Model3(data_structure)

TotalModel(cds_4 := new_data_structure(), behaviour = 0)
cds_4

{'a': 1.0,
 'b': 2.0,
 'c': 3.0,
 'd': 4.0,
 'e': 5.775606279634263,
 'f': 5.0,
 'model1_output': 1.2412888640075004,
 'model2_output': 2.775606279634263,
 'model3_output': 0.9623193120848894}

Okay great, works: but now with a different behaviour:

try:
    TotalModel(cds_5 := new_data_structure(), behaviour = 1)
    print("Works")
except:
    print("Doesn't work!")

cds_5

Works

{'a': 1.0,
 'b': 2.0,
 'c': 3.0,
 'd': 4.0,
 'e': 4.2412888640075,
 'f': 5.0,
 'model1_output': 1.2412888640075004,
 'model2_output': 2.775606279634263,
 'model3_output': 1.1961853438210768}

Great this works too!

But wait: if behaviour is 1, this actually shouldn’t! Then, model 2 depends on e, which is produced by model 4, which is run after model2, so how does it even work?

Well, it uses the value of e that was still present in the input data structure. There’s no check or guarantee there. If we remove e from the input data structure, it will break:

def new_data_structure_without_e():
    return {    'a': 1.0, 'b' : 2.0, "c" : 3.0, "d" : 4.0, "f" : 5.0}

try: 
    TotalModel(cds_6 := new_data_structure_without_e(), behaviour = 1)
    print("Works")
except KeyError as e:
    print("KeyError", e)

KeyError 'e'

The only way out is to make the ordering dependent on behaviour:

def TotalModelSwitch(data_structure, behaviour : int):

    Model1(data_structure)

    # Behaviour-dependent ordering:
    if behaviour == 0:
        AlternativeModel2(data_structure , behaviour)
        AlternativeModel4(data_structure, behaviour)
    elif behaviour == 1:
        AlternativeModel4(data_structure, behaviour)
        AlternativeModel2(data_structure , behaviour)

    Model3(data_structure)


TotalModelSwitch(cds_7 := new_data_structure_without_e(), behaviour = 0)
cds_7

{'a': 1.0,
 'b': 2.0,
 'c': 3.0,
 'd': 4.0,
 'f': 5.0,
 'model1_output': 1.2412888640075004,
 'model2_output': 2.775606279634263,
 'e': 5.775606279634263,
 'model3_output': 0.9623193120848894}

TotalModelSwitch(cds_8 := new_data_structure_without_e(), behaviour = 1)
cds_8

{'a': 1.0,
 'b': 2.0,
 'c': 3.0,
 'd': 4.0,
 'f': 5.0,
 'model1_output': 1.2412888640075004,
 'e': 4.2412888640075,
 'model2_output': 2.5563586574101667,
 'model3_output': 1.1961854882696068}

Notice now that model2_output is different than the previous ‘working’ one: 2.55 here vs 2.78 there.

For such small models, this solves the problem. However, imagine you’d have 5 models each have 3 possible alternative implementations. That means you’d have to figure out for \(3^5 = 243\) model combinations what the model computation order should be, if it even exists!. And you’d have to track what exactly all models do in their internals before you can even try to figure it out.

It’s also quite error-prone: initialising everything with 0 will work for arbitrary orders, it just won’t be correct (like we saw before!), and you’d have no way to validate it easily.

Cycles

If we have models whose inputs depend on eachothers outputs, we can’t actually write a TotalModel in a simple sequential manner:

def Model5(data_structure):
    e = data_structure['e']
    data_structure['f'] = e**0.7



def TotalModelCycle(data_structure, four_before_five : bool):

    Model1(data_structure)
    Model2(data_structure)

    # Try both!
    if four_before_five: 
        Model4(data_structure)
        Model5(data_structure)
    else:
        Model5(data_structure)
        Model4(data_structure)

    Model3(data_structure)

# Also have to seed the dict without 'f' now: model 5 produces it.
def new_data_structure_without_ef():
    return {    'a': 1.0, 'b' : 2.0, "c" : 3.0, "d" : 4.0}

try:
    TotalModelCycle(cds_9 := new_data_structure_without_ef(), True)
except Exception as e:
    print("KeyError: ", e)


try:
    TotalModelCycle(cds_10 := new_data_structure_without_ef(), False)
except Exception as e:
    print("KeyError: ", e)

KeyError:  'f'
KeyError:  'e'

Both orderings fail because Model4 and Model5 form a cyclic dependency: each requires the output of the other. Note that if e and f had been present in the dictionary — e.g. from a previous computation — both orders would have worked, but they would have been incorrect.

The natural reaction is to add a fixed point iteration here:

def TotalModelWithFixedPoint(data_structure, n_iter : int):


    Model1(data_structure)
    Model2(data_structure)

    data_structure['f']  = 1.0

    tolerance = 1e-13

    for i in range(n_iter):
        old_f = data_structure['f']
        Model4(data_structure)
        Model5(data_structure)
        if abs(data_structure['f'] - old_f) < tolerance:
            print(f"Fixed point iteration converged after {i+1} iterations.")
            break
    else:
        print(f"Fixed point iteration did not converge in {n_iter} iterations.")
    
    Model3(data_structure)


TotalModelWithFixedPoint(cds_11 := new_data_structure_without_ef(), 100)
cds_11

Fixed point iteration converged after 41 iterations.

{'a': 1.0,
 'b': 2.0,
 'c': 3.0,
 'd': 4.0,
 'model1_output': 1.2412888640075004,
 'model2_output': 30.473755232097503,
 'f': 2.5346441657675944,
 'e': 3.7759330297750244,
 'model3_output': 6.799615636706353}

This works fine. But now if you have an alternative Model4 that doesn’t actually use ‘f’ at all; you’re wasting the fixed point iteration since there’s no dependency anymore. And if you have a single alternative this might be doable, but what if you have 243 possible combinations? You’d need to find all possible fixed point iterations and at some point you’d probably give up and put the entire thing in a fixed point iteration, which is very wasteful if you only have a small cycle.

Summary

Imperative programming is widely used in scientific computing, is flexible and intuitive, and matches closely with how a computer actually executes code.

But for coupling a large number of models there are a number of problems:

• A priori, for a given set of models, we don’t actually have a clear model-computation ordering. Hand-deriving them is error-prone and inflexible with respect to alternative models.

• When models have cycles, imperative programming requires you to reason about these possible cycles and implement the correct solvers or fixed point iterations, which again is error-prone and inflexible with respect to alternative models.

---

[1]

This is not meant as a standard or even particularly realistic way to implement model couplings: nevertheless, imperative programming allows you to do this, and it shows easily why this can be a problem. Making it better (explicit inputs /outputs & no mutable state for no hidden dependencies) is possible, but moves the design toward a functional programming style. This is precisely the direction which we will argue in the next notebooks: if you fully commit to it you’ll not only make it easier to reason about, it will be automatically be validatable!