<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Facetted subgrids
We can split the image (facetting), and we can split the grid (subgrids) and do a lot of operations separately. This works out relatively straightforwardly. However, can we combine the two operations? As it turns out the answer is "yes", but the implementation has a few subtleties to it, so this notebooks walks through it step by step.
Step1:  Problem definition
To enable fast imaging we are interested in various types of locality. Specifically for images we care about locality in image space, whereas for dealing with visibilities we need locality in grid space (sub-grids). For example, we might have nodes storing pieces (facets) of the image, and now want to reconstruct a certain portion of the grid for de-gridding. Or the other way around - having gridded visibilities we want to propagate its contributions to facet nodes.
Mathematically speaking, facets and sub-grids are low-passes in time and frequency spaces. Let us express this mathematically as one-dimensional functions, and let us view everything from grid space. Let $G$ be the full grid data, assumed to be bounded / repeating past $\frac 12$ and sampled at a rate $N$
Step2:  Now let $A_i$, $B_j$ be grid and facet masks respectively, so $A_i(x) = 0$ iff $\left|x-x_i\right| > x_A$, $\sum_i A_i = 1$
and $\mathcal F_y B(y) = 0$ iff $\left|y - y_i\right| > y_B$, $\sum_i \mathcal F A_i = 1$.  Then we define
Step3:  Note that to cover the entire number space with functions of finite support requires an infinite number of sub-grids and facets. However, in practice this does not really matter, as due to finite sampling we are effectively looking at infinite (repating) equivalence classes of facet/sub-grids anyway. Due to linearity this does not change the argument, so we will use the simpler representation.
The recombination problem
Clearly we have
Step4:  After all, $A_i (B_j \ast G) \ne B_j \ast A_i G$. This has direct impact on how compact we can represent these
Step5:  That is pretty close to what we are looking for
Step6:  Close enough, this precision is more than enough for our purposes! However, we are not quite finished yet.
Correction for convolution
Our original goal was to construct $B_j \ast A_i G$. At this point we have a good approximation for $n_j \ast A_i G$. So it might look like all we need is the truncated inverse of the wave function $b_j$ such that $b_j \ast n_j = B_j$
Step7:  The function $b$ looks pretty wild in grid space, with a lot of high frequency. But that is to be expected to some degree, after all it "undoes" the well-behavedness of $n$. The second diagram shows that we indeed have $b_j\ast n_j = B_j$ as intended.
Now let us attempt to go after the grand price, $b_j \ast m_i (n_j \ast a_i G)$
Step8:  Looks okay-ish at first glance. However, let us measure the errors objectively
Step9:  A lot worse than we had hoped!
It is not too hard to trace the source of the errors back
Step10:  Note that the base accuracy performance of the wave function has improved a bit simply because we have increased $y_n$ slightly. However, the important part is that we can make much better use of the results
Step11:  Reducing the grid-size / image sampling rate has clearly not impacted the precision too much. Note that while the $(1-m)(n*AG)$ error term should be zero in the $x_m$ grid region, with the reduced size the errors get aliased back onto the subgrid. For our parameters this happens to cause an error "well" in the centre of the sub-grid.
However, limiting in grid space is only half the story. As established at the start we not only want to reduce our grid size down to $x_m$, but want to also get rid of image data beyond $y_n$
Step12:  As intended, $n_j\ast a_iG$ drops to zero past $y_N$. This is clearly true even after downsampling, so we can truncate the signal in image space as well. This does not impact signal reconstruction in any visible way.
Step13:  The errors have not changed too much even though we have now significantly reduced the number of samples in both grid and image space. The downsampling operation eliminated the incidental "well" in grid space, but $n_j\ast a_iG \approx m_i (n_j \ast a_i G)$ holds just as it did before. This is what we were after!
For the final step, let's convolve with $b_j$. Notice that
Step14:  And we made it out the other end. Seems the approximation is indeed robust against truncation and down-sampling as intended.
Sub-grids recombination
So we can obtain facet data from sub-grids effectively. But does this also work in the "opposite" direction, if we want to construct a sub-grid from facets?
\begin{align}
S_i = A_iG &= A_i \left( \left( \sum_j B_j \right) \ast  G \right) = \sum_j A_i \left( B_j \ast  G \right) \
\end{align}
Let us split $A_i$ and $B_j$ as we did before. Now what we would like to establish is that
Step15:  This is entirely dual to the original method -- we are doing everything backwards. Now $G$ needs to get convolved with $b_j$ right away, and we end with the $A_i$ mask cutting away a bit of grid space.
It is instructive to have closer look at the function of this mask operation, especially considering the reduced case
Step16:  Reproduction of the original signal quickly detoriates beyond the $x_A$ region. Especially note how $n\ast mG'$ tends to jump around violently beyond $x_A$. This is even worse for the reduced case, where the error aliases around the $x_m$ edges.
Data sizes and scaling
Let us have a closer look and check exactly how effective this has been. If $A (B \ast G) = B (A \ast G)$ was true, we would expect that we would be able to represent the dependencies between a facet and a subgrid with about $\lceil 4x_Ay_B \rceil$ samples.
Our method only allows truncating to $x_m$ in grid and $y_n$ in image space, therefore we need $\lceil 4x_my_n \rceil \approx \lceil 4(x_A+x_n)y_n \rceil$ samples
Step17:  As expected, modulo some rounding of pixels at the edges. We are not actually that far away from the limits given by information theory. However, this is clearly for a very small data size, so what would happen if we started doing this for larger grids? Would we be able to retain the reached accuracy and efficiency?
To do this quickly, let us attempt to predict how much error we expect from our method. As we have seen before, the error of
$$b_j \ast m_i (n_j \ast a_i G)$$
consists of
Step18:  This is clearly quite pessimistic even when comparing against the worst case, but it gives us a decent starting point.
Let us assume that we increase the facet size ($y_B$). To make sure that we compare like with like, let us attempt to keep $\mathcal F n_j(y_B)$ constant. This means that we need
Step19:  At the same time, the "overhead" of our solution decreases slowly depending on the sub-grid size $x_A$, eventually approximating $y_N/y_B$. This is of course because the constant grid margin eventually stops mattering
Step20:  So how much accuracy can we get for a certain overhead, optimally? If we settle on a certain overhead value $o$, we can calculate $y_n$ from a value of $x_ny_n$
Step21:  Facet Split
So far we have been rather cavalier with using arrays larger than their "true" size. This simplifies things mathematically, yet we must remember that we don't actually have the $b_j \ast F_j$ available. Instead, we will handle an FFT grid, which is actually sampled at a finite rate - which we can represent as multiplying by a comb function
Step22:  The problem here is that $\mathcal Fm$ never approaches zero - it is a sinc function that doesn't fall off very much at all. Therefore we have $\mathcal F m_i(y) \ne \left\operatorname{III}_{2y_P} \ast \mathcal F m_i \right$ no matter how much we restrict the $y$ region.
However, we can actually choose $m$
Step23:  However, now this also limits $m$ to $y_N$ in image space
Step24:  Good starting point. We just need to figure out how far we can reduce the image size (here $y_P$) without impacting our ability to reconstruct the image. What we need is
Step25:  This shows that by using the low-pass-filtered $m'$ we can indeed perform the entire computation using just $2y_P$ sample points
Step26:  Note that we are only truncating terms here that we already know to fall to zero. However, the choice of $y_P$ is still quite subtle here, as we need padding to make sure that the convolutions work out in the correct way.
Actual implementation
So to review, what we have identified are two approximations
Step27:  So let us start by reconstructing sub-grids from facets. First step is to convolve in $b_j$, derived from the PSWF. This is a cheap multiplication in image space. We then pad to $2y_P$ which yields us
$$\operatorname{III}_{2y_P} \ast  \mathcal Fb_j \mathcal FG$$
in FFT convention.
Step28:  Next we need to cut out the appropriate sub-grid for $m_i$ for the subgrid we want to construct (here
Step29:  Note that we would clearly not construct them separately for a real pipeline, as they are simply shifted. Due to the truncation in frequency space these are not quite top-hat functions any more. These terms now get used to extract the sub-grid data from each facet
Step30:  Next step is to multiply in $n_j$ in order to un-do the effects of $b$ and cut out the garbage between $y_N$ and $y_P$. This means we arrive at
Step31:  Quick mid-point accuracy check against the approximation formula using full resultion. We should be looking at only rounding errors.
Furthermore, when padded to the full density the Fourier transform should fall to the error level provided by the PSWF past $x_N+x_M$ - indicating that we can reduce the sampling rate / grid size without losing information
Step32:  So our final step is to reduce the sampling rate (sub-grid size) to $x_M$. This is the step where we actually introduce the bulk of our error, as the "tail" regions outside $x_M$ get aliased in. As established before, this especially copies the $x_M+x_N$ region inside, which doesn't hurt because we are only interested in the centre $x_A$ part.
As long as $x_M$ divides the grid well, this can simply be achieved by selecting every $\frac1{x_M}$th sample
Step33:  At this point, all that is left is to put together the sum
$$S_i = a_i \sum_j \left( n_j \ast m_i (b_j \ast F_j)\right)$$
which eliminates $b_j$. Note that we know $A$ to be a pure mask, therefore this is simply truncation in grid space.
Step34:  Note the pattern of the errors in image space - this is the position dependence of the accuracy pattern from the $b_j$ multiplication. As we stitch together 5 facets, this pattern repeats 5 times.
However, note how the "stitching" works out really well - accuracy is especially good around the $y_B$ areas where the data from two facets "overlaps". This is likely because are effectively averaging over two samples here, which boost our accuracy.
Interactive playground
Test different parameters
Step35:  2D case
Clearly this can be generalised to arbitrarily high dimensions. However, note that this makes things worse in a few ways at the same time
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
from matplotlib import pylab
import matplotlib.patches as patches
import matplotlib.path as path
from ipywidgets import interact
import numpy
import sys
import random
import itertools
import time
import scipy.special
import math
pylab.rcParams['figure.figsize'] = 16, 10
pylab.rcParams['image.cmap'] = 'viridis'
try:
    sys.path.append('../..')
    from cracodile.synthesis import *
    from util.visualize import *
    print("Crocodile mode")
except ImportError:
    print("Stand-alone mode")
    # Convolution and FFT helpers
    def conv(a, b): return ifft(fft(a) * fft(b))
    def coordinates(N):
        return numpy.fft.fftshift(numpy.fft.fftfreq(N))
    def fft(a):
        if len(a.shape) == 1: return numpy.fft.fftshift(numpy.fft.fft(numpy.fft.ifftshift(a)))
        elif len(a.shape) == 2: return numpy.fft.fftshift(numpy.fft.fft2(numpy.fft.ifftshift(a)))
    def ifft(a):
        if len(a.shape) == 1: return numpy.fft.fftshift(numpy.fft.ifft(numpy.fft.ifftshift(a)))
        elif len(a.shape) == 2: return numpy.fft.fftshift(numpy.fft.ifft2(numpy.fft.ifftshift(a)))
    def pad_mid(a, N):
        N0 = a.shape[0]
        assert N >= N0
        return numpy.pad(a, len(a.shape) * [(N//2-N0//2, (N+1)//2-(N0+1)//2)], mode='constant', constant_values=0.0)
    def extract_mid(a, N):
        assert N <= a.shape[0]
        cx = a.shape[0] // 2
        s = N // 2
        if N % 2 == 0:
            return a[len(a.shape) * [slice(cx - s, cx + s)]]
        else:
            return a[len(a.shape) * [slice(cx - s, cx + s + 1)]]
    def anti_aliasing_function(shape, m, c):
        if len(numpy.array(shape).shape) == 0:
            mult = 2 - 1/shape/4
            return scipy.special.pro_ang1(m, m, c, mult*coordinates(shape))[0]
        return numpy.outer(anti_aliasing_function(shape[0], m, c),
                           anti_aliasing_function(shape[1], m, c))
    def coordinates2(N):
        N2 = N // 2
        if N % 2 == 0:
            return numpy.mgrid[-N2:N2, -N2:N2][::-1] / N
        else:
            return numpy.mgrid[-N2:N2+1, -N2:N2+1][::-1] / N
    def _show(a, name, scale, axes):
        size = a.shape[0]
        if size % 2 == 0:
            low,high = -0.5, 0.5 * (size - 2) / size
        else:
            low,high = -0.5 * (size - 1) / size, 0.5 * (size - 1) / size
        low = (low - 1/size/2) * scale
        high = (high - 1/size/2) * scale
        cax=axes.imshow(a, extent=(low,high,low,high)); axes.set_title(name);
        axes.figure.colorbar(cax,shrink=.4,pad=0.025)
    def show_grid(grid, name, theta, axes):
        return _show(grid, name, theta, axes)
    def show_image(img, name, theta, axes):
        return _show(img, name, img.shape[0] / theta, axes)
# Helper for marking ranges in a graph
def mark_range(lbl, x0, x1, y0=None, y1=None, ax=None):
    if ax is None: ax = pylab.gca()
    if y0 is None: y0 = ax.get_ylim()[1]
    if y1 is None: y1 = ax.get_ylim()[0]
    wdt = ax.get_xlim()[1] - ax.get_xlim()[0]
    ax.add_patch(patches.PathPatch(patches.Path([(x0,y0), (x0,y1)]), linestyle="dashed"))
    ax.add_patch(patches.PathPatch(patches.Path([(x1,y0), (x1,y1)]), linestyle="dashed"))
    if pylab.gca().get_yscale() == 'linear':
        lbl_y = (y0*7+y1) / 8
    else: # Some type of log scale
        lbl_y = (y0**7*y1)**(1/8)
    ax.annotate(lbl, (x1+wdt/200, lbl_y))
Explanation: Facetted subgrids
We can split the image (facetting), and we can split the grid (subgrids) and do a lot of operations separately. This works out relatively straightforwardly. However, can we combine the two operations? As it turns out the answer is "yes", but the implementation has a few subtleties to it, so this notebooks walks through it step by step.
End of explanation
N = 2000
x = coordinates(N); fx = N * x
ticks = coordinates(10); fticks = N * coordinates(12)
G = numpy.random.rand(N)-0.5
pylab.rcParams['figure.figsize'] = 16, 4
pylab.plot(x, G); pylab.legend(["G"]); pylab.xticks(ticks); pylab.show()
pylab.plot(fx, fft(G).real); pylab.legend(["F[G]"]); pylab.xticks(fticks); pylab.show()
Explanation: Problem definition
To enable fast imaging we are interested in various types of locality. Specifically for images we care about locality in image space, whereas for dealing with visibilities we need locality in grid space (sub-grids). For example, we might have nodes storing pieces (facets) of the image, and now want to reconstruct a certain portion of the grid for de-gridding. Or the other way around - having gridded visibilities we want to propagate its contributions to facet nodes.
Mathematically speaking, facets and sub-grids are low-passes in time and frequency spaces. Let us express this mathematically as one-dimensional functions, and let us view everything from grid space. Let $G$ be the full grid data, assumed to be bounded / repeating past $\frac 12$ and sampled at a rate $N$:
End of explanation
x0 = 0.2; xA = 200/N
yB = N/10
A = numpy.where(numpy.abs(x-x0) <= xA, 1, 0)
B = ifft(numpy.where(numpy.abs(fx) <= yB, 1, 0))
pylab.plot(x, A, x, B.real, x, A*G); pylab.legend(["A", "B", "AG"]);
pylab.xticks(ticks); mark_range("xA", x0-xA,x0+xA); pylab.show();
pylab.plot(fx, fft(A).real, fx, fft(B).real,
           fx, fft(conv(B, G)).real); pylab.legend(["F[A]", "F[B]", "F[B*G]"]);
pylab.ylim(-20,30); pylab.xticks(fticks); mark_range("yB",-yB,yB); pylab.show();
Explanation: Now let $A_i$, $B_j$ be grid and facet masks respectively, so $A_i(x) = 0$ iff $\left|x-x_i\right| > x_A$, $\sum_i A_i = 1$
and $\mathcal F_y B(y) = 0$ iff $\left|y - y_i\right| > y_B$, $\sum_i \mathcal F A_i = 1$.  Then we define:
\begin{align}
S_i &= A_iG && \text{sub-grid} \
F_j &= B_j \ast G && \text{facet}
\end{align}
End of explanation
pylab.plot(x, conv(B, A*G).real); pylab.legend(["B*AG"]);
pylab.xticks(ticks); mark_range("xA", x0-xA,x0+xA); pylab.show();
pylab.plot(fx, fft(A*conv(B,G)).real); pylab.legend(["F[A(B*G)]"]);
pylab.xticks(fticks); mark_range("yB", -yB,yB); pylab.show();
Explanation: Note that to cover the entire number space with functions of finite support requires an infinite number of sub-grids and facets. However, in practice this does not really matter, as due to finite sampling we are effectively looking at infinite (repating) equivalence classes of facet/sub-grids anyway. Due to linearity this does not change the argument, so we will use the simpler representation.
The recombination problem
Clearly we have:
\begin{align}
S_i = A_iG &= A_i \left( \left( \sum_j B_j \right) \ast  G \right)  \
&= A_i \sum_j \left( B_j \ast  G \right) =\sum_j A_i F_j \
F_j = B_j \ast G &= B_j \ast \left( \left( \sum_i A_i \right) G \right)  \
&= B_j \ast \left( \sum_i  A_i   G  \right) = \sum_i \left( B_j \ast S_i \right)
\end{align}
Because of linearity properties. However, this is not particular efficient for re-distributing data. While we know that $B_j \ast G$ is limited in image space, this is not true for $A_i (B_j \ast G)$ -- and on the other side of the coin, while $A_i G$ is limited in grid space, this is not the case for $B_j \ast A_i G$:
End of explanation
alpha = 0; xN = 25 / N; yN = yB
n = ifft(pad_mid(anti_aliasing_function(int(yB*2), alpha, 2*numpy.pi*yN*xN), N)).real
pylab.semilogy(x, numpy.abs(n)); pylab.legend(["n"]);
pylab.xticks(ticks); mark_range("$x_n$", -xN,xN); pylab.show();
pylab.semilogy(fx, numpy.abs(fft(n))); pylab.legend(["F[n]"]);
pylab.xticks(fticks); mark_range("$y_n=y_B$", -yN,yN); pylab.show();
Explanation: After all, $A_i (B_j \ast G) \ne B_j \ast A_i G$. This has direct impact on how compact we can represent these: As we only know $A_i F_j$ to be limited in grid space, we cannot reduce the sampling rate, and as $B_j \ast S_i$ is only limited in image space we cannot reduce the grid size.
Therefore in order to reconstruct $S_i$ from $k$ facets we need to collect $k$ grids the same size as $S_i$! This might not sound too bad, but remember that to reconstruct all $l$ sub-grids we have $l \times k$ dependencies between facets and sub-grids to cover. So if we can only reduce the amount to transfer by a factor of $O(k)$ or $O(l)$, that still means that the amount of data to exchange scales linearly with either the number of facets or the number of sub-grids.
Facet recombination
What can we do? The masks $A_i$ and $B_j$ we have chosen leave us with little wiggle room, so let us split them into sub-masks $A_i=a_im_i$ and $B_j = b_j\ast n_j$ with $a, m$ and $b, n$ being bounded in grid and image space respectively. Now let us attempt to establish:
$$B_j \ast A_iG = b_j \ast n_j \ast m_i a_i G \approx b_j \ast m_i ( n_j \ast a_i G)$$
If we ignore $b_j$ for a moment this means that we want:
$$n_j \ast m_i a_i G \approx m_i (n_j \ast a_i G)$$
Let us assume that $m_i$ is a top-hat function mask and has a greater support than $A_i$, so $a_i m_i = a_i = A_i$. Then it follows:
\begin{align}
\Leftrightarrow n_j \ast a_i' G &\approx m_i (n_j \ast a_i G) \
\Leftrightarrow (m_i - 1) (n_j \ast a_i G) &\approx 0
\end{align}
By choosing $m$ as a top hat function, we now have that $m_i - 1 = 0$ for all $|x-x_0| < x_m$. This means we need:
\begin{align}
\Leftarrow \forall |x-x_0| \ge x_m: n_j \ast a_i G &\approx 0
\end{align}
We know that $a_iG$ is truncated at $x_A$, therefore what we need is some $x_n$ such that $n_j(x)\approx 0$ for all $|x-x_0| < x_n$. Note that we still assume $n_j$ to be truncated in image space, so it is impossible for $n_j$ to fall to zero entirely. However, if we get close enough to zero we should be able to make the original approximation work as long as $x_M \ge x_A + x_n$.
The convolution
So at this point everything hinges on us identifying a function $n$ that has finite support in image (frequency) space and falls very close to zero in grid (time) space. There are a number of options here, one of the ad-hoc choices being prolate spheroidal wave functions. Let us construct one for $y_n = y_B$:
End of explanation
xM = xA + xN
m = numpy.where(numpy.abs(x-x0) <= xM, 1, 0)
ideal = conv(n, A*G)
approx = m * conv(n, A*G)
error = (1-m) * conv(n, A*G)
pylab.plot(x, ideal.real, x, approx.real); pylab.legend(["n*AG", "m(n*AG)"]);
pylab.xticks(ticks); mark_range("$x_A$", x0-xA,x0+xA); mark_range("$x_m$", x0-xM,x0+xM); pylab.show();
pylab.semilogy(x, numpy.abs(error)); pylab.legend(["(1-m)(n*AG)"]); 
pylab.xticks(ticks); mark_range("$x_A$", x0-xA,x0+xA); mark_range("$x_m$", x0-xM,x0+xM); pylab.show();
pylab.semilogy(fx, numpy.abs(fft(error))); pylab.legend(["F[(1-m)(n*AG)]"]);
pylab.xticks(fticks); mark_range("$y_n=y_B$", -yN,yN); pylab.show();
print("RMSE:", numpy.sqrt(numpy.mean(numpy.abs(error)**2)), "(image:",numpy.sqrt(numpy.mean(numpy.abs(fft(error))**2)),")")
Explanation: That is pretty close to what we are looking for: In image space the function falls almost completely to zero, while in grid space we are left with relatively low errors. We can reduce those further by increasing either $x_n$ or the sampling rate.
Let us evaluate the error terms for this function:
End of explanation
b = ifft(pad_mid(1/anti_aliasing_function(int(yB*2), alpha, 2*numpy.pi*yN*xN), N)).real
pylab.plot(x, b); pylab.legend(["b"]); pylab.xticks(ticks); pylab.show();
pylab.semilogy(fx, numpy.abs(fft(b)), fx, numpy.abs(fft(conv(b,n)))); pylab.legend(["F[b]", "F[b*n]"]);
pylab.xticks(fticks); mark_range("$y_n=y_B$", -yN,yN); pylab.show();
Explanation: Close enough, this precision is more than enough for our purposes! However, we are not quite finished yet.
Correction for convolution
Our original goal was to construct $B_j \ast A_i G$. At this point we have a good approximation for $n_j \ast A_i G$. So it might look like all we need is the truncated inverse of the wave function $b_j$ such that $b_j \ast n_j = B_j$:
End of explanation
ideal = conv(B, A*G)
approx = conv(b, m * conv(n, A*G))
error = approx - ideal
pylab.plot(x,ideal.real, x, approx.real); pylab.legend(["B*AG", "b*m(n*AG)"]);
pylab.xticks(ticks); mark_range("$x_A$", x0-xA,x0+xA); mark_range("$x_m$", x0-xM,x0+xM); pylab.show();
Explanation: The function $b$ looks pretty wild in grid space, with a lot of high frequency. But that is to be expected to some degree, after all it "undoes" the well-behavedness of $n$. The second diagram shows that we indeed have $b_j\ast n_j = B_j$ as intended.
Now let us attempt to go after the grand price, $b_j \ast m_i (n_j \ast a_i G)$:
End of explanation
pylab.semilogy(x, numpy.abs(error)); pylab.legend(["b*(1-m)(n*AG)"]);
pylab.xticks(ticks); mark_range("$x_A$", x0-xA,x0+xA); mark_range("$x_m$", x0-xM,x0+xM); pylab.show();
pylab.semilogy(fx, numpy.abs(fft(error))); pylab.legend(["F[b*(1-m)(n*AG)]"]); 
pylab.xticks(fticks); mark_range("$y_n=y_B$", -yN,yN); pylab.show();
print("RMSE:", numpy.sqrt(numpy.mean(numpy.abs(error)**2)), "(image:",numpy.sqrt(numpy.mean(numpy.abs(fft(error))**2)),")")
Explanation: Looks okay-ish at first glance. However, let us measure the errors objectively:
End of explanation
alpha = 0; yN = yB + N * 0.02
pswf = anti_aliasing_function(int(yN*2), alpha, 2*numpy.pi*yN*xN)
n = ifft(pad_mid(pswf, N)).real
b = ifft(pad_mid(1/extract_mid(pswf, int(yB*2)), N)).real
error1 = (1-m) * conv(n, A*G)
error2 = conv(b, (1-m) * conv(n, A*G))
pylab.semilogy(fx, numpy.abs(fft(n)), fx, numpy.abs(fft(b)), fx, numpy.abs(fft(conv(b,n))));
pylab.legend(["F[n]","F[b]", "F[b*n]"]);
pylab.xticks(fticks); mark_range("$y_B$", -yB,yB); mark_range("$y_n$", -yN,yN); pylab.show();
pylab.semilogy(fx, numpy.abs(fft(error1)), fx, numpy.abs(fft(error2))); pylab.legend(["F[(1-m)(n*AG)]", "F[b*(1-m)(n*AG)]"]);
pylab.xticks(fticks); mark_range("$y_B$", -yB,yB); mark_range("$y_n$", -yN,yN); pylab.show();
print("RMSE (w/o b):", numpy.sqrt(numpy.mean(numpy.abs(error1)**2)), "(image:",numpy.sqrt(numpy.mean(numpy.abs(fft(error1))**2)),")")
print("RMSE:", numpy.sqrt(numpy.mean(numpy.abs(error2)**2)), "(image:",numpy.sqrt(numpy.mean(numpy.abs(fft(error2))**2)),")")
Explanation: A lot worse than we had hoped!
It is not too hard to trace the source of the errors back: Note that in image space the error term $(1-m)(n*AG)$ had a distinct peak at $\pm y_n$ (by about $10^2$), which is where the Fourier transform of $b$ also reaches a value about $10^5$ higher than at the lowest point (the high frequencies we noted before!).
Therefore it should be no surprise that we end up with an error that is a whopping $10^7$ higher at the edges of the facet. This error "spreads out" a bit more in grid space just because of the randomness of $G$, but we clearly have quite nasty worst case errors now.
Damage control
This is a problem inherent to prolate spheroidal wave functions - and likely all functions of the type we might select for $n$ and therefore $b$. However, note that we are constructing
$$b_j \ast m_i (n_j \ast  a_i G) \approx B_j \ast A_i G$$
and simply chose $y_n = y_B$. Clearly the support of $b$ must be the same as $B$, but we can easily make the support of $n$ larger as long as we maintain $b_j\ast n_j = B_j$. As we have just seen, the inverse wave function falls off very quickly away from the image edges, so let us simply truncate those parts away. That way, we will disconnect from the error peak of $(1-m_i) (n_j \ast G)$ (which will remain at $y_n$ and subsequentially get filtered out due to $y_b<y_n$) and substantially reduce the magnification of error due to $b$:
End of explanation
selM = (numpy.abs(x-x0) <= xM + 1e-13)
def pad_by_sel(sel, x): xp = numpy.zeros(len(sel), dtype=x.dtype); xp[sel] = x; return xp
ideal = conv(n, A*G)
approx = m * conv(n, A*G)
approx_r = pad_by_sel(selM, conv(n[numpy.abs(x) <= xM], (A*G)[selM]))
error1 = approx - ideal
error1_r = approx_r - ideal
pylab.plot(x[selM], ideal[selM].real, x[selM], approx_r[selM].real); pylab.legend(["n*AG", "[r] m(n*AG)"]);
mark_range("$x_A$", x0-xA,x0+xA); mark_range("$x_m$", x0-xM,x0+xM); pylab.show();
pylab.semilogy(fx[selM], numpy.abs(error1_r[selM])); pylab.legend(["[r] (1-m)(n*AG)"]); pylab.show();
pylab.semilogy(fx, numpy.abs(fft(error1)), fx, numpy.abs(fft(error1_r))); pylab.legend(["F[(1-m)(n*AG)]", "[r] F[(1-m)(n*AG)]"]);
mark_range("$y_n$", -yN,yN); mark_range("$y_B$", -yB,yB); pylab.show();
print("RMSE (w/o b):", numpy.sqrt(numpy.mean(numpy.abs(error1)**2)), "(image:",numpy.sqrt(numpy.mean(numpy.abs(fft(error1))**2)),")")
print("RMSE (w/o b, reduced):", numpy.sqrt(numpy.mean(numpy.abs(error1_r)**2)), "(image:",numpy.sqrt(numpy.mean(numpy.abs(fft(error1_r))**2)),")")
Explanation: Note that the base accuracy performance of the wave function has improved a bit simply because we have increased $y_n$ slightly. However, the important part is that we can make much better use of the results: Now we only lose an order of magnitude of error due to the convolution with $b$!
Downsampling
The point of this effort was to reduce the storage size of intermediate data involved in reconstructing a facet from all subgrids or a subgrid from all facets respectively. We have found a number of ways we can "compress" the signal, but have we achieved our purpose?
Let us again have a look at the approximation formula:
$$b_j \ast m_i (n_j \ast a_i G) \approx B_j \ast A_i G$$
We showed that $n_j \ast a_i' G \approx m_i (n_j \ast a_i G)$, therefore this term is approximately bound in grid space ($x_m$) as well as image space ($y_n$).
End of explanation
fx_r = N * coordinates(len(approx_r[selM]))
selN = numpy.abs(fx_r) <= yN
approx_core = ifft(fft(approx_r[selM])[selN])
approx_r2 = pad_by_sel(selM, ifft(pad_by_sel(selN, fft(approx_core))))
error1_r2 = ideal - approx_r2
pylab.semilogy(fx, numpy.abs(fft(ideal)),fx_r[selN], numpy.abs(fft(approx_core)));
mark_range("$y_n$", -yN,yN); mark_range("$y_B$", -yB,yB);
pylab.legend(["F[n*AG]", "[r²] F[m(n*AG)]"]); pylab.show();
pylab.plot(x[selM], ideal[selM].real, x[selM], approx_r2[selM].real); pylab.legend(["n*AG", "[r²] m(n*AG)"]);
mark_range("$x_A$", x0-xA,x0+xA); mark_range("$x_m$", x0-xM,x0+xM); pylab.show();
Explanation: Reducing the grid-size / image sampling rate has clearly not impacted the precision too much. Note that while the $(1-m)(n*AG)$ error term should be zero in the $x_m$ grid region, with the reduced size the errors get aliased back onto the subgrid. For our parameters this happens to cause an error "well" in the centre of the sub-grid.
However, limiting in grid space is only half the story. As established at the start we not only want to reduce our grid size down to $x_m$, but want to also get rid of image data beyond $y_n$:
End of explanation
pylab.semilogy(x[selM], numpy.abs(error1_r[selM]), x[selM], numpy.abs(error1_r2[selM]));
pylab.legend(["[r] (1-m)(n*AG)", "[r²] (1-m)(n*AG)"]); pylab.show();
pylab.semilogy(fx, numpy.abs(fft(error1)), fx, numpy.abs(fft(error1_r)), fx, numpy.abs(fft(error1_r2)));
pylab.legend(["F[(1-m)(n*AG)]", "[r] F[(1-m)(n*AG)]", "[r²] F[(1-m)(n*AG)]"]);
mark_range("$y_n$", -yN,yN); mark_range("$y_B$", -yB,yB); pylab.show();
print("RMSE (w/o b):", numpy.sqrt(numpy.mean(numpy.abs(error1)**2)), "(image:",numpy.sqrt(numpy.mean(numpy.abs(fft(error1))**2)),")")
print("RMSE (w/o b, reduced):", numpy.sqrt(numpy.mean(numpy.abs(error1_r)**2)), "(image:",numpy.sqrt(numpy.mean(numpy.abs(fft(error1_r))**2)),")")
print("RMSE (w/o b, twice reduced):", numpy.sqrt(numpy.mean(numpy.abs(error1_r2)**2)), "(image:",numpy.sqrt(numpy.mean(numpy.abs(fft(error1_r2))**2)),")")
Explanation: As intended, $n_j\ast a_iG$ drops to zero past $y_N$. This is clearly true even after downsampling, so we can truncate the signal in image space as well. This does not impact signal reconstruction in any visible way.
End of explanation
pylab.semilogy(x, numpy.abs(conv(b,error1_r)), x, numpy.abs(conv(b,error1_r2)));
mark_range("$x_A$", x0-xA,x0+xA); mark_range("$x_m$", x0-xM,x0+xM)
pylab.legend(["[r] b*(1-m)(n*AG)", "[r²] b*(1-m)(n*AG)"]); pylab.show();
pylab.semilogy(fx, numpy.abs(fft(conv(b,error1))), fx, numpy.abs(fft(conv(b,error1_r))), fx, numpy.abs(fft(conv(b,error1_r2))));
pylab.legend(["F[b*(1-m)(n*AG)]", "[r] F[b*(1-m)(n*AG)]", "[r²] F[b*(1-m)(n*AG)]"]);
mark_range("$y_n$", -yN,yN); mark_range("$y_B$", -yB,yB); pylab.show();
print("RMSE:", numpy.sqrt(numpy.mean(numpy.abs(conv(b,error1))**2)), "(image:",numpy.sqrt(numpy.mean(numpy.abs(fft(conv(b,error1)))**2)),")")
print("RMSE (reduced):", numpy.sqrt(numpy.mean(numpy.abs(conv(b,error1_r))**2)), "(image:",numpy.sqrt(numpy.mean(numpy.abs(fft(conv(b,error1_r)))**2)),")")
print("RMSE (twice reduced):", numpy.sqrt(numpy.mean(numpy.abs(conv(b,error1_r2))**2)), "(image:",numpy.sqrt(numpy.mean(numpy.abs(fft(conv(b,error1_r2)))**2)),")")
Explanation: The errors have not changed too much even though we have now significantly reduced the number of samples in both grid and image space. The downsampling operation eliminated the incidental "well" in grid space, but $n_j\ast a_iG \approx m_i (n_j \ast a_i G)$ holds just as it did before. This is what we were after!
For the final step, let's convolve with $b_j$. Notice that:
$$F_j = \sum_i \left( B_j \ast A_i G \right)
  \approx \sum_j b_j \ast n_i (m_j \ast a_i G)
   = b_j \ast \sum_j n_i (m_j \ast a_i G)$$
So we only need to do this once per facet, and can therefore easily work at full facet resolution. This doesn't matter for data transfers, but it is good to know that this can also be done efficiently from a computational point of view.
End of explanation
Gp = conv(b, G)
ideal = A * conv(n, Gp)
approx = A * conv(n, m * Gp)
error = A * conv(n, (1-m) * Gp)
pylab.plot(x, ideal.real, x, approx.real); pylab.legend(["A(n*G')", "A(n*mG')"]);
pylab.xticks(ticks); mark_range("$x_A$", x0-xA,x0+xA); mark_range("$x_m$", x0-xM,x0+xM); pylab.show();
pylab.semilogy(x, numpy.abs(error)); pylab.legend(["A(n*(1-m)G')"]); #pylab.ylim(1e-12,1e-7)
pylab.xticks(ticks); mark_range("$x_A$", x0-xA,x0+xA); mark_range("$x_m$", x0-xM,x0+xM); pylab.show();
pylab.semilogy(fx, numpy.abs(fft(error))); pylab.legend(["F[n*(1-m)G']"]);
pylab.xticks(fticks); mark_range("$y_n$", -yN,yN); mark_range("$y_B$", -yB,yB); pylab.show();
print("RMSE:", numpy.sqrt(numpy.mean(numpy.abs(error)**2)), "(image:",numpy.sqrt(numpy.mean(numpy.abs(fft(error))**2)),")")
Explanation: And we made it out the other end. Seems the approximation is indeed robust against truncation and down-sampling as intended.
Sub-grids recombination
So we can obtain facet data from sub-grids effectively. But does this also work in the "opposite" direction, if we want to construct a sub-grid from facets?
\begin{align}
S_i = A_iG &= A_i \left( \left( \sum_j B_j \right) \ast  G \right) = \sum_j A_i \left( B_j \ast  G \right) \
\end{align}
Let us split $A_i$ and $B_j$ as we did before. Now what we would like to establish is that:
$$A_i (B_j \ast G) = a_i m_i (n_j \ast b_j \ast G) \approx a_i \left(n_j \ast m_i (b_j \ast G)\right)$$
If we set $G' = b_j \ast G$ for the moment and use $a_im_i = a_i = A_i$, what we need to show is:
\begin{align}
a_i (n_j \ast G') &\approx a_i (n_j \ast m_i G') \
\Leftrightarrow a_i (n_j \ast (1 - m_i) G') &\approx 0 \
\Leftarrow \forall |x-x_0| \le x_a: n_j \ast (1 - m_i) G' &\approx 0 \
\Leftarrow \forall |x-x_0| \le x_a + x_n: (1 - m_i) G' &\approx 0
\end{align}
Which should be the case as long as $x_m \ge x_a + x_n$, like established before.
End of explanation
approx_r = conv(n[numpy.abs(x) <= xM], Gp[selM])
pylab.plot(x, conv(n, Gp).real, x, conv(n, m * Gp).real, x[selM], approx_r.real);
pylab.legend(["n*G'", "n*mG'", "[r]n*mG'"]); pylab.ylim(-1,1)
pylab.xticks(ticks); mark_range("$x_A$", x0-xA,x0+xA); mark_range("$x_m$", x0-xM,x0+xM); pylab.show();
Explanation: This is entirely dual to the original method -- we are doing everything backwards. Now $G$ needs to get convolved with $b_j$ right away, and we end with the $A_i$ mask cutting away a bit of grid space.
It is instructive to have closer look at the function of this mask operation, especially considering the reduced case:
End of explanation
print("Samples: %d (%d*%d=%d)" % (N, 1, N, N))
print("Samples (reduced): %d (%.2f*%d=%d, limit %.2f*%d=%d)" % (numpy.sum(selM), 2*xM,N,2*xM*N, 2*xA,N,2*xA*N, ))
print("Samples (twice reduced): %d (%.2f*%d=%d, limit %.2f*%d=%d)" %
      (numpy.sum(selN), 2*xM,2*yN,4*xM*yN, 2*xA,2*yB,4*xA*yB))
Explanation: Reproduction of the original signal quickly detoriates beyond the $x_A$ region. Especially note how $n\ast mG'$ tends to jump around violently beyond $x_A$. This is even worse for the reduced case, where the error aliases around the $x_m$ edges.
Data sizes and scaling
Let us have a closer look and check exactly how effective this has been. If $A (B \ast G) = B (A \ast G)$ was true, we would expect that we would be able to represent the dependencies between a facet and a subgrid with about $\lceil 4x_Ay_B \rceil$ samples.
Our method only allows truncating to $x_m$ in grid and $y_n$ in image space, therefore we need $\lceil 4x_my_n \rceil \approx \lceil 4(x_A+x_n)y_n \rceil$ samples:
End of explanation
def overhead_approx(yB, yN, xA, xN):
    return numpy.ceil(4*(xA+xN)*yN) / numpy.ceil(4*xA*yB)
def error_approx(yB, yN, xN, alpha=0, dim=1, hexagon=False):
    # gridding error
    assert yB < yN
    pswf = anti_aliasing_function(int(yN)*2, alpha, 2*numpy.pi*yN*xN)
    pswf /= numpy.prod(numpy.arange(2*alpha-1,0,-2, dtype=float)) # double factorial    
    grid_error = numpy.abs(numpy.sum(pswf[::2] - pswf[1::2]))
    # correction error
    b_error = numpy.abs(pswf[int(yN) + int(yB)])
    if dim >= 2 and hexagon:
        b_error *= numpy.abs(pswf[int(yN) + int(yB/2)])**(dim-1)
    else:
        b_error **= dim
    return numpy.abs(grid_error) / (2*xM) / b_error
print("Predicted error:", error_approx(yB, yN, xN))
print("Observed worst case:", numpy.max(numpy.abs(fft(conv(b,error1_r2)))))
Explanation: As expected, modulo some rounding of pixels at the edges. We are not actually that far away from the limits given by information theory. However, this is clearly for a very small data size, so what would happen if we started doing this for larger grids? Would we be able to retain the reached accuracy and efficiency?
To do this quickly, let us attempt to predict how much error we expect from our method. As we have seen before, the error of
$$b_j \ast m_i (n_j \ast a_i G)$$
consists of:
The error term $(1 - m_i) (n_j \ast a_i G)$, which for $x_m \ge x_A + x_n$ should be roughly equal to the PSWF grid values outside the $x_n$ region. This is mainly limited by the sampling of the wave function in image space, so we can approximate it by calculating the maximum-frequency value ($1/y_n$) in grid space.
We also have to note that this error gets "concentrated" more the smaller the sub-grid is. We should therefore divide by $2x_N$.
The $b_j$ convolution, which might magnify the error by the dynamic range of $\mathcal F b_j$. As $\mathcal F b_j(0) = 1$ for the PSWF, this is given by $\mathcal F n_j(y_B)$.
End of explanation
yNs = 2**numpy.arange(5,20); yBs = yNs*yB/yN; xNs = xN*yN/yNs
errs_approx = numpy.vectorize(error_approx)(yBs, yNs, xNs)
pylab.loglog(yBs, errs_approx); pylab.xlabel('$y_B$'); pylab.show()
Explanation: This is clearly quite pessimistic even when comparing against the worst case, but it gives us a decent starting point.
Let us assume that we increase the facet size ($y_B$). To make sure that we compare like with like, let us attempt to keep $\mathcal F n_j(y_B)$ constant. This means that we need:
$y_B / y_n$ constant so we sample the same point. This effectively means that the image margin stays constant relative to the size of the image.
$y_n x_n$ constant so the wave function parameter doesn't change. This effectively means that the grid margin stays constant in terms of the grid resolution.
This error approximation steadily improves from 100 to 100,000 (we increase size in factor-of-two steps to reduce jitter):
End of explanation
xAs = numpy.hstack([numpy.arange(0.01,0.05,0.01), numpy.arange(0.05,0.2,0.05)])
for xA_ in xAs:
    pylab.semilogx(yNs, numpy.vectorize(overhead_approx)(yNs*0.75, yNs, xA_, xN*yN/yNs));
pylab.legend(["xA=%.2f" % xA for xA in xAs]); pylab.show()
Explanation: At the same time, the "overhead" of our solution decreases slowly depending on the sub-grid size $x_A$, eventually approximating $y_N/y_B$. This is of course because the constant grid margin eventually stops mattering:
End of explanation
ov = 1.5
pylab.rcParams['figure.figsize'] = 16, 4
for yN_ in [ 256, 2048]:
    xNyNs = numpy.hstack([numpy.arange(0.125, 1, 0.125), numpy.arange(1, 5, 1), numpy.arange(5, 6, 0.125), numpy.arange(6,7,0.25)])
    xNs = xNyNs / yN_
    for xA_ in xAs:
        yBs = (yN_ + xNyNs / xA_) / ov
        sel = yBs < yN_
        if numpy.sum(sel) == 0: continue
        errs_approx = numpy.vectorize(error_approx)(yBs[sel], yN_, xNs[sel], dim=2)
        pylab.semilogy(xNyNs[sel], errs_approx)
    pylab.ylim(pylab.gca().get_ylim()[0], 1); pylab.xlabel("$x_ny_n$")
    pylab.legend(["$x_A=%.2f$" % xA for xA in xAs]); pylab.title("$o=%.1f, y_N=%.1f$" % (ov, yN_)); pylab.show()
Explanation: So how much accuracy can we get for a certain overhead, optimally? If we settle on a certain overhead value $o$, we can calculate $y_n$ from a value of $x_ny_n$:
$$o = \frac{4(x_A+x_n)y_n}{4x_Ay_B} \quad\Rightarrow\quad y_n = y_Bo - \frac{x_ny_n}{x_A}
\quad\text{, or: }
y_B = \frac{y_n + \frac{x_ny_n}{x_A}}{o}
$$
End of explanation
yP = yN + yB/2
IIIp_d = numpy.zeros(N)
for i in range(-N // int(yP) // 4+1, (N-1) // int(yP) // 4+1):
    IIIp_d[N//2+2*int(yP)*i] = 1
pylab.plot(fx, IIIp_d, fx, fft(m).real, fx, fft(m).real - conv(fft(m), IIIp_d).real)
pylab.legend(["III", "Fm", "Fm-Fm*III"], loc=2)
pylab.ylim(-10,10)
mark_range("yN", -yN, yN); mark_range("yB", -yB, yB); mark_range("yP", -yN-yB, yN+yB); mark_range("2yP", -2*yP, 2*yP)
Explanation: Facet Split
So far we have been rather cavalier with using arrays larger than their "true" size. This simplifies things mathematically, yet we must remember that we don't actually have the $b_j \ast F_j$ available. Instead, we will handle an FFT grid, which is actually sampled at a finite rate - which we can represent as multiplying by a comb function:
$$(b_j \ast G) \operatorname{III}_{\frac1{2y_P}}$$
The choice of $y_P$ corresponds to the actual memory size of the facet for the purpose of this algorithm -- the finer we need to sample, the more memory we will have to use to hold the facet in our calculations. Clearly we don't lose any information as long as $y_P \ge y_B$, so we might be tempted to choose $y_P = y_B$. Yet that might not be enough, after all if we try to find $S_i$ in the obvious way we would like:
$$n_j \ast m_i (b_j \ast G) = n_j \ast m_i (b_j \ast G) \operatorname{III}_{\frac1{2y_P}}$$
However, this is not the case. To understand why, note that if we Fourier-transform both sides and re-order slightly, we get:
\begin{align}
\mathcal F n_j (\mathcal F m_i \ast  \mathcal Fb_j \mathcal FG) &=
\mathcal F n_j (\operatorname{III}{2y_P} \ast \mathcal F m_i \ast  \mathcal Fb_j \mathcal FG ) \
\Leftrightarrow
\mathcal F n_j \left( ( \mathcal F m_i - \operatorname{III}{2y_P} \ast \mathcal F m_i) \ast  \mathcal Fb_j \mathcal FG \right) &= 0 \
\Leftarrow
\forall |y| < y_N: \left ( \mathcal F m_i - \operatorname{III}_{2y_P} \ast \mathcal F m_i) \ast  \mathcal Fb_j \mathcal FG\right &= 0 \
\Leftarrow
\forall |y| < y_N + y_B: \left  \mathcal Fm_i - \operatorname{III}_{2y_P} \ast \mathcal F m_i\right &= 0 \
\end{align}
Let us have a closer look:
End of explanation
xM = xA + 2 * xN
mp = conv(m, n).real
pylab.semilogy(x, numpy.maximum(1e-15, numpy.abs(m)));
pylab.semilogy(x, numpy.abs(mp));
pylab.ylim((1e-13, 1e1)); pylab.legend(["m", "m*n"])
mark_range("$x_0+x_A+2x_N$", x0-xM, x0+xM); mark_range("$x_0-x_A-x_N$", x0+xA+xN, x0-xA-xN);
Explanation: The problem here is that $\mathcal Fm$ never approaches zero - it is a sinc function that doesn't fall off very much at all. Therefore we have $\mathcal F m_i(y) \ne \left\operatorname{III}_{2y_P} \ast \mathcal F m_i \right$ no matter how much we restrict the $y$ region.
However, we can actually choose $m$: The two properties we need to conserve is that $m(x)=1$ for $|x| < x_A+x_N$, and we would like to fall it to zero by $x_M$ in order to limit the amount of data we need to exchange between facets and subgrids. This means that if we reserve some space between $x_A+x_N$ and $x_M$, we can fill it in whatever way we want in order to make $m$ more manageable.
In fact we can simply utilise the same prolate spheroidal wave function we used before:
$$m' = m \ast n$$
We re-use the same function just for convenience, this way we don't have to involve a second PSWF function. In this case, we know this increases the support of $m'$ in grid space to $x_{m'}=x_A+2x_N$:
End of explanation
pylab.semilogy(fx, IIIp_d, fx, numpy.abs(fft(mp)), fx, numpy.abs(fft(mp) - conv(fft(mp), IIIp_d)))
mark_range("$y_N$", -yN, yN); mark_range("$y_B$", -yB, yB); mark_range("$y_P$", -yP, yP)
mark_range("$2y_P$", -2*yP, 2*yP);  mark_range("$y_B+y_N$\n$=2y_P-y_N$", -yB-yN, yB+yN)
pylab.legend(["III", "Fm'", "Fm'-Fm'*III"], loc=2);
Explanation: However, now this also limits $m$ to $y_N$ in image space:
End of explanation
pylab.semilogy(fx, numpy.abs(conv(fft(b) * fft(G), fft(mp) - conv(fft(mp), IIIp_d))),
               fx, numpy.abs(fft(n) * conv(fft(b) * fft(G), fft(mp) - conv(fft(mp), IIIp_d))))
mark_range("$y_N$", -yN, yN); mark_range("$y_B$", -yB, yB); mark_range("yP", -yP, yP)
pylab.legend(["$FbFG*(Fm-Fm*III)$", "$Fn(FbFG*(Fm-Fm*III$))"]);
Explanation: Good starting point. We just need to figure out how far we can reduce the image size (here $y_P$) without impacting our ability to reconstruct the image. What we need is:
\begin{align}
\mathcal F n_j (\mathcal F m'i \ast  \mathcal Fb_j \mathcal FG) &=
\mathcal F n_j (\operatorname{III}{2y_P} \ast \mathcal F m'_i \ast  \mathcal Fb_j \mathcal FG) \
\Leftarrow \forall |y| < y_N + y_B:
\left \mathcal F m_i' - \operatorname{III}_{2y_P} \ast \mathcal F m_i' \right &= 0
\end{align}
As the image "loops" every $2y_P$ and the width of $\mathcal Fm'$ is $y_N$, this means:
$$2y_P - y_N \ge y_N + y_B \;\Rightarrow\; y_P \ge y_N + \frac12 y_B$$
With this sizing, the error relative to $\mathcal Fm_i\ast \mathcal Fb_j\mathcal FG$ is zero within $y_N$, and therefore close to zero entirely for $\mathcal Fn_j (\mathcal Fm_i\ast \mathcal Fb_j\mathcal FG)$:
End of explanation
def red_2yP(xs): return extract_mid(xs, int(2*yP))
ref = fft(n) * conv(fft(b) * fft(G), fft(mp))
reduced = red_2yP(fft(n)) * conv(red_2yP(fft(b) * fft(G)), red_2yP(fft(mp)))
pylab.semilogy(red_2yP(fx),  numpy.abs(reduced - red_2yP(ref)));
Explanation: This shows that by using the low-pass-filtered $m'$ we can indeed perform the entire computation using just $2y_P$ sample points:
End of explanation
nsubgrid = int(1 / (2*xA)); subgrid_size = int(N * 2*xA)
assert nsubgrid * subgrid_size == N
nfacet = int(N / (2*yB)); facet_size = int(2*yB)
assert nfacet * facet_size == N
print(nsubgrid,"subgrids,",nfacet,"facets")
subgrid = numpy.empty((nsubgrid, subgrid_size), dtype=complex)
facet = numpy.empty((nfacet, facet_size), dtype=complex)
for i in range(nsubgrid):
    subgrid[i] = extract_mid(numpy.roll(G, -i * subgrid_size), subgrid_size)
FG = fft(G)
for j in range(nfacet):
    facet[j] = extract_mid(numpy.roll(FG, -j * facet_size), facet_size)
Explanation: Note that we are only truncating terms here that we already know to fall to zero. However, the choice of $y_P$ is still quite subtle here, as we need padding to make sure that the convolutions work out in the correct way.
Actual implementation
So to review, what we have identified are two approximations:
\begin{align}
S_i &\approx
a_i \sum_j \left( n_j \ast m_i (b_j \ast F_j)\right) \
F_j &\approx b_j \ast \sum_i m_i (n_j \ast S_i)
\end{align}
Some quick observations: The slight asymmetry comes from the fact that $a_i = A_i$ and assumed to be a mask (and therefore $a_iS_i = S_i$), but $b_j \ne B_j$. Also note that $b_j$ always gets applied on the image side, whereas $n_j$ gets applied on the grid plane, therefore appearing in pairs if we consider the round-trip. This is correct, and has the same underlying reason why we can use the same type of kernels for gridding and degridding.
End of explanation
yB_size = int(yB*2)
yN_size = int(yN*2)
# Find yP that allows us to align subgrid masks easily (brute force!)
yP_size = int(yP*2)+1
while numpy.abs(yP_size * 2*xA - int(yP_size * 2*xA)) >= 1e-13 or \
      numpy.abs(yP_size * 2*xM - int(yP_size * 2*xM)) >= 1e-13:
    yP_size+=1
xM_yP_size = int(2*xM*yP_size)
xMxN_yP_size = int(2*(xM+xN)*yP_size)
xMxN_yN_size = int(2*(xM+xN)*yN_size)
pswf = anti_aliasing_function(yN_size, alpha, 2*numpy.pi*yN*xN).real
Fb = 1/extract_mid(pswf, yB_size)
b = ifft(pad_mid(Fb, N))
n = ifft(pad_mid(pswf, N))
Fn = extract_mid(fft(n), yN_size)
FBjFj = numpy.empty((nfacet, yP_size), dtype=complex)
for j in range(nfacet):
    FBjFj[j] = pad_mid(facet[j] * Fb, yP_size)
Explanation: So let us start by reconstructing sub-grids from facets. First step is to convolve in $b_j$, derived from the PSWF. This is a cheap multiplication in image space. We then pad to $2y_P$ which yields us
$$\operatorname{III}_{2y_P} \ast  \mathcal Fb_j \mathcal FG$$
in FFT convention.
End of explanation
facet_m0 = conv(n, pad_mid(numpy.ones(int(N*xM*2)), N))
facet_m0_trunc = ifft(extract_mid(fft(facet_m0), yP_size))
pylab.semilogy(coordinates(yP_size), numpy.abs(facet_m0_trunc));
mark_range("$x_M$", xM, -xM); mark_range("$x_M+x_N$", -xM-xN, xM+xN);
Explanation: Next we need to cut out the appropriate sub-grid for $m_i$ for the subgrid we want to construct (here: the centre subgrid). As explained above, to do this with a fact that has only been padded to $y_P$ we need to consider the truncated $\mathcal F m_i$:
$$\mathcal F^{-1}[\Pi_{2y_P} \mathcal F m_i \ast  \operatorname{III}{2y_P} \ast \mathcal Fb_j \mathcal FG]
= \mathcal F^{-1}[\Pi{2y_P} \mathcal F m_i] \mathcal F^{-1}[ \operatorname{III}{2y_P} \ast \mathcal Fb_j \mathcal FG]
$$
First let us construct the $m_i' = \mathcal F^{-1}[\Pi{2y_P} \mathcal F m_i]$ terms:
End of explanation
MiBjFj = numpy.empty((nsubgrid, nfacet, yP_size), dtype=complex)
assert numpy.abs(yP_size * 2*xA - int(yP_size * 2*xA)) < 1e-13
for i in range(nsubgrid):
    for j in range(nfacet):
        MiBjFj[i,j] = facet_m0_trunc * numpy.roll(ifft(FBjFj[j]), -i*int(yP_size * 2*xA))
Explanation: Note that we would clearly not construct them separately for a real pipeline, as they are simply shifted. Due to the truncation in frequency space these are not quite top-hat functions any more. These terms now get used to extract the sub-grid data from each facet:
End of explanation
Fn = extract_mid(fft(n), yN_size)
NjMiBjFj = numpy.empty((nsubgrid, nfacet, yN_size), dtype=complex)
for i in range(nsubgrid):
    for j in range(nfacet):
        NjMiBjFj[i,j] = Fn * extract_mid(fft(MiBjFj[i,j]), yN_size) * yP_size / N
Explanation: Next step is to multiply in $n_j$ in order to un-do the effects of $b$ and cut out the garbage between $y_N$ and $y_P$. This means we arrive at:
$$n_j \ast m'i (b_j \ast G) \operatorname{III}{\frac1{2y_P}} =
  \mathcal F n_j \mathcal F\left[ m'i (b_j \ast G) \operatorname{III}{\frac1{2y_P}} \right]$$
Which we truncate further to $y_N$, as we are now done with convolutions in image space.
End of explanation
fig = pylab.figure(figsize=(16, 8))
ax1, ax2 = fig.add_subplot(211), fig.add_subplot(212)
for i in range(nsubgrid):
    Gr = numpy.roll(G, -i * subgrid_size)
    for j in range(nfacet):
        Grr = ifft(numpy.roll(fft(Gr), -j * facet_size))
        ax1.semilogy(yN_size*coordinates(yN_size),
                       numpy.abs(extract_mid(fft(conv(n, facet_m0 * conv(b, Grr))), yN_size)
                                 - NjMiBjFj[i,j]))
        ax2.semilogy(x, numpy.abs(ifft(pad_mid(NjMiBjFj[i,j], N))));
ax1.set_title("Error compared with $n_j * m_i(B_j * G)$")
mark_range("yN", -yN, yN,ax=ax1); mark_range("yB", -yB, yB,ax=ax1);
ax2.set_title("Signal left in grid space")
mark_range("xA", -xA, xA, ax=ax2); mark_range("xM", -xM, xM, ax=ax2); mark_range("xN+xN", -xM-xN, xM+xN, ax=ax2)
Explanation: Quick mid-point accuracy check against the approximation formula using full resultion. We should be looking at only rounding errors.
Furthermore, when padded to the full density the Fourier transform should fall to the error level provided by the PSWF past $x_N+x_M$ - indicating that we can reduce the sampling rate / grid size without losing information:
End of explanation
assert(numpy.abs(int(yN * xM) - yN * xM) < 1e-13)
assert(numpy.abs(int(1/2/xM) - 1/2/xM) < 1e-13)
xM_yN_size = int(xM*2*yN*2)
RNjMiBjFj = numpy.empty((nsubgrid, nfacet, xM_yN_size), dtype=complex)
RNjMiBjFj[:,:] = NjMiBjFj[:,:,::int(1/2/xM)]
print("Split", nfacet,"*",facet_size, "=", N, "data points into",
      nfacet,'*',nsubgrid,'*',xM_yN_size,'=',nsubgrid*nfacet*xM_yN_size, ", overhead", nsubgrid*xM_yN_size/facet_size-1)
for i,j in itertools.product(range(nsubgrid), range(nfacet)):
    Grr = ifft(numpy.roll(fft(numpy.roll(G, -i * subgrid_size)), -j * facet_size))
    pylab.semilogy(xM*2*coordinates(int(xM*2*N)),
                   numpy.abs( extract_mid(facet_m0 * conv(conv(n,b), Grr), int(xM*2*N))  -
                              ifft(pad_mid(RNjMiBjFj[i,j], int(xM*2*N)))))
pylab.title("Error compared with $m_i (B_j * G)$")
mark_range("xA", -xA, xA); mark_range("xM", -xM, xM)
Explanation: So our final step is to reduce the sampling rate (sub-grid size) to $x_M$. This is the step where we actually introduce the bulk of our error, as the "tail" regions outside $x_M$ get aliased in. As established before, this especially copies the $x_M+x_N$ region inside, which doesn't hurt because we are only interested in the centre $x_A$ part.
As long as $x_M$ divides the grid well, this can simply be achieved by selecting every $\frac1{x_M}$th sample:
End of explanation
fig = pylab.figure(figsize=(16, 8))
ax1, ax2 = fig.add_subplot(211), fig.add_subplot(212)
err_sum = 0
for i in range(nsubgrid):
    approx = numpy.zeros(int(xM*2*N), dtype=complex)
    for j in range(nfacet):
        approx += numpy.roll(pad_mid(RNjMiBjFj[i,j], int(xM*2*N)), j * int(xM*2*yB*2))
    approx = extract_mid(ifft(approx), int(xA*2*N))
    ax1.semilogy(xA*2*coordinates(subgrid_size), numpy.abs( approx - subgrid[i] ))
    ax2.semilogy(N*coordinates(subgrid_size), numpy.abs( fft(approx - subgrid[i]) ))
    err_sum += numpy.abs(approx - subgrid[i])**2
mark_range("xA", -xA, xA, ax=ax1); ax1.set_title("Error compared with $S_i = A_i G$")
mark_range("yB", -yB, yB, ax=ax2); mark_range("yN", -yN, yN, ax=ax2); 
pylab.show()
print("MRSE:", numpy.sqrt(numpy.mean(err_sum)))
Explanation: At this point, all that is left is to put together the sum
$$S_i = a_i \sum_j \left( n_j \ast m_i (b_j \ast F_j)\right)$$
which eliminates $b_j$. Note that we know $A$ to be a pure mask, therefore this is simply truncation in grid space.
End of explanation
seed = numpy.random.randint(2**31)
@interact(N=(0,8192),
          x0=(-0.5,0.5,0.1), xA=(0,0.5,0.01), xN=(0,0.5,0.01), xM=(0,0.5,0.01),
          yN=(0,1024,25), yB=(0,1024,25), alpha=(0,20))
def test_it(N=N, x0=x0, xA=xA, xN=xN, xM=xM, yN=yN, yB=yB, alpha=alpha):
    x = coordinates(N); fx = N * x
    G = numpy.random.RandomState(seed).rand(N) - .5
    A = numpy.roll(pad_mid(numpy.ones(int(2*xA*N)), N), int(x0*N))
    m = numpy.roll(pad_mid(numpy.ones(int(2*xM*N)), N), int(x0*N))
    selM = (m == 1)
    selM0 = pad_mid(numpy.ones(int(2*xM*N),dtype=bool), N)
    selN = pad_mid(numpy.ones(int(2*xM*2*yN),dtype=bool), numpy.sum(selM0))
    pswf = anti_aliasing_function(int(yN*2), alpha, 2*numpy.pi*yN*xN)
    pswf /= numpy.prod(numpy.arange(2*alpha-1,0,-2, dtype=float)) # double factorial
    n = ifft(pad_mid(pswf, N)).real
    m = conv(n,m).real
    b = ifft(pad_mid(1/extract_mid(pswf, int(yB*2)), N)).real
    Gp = conv(b, G)
    ideal = A * conv(n, Gp)
    approx = A * conv(n, m * Gp)
    error = A * conv(n, (1-m) * Gp)
    approx_r = A[selM] * conv(n[selM0], Gp[selM])
    error_r = approx_r - ideal[selM]
    approx_core = ifft(fft(n[selM0])[selN] * fft(Gp[selM])[selN])
    approx_r2 = A[selM] * ifft(pad_by_sel(selN, fft(approx_core)))
    error_r2 = approx_r2 - ideal[selM]
    print("PSWF parameter:", 2*numpy.pi*xN*yN)
    print("Worst error:", numpy.max(numpy.abs(error_r2)), "(image:", numpy.max(numpy.abs(fft(error))),
          ", predicted:", error_approx(yB, yN, xN, alpha=alpha),")")
    print("RMSE:", numpy.sqrt(numpy.mean(numpy.abs(error)**2)), "(reduced:", numpy.sqrt(numpy.mean(numpy.abs(error_r)**2)),
              ", +downsample:", numpy.sqrt(numpy.mean(numpy.abs(error_r2)**2)), ")")
    print("RMSE image:", numpy.sqrt(numpy.mean(numpy.abs(fft(error))**2)), "(reduced:", numpy.sqrt(numpy.mean(numpy.abs(fft(error_r))**2)),
              ", +downsample:", numpy.sqrt(numpy.mean(numpy.abs(fft(error_r2))**2)), ")")
    print("Samples:", len(approx_core), "(",2*xM, "*", 2*yN,"=",4*xM*yN, ", overhead: %.2f)" % (xM*yN/xA/yB))
    ticks = coordinates(10)
    pylab.figure(figsize=(16, 18))
    pylab.subplot(4,1,1); pylab.title("Input");
    pylab.plot(x, A, x, conv(b,n).real, x, n, x, G, x, m); pylab.legend(["A", "B", "n", "G", "m"]); pylab.xticks(ticks);
    mark_range("$x_A$", x0-xA,x0+xA); mark_range("$x_m$", x0-xM,x0+xM)
    pylab.subplot(4,1,2); pylab.title("Output")
    pylab.plot(x, conv(n, Gp).real, x, conv(n, m * Gp).real, x[selM], conv(n[selM0], Gp[selM]).real);
    pylab.legend(["n*(b*G)", "n*m(b*G)", "[r]n*m(b*G)"]); pylab.ylim((-0.6,0.6)); pylab.xticks(ticks)
    mark_range("$x_A$", x0-xA,x0+xA); mark_range("$x_m$", x0-xM,x0+xM);
    pylab.subplot(4,1,3); pylab.title("Errors (Grid space)");
    pylab.semilogy(x, numpy.abs(n), x, numpy.abs(conv(n, (1-m) * Gp)), x[selM], numpy.abs(error_r), x[selM], numpy.abs(error_r2));
    pylab.legend(["n", "n*(1-m)(b*G)", "[r] A(n*(1-m)(b*G))", "[r+d] A(n*(1-m)(b*G))"]); pylab.xticks(ticks); 
    mark_range("$x_A$", x0-xA,x0+xA); mark_range("$x_n$", -xN,xN)
    pylab.subplot(4,1,4); pylab.title("Errors (Image space)")
    pylab.semilogy(fx, numpy.abs(fft(n)), fx, numpy.abs(fft(b)), fx, numpy.abs(fft(conv(n, m * Gp))), fx, numpy.abs(fft(error)),
                   N*coordinates(len(error_r)), numpy.abs(fft(conv(n[selM0], Gp[selM]))),
                   N*coordinates(len(error_r)), numpy.abs(fft(error_r)));
    pylab.legend(["F[n]", "F[b]", "F[n*m(b*G)]", "F[A(n*(1-m)(b*G))]", "[r] F[n*m(b*G)]", "[r] F[A(n*(1-m)(b*G))]"]);
    mark_range("$y_n$", -yN,yN); mark_range("$y_B$", -yB,yB);
    pylab.xticks(N*ticks); pylab.show()
seed = numpy.random.randint(2**31)
@interact(N=(0,8192),
          x0=(-0.5,0.5,0.1), xA=(0,0.5,0.01), xN=(0,0.5,0.01), xM=(0,0.5,0.01),
          yN=(0,1024,25), yB=(0,1024,25), alpha=(0,20))
def test_it(N=N, x0=x0, xA=xA, xN=xN, xM=xM, yN=yN, yB=yB, alpha=alpha):
    x = coordinates(N); fx = N * x
    G = numpy.random.RandomState(seed).rand(N) - .5
    A = numpy.roll(pad_mid(numpy.ones(int(2*xA*N)), N), int(x0*N))
    m = numpy.roll(pad_mid(numpy.ones(int(2*xM*N)), N), int(x0*N))
    selM = (m == 1)
    selM0 = pad_mid(numpy.ones(int(2*xM*N),dtype=bool), N)
    selN = pad_mid(numpy.ones(int(2*xM*2*yN),dtype=bool), numpy.sum(selM0))
    pswf = anti_aliasing_function(int(yN*2), alpha, 2*numpy.pi*yN*xN)
    pswf /= numpy.prod(numpy.arange(2*alpha-1,0,-2, dtype=float)) # double factorial
    n = ifft(pad_mid(pswf, N)).real
    b = ifft(pad_mid(1/extract_mid(pswf, int(yB*2)), N)).real
    m = conv(n,m).real
    ideal = conv(n, A*G)
    approx = m * conv(n, A*G)
    error = (1-m) * conv(n, A*G)
    error_b = conv(b, error)
    approx_r = conv(n[selM0], (A*G)[selM])
    error = ideal - approx
    error_r = ideal[selM] - approx_r
    approx_core = ifft(fft(approx_r)[selN])
    approx_r2 = ifft(pad_by_sel(selN, fft(approx_core)))
    error_r2 = ideal[selM] - approx_r2
    error_rb = conv(b, pad_by_sel(selM, error_r))
    error_r2b = conv(b, pad_by_sel(selM, error_r2))
    print("PSWF parameter:", 2*numpy.pi*xN*yN)
    print("Worst error:", numpy.max(numpy.abs(error_r2b)), "(image:", numpy.max(numpy.abs(fft(error_r2b))),
          ", predicted:", error_approx(yB, yN, xN, alpha=alpha),")")
    print("RMSE:", numpy.sqrt(numpy.mean(numpy.abs(error_b)**2)), "(reduced:", numpy.sqrt(numpy.mean(numpy.abs(error_rb)**2)),
              ", +downsample:", numpy.sqrt(numpy.mean(numpy.abs(error_r2b)**2)), ")")
    print("RMSE image:", numpy.sqrt(numpy.mean(numpy.abs(fft(error_b))**2)), "(reduced:", numpy.sqrt(numpy.mean(numpy.abs(fft(error_rb))**2)),
              ", +downsample:", numpy.sqrt(numpy.mean(numpy.abs(fft(error_r2b))**2)), ")")
    print("Samples:", len(approx_core), "(",2*xM, "*", 2*yN,"=",4*xM*yN, ", overhead: %.2f)" % (xM*yN/xA/yB))
    # Input graph
    ticks = coordinates(10)
    pylab.figure(figsize=(16, 18))
    pylab.subplot(4,1,1); pylab.title("Input");
    pylab.plot(x, A, x, conv(b,n), x, n, x, G, x, m); pylab.legend(["A", "B", "n", "G", "m"]); pylab.xticks(ticks);
    mark_range("$x_A$", x0-xA,x0+xA); mark_range("$x_m$", x0-xM,x0+xM)
    # Output graph
    pylab.subplot(4,1,2); pylab.title("Output");
    pylab.plot(x, conv(b, ideal), x, conv(b,approx), x, conv(b,pad_by_sel(selM, approx_r2)))
    pylab.ylim((-0.5,0.5)); pylab.legend(["B*AG", "b*m(n*aG)", "[r+d] b*m(n*aG)"]);
    mark_range("$x_A$", x0-xA,x0+xA); mark_range("$x_m$", x0-xM,x0+xM)
    # Error graph (image space)
    pylab.subplot(4,1,3); pylab.title("Errors (Grid space)");
    pylab.semilogy(x, numpy.abs(n), x, numpy.abs(error), x[selM], numpy.abs(error_r), x[selM], numpy.abs(error_r2), x, numpy.abs(error_r2b));
    mark_range("$x_n$", -xN,xN); mark_range("$x_m$", x0-xM,x0+xM)
    pylab.legend(["n","(1-m)(n*AG)","[r] (1-m)(n*AG)","[r+d] (1-m)(n*AG)","[r+d] b*(1-m)(n*AG)"]);
    pylab.xticks(ticks);
    # Error graph (frequency space)
    pylab.subplot(4,1,4); pylab.title("Errors (Image space)")
    pylab.semilogy(fx, numpy.abs(fft(n)), fx, numpy.abs(fft(b)), fx, numpy.abs(fft(ideal)), fx, numpy.abs(fft(error)),
                   N*coordinates(len(error_r)), numpy.abs(fft(error_r2)), fx, numpy.abs(fft(error_r2b)));
    mark_range("$y_n$", -yN,yN); mark_range("$y_B$", -yB,yB);
    pylab.legend(["n", "b", "n*AG", "(1-m)(n*AG)", "[r+d] (1-m)(n*AG)", "[r+d] b*(1-m)(n*AG)"]);
    pylab.xticks(N * ticks);
    pylab.show()
Explanation: Note the pattern of the errors in image space - this is the position dependence of the accuracy pattern from the $b_j$ multiplication. As we stitch together 5 facets, this pattern repeats 5 times.
However, note how the "stitching" works out really well - accuracy is especially good around the $y_B$ areas where the data from two facets "overlaps". This is likely because are effectively averaging over two samples here, which boost our accuracy.
Interactive playground
Test different parameters:
End of explanation
seed = numpy.random.randint(2**31)
# Subgrid region marker, assuming size to be even
def grid_patch(size):
    return patches.Rectangle((-(size+1)/theta/2, -(size+1)/theta/2), size/theta, size/theta, fill=False, linestyle="dashed")
def pad_by_sel_(shape, sel, x):
    xp = numpy.zeros(shape, dtype=x.dtype); xp[sel] = x; return xp    
@interact(N=(0,1024,128),
          x0=(-0.5,0.5,0.1), y0=(-0.5,0.5,0.1), xA=(0,0.5,0.01), xN=(0,0.5,0.01), xM=(0,0.5,0.01),
          yN=(0,1024,25), yB=(0,1024,25), alpha=(0,20))
def test_it(N=N, x0=x0, y0=-0.1, xA=xA, xN=xN, xM=xM, yN=yN, yB=yB, alpha=alpha, hexagon=False):
    x,y = coordinates2(N)
    G = numpy.random.RandomState(seed).rand(N, N) - .5
    A = numpy.roll(pad_mid(numpy.ones((int(2*xA*N), int(2*xA*N))), N), (int(y0*N), int(x0*N)), (0,1))
    m = numpy.roll(pad_mid(numpy.ones((int(2*xM*N), int(2*xM*N))), N), (int(y0*N), int(x0*N)), (0,1))
    A = numpy.where((numpy.abs(x-x0) <= xA) & (numpy.abs(y-y0) <= xA), 1, 0)
    m = numpy.where((numpy.abs(x-x0) <= xM) & (numpy.abs(y-y0) <= xM), 1, 0)
    selM_1 = pad_mid(numpy.ones(int(2*xM*N), dtype=bool), N)
    selM = numpy.ix_(numpy.roll(selM_1, int(y0*N)), numpy.roll(selM_1, int(x0*N)))
    selM0 = numpy.ix_(selM_1, selM_1)
    selN_1 = pad_mid(numpy.ones(int(2*xM*2*yN),dtype=bool), numpy.sum(selM_1))
    selN = numpy.ix_(selN_1, selN_1)
    pswf = anti_aliasing_function((int(yN*2),int(yN*2)), alpha, 2*numpy.pi*yN*xN)
    pswf /= numpy.prod(numpy.arange(2*alpha-1,0,-2, dtype=float)) # double factorial
    n = ifft(pad_mid(pswf, N)).real
    b = pad_mid(1/extract_mid(pswf, int(yB*2)), N)
    if hexagon: b[numpy.where(numpy.abs(N*x) > yB - numpy.abs(N*y) / 2)] = 0
    b = ifft(b).real
    m = conv(n,m).real
    ideal = conv(n, A*G)
    approx = m * conv(n, A*G)
    error = (1-m) * conv(n, A*G)
    error_b = conv(b, error)
    approx_r = conv(n[selM0], (A*G)[selM])
    error = ideal - approx
    error_r = ideal[selM] - approx_r
    approx_core = ifft(fft(approx_r)[selN])
    approx_r2 = ifft(pad_by_sel_(approx_r.shape, selN, fft(approx_core)))
    error_r2 = ideal[selM] - approx_r2
    error_rb = conv(b, pad_by_sel_(ideal.shape, selM, error_r))
    error_r2b = conv(b, pad_by_sel_(ideal.shape, selM, error_r2))
    print("PSWF parameter:", 2*numpy.pi*xN*yN)
    print("Worst error:", numpy.max(numpy.abs(error_r2b)), "(image:", numpy.max(numpy.abs(fft(error_r2b))),
          ", predicted:", error_approx(yB, yN, xN, alpha=alpha, dim=2, hexagon=hexagon),")")
    print("RMSE:", numpy.sqrt(numpy.mean(numpy.abs(error_b)**2)), "(reduced:", numpy.sqrt(numpy.mean(numpy.abs(error_rb)**2)),
              ", +downsample:", numpy.sqrt(numpy.mean(numpy.abs(error_r2b)**2)), ")")
    print("RMSE image:", numpy.sqrt(numpy.mean(numpy.abs(fft(error_b))**2)), "(reduced:", numpy.sqrt(numpy.mean(numpy.abs(fft(error_rb))**2)),
              ", +downsample:", numpy.sqrt(numpy.mean(numpy.abs(fft(error_r2b))**2)), ")")
    print("Samples:", approx_core.shape[0] * approx_core.shape[1], "(",2*xM, "² *", 2*yN,"² =",(4*xM*yN)**2, ", overhead: %.2f)" % (xM*yN/xA/yB)**2)
    fig = pylab.figure(figsize=(16,4))
    show_grid(numpy.log(numpy.maximum(1e-20,numpy.abs(n))) / numpy.log(10), "n", N, axes=fig.add_subplot(121))
    show_grid(numpy.log(numpy.maximum(1e-20,numpy.abs(b))) / numpy.log(10), "b", N, axes=fig.add_subplot(122))
    fig = pylab.figure(figsize=(16,4))
    show_image(numpy.log(numpy.maximum(1e-20,numpy.abs(fft(n)))) / numpy.log(10), "F[n]", N, axes=fig.add_subplot(121))
    show_image(numpy.log(numpy.maximum(1e-20,numpy.abs(fft(b)))) / numpy.log(10), "F[b]", N, axes=fig.add_subplot(122))
    fig = pylab.figure(figsize=(16,4))
    show_grid(numpy.abs(conv(b,ideal)), "b*n*AG", N, axes=fig.add_subplot(121))
    show_grid(numpy.abs(conv(b,pad_by_sel_(ideal.shape, selM, approx_r2))), "[r+d] b*m(n*aG)", N, axes=fig.add_subplot(122))
    fig = pylab.figure(figsize=(16,4))
    show_grid(numpy.log(numpy.maximum(1e-20,numpy.abs(error))) / numpy.log(10), "log_10 (1-m)(n*AG)", N, axes=fig.add_subplot(131))
    show_grid(numpy.log(numpy.maximum(1e-20,numpy.abs(pad_by_sel_(ideal.shape, selM, error_r)))) / numpy.log(10), "[r] log_10 (1-m)(n*AG)", N, axes=fig.add_subplot(132))
    show_grid(numpy.log(numpy.maximum(1e-20,numpy.abs(pad_by_sel_(ideal.shape, selM, error_r2)))) / numpy.log(10), "[r+d] log_10 (1-m)(n*AG)", N, axes=fig.add_subplot(133))
    fig = pylab.figure(figsize=(16,6))
    show_grid(numpy.log(numpy.maximum(1e-20,numpy.abs(fft(error)))) / numpy.log(10), "log_10 F[(1-m)(n*AG)]", N, axes=fig.add_subplot(121))
    show_grid(numpy.log(numpy.maximum(1e-20,numpy.abs(fft(error_r2b)))) / numpy.log(10), "[r+d] log_10 F[b*(1-m)(n*AG)]", N, axes=fig.add_subplot(122))
    pylab.show()
    #show_image(b, "b", 1)
    return
xs,ys = coordinates2(len(pswf))
pylab.rcParams['figure.figsize'] = 20, 16
pylab.contour(xs,ys, numpy.log(1/numpy.outer(pswf, pswf)) / numpy.log(10), levels=numpy.arange(0,12))
r = 0.4
pylab.gca().add_patch(patches.Circle((0,0), radius=r, fill=False))
pylab.gca().add_patch(patches.Polygon(r*numpy.transpose([numpy.cos(numpy.arange(6)/6*2*numpy.pi), numpy.sin(numpy.arange(6)/6*2*numpy.pi)]), True, fill=False))
pylab.colorbar()
Explanation: 2D case
Clearly this can be generalised to arbitrarily high dimensions. However, note that this makes things worse in a few ways at the same time:
Any overhead will be per image axis, so 2D overhead is squared 1D overhead
As we would use the outer product of the PSWF, the worst-case error is also squared (i.e. the corners are quite a bit worse).
There is unfortunately no good way around this, as the only way to significantly reduce overhead is to throw less of the image away, however this will increase exactly the error multiplier that we now have to square. As the corners are really the problem here, it might increase efficiency if we tile the plane into hexagons instead of squares...
End of explanation
<END_TASK>