Add second order partial derivatives w.r.t. parameters

AUTHORS.md

@@ -8,6 +8,11 @@
 
 # Contributors
 
+## pyMOR 2020.1
+
+* Tim Keil, tim.keil@uni-muenster.de
+    * second order derivatives for parameters
+
 ## pyMOR 2019.2
 
 * Linus Balicki, linus.balicki@ovgu.de

src/pymor/parameters/functionals.py

@@ -49,7 +49,6 @@ def d_mu(self, component, index=()):
                 return ConstantParameterFunctional(1, name=self.name + '_d_mu')
         return ConstantParameterFunctional(0, name=self.name + '_d_mu')
 
-
 class GenericParameterFunctional(ParameterFunctionalInterface):
     """A wrapper making an arbitrary Python function a |ParameterFunctional|
 
@@ -69,9 +68,12 @@ class GenericParameterFunctional(ParameterFunctionalInterface):
     derivative_mappings
         A dict containing all partial derivatives of each component and index in the
         |ParameterType| with the signature `derivative_mappings[component][index](mu)`
+    second_derivative_mappings
+        A dict containing all second order partial derivatives of each component and index in the
+        |ParameterType| with the signature `second_derivative_mappings[component_i][index_i][component_j][index_j](mu)`
     """
 
-    def __init__(self, mapping, parameter_type, name=None, derivative_mappings=None):
+    def __init__(self, mapping, parameter_type, name=None, derivative_mappings=None, second_derivative_mappings=None):
         self.__auto_init(locals())
         self.build_parameter_type(parameter_type)
 
@@ -92,10 +94,20 @@ def d_mu(self, component, index=()):
                                   partial derivatives in self.parameter_type')
             else:
                 if component in self.derivative_mappings:
-                    return GenericParameterFunctional(self.derivative_mappings[component][index],
-                                                      self.parameter_type, name=self.name + '_d_mu')
+                    if self.second_derivative_mappings is None:
+                        return GenericParameterFunctional(self.derivative_mappings[component][index],
+                                                        self.parameter_type, name=self.name + '_d_mu')
+                    else:
+                        if component in self.second_derivative_mappings:
+                            return GenericParameterFunctional(self.derivative_mappings[component][index],
+                                                    self.parameter_type, name=self.name + '_d_mu',
+                                                    derivative_mappings=self.second_derivative_mappings[component][index])
+                        else:
+                            return GenericParameterFunctional(self.derivative_mappings[component][index],
+                                                    self.parameter_type, name=self.name + '_d_mu',
+                                                    derivative_mappings={})
                 else:
-                    raise ValueError('derivative mappings does not contain item {}'.format(component))
+                    raise ValueError('derivative expressions do not contain item {}'.format(component))
         return ConstantParameterFunctional(0, name=self.name + '_d_mu')
 
 
@@ -118,10 +130,12 @@ class ExpressionParameterFunctional(GenericParameterFunctional):
         The |ParameterType| of the |Parameters| the functional expects.
     name
         The name of the functional.
-
     derivative_expressions
         A dict containing a Python expression for the partial derivatives of each
         parameter component.
+    second_derivative_expressions
+        A dict containing a list of dicts of Python expressions for all second order partial derivatives of each
+        parameter component i and j.
     """
 
     functions = {k: getattr(np, k) for k in {'sin', 'cos', 'tan', 'arcsin', 'arccos', 'arctan', 'arctan2',
@@ -133,7 +147,7 @@ class ExpressionParameterFunctional(GenericParameterFunctional):
     functions['polar'] = lambda x: (np.linalg.norm(x, axis=-1), np.arctan2(x[..., 1], x[..., 0]) % (2*np.pi))
     functions['np'] = np
 
-    def __init__(self, expression, parameter_type, name=None, derivative_expressions=None):
+    def __init__(self, expression, parameter_type, name=None, derivative_expressions=None, second_derivative_expressions=None):
         self.expression = expression
         code = compile(expression, '<expression>', 'eval')
         functions = self.functions
@@ -153,12 +167,28 @@ def get_lambda(exp_code):
                 derivative_mappings[key] = exp_array
         else:
             derivative_mappings = None
-        super().__init__(exp_mapping, parameter_type, name, derivative_mappings)
+        if second_derivative_expressions is not None:
+            second_derivative_mappings = second_derivative_expressions.copy()
+            for (key_i,key_dicts) in second_derivative_mappings.items():
+                key_dicts_array = np.array(key_dicts, dtype=object)
+                for key_dict in np.nditer(key_dicts_array, op_flags=['readwrite'], flags= ['refs_ok']):
+                    for (key_j, exp) in key_dict[()].items():
+                        exp_array = np.array(exp, dtype=object)
+                        for exp in np.nditer(exp_array, op_flags=['readwrite'], flags= ['refs_ok']):
+                            exp_code = compile(str(exp), '<expression>', 'eval')
+                            mapping = get_lambda(exp_code)
+                            exp[...] = mapping
+                        key_dict[()][key_j] = exp_array
+                second_derivative_mappings[key_i] = key_dicts_array
+        else:
+            second_derivative_mappings = None
+        super().__init__(exp_mapping, parameter_type, name, derivative_mappings, second_derivative_mappings)
         self.__auto_init(locals())
 
     def __reduce__(self):
         return (ExpressionParameterFunctional,
-                (self.expression, self.parameter_type, getattr(self, '_name', None), self.derivative_expressions))
+                (self.expression, self.parameter_type, getattr(self, '_name', None),
+                 self.derivative_expressions, self.second_derivative_expressions))
 
 
 class ProductParameterFunctional(ParameterFunctionalInterface):
@@ -212,11 +242,9 @@ def evaluate(self, mu=None):
     def d_mu(self, component, index=()):
         raise NotImplementedError
 
-
 class ConstantParameterFunctional(ParameterFunctionalInterface):
     """|ParameterFunctional| returning a constant value for each parameter.
 
-
     Parameters
     ----------
     constant_value

src/pymortests/mu_derivatives.py

@@ -19,21 +19,39 @@ def test_ProjectionParameterFunctional():
     derivative_to_first_index = pf.d_mu('mu', 0)
     derivative_to_second_index = pf.d_mu('mu', 1)
 
+    second_derivative_first_first = pf.d_mu('mu', 0).d_mu('mu', 0)
+    second_derivative_first_second = pf.d_mu('mu', 0).d_mu('mu', 1)
+    second_derivative_second_first = pf.d_mu('mu', 1).d_mu('mu', 0)
+    second_derivative_second_second = pf.d_mu('mu', 1).d_mu('mu', 1)
+
     der_mu_1 = derivative_to_first_index.evaluate(mu)
     der_mu_2 = derivative_to_second_index.evaluate(mu)
 
+    hes_mu_1_mu_1 = second_derivative_first_first.evaluate(mu)
+    hes_mu_1_mu_2 = second_derivative_first_second.evaluate(mu)
+    hes_mu_2_mu_1 = second_derivative_second_first.evaluate(mu)
+    hes_mu_2_mu_2 = second_derivative_second_second.evaluate(mu)
+
     assert pf.evaluate(mu) == 10
     assert der_mu_1 == 1
     assert der_mu_2 == 0
+    assert hes_mu_1_mu_1 == 0
+    assert hes_mu_1_mu_2 == 0
+    assert hes_mu_2_mu_1 == 0
+    assert hes_mu_2_mu_2 == 0
 
 
 def test_ExpressionParameterFunctional():
-    dict_of_d_mus = {'mu': ['100', '2'], 'nu': 'cos(nu)'}
+    dict_of_d_mus = {'mu': ['200 * mu[0]', '2 * mu[0]'], 'nu': 'cos(nu)'}
+
+    dict_of_second_derivative = {'mu': [{'mu': ['200', '2'], 'nu': '0'},
+                              {'mu': ['2', '0'], 'nu': '0'}],
+                       'nu': {'mu': ['0', '0'], 'nu': '-sin(nu)'}}
 
-    epf = ExpressionParameterFunctional('100 * mu[0] + 2 * mu[1] + sin(nu)',
+    epf = ExpressionParameterFunctional('100 * mu[0]**2 + 2 * mu[1] * mu[0] + sin(nu)',
                                         {'mu': (2,), 'nu': ()},
-                                        'functional_with_derivative',
-                                        dict_of_d_mus)
+                                        'functional_with_derivative_and_second_derivative',
+                                        dict_of_d_mus, dict_of_second_derivative)
 
     mu = {'mu': (10,2), 'nu': 0}
 
@@ -45,10 +63,39 @@ def test_ExpressionParameterFunctional():
     der_mu_2 = derivative_to_second_mu_index.evaluate(mu)
     der_nu = derivative_to_nu_index.evaluate(mu)
 
-    assert epf.evaluate(mu) == 100*10 + 2*2 + 0
-    assert der_mu_1 == 100
-    assert der_mu_2 == 2
+    second_derivative_first_mu_first_mu = epf.d_mu('mu', 0).d_mu('mu', 0)
+    second_derivative_first_mu_second_mu = epf.d_mu('mu', 0).d_mu('mu', 1)
+    second_derivative_first_mu_nu = epf.d_mu('mu', 0).d_mu('nu')
+    second_derivative_second_mu_first_mu = epf.d_mu('mu', 1).d_mu('mu', 0)
+    second_derivative_second_mu_second_mu = epf.d_mu('mu', 1).d_mu('mu', 1)
+    second_derivative_second_mu_nu = epf.d_mu('mu', 1).d_mu('nu')
+    second_derivative_nu_first_mu = epf.d_mu('nu').d_mu('mu', 0)
+    second_derivative_nu_second_mu = epf.d_mu('nu').d_mu('mu', 1)
+    second_derivative_nu_nu = epf.d_mu('nu').d_mu('nu')
+
+    hes_mu_1_mu_1 = second_derivative_first_mu_first_mu.evaluate(mu)
+    hes_mu_1_mu_2 = second_derivative_first_mu_second_mu.evaluate(mu)
+    hes_mu_1_nu = second_derivative_first_mu_nu.evaluate(mu)
+    hes_mu_2_mu_1 = second_derivative_second_mu_first_mu.evaluate(mu)
+    hes_mu_2_mu_2 = second_derivative_second_mu_second_mu.evaluate(mu)
+    hes_mu_2_nu = second_derivative_second_mu_nu.evaluate(mu)
+    hes_nu_mu_1 = second_derivative_nu_first_mu.evaluate(mu)
+    hes_nu_mu_2 = second_derivative_nu_second_mu.evaluate(mu)
+    hes_nu_nu = second_derivative_nu_nu.evaluate(mu)
+
+    assert epf.evaluate(mu) == 100*10**2 + 2*2*10 + 0
+    assert der_mu_1 == 200 * 10
+    assert der_mu_2 == 2 * 10
     assert der_nu == 1
+    assert hes_mu_1_mu_1 == 200
+    assert hes_mu_1_mu_2 == 2
+    assert hes_mu_1_nu == 0
+    assert hes_mu_2_mu_1 == 2
+    assert hes_mu_2_mu_2 == 0
+    assert hes_mu_2_nu == 0
+    assert hes_nu_mu_1 == 0
+    assert hes_nu_mu_2 == 0
+    assert hes_nu_nu == -0
 
 def test_ExpressionParameterFunctional_for_2d_array():
     dict_of_d_mus_2d = {'mu': [['10', '20'],['1', '2']], 'nu': 'cos(nu)'}
@@ -82,13 +129,17 @@ def test_ExpressionParameterFunctional_for_2d_array():
     assert der_nu == 1
 
 def test_d_mu_of_LincombOperator():
-    dict_of_d_mus = {'mu': ['100', '2'], 'nu': 'cos(nu)'}
+    dict_of_d_mus = {'mu': ['100', '2 * mu[0]'], 'nu': 'cos(nu)'}
+
+    dict_of_second_derivative = {'mu': [{'mu': ['0', '2'], 'nu': '0'},
+                              {'mu': ['2', '0'], 'nu': '0'}],
+                       'nu': {'mu': ['0', '0'], 'nu': '-sin(nu)'}}
 
     pf = ProjectionParameterFunctional('mu', (2,), (0,))
-    epf = ExpressionParameterFunctional('100 * mu[0] + 2 * mu[1] + sin(nu)',
+    epf = ExpressionParameterFunctional('100 * mu[0] + 2 * mu[1] * mu[0] + sin(nu)',
                                         {'mu': (2,), 'nu': ()},
                                         'functional_with_derivative',
-                                        dict_of_d_mus)
+                                        dict_of_d_mus, dict_of_second_derivative)
 
     mu = {'mu': (10,2), 'nu': 0}
 
@@ -107,7 +158,37 @@ def test_d_mu_of_LincombOperator():
     eval_mu_2 = op_sensitivity_to_second_mu.evaluate_coefficients(mu)
     eval_nu = op_sensitivity_to_nu.evaluate_coefficients(mu)
 
-    assert operator.evaluate_coefficients(mu) == [1., 10, 1004.]
+    second_derivative_first_mu_first_mu = operator.d_mu('mu', 0).d_mu('mu', 0)
+    second_derivative_first_mu_second_mu = operator.d_mu('mu', 0).d_mu('mu', 1)
+    second_derivative_first_mu_nu = operator.d_mu('mu', 0).d_mu('nu')
+    second_derivative_second_mu_first_mu = operator.d_mu('mu', 1).d_mu('mu', 0)
+    second_derivative_second_mu_second_mu = operator.d_mu('mu', 1).d_mu('mu', 1)
+    second_derivative_second_mu_nu = operator.d_mu('mu', 1).d_mu('nu')
+    second_derivative_nu_first_mu = operator.d_mu('nu').d_mu('mu', 0)
+    second_derivative_nu_second_mu = operator.d_mu('nu').d_mu('mu', 1)
+    second_derivative_nu_nu = operator.d_mu('nu').d_mu('nu')
+
+    hes_mu_1_mu_1 = second_derivative_first_mu_first_mu.evaluate_coefficients(mu)
+    hes_mu_1_mu_2 = second_derivative_first_mu_second_mu.evaluate_coefficients(mu)
+    hes_mu_1_nu = second_derivative_first_mu_nu.evaluate_coefficients(mu)
+    hes_mu_2_mu_1 = second_derivative_second_mu_first_mu.evaluate_coefficients(mu)
+    hes_mu_2_mu_2 = second_derivative_second_mu_second_mu.evaluate_coefficients(mu)
+    hes_mu_2_nu = second_derivative_second_mu_nu.evaluate_coefficients(mu)
+    hes_nu_mu_1 = second_derivative_nu_first_mu.evaluate_coefficients(mu)
+    hes_nu_mu_2 = second_derivative_nu_second_mu.evaluate_coefficients(mu)
+    hes_nu_nu = second_derivative_nu_nu.evaluate_coefficients(mu)
+
+    assert operator.evaluate_coefficients(mu) == [1., 10, 1040.]
     assert eval_mu_1 == [0., 1., 100.]
-    assert eval_mu_2 == [0., 0., 2.]
+    assert eval_mu_2 == [0., 0., 2. * 10]
     assert eval_nu == [0., 0., 1.]
+
+    assert hes_mu_1_mu_1 == [0. ,0. , 0.]
+    assert hes_mu_1_mu_2 == [0. ,0. ,2]
+    assert hes_mu_1_nu == [0. ,0. ,0]
+    assert hes_mu_2_mu_1 == [0. ,0. ,2]
+    assert hes_mu_2_mu_2 == [0. ,0. ,0]
+    assert hes_mu_2_nu == [0. ,0. ,0]
+    assert hes_nu_mu_1 == [0. ,0. ,0]
+    assert hes_nu_mu_2 == [0. ,0. ,0]
+    assert hes_nu_nu == [0. ,0. ,-0]