Refactor Components into parent & element

[mkoeppe]

Currently every Components object has lots of metadata attributes in addition to the actual data dictionary in ._comp.

If one has many different Components objects with the same symmetry metadata, we can reduce storage space as follows.

We create new classes CompParent, CompParentWithSym, ..., which
store the symmetry attributes and become UniqueRepresentation. We make Components objects elements of these parents.

The parents will also have index iterator methods.

Data associated with symmetries and dimensions, computed currently each time for each CompWithSym object in methods such as __init__, __add__, trace, ... can be precomputed and cached in the parent using @cached_method in the parent class).

This will make the code in #30229 (subspaces of tensor with symmetries) more elegant because it no longer needs a dummy Components object to represent the symmetry but rather a CompParent object.

See also:

• Action of a sympy TensorSymmetry #32029 Action of a sympy TensorSymmetry

Depends on #31904
Depends on #29619

CC: @mjungmath @egourgoulhon @tscrim @honglizhaobob

Component: linear algebra

Author: Michael Jung, Matthias Koeppe

Branch/Commit: public/30307 @ e5a81f2

Issue created by migration from https://trac.sagemath.org/ticket/30307

[mkoeppe]

comment:2

Help with implementing this would be very welcome!

[mjungmath]

comment:3

That's a nice idea, +1!

Data associated with symmetries, computed currently each time for each CompWithSym object in methods such as __init__, __add__, trace, ... can then be precomputed and cached in the parent (for example using @cached_method in the parent class).

I think it would be a good idea to cache results from trace, contract, ... as well.

The most annoying thing with this ticket will most likely be the docstring that has to be adapted...is there a good way to use search&replace?

[mjungmath]

comment:4

Currently, the index generators and manipulators take (mostly) lists and can therefore not be cached. To use the caching most efficiently, I suggest we switch entirely to tuples.

[mjungmath]

Branch: public/30307

[sagetrac-git]

Commit: 943f9ea

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

fec4b06	fix sage -b after 30622
dff846c	go through make/build/install
943f9ea	Trac #30307: CompParent and CompParentWithSym

[mjungmath]

comment:7

Before I go on with this ticket, could you please take a look whether this meets your rough idea?

Is anyone willing to work on the docstrings...? Perhaps there's an efficient way to do this?

[sagetrac-git]

Changed commit from 943f9ea to a2514b6

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

a2514b6

Trac #30307: no symmetries -> CompParent without symmetries

[mkoeppe]

comment:9

Replying to @mjungmath:

Before I go on with this ticket, could you please take a look whether this meets your rough idea?

Yes, this is going in the direction that I had in mind.

Some more of the normalization happening in CompParentWithSym.__init__ should probably be moved to __classcall_private__

[mjungmath]

comment:10

What would be a proper category for the parent?

[mjungmath]

comment:11

Replying to @mkoeppe:

Some more of the normalization happening in CompParentWithSym.__init__ should probably be moved to __classcall_private__

For example?

[mkoeppe]

comment:12

For example the order of the component indices does not matter

[egourgoulhon]

comment:13

Replying to @mjungmath:

What would be a proper category for the parent?

Take the following with a grain of salt (this is just a rough/naive thought): It seems to me that the current proposal is a kind of abuse of Sage's parent/element scheme, for CompParent does not correspond to a "genuine" mathematical object (hence maybe your question...). In particular, it does not know about the ring on which the components are based. I do not deny that a reorganization of Components would be a good thing, but maybe at the class level, not at the parent/element level, the issue here being mostly the storage of metadata.

[mkoeppe]

comment:14

Just use the category of sets. I don't think this is an abuse.

Think of an instance of CompParent as the set of all possible Components instances that have the specified symmetry.

By using the parent/element framework, you will get coercion for free - so elements with a coarser symmetry will be in a finer parent.

[mjungmath]

comment:15

Replying to @mkoeppe:

For example the order of the component indices does not matter

Right!

Replying to @mkoeppe:

Just use the category of sets. I don't think this is an abuse.

Think of an instance of CompParent as the set of all possible Components instances that have the specified symmetry.

By using the parent/element framework, you will get coercion for free - so elements with a coarser symmetry will be in a finer parent.

For me, it doesn't feel like an abuse either. Besides, I am sure you can make that notion of compontents mathematically rigorous. But I am not sure whether this becomes a well-defined set then... What about the category of objects?

[mjungmath]

comment:16

Alright, I make progress here. I changed only the backend such that most doctests should still pass, and the module can be used exactly as before. The elementary examples already passed.

I'll push my branch tomorrow.

I am looking forward to some benchmarks as soon as this branch is ready. :)

[tscrim]

comment:17

I can understand why this might seem like an abuse because it is a set of objects, but by extension, all of the combinatorial objects would be an abuse as well. So even though we will not be putting an extra mathematical structure on the components, this is still a valid use of the framework because there is a set of objects (the parent) and the individual objects (the elements).

Something to consider, however, is a potential speed penalty for the extra levels of initialization. This can be somewhat mitigated by using Cython, but performance regression testing is warranted here.

[mjungmath]

comment:18

Replying to @tscrim:

This can be somewhat mitigated by using Cython, but performance regression testing is warranted here.

How so?

[tscrim]

comment:19

Replying to @mjungmath:

Replying to @tscrim:

This can be somewhat mitigated by using Cython, but performance regression testing is warranted here.

How so?

Because Cython is faster than Python, including potential benefits from direct C-level function calls.

[mjungmath]

comment:20

Yes, that's clear. I meant, how to implement? Make CompParent and Components both simply extension types?

Would cdpefing the symmetry related functions cause a speedup btw?

Since this involves only integers, i.e. struct types, one could even try to Cythonize these completely.

[sagetrac-git]

Changed commit from a2514b6 to fc22d5d

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

fc22d5d

Trac 30307: Cythonize, restructuring, factory method to resemble old behavior

[egourgoulhon]

comment:22

Replying to @tscrim:

I can understand why this might seem like an abuse because it is a set of objects, but by extension, all of the combinatorial objects would be an abuse as well. So even though we will not be putting an extra mathematical structure on the components, this is still a valid use of the framework because there is a set of objects (the parent) and the individual objects (the elements).

Thank you Matthias, Michael and Travis for your replies. I'm convinced now that parent/element is a good approach here (I was bothered by the lack of algebraic structure of the parent), especially for symmetry-based coercions.

Something to consider, however, is a potential speed penalty for the extra levels of initialization. This can be somewhat mitigated by using Cython, but performance regression testing is warranted here.

Certainly!

---

New commits:

fc22d5d

Trac 30307: Cythonize, restructuring, factory method to resemble old behavior

[mkoeppe]

Changed dependencies from #31904 to #31904, #29619

[sagetrac-git]

Changed commit from 3070443 to edd1572

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

0ca670b	sage.tensor.modules: Update doctest output for trivial output changes
02d6917	VectorFieldModule.tensor_from_comp: Use CompParent*
c68734c	sage.manifolds: Update doctest output for trivial output changes
8ed01a7	CompParent, Components_dict: Fix up module action
edd1572	TensorField, ComponentsWithSym_dict: Remove comparison of _sym, _antisym with empty list, normalize to tuples

[mkoeppe]

comment:91

sage.tensor.modules passes all tests except for _test_not_implemented_methods (to be fixed by #29619)

There are still various errors in sage.manifolds.
For example, in src/sage/manifolds/differentiable/multivectorfield.py:

sage -t --random-seed=0 src/sage/manifolds/differentiable/multivectorfield.py
**********************************************************************
File "src/sage/manifolds/differentiable/multivectorfield.py", line 117, in sage.manifolds.differentiable.multivectorfield.MultivectorField
Failed example:
    a.display(eU)
Expected:
    a = (x*y^2 + 2*x) d/dx/\d/dy
Got:
    a = (x*y^2 + 2*x) d/dx*d/dy + (-x*y^2 - 2*x) d/dy*d/dx

Help in spotting where these errors are coming from would be welcome

[mjungmath]

comment:92

For some reason, the restrict method doesn't work properly:

sage: M = Manifold(2, 'M')
sage: U = M.open_subset('U') ; V = M.open_subset('V')
sage: M.declare_union(U,V)   # M is the union of U and V
sage: c_xy.<x,y> = U.chart() ; c_uv.<u,v> = V.chart()
sage: xy_to_uv = c_xy.transition_map(c_uv, (x+y, x-y),
....:                         intersection_name='W',
....:                         restrictions1= x>0, restrictions2= u+v>0)
sage: uv_to_xy = xy_to_uv.inverse()
sage: W = U.intersection(V)
sage: eU = c_xy.frame() ; eV = c_uv.frame()
sage: a = M.multivector_field(2, name='a')
sage: a[eU,0,1] = x*y^2 + 2*x
sage: a.add_comp_by_continuation(eV, W, c_uv)
sage: a.retrict(U)
Tensor field a of type (2,0) on the Open subset U of the 2-dimensional differentiable manifold M

But it should return

2-vector field a on the Open subset U of the 2-dimensional differentiable manifold M

[mjungmath]

comment:93

Okay, I spotted the error. It's caused by these lines:

+        self._sym = tuple(self._sym)
+        self._antisym = tuple(self._antisym)

When a tensor field is displayed, it must be restricted to a parallelizable subset first. If that restriction does not exists, it's constructed from scratch. This is done by the tensor method of the corresponding vector field module of this subset. In particular, the parameters sym and antisym are passed to this method. In our particular case the code will be executed until

        elif tensor_type[0] > 1 and tensor_type[1] == 0 and antisym:
            if isinstance(antisym[0], (int, Integer)):
                # a single antisymmetry is provided as a tuple or a
                # range object; it is converted to a 1-item list:
                antisym = [tuple(antisym)]
            if isinstance(antisym, list):
                antisym0 = antisym[0]
            else:
                antisym0 = antisym
            if len(antisym0) == tensor_type[0]:
                return self.alternating_contravariant_tensor(
                                              tensor_type[0], name=name,
                                              latex_name=latex_name)

With the above change, however, len(antisym0) returns 1, but it should be 2 because the code turns [(0, 1)] into ((0,1),). So, this case will be ignored and a tensor without symmetries will be created instead.

[mjungmath]

comment:94

One way to solve this could be to adapt the preprocessing of antisym0 and sym0 in the tensor method, assuming antisym and sym is a tuple.

[mkoeppe]

comment:95

Thanks for spotting this!

[sagetrac-git]

Changed commit from edd1572 to e5a81f2

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

fef81a5	FiniteRankFreeModule, VectorFieldModule: Accept both lists and tuples for sym, antisym
922b360	TensorField, PseudoRiemannianMetric: Adjust to normalization of _sym, _antisym to tuples
e5a81f2	sage.manifolds: Update doctest output for trivial output changes

[mkoeppe]

comment:97

With these changes, only failures that look like this remain:

sage -t --random-seed=0 src/sage/manifolds/differentiable/examples/sphere.py
**********************************************************************
File "src/sage/manifolds/differentiable/examples/sphere.py", line 56, in sage.manifolds.differentiable.examples.sphere
Failed example:
    h.display()
Exception raised:
    Traceback (most recent call last):
      File "/Users/mkoeppe/s/sage/sage-rebasing/worktree-gcc11/src/sage/doctest/forker.py", line 714, in _run
        self.compile_and_execute(example, compiler, test.globs)
      File "/Users/mkoeppe/s/sage/sage-rebasing/worktree-gcc11/src/sage/doctest/forker.py", line 1133, in compile_and_execute
        exec(compiled, globs)
      File "<doctest sage.manifolds.differentiable.examples.sphere[5]>", line 1, in <module>
        h.display()
      File "/Users/mkoeppe/s/sage/sage-rebasing/worktree-gcc11/src/sage/manifolds/differentiable/tensorfield.py", line 1852, in display
        return rst.display(frame, chart)
      File "/Users/mkoeppe/s/sage/sage-rebasing/worktree-gcc11/src/sage/tensor/modules/free_module_tensor.py", line 723, in display
        coef = comp[ind_arg]
      File "/Users/mkoeppe/s/sage/sage-rebasing/worktree-gcc11/src/sage/tensor/modules/comp_element_dict.py", line 2533, in __getitem__
        return self._output_formatter(self._comp[ind])
      File "/Users/mkoeppe/s/sage/sage-rebasing/worktree-gcc11/src/sage/manifolds/scalarfield.py", line 1786, in coord_function
        raise ValueError("no starting chart could be found to " +
    ValueError: no starting chart could be found to compute the expression in the Chart (E^3, (x, y, z))

[mkoeppe]

comment:98

This is likely not the final form of this ticket.

The goals of this ticket -- one-time precomputation of data related to the symmetries -- have not been fully achieved yet because the parent class is now a module over a base ring. This is fine for numerical tensors that are all over the same ring (QQ or RDF) but for the symbolic tensors that are used in sage.manifolds, there are many base rings, so the symmetries will be recomputed for each of them.

But for now it was more important to me to have separate parents for separate base rings, as this will allow us to dispatch to different implementation backends (= element classes) - as part of @honglizhaobob's project (#31991).

My solution for allowing more shared precomputation would be to introduce another class, similar to https://docs.sympy.org/latest/modules/tensor/tensor.html#sympy.tensor.tensor.TensorSymmetry, which only captures the symmetry group (with action).

Apart from this, there is still some code that should be pushed from element to parent, for example the symmetries from tensor contraction.

And there are branches in the code that can no longer be reached because the coercion system, via CompParent._element_constructor_, now guarantees that the inputs of single-underscore methods such as _add_ have the same symmetry already.

Also, I have not done any time measurements.

[mkoeppe]

comment:100

Further refactoring will go through #32029 (Action of a sympy TensorSymmetry)

[mjungmath]

comment:101

Looks good so far!

What about the idea with the morphisms in comment:28?

[mkoeppe]

comment:102

Replying to @mjungmath:

What about the idea with the morphisms in comment:28?

The next step into this direction should be to fix/complete the existing code for coercions.

sage: cp = CompParentWithSym(QQ, 4, sym=((1, 2), (3, 4)))
sage: cp2 = CompParentWithSym(QQ, 4, sym=((1, 2, 3, 4)))
sage: cp.coerce_map_from(cp2)   # returns None, should return injection to cp

Help with this is welcome!

[mkoeppe]

comment:107

Restarting this effort with a smaller step in #34497.