Source code for pybs.series.operations

from fractions import Fraction
from math import factorial

from pybs.utils import memoized2 as memoized, LinearCombination
from pybs.unordered_tree import leaf
from pybs.combinations import \
    empty_tree, \
    split, \
    subtrees, \
    antipode_ck, \
    tree_commutator as tree_commutator
from pybs.series import \
    BseriesRule, \
    VectorfieldRule, \
    ForestRule, \
    tree_tuples_of_order
# TODO: Include Sphinx in installation?


[docs]def hf_composition(a): r"""The composition :math:`a b`, where :math:`b` is the B-series representing the vector field corresponding to the exact solution. That is :math:`b = \delta_{\circ}`. """ # TODO: include some trees in the sphinx docs. if a(empty_tree) != 1: raise ValueError( 'Composition can only be performed on consistent B-series.') def new_rule(tree): if tree == empty_tree: return 0 else: result = 1 for subtree, multiplicity in tree.items(): result *= a(subtree) ** multiplicity return result return BseriesRule(new_rule)
[docs]def lie_derivative(c, b, truncate=False): """The Lie-derivative of ``c`` with respect to ``b`` as a new :class:`series.Bseries.BseriesRule`. """ if b(empty_tree) != 0: raise ValueError( 'The second argument does not satisfy b(ButcherEmptyTree()) == 0.') @memoized def new_rule(tree): result = 0 if tree == empty_tree: return result pairs = split(tree, truncate) for pair, multiplicity in pairs.items(): result += multiplicity * b(pair[0]) * c(pair[1]) return result return VectorfieldRule(new_rule)
[docs]def modified_equation(a, quadratic_vectorfield=False): """The modified equation. Mostely equivalent to :func:`log`.""" if a(empty_tree) != 1 or a(leaf) != 1: # TODO: Check last condition. raise ValueError( 'Can not calculate the modified equation for this BseriesRule.') @memoized def new_rule(tree): if tree == empty_tree: return 0 result = a(tree) c = new_rule # This is a BseriesRule. Caution: Recursive! for j in range(2, tree.order() + 1): c = lie_derivative(c, new_rule, True) # TODO: Is this memoized? result -= Fraction(c(tree), factorial(j)) return result result = VectorfieldRule(new_rule) if quadratic_vectorfield: result = remove_non_binary(result) return result
[docs]def log(a, quadratic_vectorfield=False): """The logarithm. If ``a`` is the B-series rule for a numerical method, it returns the rule for the modified equation. .. note:: Much slower than :func:`modified_equation`. """ if a(empty_tree) != 1: raise ValueError( 'Can not calculate the logarithm for this rule.') @memoized def new_rule(tree): if tree == empty_tree: return 0 a_2 = remove_empty_tree(a) result = a_2(tree) b = a_2 for n in range(2, tree.order() + 1): b = composition(b, a_2) result += ((-1)**(n+1)) * Fraction(b(tree), n) return result result = VectorfieldRule(new_rule) if quadratic_vectorfield: result = remove_non_binary(result) return result
[docs]def exp(a, quadratic_vectorfield=False): """The exponential. If ``a`` is the rule for the B-serie of some modified equation, it returns the B-series rule for the numerical method. """ if a(empty_tree) != 0: raise ValueError( 'Can not calculate the exponential for this rule.') @memoized def new_rule(tree): if tree == empty_tree: return 1 result = a(tree) b = a for n in range(2, tree.order() + 1): b = composition(b, a) result += Fraction(b(tree), factorial(n)) return result result = BseriesRule(new_rule) if quadratic_vectorfield: result = remove_non_binary(result) return result
[docs]def remove_non_binary(a): """Sets the value at all non-binary trees to zero. Used for quadratic vector fields.""" base_rule = a._call def new_rule(tree): if tree == empty_tree or tree.is_binary(): return base_rule(tree) else: return 0 return type(a)(new_rule)
[docs]def remove_empty_tree(a): """Returns :math:`a - I`. Used by :func:`log`.""" base_rule = a._call et = empty_tree def new_rule(tree): if tree == et: return 0 else: return base_rule(tree) return type(a)(new_rule)
[docs]def composition_ssa(a, b): """Same as :func:`composition`, except that it halves the stepsize afterwards. Equivalent to ``stepsize_adjustment(composition(a,b),Fraction(1, 2))``. """ if a(empty_tree) != 1: raise ValueError( 'Composition can only be performed on consistent B-series.') @memoized def new_rule(tree): sub_trees = subtrees(tree) result = 0 for pair, multiplicity in sub_trees.items(): result += a(pair[0]) * b(pair[1]) * multiplicity return Fraction(result, 2**tree.order()) return BseriesRule(new_rule) # TODO: account for ForestRule.
[docs]def stepsize_adjustment(a, A): """Corresponds to letting h -> A h.""" base_rule = a._call def new_rule(tree): return A**tree.order() * base_rule(tree) return BseriesRule(new_rule)
[docs]def composition(a, b): r"""Composition of methods, b after a. Return the composition :math:`a \circ b`. The returned object is a :class:`BseriesRule` if :math:`a(\emptyset) = 1` and a :class:`ForestRule` otherwise. """ @memoized def new_rule(arg): # arg is tree or forest. result = 0 for pair, multiplicity in subtrees(arg).items(): result += a(pair[0]) * b(pair[1]) * multiplicity return result if a(empty_tree) != 1: return ForestRule(new_rule) else: return BseriesRule(new_rule)
[docs]def inverse(a): r"""Return the inverse of *a* in the Butcher group. The returned BseriesRule is calculated as :math:`\left(a \circ S\right)(\tau)`, where :math:`S` denotes the antipode. """ # TODO: Test that a is of the right kind. @memoized def new_rule(tree): return a(antipode_ck(tree)) return BseriesRule(new_rule)
[docs]def conjugate(a, c): """The conjugate of ``a`` with change of coordinates ``c``. Calculated with compositions and inverse.""" return composition(inverse(c), composition(a, c))
[docs]def conjugate_by_commutator(a, c): """The conjugate of ``a`` with change of coordinates ``c``. Calculated using the commutator. """ def new_rule(tree): if tree == empty_tree: return 0 tmp = a result = 0 for n in range(tree.order() + 1): result += Fraction((-1)**n, factorial(n)) * tmp(tree) tmp = series_commutator(c, tmp) return result return BseriesRule(new_rule)
[docs]def adjoint(a): """The adjoint is the inverse with reversed time step. Returns a BseriesRule. """ def new_rule(tree): return (-1)**tree.order() * inverse(a)(tree) return BseriesRule(new_rule) # TODO: Adjoint of a modified equation p. 324 HLW.
[docs]def series_commutator(a, b): """Corresponds to tree commutator, just for series.""" # TODO: TEST ME! def new_rule(tree): order = tree.order() if order in new_rule.orders: return new_rule.storage[tree] * tree.symmetry() # TODO: Move the correction by symmetry to initialisation. else: result = LinearCombination() for tree1, tree2 in tree_tuples_of_order(order): result += \ Fraction(a(tree1) * b(tree2), tree1.symmetry() * tree2.symmetry()) * \ tree_commutator(tree1, tree2) new_rule.orders.add(order) new_rule.storage += result return new_rule.storage[tree] * tree.symmetry() new_rule.storage = LinearCombination() new_rule.orders = set((0,)) # the coefficient of the empty tree is always 0. return BseriesRule(new_rule)