Sage: Ticket #20330: hyperbolic_geodesic midpoint bugfix
https://trac.sagemath.org/ticket/20330
<p>
midpoint() fails to end for some geodesic with finite length
example:
</p>
<pre class="wiki">sage: g=HyperbolicPlane().UHP().get_geodesic(2*I,2*I+1)
sage: g.midpoint()'
...
RuntimeError: maximum recursion depth exceeded in __instancecheck__
</pre>en-usSagehttps://trac.sagemath.org/chrome/site/logo_sagemath_trac.png
https://trac.sagemath.org/ticket/20330
Trac 1.2Javier Honrubia GonzálezWed, 30 Mar 2016 19:30:36 GMTcomponent, type changed; keywords set
https://trac.sagemath.org/ticket/20330#comment:1
https://trac.sagemath.org/ticket/20330#comment:1
<ul>
<li><strong>keywords</strong>
<em>hyperbolic</em> <em>geometry</em> <em>geodesic</em> added
</li>
<li><strong>component</strong>
changed from <em>PLEASE CHANGE</em> to <em>geometry</em>
</li>
<li><strong>type</strong>
changed from <em>PLEASE CHANGE</em> to <em>defect</em>
</li>
</ul>
TicketJavier Honrubia GonzálezWed, 30 Mar 2016 19:32:43 GMTdescription changed
https://trac.sagemath.org/ticket/20330#comment:2
https://trac.sagemath.org/ticket/20330#comment:2
<ul>
<li><strong>description</strong>
modified (<a href="/ticket/20330?action=diff&version=2">diff</a>)
</li>
</ul>
TicketJavier Honrubia GonzálezWed, 30 Mar 2016 19:33:27 GMTdescription changed
https://trac.sagemath.org/ticket/20330#comment:3
https://trac.sagemath.org/ticket/20330#comment:3
<ul>
<li><strong>description</strong>
modified (<a href="/ticket/20330?action=diff&version=3">diff</a>)
</li>
</ul>
TicketJavier Honrubia GonzálezWed, 30 Mar 2016 20:47:28 GMT
https://trac.sagemath.org/ticket/20330#comment:4
https://trac.sagemath.org/ticket/20330#comment:4
<p>
The problem seems to be triggered when calculating the inverse of the _to_std_geod(start) matrix
</p>
TicketTravis ScrimshawSat, 09 Apr 2016 21:04:47 GMT
https://trac.sagemath.org/ticket/20330#comment:5
https://trac.sagemath.org/ticket/20330#comment:5
<p>
It seems that the matrix is sufficiently complicated that it just runs up against the Python function call limit. If I run <code>simplify_full</code> 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.
</p>
<p>
Although we might benefit from having a special case for inverting 2x2 matrices using the explicit form:
</p>
<pre class="wiki">[A B]^-1 = __1__ [ D -B]
[C D] AD-BC [-C A]
</pre><p>
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.
</p>
TicketJavier Honrubia GonzálezWed, 27 Apr 2016 15:54:18 GMT
https://trac.sagemath.org/ticket/20330#comment:6
https://trac.sagemath.org/ticket/20330#comment:6
<p>
I had tried <code>simplify_full</code> (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 <code>simplify_full</code> 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 <code>self._model.isometry_in_model(S)</code> (although the double 'full_simplified' expression successes)
Obviously works flawlessly with
</p>
<pre class="wiki">h=HyperbolicPlane().UHP().get_geodesic(2.0*I,2.0*I+1.0)
h.midpoint()
</pre>
TicketJavier Honrubia GonzálezSun, 01 May 2016 18:22:16 GMTbranch set
https://trac.sagemath.org/ticket/20330#comment:7
https://trac.sagemath.org/ticket/20330#comment:7
<ul>
<li><strong>branch</strong>
set to <em>u/jhonrubia6/hyperbolic_geodesic_midpoint_bugfix</em>
</li>
</ul>
TicketJavier Honrubia GonzálezSun, 01 May 2016 18:41:59 GMTstatus changed; commit set
https://trac.sagemath.org/ticket/20330#comment:8
https://trac.sagemath.org/ticket/20330#comment:8
<ul>
<li><strong>status</strong>
changed from <em>new</em> to <em>needs_review</em>
</li>
<li><strong>commit</strong>
set to <em>89500eb900f245d091afcc4461f4d8f19f6db5b9</em>
</li>
</ul>
<p>
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.
</p>
<hr />
<p>
New commits:
</p>
<table class="wiki">
<tr><td><a class="ext-link" href="http://git.sagemath.org/sage.git/commit/?id=89500eb900f245d091afcc4461f4d8f19f6db5b9"><span class="icon"></span>89500eb</a></td><td><code>If the S matrix is not numeric, apply full_simplify() in order to avoid the runtime error.</code>
</td></tr></table>
TicketJavier Honrubia GonzálezSun, 01 May 2016 18:42:23 GMTauthor set
https://trac.sagemath.org/ticket/20330#comment:9
https://trac.sagemath.org/ticket/20330#comment:9
<ul>
<li><strong>author</strong>
set to <em>Javier Honrubia González</em>
</li>
</ul>
TicketVincent DelecroixSun, 01 May 2016 22:53:02 GMT
https://trac.sagemath.org/ticket/20330#comment:10
https://trac.sagemath.org/ticket/20330#comment:10
<p>
Hello,
</p>
<p>
Why in all this code the type of number is not preserved
</p>
<pre class="wiki">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'>
</pre><p>
If I am using floating point number it is likely that I want to do floating point computations.
</p>
TicketVincent DelecroixSun, 01 May 2016 22:54:53 GMT
https://trac.sagemath.org/ticket/20330#comment:11
https://trac.sagemath.org/ticket/20330#comment:11
<p>
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)?
</p>
TicketJavier Honrubia GonzálezMon, 02 May 2016 06:34:39 GMTstatus changed
https://trac.sagemath.org/ticket/20330#comment:12
https://trac.sagemath.org/ticket/20330#comment:12
<ul>
<li><strong>status</strong>
changed from <em>needs_review</em> to <em>needs_work</em>
</li>
</ul>
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/20330#comment:10" title="Comment 10">vdelecroix</a>:
</p>
<blockquote class="citation">
<p>
Hello,
</p>
<p>
Why in all this code the type of number is not preserved
</p>
<pre class="wiki">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'>
</pre><p>
If I am using floating point number it is likely that I want to do floating point computations.
</p>
</blockquote>
<p>
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 <code>_crossratio_matrix()</code> as in
</p>
<pre class="wiki">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'>
</pre><p>
even
</p>
<pre class="wiki">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'>
</pre><p>
One thing we can try is to check the type of p1,p2,p3 and construct the matrix in CDF if we have <code>isinstance(p1,numbers.Complex)</code> using
</p>
<pre class="wiki">return matrix(CDF,[[p1 - p2, (p1 - p2)*(-p0)],[p1 - p0, (p1 - p0)*(-p2)]])
</pre><p>
What do you think of this approach? Any pitfalls ahead I am not aware of?
</p>
TicketJavier Honrubia GonzálezMon, 02 May 2016 06:37:19 GMT
https://trac.sagemath.org/ticket/20330#comment:13
https://trac.sagemath.org/ticket/20330#comment:13
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/20330#comment:11" title="Comment 11">vdelecroix</a>:
</p>
<blockquote class="citation">
<p>
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)?
</p>
</blockquote>
<p>
Well, it seems that all the code is working on symbolics as in
</p>
<pre class="wiki">sage: UHP = HyperbolicPlane().UHP()
sage: P = UHP.random_point().coordinates()
sage: type(P)
<type 'sage.symbolic.expression.Expression'>
</pre>
TicketJavier Honrubia GonzálezMon, 02 May 2016 06:54:11 GMT
https://trac.sagemath.org/ticket/20330#comment:14
https://trac.sagemath.org/ticket/20330#comment:14
<p>
On a related topic. If I check with number.Complex I get unwanted results as in
</p>
<pre class="wiki">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
</pre>
TicketJavier Honrubia GonzálezMon, 02 May 2016 07:40:28 GMT
https://trac.sagemath.org/ticket/20330#comment:15
https://trac.sagemath.org/ticket/20330#comment:15
<p>
I correct myself the crossratio matrix is
<code><type 'sage.matrix.matrix_generic_dense.Matrix_generic_dense'></code>
it gets converted to symbolic when multiplied by
<code>B = matrix([[1, 0], [0,I]])</code>
so if we change it to
<code>matrix([[1r,0r],[0r,CC(0,1)]])</code>
we mantain the type to
<code><type 'sage.matrix.matrix_generic_dense.Matrix_generic_dense'></code>
It is what you are requesting or it you require to get a
<code><type 'sage.matrix.matrix_complex_double_dense.Matrix_complex_double_dense'></code>
</p>
<p>
In that case, will we have to construct
</p>
<pre class="wiki">return matrix(CDF,[[p1 - p2, (p1 - p2)*(-p0)],[p1 - p0, (p1 - p0)*(-p2)]])
</pre><p>
as mentioned before when the points p1,p2,p3 are instances of <a class="wiki" href="https://trac.sagemath.org/wiki/ComplexNumber">ComplexNumber</a> and leaving as it is in the rest of cases?
</p>
TicketVincent DelecroixMon, 02 May 2016 15:10:20 GMT
https://trac.sagemath.org/ticket/20330#comment:16
https://trac.sagemath.org/ticket/20330#comment:16
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/20330#comment:14" title="Comment 14">jhonrubia6</a>:
</p>
<blockquote class="citation">
<p>
On a related topic. If I check with number.Complex I get unwanted results as in
</p>
<pre class="wiki">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
</pre></blockquote>
<p>
Which version of Sage are you using? On sage-7.2.beta6 complex numbers are perfectly recognized by Python <code>numbers.Complex</code>:
</p>
<pre class="wiki">sage: import numbers
sage: isinstance(CC(0,1), numbers.Complex)
True
sage: isinstance(CDF(0,1), numbers.Complex)
True
</pre>
TicketVincent DelecroixMon, 02 May 2016 15:21:46 GMT
https://trac.sagemath.org/ticket/20330#comment:17
https://trac.sagemath.org/ticket/20330#comment:17
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/20330#comment:15" title="Comment 15">jhonrubia6</a>:
</p>
<blockquote class="citation">
<p>
I correct myself the crossratio matrix is
<code><type 'sage.matrix.matrix_generic_dense.Matrix_generic_dense'></code>
it gets converted to symbolic when multiplied by
<code>B = matrix([[1, 0], [0,I]])</code>
so if we change it to
<code>matrix([[1r,0r],[0r,CC(0,1)]])</code>
we mantain the type to
<code><type 'sage.matrix.matrix_generic_dense.Matrix_generic_dense'></code>
It is what you are requesting or it you require to get a
<code><type 'sage.matrix.matrix_complex_double_dense.Matrix_complex_double_dense'></code>
</p>
<p>
In that case, will we have to construct
</p>
<pre class="wiki">return matrix(CDF,[[p1 - p2, (p1 - p2)*(-p0)],[p1 - p0, (p1 - p0)*(-p2)]])
</pre><p>
as mentioned before when the points p1,p2,p3 are instances of <a class="wiki" href="https://trac.sagemath.org/wiki/ComplexNumber">ComplexNumber</a> and leaving as it is in the rest of cases?
</p>
</blockquote>
<p>
It would be bad to use <code>CC(0,1)</code> by default for the complex imaginary unit. The only viable solution is to use the base ring of the coordinates
</p>
<pre class="wiki">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")
</pre><p>
It seems to work in most cases
</p>
<pre class="wiki">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
</pre>
TicketJavier Honrubia GonzálezMon, 02 May 2016 15:32:38 GMT
https://trac.sagemath.org/ticket/20330#comment:18
https://trac.sagemath.org/ticket/20330#comment:18
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/20330#comment:16" title="Comment 16">vdelecroix</a>:
</p>
<blockquote class="citation">
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/20330#comment:14" title="Comment 14">jhonrubia6</a>:
</p>
<blockquote class="citation">
<p>
On a related topic. If I check with number.Complex I get unwanted results as in
</p>
<pre class="wiki">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
</pre></blockquote>
<p>
Which version of Sage are you using? On sage-7.2.beta6 complex numbers are perfectly recognized by Python <code>numbers.Complex</code>:
</p>
<pre class="wiki">sage: import numbers
sage: isinstance(CC(0,1), numbers.Complex)
True
sage: isinstance(CDF(0,1), numbers.Complex)
True
</pre></blockquote>
<p>
I am not at my installation right now (I'll check later). I tried the previous code in cloud.sagemath.com
</p>
<pre class="wiki">sage: import numbers
sage: isinstance(CC(0,1), numbers.Complex)
False
sage: isinstance(CDF(0,1), numbers.Complex)
False
</pre><p>
But then, maybe cloud.sagemath.com is not in the appropiate version. As I said I'll check later
</p>
TicketJavier Honrubia GonzálezMon, 02 May 2016 15:34:09 GMT
https://trac.sagemath.org/ticket/20330#comment:19
https://trac.sagemath.org/ticket/20330#comment:19
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/20330#comment:17" title="Comment 17">vdelecroix</a>:
</p>
<blockquote class="citation">
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/20330#comment:15" title="Comment 15">jhonrubia6</a>:
</p>
<blockquote class="citation">
<p>
I correct myself the crossratio matrix is
<code><type 'sage.matrix.matrix_generic_dense.Matrix_generic_dense'></code>
it gets converted to symbolic when multiplied by
<code>B = matrix([[1, 0], [0,I]])</code>
so if we change it to
<code>matrix([[1r,0r],[0r,CC(0,1)]])</code>
we mantain the type to
<code><type 'sage.matrix.matrix_generic_dense.Matrix_generic_dense'></code>
It is what you are requesting or it you require to get a
<code><type 'sage.matrix.matrix_complex_double_dense.Matrix_complex_double_dense'></code>
</p>
<p>
In that case, will we have to construct
</p>
<pre class="wiki">return matrix(CDF,[[p1 - p2, (p1 - p2)*(-p0)],[p1 - p0, (p1 - p0)*(-p2)]])
</pre><p>
as mentioned before when the points p1,p2,p3 are instances of <a class="wiki" href="https://trac.sagemath.org/wiki/ComplexNumber">ComplexNumber</a> and leaving as it is in the rest of cases?
</p>
</blockquote>
<p>
It would be bad to use <code>CC(0,1)</code> by default for the complex imaginary unit. The only viable solution is to use the base ring of the coordinates
</p>
<pre class="wiki">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")
</pre><p>
It seems to work in most cases
</p>
<pre class="wiki">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
</pre></blockquote>
<p>
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
</p>
TicketgitThu, 05 May 2016 21:24:06 GMTcommit changed
https://trac.sagemath.org/ticket/20330#comment:20
https://trac.sagemath.org/ticket/20330#comment:20
<ul>
<li><strong>commit</strong>
changed from <em>89500eb900f245d091afcc4461f4d8f19f6db5b9</em> to <em>afdf46e9237f644889eb993a6236c31e91c6cf4d</em>
</li>
</ul>
<p>
Branch pushed to git repo; I updated commit sha1. New commits:
</p>
<table class="wiki">
<tr><td><a class="ext-link" href="http://git.sagemath.org/sage.git/commit/?id=afdf46e9237f644889eb993a6236c31e91c6cf4d"><span class="icon"></span>afdf46e</a></td><td><code>Modified the way it calculates the B matrix that transforms (0,1,inf)->(0,I,inf) trying to preserve the type of the calculations</code>
</td></tr></table>
TicketVincent DelecroixFri, 06 May 2016 01:41:47 GMTreviewer set
https://trac.sagemath.org/ticket/20330#comment:21
https://trac.sagemath.org/ticket/20330#comment:21
<ul>
<li><strong>reviewer</strong>
set to <em>Vincent Delecroix</em>
</li>
</ul>
<ol><li>The following syntax is wrong
<pre class="wiki">TESTS::
bla bla.
sage: one_test()
</pre>It should be
<pre class="wiki">TESTS:
bla bla.::
sage: one_test()
</pre></li></ol><ol start="2"><li>You should get rid of trailing whitespaces, that is line that only contains spaces.
</li></ol><ol start="3"><li>The <code>.is_numeric()</code> method only applies to symbolic expressions. There is something weird with your code <code>if S[0,0].is_numeric()</code> since you might also want to manipulate floating point objects.
</li></ol>
TicketgitSun, 08 May 2016 11:49:39 GMTcommit changed
https://trac.sagemath.org/ticket/20330#comment:22
https://trac.sagemath.org/ticket/20330#comment:22
<ul>
<li><strong>commit</strong>
changed from <em>afdf46e9237f644889eb993a6236c31e91c6cf4d</em> to <em>663a76cb67fb1598301ce842dccd8341ab07218b</em>
</li>
</ul>
<p>
Branch pushed to git repo; I updated commit sha1. New commits:
</p>
<table class="wiki">
<tr><td><a class="ext-link" href="http://git.sagemath.org/sage.git/commit/?id=663a76cb67fb1598301ce842dccd8341ab07218b"><span class="icon"></span>663a76c</a></td><td><code>Removed trailing blank spaces. Fixed the TEST syntax. Modified the matrix element check for symbolic matrix check</code>
</td></tr></table>
TicketJavier Honrubia GonzálezSun, 08 May 2016 16:32:52 GMT
https://trac.sagemath.org/ticket/20330#comment:23
https://trac.sagemath.org/ticket/20330#comment:23
<p>
I've also modified the doctoring in <code>_to_std_geod()</code> method since it randomly failed (depending on the random UHP points, for example:
</p>
<pre class="wiki">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
</pre>
TicketJavier Honrubia GonzálezSun, 08 May 2016 16:33:22 GMTstatus changed
https://trac.sagemath.org/ticket/20330#comment:24
https://trac.sagemath.org/ticket/20330#comment:24
<ul>
<li><strong>status</strong>
changed from <em>needs_work</em> to <em>needs_review</em>
</li>
</ul>
TicketVincent DelecroixSun, 08 May 2016 18:22:32 GMT
https://trac.sagemath.org/ticket/20330#comment:25
https://trac.sagemath.org/ticket/20330#comment:25
<ul><li>the syntax is for refering to trac tickets is <code>:trac:`number`</code>
</li></ul><ul><li>I don't understand your code. Instead of copy/paste you would better do
<pre class="wiki">if isinstance(S, Matrix_symbolic_dense):
S = S.simplify_full()
S_1 = S.inverse()
... etc ...
</pre></li></ul>
TicketVincent DelecroixSun, 08 May 2016 18:52:36 GMT
https://trac.sagemath.org/ticket/20330#comment:26
https://trac.sagemath.org/ticket/20330#comment:26
<p>
You might want to add some doctest similar to the following to check that floating points remain floating points in any method.
</p>
<pre class="wiki">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
</pre>
TicketJavier Honrubia GonzálezSun, 08 May 2016 19:23:08 GMT
https://trac.sagemath.org/ticket/20330#comment:27
https://trac.sagemath.org/ticket/20330#comment:27
<p>
Ok. I will try, but would it not be better to get its own trac ticket for that?
</p>
TicketgitSun, 08 May 2016 20:58:32 GMTcommit changed
https://trac.sagemath.org/ticket/20330#comment:28
https://trac.sagemath.org/ticket/20330#comment:28
<ul>
<li><strong>commit</strong>
changed from <em>663a76cb67fb1598301ce842dccd8341ab07218b</em> to <em>62159c386ed9ea95663b50606ae770f3e1afdbee</em>
</li>
</ul>
<p>
Branch pushed to git repo; I updated commit sha1. New commits:
</p>
<table class="wiki">
<tr><td><a class="ext-link" href="http://git.sagemath.org/sage.git/commit/?id=62159c386ed9ea95663b50606ae770f3e1afdbee"><span class="icon"></span>62159c3</a></td><td><code>Code simplified as noted by reviewer. docstring syntax error corrected. New test regarding floating points added in methods.</code>
</td></tr></table>
TicketVincent DelecroixTue, 10 May 2016 04:53:19 GMTstatus changed
https://trac.sagemath.org/ticket/20330#comment:29
https://trac.sagemath.org/ticket/20330#comment:29
<ul>
<li><strong>status</strong>
changed from <em>needs_review</em> to <em>needs_work</em>
</li>
</ul>
<p>
Thanks for incorporing some (unrelated) changes for floating point testing.
</p>
<p>
There are still some issues with documentation
</p>
<ul><li>The code in the tests/examples must be indented
<pre class="wiki">TESTS:
This is a test::
sage: 1+1
2
</pre>or
<pre class="wiki">TESTS::
sage: 1+1
2
</pre>(the same holds for EXAMPLES)
</li></ul><ul><li>before any test there should be two colons (you forgot some at one place)
</li></ul>
TicketgitTue, 10 May 2016 20:16:04 GMTcommit changed
https://trac.sagemath.org/ticket/20330#comment:30
https://trac.sagemath.org/ticket/20330#comment:30
<ul>
<li><strong>commit</strong>
changed from <em>62159c386ed9ea95663b50606ae770f3e1afdbee</em> to <em>929c7c73d9968adb10414dedf879060eb235728f</em>
</li>
</ul>
<p>
Branch pushed to git repo; I updated commit sha1. New commits:
</p>
<table class="wiki">
<tr><td><a class="ext-link" href="http://git.sagemath.org/sage.git/commit/?id=929c7c73d9968adb10414dedf879060eb235728f"><span class="icon"></span>929c7c7</a></td><td><code>Revised rst syntax. Tried to adhere to pep8 recommendations. Some long lines reaming though</code>
</td></tr></table>
TicketJavier Honrubia GonzálezTue, 10 May 2016 20:17:01 GMTstatus changed
https://trac.sagemath.org/ticket/20330#comment:31
https://trac.sagemath.org/ticket/20330#comment:31
<ul>
<li><strong>status</strong>
changed from <em>needs_work</em> to <em>needs_review</em>
</li>
</ul>
<p>
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.
</p>
TicketVincent DelecroixThu, 12 May 2016 05:01:44 GMTstatus changed
https://trac.sagemath.org/ticket/20330#comment:32
https://trac.sagemath.org/ticket/20330#comment:32
<ul>
<li><strong>status</strong>
changed from <em>needs_review</em> to <em>needs_work</em>
</li>
</ul>
<p>
Still some documentation issues
</p>
<ul><li><code>:meth:`midpoint()`</code> should be <code>:meth`midpoint`</code>.
</li></ul><ul><li>before any tests there must be two colons. There is at least one missing occurrence after <code>... floating points in :meth:`dist()`.</code>
</li></ul>
TicketgitSun, 15 May 2016 17:02:24 GMTcommit changed
https://trac.sagemath.org/ticket/20330#comment:33
https://trac.sagemath.org/ticket/20330#comment:33
<ul>
<li><strong>commit</strong>
changed from <em>929c7c73d9968adb10414dedf879060eb235728f</em> to <em>dd193839547305ed790c523e905c3e3f63ce563a</em>
</li>
</ul>
<p>
Branch pushed to git repo; I updated commit sha1. New commits:
</p>
<table class="wiki">
<tr><td><a class="ext-link" href="http://git.sagemath.org/sage.git/commit/?id=dd193839547305ed790c523e905c3e3f63ce563a"><span class="icon"></span>dd19383</a></td><td><code>Removed trailing blanks, added double colons before sage code, and revised indentation.</code>
</td></tr></table>
TicketJavier Honrubia GonzálezSun, 15 May 2016 17:04:54 GMTstatus changed
https://trac.sagemath.org/ticket/20330#comment:34
https://trac.sagemath.org/ticket/20330#comment:34
<ul>
<li><strong>status</strong>
changed from <em>needs_work</em> to <em>needs_review</em>
</li>
</ul>
<p>
I think I got them all this time. Used pep8 checkers and a rest online checker.
</p>
TicketVincent DelecroixMon, 16 May 2016 23:39:36 GMT
https://trac.sagemath.org/ticket/20330#comment:35
https://trac.sagemath.org/ticket/20330#comment:35
<p>
Why did you do these kind of changes
</p>
<div class="wiki-code"><div class="code"><pre><span class="gd">-#***********************************************************************
</span><span class="gi">+# ***********************************************************************
</span></pre></div></div>
TicketJavier Honrubia GonzálezTue, 17 May 2016 18:01:52 GMT
https://trac.sagemath.org/ticket/20330#comment:36
https://trac.sagemath.org/ticket/20330#comment:36
<p>
The line <code>#********</code> produced an error
<strong> E265 block comment should start with '# </strong>'
in the PEP8 check.
The pep8 recommendation says
"<em>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 # .</em>"
</p>
TicketVincent DelecroixTue, 17 May 2016 18:31:13 GMT
https://trac.sagemath.org/ticket/20330#comment:37
https://trac.sagemath.org/ticket/20330#comment:37
<p>
As you noticed these are <em>recommendations</em>. 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.
</p>
TicketgitWed, 18 May 2016 17:43:27 GMTcommit changed
https://trac.sagemath.org/ticket/20330#comment:38
https://trac.sagemath.org/ticket/20330#comment:38
<ul>
<li><strong>commit</strong>
changed from <em>dd193839547305ed790c523e905c3e3f63ce563a</em> to <em>8fea9c2fd8182498e593fbd2e03e151ebefdbea0</em>
</li>
</ul>
<p>
Branch pushed to git repo; I updated commit sha1. New commits:
</p>
<table class="wiki">
<tr><td><a class="ext-link" href="http://git.sagemath.org/sage.git/commit/?id=8fea9c2fd8182498e593fbd2e03e151ebefdbea0"><span class="icon"></span>8fea9c2</a></td><td><code>Removed blank spaces after #</code>
</td></tr></table>
TicketJavier Honrubia GonzálezWed, 18 May 2016 17:46:46 GMT
https://trac.sagemath.org/ticket/20330#comment:39
https://trac.sagemath.org/ticket/20330#comment:39
<p>
Ok. I reverted those changes. What do you think about adding pictures illustrating the geodesics?. If so, I think we have to rename <code>show</code> methods by <code>plot</code> for <code>sphinx_plot</code> to work.
</p>
TicketVincent DelecroixWed, 18 May 2016 17:48:35 GMT
https://trac.sagemath.org/ticket/20330#comment:40
https://trac.sagemath.org/ticket/20330#comment:40
<p>
It would be nice be open a new ticket for that purpose.
</p>
TicketJavier Honrubia GonzálezWed, 18 May 2016 17:51:27 GMT
https://trac.sagemath.org/ticket/20330#comment:41
https://trac.sagemath.org/ticket/20330#comment:41
<p>
Ok. I've opened <a class="closed ticket" href="https://trac.sagemath.org/ticket/20530" title="#20530: enhancement: Add pictures to hyperbolic_geodesic.py (closed: fixed)">#20530</a> for that.
</p>
TicketVincent DelecroixWed, 18 May 2016 18:14:16 GMT
https://trac.sagemath.org/ticket/20330#comment:42
https://trac.sagemath.org/ticket/20330#comment:42
<p>
The method <code>get_B</code> would be simpler and safer as follows
</p>
<pre class="wiki">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])
</pre>
TicketgitWed, 18 May 2016 18:31:26 GMTcommit changed
https://trac.sagemath.org/ticket/20330#comment:43
https://trac.sagemath.org/ticket/20330#comment:43
<ul>
<li><strong>commit</strong>
changed from <em>8fea9c2fd8182498e593fbd2e03e151ebefdbea0</em> to <em>2ba1c96e25d7e128adc832bc4c4e4c51e8835b98</em>
</li>
</ul>
<p>
Branch pushed to git repo; I updated commit sha1. New commits:
</p>
<table class="wiki">
<tr><td><a class="ext-link" href="http://git.sagemath.org/sage.git/commit/?id=2ba1c96e25d7e128adc832bc4c4e4c51e8835b98"><span class="icon"></span>2ba1c96</a></td><td><code>Repleaced the code for _get_B()</code>
</td></tr></table>
TicketJavier Honrubia GonzálezWed, 18 May 2016 18:33:30 GMT
https://trac.sagemath.org/ticket/20330#comment:44
https://trac.sagemath.org/ticket/20330#comment:44
<p>
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? ;-)
</p>
TicketNils BruinWed, 18 May 2016 19:20:03 GMT
https://trac.sagemath.org/ticket/20330#comment:45
https://trac.sagemath.org/ticket/20330#comment:45
<p>
Once you've established which parent P to work over, you need a way to find <code>I</code> 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 <code>P(0),P(1),P(-1)</code>, constructing <code>I</code> in P should be a similar procedure. There is one solution that works most of the time via cyclotomic fields: <code>P(CyclotomicField(4).0)</code>
In a way, if that fails for <code>P</code> then <code>P</code> is probably not a good parent for this kind of computations.
</p>
<p>
There are some numerical noise problems with <code>CC(CyclotomicField(4),0)</code> which mightbe worth fixing.
</p>
<p>
Currently, <code>GF(5)(CyclotomicField(4).0)</code> doesn't work.
</p>
<p>
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)
</p>
TicketJavier Honrubia GonzálezWed, 18 May 2016 19:48:21 GMT
https://trac.sagemath.org/ticket/20330#comment:46
https://trac.sagemath.org/ticket/20330#comment:46
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/20330#comment:45" title="Comment 45">nbruin</a>:
</p>
<blockquote class="citation">
<p>
Once you've established which parent P to work over, you need a way to find <code>I</code> 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 <code>P(0),P(1),P(-1)</code>, constructing <code>I</code> in P should be a similar procedure. There is one solution that works most of the time via cyclotomic fields: <code>P(CyclotomicField(4).0)</code>
In a way, if that fails for <code>P</code> then <code>P</code> is probably not a good parent for this kind of computations.
</p>
<p>
There are some numerical noise problems with <code>CC(CyclotomicField(4),0)</code> which mightbe worth fixing.
</p>
<p>
Currently, <code>GF(5)(CyclotomicField(4).0)</code> doesn't work.
</p>
<p>
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)
</p>
</blockquote>
<p>
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 <a class="missing wiki">CyclotomicField?</a>? Why is it more robust? aside of the problems you mention of <code>CC(CyclotomicField(4).0)</code> which may be should be accomplished on a different ticket
</p>
TicketVincent DelecroixWed, 18 May 2016 20:02:28 GMT
https://trac.sagemath.org/ticket/20330#comment:47
https://trac.sagemath.org/ticket/20330#comment:47
<p>
One way to check the sign would be to check the embedding in <code>CC</code> as <code>CC(custom_I) == CC.gen()</code>.
</p>
<p>
And an other way to find <code>I</code> would be <code>(-P.one()).sqrt()</code>... but currently there is no specification of this <code>sqrt</code> operation and sometimes you get a symbolic element where you do not want to.
</p>
<pre class="wiki">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
</pre><p>
I guess it is reasonable to have a temporary function in <code>sage/rings/</code>
</p>
<pre class="wiki">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.
"""
...
</pre><p>
We could make it better later on and possibly remove if we find some other way. At the same time I would advocate that <code>.sqrt()</code> should return an element in the same ring or, if it is not possible, in its algebraic closure.
</p>
TicketNils BruinWed, 18 May 2016 22:39:39 GMT
https://trac.sagemath.org/ticket/20330#comment:48
https://trac.sagemath.org/ticket/20330#comment:48
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/20330#comment:46" title="Comment 46">jhonrubia6</a>:
</p>
<blockquote class="citation">
<p>
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 <a class="missing wiki">CyclotomicField?</a>? Why is it more robust? aside of the problems you mention of <code>CC(CyclotomicField(4).0)</code> which may be should be accomplished on a different ticket
</p>
</blockquote>
<p>
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 <code>CyclotomicField(4)</code> 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.
</p>
TicketJavier Honrubia GonzálezWed, 25 May 2016 19:10:55 GMT
https://trac.sagemath.org/ticket/20330#comment:49
https://trac.sagemath.org/ticket/20330#comment:49
<p>
The code
</p>
<pre class="wiki">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])
</pre><p>
fails when <code>a = QQbar(1+I)</code> resulting in
</p>
<pre class="wiki">Exception raised:
Traceback (most recent call last):
...
TypeError: Illegal initializer for algebraic number
</pre><p>
The offending code being the <code>matrix(P, 2 [one, zero, zero, -I])</code>
</p>
<pre class="wiki">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
</pre><p>
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 <a class="missing wiki">CyclotomicField?</a> (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.
</p>
TicketVincent DelecroixWed, 25 May 2016 20:14:44 GMTstatus, milestone changed
https://trac.sagemath.org/ticket/20330#comment:50
https://trac.sagemath.org/ticket/20330#comment:50
<ul>
<li><strong>status</strong>
changed from <em>needs_review</em> to <em>positive_review</em>
</li>
<li><strong>milestone</strong>
changed from <em>sage-7.2</em> to <em>sage-7.3</em>
</li>
</ul>
<p>
You are right. We should think about improvements for later tickets! This one already does her part of the job.
</p>
TicketVolker BraunThu, 26 May 2016 07:15:22 GMTstatus, branch changed; resolution set
https://trac.sagemath.org/ticket/20330#comment:51
https://trac.sagemath.org/ticket/20330#comment:51
<ul>
<li><strong>status</strong>
changed from <em>positive_review</em> to <em>closed</em>
</li>
<li><strong>resolution</strong>
set to <em>fixed</em>
</li>
<li><strong>branch</strong>
changed from <em>u/jhonrubia6/hyperbolic_geodesic_midpoint_bugfix</em> to <em>2ba1c96e25d7e128adc832bc4c4e4c51e8835b98</em>
</li>
</ul>
Ticket