From c573a02f35ee6381b7d8abdb9a8c1398674a884c Mon Sep 17 00:00:00 2001 From: Jean-Christophe Giret Date: Tue, 25 Apr 2023 17:48:18 +0200 Subject: [PATCH 1/4] feat(MDA): filter the non-numeric array for assembling the MDA residual vector --- .../core/derivatives/jacobian_assembly.py | 10 ++ src/gemseo/mda/jacobi.py | 16 +- src/gemseo/mda/mda.py | 71 +++++++-- src/gemseo/mda/newton.py | 5 +- tests/mda/test_mda.py | 143 +++++++++++++++++- 5 files changed, 222 insertions(+), 23 deletions(-) diff --git a/src/gemseo/core/derivatives/jacobian_assembly.py b/src/gemseo/core/derivatives/jacobian_assembly.py index 375f7c488e..0c81da7c6a 100644 --- a/src/gemseo/core/derivatives/jacobian_assembly.py +++ b/src/gemseo/core/derivatives/jacobian_assembly.py @@ -642,6 +642,16 @@ class JacobianAssembly: self.coupling_structure, ) + # Exclude the non numeric couplings from the coupling minimal list + non_numeric_couplings = set() + for discipline in self.coupling_structure.disciplines: + inputs = discipline.get_input_data_names() + inputs_couplings = set(inputs) & set(couplings_minimal) + for coupling in inputs_couplings: + if not discipline.input_grammar.is_array(coupling, numeric_only=True): + non_numeric_couplings.add(coupling) + couplings_minimal = sorted(list(set(couplings_minimal) - non_numeric_couplings)) + couplings_and_res = couplings_minimal.copy() couplings_and_states = couplings_minimal.copy() # linearize all the disciplines diff --git a/src/gemseo/mda/jacobi.py b/src/gemseo/mda/jacobi.py index b7fda22a84..f231f38a0b 100644 --- a/src/gemseo/mda/jacobi.py +++ b/src/gemseo/mda/jacobi.py @@ -142,7 +142,9 @@ class MDAJacobi(MDA): inputs = self.get_input_data_names() strong_cpl = self.coupling_structure.all_couplings - self._input_couplings = set(strong_cpl) & set(inputs) + self._input_couplings = set(strong_cpl) & set(inputs) - set( + self._non_numeric_array_variables + ) def execute_all_disciplines( self, @@ -207,10 +209,10 @@ class MDAJacobi(MDA): new_couplings, log_normed_residual=self._log_convergence, ) - x_np1 = self._compute_nex_iterate(current_couplings, new_couplings) + x_np1 = self._compute_next_iterate(current_couplings, new_couplings) current_couplings = x_np1 - def _compute_nex_iterate( + def _compute_next_iterate( self, current_couplings: ndarray, new_couplings: ndarray, @@ -236,8 +238,14 @@ class MDAJacobi(MDA): self._dx_n.append(new_couplings - current_couplings) self._g_x_n.append(new_couplings) coupl_names = self._input_couplings + if self.sizes is None: - self.sizes = {key: self.local_data[key].size for key in coupl_names} + self.sizes = { + key: self.local_data[key].size + for key in coupl_names + if self.input_grammar.is_array(key, numeric_only=True) + } + dxn = self._dx_n[-1] dxn_1 = self._dx_n[-2] g_n = self._g_x_n[-1] diff --git a/src/gemseo/mda/mda.py b/src/gemseo/mda/mda.py index bcaeee4878..013c48de1b 100644 --- a/src/gemseo/mda/mda.py +++ b/src/gemseo/mda/mda.py @@ -16,6 +16,7 @@ # INITIAL AUTHORS - API and implementation and/or documentation # :author: Francois Gallard, Charlie Vanaret # OTHER AUTHORS - MACROSCOPIC CHANGES +# :author: Vincent Ambert """Base class for all Multi-disciplinary Design Analyses (MDA).""" from __future__ import annotations @@ -284,18 +285,32 @@ class MDA(MDODiscipline): self._input_couplings = list(input_couplings) def _current_input_couplings(self) -> ndarray: - """Return the current values of the input coupling variables.""" - input_couplings = list(self.get_outputs_by_name(self._input_couplings)) - if not input_couplings: + """Return the current values of the input coupling numeric variables.""" + numeric_input_couplings_names = [ + coupling + for coupling in self._input_couplings + if self.input_grammar.is_array(coupling, numeric_only=True) + ] + numeric_input_couplings = list( + self.get_outputs_by_name(numeric_input_couplings_names) + ) + if not numeric_input_couplings: return array([]) - return concatenate(input_couplings) + return concatenate(numeric_input_couplings) def _current_strong_couplings(self) -> ndarray: - """Return the current values of the strong coupling variables.""" - couplings = list(self.get_outputs_by_name(self.strong_couplings)) - if not couplings: + """Return the current values of the strong coupling numeric variables.""" + numeric_strong_couplings_names = [ + strong_coupling + for strong_coupling in self.strong_couplings + if self.output_grammar.is_array(strong_coupling, numeric_only=True) + ] + numeric_strong_couplings = list( + self.get_outputs_by_name(numeric_strong_couplings_names) + ) + if not numeric_strong_couplings: return array([]) - return concatenate(couplings) + return concatenate(numeric_strong_couplings) def _retrieve_diff_inouts( self, compute_all_jacobians: bool = False @@ -315,6 +330,18 @@ class MDA(MDODiscipline): if self.RESIDUALS_NORM in outputs: outputs = list(outputs) outputs.remove(self.RESIDUALS_NORM) + + # Filter the non-numeric arrays + inputs = [ + input + for input in inputs + if self.input_grammar.is_array(input, numeric_only=True) + ] + outputs = [ + output + for output in outputs + if self.output_grammar.is_array(output, numeric_only=True) + ] return inputs, outputs def _couplings_warm_start(self) -> None: @@ -330,8 +357,10 @@ class MDA(MDODiscipline): def _check_coupling_types(self) -> None: """Check that the coupling variables are of type array in the grammars. - Raises: - TypeError: When at least one of the coupling variables is not an array. + If non-array coupling variables are present, they will be filtered and not taken + into account in the MDA residual. Yet, this method warns the user that some of + the coupling variables are non-numeric array, in case of this event follows an + improper setup. """ not_arrays = [] for discipline in self.disciplines: @@ -344,9 +373,11 @@ class MDA(MDODiscipline): if not_arrays: not_arrays = sorted(set(not_arrays)) - raise TypeError( - f"The coupling variables {not_arrays} must be of type array." + LOGGER.warning( + "The coupling variable(s) %s is/are not an array of numeric values.", + not_arrays, ) + self._non_numeric_array_variables = not_arrays def reset_disciplines_statuses(self) -> None: """Reset all the statuses of the disciplines.""" @@ -388,10 +419,14 @@ class MDA(MDODiscipline): couplings_adjoint = sorted( set(self.all_couplings) + - set(self._non_numeric_array_variables) - residual_variables.keys() - set(residual_variables.values()) ) + outputs = list(set(outputs) - set(self._non_numeric_array_variables)) + inputs = list(set(inputs) - set(self._non_numeric_array_variables)) + self.jac = self.assembly.total_derivatives( self.local_data, outputs, @@ -572,6 +607,18 @@ class MDA(MDODiscipline): if self.RESIDUALS_NORM in outputs: outputs.remove(self.RESIDUALS_NORM) + # Remove non-numeric arrays that cannot be differentiated + inputs = [ + input + for input in inputs + if self.input_grammar.is_array(input, numeric_only=True) + ] + outputs = [ + output + for output in outputs + if self.output_grammar.is_array(output, numeric_only=True) + ] + return super().check_jacobian( input_data=input_data, derr_approx=derr_approx, diff --git a/src/gemseo/mda/newton.py b/src/gemseo/mda/newton.py index 666ff32b76..ad5f8af68b 100644 --- a/src/gemseo/mda/newton.py +++ b/src/gemseo/mda/newton.py @@ -211,9 +211,12 @@ class MDANewtonRaphson(MDARoot): Compute the increment :math:`-[dR/dW]^{-1}.R` and update the MDA local_data. """ + numerical_couplings = sorted( + list(set(self.all_couplings) - set(self._non_numeric_array_variables)) + ) newton_step = self.assembly.compute_newton_step( self.local_data, - self.all_couplings, + numerical_couplings, self.relax_factor, self.__newton_linear_solver_name, matrix_type=self.matrix_type, diff --git a/tests/mda/test_mda.py b/tests/mda/test_mda.py index 7e1fcdfce5..37c3bf563c 100644 --- a/tests/mda/test_mda.py +++ b/tests/mda/test_mda.py @@ -17,6 +17,7 @@ # documentation # :author: Francois Gallard # OTHER AUTHORS - MACROSCOPIC CHANGES +# :author: Vincent Ambert from __future__ import annotations import logging @@ -39,6 +40,8 @@ from gemseo.problems.scalable.linear.disciplines_generator import ( from gemseo.problems.sellar.sellar import Sellar1 from gemseo.problems.sellar.sellar import Sellar2 from gemseo.problems.sellar.sellar import SellarSystem +from numpy import zeros +from numpy.testing import assert_almost_equal DIRNAME = os.path.dirname(__file__) @@ -170,7 +173,7 @@ def test_consistency_fail(desc, log_message, caplog): @pytest.mark.parametrize( "grammar_type", [MDODiscipline.GrammarType.JSON, MDODiscipline.GrammarType.SIMPLE] ) -def test_array_couplings(mda_class, grammar_type): +def test_array_couplings(mda_class, grammar_type, caplog): disciplines = create_disciplines_from_desc( [("A", ["x", "y1"], ["y2"]), ("B", ["x", "y2"], ("y1",))], grammar_type=grammar_type, @@ -184,8 +187,10 @@ def test_array_couplings(mda_class, grammar_type): with pytest.raises(InvalidDataError): a_disc.execute({"x": 2.0}) - with pytest.raises(TypeError, match="must be of type array"): + with caplog.at_level(logging.WARNING): mda_class(disciplines, grammar_type=grammar_type) + msg = "The coupling variable(s) ['y1'] is/are not an array of numeric values." + assert msg in caplog.text def test_convergence_warning(caplog): @@ -294,10 +299,7 @@ def test_not_numeric_couplings(): sub_prop = prop.get("items") sub_prop["type"] = "string" - with pytest.raises( - TypeError, match=r"The coupling variables \['y\_1'\] must be of type array\." - ): - MDA([sellar1, sellar2]) + MDA([sellar1, sellar2]) @pytest.mark.parametrize("mda_class", [MDAJacobi, MDAGaussSeidel, MDANewtonRaphson]) @@ -310,3 +312,132 @@ def test_get_sub_disciplines(mda_class, sellar_disciplines): """ mda = mda_class(sellar_disciplines) assert mda.get_sub_disciplines() == mda.disciplines == sellar_disciplines + + +class DiscWithNonNumericInputs1(MDODiscipline): + def __init__(self): + super().__init__() + self.input_grammar.update_from_data({"x": zeros(1)}) + self.input_grammar.update_from_data({"a": zeros(1)}) + self.input_grammar.update_from_data({"b_file": "my_b_file"}) + + self.output_grammar.update_from_data({"y": zeros(1)}) + self.output_grammar.update_from_data({"b": zeros(1)}) + self.output_grammar.update_from_data({"a_file": "my_a_file"}) + + self.default_inputs["x"] = np.array([0.5]) + self.default_inputs["a"] = np.array([0.5]) + self.default_inputs["b_file"] = "my_b_file" + + def _run(self) -> None: + x = self.local_data["x"] + a = self.local_data["a"] + y = x + b = 2 * a + self.local_data.update({"y": y, "a_file": "my_a_file", "b": b}) + + def _compute_jacobian( + self, + inputs, + outputs, + ) -> None: + self.jac = {} + self.jac["y"] = {} + self.jac["b"] = {} + self.jac["y"]["x"] = np.array([[1.0]]) + self.jac["y"]["a"] = np.array([[0.0]]) + self.jac["b"]["x"] = np.array([[0.0]]) + self.jac["b"]["a"] = np.array([[2.0]]) + + +class DiscWithNonNumericInputs2(MDODiscipline): + def __init__(self): + super().__init__() + self.input_grammar.update_from_data({"y": zeros(1)}) + self.input_grammar.update_from_data({"a_file": "my_b_file"}) + self.output_grammar.update_from_data({"x": zeros(1)}) + self.output_grammar.update_from_data({"b_file": "my_a_file"}) + self.default_inputs["y"] = np.array([0.5]) + self.default_inputs["a_file"] = "my_a_file" + + def _run(self) -> None: + y = self.local_data["y"] + x = 1.0 - 0.3 * y + self.local_data.update({"x": x, "b_file": "my_b_file"}) + + def _compute_jacobian( + self, + inputs, + outputs, + ) -> None: + self.jac = {} + self.jac["x"] = {} + self.jac["x"]["y"] = np.array([[-0.3]]) + + +class DiscWithNonNumericInputs3(MDODiscipline): + def __init__(self): + super().__init__() + self.input_grammar.update_from_data({"x": zeros(1)}) + self.input_grammar.update_from_data({"y": zeros(1)}) + self.input_grammar.update_from_data({"b": zeros(1)}) + self.input_grammar.update_from_data({"a_file": "my_a_file"}) + + self.output_grammar.update_from_data({"obj": zeros(1)}) + self.output_grammar.update_from_data({"out_file_2": "my_a_file"}) + self.default_inputs["x"] = np.array([0.5]) + self.default_inputs["y"] = np.array([0.5]) + self.default_inputs["b"] = np.array([0.5]) + self.default_inputs["a_file"] = "my_a_file" + + def _run(self) -> None: + x = self.local_data["x"] + y = self.local_data["y"] + obj = x - y + self.local_data.update({"obj": obj, "out_file_2": "my_out_file"}) + + def _compute_jacobian( + self, + inputs, + outputs, + ) -> None: + x = self.local_data["x"][0] + y = self.local_data["y"][0] + self.jac = {} + self.jac["obj"] = {} + self.jac["obj"]["x"] = np.array([[1.0]]) + self.jac["obj"]["y"] = np.array([[-1.0]]) + self.jac["obj"]["b"] = np.array([[x - y]]) + + +@pytest.mark.parametrize("mda_class", [MDAJacobi, MDAGaussSeidel, MDANewtonRaphson]) +def test_sellar_mda_with_non_numeric_couplings(mda_class, sellar_disciplines): + """Test the get_sub_disciplines method. + + Args: + mda_class: The specific MDA to be tested. + sellar_disciplines: Fixture that returns the disciplines of the Sellar problem. + """ + disciplines = [ + DiscWithNonNumericInputs1(), + DiscWithNonNumericInputs2(), + DiscWithNonNumericInputs3(), + ] + mda = mda_class(disciplines) + mda.add_differentiated_inputs(["a"]) + mda.add_differentiated_outputs(["obj"]) + + inputs = { + "a_file": "test", + "b_file": "test", + } + mda_output = mda.execute(inputs) + assert_almost_equal(mda_output["obj"], 0, decimal=5) + + assert mda.check_jacobian( + input_data=inputs, + inputs=["a"], + outputs=["obj"], + linearization_mode="adjoint", + threshold=1e-3, + ) -- GitLab From 000077ee9d03d9ade066f80681caf02a710fa54c Mon Sep 17 00:00:00 2001 From: Jean-Christophe Giret Date: Tue, 9 May 2023 15:29:45 +0200 Subject: [PATCH 2/4] docs: add changelog fragment --- changelog/fragments/731.added.rst | 1 + 1 file changed, 1 insertion(+) create mode 100644 changelog/fragments/731.added.rst diff --git a/changelog/fragments/731.added.rst b/changelog/fragments/731.added.rst new file mode 100644 index 0000000000..0951ccfb34 --- /dev/null +++ b/changelog/fragments/731.added.rst @@ -0,0 +1 @@ +:class:`.MDA`s resolution and linearization now filter non-numeric coupling variables. -- GitLab From 1f4c986d18b643ba8839cb30abe23e2d7dc271bf Mon Sep 17 00:00:00 2001 From: Jean-Christophe Giret Date: Tue, 9 May 2023 15:34:57 +0200 Subject: [PATCH 3/4] fix: minor fixes in MDA tests --- tests/mda/test_mda.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/tests/mda/test_mda.py b/tests/mda/test_mda.py index 37c3bf563c..a5ba090d2a 100644 --- a/tests/mda/test_mda.py +++ b/tests/mda/test_mda.py @@ -411,12 +411,11 @@ class DiscWithNonNumericInputs3(MDODiscipline): @pytest.mark.parametrize("mda_class", [MDAJacobi, MDAGaussSeidel, MDANewtonRaphson]) -def test_sellar_mda_with_non_numeric_couplings(mda_class, sellar_disciplines): +def test_mda_with_non_numeric_couplings(mda_class): """Test the get_sub_disciplines method. Args: mda_class: The specific MDA to be tested. - sellar_disciplines: Fixture that returns the disciplines of the Sellar problem. """ disciplines = [ DiscWithNonNumericInputs1(), -- GitLab From 156b9e81a1716058d88142d7ba0bd482d568fbf6 Mon Sep 17 00:00:00 2001 From: Jean-Christophe Giret Date: Fri, 12 May 2023 11:59:59 +0200 Subject: [PATCH 4/4] fix: preserve the _input_coupling protected attribute and add the _numeric_input_coupling protected one --- src/gemseo/mda/jacobi.py | 7 ++++--- src/gemseo/mda/mda.py | 7 ++++++- 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/src/gemseo/mda/jacobi.py b/src/gemseo/mda/jacobi.py index f231f38a0b..6c166ad086 100644 --- a/src/gemseo/mda/jacobi.py +++ b/src/gemseo/mda/jacobi.py @@ -142,8 +142,9 @@ class MDAJacobi(MDA): inputs = self.get_input_data_names() strong_cpl = self.coupling_structure.all_couplings - self._input_couplings = set(strong_cpl) & set(inputs) - set( - self._non_numeric_array_variables + self._input_couplings = sorted(set(strong_cpl) & set(inputs)) + self._numeric_input_couplings = sorted( + set(self._input_couplings) - set(self._non_numeric_array_variables) ) def execute_all_disciplines( @@ -237,7 +238,7 @@ class MDAJacobi(MDA): """ self._dx_n.append(new_couplings - current_couplings) self._g_x_n.append(new_couplings) - coupl_names = self._input_couplings + coupl_names = self._numeric_input_couplings if self.sizes is None: self.sizes = { diff --git a/src/gemseo/mda/mda.py b/src/gemseo/mda/mda.py index 013c48de1b..e29c7a80ee 100644 --- a/src/gemseo/mda/mda.py +++ b/src/gemseo/mda/mda.py @@ -177,6 +177,8 @@ class MDA(MDODiscipline): self.strong_couplings = self.coupling_structure.strong_couplings self.all_couplings = self.coupling_structure.all_couplings self._input_couplings = [] + self._numeric_input_couplings = [] + self._non_numeric_array_variables = [] self.matrix_type = JacobianAssembly.JacobianType.MATRIX self.use_lu_fact = use_lu_fact # By default don't use an approximate cache for linearization @@ -282,7 +284,10 @@ class MDA(MDODiscipline): def _compute_input_couplings(self) -> None: """Compute the strong couplings that are inputs of the MDA.""" input_couplings = set(self.strong_couplings) & set(self.get_input_data_names()) - self._input_couplings = list(input_couplings) + self._input_couplings = sorted(list(input_couplings)) + self._numeric_input_couplings = sorted( + input_couplings - set(self._non_numeric_array_variables) + ) def _current_input_couplings(self) -> ndarray: """Return the current values of the input coupling numeric variables.""" -- GitLab