Automatic exclusion of non-domain points in things like arcsec

[kcrisman]

From an email exchange about plotting arcsec, slightly out of order:

The immediate context for this issue is a couple of problems that include, among other things, graphing arcsec and its derivative (in one case), arcsec and arccsc (in the other). Since the range of intended functions is known, it appears that this would work if the students follow instructions (big if). But more generally, we want students to have total control over their "graphing calculator" -- at least over the input boxes we provide -- so there's nothing stopping them from entering arcsec(x/2) or something else that we can't anticipate.

The underlying problem is that the plotting code ignores everything
between -1 and 1 (since those values aren't real), but then simply
connects (-1,3) and (1,.2) (or so) with a line segment when it's done. I
have a couple ideas for fixing this, but they won't help you in the near
term.

---

One workaround for your @interact would be to select an option for
exclude based on the function: I'm thinking of something like this:

@interact
def foo(fn=textbox(), ....):
setup stuff...

if fn == 'arcsin':
    exclude=[]
elif fn == 'arcsec':
    exclude=[-1,1]

plot(fn, ...., exclude=exclude)

The idea is just to look at the function being plotted and set a list of
points to exclude based on that. It would be tedious to code, but would
work.

---

I think that fixing this is possible. I spent about a half hour looking at http://hg.sagemath.org/sage-main/file/9ab4ab6e12d0/sage/plot/plot.py just now and I am pretty sure that one could extract bad points like nan in generate_plot_points. Currently we just do

data = [data[i] for i in range(len(data)) if i not in exception_indices]

which means we completely ignore them, but presumably we could keep this data somehow and return it as part of generate_plot_points (not sure whether that would have to be deprecated) so that we could add that somehow to the exclude data around http://hg.sagemath.org/sage-main/file/9ab4ab6e12d0/sage/plot/plot.py#l1279

The problem is that there could be a lot of these, and that could slow things down a lot to check all of them, but checking for a "consecutive series" would also take some time.

Anyway, it's at least doable in principle, but would take some time and care to implement, making sure it didn't cause any other regressions.

---

Apply to devel/sage:

1. attachment: trac_13246-rebase.patch
2. attachment: trac_13246-parametric.patch

Depends on #16439

CC: @dandrake @ppurka @jondo @vbraun

Component: graphics

Author: Punarbasu Purkayastha

Branch/Commit: 384b5a2

Reviewer: Ralf Stephan, Karl-Dieter Crisman

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

[kcrisman]

arcsec plot

[kcrisman]

comment:1

Attachment: tmp_1.png

Here is what the bad plot looks like. This function should not be defined in the middle.

[kcrisman]

comment:2

Hey ppurka, I know you like a challenge and know the plotting code pretty well... :)

[ppurka]

comment:3

Replying to @kcrisman:

Hey ppurka, I know you like a challenge and know the plotting code pretty well... :)

Hehe, I will have a look. :)

This problem had come up in sage-support or sage-devel too a couple of months back.

[dandrake]

comment:4

Here's an idea I thought of that would help avoid these ugly situations. The problem is that we build our plot out of lots of tiny line segments, and we have a line segment connecting things that it should not connect. So:

Let w = xmax - xmin. After building the list of (x,y) points for the plot, go through the xs. If the difference between any two consecutive x values is "too big", don't connect them.

Since going through the list would be expensive and rarely needed, we would only do this with a keyword, say detect_gaps. If set to True, it would default to ignoring gaps of, say, 0.01*w. Otherwise, the user could specify a proportion: detect_gaps=0.05 or whatever. This procedure would only run if detect_gaps was True or a number between 0 and 1.

Does that seem reasonable? That would fix the original problem at the cost of a little speed.

[ppurka]

comment:5

@dandrake: I guess the point of this ticket is to make the exclusion automatic. Surely, there will be a slowdown if we do this automatically for every plot; the question is what will be the magnitude of the slowdown. I suspect there won't be much slowdown with the default plot_points=200.

Your solution also sounds plausible. In case we do adopt your solution, I would like it to be a part of the exclude keyword in the plot command since that is already being used to exclude user defined points from the plot. Something like exclude='automatic' for a default proportion of 0.05 or exclude=0.02 for a user defined proportion. Hopefully, this won't break the existing code and it also won't introduce any new keywords. There are already two such keywords which deal with not plotting some points - exclude and detect_poles.

[kcrisman]

comment:6

Hmm, this idea is not bad. Particularly since we pick a default step size for the x values, one could simply choose some multiple of that (like 1 or 2) as the minimum to remove in that case. There are some issues with making sure behavior is still right for generate_plot_points in non-normal cases (no recursion, etc.) but this has some promise.

[ppurka]

comment:7

Can you check many different plots with and without this patch? It should fix two things:

1. automatically exclude invalid regions
2. fix a bug in the exclude code. This bug can be found if you try to plot the following for instance and expect the point 1.5 to be excluded.

plot(arcsec, -2, 2, exclude=[-1.5, 0, 0.5, 1.5])

[ppurka]

comment:8

Hold on.. this patch is not right. It has broken at least two plots:

plot(sin(1/x), (x, -1, 1))
plot(sin(x), (x, 0, 10), linestyle='-.')

The first one I can understand, but I don't understand the problem with the second one.

Edit: Can't reproduce this. See comment:10.

[kcrisman]

comment:9

My totally gut guess is that you're not taking adaptive recursion into account enough. Of course, one would think that would make the distance between plot points even smaller, but who knows?

I'd insert a print statement in your while loop to see what's being excluded. Also, does sin break only with the linestyle, or in general?

[ppurka]

comment:10

Actually, I am trying to reproduce it on my laptop now, and I can't reproduce the failures. I don't know what happened back then. I was manually going through each plot command in the docstrings and plotting it and verifying visually whether it was all right. For both the above examples I was getting a straight horizontal line. Now, when I am plotting it again, I don't see any anomalies. :/

I will let this patch bake for a couple more days, before I set it up for needs_review. Maybe it corrects itself automatically over time.

[ppurka]

comment:11

I manually and visually (!) tested all the plots in the doctests of sage/plot/plot.py and this patch hasn't adversely affected any of them. I also changed the test to 2 * (distance between two x-axis points).

Setting this up for review. If you know a bunch of ill-behaving functions like arcsec then please mention them here so that I can run more tests.

[jdemeyer]

comment:13

Please fill in your real name as Author.

[ppurka]

Author: Punarbasu Purkayastha

[ppurka]

Work Issues: fix doctests

[ppurka]

Attachment: trac_13246-automatic-exclusion.patch.gz

apply to devel/sage

[ppurka]

Changed work issues from fix doctests to none

[kcrisman]

comment:17

Notice this (unsurprisingly) needs a slight rebase.

Impressively, this also seems to work correctly with things like

sage: parametric_plot( (arcsec(x), x^2), (x,-2,2))

although

parametric_plot( (x,arcsec(x)), (x,-2,2))

would cause an error always - something that probably should be fixed, perhaps this is a good place to do so? Note that there is no helpful error message like with

WARNING: When plotting, failed to evaluate function at 99 points.

I'm trying to rack my brain thinking of where 2 times the difference between two consecutive plot points might cause a problem in the parametric case. Even

parametric_plot( (100000*arcsec(x), x), (x,-2,2),aspect_ratio='automatic')

doesn't seem to trigger additional exclusions. Do you see what I mean? One could imagine that in a parametric context there could be a very narrow parameter range, but a lot of room on the x-values and so everything would be excluded.

[sagetrac-git]

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

e79bdb8

fix parametric and polar plots if exclude point is outside the domain

[ppurka]

comment:58

I think I have fixed the problem with the parametric and polar plots. Needs review. :)

[kcrisman]

comment:59

Thank you for all your very hard work on this!

[ppurka]

comment:60

Replying to @kcrisman:

Thank you for all your very hard work on this!

Thanks for your extensive testing! You have uncovered many bugs in this ticket :-)

[vbraun]

comment:61

I get

sage -t --long src/sage/plot/plot.py
**********************************************************************
File "src/sage/plot/plot.py", line 1276, in sage.plot.plot.?
Failed example:
    parametric_plot((x, arcsec(x)), (x, -2, 2))
Exception raised:
    Traceback (most recent call last):
      File "/home/release/Sage/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 480, in _run
        self.execute(example, compiled, test.globs)
      File "/home/release/Sage/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 839, in execute
        exec compiled in globs
      File "<doctest sage.plot.plot.?[14]>", line 1, in <module>
        parametric_plot((x, arcsec(x)), (x, -Integer(2), Integer(2)))
      File "/home/release/Sage/local/lib/python2.7/site-packages/sage/misc/decorators.py", line 550, in wrapper
        return func(*args, **options)
      File "/home/release/Sage/local/lib/python2.7/site-packages/sage/plot/plot.py", line 1659, in parametric_plot
        return plot(funcs, *args, **kwargs)
      File "/home/release/Sage/local/lib/python2.7/site-packages/sage/misc/decorators.py", line 705, in wrapper
        return func(*args, **kwds)
      File "/home/release/Sage/local/lib/python2.7/site-packages/sage/misc/decorators.py", line 550, in wrapper
        return func(*args, **options)
      File "/home/release/Sage/local/lib/python2.7/site-packages/sage/plot/plot.py", line 1176, in plot
        G = _plot(funcs, *args, **kwds)
      File "/home/release/Sage/local/lib/python2.7/site-packages/sage/plot/plot.py", line 1383, in _plot
        newdata.append((fdata, g(x)))
      File "wrapper_rdf.pyx", line 80, in sage.ext.interpreters.wrapper_rdf.Wrapper_rdf.__call__ (build/cythonized/sage/ext/interpreters/wrapper_rdf.c:1689)
    TypeError: can't convert complex to float

[ppurka]

comment:62

This shouldn't happen at all. This is what I fixed in #16439. Does this happen if you apply this ticket after #16439?

[ppurka]

comment:63

hmm.. git is not merging the two tickets correctly. I will have to check.

[sagetrac-git]

Changed commit from e79bdb8 to 384b5a2

[sagetrac-git]

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

384b5a2

catch TypeError after #16439 is applied

[ppurka]

comment:65

@vbraun: this should fix the doctest error. This ticket by itself was working fine. Applying this after #16439 introduced a TypeError that I wasn't catching. Now it all works fine.

Edit: I think it is because of this git commit in #16439 (that is not merged into this ticket) that a TypeError started getting raised.

[ppurka]

comment:66

breaks plot in this case:

sage: var('n,a')
sage: hbar, m = 1,1
sage: psi(x,t,n)=sqrt(2/a)*sin(n*pi*x/a)*e^(-i*n^2*pi^2*hbar*t/(2*m*a^2)); 
sage: print psi
sage: Psi(x,t)=1/sqrt(2)*psi(x,t,1)+1/sqrt(2)*psi(x,t,2); 
sage: print Psi
sage: P(x,t,a) = Psi.conjugate()*Psi
sage: print P.expand()
sage: plot(P(x,1,1),x,0,1)

[ppurka]

comment:67

False alarm! It is working as intended. :-)

p = plot(P(x,1,1),x,0,1)
po = p._objects[0]
po.xdata
[0.0,0.5,1.0]

[kcrisman]

comment:68

I was happy with this, so I guess if Volker can test if the doctest error he saw is gone then we are okay.

[vbraun]

Changed branch from u/ppurka/automatic_exclusion_of_non_domain_points_in_things_like_arcsec to 384b5a2