Compute the Jacobian of a discipline with finite differences

In this example, we will compute the Jacobians of some outputs of an MDODiscipline with respect to some inputs using the finite difference method.

from __future__ import annotations

from numpy import array

from gemseo.disciplines.auto_py import AutoPyDiscipline

First, we create a discipline, e.g. an AutoPyDiscipline:

def f(a=0.0, b=0.0):
    y = a**2 + b
    z = a**3 + b**2
    return y, z


discipline = AutoPyDiscipline(f)

We can execute it with its default input values:

discipline.execute()
discipline.local_data

{'b': array([0.]), 'a': array([0.]), 'z': array([0.]), 'y': array([0.])}

or with custom ones:

discipline.execute({"a": array([1.0])})
discipline.local_data

{'a': array([1.]), 'b': array([0.]), 'z': array([1.]), 'y': array([1.])}

Then, we use the method MDODiscipline.linearize() to compute the derivatives:

jacobian_data = discipline.linearize()
jacobian_data

{}

There is no Jacobian data because we need to set the input variables against which to compute the Jacobian of the output ones. For that, we use the method add_differentiated_inputs(). We also need to set these output variables: with the method add_differentiated_outputs(). For instance, we may want to compute the derivative of "z" with respect to "a" only:

discipline.add_differentiated_inputs(["a"])
discipline.add_differentiated_outputs(["z"])
jacobian_data = discipline.linearize()
jacobian_data

{'z': {'a': array([[1.e-14]])}}

We can have a quick look at the values of these derivatives and verify that they are equal to the analytical results, up to the numerical precision.

By default, GEMSEO uses MDODiscipline.default_inputs as input data for which to compute the Jacobian one. We can change them with input_data:

jacobian_data = discipline.linearize({"a": array([1.0])})
jacobian_data

{'z': {'a': array([[3.0000003]])}}

We can also force the discipline to compute the derivatives of all the outputs with respect to all the inputs:

jacobian_data = discipline.linearize(compute_all_jacobians=True)
jacobian_data

{'y': {'a': array([[1.e-07]]), 'b': array([[1.]])}, 'z': {'a': array([[1.e-14]]), 'b': array([[1.e-07]])}}

we can change the approximation type to complex step and compare the results:

discipline.set_jacobian_approximation(
    jac_approx_type=discipline.ApproximationMode.COMPLEX_STEP
)
jacobian_data_complex_step = discipline.linearize(compute_all_jacobians=True)
jacobian_data_complex_step

/home/docs/checkouts/readthedocs.org/user_builds/gemseo/envs/5.3.2/lib/python3.9/site-packages/gemseo/disciplines/auto_py.py:384: ComplexWarning: Casting complex values to real discards the imaginary part
  processed_data[key] = float(val[0])

{'y': {'a': array([[0.]]), 'b': array([[0.]])}, 'z': {'a': array([[0.]]), 'b': array([[0.]])}}

Lastly, We can change the approximation type to centered differences and compare the results:

discipline.set_jacobian_approximation(
    jac_approx_type=discipline.ApproximationMode.CENTERED_DIFFERENCES
)
jacobian_data_centered_differences = discipline.linearize(compute_all_jacobians=True)
jacobian_data_centered_differences

{'y': {'a': array([[0.]]), 'b': array([[1.]])}, 'z': {'a': array([[1.e-14]]), 'b': array([[0.]])}}