Source code for nutils.expression

# Copyright (c) 2014 Evalf
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.

This module defines the function :func:`parse`, which parses a tensor

import re, collections, functools

# Convenience function to create a constant in ExpressionAST (details in
# docstring of `parse` below).
_ = lambda arg: (None, arg)

def _sp(count, singular, plural):
  '''format ``count``+ ``singular` or ``plural`` depending on ``count``'''
  return '{} {}'.format(count, singular if count == 1 else plural)

[docs]class ExpressionSyntaxError(Exception): pass
[docs]class AmbiguousAlignmentError(Exception): pass
class _IntermediateError(Exception): '''Intermediate exception, to be catched and converted into an ``ExpressionSyntaxError``.''' def __init__(self, msg, at=None, count=None): self.msg = msg = at self.count = count super().__init__(msg) _Token = collections.namedtuple('_Token', ['type', 'data', 'pos']) _Token.__doc__ = 'An indivisible part of an expression string.' _Token.type.__doc__ = 'The type of this token.' = 'Substring of the expression string that belongs to this token.' _Token.pos.__doc__ = 'The start position of the token in the expression string.' _Length = collections.namedtuple('_Length', ['pos']) _Length.__doc__ = 'Yet unknown length, introduced at ``pos`` in the expression string.' _Length.pos.__doc__ = 'The position where this :class:`_Length` is introduced.' class _Array: '''ExpressionAST with shape, indices. The :class:`_Array` class combines an ExpressionAST with shape and indices and maintains a list of summed indices in the expression string resulting in this :class:`_Array`. Attributes ---------- ast : :class:`tuple` The ExpressionAST (see :func:`parse`). indices : :class:`str` The indices of the array represented by the :attr:`ast`. shape : :class:`tuple` of :class:`int`\\s or :class:`_Length`\\s The shape of the array represented by the :attr:`ast`. summed : :class:`frozenset` of indices (:class:`str`) A set of indices that are summed in the expression string resulting in this :class:`_Array`. The indices are not allowed in expressions involving this :class:`_Array`. For example, index `i` in expression string ``'a_ij b_i'`` is summed and cannot be used in an expression like ``('a_ij b_i) c_i``. linked_lengths : :class:`frozenset` of :class:`frozensets` of :class:`_Length`\\s and :class:`int`\\s A set of sets of :class:`_Length`\\s and :class:`int`\\s. A :class:`_Length` is introduced if an axis of an :class:`_Array` has an unknown length, e.g. a dirac has two axes of equal, but unknown length. All class:`_Length`\\s in a set have the same length. If a set contains an :class:`int` the class:`_Length`\\s are resolved. ndim : :class:`int` The number of dimensions of this :class:`_Array`. Args ---- ast : :class:`tuple` See :attr:`_Array.ast`. indices : :class:`str` See :attr:`_Array.indices`. shape : :class:`tuple` of :class:`int`\\s or :class:`_Length`\\s See :attr:`_Array.shape`. summed : :class:`frozenset` of indices (:class:`str`) See :attr:`_Array.summed`. linked_lengths : :class:`frozenset` of :class:`frozensets` of :class:`_Length`\\s and :class:`int`\\s See :attr:`_Array.linked_lengths`. ''' @classmethod def wrap(cls, ast, indices, shape): '''Create an :class:`_Array` by wrapping ``ast``. The ``ast`` should be a constant or variable. Duplicate indices are summed and numeric indices are replaced by a getitem. ''' if len(indices) != len(shape): raise _IntermediateError('Expected {}, got {}.'.format(_sp(len(shape), 'index', 'indices'), len(indices))) return cls._apply_indices(ast, 0, indices, shape, frozenset(), {}) @classmethod def _apply_indices(cls, ast, offset, indices, shape, summed, linked_lengths): '''Wrap ``ast`` in an :class:`_Array`, thereby summing indices occuring twice and applying numeric indices. When wrapping a variable or gradient the indices of may appear twice, indicating summation, or numeric, indicating a getitem. This method wraps ``ast`` and applies summation and getitem if needed. Args ---- ast : :class:`tuple` See :attr:`_Array.ast`. offset : :class:`int` Start at index ``offset`` when looking for indices occuring twice (in the entire list of ``indices``, not only those in ``indices[offset:]``) or numeric indices. The list ``indices[offset:]`` is assumed to be already processed. indices : :class:`str` See :attr:`_Array.indices`. shape : :class:`tuple` of :class:`int`\\s or :class:`_Length`\\s See :attr:`_Array.shape`. ``indices``. summed : :class:`frozenset` of indices (:class:`str`) See :attr:`_Array.summed`. linked_lengths : :class:`frozenset` of :class:`frozensets` of :class:`_Length`\\s and :class:`int`\\s See :attr:`_Array.linked_lengths`. Returns ------- wrapped_ast : :class:`_Array foo : bar` ''' summed = set(summed) linked_lengths = set(linked_lengths) i = offset dims = tuple(range(len(indices))) while i < len(indices): index = indices[i] j = indices.index(index) if '0' <= index <= '9': index = int(index) if isinstance(shape[i], int) and index >= shape[i]: raise _IntermediateError('Index of dimension {} with length {} out of range.'.format(dims[i], shape[i])) ast = 'getitem', ast, _(i), _(index) indices = indices[:i] + indices[i+1:] shape = shape[:i] + shape[i+1:] dims = dims[:i] + dims[i+1:] elif index in summed: raise _IntermediateError('Index {!r} occurs more than twice.'.format(index)) elif j < i: linked_lengths = set(cls._update_lengths(linked_lengths, index, shape[j], shape[i])) ast = 'trace', ast, _(j), _(i) indices = indices[:j] + indices[j+1:i] + indices[i+1:] shape = shape[:j] + shape[j+1:i] + shape[i+1:] dims = dims[:j] + dims[j+1:i] + dims[i+1:] summed.add(index) i -= 1 else: if isinstance(shape[i], _Length) and not any(shape[i] in g for g in linked_lengths): linked_lengths.add(frozenset([shape[i]])) i += 1 return cls(ast, indices, shape, summed, linked_lengths) @classmethod def stack(cls, arrays, index): '''Stack ``arrays`` along axis ``index``. The arrays are stacked in given order. All arrays should have matching shapes, except for the axis labeled ``index``. If an array does not have the supplied ``index``, the array is expanded with an axis of length one before stacking. For example, stacking a scalar and an array with shape ``{i: 2}`` along ``i`` gives an array with shape ``{i: 3}``. Args ---- arrays : a :class:`` of :class:`_Array` objects The arrays to stack. index : :class:`str` The index along which to stack the ``arrays``. Returns ------- array : :class:`_Array` The stacked array. ''' # TODO: assert is_valid_lhs_indices(index) if len(arrays) == 0: raise _IntermediateError('Cannot stack 0 arrays.') if len(set(frozenset(array.indices) - {index} for array in arrays)) != 1: raise _IntermediateError( 'Cannot stack arrays with unmatched indices (excluding the stack index {!r}): {}.' .format(index, ', '.join(array.indices for array in arrays))) indices = index + ''.join(i for i in arrays[0].indices if i != index) arrays = [(array.append_axis(index, 1) if index not in array.indices else array).transpose(indices) for array in arrays] if len(arrays) == 1: return arrays[0] helper = arrays[0].replace(indices=arrays[0].indices[1:], shape=arrays[0].shape[1:]) for other in arrays[1:]: other = other.replace(indices=other.indices[1:], shape=other.shape[1:]) shape, linked_lengths = helper._join_shapes(other) helper = helper.replace(shape=shape, linked_lengths=linked_lengths, summed=helper.summed | other.summed) # Apply `helper.linked_lengths` to all `arrays`. If the lengths at # `index` is not known at this point, we won't be able to resolve this # ever, so raise an exception here. length = 0 for array in arrays: shape = array._simplify_shape(helper.linked_lengths) if isinstance(shape[0], _Length): raise _IntermediateError('Cannot determine the length of the stack axis, because the length at {} is unknown.'.format(shape[0].pos), at=shape[0].pos) length += shape[0] ast = ('concatenate',) + tuple(array.ast for array in arrays) return helper.replace(ast=ast, indices=indices, shape=(length,)+helper.shape) @staticmethod def align(*arrays): '''Align ``arrays`` to the first array. Args ---- arrays : :class:`_Array` The arrays to align. Returns ------- aligned_arrays : :class:`tuple` of :class:`_Array` objects The aligned arrays. ''' assert len(arrays) > 0 if len(set(frozenset(array.indices) for array in arrays)) != 1: raise _IntermediateError( 'Cannot align arrays with unmatched indices: {}.' .format(', '.join(array.indices for array in arrays))) arrays = [array.transpose(arrays[0].indices) for array in arrays] helper = arrays[0] for other in arrays[1:]: shape, linked_lengths = helper._join_shapes(other) helper = helper.replace(shape=shape, linked_lengths=linked_lengths, summed=helper.summed | other.summed) return tuple(array.replace(shape=helper.shape, linked_lengths=helper.linked_lengths, summed=helper.summed) for array in arrays) def __init__(self, ast, indices, shape, summed, linked_lengths): assert isinstance(indices, str) self.ast = tuple(ast) self.indices = indices self.shape = tuple(shape) self.summed = frozenset(summed) self.linked_lengths = frozenset(linked_lengths) self.ndim = len(self.indices) def _join_shapes(self, other): '''Verify ``self + other`` is valid and return the resulting shape and linked lengths. Args ---- other : :class:`_Array` Should have the same (order of) indices as this array. Returns ------- shape : :class:`tuple` The simplified shape of ``self + other``. linked_lengths : :class:`frozenset` of :class:`frozensets` of :class:`_Length`\\s and :class:`int`\\s See :attr:`_Array.linked_lengths`. Updated with links resulting from applying ``self + other``. ''' assert self.indices == other.indices, 'unaligned' groups = set(self.linked_lengths | other.linked_lengths) for index, a, b in zip(self.indices, self.shape, other.shape): if a == b: continue if not isinstance(a, _Length) and not isinstance(b, _Length): raise _IntermediateError('Shapes at index {!r} differ: {}, {}.'.format(index, a, b)) groups.add(frozenset({a, b})) linked_lengths = self._join_lengths(other, groups) return self._simplify_shape(linked_lengths), linked_lengths def _simplify_shape(self, linked_lengths): '''Return simplified shape by replacing :class:`_Length`\\s with :class:`int`\\s according to the ``linked_lengths``.''' shape = [] cache = {k: v for v in linked_lengths for k in v} for length in self.shape: if isinstance(length, _Length): for l in cache[length]: if not isinstance(l, _Length): length = l break shape.append(length) return shape def _join_lengths(*args): '''Return updated linked lengths resulting from ``self + other``.''' groups = set() for arg in args: groups |= arg.linked_lengths if isinstance(arg, _Array) else arg cache = {} for g in groups: # g = frozenset(itertools.chain.from_iterable(map(linked_lenghts.get, g))) new_g = set() for k in g: new_g |= cache.get(k, frozenset([k])) new_g = frozenset(new_g) cache.update((k, new_g) for k in new_g) linked_lengths = frozenset(cache.values()) # Verify. for g in linked_lengths: known = tuple(sorted(set(k for k in g if not isinstance(k, _Length)))) if len(known) > 1: raise _IntermediateError('Axes have different lengths: {}.'.format(', '.join(map(str, known)))) return linked_lengths @staticmethod def _update_lengths(linked_lengths, index, a, b): '''Add link ``a``, ``b`` to ``linked_lengths``.''' cache = {l: g for g in linked_lengths for l in g} if a != b: if not isinstance(a, _Length) and not isinstance(b, _Length): raise _IntermediateError('Shapes at index {!r} differ: {}, {}.'.format(index, a, b)) g = cache.get(a, frozenset([a])) | cache.get(b, frozenset([b])) cache.update((k, g) for k in g) # Verify. known = tuple(sorted(set(k for k in g if not isinstance(k, _Length)))) if len(known) > 1: raise _IntermediateError('Shapes at index {!r} differ: {}.'.format(index, ', '.join(map(str, known)))) elif isinstance(a, _Length): cache.setdefault(a, frozenset([a])) return frozenset(cache.values()) def __neg__(self): '''Return -self.''' return self.replace(ast=('neg', self.ast)) def _add_sub(self, other, op, name): '''Return op(self, other).''' if frozenset(self.indices) != frozenset(other.indices): raise _IntermediateError('Cannot {} arrays with unmatched indices: {!r}, {!r}.'.format(name, self.indices, other.indices)) other = other.transpose(self.indices) shape, linked_lengths = self._join_shapes(other) return _Array((op, self.ast, other.ast), self.indices, shape, self.summed, linked_lengths) def __add__(self, other): '''Return self+other.''' return self._add_sub(other, 'add', 'add') def __sub__(self, other): '''Return self-other.''' return self._add_sub(other, 'sub', 'subtract') def __mul__(self, other): '''Return self*other.''' for a, b in ((self, other), (other, self)): for index in sorted(frozenset(a.indices) | a.summed): if index in b.summed: raise _IntermediateError('Index {!r} occurs more than twice.'.format(index)) common = [] for index, length in zip(self.indices, self.shape): if index in other.indices: common.append(index) else: other = other.append_axis(index, length) for index, length in zip(other.indices, other.shape): if index not in self.indices: self = self.append_axis(index, length) indices = self.indices other = other.transpose(indices) shape, linked_lengths = self._join_shapes(other) ast = 'mul', self.ast, other.ast for index in reversed(common): i = self.indices.index(index) ast = 'sum', ast, _(i) indices = indices[:i] + indices[i+1:] shape = shape[:i] + shape[i+1:] return _Array(ast, indices, shape, self.summed | other.summed | frozenset(common), linked_lengths) def __truediv__(self, other): '''Return self/value.''' if other.ndim > 0: raise _IntermediateError('A denominator must have dimension 0.') for index in sorted((self.summed | set(self.indices)) & other.summed): raise _IntermediateError('Index {!r} occurs more than twice.'.format(index)) return _Array(('truediv', self.ast, other.ast), self.indices, self.shape, self.summed | other.summed, self._join_lengths(other)) def __pow__(self, other): '''Return self**value.''' if other.ndim > 0: raise _IntermediateError('An exponent must have dimension 0.') for index in sorted((self.summed | set(self.indices)) & other.summed): raise _IntermediateError('Index {!r} occurs more than twice.'.format(index)) return _Array(('pow', self.ast, other.ast), self.indices, self.shape, self.summed | other.summed, self._join_lengths(other)) def grad(self, index, geom, type): '''Return the gradient w.r.t. ``geom``.''' assert geom.ndim == 1 assert not isinstance(geom.shape[0], _Length) assert type in ('grad','surfgrad') ast = type, self.ast, _(geom) return _Array._apply_indices(ast, self.ndim, self.indices+index, self.shape+geom.shape, self.summed, self.linked_lengths) def derivative(self, arg, indices): 'Return the derivative to ``arg``.' return _Array._apply_indices(('derivative', self.ast, arg.ast), self.ndim, self.indices+indices, self.shape+arg.shape, self.summed, self.linked_lengths) def append_axis(self, index, length): '''Return an :class:`_Array` with one additional axis.''' if index in self.indices or index in self.summed: raise _IntermediateError('Duplicate index: {!r}.'.format(index)) linked_lengths = self.linked_lengths if isinstance(length, _Length): for group in linked_lengths: if length in group: break else: linked_lengths |= frozenset({frozenset({length})}) return _Array(('append_axis', self.ast, _(length)), self.indices+index, self.shape+(length,), self.summed, linked_lengths) def transpose(self, indices): '''Return an :class:`_Array` transposed according to ``indices``.''' if len(indices) != len(set(indices)): raise _IntermediateError('Cannot transpose from {!r} to {!r}: duplicate indices.'.format(self.indices, indices)) elif set(self.indices) != set(indices): raise _IntermediateError('Cannot transpose from {!r} to {!r}: indices differ.'.format(self.indices, indices)) if self.indices == indices: return self else: transpose = tuple(map(self.indices.index, indices)) shape = tuple(map(self.shape.__getitem__, transpose)) return _Array(('transpose', self.ast, _(transpose)), indices, shape, self.summed, self.linked_lengths) def replace(self, **updates): '''Return a copy of this :class:`_Array` with attributes replaced by ``updates``.''' kwargs = dict(ast=self.ast, indices=self.indices, shape=self.shape, summed=self.summed, linked_lengths=self.linked_lengths) kwargs.update(updates) return _Array(**kwargs) class _ExpressionParser: '''Expression parser Args ---- expression : :class:`str` See argument ``expression`` of :func:`parse`. variables : :class:`dict` of :class:`str` and :class:`nutils.function.Array` pairs See argument ``variables`` of :func:`parse`. functions : :class:`dict` of :class:`str` and :class:`int` pairs See argument ``functions`` of :func:`parse`. arg_shapes : :class:`dict` of :class:`str` and :class:`tuple` or :class:`int`\\s pairs See argument ``arg_shapes`` of :func:`parse`. default_geometry_name : class:`str` See argument ``default_geometry_name`` of :func:`parse`. ''' eye_symbols = '$', 'δ' normal_symbols = 'n', def __init__(self, expression, variables, functions, arg_shapes, default_geometry_name): self.expression = expression self.variables = variables self.functions = functions self.arg_shapes = dict(arg_shapes) self.default_geometry_name = default_geometry_name def highlight(f): 'wrap ``f`` in a function that converts ``_IntermediateError`` objects' def wrapper(self, *args, **kwargs): if hasattr(self, '_tokens'): pos = self._next.pos else: pos = 0 try: return f(self, *args, **kwargs) except _IntermediateError as e: if is None: at = pos count = self._next.pos - pos if self._next.pos > pos else len( else: at = count = 1 if e.count is None else e.count raise ExpressionSyntaxError(e.msg + '\n' + self.expression + '\n' + ' '*at + '^'*count) from e return wrapper def _consume(self): 'advance to next token' self._index += 1 if self._index >= len(self._tokens): raise _IntermediateError('Unexpected end of expression.', at=len(self.expression)) return self._current def _consume_if_whitespace(self): 'advance to next token if it is a whitespace' if self._next.type == 'whitespace': self._consume() @highlight def _consume_assert_whitespace(self): 'assert the next token is whitespace, skip it, and advance to next token' if self._consume().type != 'whitespace': raise _IntermediateError('Missing whitespace.', at=self._current.pos) @highlight def _consume_assert_equal(self, value, msg=None): 'assert the next token is equal to ``value``' token = self._consume() if token.type != value: if msg is None: msg = 'Expected {!r}.'.format(value) raise _IntermediateError(msg, at=token.pos) return token @property def _current(self): 'the current token' return self._tokens[self._index] @property def _next(self): 'the next token' return self._tokens[min(len(self._tokens)-1, self._index+1)] @property def _next_non_whitespace(self): 'the next non-whitespace token' return self._tokens[self._index+2] if self._next.type == 'whitespace' else self._next def _get_variable(self, name): 'get variable by ``name`` or raise an error' value = self.variables.get(name, None) if value is None: raise _IntermediateError('Unknown variable: {!r}.'.format(name)) return value def _get_geometry(self, name): 'get geometry by ``name`` or raise an error' geom = self._get_variable(name) if geom.ndim != 1: raise _IntermediateError('Invalid geometry: expected 1 dimension, but {!r} has {}.'.format(name, geom.ndim)) return geom def _get_arg(self, name, indices, indices_start): 'get arg by ``name`` or raise an error' if name in self.arg_shapes: shape = self.arg_shapes[name] if len(shape) != len(indices): raise _IntermediateError('Argument {!r} previously defined with {} instead of {}.'.format(name, _sp(len(shape), 'axis', 'axes'), len(indices))) else: shape = tuple(_Length(indices_start+i) for i, j in enumerate(indices)) self.arg_shapes[name] = shape return _Array.wrap(('arg', _(name)) + tuple(map(_, shape)), indices, shape) @highlight def parse_lhs_arg(self, seen_lhs): 'parse lhs arg, e.g. the "x_ij" in "x_kk(x_ij=a_ij)"' token = self._consume() if token.type != 'variable': raise _IntermediateError("Expected an argument, e.g. 'argname'.") if'?'): raise _IntermediateError("The argument name at the left hand side of a substitution must not be prefixed by a '?'.") name = if name in seen_lhs: raise _IntermediateError("Argument {!r} occurs more than once.".format(name)) seen_lhs[name] = token indices = self._consume().data if self._next.type == 'indices' else '' for i, index in enumerate(indices): if index in indices[i+1:]: raise _IntermediateError('Repeated indices are not allowed on the left hand side.') elif '0' <= index <= '9': raise _IntermediateError('Numeric indices are not allowed on the left hand side.') return self._get_arg(name, indices, self._current.pos) @highlight def parse_var(self): 'parse a component of a term, e.g. "1", "a_i", "(2 a_i)", "a_i^2", "abs(x)"' if self._next.type == '(': self._consume() value = self.parse_subexpression() self._consume_assert_equal(')') value = value.replace(ast=('group', value.ast)) elif self._next.type == '[': self._consume() value = self.parse_subexpression() self._consume_assert_equal(']') value = value.replace(ast=('jump', value.ast)) if self._next.type == 'geometry': geometry_name = self._consume().data else: geometry_name = self.default_geometry_name geom = self._get_geometry(geometry_name) if self._next.type == 'indices': indices = self._consume().data value *= _Array.wrap(('normal', _(geom)), indices, geom.shape) elif self._next.type == '{': self._consume() value = self.parse_subexpression() self._consume_assert_equal('}') value = value.replace(ast=('mean', value.ast)) elif self._next.type == '<': self._consume() args = self.parse_comma_separated(end='>', parse_item=self.parse_subexpression) indices = self._consume() if indices.type != 'indices': raise _IntermediateError('Expected 1 index.', at=indices.pos, count=len( if len( != 1: raise _IntermediateError('Expected 1 index, got {}.'.format(len(, at=indices.pos, count=len( if '0' <= <= '9': raise _IntermediateError('Expected a non-numeric index, got {!r}.'.format(, at=indices.pos, count=len( value = _Array.stack(args, elif self._next.type == 'eye': self._consume() if self._next.type == 'indices': indices = self._consume().data else: indices = '' length = _Length(self._current.pos) value = _Array.wrap(('eye', _(length)), indices, (length, length)) elif self._next.type == 'normal': self._consume() if self._next.type == 'geometry': geometry_name = self._consume().data else: geometry_name = self.default_geometry_name geom = self._get_geometry(geometry_name) if self._next.type == 'indices': indices = self._consume().data else: indices = '' value = _Array.wrap(('normal', _(geom)), indices, geom.shape) elif self._next.type == 'variable': token = self._consume() name = if name in self.functions and name not in self.variables: # function (and not overriden as variable) self._consume_assert_equal('(', msg="Expected '(' for function {}.".format(name)) args = self.parse_comma_separated(end=')', parse_item=self.parse_subexpression) nargs = self.functions[name] if len(args) != nargs: raise _IntermediateError('Function {!r} takes {}, got {}.'.format(name, _sp(nargs, 'argument', 'arguments'), len(args))) args = _Array.align(*args) value = args[0].replace(ast=('call', _(name))+tuple(arg.ast for arg in args)) elif name.startswith('?'): indices = self._consume().data if self._next.type == 'indices' else '' value = self._get_arg(name[1:], indices, self._current.pos) else: raw = self._get_variable(name) indices = self._consume().data if self._next.type == 'indices' else '' value = _Array.wrap(_(raw), indices, raw.shape) else: raise _IntermediateError('Expected a variable, group or function call.') if self._next.type == 'gradient': token = self._consume() if',?'): name =[2:] if '_' in name: name, indices = name.split('_', 1) indices_start = token.pos+3+len(name) else: indices = '' indices_start = 0 arg = self._get_arg(name, indices, indices_start) value = value.derivative(arg, indices) else: gradtype = {',': 'grad', ';': 'surfgrad'}[[0]] if '_' in[1:]: geometry_name, indices =[1:].split('_', 1) else: geometry_name = self.default_geometry_name indices =[1:] geom = self._get_geometry(geometry_name) for i, index in enumerate(indices): value = value.grad(index, geom, gradtype) elif self._next.type == 'indices': raise _IntermediateError("Indices can only be specified for variables, e.g. 'a_ij', not for groups, e.g. '(a+b)_ij'.", at=self._next.pos, count=len( if self._next.type == '(': self._consume() subs = self.parse_comma_separated(end=')', parse_item=functools.partial(self.parse_substitution, seen_lhs={})) if not subs: raise _IntermediateError("Zero substitutions are not allowed.") ast = ['substitute', value.ast] links = [] for lhs, rhs in subs: ast += [lhs.ast, rhs.ast] links += [rhs.linked_lengths, frozenset(zip(lhs.shape, rhs.shape))] value = value.replace(ast=ast, linked_lengths=value._join_lengths(*links)) if self._next.type == '^': token = self._consume() if self._next.type == '(': self._consume() exponent = self.parse_subexpression() self._consume_assert_equal(')') else: if self._next.type == '-': self._consume() negate = True else: negate = False exponent = self.parse_const_scalar() if negate: exponent = -exponent value = value**exponent return value @highlight def parse_const_scalar(self): 'parse a constant scalar, e.g. "1", "1.0", "0.1"' token = self._consume() if token.type == 'int': value = _Array.wrap(_(int(, '', []) elif token.type == 'float': value = _Array.wrap(_(float(, '', []) else: raise _IntermediateError('Expected a number.') if self._next.type == 'gradient': self._consume() raise _IntermediateError('Taking a derivative of a constant is not allowed.') return value @highlight def parse_const(self): 'parse a const, possibly with indices, e.g. "1_j"' value = self.parse_const_scalar() if self._next.type == 'indices': token = self._consume() indices = for i, index in enumerate(indices): if '0' <= index <= '9': raise _IntermediateError('Numeric indices are not allowed on constant values.') if index in indices[1+i:]: raise _IntermediateError('Indices of a constant value may not be repeated.') value = value.append_axis(index, _Length(pos=token.pos+i)) if self._next.type == 'gradient': self._consume() raise _IntermediateError('Taking a derivative of a constant is not allowed.') return value @highlight def parse_numerator(self): 'parse the numerator part of a fraction' if self._next.type in ('int', 'float'): value = self.parse_const() else: value = self.parse_var() while True: stop = self._next.pos if self._next_non_whitespace.type in (')', ']', '}', '>', 'EOF', '+', '-', '/', '|', ','): break self._consume_assert_whitespace() value *= self.parse_var() return value @highlight def parse_denominator(self): 'parse the denominator part of a fraction' value = self.parse_numerator() if value.ndim > 0: raise _IntermediateError('A denominator must have dimension 0.') return value def parse_comma_separated(self, end, parse_item): 'parse comma separated values until end token, e.g. "1, 2 (a_ij b_j + 3))" with end token ")"' items = [] self._consume_if_whitespace() if self._next.type != end: while True: items.append(parse_item()) self._consume_if_whitespace() if self._next.type != ',': break self._consume_assert_equal(',') self._consume_assert_whitespace() self._consume_assert_equal(end) return items @highlight def parse_substitution(self, seen_lhs): 'parse a substitution, e.g. "x_ij=a_ij" in "?x_kk(x_ij=a_ij)"' lhs = self.parse_lhs_arg(seen_lhs) self._consume_if_whitespace() self._consume_assert_equal('=') self._consume_if_whitespace() rhs = self.parse_subexpression() if set(lhs.indices) != set(rhs.indices): raise _IntermediateError('Left and right hand side should have the same indices, got {!r} and {!r}.'.format(lhs.indices, rhs.indices)) rhs = rhs.transpose(lhs.indices) return lhs, rhs @highlight def parse_term(self): 'parse a term, e.g. "a b_i (2 c_i + 1)"' value = self.parse_numerator() if self._next_non_whitespace.type == '/': self._consume_assert_whitespace() token = self._consume() assert token.type == '/' self._consume_assert_whitespace() denominator = self.parse_denominator() value /= denominator return value @highlight def parse_subexpression(self): 'parse a scope: the entire expression or a subexpression between parentheses' self._consume_if_whitespace() negate = self._next.type == '-' if negate: self._consume() self._consume_if_whitespace() value = self.parse_term() if negate: value = -value while self._next_non_whitespace.type not in ('|', 'EOF', '_', ')', ']', '}', '>', ','): self._consume_assert_whitespace() op_token = self._consume() if op_token.type not in '+-': raise _IntermediateError('Expected {!r} or {!r}.'.format('+', '-'), at=op_token.pos, count=len( self._consume_assert_whitespace() r_value = self.parse_term() value = {'+': value.__add__, '-': value.__sub__}[op_token.type](r_value) self._consume_if_whitespace() return value @highlight def tokenize(self): 'subdivide :attr:`expression` in indivisible tokens' pos = 0 tokens = [_Token('BOF', '', pos)] while pos < len(self.expression): m = re.match(r'\s+', self.expression[pos:]) if m: tokens.append(_Token('whitespace',, pos)) pos += m.end() continue if self.expression[pos] in '+-^/|=[]{}()<>,': tokens.append(_Token(self.expression[pos], self.expression[pos], pos)) pos += 1 continue m_variable = re.match(r'[?]?[a-zA-Zα-ωΑ-Ω][a-zA-Zα-ωΑ-Ω0-9]*', self.expression[pos:]) m_variable = if m_variable else '' m_eye = _string_startswith(self.expression, self.eye_symbols, start=pos) # Insert eye or normal symbols only if we can't match a longer variable name. if m_eye and len(m_variable) <= len(m_eye): tokens.append(_Token('eye', m_eye, pos)) pos += len(m_eye) continue m_normal = _string_startswith(self.expression, self.normal_symbols, start=pos) if m_normal and len(m_variable) <= len(m_normal): tokens.append(_Token('normal', m_normal, pos)) pos += len(m_normal) continue if m_variable: tokens.append(_Token('variable', m_variable, pos)) pos += len(m_variable) continue m = re.match(r'[0-9]*[.][0-9]*', self.expression[pos:]) if m: if'0') and not'0.'): raise _IntermediateError('Leading zeros are forbidden.', at=pos, count=len( tokens.append(_Token('float',, pos)) pos += m.end() continue m = re.match(r'[0-9]+', self.expression[pos:]) if m: if'0') and not == '0': raise _IntermediateError('Leading zeros are forbidden.', at=pos, count=len( tokens.append(_Token('int',, pos)) pos += m.end() continue if self.expression[pos] == '_': pos += 1 parts = 0 m = re.match(r'[a-zA-Zα-ωΑ-Ω][a-zA-Zα-ωΑ-Ω0-9]*_', self.expression[pos:]) if m: withgeom =[:-1] tokens.append(_Token('geometry',[:-1], pos)) pos += m.end() else: withgeom = None m = re.match(r'[a-zA-Z0-9]+', self.expression[pos:]) if m: tokens.append(_Token('indices',, pos)) pos += m.end() parts += 1 m_arg = re.match(r',[?][a-zA-Zα-ωΑ-Ω][a-zA-Zα-ωΑ-Ω0-9]*(_[a-zA-Z0-9]+)?', self.expression[pos:]) m_geom = re.match(r'[,;]([a-zA-Zα-ωΑ-Ω][a-zA-Zα-ωΑ-Ω0-9]*_)?([a-zA-Z0-9]+)', self.expression[pos:]) if m_arg: tokens.append(_Token('gradient',, pos)) pos += m_arg.end() parts += 1 elif m_geom: if withgeom is not None and not variant_geom =[0] + withgeom + '_' + variant_default =[0] + self.default_geometry_name + '_' + raise _IntermediateError('Missing geometry, e.g. {!r} or {!r}.'.format(variant_geom, variant_default), at=pos) tokens.append(_Token('gradient',, pos)) pos += m_geom.end() parts += 1 if parts == 0: raise _IntermediateError('Missing indices.', at=pos) continue raise _IntermediateError('Unknown symbol: {!r}.'.format(self.expression[pos]), at=pos) tokens.append(_Token('EOF', '', pos)) self._tokens = tokens self._index = 0 def _string_startswith(string, prefixes, start=0): assert not isinstance(prefixes, str) for prefix in prefixes: if string.startswith(prefix, start): return prefix def _replace_lengths(ast, lengths): 'replace all :class:`_Length` objects in ``ast`` with the lengths in ``lengths``' if ast[0] is not None: return (ast[0],) + tuple(_replace_lengths(arg, lengths) for arg in ast[1:]) elif isinstance(ast[1], _Length): return _(lengths[ast[1]]) else: return ast
[docs]def parse(expression, variables, functions, indices, arg_shapes={}, default_geometry_name='x'): '''Parse ``expression`` and return AST. This function parses a tensor expression with `Einstein Summation Convection`_ stored in a :class:`str` and returns an Abstract Syntax Tree (AST). The syntax of ``expression`` is as follows: * **Integers** or **decimal numbers** are denoted in the usual way. Examples: ``1``, ``1.2``, ``.2``. A number may not start with a zero, except when followed by a dot: ``0.1`` is valid, but ``01`` is not. * **Variables** are denoted with a string of alphanumeric characters. The first character may not be a numeral. Unlike Python variables, underscores are not allowed, as they have a special meaning. If the variable is an array with one or more axes, all those axes should be labeled with a latin character, the index, and appended to the variable with an underscore. For example an array ``a`` with two axes can be denoted with ``a_ij``. Optionally, a single numeral may be used to select an item at the concerning axis. Example: in ``a_i0`` the first axis of ``a`` is labeled ``i`` and the first element of the second axis is selected. If the same index occurs twice, the trace is taken along the concerning axes. Example: the trace of the first and third axes of ``b`` is denoted by ``b_iji``. It is invalid to specify an index more than twice. The following names cannot be used as variables: ``n``, ``δ``, ``$``. The variable named ``x``, or the value of argument ``default_geometry_name``, has a special meaning, detailed below. * A term, the **product** of two or more arrays or scalars, is denoted by space-separated variables, constants or compound expressions. Example: ``a b c`` denotes the product of the scalars ``a``, ``b`` and ``c``. A term may start with a number, but a number is not allowed in other parts of the term. Example: ``2 a`` denotes two times ``a``; ``2 2 a`` and ``2 a 2``` are invalid. When two arrays in a term have the same index, this index is summed. Example: ``a_i b_i`` denotes the inner product of ``a`` and ``b`` and ``A_ij b_j``` a matrix vector product. It is not allowed to use an index more than twice in a term. * The operator ``/`` denotes a **fraction**. Example: in ``a b / c d`` ``a b`` is the numerator and ``c d`` the denominator. Both the numerator and the denominator may start with a number. Example: ``2 a / 3 b``. The denominator must be a scalar. Example: ``2 / a_i b_i`` is valid, but ``2 a_i / b_i`` is not. .. warning:: This syntax is different from the Python syntax. In Python ``a*b / c*d`` is mathematically equivalent to ``a*b*d/c``. * The operators ``+`` and ``-`` denote **add** and **subtract**. Both operators should be surrounded by whitespace, e.g. ``a + b``. Both operands should have the same shape. Example: ``a_ij + b_i c_j`` is a valid, provided that the lengths of the axes with the same indices match, but ``a_ij + b_i`` is invalid. At the beginning of an expression or a compound ``-`` may be used to negate the following term. Example: in ``-a b + c`` the term ``a b`` is negated before adding ``c``. It is not allowed to negate other terms: ``a + -b`` is invalid, so is ``a -b``. * An expression surrounded by parentheses is a **compound expression** and can be used as single entity in a term. Example: ``(a_i + b_i) c_i`` denotes the inner product of ``a_i + b_i`` with ``c_i``. * **Exponentiation** is denoted by a ``^``, where the left and right operands should be a number, variable or compound expression and the right operand should be a scalar. Example: ``a^2`` denotes the square of ``a``, ``a^-2`` denotes ``a`` to the power ``-2`` and ``a^(1 / 2)`` the square root of ``a``. * An **argument** is denoted by a name — following the same rules as a variable name — prefixed with a question mark. An argument is a scalar or array with a yet unknown value. Example: ``basis_i ?coeffs_i`` denotes the inner product of a basis with unknown coefficient vector ``?coeffs``. If possible the shape of the argument is deduced from the expression. In the previous example the shape of ``?coeffs`` is equal to the shape of ``basis``. If the shape cannot be deduced from the expression the shape should be defined manually (see :func:`parse`). Arguments and variables live in separate namespaces: ``?x`` and ``x`` are different entities. * An argument may be **substituted** by appending without whitespace ``(arg = value)`` to a variable of compound expression, where ``arg`` is an argument and ``value`` the substitution. The substitution applies to the variable of compound expression only. The value may be an expression. Example: ``2 ?x(x = 3 + y)`` is equivalent to ``2 (3 + y)`` and ``2 ?x(x=y) + 3`` is equivalent to ``2 (y) + 3``. It is possible to apply multiple substitutions. Example: ``(?x + ?y)(x = 1, y = )2`` is equivalent to ``1 + 2``. * The **gradient** of a variable to the default geometry — the default geometry is variable ``x`` unless overriden by the argument ``default_geometry_name`` — is denoted by an underscore, a comma and an index. If the variable is an array with more than one axis, the underscore is omitted. Example: ``a_,i`` denotes the gradient of the scalar ``a`` to the geometry and ``b_i,j`` the gradient of vector ``b``. The gradient of a compound expression is denoted by an underscore, a comma and an index. Example: ``(a_i + b_j)_,k`` denotes the gradient of ``a_i + b_j``. The usual summation rules apply and it is allowed to use a numeral as index. The **surface gradient** is denoted with a semicolon instead of a comma, but follows the same rules as the gradient otherwise. Example: ``a_i;j`` is the sufrace gradient of ``a_i`` to the geometry. It is also possible to take the gradient to another geometry by appending the name of the geometry, which should exist as a variable, and an underscore directly after the comma of semicolon. Example: ``a_i,altgeom_j`` denotes the gradient of ``a_i`` to ``altgeom`` and the gradient axis has index ``j``. Futhermore, it is possible to take the **derivative** to an argument by adding the argument with appropriate indices after the comma. Example: ``(?x^2)_,?x`` denotes the derivative of ``?x^2`` to ``?x``, which is equivalent to ``2 ?x``, and ``(?y_i ?y_i),?y_j`` is the derivative of ``?y_i ?y_i`` to ``?y_j``, which is equivalent to ``2 ?y_j``. * The **normal** of the default geometry is denoted by ``n_i``, where the index ``i`` may be replaced with an index of choice. The normal with respect to different geometry is denoted by appending an underscore with the name of the geometry right after ``n``. Example: ``n_altgeom_j`` is the normal with respect to geometry ``altgeom``. * A **dirac** is denoted by ``δ`` or ``$`` and takes two indices. The shape of the dirac is deduced from the expression. Example: let ``A`` be a square matrix with three rows and columns, then ``δ_ij`` in ``(A_ij - λ δ_ij) x_j`` has three rows and columns as well. * An expression surrounded by square brackets or curly braces denotes the **jump** or **mean**, respectively, of the enclosed expression. Example: ``[ a_i ]`` denotes the jump of ``a_i`` and ``{ a_i + b_i }`` denotes the mean of ``a_i + b_i``. * A **function call** is denoted by a name — following the same rules as for a variable name — directly followed by the left parenthesis ``(``, without a space. The arguments to the function are separated by a comma and at least one space. The function is applied pointwise to the arguments and all arguments should have the same shape. Example: ``f(x_i, y_i)``.denotes the call to function ``f`` with arguments ``x_i`` and ``y_i``. Functions and variables share a namespace: defining a variable with the same name as a function renders the function inaccessible. * A **stack** of two or more arrays along an axis is denoted by a ``<`` followed by comma and space separated arrays followed by ``>`` and an index. If an argument does not have an axis with the specified stack index, the argument is expanded with an axis of length one. Beside the stack axis, all arguments should have the same shape. Example: ``<1, x_i>_i``, with ``x`` a vector of length three, creates an array with components ``1``, ``x_0``, ``x_1``, ``x_2``. .. _`Einstein Summation Convection`: Args ---- expression : :class:`str` The expression to parse. See :mod:`~nutils.expression` for the expression syntax. variables : :class:`dict` of :class:`str` and :class:`nutils.function.Array` pairs A :class:`dict` of variable names and array pairs. All variables used in the ``expression`` should exist in ``variables``. functions : :class:`dict` of :class:`str` and :class:`int` pairs A :class:`dict` of function names and number of arguments pairs. All functions used in the ``expression`` should exist in ``functions``. indices : :class:`str` The indices used for aligning the resulting array. For example, let ``expression`` be ``'a_ij'``. If ``indices`` is ``'ij'``, then the returned array is simply ``variables['a']``, but if ``indices`` is ``'ji'`` the transpose of ``variables['a']`` is returned. All indices of the ``expression`` should be listed precisely once. arg_shapes : :class:`dict` of :class:`str` and :class:`tuple` or :class:`int`\\s pairs A :class:`dict` of argument names and shapes. If ``expression`` contains an argument not present in ``arg_shapes`` the shape will be decuded from the expression and added to a copy of ``arg_shapes``. default_geometry_name : :class:`str` The name of the default geometry variable. When computing a gradient or the normal, e.g. ``'f_,i'`` or ``'n_i'``, this variable is used as the geometry, unless the geometry is explicitly mentioned in the expression. Default: ``'x'``. Returns ------- ast : :class:`tuple` The parsed ``expression`` as an abstract syntax tree (AST). The AST is a :class:`tuple` of an opcode and arguments. The special opcode ``None`` indicates that the single argument is used verbatim. All other opcodes have AST as arguments. The following opcodes exist:: (None, const) ('group', group) ('arg', name, *shape) ('substitute', array, arg, value) ('call', func, arg) ('eye', length) ('normal', geom) ('getitem', array, dim, index) ('trace', array, n1, n2) ('sum', array, axis) ('concatenate', *args) ('grad', array, geom) ('surfgrad', array, geom) ('derivative', func, target) ('append_axis', array, length) ('transpose', array, trans) ('jump', array) ('mean', array) ('neg', array) ('add', left, right) ('sub', left, right) ('mul', left, right) ('truediv', left, right) ('pow', left, right) arg_shapes : :class:`dict` of :class:`str` and :class:`tuple` of :class:`int`\\s pairs A copy of ``arg_shapes`` updated with shapes of arguments present in this ``expression``. ''' parser = _ExpressionParser(expression, variables, functions, arg_shapes, default_geometry_name) parser.tokenize() value = parser.parse_subexpression() parser._consume_assert_equal('EOF', msg='Unexpected symbol at end of expression.') if indices is None: if value.ndim > 1: raise AmbiguousAlignmentError( 'Cannot unambiguously align the array because the array has more than one dimension.\n' + expression + '\n' + '^'*len(expression)) ast = value.ast else: try: ast = value.transpose(indices).ast except _IntermediateError as e: raise ExpressionSyntaxError(e.msg + '\n' + expression + '\n' + '^'*len(expression)) from e lengths = {} undetermined = set() for group in value.linked_lengths: val = None for i in group: if not isinstance(i, _Length): assert val is None val = i if val is None: undetermined.update(i.pos for i in group) else: lengths.update((length, val) for length in group) for pos in sorted(undetermined): raise ExpressionSyntaxError('Length of axis cannot be determined from the expression.' + '\n' + expression + '\n' + ' '*pos + '^') arg_shapes = dict(arg_shapes) for arg, shape in parser.arg_shapes.items(): arg_shapes[arg] = tuple(lengths.get(i, i) for i in shape) return _replace_lengths(ast, lengths), arg_shapes
# vim:shiftwidth=2:softtabstop=2:expandtab:foldmethod=indent:foldnestmax=2