3356190533

Build # 3356190533

Build Type

github

Committed by GitHub

Commit Message

<a href="https://github.com/tBuLi/kingdon/commit/<a class=hub.com/tBuLi/kingdon/commit/8cf9b624d4b40f8a3a912abd2c391e72b2176dc5">8cf9b624d<a href="https://github.com/tBuLi/kingdon/commit/8cf9b624d4b40f8a3a912abd2c391e72b2176dc5">">Merge </a><a class="double-link" href="https://github.com/tBuLi/kingdon/commit/<a class="double-link" href="https://github.com/tBuLi/kingdon/commit/a10360b727c40e6b80e8caa5190320aa0c5e21ac">a10360b72</a>">a10360b72</a><a href="https://github.com/tBuLi/kingdon/commit/8cf9b624d4b40f8a3a912abd2c391e72b2176dc5"> into 238ae7ff6">238ae7ff6</a>

Pull Request Pull Request #9: Introduced OperatorDict, cleaned-up code significantly

Run Details

282 of 282 new or added lines in 5 files covered. (100.0%)

530 of 559 relevant lines covered (94.81%)

0.95 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

93.71

/kingdon/codegen.py

from itertools import product, combinations, groupby
from collections import namedtuple
from functools import reduce
import operator
import linecache
import warnings

from sympy import simplify, sympify, Add, Mul
from sympy.utilities.lambdify import lambdify
from numba import njit


TermTuple = namedtuple('TermTuple', ['key_out', 'keys_in', 'sign', 'values_in', 'termstr'])
TermTuple.__doc__ = """
TermTuple represents a single monomial in a product of multivectors.

:param key_out: is the basis blade to which this monomial belongs.
:param keys_in: are the input basis blades in this monomial.
:param sign: Sign of the monomial.
:param values_in: Input values. Typically sympy symbols.
:param termstr: The string representation of this monomial, e.g. '-x*y'.
"""


def term_tuple(items, sign_func, keyout_func=operator.xor):
    """
    Create a single term in a multivector product between the basis blades present in `items`.
    """
    keys_in, values_in = zip(*items)
    sign = reduce(operator.mul, (sign_func(pair) for pair in combinations(keys_in, r=2)))

    if not sign:
        return TermTuple(key_out=0, keys_in=keys_in, sign=sign, values_in=values_in, termstr='')
    # Values could also have signs, e.g. -x1. Multiply signs by these.
    sign_mod, values_in = zip(*(v.args if isinstance(v, Mul) else (1, v) for v in values_in))
    sign = reduce(operator.mul, sign_mod, sign)
    key_out = reduce(keyout_func, keys_in)
    return TermTuple(key_out=key_out,
                     keys_in=keys_in,
                     sign=sign,
                     values_in=values_in,
                     termstr=f'{"+" if sign > 0 else "-"}{"*".join(v.name for v in values_in)}')


def codegen_product(*mvs, name_base, filter_func=lambda tt: tt.sign, sign_func=None, keyout_func=operator.xor, asdict=False, sympy=False):
    """
    Helper function for the codegen of all product-type functions.

    :param mvs: Positional-argument :class:`~kingdon.multivector.MultiVector`.
    :param filter_func: A condition which should be true in the preprocessing of terms.
        For example, for the geometric product we filter out all values for which
        ei*ej = 0 since these do not have to be considered anyway.
        Input is a TermTuple.
    :param sign_func: function to compute sign between terms. E.g. algebra.signs[ei, ej]
        for metric dependent products. Input: 2-tuple of blade indices, e.g. (ei, ej).
    :param name_base: base name for the generated code.
    :param asdict: If true, return the dict of strings before converting to a function.
    """
    sortfunc = lambda x: x.key_out
    algebra = mvs[0].algebra
    if sign_func is None:
        sign_func = lambda pair: algebra.signs[pair]

    terms = filter(filter_func, (term_tuple(items, sign_func, keyout_func=keyout_func)
                                 for items in product(*(mv.items() for mv in mvs))))
    # TODO: Can we loop over the basis blades in such a way that no sort is needed?
    sorted_terms = sorted(terms, key=sortfunc)
    if not sympy:
        res = {k: "".join(term.termstr for term in group)
               for k, group in groupby(sorted_terms, key=sortfunc)}
    else:
        res = {k: Add(*(Mul(*(term.values_in if term.sign == 1 else (term.sign, *term.values_in)), evaluate=False)
                        for term in group), evaluate=False)
               for k, group in groupby(sorted_terms, key=sortfunc)}

    if asdict:
        return res
    elif not sympy:
        return _func_builder(res, *mvs, name_base=name_base)
    else:
        return _lambdify_binary(*mvs, res)


def codegen_gp(x, y, symbolic=False):
    """
    Generate the geometric product between :code:`x` and :code:`y`.

    :param x: Fully symbolic :class:`~kingdon.multivector.MultiVector`.
    :param y: Fully symbolic :class:`~kingdon.multivector.MultiVector`.
    :param symbolic: If true, return a symbolic :class:`~kingdon.multivector.MultiVector` instead of a lambda function.
    :return: tuple with integers indicating the basis blades present in the
        product in binary convention, and a lambda function that perform the product.
    """
    keys_out, func = codegen_product(x, y, name_base='gp')
    if symbolic:
        return x.fromkeysvalues(x.algebra, keys_out, func(x.values(), y.values()))

    return keys_out, func


def codegen_conj(x, y):
    if x.algebra.simplify:
        res = codegen_product(x, y, ~x, name_base='gp', asdict=True, sympy=True)
        res = {k: str(simp_expr) for k, expr in res.items() if (simp_expr := simplify(expr))}
    else:
        res = codegen_product(x, y, ~x, name_base='gp', asdict=True)
    return _func_builder(res, x, y, name_base="conj")


def codegen_cp(x, y, symbolic=False):
    """
    Generate the commutator product of `x := self` and `y := other`: `x.cp(y) = 0.5*(x*y-y*x)`.

    :return: tuple of keys in binary representation and a lambda function.
    """
    algebra = x.algebra
    filter_func = lambda tt: (algebra.signs[tt.keys_in] - algebra.signs[tt.keys_in[::-1]])
    keys_out, func = codegen_product(x, y, filter_func=filter_func, name_base='cp')
    if symbolic:
        return x.fromkeysvalues(algebra, keys_out, func(x.values(), y.values()))

    return keys_out, func


def codegen_acp(x, y, symbolic=False):
    """
    Generate the anti-commutator product of `x := self` and `y := other`: `x.acp(y) = 0.5*(x*y+y*x)`.

    :return: tuple of keys in binary representation and a lambda function.
    """
    algebra = x.algebra
    filter_func = lambda tt: (algebra.signs[tt.keys_in] + algebra.signs[tt.keys_in[::-1]])
    keys_out, func = codegen_product(x, y, filter_func=filter_func, name_base='acp')
    if symbolic:
        return x.fromkeysvalues(algebra, keys_out, func(x.values(), y.values()))

    return keys_out, func


def codegen_ip(x, y, diff_func=abs, symbolic=False):
    """
    Generate the inner product of `x := self` and `y := other`.

    :param diff_func: How to treat the difference between the binary reps of the basis blades.
        if :code:`abs`, compute the symmetric inner product. When :code:`lambda x: -x` this
        function generates left-contraction, and when :code:`lambda x: x`, right-contraction.
    :return: tuple of keys in binary representation and a lambda function.
    """
    algebra = x.algebra
    filter_func = lambda tt: tt.key_out == diff_func(tt.keys_in[0] - tt.keys_in[1])
    keys_out, func = codegen_product(x, y, filter_func=filter_func, name_base='ip')
    if symbolic:
        return x.fromkeysvalues(algebra, keys_out, func(x.values(), y.values()))

    return keys_out, func


def codegen_lc(x, y):
    """
    Generate the left-contraction of `x := self` and `y := other`.

    :return: tuple of keys in binary representation and a lambda function.
    """
    return codegen_ip(x, y, diff_func=lambda x: -x)


def codegen_rc(x, y):
    """
    Generate the right-contraction of `x := self` and `y := other`.

    :return: tuple of keys in binary representation and a lambda function.
    """
    return codegen_ip(x, y, diff_func=lambda x: x)


def codegen_sp(x, y):
    """
    Generate the scalar product of `x := self` and `y := other`.

    :return: tuple of keys in binary representation and a lambda function.
    """
    return codegen_ip(x, y, diff_func=lambda x: 0)


def codegen_proj(x, y):
    """
    Generate the projection of `x := self` onto `y := other`: :math:`(x \cdot y) / y`.

    :return: tuple of keys in binary representation and a lambda function.
    """
    filter_func = lambda tt: tt.key_out == abs(tt.keys_in[0] - tt.keys_in[1]) ^ tt.keys_in[2]
    res = codegen_product(x, y, ~y, name_base='gp', filter_func=filter_func, asdict=True)
    return _func_builder(res, x, y, name_base='proj')


def codegen_op(x, y, symbolic=False):
    """
    Generate the outer product of `x := self` and `y := other`: `x.op(y) = x ^ y`.

    :x: MultiVector
    :y: MultiVector
    :return: dictionary with integer keys indicating the corresponding basis blade in binary convention,
        and values which are a 3-tuple of indices in `x`, indices in `y`, and a lambda function.
    """
    algebra = x.algebra
    filter_func = lambda tt: tt.key_out == sum(tt.keys_in)
    sign_func = lambda pair: (-1)**algebra.swaps[pair]
    keys_out, func = codegen_product(x, y, filter_func=filter_func, sign_func=sign_func, name_base='op')
    if symbolic:
        return x.fromkeysvalues(algebra, keys_out, func(x.values(), y.values()))

    return keys_out, func


def codegen_rp(x, y):
    """
    Generate the regressive product of :code:`x` and :code:`y`:,
    :math:`x \\vee y`.

    :param x:
    :param y:
    :return: tuple of keys in binary representation and a lambda function.
    """
    algebra = x.algebra
    keyout_func = lambda tot, key_in: len(algebra) - 1 - (key_in ^ tot)
    filter_func = lambda tt: len(algebra) - 1 == sum(tt.keys_in) - tt.key_out
    # Sign is composed of dualization of each blade, exterior product, and undual.
    sign_func = lambda pair: (
        algebra.signs[pair[0], len(algebra) - 1 - pair[0]] *
        algebra.signs[pair[1], len(algebra) - 1 - pair[1]] *
        (-1)**algebra.swaps[len(algebra) - 1 - pair[0], len(algebra) - 1 - pair[1]] *
        algebra.signs[len(algebra) - 1 - (pair[0] ^ pair[1]), pair[0] ^ pair[1]]
    )

    return codegen_product(
        x, y,
        filter_func=filter_func,
        keyout_func=keyout_func,
        sign_func=sign_func,
        name_base='rp'
    )


Fraction = namedtuple('Fraction', ['numer', 'denom'])
Fraction.__doc__ = """
Tuple representing a fraction.
"""

def codegen_inv(x, symbolic=False):
    """
    Generate code for the inverse of :code:`x`.
    Currently, this always uses the Shirokov inverse, which is works in any algebra,
    but it can be expensive to compute.
    In the future this should be extended to use dedicated solutions for known cases.
    """
    k = 2 ** ((x.algebra.d + 1) // 2)
    x_i = x
    if x.grades == (0,):
        adj_x = 1
    else:
        for i in range(1, k + 1):
            # Sympify ratio to keep the ratios exact and avoid floating point errors.
            c_i = (sympify(k) / i) * x_i[0] if x_i[0] else x_i[0]
            adj_x = (x_i - c_i)
            if x.algebra.simplify:
                keys, values = zip(*((k, simp_expr) for k, expr in adj_x.items() if (simp_expr := simplify(expr))))
                adj_x = adj_x.fromkeysvalues(adj_x.algebra, keys=keys, values=values)
            x_i = x * adj_x
            if x_i:
                if x.algebra.simplify:
                    keys, values = zip(*((k, simp_expr) for k, expr in x_i.items() if (simp_expr := simplify(expr))))
                    x_i = x_i.fromkeysvalues(x_i.algebra, keys=keys, values=values)

                if x_i.grades == (0,):
                    break
            else:
                break

    if symbolic:
        return Fraction(adj_x, x_i[0])
    xinv = x.algebra.multivector({k: v / x_i[0] for k, v in adj_x.items()})
    return _lambdify_unary(x, xinv)


def codegen_div(x, y):
    """
    Generate code for :math:`x y^{-1}`.
    """
    adjy, normsqy = codegen_inv(y, symbolic=True)
    if not normsqy:
        raise ZeroDivisionError
    x_adjy = x * adjy
    if x.algebra.simplify:
        xdivy = x.algebra.multivector({k: simplify(v / normsqy) for k, v in x_adjy.items()})
    else:
        xdivy = x.algebra.multivector({k: v / normsqy for k, v in x_adjy.items()})
    return _lambdify_binary(x, y, xdivy)

def codegen_normsq(x):
    if x.algebra.simplify:
        res = codegen_product(x, ~x, name_base='gp', asdict=True, sympy=True)
        res = {k: str(simp_expr) for k, expr in res.items() if (simp_expr := simplify(expr))}
    else:
        res = codegen_product(x, ~x, name_base='gp', asdict=True)
    return _func_builder(res, x, name_base="normsq")


def codegen_outerexp(x, asterms=False):
    if len(x.grades) != 1:
        warnings.warn('Outer exponential might not converge for mixed-grade multivectors.', RuntimeWarning)
    k = x.algebra.d

    Ws = [x.algebra.scalar([1]), x]
    j = 2
    while j <= k:
        Wj = Ws[-1] ^ x
        # Dividing like this avoids floating point numbers, which is excellent.
        Wj._values = tuple(v / j for v in Wj._values)
        if Wj:
            Ws.append(Wj)
            j += 1
        else:
            break

    if asterms:
        return Ws
    L = reduce(operator.add, Ws)
    return _func_builder(dict(L.items()), x, name_base='outerexp')


def _lambdify_binary(x, y, x_bin_y):
    xy_symbols = [list(x.values()), list(y.values())]
    func = lambdify(xy_symbols, list(x_bin_y.values()), cse=x.algebra.cse)
    return tuple(x_bin_y.keys()), njit(func) if x.algebra.numba else func


def _lambdify_unary(x, x_unary):
    func = lambdify([list(x.values())], list(x_unary.values()), cse=x.algebra.cse)
    return tuple(x_unary.keys()), njit(func) if x.algebra.numba else func


def _lambdify_mv(free_symbols, mv):
    # TODO: Numba wants a tuple in the line below, but simpy only produces a
    #  list as output if this is a list, not a tuple. See if we can solve this.
    func = lambdify(free_symbols, list(mv.values()), cse=mv.algebra.cse)
    return tuple(mv.keys()), njit(func) if mv.algebra.numba else func

def _func_builder(res_vals: dict, *mvs, name_base: str):
    """
    Build a Python functions for the product between given multivectors.

    :param res_vals: Dict to be converted into a function. The keys correspond to the basis blades in binary,
        while the values are strings to be converted into source code.
    :param mvs: all the multivectors that the resulting function is a product of.
    :param name_base: Base for the name of the function. E.g. 'gp'.
    :return: tuple of output keys of the callable, and the callable.
    """
    args = [f'arg{i}' for i in range(1, len(mvs) + 1)]
    keys_str_per_arg = ["_".join(str(k) for k in mv.keys()) for mv in mvs]
    func_name = f'{name_base}_{"_x_".join(keys_str_per_arg)}'
    header = f'def {func_name}({", ".join(args)}):'
    if res_vals:
        body = "\n".join(f'    {",".join(v.name for v in mv.values())}, = {arg}' for mv, arg in zip(mvs, args))
        return_val = f'    return ({", ".join(v if isinstance(v, str) else str(v) for v in res_vals.values())},)'
    else:
        body = ''
        return_val = f'    return tuple()'
    func_source = f'{header}\n{body}\n{return_val}'

    # Dynamically build a function
    func_locals = {}
    c = compile(func_source, func_name, 'exec')
    exec(c, {}, func_locals)

    # Add the generated code to linecache such that it is inspect-safe.
    linecache.cache[func_name] = (len(func_source), None, func_source.splitlines(True), func_name)
    func = func_locals[func_name]
    return tuple(res_vals.keys()), njit(func) if mvs[0].algebra.numba else func

1	from itertools import product, combinations, groupby	1✔
2	from collections import namedtuple	1✔
3	from functools import reduce	1✔
4	import operator	1✔
5	import linecache	1✔
6	import warnings	1✔
7
8	from sympy import simplify, sympify, Add, Mul	1✔
9	from sympy.utilities.lambdify import lambdify	1✔
10	from numba import njit	1✔
11
12
13	TermTuple = namedtuple('TermTuple', ['key_out', 'keys_in', 'sign', 'values_in', 'termstr'])	1✔
14	TermTuple.__doc__ = """	1✔
15	TermTuple represents a single monomial in a product of multivectors.
16
17	:param key_out: is the basis blade to which this monomial belongs.
18	:param keys_in: are the input basis blades in this monomial.
19	:param sign: Sign of the monomial.
20	:param values_in: Input values. Typically sympy symbols.
21	:param termstr: The string representation of this monomial, e.g. '-x*y'.
22	"""
23
24
25	def term_tuple(items, sign_func, keyout_func=operator.xor):	1✔
26	"""
27	Create a single term in a multivector product between the basis blades present in `items`.
28	"""
29	keys_in, values_in = zip(*items)	1✔
30	sign = reduce(operator.mul, (sign_func(pair) for pair in combinations(keys_in, r=2)))	1✔
31
32	if not sign:	1✔
33	return TermTuple(key_out=0, keys_in=keys_in, sign=sign, values_in=values_in, termstr='')	1✔
34	# Values could also have signs, e.g. -x1. Multiply signs by these.
35	sign_mod, values_in = zip(*(v.args if isinstance(v, Mul) else (1, v) for v in values_in))	1✔
36	sign = reduce(operator.mul, sign_mod, sign)	1✔
37	key_out = reduce(keyout_func, keys_in)	1✔
38	return TermTuple(key_out=key_out,	1✔
39	keys_in=keys_in,
40	sign=sign,
41	values_in=values_in,
42	termstr=f'{"+" if sign > 0 else "-"}{"*".join(v.name for v in values_in)}')
43
44
45	def codegen_product(*mvs, name_base, filter_func=lambda tt: tt.sign, sign_func=None, keyout_func=operator.xor, asdict=False, sympy=False):	1✔
46	"""
47	Helper function for the codegen of all product-type functions.
48
49	:param mvs: Positional-argument :class:`~kingdon.multivector.MultiVector`.
50	:param filter_func: A condition which should be true in the preprocessing of terms.
51	For example, for the geometric product we filter out all values for which
52	ei*ej = 0 since these do not have to be considered anyway.
53	Input is a TermTuple.
54	:param sign_func: function to compute sign between terms. E.g. algebra.signs[ei, ej]
55	for metric dependent products. Input: 2-tuple of blade indices, e.g. (ei, ej).
56	:param name_base: base name for the generated code.
57	:param asdict: If true, return the dict of strings before converting to a function.
58	"""
59	sortfunc = lambda x: x.key_out	1✔
60	algebra = mvs[0].algebra	1✔
61	if sign_func is None:	1✔
62	sign_func = lambda pair: algebra.signs[pair]	1✔
63
64	terms = filter(filter_func, (term_tuple(items, sign_func, keyout_func=keyout_func)	1✔
65	for items in product(*(mv.items() for mv in mvs))))
66	# TODO: Can we loop over the basis blades in such a way that no sort is needed?
67	sorted_terms = sorted(terms, key=sortfunc)	1✔
68	if not sympy:	1✔
69	res = {k: "".join(term.termstr for term in group)	1✔
70	for k, group in groupby(sorted_terms, key=sortfunc)}
71	else:
72	res = {k: Add((Mul((term.values_in if term.sign == 1 else (term.sign, *term.values_in)), evaluate=False)	1✔
73	for term in group), evaluate=False)
74	for k, group in groupby(sorted_terms, key=sortfunc)}
75
76	if asdict:	1✔
77	return res	1✔
78	elif not sympy:	1✔
79	return _func_builder(res, *mvs, name_base=name_base)	1✔
80	else:
81	return _lambdify_binary(*mvs, res)	×
82
83
84	def codegen_gp(x, y, symbolic=False):	1✔
85	"""
86	Generate the geometric product between :code:`x` and :code:`y`.
87
88	:param x: Fully symbolic :class:`~kingdon.multivector.MultiVector`.
89	:param y: Fully symbolic :class:`~kingdon.multivector.MultiVector`.
90	:param symbolic: If true, return a symbolic :class:`~kingdon.multivector.MultiVector` instead of a lambda function.
91	:return: tuple with integers indicating the basis blades present in the
92	product in binary convention, and a lambda function that perform the product.
93	"""
94	keys_out, func = codegen_product(x, y, name_base='gp')	1✔
95	if symbolic:	1✔
96	return x.fromkeysvalues(x.algebra, keys_out, func(x.values(), y.values()))	×
97
98	return keys_out, func	1✔
99
100
101	def codegen_conj(x, y):	1✔
102	if x.algebra.simplify:	1✔
103	res = codegen_product(x, y, ~x, name_base='gp', asdict=True, sympy=True)	1✔
104	res = {k: str(simp_expr) for k, expr in res.items() if (simp_expr := simplify(expr))}	1✔
105	else:
106	res = codegen_product(x, y, ~x, name_base='gp', asdict=True)	×
107	return _func_builder(res, x, y, name_base="conj")	1✔
108
109
110	def codegen_cp(x, y, symbolic=False):	1✔
111	"""
112	Generate the commutator product of `x := self` and `y := other`: `x.cp(y) = 0.5(xy-y*x)`.
113
114	:return: tuple of keys in binary representation and a lambda function.
115	"""
116	algebra = x.algebra	1✔
117	filter_func = lambda tt: (algebra.signs[tt.keys_in] - algebra.signs[tt.keys_in[::-1]])	1✔
118	keys_out, func = codegen_product(x, y, filter_func=filter_func, name_base='cp')	1✔
119	if symbolic:	1✔
120	return x.fromkeysvalues(algebra, keys_out, func(x.values(), y.values()))	×
121
122	return keys_out, func	1✔
123
124
125	def codegen_acp(x, y, symbolic=False):	1✔
126	"""
127	Generate the anti-commutator product of `x := self` and `y := other`: `x.acp(y) = 0.5(xy+y*x)`.
128
129	:return: tuple of keys in binary representation and a lambda function.
130	"""
131	algebra = x.algebra	1✔
132	filter_func = lambda tt: (algebra.signs[tt.keys_in] + algebra.signs[tt.keys_in[::-1]])	1✔
133	keys_out, func = codegen_product(x, y, filter_func=filter_func, name_base='acp')	1✔
134	if symbolic:	1✔
135	return x.fromkeysvalues(algebra, keys_out, func(x.values(), y.values()))	×
136
137	return keys_out, func	1✔
138
139
140	def codegen_ip(x, y, diff_func=abs, symbolic=False):	1✔
141	"""
142	Generate the inner product of `x := self` and `y := other`.
143
144	:param diff_func: How to treat the difference between the binary reps of the basis blades.
145	if :code:`abs`, compute the symmetric inner product. When :code:`lambda x: -x` this
146	function generates left-contraction, and when :code:`lambda x: x`, right-contraction.
147	:return: tuple of keys in binary representation and a lambda function.
148	"""
149	algebra = x.algebra	1✔
150	filter_func = lambda tt: tt.key_out == diff_func(tt.keys_in[0] - tt.keys_in[1])	1✔
151	keys_out, func = codegen_product(x, y, filter_func=filter_func, name_base='ip')	1✔
152	if symbolic:	1✔
153	return x.fromkeysvalues(algebra, keys_out, func(x.values(), y.values()))	×
154
155	return keys_out, func	1✔
156
157
158	def codegen_lc(x, y):	1✔
159	"""
160	Generate the left-contraction of `x := self` and `y := other`.
161
162	:return: tuple of keys in binary representation and a lambda function.
163	"""
164	return codegen_ip(x, y, diff_func=lambda x: -x)	1✔
165
166
167	def codegen_rc(x, y):	1✔
168	"""
169	Generate the right-contraction of `x := self` and `y := other`.
170
171	:return: tuple of keys in binary representation and a lambda function.
172	"""
173	return codegen_ip(x, y, diff_func=lambda x: x)	1✔
174
175
176	def codegen_sp(x, y):	1✔
177	"""
178	Generate the scalar product of `x := self` and `y := other`.
179
180	:return: tuple of keys in binary representation and a lambda function.
181	"""
182	return codegen_ip(x, y, diff_func=lambda x: 0)	1✔
183
184
185	def codegen_proj(x, y):	1✔
186	"""
187	Generate the projection of `x := self` onto `y := other`: :math:`(x \cdot y) / y`.
188
189	:return: tuple of keys in binary representation and a lambda function.
190	"""
191	filter_func = lambda tt: tt.key_out == abs(tt.keys_in[0] - tt.keys_in[1]) ^ tt.keys_in[2]	1✔
192	res = codegen_product(x, y, ~y, name_base='gp', filter_func=filter_func, asdict=True)	1✔
193	return _func_builder(res, x, y, name_base='proj')	1✔
194
195
196	def codegen_op(x, y, symbolic=False):	1✔
197	"""
198	Generate the outer product of `x := self` and `y := other`: `x.op(y) = x ^ y`.
199
200	:x: MultiVector
201	:y: MultiVector
202	:return: dictionary with integer keys indicating the corresponding basis blade in binary convention,
203	and values which are a 3-tuple of indices in `x`, indices in `y`, and a lambda function.
204	"""
205	algebra = x.algebra	1✔
206	filter_func = lambda tt: tt.key_out == sum(tt.keys_in)	1✔
207	sign_func = lambda pair: (-1)**algebra.swaps[pair]	1✔
208	keys_out, func = codegen_product(x, y, filter_func=filter_func, sign_func=sign_func, name_base='op')	1✔
209	if symbolic:	1✔
210	return x.fromkeysvalues(algebra, keys_out, func(x.values(), y.values()))	×
211
212	return keys_out, func	1✔
213
214
215	def codegen_rp(x, y):	1✔
216	"""
217	Generate the regressive product of :code:`x` and :code:`y`:,
218	:math:`x \\vee y`.
219
220	:param x:
221	:param y:
222	:return: tuple of keys in binary representation and a lambda function.
223	"""
224	algebra = x.algebra	1✔
225	keyout_func = lambda tot, key_in: len(algebra) - 1 - (key_in ^ tot)	1✔
226	filter_func = lambda tt: len(algebra) - 1 == sum(tt.keys_in) - tt.key_out	1✔
227	# Sign is composed of dualization of each blade, exterior product, and undual.
228	sign_func = lambda pair: (	1✔
229	algebra.signs[pair[0], len(algebra) - 1 - pair[0]] *
230	algebra.signs[pair[1], len(algebra) - 1 - pair[1]] *
231	(-1)*algebra.swaps[len(algebra) - 1 - pair[0], len(algebra) - 1 - pair[1]]
232	algebra.signs[len(algebra) - 1 - (pair[0] ^ pair[1]), pair[0] ^ pair[1]]
233	)
234
235	return codegen_product(	1✔
236	x, y,
237	filter_func=filter_func,
238	keyout_func=keyout_func,
239	sign_func=sign_func,
240	name_base='rp'
241	)
242
243
244	Fraction = namedtuple('Fraction', ['numer', 'denom'])	1✔
245	Fraction.__doc__ = """	1✔
246	Tuple representing a fraction.
247	"""
248
249	def codegen_inv(x, symbolic=False):	1✔
250	"""
251	Generate code for the inverse of :code:`x`.
252	Currently, this always uses the Shirokov inverse, which is works in any algebra,
253	but it can be expensive to compute.
254	In the future this should be extended to use dedicated solutions for known cases.
255	"""
256	k = 2 ** ((x.algebra.d + 1) // 2)	1✔
257	x_i = x	1✔
258	if x.grades == (0,):	1✔
259	adj_x = 1	1✔
260	else:
261	for i in range(1, k + 1):	1✔
262	# Sympify ratio to keep the ratios exact and avoid floating point errors.
263	c_i = (sympify(k) / i) * x_i[0] if x_i[0] else x_i[0]	1✔
264	adj_x = (x_i - c_i)	1✔
265	if x.algebra.simplify:	1✔
266	keys, values = zip(*((k, simp_expr) for k, expr in adj_x.items() if (simp_expr := simplify(expr))))	1✔
267	adj_x = adj_x.fromkeysvalues(adj_x.algebra, keys=keys, values=values)	1✔
268	x_i = x * adj_x	1✔
269	if x_i:	1✔
270	if x.algebra.simplify:	1✔
271	keys, values = zip(*((k, simp_expr) for k, expr in x_i.items() if (simp_expr := simplify(expr))))	1✔
272	x_i = x_i.fromkeysvalues(x_i.algebra, keys=keys, values=values)	1✔
273
274	if x_i.grades == (0,):	1✔
275	break	1✔
276	else:
277	break	1✔
278
279	if symbolic:	1✔
280	return Fraction(adj_x, x_i[0])	1✔
281	xinv = x.algebra.multivector({k: v / x_i[0] for k, v in adj_x.items()})	1✔
282	return _lambdify_unary(x, xinv)	1✔
283
284
285	def codegen_div(x, y):	1✔
286	"""
287	Generate code for :math:`x y^{-1}`.
288	"""
289	adjy, normsqy = codegen_inv(y, symbolic=True)	1✔
290	if not normsqy:	1✔
291	raise ZeroDivisionError	1✔
292	x_adjy = x * adjy	1✔
293	if x.algebra.simplify:	1✔
294	xdivy = x.algebra.multivector({k: simplify(v / normsqy) for k, v in x_adjy.items()})	1✔
295	else:
296	xdivy = x.algebra.multivector({k: v / normsqy for k, v in x_adjy.items()})	×
297	return _lambdify_binary(x, y, xdivy)	1✔
298
299	def codegen_normsq(x):	1✔
300	if x.algebra.simplify:	1✔
301	res = codegen_product(x, ~x, name_base='gp', asdict=True, sympy=True)	1✔
302	res = {k: str(simp_expr) for k, expr in res.items() if (simp_expr := simplify(expr))}	1✔
303	else:
304	res = codegen_product(x, ~x, name_base='gp', asdict=True)	×
305	return _func_builder(res, x, name_base="normsq")	1✔
306
307
308	def codegen_outerexp(x, asterms=False):	1✔
309	if len(x.grades) != 1:	1✔
310	warnings.warn('Outer exponential might not converge for mixed-grade multivectors.', RuntimeWarning)	×
311	k = x.algebra.d	1✔
312
313	Ws = [x.algebra.scalar([1]), x]	1✔
314	j = 2	1✔
315	while j <= k:	1✔
316	Wj = Ws[-1] ^ x	1✔
317	# Dividing like this avoids floating point numbers, which is excellent.
318	Wj._values = tuple(v / j for v in Wj._values)	1✔
319	if Wj:	1✔
320	Ws.append(Wj)	1✔
321	j += 1	1✔
322	else:
323	break	1✔
324
325	if asterms:	1✔
326	return Ws	×
327	L = reduce(operator.add, Ws)	1✔
328	return _func_builder(dict(L.items()), x, name_base='outerexp')	1✔
329
330
331	def _lambdify_binary(x, y, x_bin_y):	1✔
332	xy_symbols = [list(x.values()), list(y.values())]	1✔
333	func = lambdify(xy_symbols, list(x_bin_y.values()), cse=x.algebra.cse)	1✔
334	return tuple(x_bin_y.keys()), njit(func) if x.algebra.numba else func	1✔
335
336
337	def _lambdify_unary(x, x_unary):	1✔
338	func = lambdify([list(x.values())], list(x_unary.values()), cse=x.algebra.cse)	1✔
339	return tuple(x_unary.keys()), njit(func) if x.algebra.numba else func	1✔
340
341
342	def _lambdify_mv(free_symbols, mv):	1✔
343	# TODO: Numba wants a tuple in the line below, but simpy only produces a
344	# list as output if this is a list, not a tuple. See if we can solve this.
345	func = lambdify(free_symbols, list(mv.values()), cse=mv.algebra.cse)	1✔
346	return tuple(mv.keys()), njit(func) if mv.algebra.numba else func	1✔
347
348	def _func_builder(res_vals: dict, *mvs, name_base: str):	1✔
349	"""
350	Build a Python functions for the product between given multivectors.
351
352	:param res_vals: Dict to be converted into a function. The keys correspond to the basis blades in binary,
353	while the values are strings to be converted into source code.
354	:param mvs: all the multivectors that the resulting function is a product of.
355	:param name_base: Base for the name of the function. E.g. 'gp'.
356	:return: tuple of output keys of the callable, and the callable.
357	"""
358	args = [f'arg{i}' for i in range(1, len(mvs) + 1)]	1✔
359	keys_str_per_arg = ["_".join(str(k) for k in mv.keys()) for mv in mvs]	1✔
360	func_name = f'{name_base}_{"_x_".join(keys_str_per_arg)}'	1✔
361	header = f'def {func_name}({", ".join(args)}):'	1✔
362	if res_vals:	1✔
363	body = "\n".join(f' {",".join(v.name for v in mv.values())}, = {arg}' for mv, arg in zip(mvs, args))	1✔
364	return_val = f' return ({", ".join(v if isinstance(v, str) else str(v) for v in res_vals.values())},)'	1✔
365	else:
366	body = ''	1✔
367	return_val = f' return tuple()'	1✔
368	func_source = f'{header}\n{body}\n{return_val}'	1✔
369
370	# Dynamically build a function
371	func_locals = {}	1✔
372	c = compile(func_source, func_name, 'exec')	1✔
373	exec(c, {}, func_locals)	1✔
374
375	# Add the generated code to linecache such that it is inspect-safe.
376	linecache.cache[func_name] = (len(func_source), None, func_source.splitlines(True), func_name)	1✔
377	func = func_locals[func_name]	1✔
378	return tuple(res_vals.keys()), njit(func) if mvs[0].algebra.numba else func	1✔

tBuLi / kingdon / 3356190533

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous