hyperbolic_geodesic midpoint bugfix

[Javier Honrubia González]

midpoint() fails to end for some geodesic with finite length example:

sage: g=HyperbolicPlane().UHP().get_geodesic(2*I,2*I+1)
sage: g.midpoint()'
...
RuntimeError: maximum recursion depth exceeded in __instancecheck__

[Javier Honrubia González]

The problem seems to be triggered when calculating the inverse of the _to_std_geod(start) matrix

[Travis Scrimshaw]

It seems that the matrix is sufficiently complicated that it just runs up against the Python function call limit. If I run simplify_full on the matrix first, then it works (and I can even run in again to simplify it further). So this would be one way to hack around this.

Although we might benefit from having a special case for inverting 2x2 matrices using the explicit form:

[A B]^-1 = __1__ [ D -B]
[C D]      AD-BC [-C  A]

This would give another way around this (it should be faster and in some ways, more robust) as this code only deals with 2x2 matrices.

[Javier Honrubia González]

I had tried simplify_full (although I hadn't thought of applying it multiple times) but I am somehow disappointed with the solution even that it gets since the result is still too complicated (in the motivating example should have Re(z)=1/2 instead of the monster that it shows. I am trying to understand the code to see how to get a simpler expression, if possible without resorting to simplify_full the result. The transformation S it uses to carry the geodesic to the imaginary axis it is not even recognized as an isometry by the model, fails the test self._model.isometry_in_model(S) (although the 'full_simplified' expression successes) Obviously works flawlessly with

h=HyperbolicPlane().UHP().get_geodesic(2.0*I,2.0*I+1.0)
h.midpoint()

[Javier Honrubia González]

The matrix of the isometry that sends the geodesic to the imaginary axis is a symbolic array, when this matrix is numeric the algorithm evaluates efficiently but when is purely symbolic it is necessary a double full_simplify() operation so that the posterior operation of inverting the matrix does not get a runtime error. After profiling the code, I do not see (surprisingly) any speed advantage in coding explicitly the inverse matrix, so I decided to left that unchanged.

---

New commits:

​89500eb

If the S matrix is not numeric, apply full_simplify() in order to avoid the runtime error.

[Vincent Delecroix]

Hello,

Why in all this code the type of number is not preserved

sage: g = HyperbolicPlane().UHP().get_geodesic(CC(0,2),CC(1,2))
sage: type(g.endpoints()[0].coordinates())
<type 'sage.rings.complex_number.ComplexNumber'>
sage: type(g.midpoint().coordinates())
<type 'sage.symbolic.expression.Expression'>

If I am using floating point number it is likely that I want to do floating point computations.

[Vincent Delecroix]

Related question: in your code your assume that the coordinates are symbolic. Why can you make such assumption (since I am able to construct geodesic with any kind of coordinates)?

[Javier Honrubia González]

Replying to vdelecroix:

Hello,

Why in all this code the type of number is not preserved

sage: g = HyperbolicPlane().UHP().get_geodesic(CC(0,2),CC(1,2))
sage: type(g.endpoints()[0].coordinates())
<type 'sage.rings.complex_number.ComplexNumber'>
sage: type(g.midpoint().coordinates())
<type 'sage.symbolic.expression.Expression'>

If I am using floating point number it is likely that I want to do floating point computations.

Since I was trying to fix a bug, I did not question the way the module was implemented in the beginning (I assumed it was discussed, and reached a consensus). Taking a closer look, the effect you mention is caused by _crossratio_matrix() as in

sage: (p1, p2, p3) = [UHP.random_point().coordinates() for k in range(3)]
sage: A = HyperbolicGeodesicUHP._crossratio_matrix(p1, p2, p3)
sage: type(A)
<type 'sage.matrix.matrix_symbolic_dense.Matrix_symbolic_dense'>

even

sage: B = HyperbolicGeodesicUHP._crossratio_matrix(CC(1,1), CC(3,2), CC(2,2))
sage: type(B)
<type 'sage.matrix.matrix_generic_dense.Matrix_generic_dense'>

One thing we can try is to check the type of p1,p2,p3 and construct the matrix in CDF if we have isinstance(p1,numbers.Complex) using

return matrix(CDF,[[p1 - p2, (p1 - p2)*(-p0)],[p1 - p0, (p1 - p0)*(-p2)]])

What do you think of this approach? Any pitfalls ahead I am not aware of?

[Javier Honrubia González]

Replying to vdelecroix:

Related question: in your code your assume that the coordinates are symbolic. Why can you make such assumption (since I am able to construct geodesic with any kind of coordinates)?

Well, it seems that all the code is working on symbolics as in

sage: UHP = HyperbolicPlane().UHP()
sage: P = UHP.random_point().coordinates()
sage: type(P)
<type 'sage.symbolic.expression.Expression'>

[Javier Honrubia González]

On a related topic. If I check with number.Complex I get unwanted results as in

sage: import numbers
sage: from sage.rings.complex_number import ComplexNumber
sage: Q = UHP.get_point(CC(sqrt(2)*(1+I)))
sage: Q.coordinates()
1.41421356237309 + 1.41421356237309*I
sage: isinstance(Q.coordinates(),numbers.Complex)
False
sage: isinstance(Q.coordinates(),ComplexNumber)
True

[Javier Honrubia González]

I correct myself the crossratio matrix is <type 'sage.matrix.matrix_generic_dense.Matrix_generic_dense'> it gets converted to symbolic when multiplied by B = matrix([[1, 0], [0,I]]) so if we change it to matrix([[1r,0r],[0r,CC(0,1)]]) we mantain the type to <type 'sage.matrix.matrix_generic_dense.Matrix_generic_dense'> It is what you are requesting or it you require to get a <type 'sage.matrix.matrix_complex_double_dense.Matrix_complex_double_dense'>

In that case, will we have to construct

return matrix(CDF,[[p1 - p2, (p1 - p2)*(-p0)],[p1 - p0, (p1 - p0)*(-p2)]])

as mentioned before when the points p1,p2,p3 are instances of ComplexNumber and leaving as it is in the rest of cases?

[Vincent Delecroix]

Replying to jhonrubia6:

On a related topic. If I check with number.Complex I get unwanted results as in

sage: import numbers
sage: from sage.rings.complex_number import ComplexNumber
sage: Q = UHP.get_point(CC(sqrt(2)*(1+I)))
sage: Q.coordinates()
1.41421356237309 + 1.41421356237309*I
sage: isinstance(Q.coordinates(),numbers.Complex)
False
sage: isinstance(Q.coordinates(),ComplexNumber)
True

Which version of Sage are you using? On sage-7.2.beta6 complex numbers are perfectly recognized by Python numbers.Complex:

sage: import numbers
sage: isinstance(CC(0,1), numbers.Complex)
True
sage: isinstance(CDF(0,1), numbers.Complex)
True

[Vincent Delecroix]

Replying to jhonrubia6:

I correct myself the crossratio matrix is <type 'sage.matrix.matrix_generic_dense.Matrix_generic_dense'> it gets converted to symbolic when multiplied by B = matrix([[1, 0], [0,I]]) so if we change it to matrix([[1r,0r],[0r,CC(0,1)]]) we mantain the type to <type 'sage.matrix.matrix_generic_dense.Matrix_generic_dense'> It is what you are requesting or it you require to get a <type 'sage.matrix.matrix_complex_double_dense.Matrix_complex_double_dense'>

In that case, will we have to construct

return matrix(CDF,[[p1 - p2, (p1 - p2)*(-p0)],[p1 - p0, (p1 - p0)*(-p2)]])

as mentioned before when the points p1,p2,p3 are instances of ComplexNumber and leaving as it is in the rest of cases?

It would be bad to use CC(0,1) by default for the complex imaginary unit. The only viable solution is to use the base ring of the coordinates

def get_I(a):
    from sage.structure.element import Element
    if isinstance(a, complex):
        return complex("j")
    elif isinstance(a, Expression):
        return SR("I")
    elif isinstance(a, Element):
        return a.parent().gen()
    else:
        raise ValueError("not a complex number")

It seems to work in most cases

sage: get_I(SR(1))
I
sage: parent(_)
Symbolic Ring

sage: get_I(QQbar(1+I))
I
sage: parent(_)
Algebraic Field

sage: get_I(CDF.an_element())
1.0*I
sage: parent(_)
Complex Double Field

sage: get_I(ComplexField(256).an_element())
1.000000000000000000000000000000000000000000000000000000000000000000000000000*I
sage: parent(_)
Complex Field with 256 bits of precision

[Javier Honrubia González]

Replying to vdelecroix:

Replying to jhonrubia6:

On a related topic. If I check with number.Complex I get unwanted results as in

sage: import numbers
sage: from sage.rings.complex_number import ComplexNumber
sage: Q = UHP.get_point(CC(sqrt(2)*(1+I)))
sage: Q.coordinates()
1.41421356237309 + 1.41421356237309*I
sage: isinstance(Q.coordinates(),numbers.Complex)
False
sage: isinstance(Q.coordinates(),ComplexNumber)
True

Which version of Sage are you using? On sage-7.2.beta6 complex numbers are perfectly recognized by Python numbers.Complex:

sage: import numbers
sage: isinstance(CC(0,1), numbers.Complex)
True
sage: isinstance(CDF(0,1), numbers.Complex)
True

I am not at my installation right now (I'll check later). I tried the previous code in cloud.sagemath.com

sage: import numbers
sage: isinstance(CC(0,1), numbers.Complex)
False
sage: isinstance(CDF(0,1), numbers.Complex)
False

But then, maybe cloud.sagemath.com is not in the appropiate version. As I said I'll check later

[Javier Honrubia González]

Replying to vdelecroix:

Replying to jhonrubia6:

I correct myself the crossratio matrix is <type 'sage.matrix.matrix_generic_dense.Matrix_generic_dense'> it gets converted to symbolic when multiplied by B = matrix([[1, 0], [0,I]]) so if we change it to matrix([[1r,0r],[0r,CC(0,1)]]) we mantain the type to <type 'sage.matrix.matrix_generic_dense.Matrix_generic_dense'> It is what you are requesting or it you require to get a <type 'sage.matrix.matrix_complex_double_dense.Matrix_complex_double_dense'>

In that case, will we have to construct

return matrix(CDF,[[p1 - p2, (p1 - p2)*(-p0)],[p1 - p0, (p1 - p0)*(-p2)]])

as mentioned before when the points p1,p2,p3 are instances of ComplexNumber and leaving as it is in the rest of cases?

It would be bad to use CC(0,1) by default for the complex imaginary unit. The only viable solution is to use the base ring of the coordinates

def get_I(a):
    from sage.structure.element import Element
    if isinstance(a, complex):
        return complex("j")
    elif isinstance(a, Expression):
        return SR("I")
    elif isinstance(a, Element):
        return a.parent().gen()
    else:
        raise ValueError("not a complex number")

It seems to work in most cases

sage: get_I(SR(1))
I
sage: parent(_)
Symbolic Ring

sage: get_I(QQbar(1+I))
I
sage: parent(_)
Algebraic Field

sage: get_I(CDF.an_element())
1.0*I
sage: parent(_)
Complex Double Field

sage: get_I(ComplexField(256).an_element())
1.000000000000000000000000000000000000000000000000000000000000000000000000000*I
sage: parent(_)
Complex Field with 256 bits of precision

Understood. I appreciate your insight, I have not the experience in Sage programming to foresee these subleties. I'll introduce your code in the next patch

[git]

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

​afdf46e

Modified the way it calculates the B matrix that transforms (0,1,inf)->(0,I,inf) trying to preserve the type of the calculations

[Vincent Delecroix]

1. The following syntax is wrong

   TESTS::
   
   bla bla.
   
       sage: one_test()
   

   It should be

   TESTS:
   
   bla bla.::
   
       sage: one_test()   
   

2. You should get rid of trailing whitespaces, that is line that only contains spaces.

3. The .is_numeric() method only applies to symbolic expressions. There is something weird with your code if S[0,0].is_numeric() since you might also want to manipulate floating point objects.

[git]

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

​663a76c

Removed trailing blank spaces. Fixed the TEST syntax. Modified the matrix element check for symbolic matrix check

[Javier Honrubia González]

I've also modified the doctoring in _to_std_geod() method since it randomly failed (depending on the random UHP points, for example:

sage: UHP = HyperbolicPlane().UHP()
sage: p1=0.869685163129173 + 0.897487924526729*I
sage: p2=5.71917016919078 + 7.09666543326670*I
sage: g = UHP.get_geodesic(p1, p2)
sage: p= g.midpoint().coordinates()
sage: A = g._to_std_geod(p);A # Send midpoint to I.
sage: A = UHP.get_isometry(A)
sage: [s, e]= g.complete().endpoints();s;e;p
sage: bool(abs(A(s).coordinates()) < 10**-9)
True
sage: bool(abs(A(g.midpoint()).coordinates() - I) < 10**-9)
True
sage: bool(A(e).coordinates() == infinity)
False
sage: bool(abs(A(e).coordinates())>10**9)
True

[Vincent Delecroix]

• the syntax is for refering to trac tickets is :trac:`number`

• I don't understand your code. Instead of copy/paste you would better do

  if isinstance(S, Matrix_symbolic_dense):
      S = S.simplify_full()
  S_1 = S.inverse()
  ... etc ...
  

[Vincent Delecroix]

You might want to add some doctest similar to the following to check that floating points remain floating points in any method.

sage: g = UHP.get_geodesic(CC(0,1), CC(2,2))
sage: UHP.dist(g.start(), g.end())
1.45057451382258
sage: parent(_)
Real Field with 53 bits of precision
sage: g.midpoint()
Point in UHP 0.666666666666666 + 1.69967317119759*I
sage: parent(g.midpoint().coordinates())
Complex Field with 53 bits of precision
sage: gc = g.complete()
Geodesic in UHP from -0.265564437074637 to 3.76556443707464
sage: parent(gc.start().coordinates())
Real Field with 53 bits of precision
sage: parent(gc._to_std_geod(g.start().coordinates()))
Full MatrixSpace of 2 by 2 dense matrices over Complex Field with 53 bits of precision

[Javier Honrubia González]

Ok. I will try, but would it not be better to get its own trac ticket for that?

[git]

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

​62159c3

Code simplified as noted by reviewer. docstring syntax error corrected. New test regarding floating points added in methods.

[Vincent Delecroix]

Thanks for incorporing some (unrelated) changes for floating point testing.

There are still some issues with documentation

• The code in the tests/examples must be indented

  TESTS:
  
  This is a test::
  
      sage: 1+1
      2
  

  or

  TESTS::
  
      sage: 1+1
      2
  

  (the same holds for EXAMPLES)

• before any test there should be two colons (you forgot some at one place)

[git]

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

​929c7c7

Revised rst syntax. Tried to adhere to pep8 recommendations. Some long lines reaming though

[Javier Honrubia González]

I tried to fix long lines too (as noted by flake8) , but there are some I do not find a way to reduce, so I left them.

[Vincent Delecroix]

Still some documentation issues

• :meth:`midpoint()` should be :meth`midpoint`.

• before any tests there must be two colons. There is at least one missing occurrence after ... floating points in :meth:`dist()`.

[git]

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

​dd19383

Removed trailing blanks, added double colons before sage code, and revised indentation.

[Javier Honrubia González]

I think I got them all this time. Used pep8 checkers and a rest online checker.

[Vincent Delecroix]

Why did you do these kind of changes

-#***********************************************************************
+# ***********************************************************************

[Javier Honrubia González]

The line #******** produced an error E265 block comment should start with '# ' in the PEP8 check. The pep8 recommendation says "Each line of a block comment starts with a # and a single space (unless it is indented text inside the comment). Paragraphs inside a block comment are separated by a line containing a single # ."

[Vincent Delecroix]

As you noticed these are recommendations. Everywhere else in Sage, there is no space after the comment (especially in banners). It has to be consistent within Sage source code before being pep8 complient.

[git]

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

​8fea9c2

Removed blank spaces after #

[Javier Honrubia González]

Ok. I reverted those changes. What do you think about adding pictures illustrating the geodesics?. If so, I think we have to rename show methods by plot for sphinx_plot to work.

[Vincent Delecroix]

It would be nice be open a new ticket for that purpose.

[Javier Honrubia González]

Ok. I've opened #20530 for that.

[Vincent Delecroix]

The method get_B would be simpler and safer as follows

from sage.structure.element import Element
from sage.symbolic.expression import Expression

if isinstance(a, (int,float,complex)):  # Python number
    a = CDF(a)

if isinstance(a, Expression):           # symbolic
    P = SR
    zero = SR.zero()
    one = SR.one()
    I = SR("I")
elif isinstance(a, Element):            # Sage number
    P = a.parent()
    zero = P.zero()
    one = P.one()
    I = P.gen()
    if I.is_one() or (I*I).is_one() or not (-I*I).is_one():
        raise ValueError("invalid number")
else:
    raise ValueError("not a complex number")

return matrix(P, 2, [one, zero, zero, -I])

[git]

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

​2ba1c96

Repleaced the code for _get_B()

[Javier Honrubia González]

Ok. It seems nicer coding. I think you have contributed more code than me to this ticket Maybe we should exchange roles, put you as author and me as reviewer? ;-)

[Nils Bruin]

Once you've established which parent P to work over, you need a way to find I in there. Just as getting 0, 1, -1 into P is to construct them in the simplest parent (ZZ in this case) and converting them into P via P(0),P(1),P(-1), constructing I in P should be a similar procedure. There is one solution that works most of the time via cyclotomic fields: P(CyclotomicField(4).0) In a way, if that fails for P then P is probably not a good parent for this kind of computations.

There are some numerical noise problems with CC(CyclotomicField(4),0) which mightbe worth fixing.

Currently, GF(5)(CyclotomicField(4).0) doesn't work.

Code along these lines could be a lot more robust for new parents than the special-casing approach taken now. I'm not quite sure we have suitable infrastructure for it to use it now, but perhaps we should (the main issue is of course that the choice of I isn't unique, whereas the integers convert uniquely into any integral domain)

[Javier Honrubia González]

Replying to nbruin:

Once you've established which parent P to work over, you need a way to find I in there. Just as getting 0, 1, -1 into P is to construct them in the simplest parent (ZZ in this case) and converting them into P via P(0),P(1),P(-1), constructing I in P should be a similar procedure. There is one solution that works most of the time via cyclotomic fields: P(CyclotomicField(4).0) In a way, if that fails for P then P is probably not a good parent for this kind of computations.

There are some numerical noise problems with CC(CyclotomicField(4),0) which mightbe worth fixing.

Currently, GF(5)(CyclotomicField(4).0) doesn't work.

Code along these lines could be a lot more robust for new parents than the special-casing approach taken now. I'm not quite sure we have suitable infrastructure for it to use it now, but perhaps we should (the main issue is of course that the choice of I isn't unique, whereas the integers convert uniquely into any integral domain)

I am out of my field here, but if I understand you right, are you proposing to replace the elif block proposed by Vincent by something based on the CyclotomicField?? Why is it more robust? aside of the problems you mention of CC(CyclotomicField(4).0) which may be should be accomplished on a different ticket

[Vincent Delecroix]

One way to check the sign would be to check the embedding in CC as CC(custom_I) == CC.gen().

And an other way to find I would be (-P.one()).sqrt()... but currently there is no specification of this sqrt operation and sometimes you get a symbolic element where you do not want to.

sage: (-QQ.one()).sqrt().parent()
Symbolic Ring
sage: (-RR.one()).sqrt().parent()
Complex Field with 53 bits of precision
sage: (-AA.one()).sqrt().parent()
Algebraic Field
sage: K = QuadraticField(-1)
sage: (-K.one()).sqrt().parent()
Number Field in a with defining polynomial x^2 + 1
sage: K = QuadraticField(2)
sage: (-K.one()).sqrt().parent()
Symbolic Ring

I guess it is reasonable to have a temporary function in sage/rings/

def imaginary_unit(ring):
    r"""
    Return the imaginary unit in ``ring`` or in its algebraic closure if
    `-1` has no square root.

    If ``ring`` has a defined embedding in some complex field, then the
    imaginary unit is normalized so that its image in this complex field
    is its canonical imaginary unit.
    """
    ...

We could make it better later on and possibly remove if we find some other way. At the same time I would advocate that .sqrt() should return an element in the same ring or, if it is not possible, in its algebraic closure.

[Nils Bruin]

Replying to jhonrubia6:

I am out of my field here, but if I understand you right, are you proposing to replace the elif block proposed by Vincent by something based on the CyclotomicField?? Why is it more robust? aside of the problems you mention of CC(CyclotomicField(4).0) which may be should be accomplished on a different ticket

There is no particular reason for P.gen(0) to be a fourth root of unity, even if P contains a fourth root of unity. Therefore, the proposed code can easily fail if it doesn't really have to. Hence, it would be nice to express intent a little more explicitly here. For CyclotomicField(4) it's clear its generator is a fourth root of unity, so that's a candidate for a "universal" place to find such an element. I'm not so sure it's the best place, but it's at least a place to hook into the general conversion infrastructure.

[Javier Honrubia González]

The code

if isinstance(a, (int,float,complex)): # Python number
            a = CDF(a)

        if isinstance(a, Expression):          # symbolic
            P = SR
            zero = SR.zero()
            one = SR.one()
            I = SR("I")
        elif isinstance(a, Element):           # Sage number
            P = a.parent()
            zero = P(0)                        # Get 0 in Parent ring
            one = P(1)                         # Get 1 in Parent ring
            CF = CyclotomicField(4)            # CF.gen() is guaranteed to be I
            I = P(CF.gen())                    # Get I in Parent ring
            # In case there is some numerical 'noise'
            # e.g. P=CC ; P(CyclotomicField) == 6.123233995736766e-17 + 1.0*I
            I = 0.5*(I- I.conjugate())
            if I.is_one() or (I*I).is_one() or not (-I*I).is_one():
                raise ValueError("invalid number")
        else:
            raise ValueError("not a complex number")

        return matrix(P, 2, [one, zero, zero, -I])

fails when a = QQbar(1+I) resulting in

Exception raised:
    Traceback (most recent call last):
...
    TypeError: Illegal initializer for algebraic number

The offending code being the matrix(P, 2 [one, zero, zero, -I])

sage: P = QQbar(1+I).parent()
sage: zero = P(0); one = P(1)
sage: matrix(P, 2, [one, zero, zero, -P(CyclotomicField(4).gen())])
Exception raised:
    Traceback (most recent call last):
...
    TypeError: Illegal initializer for algebraic number
sage: matrix(P,2,[one,zero,zero,-P.gen()])
[ 1  0]
[ 0 -I]
sage: type(P(CyclotomicField(4).gen()));P(CyclotomicField(4).gen())
<class 'sage.rings.qqbar.AlgebraicNumber'>
I
sage: type(P.gen());P.gen()
<class 'sage.rings.qqbar.AlgebraicNumber'>
I

Do we have to open a ticket for that and leave this stopped or we leave the code as it is in the latest commit make it to the release and wait until the problems with CyclotomicField? (numerical noise, and this particular matrix generation) get solved to open a new ticket on this module? IMHO , the commit as it is solves the bug of the ticket (not introducing any new one) while both of your ideas go into the direction of a improvement of the code.

[Vincent Delecroix]

You are right. We should think about improvements for later tickets! This one already does her part of the job.