Tutorial
The basic goal of this tutorial is to get up and running with the amisc package. It assumes you have a working Python 3.11+ environment with amisc >= v0.5.0 properly installed.
Contents:
- Introduce the basic classes
- Interfacing an external model
- Connecting models together
- Working with model fidelities
- Training a surrogate model
- Evaluate surrogate performance
- Saving and loading from file
Estimated time to complete: 1 hour Download Notebook
The basic classes¶
Multidisciplinary (MD) systems are constructed from 3 objects:
- Variables - These are inputs, outputs, quantities of interest (QoIs), etc. They are the most basic element of any MD system and serve as the datapaths or connections between models. They can be random variables, scalars, field quantities, etc.
- Components - These are the individual discipline models. They map a set of inputs to a set of outputs.
- System - This is the top-level MD system. It connects multiple component models and manages the flow of data between them.
So let's start by defining some variables.
from amisc import Variable
x = Variable()
y = Variable(description='My variable')
z = Variable(description='Another variable', units='m/s')
print(x, y, z) # will show "x, y, z"
x y z
As you can see, variables are nothing more than a placeholder, just like in any equation. You can treat them mostly as a string with some extra optional properties (like units or a description). In fact, comparing variables to their string counterparts will return True!
print(x == 'x')
print(x == x.name) # alternatively
True True
Of course, these variables can do more than just plain strings -- we will see some more examples later.
For now, let's move to constructing a component model. A model is simply a function that takes a set of inputs and computes a set of outputs. For example, say you have an algebraic "model" (i.e. function) that predicts the force of gravity between two objects:
$$ F_g = G\frac{m_1 m_2}{r^2} $$
We can write this as a Python function fairly simply:
def force_of_gravity(m1, m2, r):
G = 6.6743e-11
Fg = G * (m1 * m2) / r ** 2
return Fg
And now, we can integrate this model into amisc by wrapping it with the Component:
from amisc import Component
gravity_model = Component(force_of_gravity)
We can view the input and output variables by printing our new component:
print(gravity_model)
---- gravity_model ---- Inputs: [m1, m2, r] Outputs: [Fg] Model: <function force_of_gravity at 0x7f55a98cdee0>
Note how the m1, m2, r inputs were inspected directly from the function signature, as well as the Fg output from the return statement -- this highlights that Component models are nothing more than a callable function with a set of Variable inputs and Variable outputs.
Let's now say we have a second model that predicts the acceleration of our second object under only the influence of gravity:
$$ A_g = F_g / m_2 $$
And the corresponding Component:
def acceleration(Fg, m2):
Ag = Fg / m2
return Ag
accel_model = Component(acceleration)
We would now like to predict both models in concert, i.e. first compute $F_g$ using the gravity model and then $A_g$ using the acceleration model. We construct this MD system as:
from amisc import System
md_system = System(gravity_model, accel_model)
We can now make all of our predictions with a single call on the System:
pred = md_system.predict({'m1': 5.972e24, # kg (the Earth)
'm2': 68, # kg (and you)
'r' : 6.37e6 # m (on the surface of Earth)
}, use_model='best')
Now, printing the outputs, we see the value of the force as Fg in Newtons and the expected acceleration due to gravity of $A_g = 9.8\ m/s^2$:
print(pred)
{'Fg': array([667.96786664]), 'Ag': array([9.82305686])}
Interfacing a model¶
The last example was fairly simplified, since we only had three inputs $m_1, m_2$ and $r$, and the model had a known, algebraic expression. In general, we might have many more inputs/outputs, each with more details like units, domains, nominal values, etc. and our model might be a more complex finite-element simulation, for example, that maybe runs on a remote machine or uses some software outside of Python.
From the Component's perspective, all of these details are irrelevant -- the interface that we used for the simple gravity model applies to any model. That is to say, the Component is a black-box wrapper that interfaces any external model into the amisc framework. The beauty of Python as a "glue" language means we can make any external calls we need right from within our Component model.
As an example, let's say we have some binary on our local machine that reads from an input.ini config file, runs a finite-element thermal analysis of a heater component, and writes to a results.out file. If we want to study the impact of the ambient temperature $T_a$ and the heating element diameter $d$ on the maximum thermal stress $\sigma$ in the component, then we would interface the model with amisc like this:
from amisc import Component
def thermal_analysis(Ta, d):
"""Run FEA thermal study of a heater component.
:param Ta: the ambient temperature (K)
:param d: the heating component diameter (m)
:returns: sigma, the maximum thermal stress (kPA)
"""
import subprocess
with open('input.ini', 'a') as fd:
fd.writelines([f'Ambient temperature = {Ta}',
f'Component diameter = {d}'
])
subprocess.run(['my_fea_binary', '--input_file', 'input.ini', '--out', 'results.out'])
with open('results.out', 'r') as fd:
results = fd.readlines()
sigma = float(results[-1]) # hypothetical, if the result was stored on the last line
return sigma
thermal_model = Component(thermal_analysis)
print(thermal_model)
---- thermal_model ---- Inputs: [Ta, d] Outputs: [sigma] Model: <function thermal_analysis at 0x7f55a97e7740>
As the variables grow in number and complexity, it will often be more convenient to batch all inputs or outputs together and unpack them inside the model where needed. For example, if the thermal model above had many more inputs or outputs, we would define the Component to take a single set of inputs and return a single set of outputs:
def thermal_analysis(inputs): # single `dict` of inputs
Ta = inputs['Ta']
d = inputs['d']
other_inputs = ...
# Compute the model
sigma = float(results[-1])
other = ...
return {'sigma': sigma, 'other': ...} # single `dict` of outputs
inputs = ['Ta', 'd', 'x1', 'x2'] # ...
outputs = ['sigma', 'other'] # ...
thermal_model = Component(thermal_analysis, inputs, outputs)
print(thermal_model)
---- thermal_model ---- Inputs: [Ta, d, x1, x2] Outputs: [sigma, other] Model: <function thermal_analysis at 0x7f55a9aafe20>
As one last note about models, amisc will always treat positional arguments as numeric input Variable's. If you need to pass extra arguments or configs to your model, but don't necessarily want them to be considered as numeric variables, then you can pass them as keyword arguments to Component. For example, say your model needs to load additional settings from a config file, then you would pass config_file as a kwarg to Component:
def thermal_analysis(numeric_inputs, config_file=None, **kwargs):
Ta = numeric_inputs['Ta']
d = ...
# etc.
with open(config_file, 'r') as fd:
# Read in additional configs for the model
configs = fd.readlines()
return numeric_outputs
inputs = ['Ta', 'd']
outputs = ['sigma']
kwargs = {'extra': 'configs', 'go': 'here'}
thermal_model = Component(thermal_analysis, inputs, outputs, config_file='my_config.json', **kwargs)
Connecting models together¶
We have already seen that connecting models together is as easy as wrapping them in a System(comp1, comp2, ...). Since each individual Component specifies its inputs and outputs, the System automatically connects the components into a graph-like data structure by drawing edges between one component's outputs and another component's inputs.
Let's look at an example of a three-component, purely feed-forward system:
\begin{align} y_1 &= f_1(x) = x \sin(\pi x)\\ y_2 &= f_2(y_1) = \frac{1}{1 + 25y_1^2}\\ y_3 &= f_3(x, y_2) = x \cos(\pi y_2) \end{align}
To implement this MD system in amisc, we define each component model and wrap them in a System. Since this is a simple example and there are no extra kwargs, we'll skip the definition of each Component and pass the models directly to System:
import numpy as np
def f1(x):
y1 = x * np.sin(np.pi * x)
return y1
def f2(y1):
y2 = 1 / (1 + 25*y1**2)
return y2
def f3(x, y2):
y3 = x * np.cos(np.pi * y2)
return y3
md_system = System(f1, f2, f3)
We can inspect the System graph to show that we do indeed get the $f_1\rightarrow f_2\rightarrow f_3$ coupling that we expected:
import networkx as nx
nx.draw(md_system.graph(), with_labels=True)
print(md_system)
---- System_157 ---- amisc version: 0.5.2 Refinement level: 0 Components: f1, f2, f3 Inputs: x Outputs: y3, y2, y1
We can now view the model outputs over the range $x\in(0, 1)$.
import matplotlib.pyplot as plt
x_grid = {'x': np.linspace(0, 1, 100)}
y_grid = md_system.predict(x_grid, use_model='best')
fig, ax = plt.subplots(1, 3, figsize=(8, 3), layout='tight', sharey='row')
ax[0].plot(x_grid['x'], y_grid['y1'], '-r')
ax[1].plot(x_grid['x'], y_grid['y2'], '-r')
ax[2].plot(x_grid['x'], y_grid['y3'], '-r')
ax[0].set_xlabel('$x$')
ax[1].set_xlabel('$x$')
ax[2].set_xlabel('$x$')
ax[0].set_ylabel('$f(x)$')
plt.show()
One of the primary motivations for studying models in a "decoupled" way like this is that individual component models are typically much simpler when treated independently than when they are convolved with other models. As we will see, approximating the models in this decoupled fashion has many advantages over the traditional "black-box" approach.
Before we get to building surrogates, we must first understand one of the core features of amisc, and that is multi-fidelity hierarchies.
Model fidelities¶
The AMISC algorithm builds on the principle that we can approximate a high-fidelity model by summing or "averaging" over many lower-fidelity versions of the model. As such, the primary way to build surrogates in amisc is by specifying some set of fidelity levels.
The fidelity "level" of a model can very simply be represented as a whole number, with 0 being the "worst" fidelity and counting up 1, 2, 3, ... to higher levels. For example, if you are building a model of a fluid, level 0 might be an analytical potential flow model, level 1 might be incompressible flow, level 2 might be a RANS solver, and so on up the scale, with each level bringing greater accuracy.
More generally, you might have more than one knob by which to tune the fidelity of your model. So instead of one number associated with one fidelity "level", you would have a set of numbers associated with each level, i.e. $(0, 0, ...)$, or $(0, 1, ...)$ etc. We call each set of these numbers a multi-index. AMISC itself is an algorithm for adaptively building up groups of these multi-indices to approximate a model.
For our purposes, overall fidelity is divided into two multi-indices: physical model fidelity and parametric model fidelity. Parametric fidelity is further grouped into training data and surrogate fidelities. Each are summarized below:
- Physical fidelity - denoted by the multi-index $\alpha = (\alpha_1, \alpha_2, \dots)$ → controls the approximation accuracy of the actual physics of the underlying model. This is what one typically thinks of when they hear model "fidelity". This includes PDE mesh refinement, time step size, simplifying physics, etc.
- Parametric fidelity - denoted by the multi-index $\beta = (\beta_1, \beta_2, \dots)$ → controls the approximation accuracy of the surrogate or metamodel itself. We further divide $\beta$ into indices that control the amount of training data (data fidelity) and indices that control the complexity of the surrogate (surrogate fidelity). As $\beta$ increases, both the amount of training data in the approximation increases, as well as the surrogate complexity (i.e. number of hidden layers, nodes, etc.), in both cases resulting in a more accurate surrogate model.
As an example, consider the model:
\begin{align} y &= f(x) = \cos\left(\frac{\pi}{2}(x + \frac{4}{5} + \frac{\epsilon}{5})\right)\\ \epsilon &= 2^{-\alpha_0}\\ \text{for}\ \alpha_0 &= (0, 1, 2, \dots)\\ \end{align}
As the $\alpha$ index increases, $\epsilon\rightarrow 0$ and the "fidelity" of the model increases. Let's build this model and show some predictions for varying $\alpha$ fidelity. Note that we should request the model_fidelity keyword in our component model.
def multilevel_model(inputs, model_fidelity=(0,)):
alpha0 = model_fidelity[0]
eps = 2**(-alpha0)
return {'y': np.cos(np.pi/2 * (inputs['x'] + 4/5 + (1/5)*eps))}
mf_comp = Component(multilevel_model, Variable('x'), Variable('y'), model_fidelity=(2,))
# Plot for varying fidelity
x_grid = {'x': np.linspace(-1, 1, 100)}
y_truth = multilevel_model(x_grid, (10,)) # Use a really high alpha for the "truth"
fig, ax = plt.subplots()
ax.plot(x_grid['x'], y_truth['y'], '-k', label='Truth')
for i in range(3):
y_grid = mf_comp.call_model(x_grid, model_fidelity=(i,))
ax.plot(x_grid['x'], y_grid['y'], label=f'$\\alpha$={i}')
ax.set_xlabel('Input ($x$)')
ax.set_ylabel('Output ($y$)')
ax.legend()
plt.show()
Note how we chose to use the internal Component.call_model function, which is a lightweight wrapper around the actual model itself -- we will see later some of the benefits of using call_model over the raw model function.
Training a surrogate¶
Now that we have an understanding of the component models and their fidelity multi-indices, let's get on with the main purpose of amisc -- building surrogate approximations of the models.
As an example, let's say you have a model that can be refined up to $\alpha=(2,)$ and $\beta=(2,)$. We would start by building the $(\alpha, \beta)=(0, 0)$ multi-index, then the $(0, 1)$ index, then $(1, 0)$, and so on.
Here is a graph that illustrates this surrogate building process:
import matplotlib.pyplot as plt
import matplotlib.patches as patches
fig, ax = plt.subplots()
ax.set_xlabel(r'Parametric fidelity ($\beta$)')
ax.set_ylabel(r'Model fidelity ($\alpha$)')
ax.set_xlim(0, 3)
ax.set_ylim(0, 3)
ax.set_xticks([0.5, 1.5, 2.5])
ax.set_xticklabels(['0', '1', '2'])
ax.set_yticks([0.5, 1.5, 2.5])
ax.set_yticklabels(['0', '1', '2'])
# Draw a box for each multi-index
boxes = [
(0, 0, 'red'),
(1, 0, 'red'),
(0, 1, 'red'),
(1, 1, 'gray'),
(0, 2, 'gray'),
(2, 0, 'gray')
]
width = 1
height = 1
alpha = 0.5
# Add boxes to the plot
for i, (x, y, color) in enumerate(boxes):
label = 'Activated' if color == 'red' else 'Candidate'
rect = patches.Rectangle((x, y), width, height, linewidth=1,
label=label if i in [0, 3] else '',
edgecolor=color, facecolor=color, alpha=alpha)
ax.add_patch(rect)
ax.legend()
plt.show()
Note how the red boxes are the "activated" multi-indices (i.e. the ones we have built so far) and the gray boxes are the next "candidate" multi-indices, or the nearest neighbors that we may choose to build next.
In higher dimensions, this diagram would get harder to visualize, but the premise is the same: we approximate a model by stacking up lower-fidelity approximations in an iterative fashion.
Recall the multi-level model from earlier:
\begin{align} y &= f(x) = \cos\left(\frac{\pi}{2}(x + \frac{4}{5} + \frac{\epsilon}{5})\right)\\ \epsilon &= 2^{-\alpha_0}\\ \text{for}\ \alpha_0 &= (0, 1, 2, \dots)\\ \end{align}
Let's iteratively construct a surrogate for this model, with $\beta$ representing the amount of training data (note the data_fidelity keyword). We'll also train the surrogate over the domain $x\in (-1, 1)$.
comp = Component(multilevel_model,
Variable('x', domain=(-1, 1)),
Variable('y'),
model_fidelity=(2,),
data_fidelity=(2,))
import itertools
print('alpha beta')
for multi_index in itertools.product(range(3), range(3)):
alpha = multi_index[:1]
beta = multi_index[1:]
print(f'{str(alpha):5s} {str(beta):4s}')
comp.activate_index(alpha, beta) # 'training' happens here as gray boxes -> red boxes
alpha beta (0,) (0,) (0,) (1,) (0,) (2,) (1,) (0,) (1,) (1,) (1,) (2,) (2,) (0,) (2,) (1,) (2,) (2,)
Now that we've "activated" all of the multi-indices, let's take a look at the results. In the following grid, we'll plot increasing $\beta$ fidelity along the x-axis and increasing $\alpha$ fidelity along the y-axis.
xg = {'x': np.linspace(-1, 1, 100)}
y_truth = comp.call_model(xg, model_fidelity=(15,))
fig, axs = plt.subplots(3, 3, sharey='row', sharex='col')
for alpha in range(3):
for beta in range(3):
ax = axs[2-alpha, beta]
xi, yi = comp.training_data.get((alpha,), (beta,))
y_surr = comp.interpolator.predict(xg, comp.misc_states[(alpha,), (beta,)], (xi, yi))
s = rf'$\hat{{f}}_{{{alpha}, {beta}}}$'
ax.plot(xg['x'], y_surr['y'], '--k', label=r'{}'.format(s), linewidth=1.5)
s = rf'$\hat{{f}}_{alpha}$'
ax.plot(xg['x'], comp.call_model(xg, (alpha,))['y'], '--b', label=r'{}'.format(s), linewidth=2)
ax.plot(xg['x'], y_truth['y'], '-r', label=r'$f$', linewidth=2)
ax.plot(xi['x'], yi['y'], 'ok')
ax.set_xlabel('$x$' if alpha == 0 else '')
ax.set_ylabel('$f(x)$' if beta == 0 else '')
ax.legend()
fig.text(0.5, 0.02, r'Increasing parametric fidelity ($\beta$) $\rightarrow$', ha='center', fontweight='bold')
fig.text(0.02, 0.5, r'Increasing model fidelity ($\alpha$) $\rightarrow$', va='center', fontweight='bold',
rotation='vertical')
fig.set_size_inches(9, 9)
fig.tight_layout(pad=3, w_pad=1, h_pad=1)
plt.show()
In the grid above, the true function that we wish to approximate is the red $f(x)$. In each row, the dashed blue $\hat{f}_{\alpha}$ gives the various physical model fidelity approximations. The black markers give the surrogate training data. As we move from left to right, $\beta$ increases and so more training data is provided for a given $\hat{f}_{\alpha}$ and the surrogate approximation (black dashed line) improves. Finally, as we move from bottom to top, the surrogate matches the highest physical model fidelity better and better.
While we have shown that you can train the surrogate by iteratively calling activate_index() on the Component, the AMISC algorithm provides a globally adaptive procedure for doing this automatically. While we won't go into the details for the purposes of this tutorial, the general idea is that we search over all the possible combinations of $(\alpha, \beta)$ for every component and select the indices with the most potential for improvement.
To do this, we'll wrap our component in the System and then call System.fit(). That's it! We'll turn on the logger so you can see how this adaptive procedure unfolds.
comp.clear() # Reset
system = System(comp, name='MF tutorial')
system.set_logger(stdout=True)
system.fit()
2024-11-05 19:26:28,714 — [INFO] — MF tutorial — -------------------Refining system surrogate: iteration 1-------------------
2024-11-05 19:26:28,714 — [INFO] — MF tutorial — Initializing component multilevel_model: adding ((0,), (0,)) to active set
2024-11-05 19:26:28,720 — [INFO] — multilevel_model — Running 4 total model evaluations for component 'multilevel_model' new candidate indices: [((0,), (0,)), ((1,), (0,)), ((0,), (1,))]...
2024-11-05 19:26:28,722 — [INFO] — MF tutorial — -------------------Refining system surrogate: iteration 2-------------------
2024-11-05 19:26:28,723 — [INFO] — MF tutorial — Estimating error for component 'multilevel_model'...
2024-11-05 19:26:28,725 — [INFO] — MF tutorial — Candidate multi-index: ((1,), (0,)). Relative error: 2554768691661083.5. Error indicator: 2554768691661083.5.
2024-11-05 19:26:28,726 — [INFO] — MF tutorial — Candidate multi-index: ((0,), (1,)). Relative error: 9618816470364806.0. Error indicator: 4809408235182403.0.
2024-11-05 19:26:28,726 — [INFO] — MF tutorial — Candidate multi-index ((0,), (1,)) chosen for component 'multilevel_model'.
2024-11-05 19:26:28,733 — [INFO] — multilevel_model — Running 2 total model evaluations for component 'multilevel_model' new candidate indices: [((0,), (2,))]...
2024-11-05 19:26:28,735 — [INFO] — MF tutorial — -------------------Refining system surrogate: iteration 3-------------------
2024-11-05 19:26:28,736 — [INFO] — MF tutorial — Estimating error for component 'multilevel_model'...
2024-11-05 19:26:28,739 — [INFO] — MF tutorial — Candidate multi-index: ((0,), (2,)). Relative error: 0.2641454462385164. Error indicator: 0.1320727231192582.
2024-11-05 19:26:28,740 — [INFO] — MF tutorial — Candidate multi-index: ((1,), (0,)). Relative error: 0.2967144262656313. Error indicator: 0.2967144262656313.
2024-11-05 19:26:28,740 — [INFO] — MF tutorial — Candidate multi-index ((1,), (0,)) chosen for component 'multilevel_model'.
2024-11-05 19:26:28,741 — [INFO] — multilevel_model — Running 3 total model evaluations for component 'multilevel_model' new candidate indices: [((1,), (1,)), ((2,), (0,))]...
2024-11-05 19:26:28,743 — [INFO] — MF tutorial — -------------------Refining system surrogate: iteration 4-------------------
2024-11-05 19:26:28,745 — [INFO] — MF tutorial — Estimating error for component 'multilevel_model'...
2024-11-05 19:26:28,749 — [INFO] — MF tutorial — Candidate multi-index: ((1,), (1,)). Relative error: 0.1150682011384518. Error indicator: 0.0575341005692259.
2024-11-05 19:26:28,750 — [INFO] — MF tutorial — Candidate multi-index: ((0,), (2,)). Relative error: 0.2756308582784701. Error indicator: 0.13781542913923506.
2024-11-05 19:26:28,750 — [INFO] — MF tutorial — Candidate multi-index: ((2,), (0,)). Relative error: 0.13703027278834307. Error indicator: 0.13703027278834307.
2024-11-05 19:26:28,753 — [INFO] — MF tutorial — Candidate multi-index ((0,), (2,)) chosen for component 'multilevel_model'.
2024-11-05 19:26:28,754 — [INFO] — MF tutorial — -------------------Refining system surrogate: iteration 5-------------------
2024-11-05 19:26:28,756 — [INFO] — MF tutorial — Estimating error for component 'multilevel_model'...
2024-11-05 19:26:28,760 — [INFO] — MF tutorial — Candidate multi-index: ((1,), (1,)). Relative error: 0.09912147321186997. Error indicator: 0.049560736605934984.
2024-11-05 19:26:28,760 — [INFO] — MF tutorial — Candidate multi-index: ((2,), (0,)). Relative error: 0.09999076332676624. Error indicator: 0.09999076332676624.
2024-11-05 19:26:28,760 — [INFO] — MF tutorial — Candidate multi-index ((2,), (0,)) chosen for component 'multilevel_model'.
2024-11-05 19:26:28,761 — [INFO] — MF tutorial — -------------------Refining system surrogate: iteration 6-------------------
2024-11-05 19:26:28,763 — [INFO] — MF tutorial — Estimating error for component 'multilevel_model'...
2024-11-05 19:26:28,766 — [INFO] — MF tutorial — Candidate multi-index: ((1,), (1,)). Relative error: 0.0974652862212737. Error indicator: 0.04873264311063685.
2024-11-05 19:26:28,766 — [INFO] — MF tutorial — Candidate multi-index ((1,), (1,)) chosen for component 'multilevel_model'.
2024-11-05 19:26:28,767 — [INFO] — multilevel_model — Running 4 total model evaluations for component 'multilevel_model' new candidate indices: [((2,), (1,)), ((1,), (2,))]...
2024-11-05 19:26:28,769 — [INFO] — MF tutorial — -------------------Refining system surrogate: iteration 7-------------------
2024-11-05 19:26:28,772 — [INFO] — MF tutorial — Estimating error for component 'multilevel_model'...
2024-11-05 19:26:28,776 — [INFO] — MF tutorial — Candidate multi-index: ((1,), (2,)). Relative error: 0.0077765588767778485. Error indicator: 0.0038882794383889242.
2024-11-05 19:26:28,776 — [INFO] — MF tutorial — Candidate multi-index: ((2,), (1,)). Relative error: 0.052929814502854895. Error indicator: 0.026464907251427448.
2024-11-05 19:26:28,777 — [INFO] — MF tutorial — Candidate multi-index ((2,), (1,)) chosen for component 'multilevel_model'.
2024-11-05 19:26:28,778 — [INFO] — MF tutorial — -------------------Refining system surrogate: iteration 8-------------------
2024-11-05 19:26:28,780 — [INFO] — MF tutorial — Estimating error for component 'multilevel_model'...
2024-11-05 19:26:28,782 — [INFO] — MF tutorial — Candidate multi-index: ((1,), (2,)). Relative error: 0.008452366215639355. Error indicator: 0.004226183107819677.
2024-11-05 19:26:28,783 — [INFO] — MF tutorial — Candidate multi-index ((1,), (2,)) chosen for component 'multilevel_model'.
2024-11-05 19:26:28,783 — [INFO] — multilevel_model — Running 2 total model evaluations for component 'multilevel_model' new candidate indices: [((2,), (2,))]...
2024-11-05 19:26:28,785 — [INFO] — MF tutorial — -------------------Refining system surrogate: iteration 9-------------------
2024-11-05 19:26:28,787 — [INFO] — MF tutorial — Estimating error for component 'multilevel_model'...
2024-11-05 19:26:28,789 — [INFO] — MF tutorial — Candidate multi-index: ((2,), (2,)). Relative error: 0.004665236286971056. Error indicator: 0.002332618143485528.
2024-11-05 19:26:28,790 — [INFO] — MF tutorial — Candidate multi-index ((2,), (2,)) chosen for component 'multilevel_model'.
2024-11-05 19:26:28,791 — [INFO] — MF tutorial — -------------------Refining system surrogate: iteration 10-------------------
2024-11-05 19:26:28,792 — [INFO] — MF tutorial — Estimating error for component 'multilevel_model'...
2024-11-05 19:26:28,792 — [INFO] — MF tutorial — Component 'multilevel_model' has no available candidates left!
2024-11-05 19:26:28,793 — [INFO] — MF tutorial — No candidates left for refinement, iteration: 9
2024-11-05 19:26:28,793 — [INFO] — MF tutorial — -----------------------------Termination criteria reached: No candidates left to refine-----------------------------
2024-11-05 19:26:28,794 — [INFO] — MF tutorial — Final system surrogate: ---- MF tutorial ---- amisc version: 0.5.2 Refinement level: 9 Components: multilevel_model Inputs: x Outputs: y
At this point, you may be wondering "what exactly is the surrogate and how exactly is it training?"
By default, we use the original multivariate Lagrange polynomials described in the AMISC journal paper. You can view this as the Component.interpolator property. Lagrange polynomials work well to interpolate smooth response functions up to an input dimension of around 12-15 inputs. The training data is selected by the Leja objective function in a sparse grid format, which is generally a "space-filling" design. You can view this as the Component.training_data property. While it is beyond the scope of this tutorial, these properties are designed to be extensible, so that different interpolation or experimental design strategies could be implemented on top of the amisc framework, such as kriging or latin hypercube sampling.
The main idea that amisc brings to the proverbial surrogate table is the ability to build multifidelity surrogates of multiple models linked in a multidisciplinary fashion -- the specific mathematical interpolation/surrogate technique is abstracted out and left up to the user.
Evaluate surrogate performance¶
Now that we know how to train a surrogate, we'd like to know how good the fit is. As an example, let's go back to our three-component system:
\begin{align} y_1 &= f_1(x) = x \sin(\pi x)\\ y_2 &= f_2(y_1) = \frac{1}{1 + 25y_1^2}\\ y_3 &= f_3(x, y_2) = x \cos(\pi y_2) \end{align}
def f1(inputs):
x = inputs['x']
return {'y1': x * np.sin(np.pi * x)}
def f2(inputs):
y1 = inputs['y1']
return {'y2': 1 / (1 + 25*y1**2)}
def f3(inputs):
x = inputs['x']
y2 = inputs['y2']
return {'y3': x * np.cos(np.pi * y2)}
md_system = System(Component(f1, Variable('x', domain=(0, 1)), Variable('y1'), data_fidelity=(2,)),
Component(f2, Variable('y1', domain=(0, 1)), Variable('y2'), data_fidelity=(2,)),
Component(f3, ['x', Variable('y2', domain=(0, 1))], Variable('y3'), data_fidelity=(2, 2)),
name='3-component system'
)
Note how, for the default Component.training_data = SparseGrid, the data fidelity must be an indication of the amount of training data in each input dimension. As a result, len(data_fidelity) = len(inputs) for each component. For details, look into the amisc.training module.
Also note that models may not have any $\alpha$ fidelity indices -- this is completely fine! The surrogate will be built over any $\beta$ indices instead. However, if you do not provide $\alpha$ or $\beta$ (i.e. model_fidelity, data_fidelity, or surrogate_fidelity) then system.predict() will just call the underlying models directly, since no surrogate can be built in the absence of fidelity indices.
Now, let's run 10 training iterations.
md_system.set_logger(stdout=True)
md_system.fit(max_iter=10)
2024-11-05 19:26:28,892 — [INFO] — 3-component system — -------------------Refining system surrogate: iteration 1-------------------
2024-11-05 19:26:28,893 — [INFO] — 3-component system — Initializing component f1: adding ((), (0,)) to active set
2024-11-05 19:26:28,899 — [INFO] — f1 — Running 3 total model evaluations for component 'f1' new candidate indices: [((), (0,)), ((), (1,))]...
2024-11-05 19:26:28,900 — [INFO] — 3-component system — -------------------Refining system surrogate: iteration 2-------------------
2024-11-05 19:26:28,901 — [INFO] — 3-component system — Initializing component f2: adding ((), (0,)) to active set
2024-11-05 19:26:28,906 — [INFO] — f2 — Running 3 total model evaluations for component 'f2' new candidate indices: [((), (0,)), ((), (1,))]...
2024-11-05 19:26:28,908 — [INFO] — 3-component system — -------------------Refining system surrogate: iteration 3-------------------
2024-11-05 19:26:28,908 — [INFO] — 3-component system — Initializing component f3: adding ((), (0, 0)) to active set
2024-11-05 19:26:28,920 — [INFO] — f3 — Running 5 total model evaluations for component 'f3' new candidate indices: [((), (0, 0)), ((), (0, 1)), ((), (1, 0))]...
2024-11-05 19:26:28,922 — [INFO] — 3-component system — -------------------Refining system surrogate: iteration 4-------------------
2024-11-05 19:26:28,924 — [INFO] — 3-component system — Estimating error for component 'f1'...
2024-11-05 19:26:28,926 — [INFO] — 3-component system — Candidate multi-index: ((), (1,)). Relative error: 0.48242756879392457. Error indicator: 0.24121378439696228.
2024-11-05 19:26:28,927 — [INFO] — 3-component system — Estimating error for component 'f2'...
2024-11-05 19:26:28,929 — [INFO] — 3-component system — Candidate multi-index: ((), (1,)). Relative error: 0.0. Error indicator: 0.0.
2024-11-05 19:26:28,930 — [INFO] — 3-component system — Estimating error for component 'f3'...
2024-11-05 19:26:28,934 — [INFO] — 3-component system — Candidate multi-index: ((), (0, 1)). Relative error: 1.182613049627404e+16. Error indicator: 5913065248137020.0.
2024-11-05 19:26:28,934 — [INFO] — 3-component system — Candidate multi-index: ((), (1, 0)). Relative error: 0.619913905352461. Error indicator: 0.3099569526762305.
2024-11-05 19:26:28,936 — [INFO] — 3-component system — Candidate multi-index ((), (0, 1)) chosen for component 'f3'.
2024-11-05 19:26:28,943 — [INFO] — f3 — Running 2 total model evaluations for component 'f3' new candidate indices: [((), (0, 2))]...
2024-11-05 19:26:28,944 — [INFO] — 3-component system — -------------------Refining system surrogate: iteration 5-------------------
2024-11-05 19:26:28,947 — [INFO] — 3-component system — Estimating error for component 'f1'...
2024-11-05 19:26:28,949 — [INFO] — 3-component system — Candidate multi-index: ((), (1,)). Relative error: 0.4436664467953997. Error indicator: 0.22183322339769984.
2024-11-05 19:26:28,949 — [INFO] — 3-component system — Estimating error for component 'f2'...
2024-11-05 19:26:28,952 — [INFO] — 3-component system — Candidate multi-index: ((), (1,)). Relative error: 0.0. Error indicator: 0.0.
2024-11-05 19:26:28,952 — [INFO] — 3-component system — Estimating error for component 'f3'...
2024-11-05 19:26:28,957 — [INFO] — 3-component system — Candidate multi-index: ((), (1, 0)). Relative error: 1.0510796693677949e-16. Error indicator: 5.2553983468389743e-17.
2024-11-05 19:26:28,958 — [INFO] — 3-component system — Candidate multi-index: ((), (0, 2)). Relative error: 0.2561815692839554. Error indicator: 0.1280907846419777.
2024-11-05 19:26:28,958 — [INFO] — 3-component system — Candidate multi-index ((), (1,)) chosen for component 'f1'.
2024-11-05 19:26:28,965 — [INFO] — f1 — Running 2 total model evaluations for component 'f1' new candidate indices: [((), (2,))]...
2024-11-05 19:26:28,967 — [INFO] — 3-component system — -------------------Refining system surrogate: iteration 6-------------------
2024-11-05 19:26:28,969 — [INFO] — 3-component system — Estimating error for component 'f1'...
2024-11-05 19:26:28,973 — [INFO] — 3-component system — Candidate multi-index: ((), (2,)). Relative error: 0.33701429103278285. Error indicator: 0.16850714551639143.
2024-11-05 19:26:28,973 — [INFO] — 3-component system — Estimating error for component 'f2'...
2024-11-05 19:26:28,976 — [INFO] — 3-component system — Candidate multi-index: ((), (1,)). Relative error: 2.1688310838232674. Error indicator: 1.0844155419116337.
2024-11-05 19:26:28,976 — [INFO] — 3-component system — Estimating error for component 'f3'...
2024-11-05 19:26:28,981 — [INFO] — 3-component system — Candidate multi-index: ((), (1, 0)). Relative error: 1.1055754303262678e-16. Error indicator: 5.527877151631339e-17.
2024-11-05 19:26:28,982 — [INFO] — 3-component system — Candidate multi-index: ((), (0, 2)). Relative error: 0.2561815692839554. Error indicator: 0.1280907846419777.
2024-11-05 19:26:28,983 — [INFO] — 3-component system — Candidate multi-index ((), (1,)) chosen for component 'f2'.
2024-11-05 19:26:28,990 — [INFO] — f2 — Running 2 total model evaluations for component 'f2' new candidate indices: [((), (2,))]...
2024-11-05 19:26:28,991 — [INFO] — 3-component system — -------------------Refining system surrogate: iteration 7-------------------
2024-11-05 19:26:28,994 — [INFO] — 3-component system — Estimating error for component 'f1'...
2024-11-05 19:26:28,996 — [INFO] — 3-component system — Candidate multi-index: ((), (2,)). Relative error: 0.7939429646333563. Error indicator: 0.39697148231667817.
2024-11-05 19:26:28,997 — [INFO] — 3-component system — Estimating error for component 'f2'...
2024-11-05 19:26:29,000 — [INFO] — 3-component system — Candidate multi-index: ((), (2,)). Relative error: 0.10486699316077953. Error indicator: 0.05243349658038977.
2024-11-05 19:26:29,000 — [INFO] — 3-component system — Estimating error for component 'f3'...
2024-11-05 19:26:29,005 — [INFO] — 3-component system — Candidate multi-index: ((), (1, 0)). Relative error: 1.20439733802513e-16. Error indicator: 6.02198669012565e-17.
2024-11-05 19:26:29,006 — [INFO] — 3-component system — Candidate multi-index: ((), (0, 2)). Relative error: 0.3237545298072328. Error indicator: 0.1618772649036164.
2024-11-05 19:26:29,007 — [INFO] — 3-component system — Candidate multi-index ((), (2,)) chosen for component 'f1'.
2024-11-05 19:26:29,007 — [INFO] — 3-component system — -------------------Refining system surrogate: iteration 8-------------------
2024-11-05 19:26:29,010 — [INFO] — 3-component system — Estimating error for component 'f1'...
2024-11-05 19:26:29,011 — [INFO] — 3-component system — Component 'f1' has no available candidates left!
2024-11-05 19:26:29,011 — [INFO] — 3-component system — Estimating error for component 'f2'...
2024-11-05 19:26:29,014 — [INFO] — 3-component system — Candidate multi-index: ((), (2,)). Relative error: 0.07781981328850017. Error indicator: 0.038909906644250085.
2024-11-05 19:26:29,015 — [INFO] — 3-component system — Estimating error for component 'f3'...
2024-11-05 19:26:29,020 — [INFO] — 3-component system — Candidate multi-index: ((), (1, 0)). Relative error: 7.494890520653828e-17. Error indicator: 3.747445260326914e-17.
2024-11-05 19:26:29,021 — [INFO] — 3-component system — Candidate multi-index: ((), (0, 2)). Relative error: 0.22720071642637238. Error indicator: 0.11360035821318619.
2024-11-05 19:26:29,021 — [INFO] — 3-component system — Candidate multi-index ((), (0, 2)) chosen for component 'f3'.
2024-11-05 19:26:29,022 — [INFO] — 3-component system — -------------------Refining system surrogate: iteration 9-------------------
2024-11-05 19:26:29,025 — [INFO] — 3-component system — Estimating error for component 'f1'...
2024-11-05 19:26:29,025 — [INFO] — 3-component system — Component 'f1' has no available candidates left!
2024-11-05 19:26:29,026 — [INFO] — 3-component system — Estimating error for component 'f2'...
2024-11-05 19:26:29,029 — [INFO] — 3-component system — Candidate multi-index: ((), (2,)). Relative error: 0.08268189244193577. Error indicator: 0.041340946220967886.
2024-11-05 19:26:29,029 — [INFO] — 3-component system — Estimating error for component 'f3'...
2024-11-05 19:26:29,033 — [INFO] — 3-component system — Candidate multi-index: ((), (1, 0)). Relative error: 4.3898253753650394e-17. Error indicator: 2.1949126876825197e-17.
2024-11-05 19:26:29,034 — [INFO] — 3-component system — Candidate multi-index ((), (2,)) chosen for component 'f2'.
2024-11-05 19:26:29,034 — [INFO] — 3-component system — -------------------Refining system surrogate: iteration 10-------------------
2024-11-05 19:26:29,038 — [INFO] — 3-component system — Estimating error for component 'f1'...
2024-11-05 19:26:29,038 — [INFO] — 3-component system — Component 'f1' has no available candidates left!
2024-11-05 19:26:29,039 — [INFO] — 3-component system — Estimating error for component 'f2'...
2024-11-05 19:26:29,039 — [INFO] — 3-component system — Component 'f2' has no available candidates left!
2024-11-05 19:26:29,040 — [INFO] — 3-component system — Estimating error for component 'f3'...
2024-11-05 19:26:29,043 — [INFO] — 3-component system — Candidate multi-index: ((), (1, 0)). Relative error: 4.346287365165632e-17. Error indicator: 2.173143682582816e-17.
2024-11-05 19:26:29,044 — [INFO] — 3-component system — Candidate multi-index ((), (1, 0)) chosen for component 'f3'.
2024-11-05 19:26:29,051 — [INFO] — f3 — Running 6 total model evaluations for component 'f3' new candidate indices: [((), (1, 1)), ((), (2, 0))]...
2024-11-05 19:26:29,053 — [INFO] — 3-component system — ------------------------Termination criteria reached: Max iteration 10/10------------------------
2024-11-05 19:26:29,054 — [INFO] — 3-component system — Final system surrogate: ---- 3-component system ---- amisc version: 0.5.2 Refinement level: 10 Components: f1, f2, f3 Inputs: x Outputs: y3, y2, y1
To determine the performance of the surrogate, let's randomly sample the inputs and compare the surrogate predictions to the ground truth model. For convenience, you can evaluate the underlying model(s) using the same predict() function with the use_model flag.
xtest = md_system.sample_inputs(100)
ytest_surr = md_system.predict(xtest)
ytest_model = md_system.predict(xtest, use_model='best')
from amisc.utils import relative_error
for output in ytest_surr:
error = relative_error(ytest_surr[output], ytest_model[output])
print(f'Relative error for {output} = {error:.3f}')
Relative error for y1 = 0.017 Relative error for y2 = 0.059 Relative error for y3 = 0.109
Not bad! We get around 2% and 5% error for $y_1$ and $y_2$. Since $y_3$ is fit over both $x$ and $y_2$, it is a bit more complicated and may require more training to improve.
Another useful diagnostic is to plot the system outputs as functions of the inputs:
fig, ax = md_system.plot_slice(show_model='best', random_walk=True)
plt.show()
These plots can only visualize 1d "slices", which provides limited information in higher input dimensions. However, it does provide a good indication of how "smooth" the output response is, which may help decide if Lagrange polynomials are sufficient or if some other surrogate method should be used.
Saving and loading from file¶
After fitting a surrogate, you'll want to save the results for later use (especially if you fit a fairly expensive model). We provide two utilities for data persistence: save_to_file and load_from_file.
# Save to a temporary directory
import tempfile
import pathlib
with tempfile.TemporaryDirectory() as tmp_path:
md_system.save_to_file('my_system.yml', save_dir=tmp_path)
# Later on
loaded_system = System.load_from_file(pathlib.Path(tmp_path) / 'my_system.yml')
These utility functions use the amisc.YamlLoader class by default to save/load the system from a YAML (.yml) file. Yaml files are very similar text-based formats as .json, but they integrate better with Python objects. You can open a save.yml file in any text-editor to view (and even edit) the amisc objects that are saved there.
Important: To consistently save/load callable model functions, the model should be located in some global scope (i.e. no lambda or local functions). This follows the same rules as Python's built-in pickle module -- it only saves a reference to the function's import path, which must be available when you load back from file. If you move or rename the model functions, you can still use the same .yml save file! You'll just have to manually edit the save file to point the models to their new locations.
One last note about files -- you can also completely configure an entire amisc.System object from a .yml config file, without having to define anything manually in Python. This is achievable through the use of the three built-in YAML "tags" for variables, components, and systems. A "tag" is just an exclamation followed by the name of the object, (e.g. !Variable).
Our three-component system could be completely defined in a YAML file like this:
# Loading from file
yaml_config = """
!System
name: 3-component system
components:
!Component
- model: !!python/name:amisc.examples.models.f1
data_fidelity: (2,)
inputs: !Variable
- name: x
domain: (0, 1)
outputs: !Variable
- name: y1
call_unpacked: true # not reading from a `dict` input (just to silence some warnings)
ret_unpacked: true # not returning a `dict` output
- model: !!python/name:amisc.examples.models.f2
data_fidelity: (2,)
inputs: !Variable
- name: y1
domain: (0, 1)
outputs: !Variable
- name: y2
call_unpacked: true
ret_unpacked: true
- model: !!python/name:amisc.examples.models.f3
data_fidelity: (2,)
inputs: !Variable
- name: x
- name: y2
domain: (0, 1)
outputs: !Variable
- name: y3
call_unpacked: true
ret_unpacked: true
"""
Note the use of the !!python/name tag to indicate the importable path at which to find the respective component models -- we provide these three functions in the amisc.examples.models module for convenience. For you, one could make a models.py file in the current directory, for example, and then write some function def my_model(inputs): .... You would then reference this function via yaml as:
!!python/name:models.my_model
which will work so long as Python can find the models.py file on its search path (which always includes the current directory).
To verify our config file, let's load it using YamlLoader:
from amisc import YamlLoader
md_system = YamlLoader.load(yaml_config)
print(md_system)
---- 3-component system ---- amisc version: 0.5.2 Refinement level: 0 Components: f1, f2, f3 Inputs: x Outputs: y3, y2, y1
Advanced features¶
This completes the intro tutorial! If you have decided amisc might be useful for you, then you can view the online docs for detailed API reference and other advanced features.
These include:
- Using random variables
- Normalizing inputs/outputs
- Fitting a surrogate for high-dimensional field-quantities
- Handling feedback between models
- Parallelizing the models and training
- Extending
amisc