diff --git a/changelog/fragments/1166.feat.rst b/changelog/fragments/1166.feat.rst new file mode 100644 index 0000000000000000000000000000000000000000..a795fed06bb88affc33a536d55ef80479ad58989 --- /dev/null +++ b/changelog/fragments/1166.feat.rst @@ -0,0 +1 @@ +Jacobian can be approximated hybridly by defining certain elements of the Jacobian within the discipline and then calling the ``linearize`` method and passing one of the hybrid Approximation modes: ``HYBRID_FINITE_DIFFERENCES``, ``HYBRID_CENTERED_DIFFERENCES``, ``HYBRID_COMPLEX_STEP`` diff --git a/src/gemseo/core/discipline/discipline.py b/src/gemseo/core/discipline/discipline.py index b666ffc8bbb59bc62084ee6128011356237a0647..4631af665d9d1392698edc2f2fa0570a2bc555d1 100644 --- a/src/gemseo/core/discipline/discipline.py +++ b/src/gemseo/core/discipline/discipline.py @@ -32,6 +32,7 @@ from gemseo.core.derivatives.derivation_modes import DerivationMode from gemseo.core.discipline.base_discipline import BaseDiscipline from gemseo.utils.constants import READ_ONLY_EMPTY_DICT from gemseo.utils.derivatives.approximation_modes import ApproximationMode +from gemseo.utils.derivatives.approximation_modes import HybridApproximationMode from gemseo.utils.derivatives.derivatives_approx import DisciplineJacApprox from gemseo.utils.derivatives.error_estimators import EPSILON from gemseo.utils.enumeration import merge_enums @@ -85,7 +86,9 @@ class Discipline(BaseDiscipline, metaclass=ClassInjector): SPARSE = "sparse" """Initialized as SciPy CSR arrays filled with zeros.""" - ApproximationMode: EnumType = ApproximationMode + ApproximationMode: EnumType = merge_enums( + "ApproximationMode", StrEnum, ApproximationMode, HybridApproximationMode + ) LinearizationMode: EnumType = merge_enums( "LinearizationMode", @@ -124,6 +127,14 @@ class Discipline(BaseDiscipline, metaclass=ClassInjector): __output_names: Iterable[str] """The output names used for handling execution status and statistics.""" + __hybrid_approximation_name_to_mode: ClassVar[ + Mapping[ApproximationMode, ApproximationMode] + ] = { + ApproximationMode.HYBRID_FINITE_DIFFERENCES: ApproximationMode.FINITE_DIFFERENCES, # noqa: E501 + ApproximationMode.HYBRID_CENTERED_DIFFERENCES: ApproximationMode.CENTERED_DIFFERENCES, # noqa: E501 + ApproximationMode.HYBRID_COMPLEX_STEP: ApproximationMode.COMPLEX_STEP, + } + def __init__( # noqa: D107 self, name: str = "", @@ -252,9 +263,12 @@ class Discipline(BaseDiscipline, metaclass=ClassInjector): def __compute_jacobian(self): """Callable used for handling execution status and statistics.""" if self._linearization_mode in set(self.ApproximationMode): - self.jac = self._jac_approx.compute_approx_jac( - self.__output_names, self.__input_names - ) + if self._linearization_mode in set(HybridApproximationMode): + self._compose_hybrid_jacobian(self.__output_names, self.__input_names) + else: + self.jac = self._jac_approx.compute_approx_jac( + self.__output_names, self.__input_names + ) else: self._compute_jacobian(self.__input_names, self.__output_names) @@ -289,10 +303,15 @@ class Discipline(BaseDiscipline, metaclass=ClassInjector): jac_approx_wait_time: The time waited between two forks of the process / thread. """ + approx_method = ( + self.__hybrid_approximation_name_to_mode[jac_approx_type] + if jac_approx_type in set(HybridApproximationMode) + else jac_approx_type + ) self._jac_approx = DisciplineJacApprox( # TODO: pass the bare minimum instead of self. self, - approx_method=jac_approx_type, + approx_method=approx_method, step=jax_approx_step, parallel=jac_approx_n_processes > 1, n_processes=jac_approx_n_processes, @@ -807,3 +826,60 @@ class Discipline(BaseDiscipline, metaclass=ClassInjector): output_names: The names of the outputs to be differentiated. If empty, use all the outputs. """ + + def _compose_hybrid_jacobian( + self, + output_names: Iterable[str] = (), + input_names: Iterable[str] = (), + ) -> None: + """Compose a hybrid Jacobian using analytical and approximated expressions. + + This method allows to complete a given Jacobian for which not all + inputs-outputs have been defined by using approximation methods. + + Args: + input_names: The names of the input wrt the ``output_names`` are + linearized. + output_names: The names of the output to be linearized. + """ + analytical_jacobian = self.jac + + self.set_jacobian_approximation(self._linearization_mode) + + approximated_jac = {} + outputs_names_to_approximate = [] + input_names_to_approximate = [] + + for output_name in output_names: + if output_name not in analytical_jacobian: + analytical_jacobian[output_name] = {} + for input_name in input_names: + if input_name not in analytical_jacobian[output_name]: + approximated_jac[output_name] = {} + # Map the outputs to be differentiated wrt the corresponding inputs. + outputs_names_to_approximate.append(output_name) + input_names_to_approximate.append(input_name) + + # Compute approximated Jacobian elements. + for output_name, input_name in zip( + outputs_names_to_approximate, input_names_to_approximate + ): + jac_input_output = self._jac_approx.compute_approx_jac( + [output_name], [input_name] + ) + approximated_jac[output_name][input_name] = jac_input_output[output_name][ + input_name + ] + + # Fill in missing inputs of the Jacobian. + for output_name in output_names: + analytical_jacobian_out = analytical_jacobian[output_name] + approximated_jac_out = approximated_jac[output_name] + for input_name in input_names: + if input_name not in analytical_jacobian_out: + analytical_jacobian_out[input_name] = approximated_jac_out[ + input_name + ] + + # Recover analytical Jacobian data. + self.jac = analytical_jacobian diff --git a/src/gemseo/utils/derivatives/approximation_modes.py b/src/gemseo/utils/derivatives/approximation_modes.py index c00a3f314d8bcfb200d5b6ced2160b46c751b3cb..298c224deaddc357d69f7116f6c7a59880e794dc 100644 --- a/src/gemseo/utils/derivatives/approximation_modes.py +++ b/src/gemseo/utils/derivatives/approximation_modes.py @@ -33,3 +33,19 @@ class ApproximationMode(StrEnum): CENTERED_DIFFERENCES = "centered_differences" """The centered differences method used to approximate the Jacobians by perturbing each variable with a small real number.""" + + +class HybridApproximationMode(StrEnum): + """The approximation derivation modes for semi-analytical computations.""" + + HYBRID_COMPLEX_STEP = "hybrid_complex_step" + """"The complex step method used to approximiate the Jacobian not available + analytically by perturbing those variables with a small complex number.""" + + HYBRID_FINITE_DIFFERENCES = "hybrid_finite_differences" + """"The finite differences method used to approxiamate the Jacobian not available + analytically by perturbing those variables with a small real number.""" + + HYBRID_CENTERED_DIFFERENCES = "hybrid_centered_differences" + """The centered differences method used to approxiamate the Jacobian not available + analytically by perturbing those variables with a small real number.""" diff --git a/tests/core/test_discipline.py b/tests/core/test_discipline.py index cb08b1d947ad3dedaa7f07dcba549b825354453a..38601c6a41f0d4581e427ff7499d058b3b3e273e 100644 --- a/tests/core/test_discipline.py +++ b/tests/core/test_discipline.py @@ -68,6 +68,8 @@ from gemseo.utils.pickle import to_pickle from gemseo.utils.repr_html import REPR_HTML_WRAPPER if TYPE_CHECKING: + from collections.abc import Iterable + from gemseo.typing import StrKeyMapping Status = ExecutionStatus.Status @@ -117,6 +119,44 @@ def sobieski_chain() -> tuple[MDOChain, dict[str, ndarray]]: return chain, indata +@pytest.fixture +def hybrid_jacobian_discipline() -> Discipline: + class HybridDiscipline(Discipline): + def __init__(self) -> None: + super().__init__() + self.input_grammar.update_from_names(["x_1"]) + self.input_grammar.update_from_names(["x_2"]) + self.input_grammar.update_from_names(["x_3"]) + self.output_grammar.update_from_names(["y_1"]) + self.output_grammar.update_from_names(["y_2"]) + self.output_grammar.update_from_names(["y_3"]) + self.default_input_data = { + "x_1": array([1.0]), + "x_2": array([2.0]), + "x_3": array([1.0]), + } + + def _run(self, input_data: StrKeyMapping) -> StrKeyMapping | None: + self.io.data["y_1"] = input_data["x_1"] * input_data["x_2"] + self.io.data["y_2"] = ( + input_data["x_1"] * input_data["x_2"] * input_data["x_3"] + ) + self.io.data["y_3"] = input_data["x_1"] + + def _compute_jacobian( + self, + input_names: Iterable[str] = (), + output_names: Iterable[str] = (), + ) -> None: + self._init_jacobian() + x1 = array([self.get_input_data(with_namespaces=False)["x_1"]]) + x2 = array([self.get_input_data(with_namespaces=False)["x_2"]]) + x3 = array([self.get_input_data(with_namespaces=False)["x_3"]]) + self.jac = {"y_1": {"x_1": x2}, "y_2": {"x_2": x1 * x3}} + + return HybridDiscipline() + + @pytest.mark.xfail def test_instantiate_grammars() -> None: """Test the instantiation of the grammars.""" @@ -198,7 +238,6 @@ def test_check_jac_fdapprox() -> None: aero.linearization_mode = aero.ApproximationMode.FINITE_DIFFERENCES aero.linearize(inpts, compute_all_jacobians=True) aero.check_jacobian(inpts) - aero.linearization_mode = "auto" aero.check_jacobian(inpts) @@ -211,6 +250,25 @@ def test_check_jac_csapprox() -> None: aero.check_jacobian() +@pytest.mark.parametrize( + "hybrid_approximation_mode", + ["hybrid_complex_step", "hybrid_finite_differences", "hybrid_centered_differences"], +) +def test_check_jac_hybrid_approx( + hybrid_jacobian_discipline, hybrid_approximation_mode +) -> None: + """Test the hybrid finite difference approximation.""" + + disc = hybrid_jacobian_discipline + + inputs = disc.default_input_data + disc.linearization_mode = hybrid_approximation_mode + disc.linearize(inputs, compute_all_jacobians=True) + disc.check_jacobian( + inputs, linearization_mode=disc.ApproximationMode.FINITE_DIFFERENCES + ) + + def test_check_jac_approx_plot(tmp_wd) -> None: """Test the generation of the gradient plot.""" aero = SobieskiAerodynamics()