Opened 14 months ago
Closed 11 months ago
#32953 closed defect (fixed)
Sphere: improve handling of default charts
Reported by:  ghtobiasdiez  Owned by:  

Priority:  major  Milestone:  sage9.6 
Component:  manifolds  Keywords:  sphere 
Cc:  tscrim, nthiery, ghmjungmath, egourgoulhon  Merged in:  
Authors:  Eric Gourgoulhon  Reviewers:  Travis Scrimshaw 
Report Upstream:  N/A  Work issues:  
Branch:  20b9cb0 (Commits, GitHub, GitLab)  Commit:  20b9cb0cf890c1a717605b9e043ba398203a05e0 
Dependencies:  Stopgaps: 
Description (last modified by )
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 2sphere 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}. }
Change History (30)
comment:1 followup: 2 Changed 14 months ago by
comment:2 Changed 14 months ago by
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 2sphere 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 followup: 5 Changed 14 months ago by
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 orientationpreserving. 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 orientationpreserving and otherwise they have to differ (thatt is, a 0cochain trivializating the Z_2valued orientation 1cocycle). Using this extension, one could simply say that the chart on SPNP gets a minus sign instead of needing to introduce a new frame.
What do you think?
comment:4 Changed 14 months ago by
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 userfriendly?)
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 Changed 14 months ago by
Replying to ghtobiasdiez:
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) ....: 2form eps_g on the Open subset S^2{SP} of the 2sphere S^2 of radius 1 smoothly embedded in the Euclidean space E^3 Coordinate frame (S^2{SP}, (∂/∂yp1,∂/∂yp2)) 2form eps_g on the Open subset S^2{NP} of the 2sphere 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 2indices 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) ....: 2form eps_g on the Open subset S^2{SP} of the 2sphere S^2 of radius 1 smoothly embedded in the Euclidean space E^3 Coordinate frame (S^2{SP}, (∂/∂yp1,∂/∂yp2)) 2form eps_g on the Open subset S^2{NP} of the 2sphere 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 14 months ago by
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 followups: 8 9 Changed 14 months ago by
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 followup: 10 Changed 14 months ago by
Replying to ghtobiasdiez:
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 usingcomp
, 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 followup: 11 Changed 14 months ago by
Replying to ghtobiasdiez:
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.
comment:10 Changed 14 months ago by
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 followups: 13 15 Changed 14 months ago by
Replying to egourgoulhon:
Replying to ghtobiasdiez:
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 expressionH(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 betweenfunction('H')
and the scalar fieldH
. For the sake of clarity, different symbols shall be used. Keep in mind thatfunction('H')
injectsH
in the global namespace. By running subsequentlyH = M.scalar_field(...)
, you overwrite the Python nameH
. This can be unfortunate when you want to further manipulatefunction('H')
, like insubstitute_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 14 months ago by
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 followups: 14 16 Changed 14 months ago by
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.
comment:14 Changed 14 months ago by
Replying to ghmjungmath:
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 followup: 21 Changed 14 months ago by
Replying to ghtobiasdiez:
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, thenH(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 followup: 17 Changed 14 months ago by
Replying to ghmjungmath:
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 followup: 19 Changed 14 months ago by
Replying to egourgoulhon:
Replying to ghmjungmath:
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 beXN
.
Why is this output unexpected? Wouldn't you expect open subsets U
to have the same default chart as their ambient space M
?
comment:18 Changed 14 months ago by
I find this default chart concept in case of nonparallelizable 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 followup: 22 Changed 14 months ago by
Replying to ghmjungmath:
One would rather expect
XN.domain().default_chart()
to beXN
.Why is this output unexpected? Wouldn't you expect open subsets
U
to have the same default chart as their ambient spaceM
?
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 14 months ago by
Summary:  Sphere: vector vs coordinate frame for volume form → Sphere: improve handling of default charts 

comment:21 Changed 14 months ago by
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, thenH(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 give 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.
comment:22 followup: 23 Changed 14 months ago by
Replying to egourgoulhon:
No, not in general. Especially in the present case, where the spherical chart,
Xspher
says, does not coverU = XN.domain()
, but onlyA
, which is a strict subset ofU
. IfU
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 isXN
. 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 ofSphere.__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 UIf one creates a new open subset of
M
, it does not inherit ofXspher
as a default chart:sage: V = M.open_subset('V') sage: V.default_chart() is None Truewhich is fortunate, since
V
might not containA
.IMHO, the solution would be that
Sphere._init_stereographic
setsXN
as the default chart onU
, just after having created it. Similarly,Sphere._init_spherical
should setXspher
as the default chart onA
.
Don't you think this should be attacked on the level of the general manifolds implementation? Those (nonobvious) things can still happen if you define manifolds on the fly.
comment:23 Changed 14 months ago by
Replying to ghmjungmath:
IMHO, the solution would be that
Sphere._init_stereographic
setsXN
as the default chart onU
, just after having created it. Similarly,Sphere._init_spherical
should setXspher
as the default chart onA
.Don't you think this should be attacked on the level of the general manifolds implementation? Those (nonobvious) 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 followup: 27 Changed 14 months ago by
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 14 months ago by
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 14 months ago by
Milestone:  sage9.5 → sage9.6 

comment:27 Changed 11 months ago by
Authors:  → Eric Gourgoulhon 

Branch:  → public/manifolds/sphere_default_charts32953 
Commit:  → 20b9cb0cf890c1a717605b9e043ba398203a05e0 
Keywords:  sphere added 
Status:  new → needs_review 
Replying to ghmjungmath:
I agree, it wouldn't hurt to add this for spheres.
Here we go...
New commits:
20b9cb0  Improve handling of default charts on the sphere (#32953)

comment:28 Changed 11 months ago by
Reviewers:  → Travis Scrimshaw 

Status:  needs_review → positive_review 
Green bot (morally). I am setting this to a positive review. Feel free to revert if there are any objections.
comment:30 Changed 11 months ago by
Branch:  public/manifolds/sphere_default_charts32953 → 20b9cb0cf890c1a717605b9e043ba398203a05e0 

Resolution:  → fixed 
Status:  positive_review → closed 
Replying to ghtobiasdiez:
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 ofsrc/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:
can be shorten to
manifolds
providing a direct access to the manifold catalog.Note that the class
Sphere
is not extensively tested and might have some bugs.