Skip to content

Commit

Permalink
Updated docstrings, improved scalar division, more robust against emp…
Browse files Browse the repository at this point in the history
…ty multivectors.
  • Loading branch information
tBuLi committed Jun 14, 2023
1 parent a1eb5fa commit 9e78112
Show file tree
Hide file tree
Showing 4 changed files with 55 additions and 48 deletions.
30 changes: 15 additions & 15 deletions docs/usage.rst
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
Basic Usage
===========

The most important object in all of :mod:`kingdon` is :code:`Algebra`::
The most important object in all of :code:`kingdon` is :code:`Algebra`::

from kingdon import Algebra

Expand Down Expand Up @@ -46,17 +46,17 @@ For example, let us create two symbolic vectors :code:`u` and :code:`v` in this
>>> v
(v1) * e1 + (v2) * e2
The return type of :code:`Algebra.multivector` is an instance of :class:`~kingdon.MultiVector`.
The return type of :code:`Algebra.multivector` is an instance of :class:`~kingdon.multivector.MultiVector`.

.. note::
:code:`kingdon` offers convenience methods for common types of multivectors, such as the vectors above.
For example, the vectors above can also be created using :code:`u = alg.vector(name='u')`.
Moreover, a scalar is created by :meth:`~kingdon.Algebra.scalar`, a bivector by :meth:`~kingdon.Algebra.bivector`,
a pseudoscalar by :meth:`~kingdon.Algebra.pseudoscalar`, and so on.
Moreover, a scalar is created by :meth:`~kingdon.algebra.Algebra.scalar`, a bivector by :meth:`~kingdon.algebra.Algebra.bivector`,
a pseudoscalar by :meth:`~kingdon.algebra.Algebra.pseudoscalar`, and so on.
However, all of these merely add the corresponding :code:`grades` argument to your input and
then call :code:`alg.multivector`, so :code:`alg.multivector` is what we need to understand.

:class:`~kingdon.MultiVector`'s support common math operators:
:class:`~kingdon.multivector.MultiVector`'s support common math operators:

.. code-block::
Expand Down Expand Up @@ -87,7 +87,7 @@ line :code:`v` in the line :code:`u` by using conjugation:
we see that the result is again a vector, as it should be.

These examples should show that the symbolic multivectors of :mod:`kingdon`
These examples should show that the symbolic multivectors of :code:`kingdon`
make it easy to do symbolic computations. Moreover, we can also use :mod:`sympy` expressions
as values for the multivector:

Expand Down Expand Up @@ -132,8 +132,8 @@ For example, to repeat some of the examples above with numerical values, we coul
(0.1541) + (0.0886) * e12
A big performance bottleneck that we suffer from in Python, is that arrays over objects are very slow.
So while we could make a numpy array filled with :code:`~kingdon.MultiVector`'s, this would tank our performance.
:mod:`kingdon` gets around this problem by instead accepting numpy arrays as input. So to make a collection of
So while we could make a numpy array filled with :code:`~kingdon.multivector.MultiVector`'s, this would tank our performance.
:code:`kingdon` gets around this problem by instead accepting numpy arrays as input. So to make a collection of
3 lines, we do

.. code-block::
Expand All @@ -145,7 +145,7 @@ So while we could make a numpy array filled with :code:`~kingdon.MultiVector`'s,
([0.82499172 0.71181276 0.98052928]) * e1 + ([0.53395072 0.07312351 0.42464341]) * e2
what is important here is that the first dimension of the array has to have the expected length: 2 for a vector.
All other dimensions are not used by :mod:`kingdon`. Now we can reflect this multivector in the :code:`e1` line:
All other dimensions are not used by :code:`kingdon`. Now we can reflect this multivector in the :code:`e1` line:

.. code-block::
Expand All @@ -159,7 +159,7 @@ and with only minor performance penalties.
Operators
---------

Instances of :mod:`~kingdon.MultiVector` overload all common Geometric Algebra operators.
Instances of :mod:`~kingdon.multivector.MultiVector` overload all common Geometric Algebra operators.
Below is an overview:

.. list-table:: Operators
Expand Down Expand Up @@ -243,7 +243,7 @@ when codegen is performed for the given input, which is why this isn't the defau
Graphing using :code:`ganja.js`
-------------------------------

:mod:`kingdon` supports the :code:`ganja.js` graphing syntax. For those already familiar with
:code:`kingdon` supports the :code:`ganja.js` graphing syntax. For those already familiar with
:code:`ganja.js`, the API will feel very similar:

.. code-block::
Expand All @@ -255,7 +255,7 @@ elements to graph, whereas keyword arguments are passed to :code:`ganja.js` as o
Hence, the example above will graph the line :code:`u` with :code:`lineWidth = 3`,
and will attach the label "u" to it, and all of this will be red.
Identical to :code:`ganja.js`, valid inputs to :code:`alg.graph` are (lists of) instances
of :class:`~kingdon.MultiVector`, strings, and hexadecimal numbers to indicate colors.
of :class:`~kingdon.multivector.MultiVector`, strings, and hexadecimal numbers to indicate colors.
These strings can be simple labels, or valid SVG syntax.

.. note::
Expand All @@ -265,10 +265,10 @@ These strings can be simple labels, or valid SVG syntax.

Performance
-----------
Because :mod:`kingdon` attempts to symbolically optimize expressions
Because :code:`kingdon` attempts to symbolically optimize expressions
using :mod:`sympy` the first time they are called, the first call to any operation is always slow,
whereas subsequent calls have extremely good performance.
This is because :mod:`kingdon` first leverages the sparseness of the input,
This is because :code:`kingdon` first leverages the sparseness of the input,
*and* subsequently uses symbolic optimization to eliminate any terms that are always zero
regardless of the input.
For example, the product :math:`\mathbf{e}_{1} \mathbf{e}_{12}` of the vector :math:`\mathbf{e}_1`
Expand All @@ -277,7 +277,7 @@ and the bivector :math:`\mathbf{e}_{12}` in :math:`\mathbb{R}_{2+p',q,r}` always
In :code:`kingdon`, it will also be equally fast to compute this product in all of these algebras,
regardless of the total dimension.

Because the precomputation can get expensive, :mod:`kingdon` predefines all the popular algebras
Because the precomputation can get expensive, :code:`kingdon` predefines all the popular algebras
of :math:`d = p+q+r < 6`.
For example, a precomputed version of 3DPGA can be imported as

Expand Down
2 changes: 1 addition & 1 deletion docs/workings.rst
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
==============
Inner Workings
==============
This chapter explains how :mod:`kingdon` works internally.
This chapter explains how :code:`kingdon` works internally.
67 changes: 36 additions & 31 deletions kingdon/codegen.py
Original file line number Diff line number Diff line change
Expand Up @@ -82,13 +82,13 @@ class TermTuple(NamedTuple):
: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, tuple of :class:`~sympy.core.Symbol`.
:param values_in: Input values. Typically, tuple of :class:`sympy.core.symbol.Symbol`.
:param termstr: The string representation of this monomial, e.g. '-x*y'.
"""
key_out: int
keys_in: Tuple[int]
sign: int
values_in: Tuple[Symbol]
values_in: Tuple["sympy.core.symbol.Symbol"]
termstr: str


Expand Down Expand Up @@ -310,32 +310,35 @@ def codegen_rp(x, y):
def codegen_inv(y, x=None, symbolic=False):
alg = y.algebra
# As preprocessing we invert y*~y since y^{-1} = ~y / (y*~y)
ynormsq = y.normsq()

if ynormsq.grades == tuple():
raise ZeroDivisionError
elif ynormsq.grades == (0,):
adj_y, denom_y = ~y if x is None else x * ~y, ynormsq[0]
if len(y) == 1 and y.grades == (0,):
adj_y, denom_y = x or 1, y.e
else:
# Make a mv with the same components as ynormsq, and invert that instead.
# Although this requires some more bookkeeping, it is much more performant.
z = alg.multivector(name='z', keys=ynormsq.keys())
adj_z, denom_z = codegen_shirokov_inv(z, symbolic=True)

# Same trick, make a mv that mimicks adj_z and multiply it with ~y.
A_z = alg.multivector(name='A_z', keys=adj_z.keys())
# Faster to this `manually` instead of with ~a * A_z
res = codegen_product(~y, A_z, symbolic=True)
A_y = res if x is None else x * res

# Replace all the dummy A_z symbols by the expressions in adj_z.
subs_dict = dict(zip(A_z.values(), adj_z.values()))
adj_y = A_y.subs(subs_dict)
# Replace all the dummy b symbols by the expressions in anormsq to
# recover the actual adjoint of a and corresponding denominator.
subs_dict = dict(zip(z.values(), ynormsq.values()))
denom_y = denom_z.subs(subs_dict)
adj_y = adj_y.subs(subs_dict)
ynormsq = y.normsq()

if ynormsq.grades == tuple():
raise ZeroDivisionError
elif ynormsq.grades == (0,):
adj_y, denom_y = ~y if x is None else x * ~y, ynormsq.e
else:
# Make a mv with the same components as ynormsq, and invert that instead.
# Although this requires some more bookkeeping, it is much more performant.
z = alg.multivector(name='z', keys=ynormsq.keys())
adj_z, denom_z = codegen_shirokov_inv(z, symbolic=True)

# Same trick, make a mv that mimicks adj_z and multiply it with ~y.
A_z = alg.multivector(name='A_z', keys=adj_z.keys())
# Faster to this `manually` instead of with ~a * A_z
res = codegen_product(~y, A_z, symbolic=True)
A_y = res if x is None else x * res

# Replace all the dummy A_z symbols by the expressions in adj_z.
subs_dict = dict(zip(A_z.values(), adj_z.values()))
adj_y = A_y.subs(subs_dict)
# Replace all the dummy b symbols by the expressions in anormsq to
# recover the actual adjoint of a and corresponding denominator.
subs_dict = dict(zip(z.values(), ynormsq.values()))
denom_y = denom_z.subs(subs_dict)
adj_y = adj_y.subs(subs_dict)

if symbolic:
return Fraction(adj_y, denom_y)
Expand Down Expand Up @@ -476,8 +479,8 @@ def _lambdify_mv(free_symbols, mv):
def do_codegen(codegen, *mvs) -> CodegenOutput:
"""
:param codegen: callable that performs codegen for the given :code:`mvs`. This can be any callable
that returns either a :class:`MultiVector`, a dictionary, or an instance of :class:`CodegenOutput`.
:param *mvs: Any remaining positional arguments are taken to be symbolic :class:`Multivector`'s.
that returns either a :class:`~kingdon.multivector.MultiVector`, a dictionary, or an instance of :class:`CodegenOutput`.
:param mvs: Any remaining positional arguments are taken to be symbolic :class:`~kingdon.multivector.MultiVector`'s.
:return: Instance of :class:`CodegenOutput`.
"""
res = codegen(*mvs)
Expand All @@ -495,7 +498,7 @@ def do_codegen(codegen, *mvs) -> CodegenOutput:
def lambdify(args: dict, exprs: tuple, funcname: str, dependencies: tuple = None, printer=NumPyPrinter, dummify=False, cse=False):
"""
Function that turns symbolic expressions into Python functions. Heavily inspired by
:mod:`sympy`'s function by the same name, but adapted for the needs of :mod:`kingdon`.
:mod:`sympy`'s function by the same name, but adapted for the needs of :code:`kingdon`.
Particularly, this version gives us more control over the names of the function and its
arguments, and is more performant, particularly when the given expressions are strings.
Expand Down Expand Up @@ -532,7 +535,7 @@ def cp(A, B):
included in the CSE process.
:param cse: If :code:`True` (default), CSE is applied to the expressions and dependencies.
This typically greatly improves performance and reduces numba's initialization time.
return: Function that represents that can be used to calculate the values of exprs.
:return: Function that represents that can be used to calculate the values of exprs.
"""
if printer is NumPyPrinter:
printer = NumPyPrinter(
Expand Down Expand Up @@ -561,6 +564,8 @@ def cp(A, B):
else:
cses, _exprs = list(zip(flatten(lhsides), flatten(rhsides))), exprs

if not any(_exprs):
_exprs = tuple('0' for expr in _exprs)
funcstr = funcprinter.doprint(funcname, iterable_args, names, _exprs, cses=cses)

# Provide lambda expression with builtins, and compatible implementation of range
Expand Down
4 changes: 3 additions & 1 deletion kingdon/operator_dict.py
Original file line number Diff line number Diff line change
Expand Up @@ -93,7 +93,9 @@ def _call_binary(self, mv1, mv2):

keys_out, func, numba_func = self[mv1.keys(), mv2.keys()]
issymbolic = (mv1.issymbolic or mv2.issymbolic)
if issymbolic or not mv1.algebra.numba:
if not keys_out:
values_out = tuple()
elif issymbolic or not mv1.algebra.numba:
values_out = func(mv1.values(), mv2.values())
else:
values_out = numba_func(mv1.values(), mv2.values())
Expand Down

0 comments on commit 9e78112

Please sign in to comment.