Add floodfill

[oscarxblanco]

Dear all,
this PR adds a function called 'atfloodfill' to AT matlab and 'floodfill' to pyat in order to explore the region of acceptance from the unstable to the stable area. It has been discussed previously here #1011 .

In MATLAB it can be used in this way.

[nl, l] = atfloodfill(THERING, ...
    nturns=1000, ...
    gridsize=[10,10], ...
    axes = [1,3], ...
    window=[-10e-3,10e-3,-10e-3,10e-3], ...
    sixdoffset=findorbit6(THERING) ...
    );

and the acceptance region can be plot with the following commands

    xscale = 1e3;
    yscale = 1e3;
    figure;plot(xscale*nl(1,:),yscale*nl(2,:),'o');
    xlabel('x [mm]');ylabel('y [mm]');title('Not lost');grid;
    figure; scatter(xscale*l(1,:),yscale*l(2,:),20,log2(l(3,:)+1),'filled');
    xlabel('x [mm]');ylabel('y [mm]');title('Lost');
    a = colorbar;
    clim([0 10]);
    a.Label.String = 'log_2(turns+1)';
    colormap jet;

The boundary of the acceptance region is not provided but could be calculated from the output. Typically it could be done by filtering, interpolation over angular sectors from the mid-point of the acceptance, or with a rule similar to the game of life and joining the remaining dots. The user could decide which suits the best the needs depending on calculation time, precision and so on.

Similarly, in pyat it could be called like:

pnl_, pl_ = floodfill(ring,
                         window=(-10e-3,10e-3,-3e-3,3e-3),
                         axes=(0,2),
                         nturns=1000,
                         grid_size=(40,40),
                         )

pyat allows to parallelize the calculation with a custom tracking that chooses what particles to track next depending on the tracking results of previous particles.

[swhite2401]

Ah this is great! I don't know when I will have time to try it but I will do for sure!

[swhite2401]

@oscarxblanco I have tested it and it works very nicely!

Here is the script I used:

import at
import matplotlib.pyplot as plt
import time

ring = at.load_lattice("../lattice/hmba.mat")

t0 = time.time()
pnl, pl = at.floodfill(ring, window=(-20e-3,20e-3,-10e-3,10e-3),axes=(0,2), nturns=1000, grid_size=(40,40), use_mp=False)
t1 = time.time()
bf, sf, gf = ring.get_acceptance(['x', 'y'], [40, 40], [20.0e-3, 10.0e-3], bounds=((-1, 1), (-1, 1)), grid_mode=at.GridMode.CARTESIAN, use_mp= False)
t2 = time.time()
print(t1-t0, t2-t1)

print(pnl.shape)
print(pl.shape)

plt.plot(pnl[0,:], pnl[1, :], '.')
plt.plot(*bf)
plt.show()

Single CPU: FF: ~10s, GA: ~25s
Multi-CUP: FF: ~3s, GA: ~2s

It does present an advantage for one or few CPUs to use this method

Few differences observed wrt to the standard method, but overall the results are very consistent.

[swhite2401]

Some feedback on the implementation:
1- It would be nice to order the output such that a line can be plotted to show the boundary
2- It could be interesting to match the input of get_acceptance(), this would eventually allow to integrate this method into it using a flag (use_floodfill or grid_mode=GridMode.FLOODFILL?) instead of adding a new function (you would also profit from the re-ordering etc... already implemented)

[oscarxblanco]

Some feedback on the implementation: 1- It would be nice to order the output such that a line can be plotted to show the boundary 2- It could be interesting to match the input of get_acceptance(), this would eventually allow to integrate this method into it using a flag (use_floodfill or grid_mode=GridMode.FLOODFILL?) instead of adding a new function (you would also profit from the re-ordering etc... already implemented)

OK, I'll have a look.

[oscarxblanco]

@swhite2401 , I have modified floodfill to provide all the tracked points, a bool set to one if the particle is lost and the number of turns the particle did.

Also, floodfill is partially integrated in get_acceptance with a flag. I have to discard the number of turns in order to make it compatible with the current output. The boundary and the order of the particles is still to be defined.

Is this OK for you ? And if so, what is the method that calculates the boundary ? Floodfll tracks from unstable to stable, so, the boundary implementation probably needs to be modified.

[swhite2401]

Yes I think this is fine, you could even have created a GridMode.FLOODFILL to avoid too many test cases.

Concerning the boundary, the floodfill method returns it already, no need to compute it! The only point is to properly order the points such that you can draw a line through them. I believe this can be done by ordering them by ascending angle.

[oscarxblanco]

@swhite2401 , I have done the changes you requested. Please, let me know if there is anything I should modify.

[swhite2401]

If I plot the floodfill survived particles I still get this (in blue):

While with correct ordering I should get the orange.

[swhite2401]

In green I now run it within get_acceptance() and it looks good, nice!

[swhite2401]

All looks good to me!

[swhite2401]

The ordering in the floodfill() function could be improved in order to plot a boundary line easily but I live it up to you since this is fixed in get_acceptance().

[lfarv]

I am very much in favour of the best integration of get_acceptance and floodfill. In particular:

• using GridMode.FLOODFILL to select the method, as proposed by @swhite2401,
• using compatible names for the argument selecting the axes: floodfill uses axes, with numerical values while get_acceptance uses planes with string values. The recommended practice is to use axis or axes for 6-valued quantities (x, px…) and plane or planes for 3-valued quantities (horizontal, vertical, longitudinal). By using the auxiliary functions axis_ and plane_, the input may be numbers or strings and can be converted to the desired format,
• adding some documentation on the different methods to guide the user's choice.

@oscarxblanco and @swhite2401 : any of you interested in looking at that?

[oscarxblanco]

Dear @lfarv , sorry, I hurried and didn't wait for your comments.
Is there a way to undo the merge ? or would you prefer that I open another pull request for this ?
I think I have a few hours next week to do this.

[lfarv]

@oscarxblanco : don't worry, I reacted very late and it's easy to create a new PR. My comments about argument naming were also concerning get_acceptance and we need @swhite's advice for that.

I can add other points:

• the definition of the search range: in get_acceptance, we have the amplitudes and bounds arguments, in floodfill it's a window argument.
• for the sampling: npoints versus grid_size,
• the default number of turns is 1024 in get_acceptance and 1000 in floodfill.

Those are details but I think that consistency makes the user's life easier!

[swhite2401]

@oscarxblanco : don't worry, I reacted very late and it's easy to create a new PR. My comments about argument naming were also concerning get_acceptance and we need @swhite's advice for that.

I can add other points:

* the definition of the search range: in `get_acceptance`, we have the `amplitudes` and `bounds` arguments, in `floodfill` it's a `window` argument.

* for the sampling: `npoints` versus `grid_size`,

* the default number of turns is 1024 in `get_acceptance` and 1000 in `floodfill`.

Those are details but I think that consistency makes the user's life easier!

bounds and amplitudes are convenient to handle the different types of grids (RADIAL or CARTESIAN) to prevent mixing angles and amplitudes but it can be changed if you think a 2D array is more user friendly

npoints vs grid_size is just semantics for me, I prefer npoints and I also think it should be a positional input argument, concerning the default number of turns, I don't think this matters, I just have this habit of using powers of 2... which is really not needed in here.

Overall I would be in favor of keeping the interface of get_acceptance() unchanged (also for backward compatibility) and adapt the flood_fill option to it which does not look to complicated.

I agree that another PR is fine.