sagemath / sage

Main repository of SageMath
https://www.sagemath.org
Other
1.45k stars 482 forks source link

Sphere: improve handling of default charts #32953

Closed tobiasdiez closed 2 years ago

tobiasdiez commented 2 years ago

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 = S^2-{NP}.

CC: @tscrim @nthiery @mjungmath @egourgoulhon

Component: manifolds

Keywords: sphere

Author: Eric Gourgoulhon

Branch/Commit: 20b9cb0

Reviewer: Travis Scrimshaw

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

egourgoulhon commented 2 years ago

Replying to @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.

egourgoulhon commented 2 years ago
comment:2

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}.

tobiasdiez commented 2 years ago
comment:3

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?

mjungmath commented 2 years ago
comment:4

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.

egourgoulhon commented 2 years ago
comment:5

Replying to @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
mjungmath commented 2 years ago
comment:6

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.

tobiasdiez commented 2 years ago
comment:7

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).

egourgoulhon commented 2 years ago
comment:8

Replying to @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.

egourgoulhon commented 2 years ago
comment:9

Replying to @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.

tobiasdiez commented 2 years ago
comment:10

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.

tobiasdiez commented 2 years ago

Description changed:

--- 
+++ 
@@ -1,59 +1,16 @@
-Calculating the volume form on the sphere uses a 'strange' vector frame on the lower hemisphere.
+If one uses the method `stereographic_coordinates` to initialize steographic coordinates (instead of as via the constructor argument), then there are some inconsistencies.

-In more detail, the following test script
+For example, 

-import sage.all -from sage.manifolds.differentiable.examples.sphere import Sphere

-M = Sphere(2) -M.stereographic_coordinates()

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

tobiasdiez commented 2 years ago
comment:11

Replying to @egourgoulhon:

Replying to @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).

tobiasdiez commented 2 years ago
comment:12

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!

mjungmath commented 2 years ago
comment:13

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.

egourgoulhon commented 2 years ago
comment:14

Replying to @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!

egourgoulhon commented 2 years ago
comment:15

Replying to @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.

egourgoulhon commented 2 years ago
comment:16

Replying to @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.

mjungmath commented 2 years ago
comment:17

Replying to @egourgoulhon:

Replying to @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?

mjungmath commented 2 years ago
comment:18

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.

egourgoulhon commented 2 years ago
comment:19

Replying to @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.

egourgoulhon commented 2 years ago
comment:21

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.

mjungmath commented 2 years ago
comment:22

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.

egourgoulhon commented 2 years ago
comment:23

Replying to @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.

mjungmath commented 2 years ago
comment:24

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

What about my idea before, using dictionaries instead?

Perhaps Travis has an idea?

tscrim commented 2 years ago
comment:25

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.

egourgoulhon commented 2 years ago

Commit: 20b9cb0

egourgoulhon commented 2 years ago

Author: Eric Gourgoulhon

egourgoulhon commented 2 years ago
comment:27

Replying to @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)
egourgoulhon commented 2 years ago

Branch: public/manifolds/sphere_default_charts-32953

egourgoulhon commented 2 years ago

Changed keywords from none to sphere

tscrim commented 2 years ago
comment:28

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

tscrim commented 2 years ago

Reviewer: Travis Scrimshaw

egourgoulhon commented 2 years ago
comment:29

Thank you for the review!

vbraun commented 2 years ago

Changed branch from public/manifolds/sphere_default_charts-32953 to 20b9cb0