streamline plots

[jason]

Recently there was a thread on the matplotlib list where someone posted matplotlib code for a streamline plot. It looks very nice. If the code doesn't get included in matplotlib, we could still ask about including it in Sage to do streamline plots.

Mailing list message: ​http://permalink.gmane.org/gmane.comp.python.matplotlib.general/26362

Mailing list thread: ​http://comments.gmane.org/gmane.comp.python.matplotlib.general/26354

Example code: ​http://web.mit.edu/speth/Public/streamlines.py

Example plot: ​http://web.mit.edu/speth/Public/streamlines.png

[jason]

And here is another file from the same message with some fantastically nice plots:

Code: ​http://www.atm.damtp.cam.ac.uk/people/tjf37/streamplot.py

Result: ​http://www.atm.damtp.cam.ac.uk/people/tjf37/streamlines1.png

Result 2: ​http://www.atm.damtp.cam.ac.uk/people/tjf37/streamlines2.png

[jason]

I've attached both code files so that we have them even if the above links go down for whatever reason.

[jason]

See ​http://sagenb.org/home/pub/3237/ for a short example and rough Sage code (posted at ​http://sage.math.washington.edu/home/jason/sage-field-plots/)

[kcrisman]

Yo, this is now in mpl!

• ​the file
• ​the commit

I'm not sure how to figure out which release it's in, if it's in a release yet. It looks like it will probably be in 1.1.1, as 1.1.0 has been out for longer than this was in.

[jason]

See also ​http://sage.cs.drake.edu/home/pub/60/ for another example.

[jason]

I attached the example I pointed out in the sage.cs.drake.edu server.

[kcrisman]

#13693 is for the appropriate mpl upgrade for this.

[kcrisman]

Here are some ​upstream examples.

---

Here is another example from somewhere in this discussion, allowing phase plots. Not directly related but looks nice with such things.

# based on code from oddrobot, http://sagenb.org/home/pub/1532/
from sage.calculus.desolvers import desolve_system_rk4
x,y,t=var('x y t')
class DESolution:
    def __init__(self,system,time_range,initial,stepsize=0.05):
        
        self.tvar=time_range[0]
        self._times=srange(time_range[1],time_range[2],stepsize)
        self.vars=[v for v,_ in initial]
        self.dim=len(self.vars)
        self._soln=desolve_odeint(system, ics=[v for _,v in initial], times=self._times, dvars=self.vars, ivar=self.tvar)
        
    def phase_plot(self,vars=None,color='blue',**kwargs):
        # find which indices the specified variables are
        if vars is not None:
            vars_index=[self.vars.index(v) for v in vars]
        elif self.dim<=3:
            vars_index=range(self.dim)
        else:
            vars_index=range(2)
            
        p=line(self._soln[:,vars_index],color=color,**kwargs)
        # add an arrow head showing which way we are going around the phase line
        half=int(self._soln.shape[0]/2)
        p+=arrow(self._soln[half,vars_index], self._soln[half+1,vars_index],color=color)
        if len(vars_index)==2:
            p.axes_labels([str(self.vars[v]) for v in vars_index])
        return p
        
    def coordinates(self,colors=None,**kwargs):
        if colors is None:
            colors=rainbow(len(self.vars))
        p=Graphics()
        p+=sum(line(zip(self._times,self._soln[:,i]), color=colors[i], legend_label=str(self.vars[i]),**kwargs) for i in range(self.dim))
        return p

[tscrim]

I was hoping we could resurrect this ticket as I am teaching a class on ODE's this semester and would like to plot streamlines on top of the slope fields. Does someone have some time to work on this? Alas, I don't know how much time I might have to do so.

So streamlines are included in matplotlib: ​http://matplotlib.org/examples/images_contours_and_fields/streamplot_demo_features.html

There is also another plot package I came across that is open-source python that does streamlines: ​https://plot.ly/python/streamline-plots/

[paulmasson]

Since I spent too much time today reading graphics code, I can take a crack at this. It looks like a new class is needed, modeled on PlotField but calling the pyplot method streamplot instead of quiver. It shouldn't be too difficult but won't happen overnight.

Travis, if you see a way to avoid a new class, then let me know before I dive in tomorrow.

[paulmasson]

To clarify: by new class I mean an individual file like scatter_plot.py, which calls pyplot.scatter().

[kcrisman]

Yes, this is exactly how you should do it. You could even keep it in the same module as the vector field plots, you can see that e.g. contour and density plots are in the same file.

[tscrim]

Thank you for taking point on this Paul!

I agree with Karl-Dieter, you should definitely create a new class. However, contour and density plots are in separate files. Although, I have no preference whether or not you put in a separate file; it fits very well with plot_field.py.

[paulmasson]

This was surprisingly easy to get working! I copied over plot_field.py and left variable names the same for consistency. Then it was a matter of rearranging how the data arrays were assembled and changing the pyplot function name.

Doctests and useful examples have been updated. I left the _allowed_options untouched because they might be relevant.

---

New commits:

​852c719	Initial commit
​e6094a3	Initial working commit
​188b3c7	Update doctests
​4c2ff90	Add streamline_plot
​9576ba7	Update documentation

[tscrim]

Sorry to ask for a bit more. Could we add a way to input an initial point, allow input as a single function to treat it as a slope field, and an option to display the vector/slope field? If this is too much trouble, I am okay with positive reviewing this ticket as-is.

[paulmasson]

Travis, if by "initial point" you mean initial conditions, aren't those just part of the function arguments given to the plot? What more do you intend with this feature?

The idea of a convenience function that switches between streamline and vector/slope plots is interesting, but shouldn't be part of streamline_plot itself. That one should stand alone, partly because the streamplot function it uses behaves differently from the quiver function used by the others. Some of the examples from plot_field.py had to be left out of this new file because streamplot returned errors for them.

Would this convenience function be used regularly by others, or just you for one class? What's the rule of thumb for adding those sorts of features to Sage?

If you still want an additional convenience function added to this file, give me a name for it and its options.

[tscrim]

Replying to paulmasson:

Travis, if by "initial point" you mean initial conditions, aren't those just part of the function arguments given to the plot?

Yes, but how I am reading the code, it looks like it does streamlines starting at each integer point.

What more do you intend with this feature?

Plot a solution for a given IVP.

The idea of a convenience function that switches between streamline and vector/slope plots is interesting, but shouldn't be part of streamline_plot itself. That one should stand alone, partly because the streamplot function it uses behaves differently from the quiver function used by the others. Some of the examples from plot_field.py had to be left out of this new file because streamplot returned errors for them.

If you look at the slope plot, it redirects to the vector plot with special options. We should be able to automate this by looking at the input, a tuple for a vector plot and otherwise consider it a slope plot.

Would this convenience function be used regularly by others, or just you for one class? What's the rule of thumb for adding those sorts of features to Sage?

I would be surprised if this did not get significant uptake in a variety of other classes. One of my first forays into Python was plotting a single streamline for the Lorenz attractor. However, there are a number of things in Sage that are likely only used by a few people, so the rule is it gets added if someone is willing to positively review it.

If you still want an additional convenience function added to this file, give me a name for it and its options.

Granted, I would eventually like a 3d version, but that can most certainly be a followup ticket.

I should have time next week to actually directly work on this (i.e. write code).

[git]

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

​77651df	Remove unused keyword arguments
​9201066	Add density option

[paulmasson]

Travis, I went through the existing options, removed those that are not accepted by streamplot and added one that is, density. This value must be a single float or a pair of floats in a list, although I haven't added any validation for input values. It will at least allow you to see more streamlines and verify that they do not always start at integers.

There is another streamplot option I though might be interesting called start_points, but all it does is restrict streamlines to only those passing through an input list of points. You can test it by adding an option start_points=[[0,.1],[0,-.1]] or the like. Doesn't look too interesting to me.

[tscrim]

Replying to paulmasson:

Travis, I went through the existing options, removed those that are not accepted by streamplot and added one that is, density. This value must be a single float or a pair of floats in a list, although I haven't added any validation for input values. It will at least allow you to see more streamlines and verify that they do not always start at integers.

I added some input checking, but density will still need to be documented.

There is another streamplot option I though might be interesting called start_points, but all it does is restrict streamlines to only those passing through an input list of points. You can test it by adding an option start_points=[[0,.1],[0,-.1]] or the like. Doesn't look too interesting to me.

It is useful for showing numerical solutions to initial value problems (IVPs) of (nonlinear) ODEs. Also useful for seeing how stable solutions are. Although it is very annoying that it can't plot start points that are too close together:

sage: x,y = var('x,y')
sage: f = (x + y) / sqrt(x^2 + y^2)
sage: pts = [[-1., -2.+val/8] for val in range(10)]
sage: streamline_plot((cos(x), sin(x+y)), (x,-3,3), (y,-3,3), start_points=pts)
/home/travis/sage-build/local/lib/python2.7/site-packages/sage/repl/rich_output/display_manager.py:590: RichReprWarning: Exception in _rich_repr_ while displaying object: 
  RichReprWarning,
Graphics object consisting of 1 graphics primitive

Although this is not something that would (nor should) block this getting a positive review.

---

New commits:

​c443cc4	Merge branch 'u/paulmasson/10775' of git://trac.sagemath.org/sage into u/paulmasson/10775
​584486b	Adding start points option and slope field input.

[paulmasson]

Ah, so that's what you had in mind for the additions.

The documentation would not build without some fixes. Sphinx plots do not allow ^ for exponentiation, needs to be **, and variable declarations do not carry over between them. Documentation now builds.

---

New commits:

​9f3a5c2

Fix Sphinx plots

[paulmasson]

Don't know if you noticed, Travis, but in the example with start_points the streamlines come close to the start points but do not pass through them. I thought it might be because you didn't process the start points as in this example

​http://matplotlib.org/examples/images_contours_and_fields/streamplot_demo_start_points.html

but even with setting up a numpy array as indicated and transposing it, the streamlines still do not pass through the start points. In fact the resulting graph is exactly the same as not bothering to process the points, so apparently that isn't necessary.

Problem with streamplot? Thought you should know before any students ask what's going on.

[tscrim]

Thanks for fixing that Paul. I guess we should add a note and example about the behavior of start_points. Also, do you want me to document the density option?

Actually, I think I should not be normalizing the input when we want to consider it as a slope field, which makes sense when considering plotting the tangent lines. I believe it just adds a unnecessary computation that cancels when considering splitting dy / dx to (dx, dy).

[git]

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

​78c53a0	Improve documentation and example
​478b1e1	Code block for start_points

[paulmasson]

Documentation for implemented plot options now included. I also updated your language for start_points using that from the matplotlib website and added the note.

I'm including a commented code block for the official way to handle start_points in the event it is useful later.

[tscrim]

I disagree with adding start_points to the @options, where if that is changed, it is perpetuated in all future plots, which is not what we want. I will revert this back when I make my next round of changes.

Does the "official way" for start_points make any difference?

[paulmasson]

Replying to tscrim:

Does the "official way" for start_points make any difference?

Not that I can tell. xpos_array and ypos_array need to be explicit numpy arrays to avoid errors. The documentation says that start_points is also supposed to be a numpy array, so I added this block of commented code in case that is ever enforced, but you can remove it if you feel that is more appropriate.

[tscrim]

Let's use the official way. Also (not that it makes any real difference), you don't need to pop off the option:

    if 'start_points' in options:
        xstart_array, ystart_array = [], []
        for point in options['start_points']:
            xstart_array.append(point[0])
            ystart_array.append(point[1])
        options['start_points'] = numpy.array([xstart_array, ystart_array]).T

[git]

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

​79e8795

Update start_points handling

[tscrim]

Thank you for all your work on this.

[fbissey]

Not sure what Volker will find in his current run but I am failing to build the documentation in sage-on-gentoo

[plotting ] /scratch2/portage/sci-mathematics/sage-9999/work/sage-9999/src-python2_7/build/lib/sage/plot/streamline_plot.py:docstring of sage.plot.streamline_plot.streamline_plot:35: WARNING: Exception occurred in plotting streamline_plot-1
[plotting ] from /scratch2/portage/sci-mathematics/sage-9999/work/sage-9999/src-python2_7/doc/en/reference/plotting/sage/plot/streamline_plot.rst:
[plotting ] Traceback (most recent call last):
[plotting ] File "/usr/lib64/python2.7/site-packages/matplotlib/sphinxext/plot_directive.py", line 517, in run_code
[plotting ] six.exec_(code, ns)
[plotting ] File "/usr/lib64/python2.7/site-packages/matplotlib/externals/six.py", line 672, in exec_
[plotting ] exec("""exec _code_ in _globs_, _locs_""")
[plotting ] File "<string>", line 1, in <module>
[plotting ] File "<string>", line 2, in <module>
[plotting ] File "/scratch2/portage/sci-mathematics/sage-9999/work/sage-9999/src-python2_7/build/lib/sage/misc/decorators.py", line 554, in wrapper
[plotting ] return func(*args, **options)
[plotting ] File "/scratch2/portage/sci-mathematics/sage-9999/work/sage-9999/src-python2_7/build/lib/sage/plot/streamline_plot.py", line 316, in streamline_plot
[plotting ] for point in options['start_points']:
[plotting ] TypeError: 'NoneType' object is not iterable

I have this several times. I note that the bots have a similar problem in their logs.

[tscrim]

I had thought the first part of comment:32 was done by Paul. I removed start_points from the @options and did some little doc tweaking. Paul, can you make a quick check of my changes?

---

New commits:

​1e2b2d0	Merge branch 'u/paulmasson/10775' of git://trac.sagemath.org/sage into u/tscrim/streamline_plot-10775
​fa8cda9	Merge branch 'u/paulmasson/streamline_plot-10775' of git://trac.sagemath.org/sage into u/tscrim/streamline_plot-10775
​640cb6a	Some last little touchups.

[paulmasson]

Documentation now builds. Setting back to positive review.