Commit a5dfe96a by Stephan Rave Committed by GitHub

### Merge pull request #1277 from pymor/symbolic_expressions

Use symbolic expressions in ExpressionFunction
parents 12a8f0a6 573889d8
Pipeline #102031 passed with stages
in 41 minutes and 52 seconds
 ... ... @@ -2,6 +2,7 @@ import re MANUAL_SKIPS = ('pymor.analyticalproblems.domaindescriptions.DomainDescription.dim', 'pymor.analyticalproblems.domaindescriptions.DomainDescription.boundary_types', 'pymor.analyticalproblems.expressions.Expression.shape', 'pymor.core.base.BasicObject.name', 'pymor.core.base.BasicObject.logging_disabled', 'pymor.core.base.BasicObject.logger', ... ...
 ... ... @@ -94,32 +94,15 @@ For the definition of the source term {math}f we use an {{ ExpressionFunction }} which is given an arbitrary Python expression used to evaluate the function. In this expression, the coordinates at which the function shall be evaluated are given as the variable x. Many NumPy functions can be used directly. The entire NumPy module is available under the name np. Thus, to define {math}f we could write  (sqrt((x[0]-0.5)**2 + (x[1]-0.5)**2) <= 0.3) * 1.  However, pyMOR {{ Functions }} are required to be vectorized with respect to the coordinate x. In the case of {{ ExpressionFunction }} this means that x can be an arbitrary dimensional NumPy array of coordinates where the last array index specifies the spacial dimension. Therefore, the correct definition of {math}f is: Many NumPy functions can be used directly. Thus, to define {math}f we can write {code-cell} rhs = ExpressionFunction('(sqrt( (x[...,0]-0.5)**2 + (x[...,1]-0.5)**2) <= 0.3) * 1.', 2, ()) rhs = ExpressionFunction('(sqrt( (x[0]-0.5)**2 + (x[1]-0.5)**2) <= 0.3) * 1.', 2)  Similarly to {{ ConstantFunction }}, the second argument is the dimension of the computational domain. As the shape of the return value cannot be easily inferred from the given string expression, it has to be provided as a third argument to {{ ConstantFunction }}. For scalar functions we provide the empty tuple (), for functions returning three-dimensional vectors we would specify (3,), and for functions returning {math}2\times 2 matrices we would specify (2,2). of the computational domain. Finally, the computational domain and all data functions are collected in a {{ StationaryProblem }}: ... ... @@ -228,7 +211,7 @@ The diffusivity can be defined similarly as above: {code-cell} neumann_data = ConstantFunction(-1., 2) diffusion = ExpressionFunction('1. - (sqrt( (x[...,0]-0.5)**2 + (x[...,1]-0.5)**2) <= 0.3) * 0.999' , 2, ()) diffusion = ExpressionFunction('1. - (sqrt( (x[0]-0.5)**2 + (x[1]-0.5)**2) <= 0.3) * 0.999' , 2) problem = StationaryProblem( domain=domain, ... ... @@ -253,9 +236,8 @@ following definition: {code-cell} diffusion = ExpressionFunction( '1. - (sqrt( (np.mod(x[...,0],1./K)-0.5/K)**2 + (np.mod(x[...,1],1./K)-0.5/K)**2) <= 0.3/K) * 0.999', 2, (), values={'K': 10} '1. - (sqrt( (np.mod(x[0],1./K)-0.5/K)**2 + (np.mod(x[1],1./K)-0.5/K)**2) <= 0.3/K) * 0.999', 2, values={'K': 10} )  ... ... @@ -332,7 +314,7 @@ will be We can then make the following definition of the Neumann data: {code-cell} neumann_data = ExpressionFunction('-cos(pi*x[...,0])**2*neum[0]', 2, (), parameters= {'neum': 1}) neumann_data = ExpressionFunction('-cos(pi*x[0])**2*neum[0]', 2, parameters= {'neum': 1})  Similar to the range of the function, pyMOR cannot infer from the given ... ... @@ -346,9 +328,8 @@ We can then proceed as usual and automatically obtain a parametric {code-cell} diffusion = ExpressionFunction( '1. - (sqrt( (np.mod(x[...,0],1./K)-0.5/K)**2 + (np.mod(x[...,1],1./K)-0.5/K)**2) <= 0.3/K) * 0.999', 2, (), values={'K': 10} '1. - (sqrt( (np.mod(x[0],1./K)-0.5/K)**2 + (np.mod(x[1],1./K)-0.5/K)**2) <= 0.3/K) * 0.999', 2, values={'K': 10} ) problem = StationaryProblem( domain=domain, ... ... @@ -382,9 +363,8 @@ Next we also want to parameterize the diffusivity in the {code-cell} diffusion = ExpressionFunction( '1. - (sqrt( (np.mod(x[...,0],1./K)-0.5/K)**2 + (np.mod(x[...,1],1./K)-0.5/K)**2) <= 0.3/K) * (1 - diffu[0])', 2, (), values={'K': 10}, '1. - (sqrt( (np.mod(x[0],1./K)-0.5/K)**2 + (np.mod(x[1],1./K)-0.5/K)**2) <= 0.3/K) * (1 - diffu[0])', 2, values={'K': 10}, parameters= {'diffu': 1} )  ... ...
 ... ... @@ -114,15 +114,15 @@ problem = StationaryProblem( domain=RectDomain(), rhs=LincombFunction( [ExpressionFunction('ones(x.shape[:-1]) * 10', 2, ()), ConstantFunction(1., 2)], [ExpressionFunction('10', 2), ConstantFunction(1., 2)], [ProjectionParameterFunctional('mu'), 0.1]), diffusion=LincombFunction( [ExpressionFunction('1 - x[..., 0]', 2, ()), ExpressionFunction('x[..., 0]', 2, ())], [ExpressionFunction('1 - x[0]', 2), ExpressionFunction('x[0]', 2)], [ProjectionParameterFunctional('mu'), 1]), dirichlet_data=LincombFunction( [ExpressionFunction('2 * x[..., 0]', 2, ()), ConstantFunction(1., 2)], [ExpressionFunction('2 * x[0]', 2), ConstantFunction(1., 2)], [ProjectionParameterFunctional('mu'), 0.5]), name='2DProblem' ... ...
 ... ... @@ -131,12 +131,12 @@ import numpy as np domain = RectDomain(([-1,-1], [1,1])) indicator_domain = ExpressionFunction( '(-2/3. <= x[..., 0]) * (x[..., 0] <= -1/3.) * (-2/3. <= x[..., 1]) * (x[..., 1] <= -1/3.) * 1. \ + (-2/3. <= x[..., 0]) * (x[..., 0] <= -1/3.) * (1/3. <= x[..., 1]) * (x[..., 1] <= 2/3.) * 1.', dim_domain=2, shape_range=()) '(-2/3. <= x[0]) * (x[0] <= -1/3.) * (-2/3. <= x[1]) * (x[1] <= -1/3.) * 1. \ + (-2/3. <= x[0]) * (x[0] <= -1/3.) * (1/3. <= x[1]) * (x[1] <= 2/3.) * 1.', dim_domain=2) rest_of_domain = ConstantFunction(1, 2) - indicator_domain f = ExpressionFunction('0.5*pi*pi*cos(0.5*pi*x[..., 0])*cos(0.5*pi*x[..., 1])', dim_domain=2, shape_range=()) f = ExpressionFunction('0.5*pi*pi*cos(0.5*pi*x[0])*cos(0.5*pi*x[1])', dim_domain=2) parameters = {'diffusion': 2} thetas = [ExpressionParameterFunctional('1.1 + sin(diffusion[0])*diffusion[1]', parameters, ... ...
 ... ... @@ -35,10 +35,10 @@ def burgers_problem(v=1., circle=True, initial_data_type='sin', parameter_range= assert initial_data_type in ('sin', 'bump') if initial_data_type == 'sin': initial_data = ExpressionFunction('0.5 * (sin(2 * pi * x) + 1.)', 1, ()) initial_data = ExpressionFunction('0.5 * (sin(2 * pi * x[0]) + 1.)', 1) dirichlet_data = ConstantFunction(dim_domain=1, value=0.5) else: initial_data = ExpressionFunction('(x >= 0.5) * (x <= 1) * 1.', 1, ()) initial_data = ExpressionFunction('(x[0] >= 0.5) * (x[0] <= 1) * 1.', 1) dirichlet_data = ConstantFunction(dim_domain=1, value=0.) return InstationaryProblem( ... ... @@ -51,10 +51,10 @@ def burgers_problem(v=1., circle=True, initial_data_type='sin', parameter_range= rhs=None, nonlinear_advection=ExpressionFunction('abs(x)**exponent[0] * v', 1, (1,), {'exponent': 1}, {'v': v}), 1, {'exponent': 1}, {'v': v}), nonlinear_advection_derivative=ExpressionFunction('exponent * abs(x)**(exponent[0]-1) * sign(x) * v', 1, (1,), {'exponent': 1}, {'v': v}), 1, {'exponent': 1}, {'v': v}), ), T=0.3, ... ... @@ -94,10 +94,10 @@ def burgers_problem_2d(vx=1., vy=1., torus=True, initial_data_type='sin', parame assert initial_data_type in ('sin', 'bump') if initial_data_type == 'sin': initial_data = ExpressionFunction("0.5 * (sin(2 * pi * x[..., 0]) * sin(2 * pi * x[..., 1]) + 1.)", 2, ()) initial_data = ExpressionFunction("0.5 * (sin(2 * pi * x[0]) * sin(2 * pi * x[1]) + 1.)", 2) dirichlet_data = ConstantFunction(dim_domain=2, value=0.5) else: initial_data = ExpressionFunction("(x[..., 0] >= 0.5) * (x[..., 0] <= 1) * 1", 2, ()) initial_data = ExpressionFunction("(x[0] >= 0.5) * (x[0] <= 1) * 1", 2) dirichlet_data = ConstantFunction(dim_domain=2, value=0.) return InstationaryProblem( ... ... @@ -109,11 +109,11 @@ def burgers_problem_2d(vx=1., vy=1., torus=True, initial_data_type='sin', parame rhs=None, nonlinear_advection=ExpressionFunction("abs(x)**exponent * v", 1, (2,), {'exponent': 1}, {'v': np.array([vx, vy])}), nonlinear_advection=ExpressionFunction("abs(x[0])**exponent * v", 1, {'exponent': 1}, {'v': [vx, vy]}), nonlinear_advection_derivative=ExpressionFunction("exponent * abs(x)**(exponent-1) * sign(x) * v", 1, (2,), {'exponent': 1}, {'v': np.array([vx, vy])}), nonlinear_advection_derivative=ExpressionFunction("exponent * abs(x[0])**(exponent-1) * sign(x[0]) * v", 1, {'exponent': 1}, {'v': [vx, vy]}), ), initial_data=initial_data, ... ...
 # 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) """This module contains a basic symbolic expression library. The library is used by |ExpressionFunction| and |ExpressionParameterFunctional| by calling :func:parse_expression, which parses str expressions by replacing the names in the string with objects from this module. The result is an :class:Expression object, which can be converted to a |NumPy|-vectorized function using :meth:Expression.to_numpy. In the future, more conversion routines will be added to make the same :class:Expression usable for :mod:pymor.discretizers that use external PDE solvers. Further advantages of using this expression library: - meaningful error messages are generated at parse time of the str expression, instead of hard-to-debug errors in lambda functions at evaluation time, - expressions are automatically correctly vectorized. In particular, there is no longer a need to add ... to indexing expressions, - the shape of the resulting expressions is automatically determined. In the future, we will also provide support for symbolic differentiation of the given :class:Expressions. Every :class:Expression is built from the following atoms: - a :class:Constant, which is a fixed value of arbitrary shape, - a :class:Parameter, which is a variable of a fixed one-dimensional shape. Note that both |Parameters| and input variables are treated as a :class:Parameter in the expression. Only when calling, e.g., :meth:~Expression.to_numpy, it is determined which :class:Parameter belongs to the resulting function's |Parameters| and which :class:Parameter is treated as an input variable. More complex expressions can be built using: - :class:binary operators , - :class:negation , - :class:function calls , - :class:indexing , - :class:array construction . For binary operations of :class:Expressions  of different shape, the usual broadcasting rules apply. """ from numbers import Number from itertools import zip_longest import numpy as np from pymor.parameters.base import ParametricObject builtin_max = max def parse_expression(expression, parameters={}, values={}): if isinstance(expression, Expression): return expression locals_dict = {name: Parameter(name, dim) for name, dim in parameters.items()} return _convert_to_expression(eval(expression, dict(globals(), **values), locals_dict)) class Expression(ParametricObject): """A symbolic math expression Attributes ---------- shape The shape of the object this expression evaluates to in the sense of |NumPy|. """ shape = None def to_numpy(self, variables): """Convert expression to a |NumPy|-vectorized callable. Parameters ---------- variables List of names of :class:~Parameters  in the expression which shall be treated as input variables. """ expression = self.numpy_expr() code = compile(expression, '', 'eval') def wrapper(*args, mu={}): if not variables and args: assert len(args) == 1 mu = args[0] args = [] assert all(_broadcastable_shapes(args[0].shape[:-1], a.shape[:-1]) for a in args[1:]) if len(args) == 0: input_shape = () elif len(args) == 1: input_shape = args[0].shape[:-1] else: input_shape = (tuple(builtin_max(*s) for s in zip_longest(*(a.shape[-2::-1] for a in args), fillvalue=1)))[::-1] all_args = dict(mu) if mu else {} all_args.update({k: v for k, v in zip(variables, args)}) result = np.broadcast_to(eval(code, _numpy_functions, all_args), input_shape + self.shape) return result return wrapper def numpy_expr(self): """Called by :meth:~Expression.to_numpy.""" raise NotImplementedError def __getitem__(self, index): return Indexed(self, index) def __add__(self, other): return Sum(self, _convert_to_expression(other)) def __radd__(self, other): return Sum(_convert_to_expression(other), self) def __sub__(self, other): return Diff(self, _convert_to_expression(other)) def __rsub__(self, other): return Diff(_convert_to_expression(other), self) def __mul__(self, other): return Prod(self, _convert_to_expression(other)) def __rmul__(self, other): return Prod(_convert_to_expression(other), self) def __truediv__(self, other): return Div(self, _convert_to_expression(other)) def __rtruediv__(self, other): return Div(_convert_to_expression(other), self) def __pow__(self, other): return Pow(self, _convert_to_expression(other)) def __rpow__(self, other): return Pow(_convert_to_expression(other), self) def __mod__(self, other): return Mod(self, _convert_to_expression(other)) def __rmod__(self, other): return Mod(_convert_to_expression(other), self) def __neg__(self): return Neg(self) def __le__(self, other): return LE(self, _convert_to_expression(other)) def __ge__(self, other): return GE(self, _convert_to_expression(other)) def __lt__(self, other): return LT(self, _convert_to_expression(other)) def __gt__(self, other): return GT(self, _convert_to_expression(other)) class BaseConstant(Expression): """A constant value.""" numpy_symbol = None shape = () def numpy_expr(self): return f'array({self.numpy_symbol}, ndmin={len(self.shape)}, copy=False)' def __str__(self): return str(self.numpy_symbol) class Constant(BaseConstant): """A constant value given by a |NumPy| array.""" def __init__(self, value): self.value = value = np.array(value) self.numpy_symbol = repr(value.tolist()) self.shape = value.shape def __str__(self): return str(self.value) class Parameter(Expression): """A free parameter in an :class:Expression. Parameters represent both pyMOR |Parameters| as well as input variables. Each parameter is a vector of shape (dim,). Parameters ---------- name The name of the parameter. dim The dimension of the parameter. """ def __init__(self, name, dim): self.name, self.dim = name, dim self.numpy_symbol = name self.shape = (dim,) self.parameters_own = {name: dim} def numpy_expr(self): return str(self.numpy_symbol) def __str__(self): return str(self.numpy_symbol) class Array(Expression): """An array of scalar-valued :class:Expressions .""" def __init__(self, array): array = np.array(array) for i, v in np.ndenumerate(array): if not isinstance(v, Expression): raise ValueError(f'Array entry {v} at index {i} is not an Expression.') if v.shape not in ((), (1,)): raise ValueError(f'Array entry {v} at index {i} is not scalar valued (shape: {v.shape}).') self.array = np.vectorize(lambda x: x[0] if x.shape else x)(array) self.shape = array.shape def numpy_expr(self): entries = [v.numpy_expr() for v in self.array.flat] return f'(lambda a: array(a).T.reshape(a[0].shape + {self.shape}))(broadcast_arrays({", ".join(entries)}))' def __str__(self): expr_array = np.vectorize(str)(self.array) return str(expr_array) class BinaryOp(Expression): """Compound :class:Expression of a binary operator acting on two sub-expressions.""" numpy_symbol = None def __init__(self, first, second): if not _broadcastable_shapes(first.shape, second.shape): raise ValueError(f'Incompatible shapes of expressions "{first}" and "{second}" with shapes ' f'{first.shape} and {second.shape} for binary operator {self.numpy_symbol}') self.first, self.second = first, second self.shape = tuple(builtin_max(f, s) for f, s in zip_longest(first.shape[::-1], second.shape[::-1], fillvalue=1))[::-1] def numpy_expr(self): first_ind = second_ind = '' if len(self.first.shape) > len(self.second.shape): second_ind = (['...'] + ['newaxis'] * (len(self.first.shape) - len(self.second.shape)) + [':'] * len(self.second.shape)) second_ind = '[' + ','.join(second_ind) + ']' if len(self.first.shape) < len(self.second.shape): first_ind = (['...'] + ['newaxis'] * (len(self.second.shape) - len(self.first.shape)) + [':'] * len(self.first.shape)) first_ind = '[' + ','.join(first_ind) + ']' return f'({self.first.numpy_expr()}{first_ind} {self.numpy_symbol} {self.second.numpy_expr()}{second_ind})' def __str__(self): return f'({self.first} {self.numpy_symbol} {self.second})' class Neg(Expression): """Negated :class:Expression.""" def __init__(self, operand): self.operand = operand self.shape = operand.shape def numpy_expr(self): return f'(- {self.operand.numpy_expr()})' def __str__(self): return f'(- {self.operand})' class Indexed(Expression): """Indexed :class:Expression.""" def __init__(self, base, index): if not isinstance(index, int) and \ not (isinstance(index, tuple) and all(isinstance(i, int) for i in index)): raise ValueError(f'Indices must be ints or tuples of ints (given: {index})') if isinstance(index, int): index = (index,) if not len(index) == len(base.shape): raise ValueError(f'Wrong number of indices (given: {index} for expression "{base}" of shape {base.shape})') if not all(0 <= i < s for i, s in zip(index, base.shape)): raise ValueError(f'Invalid index (given: {index} for expression "{base}" of shape {base.shape})') self.base, self.index = base, index self.shape = base.shape[len(index):] def numpy_expr(self): index = ['...'] + [repr(i) for i in self.index] return f'{self.base.numpy_expr()}[{",".join(index)}]' def __str__(self): index = [str(i) for i in self.index] return f'{self.base}[{",".join(index)}]' class UnaryFunctionCall(Expression): """Compound :class:Expression of an unary function applied to a sub-expression. The function is applied component-wise. """ numpy_symbol = None def __init__(self, *arg): if len(arg) != 1: raise ValueError(f'{self.numpy_symbol} takes a single argument (given: {arg})') self.arg = _convert_to_expression(arg[0]) self.shape = self.arg.shape def numpy_expr(self): return f'{self.numpy_symbol}({self.arg.numpy_expr()})' def __str__(self): return f'{self.numpy_symbol}({self.arg})' class UnaryReductionCall(Expression): """Compound :class:Expression of an unary function applied to a sub-expression. The function is applied to the entire vector/matrix/tensor the sub-expression evaluates to, returning a single number. """ def __init__(self, *arg): if len(arg) != 1: raise ValueError(f'{self.numpy_symbol} takes a single argument (given {arg})') self.arg = _convert_to_expression(arg[0]) self.shape = () def numpy_expr(self): return (f'(lambda _a: {self.numpy_symbol}(_a.reshape(_a.shape[:-{len(self.arg.shape)}] + (-1,)), ' f'axis=-1))({self.arg.numpy_expr()})') def __str__(self): return f'{self.numpy_symbol}({self.arg})' def _convert_to_expression(obj): """Used internally to convert literals/list constructors to an :class:Expression.""" if isinstance(obj, Expression): return obj if isinstance(obj, Number): return Constant(obj) obj = np.array(obj) if obj.dtype == object: if obj.shape: obj = np.vectorize(_convert_to_expression)(obj) return Array(obj) else: raise ValueError(f'Cannot convert {obj} to expression') else: return Constant(obj) def _broadcastable_shapes(first, second): return all(f == s or f == 1 or s == 1 for f, s in zip(first[::-1], second[::-1])) class Sum(BinaryOp): numpy_symbol = '+' # NOQA class Diff(BinaryOp): numpy_symbol = '-' # NOQA class Prod(BinaryOp): numpy_symbol = '*' # NOQA class Div(BinaryOp): numpy_symbol = '/' # NOQA class Pow(BinaryOp): numpy_symbol = '**' # NOQA class Mod(BinaryOp): numpy_symbol = '%' # NOQA class LE(BinaryOp): numpy_symbol = '<=' # NOQA class GE(BinaryOp): numpy_symbol = '>=' # NOQA class LT(BinaryOp): numpy_symbol = '<' # NOQA class GT(BinaryOp): numpy_symbol = '>' # NOQA class sin(UnaryFunctionCall): numpy_symbol = 'sin' # NOQA class cos(UnaryFunctionCall): numpy_symbol = 'cos' # NOQA class tan(UnaryFunctionCall): numpy_symbol = 'tan' # NOQA class arcsin(UnaryFunctionCall): numpy_symbol = 'arcsin' # NOQA class arccos(UnaryFunctionCall): numpy_symbol = 'arccos' # NOQA class arctan(UnaryFunctionCall): numpy_symbol = 'arctan' # NOQA class sinh(UnaryFunctionCall): numpy_symbol = 'sinh' # NOQA class cosh(UnaryFunctionCall): numpy_symbol = 'cosh' # NOQA class tanh(UnaryFunctionCall): numpy_symbol = 'tanh' # NOQA class arcsinh(UnaryFunctionCall): numpy_symbol = 'arcsinh' # NOQA class arccosh(UnaryFunctionCall): numpy_symbol = 'arccosh' # NOQA class arctanh(UnaryFunctionCall): numpy_symbol = 'arctanh' # NOQA class exp(UnaryFunctionCall): numpy_symbol = 'exp' # NOQA