Source code for pymor.parameters.functionals

# This file is part of the pyMOR project (http://www.pymor.org).
# Copyright 2013-2020 pyMOR developers and contributors. All rights reserved.
# License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)

from numbers import Number

import numpy as np

from pymor.core.base import abstractmethod
from pymor.parameters.base import Mu, ParametricObject, Parameters
from pymor.tools.floatcmp import float_cmp


[docs]class ParameterFunctional(ParametricObject): """Interface for |Parameter| functionals. A parameter functional is simply a function mapping |Parameters| to a number. """
[docs] @abstractmethod def evaluate(self, mu=None): """Evaluate the functional for given |parameter values| `mu`.""" pass
[docs] def d_mu(self, parameter, index=0): """Return the functionals's derivative with respect to a given parameter. Parameters ---------- parameter The |Parameter| w.r.t. which to return the derivative. index Index of the |Parameter|'s component w.r.t which to return the derivative. Returns ------- New |ParameterFunctional| representing the partial derivative. """ if parameter not in self.parameters: return ConstantParameterFunctional(0, name=f'{self.name}_d_{parameter}_{index}') else: raise NotImplementedError
[docs] def __call__(self, mu=None): return self.evaluate(mu)
def _add_sub(self, other, sign): if not isinstance(other, (ParameterFunctional, Number)): return NotImplemented if isinstance(other, Number): if other == 0: return self other = ConstantParameterFunctional(other) if self.name != 'LincombParameterFunctional' or not isinstance(self, LincombParameterFunctional): if other.name == 'LincombParameterFunctional' and isinstance(other, LincombParameterFunctional): functionals = (self,) + other.functionals coefficients = (1.,) + (other.coefficients if sign == 1. else tuple(-c for c in other.coefficients)) else: functionals, coefficients = (self, other), (1., sign) elif other.name == 'LincombParameterFunctional' and isinstance(other, LincombParameterFunctional): functionals = self.functionals + other.functionals coefficients = self.coefficients + (other.coefficients if sign == 1. else tuple(-c for c in other.coefficients)) else: functionals, coefficients = self.functionals + (other,), self.coefficients + (sign,) return LincombParameterFunctional(functionals, coefficients) def _radd_sub(self, other, sign): assert not isinstance(other, ParameterFunctional) # handled by __add__/__sub__ if not isinstance(other, Number): return NotImplemented if other == 0: return self other = ConstantParameterFunctional(other) if self.name != 'LincombParameterFunctional' or not isinstance(self, LincombParameterFunctional): functionals, coefficients = (other, self), (1., sign) else: functionals = (other,) + self.functionals coefficients = (1.,) + (self.coefficients if sign == 1. else tuple(-c for c in self.coefficients)) return LincombParameterFunctional(functionals, coefficients) def __add__(self, other): return self._add_sub(other, 1.) def __sub__(self, other): return self._add_sub(other, -1.) def __radd__(self, other): return self._radd_sub(other, 1.) def __rsub__(self, other): return self._radd_sub(other, -1.) def __mul__(self, other): if not isinstance(other, (Number, ParameterFunctional)): return NotImplemented if self.name != 'ProductParameterFunctional' or not isinstance(self, ProductParameterFunctional): if isinstance(other, ProductParameterFunctional) and other.name == 'ProductParameterFunctional': return other.with_(factors=(self,) + other.factors) else: return ProductParameterFunctional((self, other)) elif isinstance(other, ProductParameterFunctional) and other.name == 'ProductParameterFunctional': factors = self.factors + other.factors return ProductParameterFunctional(factors) else: return self.with_(factors=self.factors + (other,)) __rmul__ = __mul__ def __neg__(self): return self * (-1.)
[docs]class ProjectionParameterFunctional(ParameterFunctional): """|ParameterFunctional| returning a component value of the given parameter. For given parameter map `mu`, this functional evaluates to :: mu[parameter][index] Parameters ---------- parameter The name of the parameter to return. size The size of the parameter. index See above. name Name of the functional. """ def __init__(self, parameter, size=1, index=None, name=None): assert isinstance(size, Number) if index is None and size == 1: index = 0 assert isinstance(index, Number) assert 0 <= index < size self.__auto_init(locals()) self.parameters_own = {parameter: size}
[docs] def evaluate(self, mu=None): assert self.parameters.assert_compatible(mu) return mu[self.parameter].item(self.index)
[docs] def d_mu(self, parameter, index=0): if parameter == self.parameter: assert 0 <= index < self.size if index == self.index: return ConstantParameterFunctional(1, name=f'{self.name}_d_{parameter}_{index}') return ConstantParameterFunctional(0, name=f'{self.name}_d_{parameter}_{index}')
[docs]class GenericParameterFunctional(ParameterFunctional): """A wrapper making an arbitrary Python function a |ParameterFunctional| Note that a GenericParameterFunctional can only be :mod:`pickled <pymor.core.pickle>` if the function it is wrapping can be pickled. For this reason, it is usually preferable to use :class:`ExpressionParameterFunctional` instead of :class:`GenericParameterFunctional`. Parameters ---------- mapping The function to wrap. The function has signature `mapping(mu)`. parameters The |Parameters| the functional depends on. name The name of the functional. derivative_mappings A dict containing all partial derivatives of each |Parameter| and index with the signature `derivative_mappings[parameter][index](mu)` second_derivative_mappings A dict containing all second order partial derivatives of each |Parameter| and index with the signature `second_derivative_mappings[parameter_i][index_i][parameter_j][index_j](mu)` """ def __init__(self, mapping, parameters, name=None, derivative_mappings=None, second_derivative_mappings=None): self.__auto_init(locals()) self.parameters_own = parameters
[docs] def evaluate(self, mu=None): assert self.parameters.assert_compatible(mu) value = self.mapping(mu) # ensure that we return a number not an array if isinstance(value, np.ndarray): return value.item() else: return value
[docs] def d_mu(self, parameter, index=0): if parameter in self.parameters: assert 0 <= index < self.parameters[parameter] if self.derivative_mappings is None: raise ValueError('You must provide a dict of expressions for all \ partial derivatives in self.parameters') else: if parameter in self.derivative_mappings: if self.second_derivative_mappings is None: return GenericParameterFunctional( self.derivative_mappings[parameter][index], self.parameters, name=f'{self.name}_d_{parameter}_{index}' ) else: if parameter in self.second_derivative_mappings: return GenericParameterFunctional( self.derivative_mappings[parameter][index], self.parameters, name=f'{self.name}_d_{parameter}_{index}', derivative_mappings=self.second_derivative_mappings[parameter][index] ) else: return GenericParameterFunctional( self.derivative_mappings[parameter][index], self.parameters, name=f'{self.name}_d_{parameter}_{index}', derivative_mappings={} ) else: raise ValueError('derivative expressions do not contain item {}'.format(parameter)) return ConstantParameterFunctional(0, name=f'{self.name}_d_{parameter}_{index}')
[docs]class ExpressionParameterFunctional(GenericParameterFunctional): """Turns a Python expression given as a string into a |ParameterFunctional|. Some |NumPy| arithmetic functions like `sin`, `log`, `min` are supported. For a full list see the `functions` class attribute. .. warning:: :meth:`eval` is used to evaluate the given expression. Using this class with expression strings from untrusted sources will cause mayhem and destruction! Parameters ---------- expression A Python expression in the parameter components of the given `parameters`. parameters The |Parameters| the functional depends on. 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', 'sinh', 'cosh', 'tanh', 'arcsinh', 'arccosh', 'arctanh', 'exp', 'exp2', 'log', 'log2', 'log10', 'sqrt', 'array', 'min', 'minimum', 'max', 'maximum', 'pi', 'e', 'sum', 'prod', 'abs', 'sign', 'zeros', 'ones'}} functions['norm'] = np.linalg.norm 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, parameters, name=None, derivative_expressions=None, second_derivative_expressions=None): self.expression = expression code = compile(expression, '<expression>', 'eval') functions = self.functions def get_lambda(exp_code): return lambda mu: eval(exp_code, functions, mu) exp_mapping = get_lambda(code) if derivative_expressions is not None: derivative_mappings = derivative_expressions.copy() for (key,exp) in derivative_mappings.items(): if isinstance(exp, str): exp = [exp] 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 derivative_mappings[key] = exp_array else: derivative_mappings = None if second_derivative_expressions is not None: second_derivative_mappings = second_derivative_expressions.copy() for (key_i,key_dicts) in second_derivative_mappings.items(): if isinstance(key_dicts, dict): key_dicts = [key_dicts] 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(): if isinstance(exp, str): exp = [exp] 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, parameters, name, derivative_mappings, second_derivative_mappings) self.__auto_init(locals())
[docs] def __reduce__(self): return (ExpressionParameterFunctional, (self.expression, self.parameters, getattr(self, '_name', None), self.derivative_expressions, self.second_derivative_expressions))
[docs]class ProductParameterFunctional(ParameterFunctional): """Forms the product of a list of |ParameterFunctionals| or numbers. Parameters ---------- factors A list of |ParameterFunctionals| or numbers. name Name of the functional. """ def __init__(self, factors, name=None): assert len(factors) > 0 assert all(isinstance(f, (ParameterFunctional, Number)) for f in factors) factors = tuple(factors) self.__auto_init(locals())
[docs] def evaluate(self, mu=None): assert self.parameters.assert_compatible(mu) return np.array([f.evaluate(mu) if hasattr(f, 'evaluate') else f for f in self.factors]).prod()
[docs] def d_mu(self, parameter, index=0): summands = [] for i, f in enumerate(self.factors): if hasattr(f, 'evaluate'): f_d_mu = f.d_mu(parameter, index) if isinstance(f_d_mu, ConstantParameterFunctional) and f_d_mu() == 0: continue else: continue summands.append( ProductParameterFunctional([f_d_mu if j == i else g for j, g in enumerate(self.factors)]) ) if not summands: return ConstantParameterFunctional(0, name=f'{self.name}_d_{parameter}_{index}') else: return LincombParameterFunctional(summands, [1] * len(summands), name=f'{self.name}_d_{parameter}_{index}')
[docs]class ConjugateParameterFunctional(ParameterFunctional): """Conjugate of a given |ParameterFunctional| Evaluates a given |ParameterFunctional| and returns the complex conjugate of the value. Parameters ---------- functional The |ParameterFunctional| of which the complex conjuate is taken. name Name of the functional. """ def __init__(self, functional, name=None): self.functional = functional self.name = name or f'{functional.name}_conj'
[docs] def evaluate(self, mu=None): assert self.parameters.assert_compatible(mu) return np.conj(self.functional.evaluate(mu))
[docs]class ConstantParameterFunctional(ParameterFunctional): """|ParameterFunctional| returning a constant value for each parameter. Parameters ---------- constant_value value of the functional name Name of the functional. """ def __init__(self, constant_value, name=None): self.constant_value = constant_value self.__auto_init(locals())
[docs] def evaluate(self, mu=None): return self.constant_value
[docs] def d_mu(self, parameter, index=0): return self.with_(constant_value=0, name=f'{self.name}_d_{parameter}_{index}')
[docs]class LincombParameterFunctional(ParameterFunctional): """A |ParameterFunctional| representing a linear combination of |ParameterFunctionals|. The coefficients must be provided as scalars. Parameters ---------- functionals List of |ParameterFunctionals| whose linear combination is formed. coefficients A list of scalar coefficients. name Name of the functional. Attributes ---------- functionals coefficients """ def __init__(self, functionals, coefficients, name=None): assert len(functionals) > 0 assert len(functionals) == len(coefficients) assert all(isinstance(f, ParameterFunctional) for f in functionals) assert all(isinstance(c, Number) for c in coefficients) functionals = tuple(functionals) coefficients = tuple(coefficients) self.__auto_init(locals())
[docs] def evaluate(self, mu=None): assert self.parameters.assert_compatible(mu) return sum(c * f(mu) for c, f in zip(self.coefficients, self.functionals))
[docs] def d_mu(self, parameter, index=0): functionals_d_mu = [f.d_mu(parameter, index) for f in self.functionals] return self.with_(functionals=functionals_d_mu, name=f'{self.name}_d_{parameter}_{index}')
[docs]class MinThetaParameterFunctional(ParameterFunctional): """|ParameterFunctional| implementing the min-theta approach from [Haa17]_ (Proposition 2.35). Let V denote a Hilbert space and let a: V x V -> K denote a parametric coercive bilinear form with affine decomposition :: a(u, v, mu) = sum_{q = 1}^Q theta_q(mu) a_q(u, v), for Q positive coefficient |ParameterFunctional| theta_1, ..., theta_Q and positive semi-definite component bilinear forms a_1, ..., a_Q: V x V -> K. Let mu_bar be a parameter with respect to which the coercivity constant of a(., ., mu_bar) is known, i.e. we known alpha_mu_bar > 0, s.t. :: alpha_mu_bar |u|_V^2 <= a(u, u, mu=mu_bar). The min-theta approach from [Haa17]_ (Proposition 2.35) allows to obtain a computable bound for the coercivity constant of a(., ., mu) for arbitrary parameters mu, since :: a(u, u, mu=mu) >= min_{q = 1}^Q theta_q(mu)/theta_q(mu_bar) a(u, u, mu=mu_bar). Given a list of the thetas, the |parameter values| mu_bar and the constant alpha_mu_bar, this functional thus evaluates to :: alpha_mu_bar * min_{q = 1}^Q theta_q(mu)/theta_q(mu_bar) Parameters ---------- thetas List or tuple of |ParameterFunctional| mu_bar Parameter associated with alpha_mu_bar. alpha_mu_bar Known coercivity constant. name Name of the functional. """ def __init__(self, thetas, mu_bar, alpha_mu_bar=1., name=None): assert isinstance(thetas, (list, tuple)) assert len(thetas) > 0 assert all([isinstance(theta, (Number, ParameterFunctional)) for theta in thetas]) thetas = tuple(ConstantParameterFunctional(theta) if not isinstance(theta, ParameterFunctional) else theta for theta in thetas) if not isinstance(mu_bar, Mu): mu_bar = Parameters.of(thetas).parse(mu_bar) assert Parameters.of(thetas).assert_compatible(mu_bar) thetas_mu_bar = np.array([theta(mu_bar) for theta in thetas]) assert np.all(thetas_mu_bar > 0) assert isinstance(alpha_mu_bar, Number) assert alpha_mu_bar > 0 self.__auto_init(locals()) self.thetas_mu_bar = thetas_mu_bar
[docs] def evaluate(self, mu=None): assert self.parameters.assert_compatible(mu) thetas_mu = np.array([theta(mu) for theta in self.thetas]) assert np.all(thetas_mu > 0) return self.alpha_mu_bar * np.min(thetas_mu / self.thetas_mu_bar)
[docs]class BaseMaxThetaParameterFunctional(ParameterFunctional): """|ParameterFunctional| implementing a generalization of the max-theta approach from [Haa17]_ (Exercise 5.12). Let V denote a Hilbert space and let a: V x V -> K denote a continuous bilinear form or l: V -> K a continuous linear functional, either with affine decomposition :: a(u, v, mu) = sum_{q = 1}^Q theta_q(mu) a_q(u, v) or l(v, mu) = sum_{q = 1}^Q theta_q(mu) l_q(v) for Q coefficient |ParameterFunctional| theta_1, ..., theta_Q and continuous bilinear forms a_1, ..., a_Q: V x V -> K or continuous linear functionals l_q: V -> K. Let mu_bar be a parameter with respect to which the continuity constant of a(., ., mu_bar) or l(., mu_bar) is known, i.e. we known gamma_mu_bar > 0, s.t. :: a(u, v, mu_bar) <= gamma_mu_bar |u|_V |v|_V or l(v, mu_bar) <= gamma_mu_bar |v|_V. The max-theta approach (in its generalized form) from [Haa17]_ (Exercise 5.12) allows to obtain a computable bound for the continuity constant of another bilinear form a_prime(., ., mu) or linear form l_prime(., mu) with the same affine decomposition but different theta_prime_q for arbitrary parameters mu, since :: a_prime(u, v, mu=mu) <= |max_{q = 1}^Q theta_prime_q(mu)/theta_q(mu_bar)| |a(u, v, mu=mu_bar)| or :: l_prime(v, mu=mu) <= |max_{q = 1}^Q theta_prime_q(mu)/theta_q(mu_bar)| |l(v, mu=mu_bar)|, if all theta_q(mu_bar) != 0. Given a list of the thetas, the |parameter values| mu_bar and the constant gamma_mu_bar, this functional thus evaluates to :: |gamma_mu_bar * max_{q = 1}^Q theta_prime_q(mu)/theta_q(mu_bar)| Note that we also get an upper bound if theta_prime_q(mu) == 0 for any q. However, if theta_prime_q(mu) == 0 for at least one q, we need to use the absolute value in the denominator, i.e. :: |gamma_mu_bar * max_{q = 1}^Q theta_prime_q(mu)/|theta_q(mu_bar)| | Parameters ---------- thetas List or tuple of |ParameterFunctional| of the base bilinear form a which is used for estimation. thetas_prime List or tuple of |ParameterFunctional| of the bilinear form a_prime for the numerator of the MaxThetaParameterFunctional. mu_bar Parameter associated with gamma_mu_bar. gamma_mu_bar Known continuity constant of the base bilinear form a. name Name of the functional. """ def __init__(self, thetas_prime, thetas, mu_bar, gamma_mu_bar=1., name=None): assert isinstance(thetas_prime, (list, tuple)) assert isinstance(thetas, (list, tuple)) assert len(thetas) > 0 assert len(thetas) == len(thetas_prime) assert all([isinstance(theta, (Number, ParameterFunctional)) for theta in thetas]) assert all([isinstance(theta, (Number, ParameterFunctional)) for theta in thetas_prime]) thetas = tuple(ConstantParameterFunctional(f) if not isinstance(f, ParameterFunctional) else f for f in thetas) thetas_prime = tuple(ConstantParameterFunctional(f) if not isinstance(f, ParameterFunctional) else f for f in thetas_prime) if not isinstance(mu_bar, Mu): mu_bar = Parameters.of(thetas).parse(mu_bar) assert Parameters.of(thetas).assert_compatible(mu_bar) thetas_mu_bar = np.array([theta(mu_bar) for theta in thetas]) assert not np.any(float_cmp(thetas_mu_bar, 0)) assert isinstance(gamma_mu_bar, Number) assert gamma_mu_bar > 0 self.__auto_init(locals()) self.thetas_mu_bar = thetas_mu_bar self.theta_mu_bar_has_negative = True if np.any(thetas_mu_bar < 0) else False if self.theta_mu_bar_has_negative: # If 0 is in theta_prime(mu), we need to use the absolute value to ensure # that the bound is still valid (and not zero) self.abs_thetas_mu_bar = np.array([np.abs(theta(mu_bar)) for theta in thetas])
[docs] def evaluate(self, mu=None): assert self.parameters.assert_compatible(mu) thetas_prime_mu = np.array([theta(mu) for theta in self.thetas_prime]) if np.all(np.logical_or(thetas_prime_mu < 0, thetas_prime_mu > 0)) or not self.theta_mu_bar_has_negative: return self.gamma_mu_bar * np.abs(np.max(thetas_prime_mu / self.thetas_mu_bar)) else: # special case return self.gamma_mu_bar * np.abs(np.max(thetas_prime_mu / self.abs_thetas_mu_bar))
[docs] def d_mu(self, component, index=()): raise NotImplementedError
[docs]class MaxThetaParameterFunctional(BaseMaxThetaParameterFunctional): """|ParameterFunctional| implementing the max-theta approach from [Haa17]_ (Exercise 5.12). This is a specialized version of BaseMaxThetaParameterFunctional which allows to obtain a computable bound for the continuity constant of the actual a(., ., mu) or l(., mu) for arbitrary parameters mu, since :: a(u, v, mu=mu) <= |max_{q = 1}^Q theta_q(mu)/theta_q(mu_bar)| |a(u, v, mu=mu_bar)| or :: l(v, mu=mu) <= |max_{q = 1}^Q theta_q(mu)/theta_q(mu_bar)| |l(v, mu=mu_bar)|, if all theta_q(mu_bar) != 0. Given a list of the thetas, the |parameter values| mu_bar and the constant gamma_mu_bar, this functional thus evaluates to :: |gamma_mu_bar * max{q = 1}^Q theta_q(mu)/theta_q(mu_bar)| Parameters ---------- thetas List or tuple of |ParameterFunctional| mu_bar Parameter associated with gamma_mu_bar. gamma_mu_bar Known continuity constant. name Name of the functional. """ def __init__(self, thetas, mu_bar, gamma_mu_bar=1., name=None): super().__init__(thetas, thetas, mu_bar, gamma_mu_bar, name)