Opened 10 months ago

Closed 9 months ago

#24884 closed defect (fixed)

Matrix-related fixes in differential geometry

Reported by: jdemeyer Owned by:
Priority: major Milestone: sage-8.2
Component: geometry Keywords:
Cc: egourgoulhon, tscrim Merged in:
Authors: Jeroen Demeyer, Eric Gourgoulhon Reviewers: Travis Scrimshaw
Report Upstream: N/A Work issues:
Branch: 6ae84f1 (Commits) Commit: 6ae84f147a2209e0396d508c77b4c86e706f3227
Dependencies: Stopgaps:

Description (last modified by jdemeyer)

These are changes which are needed because of #24742.

  1. src/sage/schemes/riemann_surfaces/riemann_surface.py: #24742 will disallow passing strings to the matrix constructor. It's a rarely used feature and hard to get right properly. The fix is easy: pass actual elements of the ring instead of the generator names. This is clearly the right thing to do.
  1. src/sage/tensor/modules/comp.py: I needed to change _get_list() but I have to say that I don't understand what the code is supposed to do. I wrote the code the way I did mainly to pass doctests. I also added some documentation to reflect better what the _get_list() method does. I also got rid of the evil try/except block trying to construct a matrix but happily returning a list if that failed.
  1. src/sage/manifolds/differentiable/metric.py: this is a consequence of the change in 2. and it looks like a sensible change.

Change History (22)

comment:1 Changed 10 months ago by jdemeyer

  • Branch set to u/jdemeyer/matrix_related_fixes_in_differential_geometry

comment:2 Changed 10 months ago by jdemeyer

  • Commit set to bd532ab6734f5a87df0c237f75105c92c1158c3d
  • Description modified (diff)
  • Status changed from new to needs_review

New commits:

bd532abMatrix-related fixes in differential geometry

comment:3 Changed 10 months ago by jdemeyer

  • Description modified (diff)

comment:4 Changed 10 months ago by jdemeyer

  • Description modified (diff)

comment:5 follow-up: Changed 10 months ago by tscrim

  • Reviewers set to Travis Scrimshaw
  • Status changed from needs_review to positive_review

Changes LGTM (Eric, if you see something you do not like, feel free to revert).

comment:6 Changed 10 months ago by egourgoulhon

  • Status changed from positive_review to needs_info

comment:7 in reply to: ↑ 5 ; follow-up: Changed 10 months ago by egourgoulhon

Replying to tscrim:

Changes LGTM (Eric, if you see something you do not like, feel free to revert).

Hi Travis, I was just looking to the ticket, at the same time as you apparently! The change in the doctest of src/sage/manifold/differentiable/metric.py (point 3 in the ticket description) is not desirable: it results from a forced conversion of SymPy to SR due to the change in src/sage/tensor/modules/comp.py (point 2). Jeroen, the line

return matrix(SR, resu)

forces the conversion of elements of the list resu to SR elements, while they can actually belong to something else, like SymPy. Is there a good reason for this? Do we have now to explicitly specify the ring as the first argument of matrix? Since SymPy elements do not constitute a ring in Sage, one cannot construct a matrix for them. I plan to improve this in the short future, by implementing a ring of SymPy coordinate functions. But for the time being, I would say it is preferable to return a list of SymPy elements instead of a matrix of SR elements.

comment:8 in reply to: ↑ 7 ; follow-up: Changed 9 months ago by jdemeyer

Replying to egourgoulhon:

Jeroen, the line

return matrix(SR, resu)

forces the conversion of elements of the list resu to SR elements, while they can actually belong to something else, like SymPy. Is there a good reason for this?

I just did that to be consistent with the documentation. The documentation never said "return a matrix, except when using the SymPy backend". I did make an exception for the str formatter, because that was documented (not in this method but elsewhere).

Do we have now to explicitly specify the ring as the first argument of matrix?

No, that won't be required. But again, it's good for consistency that you always get the same kind of object back (a symbolic matrix).

I would say it is preferable to return a list of SymPy elements instead of a matrix of SR elements.

That is fine for me, but then SymPy elements should always return a list. The main thing that I do not like about the current code is this try/except block:

            try:
                for i in range(self._dim):
                    for j in range(self._dim):
                        a = resu[i][j]
                        if hasattr(a, '_express'):
                            resu[i][j] = a.expr()
                resu = matrix(resu)  # for a nicer output
            except TypeError:
                pass

With this code, you really have no idea whether a matrix or a list will be returned. And if it's a matrix, you don't know over which base ring.

Travis, what do you think?

comment:9 Changed 9 months ago by egourgoulhon

  • Branch changed from u/jdemeyer/matrix_related_fixes_in_differential_geometry to public/manifolds/matrix_related_fixes_in_differential_geometry
  • Commit changed from bd532ab6734f5a87df0c237f75105c92c1158c3d to 6ae84f147a2209e0396d508c77b4c86e706f3227
  • Status changed from needs_info to needs_review

New commits:

6ae84f1Improve treatment of matrix output in sage.tensor.modules.comp.Components._get_list.

comment:10 in reply to: ↑ 8 ; follow-up: Changed 9 months ago by egourgoulhon

Replying to jdemeyer:

I just did that to be consistent with the documentation. The documentation never said "return a matrix, except when using the SymPy backend". I did make an exception for the str formatter, because that was documented (not in this method but elsewhere).

I agree the documentation of the method _get_list is pretty poor (but this is a "private" (underscored) method). Sorry about this. Note that the class Components to which it belongs does not know anything about SR or SymPy: it manipulates only elements of a ring. The SymPy or str objects appear only through the output formatter.

Do we have now to explicitly specify the ring as the first argument of matrix?

No, that won't be required. But again, it's good for consistency that you always get the same kind of object back (a symbolic matrix).

OK, thanks.

That is fine for me, but then SymPy elements should always return a list. The main thing that I do not like about the current code is this try/except block: With this code, you really have no idea whether a matrix or a list will be returned. And if it's a matrix, you don't know over which base ring.

I agree, this piece of code was not very good. Actually it pertains to the early time of the SageManifolds project, where the elements stored in Components where not guaranteed to belong to any ring. This is no longer the case for the Components used in tensor calculus: they all belong to the algebra of scalar fields and, after being processed by self._output_formatter, they belong to the ring of chart functions, sage.manifolds.chart_func.ChartFunctionRing whatever their internal representation (SR or Sympy). I've therefore modified the code of Components._get_list to

        if self._nid == 2:
            # 2-dim case: convert to matrix for a nicer output
            from sage.matrix.constructor import matrix
            from sage.structure.element import parent
            from sage.categories.rings import Rings
            if parent(resu[0][0]) in Rings():
                return matrix(resu)
        return resu

In this way, there is no special treatment for strings or whatever else that is not a ring, since the test if parent(resu[0][0]) in Rings() will return false in this case. For tensor field components, the treatment is now at the level of ChartFunctionRing (there is no longer any reference to the internal representation via _express and expr()), so it works for both SR and SymPy. In particular the doctest in src/sage/manifolds/differentiable/metric.py now displays a matrix of SymPy objects, and no longer a list, which is an improvement! Thanks for having triggered it!

I've pushed these changes in the new branch public/manifolds/matrix_related_fixes_in_differential_geometry, which is based on your branch (simply one commit ahead). I am attaching it to the ticket.

Last edited 9 months ago by egourgoulhon (previous) (diff)

comment:11 in reply to: ↑ 10 Changed 9 months ago by egourgoulhon

Replying to egourgoulhon:

An alternative to

            if parent(resu[0][0]) in Rings():
                return matrix(resu)

is

            R = parent(resu[0][0])
            if R in Rings():
                return matrix(R, resu)

I don't know which form is preferable. Note that in the latter, the ring is passed as the first argument to matrix.

comment:12 follow-ups: Changed 9 months ago by tscrim

IMO, a more ideal form is

            R = parent(resu[0][0])
            if R in Rings():
                return MatrixSpace(R,2)(resu)

Given comment:7

I plan to improve this in the short future, by implementing a ring of SymPy coordinate functions

I would also add a comment to the code saying to remove the Rings check once SymPy coordinate functions are a ring.

I think this also addresses Jeroen's comment about code flow clarity and the code LBTM (B for better). Any other comments Jeroen?

comment:13 Changed 9 months ago by jdemeyer

Instead of asking whether the parent is a ring, it seems more sensible to check for a Sage Element. i.e. I would suggest to replace

parent(x) in Rings()

by

isinstance(x, Element)

comment:14 in reply to: ↑ 12 ; follow-up: Changed 9 months ago by jdemeyer

Replying to tscrim:

IMO, a more ideal form is

            R = parent(resu[0][0])
            if R in Rings():
                return MatrixSpace(R,2)(resu)

What if the different elements have different parents? Is there a guarantee that the parent of resu[0][0] is the same as the parent of resu[0][1]? So matrix(resu) will auto-detect the correct parent (using the coercion model) and that looks better.

Of course the same argument could be made for the isinstance(resu[0][0], Element) check. But it seems less likely to mix Elements and non-Elements than to mix different parents.

comment:15 Changed 9 months ago by jdemeyer

  • Status changed from needs_review to needs_work

comment:16 in reply to: ↑ 12 Changed 9 months ago by egourgoulhon

Replying to tscrim:

I would also add a comment to the code saying to remove the Rings check once SymPy coordinate functions are a ring.

Actually, with the latest version, it works with SymPy, i.e. it does generate a matrix, the elements of which are chart functions, which are displayed as SymPy objects when SymPy is used as the symbolic engine (see the matrix containing x**2 and y**2 instead of x^2 and y^2 in the doctest in metric.py). The Rings check is then for other possible cases, not involving chart functions (I don't have any precise use case at the moment, except for the str output formatter used as illustration in the doctests). SymPy was a problem in the previous version of the code, where the matrix construction was attempted directly at the level of symbolic representations, either SR or SymPy. Now the matrix construction is at the (higher) level of chart functions, and, thanks to you IIRC ;-), chart functions now forms a ring, so that the matrix construction always succeeds.

comment:17 in reply to: ↑ 14 ; follow-up: Changed 9 months ago by egourgoulhon

Replying to jdemeyer:

What if the different elements have different parents? Is there a guarantee that the parent of resu[0][0] is the same as the parent of resu[0][1]? So matrix(resu) will auto-detect the correct parent (using the coercion model) and that looks better.

OK.

Of course the same argument could be made for the isinstance(resu[0][0], Element) check. But it seems less likely to mix Elements and non-Elements than to mix different parents.

A priori, all elements of resu should have the same parent: this is guaranted if no output formatter is used: the parent is nothing but the ring on which the Components object is constructed. If some output formatter is used, it is likely that its output is of a fixed type, hence giving rise to the same parent (in the extended sense, i.e. the parent default to the Python type if the object has no Sage Parent), like for instance the str output formatter. My concern with the Element check is that it does not guarantee that the parent is a ring. If it is not, then the matrix construction will fail.

comment:18 in reply to: ↑ 17 Changed 9 months ago by jdemeyer

  • Authors changed from Jeroen Demeyer to Jeroen Demeyer, Eric Gourgoulhon
  • Status changed from needs_work to needs_review

Replying to egourgoulhon:

If it is not, then the matrix construction will fail.

I didn't know that, but you are right. If you ask me, that's a silly artificial restriction.

In that case, maybe Eric's code is fine. Travis, what do you think?

comment:19 follow-up: Changed 9 months ago by tscrim

My understanding is that they should all be of the same parent. Yet, I don't hold a stake in this code, so I don't really care either way. I doubt this will ever make a difference in the wild (barring very esoteric uses).

comment:20 in reply to: ↑ 19 Changed 9 months ago by egourgoulhon

Replying to tscrim:

My understanding is that they should all be of the same parent. Yet, I don't hold a stake in this code, so I don't really care either way. I doubt this will ever make a difference in the wild (barring very esoteric uses).

Yes, I think it's pretty safe to leave the code in the current state.

comment:21 Changed 9 months ago by jdemeyer

  • Status changed from needs_review to positive_review

I'm taking that as positive review :-)

comment:22 Changed 9 months ago by vbraun

  • Branch changed from public/manifolds/matrix_related_fixes_in_differential_geometry to 6ae84f147a2209e0396d508c77b4c86e706f3227
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.