Plot preparation

[lfarv]

Here is the first set of plotting functions:

Plotting is split in 2 parts:

1. a generic plotting function (baseplot) which takes care of slicing the lattice into shorter elements inside the plot region, getting the data from a plot-specific function, and calling all necessary Matplotlib functions to generate the plot,
2. plot-specific functions generating the data for any kind of plot. These functions are completely independent from any graphics module, so they can be located in any at sub-package.

Only the generic function relies on Matplotlib. It is set so that Matplotlib is not required for building and installing AT. If it is unavailable, any plot attempt will generate an error at run time, but the rest of AT still runs.

By default, baseplot creates a new figure window, but if it is supplied with existing axes objects, it can be used to plot inside an existing GUI.

For a simple test, create a Lattice object and run ring.plot_beta().

Plotting still misses the synoptic of lattice elements along the s-axis, this will be the next step.
Let me know your comments on the organisation, naming, etc. of the various functions.

Note: there is still a bug in the Lattice slicing function which gives wrong results for the err.mat lattice, I am working on that.

[lfarv]

There are still problems with Matplotlib 2.x and consequently python2.7. Please test plots on python3 and Matplotlib 3.x

[lfarv]

This version now works for:

• python 2.7 and Matplotlib 2.2.4
• python 3.6 and Matplotlib 3.0.3

[lfarv]

Everything is now working for plotting optical functions, trajectory or for building new plots. To test it, install:
Matplotlib >= 2.1.2 for python 2.7
Matplotlib >= 3.0.0 for python 3.x

Then create a Lattice object ring and run:

>> ring.plot_beta()
>> ring.plot_trajectory(r_in)

and give me your comments

[T-Nicholls]

My comments on the code and my use of it so far.

• I am very happy with the plotting groundwork and it mostly functions well for me.
• We only support the plotting of up to two different 'keys' at a time for plot_linear, one for left axis and one for the right, this is fine and makes sense; however, I don't think it is explicitly mentioned in the docstrings that this is the case, and when I entered three no error was raised the third one was simply not plotted.
• I think the indexing for multi-dimensional arrays has some issues, for example:

>>> at.plot.plot_linear(lattice, ('m44', [0, 1], [1, 2]))  # works as expected
[(0, 1), (1, 2)]
>>> at.plot.plot_linear(lattice, ('m44', [0, 1], [0, 2]))  # doesn't
[(0, 0), (1, 2)]

• Also for plot_linear:
  • in the docstring, surely 2, slice(4) corresponds to T31, T32, T33, T34?
  • they are plotted as one colour on the figure and the legend, for example, says T_1234,3; I think they should be plotted as separate lines?

[T-Nicholls]

I will ask Will about what our plotting needs are. However, bearing in mind that I am not the intended user, my vision for pyAT's plotting functionality was that most physics data would be plottable through general methods on the lattice object like plot_linear and more specific functions like pltdata_beta_disp would be an example for users to base their own plotting functions off of. In summary, I think most of the physics data we calculate should be able to be plotted through a few lattice methods and the more specific plotting functions should be the example in the plotting package. Does this align with your intentions for plotting?

[lfarv]

From my experience, I think that we should have:

• very simple methods for plotting the most usual data. For instance, plot_beta will represent more than 90% of the plots: it's how you see immediately what the lattice looks like. We probably need a few other ones like plot_phase_advances and plot_closed_orbit.
• a few generic methods for plotting more exotic data from line_pass, linopt and ohmi_envelope. plot_linear and plot_trajectory are such ones, I am thinking of plot_emittance. But these will be rarely used…
• everything else will be done with user-written functions, and there is an unlimited potential for that!

Of course, all suggestions are welcome!

[lfarv]

Rebased on the updated master branch, and added the get_elements method in the Lattice object.

Here in is also a solution for the "get_ring_energy" factorisation:

• I moved the full_scan iterator into a static method _scanner of Lattice (set private because its use is tricky)
• Using that, I create a very simple get_ring_properties function, more general since it returns a dictionary with at least name, energy and periodicity keys. So get_ring_energy(ring) can be replaced by get_ring_properties(ring)['energy'], and could then be removed (or not…)

I still have to tune the search algorithm.
But I have a serious problem of cyclic imports: in radiation.py, I cannot use
from at.lattice import get_properties, but instead import at and then use at.get_properties. This may indicate a wrong module organisation…

[lfarv]

For me this is now ready for merging.

[lfarv]

Last minute: fix for Nan in emittance results reported in #104

[T-Nicholls]

@lfarv I discussed the problem of cyclical imports with Will and we think that get_ring_properties should be in lattice.utils as it works on any ring, and so _scanner should also be moved there. This would probably require a few changes, but should resolve the import problem as both lattice_object and radiation should import from lattice.utils. Side note I don't think get_ring_energy is needed and agree it could be removed.

[T-Nicholls]

This doesn't fix the warning, we need to do the same for emit3:

emit3sq = numpy.array([det(r66[s, s]) for s in _submat])
# Prevent from unrealistic negative values of the determinant
emit3 = numpy.sqrt(numpy.maximum(emit3sq, 0.0))

If I understand the code correctly, I don't think we need to do it for emit2 at all because accidentally negative values only occur for uncoupled machines, and now emit3[1] is not nan it should be caught by the eilf on line 97, please correct me if I am wrong.

[lfarv]

• You are right, the problem occurs only for uncoupled machines only, otherwise the vertical emittance is significantly positive.
• I do not think that the problem may occurs on emit3 which comes directly from r66 (at least this is not what happens with @carmignani ). It may come from the double matrix inversion at lines 98 and 100.
  It's difficult to say since I cannot reproduce the problem, but I would keep the test as it is before adding another one on emit3, which is still possible later, but only if it is necessary.

[T-Nicholls]

In @carmignani's test output you can see two warnings are raised:

test/test_physics.py::test_ohmi_envelope[refpts0]
test/test_physics.py::test_ohmi_envelope[refpts1]
  /mntdirect/_machfs/carmigna/python/at3.4/pyat/at/physics/radiation.py:87: RuntimeWarning: invalid value encountered in sqrt
    emit3 = numpy.sqrt(numpy.array([det(r66[s, s]) for s in _submat]))

test/test_physics.py::test_ohmi_envelope[refpts0]
test/test_physics.py::test_ohmi_envelope[refpts1]
  /mntdirect/_machfs/carmigna/python/at3.4/pyat/at/physics/radiation.py:100: RuntimeWarning: invalid value encountered in sqrt
    [det(r44[s, s], check_finite=False) for s in _submat[:2]]))

The first is for emit3 the second is for emit2. I also had two warnings in my environment, after applying your fix I only had one, the emit3 warning. The point I was trying to make above was that I thought that the only reason emit2 was warning was because of the failure of emit3. I developed this theory after applying the fix to the emit3 calculation but not the emit2 calculation and both warnings disappeared for me. Perhaps I was a little overzealous in saying the check for emit2 could be removed, but I definitely think that we need to add a check for emit3.

[T-Nicholls]

However, I do not know which is the cause of the test failures as I am unable to run the tests in my environment that raises warnings. So patching just emit2 my fix the tests but it certainly doesn't remove both warnings.

[lfarv]

You are right, I missed the warnings. But the 2 errors come from emit2.
emit2 and emit3 are two independent derivations from r66. Fixing one has no effect on the other. So I will add the fix on emit3 also.

A word on the meaning of all this:
r66 is a 6x6 symmetrical matrix which is the mathematical representation of a 6D ellipsoid with axes x, x', z, z', l, delta. To help in visualising this thing, we project the ellipsoid on 3 planes: x-x', z-z', l-delta. This gives 3 ellipses and their surfaces are the 3 values in emit3. For emit2, we first cut the ellipsoid along l=0 and delta=0. We are left with a 4D ellipsoid (r44), which is projected on x-x' and z-z', giving 2 ellipses (contained in the ones of emit3, of course), and their surfaces in emit2.
In the case of an uncoupled machine, the ellipsoid is "flat", its volume is zero. When projected, we get on some planes (z-z' in our example) a segment of a line instead of an ellipse (an ellipse with surface 0). Because of rounding errors, this is even "thinner" than a line and the surface is imaginary. But for me, 0 or Nan mean the same: no vertical emittance. Clipping the emittance to zero is the best fix we can apply.

What worries me is the fact that I couldn't reproduce the bug on any machine, MacOS, Windows of Linux…

[T-Nicholls]

Ahh ok, thank you for the explanation.
I agree, I found it very difficult to pin down even though I was able to replicate it in one environment.

[T-Nicholls]

@lfarv I can now repeatably get warnings in any virtual environment by pip uninstalling numpy and scipy and then adding the locations of our Diamond versions to the python path. Will says that the only difference between them is that pip install uses wheels but the Diamond versions are source installed and they should be otherwise identical. Another thing worth noting is that hmba.mat, and so the tests, work fine for me but Diamond's diad.mat lattice raised warnings and had nan for emitXY[1]. Everything about this makes it look like it is going to be very hard to pin down a single root cause for the numerical errors. I think it is best to just settle with having fixed the problem on our end, that said I would be happy to test any theories you might have, now I can more easily replicate the conditions.

[lfarv]

Now for the cyclic imports, I agree that logically get_ring_properties should belong to lattice.utils, but then lattice.utils will need to import Lattice, which introduces another cycle… Unless we replace testing isinstance(ring, Lattice) bytry: energy=ring.energy.... I had the import problem long before get_ring_properties was introduced, but until now I could manage by modifying the order of imports here and there. There should be a way to organise the modules from low-level to high level but this would need a deep reorganisation. And I think that there is a real cycle since Lattice depends on physics to implement many of its methods, and physics may depend on the properties of a Lattice.
I propose that we look at that in a separate branch, and that we terminate this pull request as soon as:

• the nan in emittance is fixed (at least for emit2)
• I have removed get_ring_energy

The present situation should not be a problem for users since they should reach everything directly from the at package and not from sub-packages.

[T-Nicholls]

@lfarv ok I am happy to postpone the rearrange until we have thought more deeply on the topic.

[lfarv]

Last modification? Tell me!

[T-Nicholls]

@lfarv as far as I'm concerned yes, now both warnings are fixed I am happy.

[lfarv]

Any remarks on this pull request before merging?

[willrogers]

There's a lot here, I haven't looked in a lot of detail at the plotting itself, but it generally looks sensible.

The only thing I would say is that we could use a readme (e.g. to specify the versions of Matplotlib and just tell me how to do a basic plot) either in the plot package or added to the pyat readme.

We should start reporting test coverage as well, but that's not to do with this pull request particularly.

[T-Nicholls]

@willrogers I added coverage to the test calls in the appveyor and travis .yml files on my wheelBuilding branch.