marrink-lab / bentopy

Packs stuff in boxes
4 stars 0 forks source link

Periodic boundary conditions #17

Open ma3ke opened 8 months ago

ma3ke commented 8 months ago

Currently, bentopy does not support periodic boundary conditions. In order to judge whether this is feasible, we must answer the following questions:

  1. Can the convolution be carried out in such a way that the periodic boundaries are taken into account?
  2. Can the placement of the segments onto the background be done in a periodic manner.

The second question is not an issue at all. That is a mere matter of implementation and a placement modulo the background size. The biggest challenge here would be to be able to take into account semi-periodic boundaries (e.g., periodic over x and y axes, but not over z axis).

But the first question is where I foresee some trouble.

Experiment

This afternoon, Bart and I messed around with implementing the scipy.signal.fftconvolve function ourselves. We partially based our experimentation on a seemingly relevant SO entry. In that SO answer, it is suggested one uses np.roll in order to capture the periodicity of the convolution. Through our experimentation, we came to the conclusion (for now) that this is an accidentally correct answer, and that it does not hold over all inputs.

From there, we moved on to studying how the periodicity is treated within the scipy.signal.convolve2d function, which has a boundary="wrap" argument. To verify our this function and to check understanding, we decided to carry out a cookie-cutter experiment with this function. In such an experiment, we take a fully occupied background, and remove only the cells that reflect the exact segment that we intend to place, setting them to be available. Therefore, we know that there is exactly one correct way to insert the segment into the background.

Supposedly, convolving over the cookie-cutter background with the matching segment should return the position of the segment placement as the location of a single 0-valued cell. The position of that cell (with appropriate translation given the segment's center) should be the exact position of the cut-out segment.

From this, we can conclude that it seems that scipy.signal.convolve2d with the boundary="wrap" argument does not produce the periodic results that we expect either. This leaves us more weary of this terrain, and suggests we should try to revisit this issue at a later stage.