Opened 12 months ago

Closed 9 months ago

#32953 closed defect (fixed)

Sphere: improve handling of default charts

Reported by: Tobias Diez Owned by:
Priority: major Milestone: sage-9.6
Component: manifolds Keywords: sphere
Cc: Travis Scrimshaw, Nicolas M. Thiéry, Michael Jung, Eric Gourgoulhon Merged in:
Authors: Eric Gourgoulhon Reviewers: Travis Scrimshaw
Report Upstream: N/A Work issues:
Branch: 20b9cb0 (Commits, GitHub, GitLab) Commit: 20b9cb0cf890c1a717605b9e043ba398203a05e0
Dependencies: Stopgaps:

Status badges

Description (last modified by Tobias Diez)

If one uses the method stereographic_coordinates to initialize steographic coordinates (instead of as via the constructor argument), then there are some inconsistencies.

For example,

sage: M = manifolds.Sphere(2)
sage: XN = M.stereographic_coordinates(pole='north')
sage: g = M.metric()
sage: U = XN.domain(); U
Open subset S^2-{NP} of the 2-sphere S^2 of radius 1 smoothly embedded in the 
 Euclidean space E^3
sage: g.restrict(U).display()
g = (cos(theta)^2 - 2*cos(theta) + 1) dy1⊗dy1
  + (cos(theta)^2 - 2*cos(theta) + 1) dy2⊗dy2

One would certainly expect the components of g to be expressed in terms of (y1,y2), i.e. the chart XN, which is the only chart that covers entirely U = S2-{NP}.

Change History (30)

comment:1 in reply to:  description ; Changed 12 months ago by Eric Gourgoulhon

Replying to gh-tobiasdiez:

As one can see, the volume form has components with respect to a coordinate frame on S2-{SP} but with respect to a vector frame on S2-{NP}.

This is because the volume form is defined with respect to the orientation provided by the frame of stereographic coordinates on S2-{SP}. The frame of stereographic coordinates on S2-{NP} having the opposite orientation, a vector frame (f1, f2) is introduced there to keep the orientation. As you can read on line 954 of src/sage/manifolds/differentiable/examples/sphere.py, the relation to the coordinate frame is very simple: f1 = -d/dy1, f2 = d/dy2. See also the comment on the orientation in the reference manual for Sphere.

A side remark:

from sage.manifolds.differentiable.examples.sphere import Sphere
M = Sphere(2)

can be shorten to

M = manifolds.Sphere(2)

manifolds providing a direct access to the manifold catalog.

Note that the class Sphere is not extensively tested and might have some bugs.

comment:2 in reply to:  1 Changed 12 months ago by Eric Gourgoulhon

Replying to egourgoulhon:

Note that the class Sphere is not extensively tested and might have some bugs.

An example of bug is

sage: M = manifolds.Sphere(2)
sage: XN = M.stereographic_coordinates(pole='north')
sage: g = M.metric()
sage: U = XN.domain(); U
Open subset S^2-{NP} of the 2-sphere S^2 of radius 1 smoothly embedded in the 
 Euclidean space E^3
sage: g.restrict(U).display()
g = (cos(theta)^2 - 2*cos(theta) + 1) dy1⊗dy1
  + (cos(theta)^2 - 2*cos(theta) + 1) dy2⊗dy2

One would certainly expect the components of g to be expressed in terms of (y1,y2), i.e. the chart XN, which is the only chart that covers entirely U = S^2-{NP}.

comment:3 Changed 12 months ago by Tobias Diez

Thanks for having a look at this.

So the vector frame f is only needed to specify the orienation? I guess then its easy to also calculate the components of the volumen form in the coordinate frame, since it's only a overall factor of minus 1.

I had a look at the implementation and would propose to circumvent this problem completely by slightly generalizing the way an orientation is represented in sage. Currently, it is as a list of frames for which the chart transitions are orientation-preserving. More generally one can assign to each chart (or frame) a sign, where two such signs have to be equal if the chart transition is orientation-preserving and otherwise they have to differ (thatt is, a 0-cochain trivializating the Z_2-valued orientation 1-cocycle). Using this extension, one could simply say that the chart on SP-NP gets a minus sign instead of needing to introduce a new frame.

What do you think?

Last edited 12 months ago by Tobias Diez (previous) (diff)

comment:4 Changed 12 months ago by Michael Jung

In general, this is a good idea. Then the question is, how do we assign these values when we introduce new charts? Currently, this is not necessary because we have reference frames/charts in the aforementioned list. Instead of giving every frame/chart a sign, I'd suggest to keep such a list of references but assign to each frame/chart only in that list a sign as you propose. (Is this user-friendly?)

I presume, however, a much better way would incorporate Chech cohomology, which goes in the direction of what you suggested in terms of cohomology. As soon as we have implemented Chech cohomology, the notion of "orientation" can easily be generalized. I think this is the proper way of doing it.

comment:5 in reply to:  3 Changed 12 months ago by Eric Gourgoulhon

Replying to gh-tobiasdiez:

So the vector frame f is only needed to specify the orienation? I guess then its easy to also calculate the components of the volumen form in the coordinate frame, since it's only a overall factor of minus 1.

Indeed, and Sage will do it automatically if required, because it knows all the automorphisms connecting the various frames. In particular, it knows the automorphism linking f to the coordinate frame on S^2-{NP}.

I had a look at the implementation and would propose to circumvent this problem completely by slightly generalizing the way an orientation is represented in sage.

I am not sure to understand what your problem is. At the initialization of the volume form, only the minimal amount of information is stored in the attribute _components. The user shall not access this private attribute directly, but should instead ask for components via the method comp() or display(). The required components are then automatically computed if they need to be and cached in _components. For instance:

sage: M = manifolds.Sphere(2, coordinates='stereographic')
sage: eps = M.metric().volume_form()
sage: for restr in eps._restrictions.values():
....:     print(restr)
....:     for comp in restr._components:
....:         print(comp)
....:
2-form eps_g on the Open subset S^2-{SP} of the 2-sphere S^2 of radius 1 smoothly embedded in the Euclidean space E^3
Coordinate frame (S^2-{SP}, (∂/∂yp1,∂/∂yp2))
2-form eps_g on the Open subset S^2-{NP} of the 2-sphere S^2 of radius 1 smoothly embedded in the Euclidean space E^3
Vector frame (S^2-{NP}, (f_1,f_2))

As you noticed, eps is initialized only w.r.t. f on S^2-{NP}. But you can access to the components in the coordinate frame by comp():

sage: Nf = M.stereographic_coordinates().frame(); Nf
Coordinate frame (S^2-{NP}, (∂/∂y1,∂/∂y2))
sage: eps.comp(Nf)
Fully antisymmetric 2-indices components w.r.t. Coordinate frame (S^2-{NP}, (∂/∂y1,∂/∂y2))

Note that the attribute _components has been updated:

sage: for restr in eps._restrictions.values():
....:     print(restr)
....:     for comp in restr._components:
....:         print(comp)
....:         
2-form eps_g on the Open subset S^2-{SP} of the 2-sphere S^2 of radius 1 smoothly embedded in the Euclidean space E^3
Coordinate frame (S^2-{SP}, (∂/∂yp1,∂/∂yp2))
2-form eps_g on the Open subset S^2-{NP} of the 2-sphere S^2 of radius 1 smoothly embedded in the Euclidean space E^3
Vector frame (S^2-{NP}, (f_1,f_2))
Coordinate frame (S^2-{NP}, (∂/∂y1,∂/∂y2))

Invoking display() instead of comp() would have produced the same result:

sage: eps.display(Nf)
eps_g = -4/(y1^4 + y2^4 + 2*(y1^2 + 1)*y2^2 + 2*y1^2 + 1) dy1∧dy2

comment:6 Changed 12 months ago by Michael Jung

I agree with Eric, there is no major problem here. Though I also agree with Tobias in the sense that introducing a dummy frame in order to get the orientation is not "nice" either.

comment:7 Changed 12 months ago by Tobias Diez

My (originial) problem was that the following code

import sage.all
from sage.manifolds.differentiable.examples.sphere import Sphere
from sage.symbolic.function_factory import function

M = Sphere(2)
stereoN = M.stereographic_coordinates(pole='north')
stereoS = M.stereographic_coordinates(pole='south')

print("Metric")
for restr in M.metric()._restrictions.values():
    print(restr)
    for comp in restr._components:
        print(comp)

print("Volume form")
for restr in M.metric().volume_form()._restrictions.values():
    print(restr)
    for comp in restr._components:
        print(comp)

chart_functions = {chart: function('H')(*chart[:]) for chart in M.atlas()}
H = M.scalar_field(chart_functions, name='H')
print(H.display())
X = M.metric().inverse().contract(H.differential())
print(X.display())

print("X")
for restr in X._restrictions.values():
    print(restr)
    for comp in restr._components:
        print(comp)

M.metric().volume_form().contract(X)

fails with the error "ValueError?: no common chart for the multiplication" (in the last line). I thought that this might be a result of the difference in vector vs coordinate frame. But I guess based on your comments, the issue is more that FreeModuleTensor#contract is directly using _components instead of trying to compute the components in the other frame using comp, right?

I've opened #32974 for the idea to improve the orientation of manifolds (using a general framework).

comment:8 in reply to:  7 ; Changed 12 months ago by Eric Gourgoulhon

Replying to gh-tobiasdiez:

My (originial) problem was that the following code ... fails with the error "ValueError?: no common chart for the multiplication" (in the last line). I thought that this might be a result of the difference in vector vs coordinate frame. But I guess based on your comments, the issue is more that FreeModuleTensor#contract is directly using _components instead of trying to compute the components in the other frame using comp, right?

No, FreeModuleTensor.contract is doing things correctly by invoking

basis = self.common_basis(other)

(cf. line 2826 of src/sage/tensor/modules/free_module_tensor.py), which computes the components in a common frame by calling comp if necessary. The issue you are facing is rather due to some flaw in the current implementation of Sphere, as already mentioned in comment:2. To circumvent it, you should use

M = manifolds.Sphere(2, coordinates='stereographic')

if you plan to work with stereographic coordinates. Then M.metric().volume_form().contract(X) returns no error.

comment:9 in reply to:  7 ; Changed 12 months ago by Eric Gourgoulhon

Replying to gh-tobiasdiez:

My (originial) problem was that the following code

...
chart_functions = {chart: function('H')(*chart[:]) for chart in M.atlas()}
H = M.scalar_field(chart_functions, name='H')

This declaration of the scalar field H is not mathematically correct: it enforces the same coordinate expression H(x,y) for all the coordinates (x,y) in the manifold's atlas, which cannot be, except for a constant scalar field. Furthermore, it introduces a confusion between function('H') and the scalar field H. For the sake of clarity, different symbols shall be used. Keep in mind that function('H') injects H in the global namespace. By running subsequently H = M.scalar_field(...), you overwrite the Python name H. This can be unfortunate when you want to further manipulate function('H'), like in substitute_function for instance.

Last edited 12 months ago by Eric Gourgoulhon (previous) (diff)

comment:10 in reply to:  8 Changed 12 months ago by Tobias Diez

Description: modified (diff)

Oh, then I misunderstood your comment above.

To circumvent it, you should use M = manifolds.Sphere(2, coordinates='stereographic')

This works perfectly, thanks!

I've updated the ticket description accordingly.

comment:11 in reply to:  9 ; Changed 12 months ago by Tobias Diez

Replying to egourgoulhon:

Replying to gh-tobiasdiez:

My (originial) problem was that the following code

...
chart_functions = {chart: function('H')(*chart[:]) for chart in M.atlas()}
H = M.scalar_field(chart_functions, name='H')

This declaration of the scalar field H is not mathematically correct: it enforces the same coordinate expression H(x,y) for all the coordinates (x,y) in the manifold's atlas, which cannot be, except for a constant scalar field. Furthermore, it introduces a confusion between function('H') and the scalar field H. For the sake of clarity, different symbols shall be used. Keep in mind that function('H') injects H in the global namespace. By running subsequently H = M.scalar_field(...), you overwrite the Python name H. This can be unfortunate when you want to further manipulate function('H'), like in substitute_function for instance.

What I intended to do was to create a new function in each chart, i.e. just represent a general function on a manifold. What would be the correct way to do this? Simply use a different name for each chart function, or can one disable the injection of the local function in the global namespace?

Side remark: I agree that using the same name might be confusing, but at least in the area I work in it seems standard to denote the chart representation of a function by the same symbol and make the distinction only by the arguments, i.e. if H: M \to \R is function on a manifold and \kappa: M \to \R^n a local chart, then H(x_1, ..., x_n) := H \circ \kappa^{-1}(x_1, ..., x_n).

comment:12 Changed 12 months ago by Tobias Diez

By the way, I really appreciate that you guys take the time to clarify and explain to me the inner workings of the manifold code!

comment:13 in reply to:  11 ; Changed 12 months ago by Michael Jung

Spheres are automatically initialized with spherical coordinates if not stated otherwise. That means, at the point of initialization, the default chart is set to the one given by spherical coordinates. I bet this to be the culprit of that "inconsistency".

If you want to initialize the sphere with stereographic coordinates right from the start, you can use manifolds.Sphere(2, coordinates='stereographic').

One would probably have to transcribe the display method in such a way that the chart used for displaying is tried to be obtained from (the restriction) of the chart which the frame is induced from before the default chart is used. This is just my guess from not looking into the code.

Incidental remark 1: Ticket #31894 attempts to solve the covering issue of spherical coordinates.

Incidental remark 2: remark 1 is also related to the orientation #31324 issue.

Last edited 12 months ago by Michael Jung (previous) (diff)

comment:14 in reply to:  13 Changed 12 months ago by Eric Gourgoulhon

Replying to gh-mjungmath:

One would probably has to transcribe the display method in such a way that the chart used for displaying is tried to be obtained from (the restriction) of the chart which the frame is induced from before the default chart is used. This is just my guess from not looking into the code.

This sounds a good idea!

comment:15 in reply to:  11 ; Changed 12 months ago by Eric Gourgoulhon

Replying to gh-tobiasdiez:

What I intended to do was to create a new function in each chart, i.e. just represent a general function on a manifold. What would be the correct way to do this?

I would say like this:

sage: M = manifolds.Sphere(2, coordinates='stereographic')
sage: XN.<y1, y2> = M.stereographic_coordinates(pole='north')
sage: XS = M.stereographic_coordinates(pole='south')
sage: H = M.scalar_field({XN: function('h')(y1, y2)}, name='H')
sage: H.add_expr_by_continuation(XS, XN.domain().intersection(XS.domain()))
sage: H.display()
H: S^2 → ℝ
on S^2-{NP}: (y1, y2) ↦ h(y1, y2)
on S^2-{SP}: (yp1, yp2) ↦ h(yp1/(yp1^2 + yp2^2), yp2/(yp1^2 + yp2^2))

Thanks to add_expr_by_continuation, the scalar field is defined in a consistent way on all the manifold.

Side remark: I agree that using the same name might be confusing, but at least in the area I work in it seems standard to denote the chart representation of a function by the same symbol and make the distinction only by the arguments, i.e. if H: M \to \R is function on a manifold and \kappa: M \to \R^n a local chart, then H(x_1, ..., x_n) := H \circ \kappa^{-1}(x_1, ..., x_n).

I see your point, but I am afraid that using this in a computer algebra system may lead to some false results.

comment:16 in reply to:  13 ; Changed 12 months ago by Eric Gourgoulhon

Replying to gh-mjungmath:

Spheres are automatically initialized with spherical coordinates if not stated otherwise. That means, at the point of initialization, the default chart is set to the one given by spherical coordinates. I bet this to be the culprit of that "inconsistency".

The issue is rather that the spherical chart appears as well as the default one on the domain of stereographic coordinates:

sage: M = manifolds.Sphere(2)
sage: XN = M.stereographic_coordinates(pole='north')
sage: M.default_chart()  # output is OK, given the initialization of M
Chart (A, (theta, phi))
sage: XN.domain().default_chart()  # output is unexpected
Chart (A, (theta, phi))

One would rather expect XN.domain().default_chart() to be XN.

comment:17 in reply to:  16 ; Changed 12 months ago by Michael Jung

Replying to egourgoulhon:

Replying to gh-mjungmath:

Spheres are automatically initialized with spherical coordinates if not stated otherwise. That means, at the point of initialization, the default chart is set to the one given by spherical coordinates. I bet this to be the culprit of that "inconsistency".

The issue is rather that the spherical chart appears as well as the default one on the domain of stereographic coordinates:

sage: M = manifolds.Sphere(2)
sage: XN = M.stereographic_coordinates(pole='north')
sage: M.default_chart()  # output is OK, given the initialization of M
Chart (A, (theta, phi))
sage: XN.domain().default_chart()  # output is unexpected
Chart (A, (theta, phi))

One would rather expect XN.domain().default_chart() to be XN.

Why is this output unexpected? Wouldn't you expect open subsets U to have the same default chart as their ambient space M?

Last edited 12 months ago by Michael Jung (previous) (diff)

comment:18 Changed 12 months ago by Michael Jung

I find this default chart concept in case of non-parallelizable manifolds a bit cumbersome, anyways.

What if we replace default_chart by default_charts, and use dictionaries instead? That is, the keys are subsets and values are default charts on that domain.

In addition, I'd opt for what I said earlier in comment:13.

comment:19 in reply to:  17 ; Changed 12 months ago by Eric Gourgoulhon

Replying to gh-mjungmath:

One would rather expect XN.domain().default_chart() to be XN.

Why is this output unexpected? Wouldn't you expect open subsets U to have the same default chart as their ambient space M?

No, not in general. Especially in the present case, where the spherical chart, Xspher says, does not cover U = XN.domain(), but only A, which is a strict subset of U. If U can be covered by a single chart, one would certainly expect that its default chart is such a covering chart. Moreover, in the current context, U is precisely defined as the domain of stereographic coordinates from the North pole (XN), so one certainly expects that its default chart is XN. That it is instead the chart of spherical coordinates happens only as an artifact of the order of creation of the open subsets and the spherical coordinate chart, in this part of Sphere.__init__:

self._init_chart_domains()   <- creates U, but not XN
self._init_embedding()
self._init_coordinates[coordinates](names)  <- creates Xspher and set it as the 
                                               default on all the preexisting 
                                               supersets of A, including U

If one creates a new open subset of M, it does not inherit of Xspher as a default chart:

sage: V = M.open_subset('V')
sage: V.default_chart() is None
True

which is fortunate, since V might not contain A.

IMHO, the solution would be that Sphere._init_stereographic sets XN as the default chart on U, just after having created it. Similarly, Sphere._init_spherical should set Xspher as the default chart on A.

comment:20 Changed 12 months ago by Eric Gourgoulhon

Summary: Sphere: vector vs coordinate frame for volume formSphere: improve handling of default charts

comment:21 in reply to:  15 Changed 12 months ago by Eric Gourgoulhon

Replying to egourgoulhon:

Thanks to add_expr_by_continuation, the scalar field is defined in a consistent way on all the manifold.

Side remark: I agree that using the same name might be confusing, but at least in the area I work in it seems standard to denote the chart representation of a function by the same symbol and make the distinction only by the arguments, i.e. if H: M \to \R is function on a manifold and \kappa: M \to \R^n a local chart, then H(x_1, ..., x_n) := H \circ \kappa^{-1}(x_1, ..., x_n).

I see your point, but I am afraid that using this in a computer algebra system may lead to some false results.

Here is an example of such false result. Suppose you define a scalar field the way you would like:

sage: M = manifolds.Sphere(2, coordinates='stereographic')
sage: XN.<y1, y2> = M.stereographic_coordinates(pole='north')
sage: XS.<yp1, yp2> = M.stereographic_coordinates(pole='south')
sage: H = M.scalar_field({XN: function('h')(y1, y2), 
....:                     XS: function('h')(yp1, yp2)}, name='H')

Then it will always take identical values at the North and South poles:

sage: NP = M((0,0), chart=XS)
sage: SP = M((0,0), chart=XN)
sage: bool(H(NP) == H(SP))
True

which is certainly not what you want for a generic scalar field on the sphere.

Last edited 12 months ago by Eric Gourgoulhon (previous) (diff)

comment:22 in reply to:  19 ; Changed 12 months ago by Michael Jung

Replying to egourgoulhon:

No, not in general. Especially in the present case, where the spherical chart, Xspher says, does not cover U = XN.domain(), but only A, which is a strict subset of U. If U can be covered by a single chart, one would certainly expect that its default chart is such a covering chart. Moreover, in the current context, U is precisely defined as the domain of stereographic coordinates from the North pole (XN), so one certainly expects that its default chart is XN. That it is instead the chart of spherical coordinates happens only as an artifact of the order of creation of the open subsets and the spherical coordinate chart, in this part of Sphere.__init__:

self._init_chart_domains()   <- creates U, but not XN
self._init_embedding()
self._init_coordinates[coordinates](names)  <- creates Xspher and set it as the 
                                               default on all the preexisting 
                                               supersets of A, including U

If one creates a new open subset of M, it does not inherit of Xspher as a default chart:

sage: V = M.open_subset('V')
sage: V.default_chart() is None
True

which is fortunate, since V might not contain A.

IMHO, the solution would be that Sphere._init_stereographic sets XN as the default chart on U, just after having created it. Similarly, Sphere._init_spherical should set Xspher as the default chart on A.

Don't you think this should be attacked on the level of the general manifolds implementation? Those (non-obvious) things can still happen if you define manifolds on the fly.

comment:23 in reply to:  22 Changed 12 months ago by Eric Gourgoulhon

Replying to gh-mjungmath:

IMHO, the solution would be that Sphere._init_stereographic sets XN as the default chart on U, just after having created it. Similarly, Sphere._init_spherical should set Xspher as the default chart on A.

Don't you think this should be attacked on the level of the general manifolds implementation? Those (non-obvious) things can still happen if you define manifolds on the fly.

I don't see any general strategy here. Probably we should leave this to the user defining the manifold on the fly. But in the specific case of the sphere, we should implement the standard charts (spherical, stereographic) in such a way that they are the default ones on their domains.

comment:24 Changed 12 months ago by Michael Jung

I agree, it wouldn't hurt to add this for spheres.

What about my idea before, using dictionaries instead?

Perhaps Travis has an idea?

comment:25 Changed 12 months ago by Travis Scrimshaw

I don't really like the idea of using dictionaries because you further increase the burden on keeping all of the information consistent when the information you want can just be directly attached to the subsets in question. So you have nC2 ~ n^2 pieces of information stored instead of n for a chain of n subsets.

comment:26 Changed 11 months ago by Matthias Köppe

Milestone: sage-9.5sage-9.6

comment:27 in reply to:  24 Changed 9 months ago by Eric Gourgoulhon

Authors: Eric Gourgoulhon
Branch: public/manifolds/sphere_default_charts-32953
Commit: 20b9cb0cf890c1a717605b9e043ba398203a05e0
Keywords: sphere added
Status: newneeds_review

Replying to gh-mjungmath:

I agree, it wouldn't hurt to add this for spheres.

Here we go...


New commits:

20b9cb0Improve handling of default charts on the sphere (#32953)

comment:28 Changed 9 months ago by Travis Scrimshaw

Reviewers: Travis Scrimshaw
Status: needs_reviewpositive_review

Green bot (morally). I am setting this to a positive review. Feel free to revert if there are any objections.

comment:29 Changed 9 months ago by Eric Gourgoulhon

Thank you for the review!

comment:30 Changed 9 months ago by Volker Braun

Branch: public/manifolds/sphere_default_charts-3295320b9cb0cf890c1a717605b9e043ba398203a05e0
Resolution: fixed
Status: positive_reviewclosed
Note: See TracTickets for help on using tickets.