tests of interpolation in lensing engine

[rmandelb]

We need to test more quantitatively the effects of interpolation on the shears/convergences that come out of the lensing engine for non-gridded positions. This is low priority from my perspective because it will not affect GREAT3, and that is where my attention is focused right now. I did assign this to myself since I'm willing to do it eventually, but if someone else wants to do this sooner than I am able, then by all means feel free to take this over.

[rmandelb]

Hi @barnabytprowe , @rmjarvis , and @gbernstein -

First, just like my pull request from yesterday, I completely understand if you can't get to this for a bit because of the holidays. I'm just posting it now while it's fresh in my mind.

I was going over the list of GalSim issues to see if there were any that might be good to do before we write up the GalSim paper, and I thought this one seemed like a good idea to address. I just noticed that I forgot to enable the commit-msg script on my new laptop, so the commits aren't showing up on here, sorry about that... but the summary is that I put together a piece of code that generates shears on a coarse grid and interpolates to a finer grid, generates shears on the finer grid, and compares the 2-point functions for the interpolated vs. direct generation of shears on the fine grid. You can check it out on this branch. I've done basic sanity checks and then started a big run (averaging over many random realizations) for 4 different interpolants. There is something about the results that I could use your advice on.

First I should say outright what my expectations were:

(1) For the power spectrum: based on discussions from a while ago, the picture I had in my head was that the interpolation acts like a smoothing function on the shear field, so the power spectrum should be multiplied by the square of the Fourier transform of the interpolant. This is something that we know analytically for our interpolants, so we can easily compare with that prediction, which (qualitatively) amounts to multiplying by 1 for small k and 0 for large k, where the turnover point depends on the kmax values for the coarse and fine grids and the specific interpolant choice.

(2) For the correlation function: no simple analytic prediction (it would come from doing the numerical integration over the modified PS). But it should be changed on small scales and then match up smoothly on large scales.

(3) 2-point functions can have finite-binning effects, but those will be the same for direct calculation vs. interpolated calculation of shears, so when we take ratios of interpolated vs. direct PS or correlation function, the finite-binning effects are largely irrelevant.

I realize now that this is a slightly naive picture, since there isn't a single PS or correlation function. e.g., if I put in EE power, then I could imagine some B-mode leakage (i.e., non-zero BB) that is not accounted for in this simple picture. Likewise, there could be some difference in what happens to xi+ vs. xi-.

So, here is an example of results for the power spectrum for EE, when I put in a cosmological shear power spectrum (I'm plotting dimensionless power, k^2 P/(2pi)) and use a linear interpolant.

In the top panel, I've plotted EE power, with blue being what I get with direct computation and red being the result with interpolation. The vertical lines show kmax for the grids with coarse and fine spacing. The bottom panel shows the ratio of (interpolated)/(direct) PS, where the red dashed line is the ideal ratio of 1 and the green dashed line is my naive prediction based on the FT of a linear interpolant. This is averaged over 2000 random realizations of the shear field.

Key points: (1) The power in the interpolated case is dropping above kmax for the coarse grid, though the scaling is not quite the naively predicted one. My initial guess is that this might relate to the fact that the naive prediction was calculated for finely spaced k, but we really have rather coarse bins, so the green curve sort of gets averaged over. Or my naive guess could be dumb, I suppose.

(2) There is a bump in the power below kmax. The power only slowly approaches the result from direct computation at the very lowest k / largest scales probed by our grid. This is a ~20% effect.

I wasn't very worried about (1), but I was worried about (2). So... that's when I started wondering about B-mode power. And indeed, there is some:

The bump in the red curve is massively significant, in case you're wondering about noise. (And i checked, there's no EB cross-power.) What's interesting to me is that if you consider the amplitude of that bump, it's comparable in order of magnitude to the EE power at 300<~ell<~2000 to within a factor of 2. These effects all depend on the interpolant - e.g., the enhanced EE power and the BB bump are about twice as pronounced for quintic interpolant than for linear.

As this message is rather long, I will spare you the figures of the correlation function, though I'm happy to post them later if anyone wants. For now I will summarize my conclusions: xi- acts like I expected, with suppression of the power on small scales (though the effect goes to larger scales than I expected). xi+ is actually enhanced on small scales by tens of percent.

My question is whether anyone can think of anything I've done wrong here, or if this seems correct? And if the latter, is there some simple physical explanation for why my naive analysis was wrong?

Note that I've not investigated the case of interpolation to non-gridded positions because then our handy power spectrum estimator does not work. I could do a correlation function analysis for that case, however. But this whole exercise has made me realize that the effect of interpolation on shears from the lensing engine is more complicated than I had previously realized.

[rmjarvis]

Interesting. I can't think of something specific that might be going on, but it does look like the B-mode plot is approximately equal to the excess power in the E-mode. (Is this true? It would be worth plotting them with the same axes to confirm.)

So my guess is that there is some extra shot noise coming into play, which adds to both E and B equally. Nothing jumps to mind where exactly this would be coming from, but that seems like the model for the effect you are seeing.

[rmandelb]

Meh, turns out they are not quite equal, though they are the same order. Here's a plot of the interpolated P_EE and P_BB with respect to the P_EE from the direct calculation: I was right in thinking that they were of similar orders of magnitude, but (interpolated P_EE/direct P_EE)-1 is about 1/3-1/2 the size of (interpolated P_BB/direct P_EE).

[rmandelb]

Sorry for the long delay since I updated this issue. Barney and I decided to work on it after the GREAT3 mid-challenge meeting which is now over.

So, @rmjarvis , I have an update and a question for you:

(1) Update: Based on discussion with Barney, it became apparent that the reason things weren't working is that it was a stupidly designed test. (He was nicer than that, but that's what it actually comes down to.) First, my input power spectrum had aliasing due to not being cut off at high k, and that's why I saw similar weird features for all interpolants - it was due to aliasing, not an interpolation effect at all. Second, the test is really an abuse of interpolation. One cannot expect to recover decent shear fields that are subsampled that much. We concluded that I should be doing two tests: one with an input grid, interpolated to an offset grid with the same size and spacing and hence the same min/max k. In this case, I can use the output power spectrum and correlation function as a test. For the second test, I would distribute points randomly, get shears on a grid with spacing determined by the typical separation between points, and interpolate to the random positions. In this case, we can compare the correlation function from the gridded points vs. for the randomly-distributed points. (I guess we will need to account for finite-gridding effects in the first instance.)

(2) Question: Is corr2 supposed to work with clang? I am trying to do this test on my new computer, and when I went to compile corr2 without OPENMP (since I thought that was a problem with clang++ when I tried recompiling GalSim), I got this kind of junk:

clang++ -I. -I/Users/rmandelb/software/cfitsio/include  -O3 -fno-strict-aliasing -Wall -Werror -g3 -c Corr2.cpp -o .odir/Corr2.o
In file included from Corr2.cpp:7:
In file included from ./BinnedCorr2.h:9:
./Cell.h:34:1: error: 'CellData' defined as a class template here but previously declared as a struct template [-Werror,-Wmismatched-tags]
class CellData<NData,M>
^
./Cell.h:31:1: note: did you mean class here?
struct CellData;
^~~~~~
class
1 error generated.
make: *** [.odir/Corr2.o] Error 1

Thanks, Rachel

[rmjarvis]

I probably haven't tested it with clang. So it doesn't surprise me that it balks at various things. It warns about things that other compilers don't care about.

You could try it without the -Wall -Werror part to see if that works. I should still fix the warnings, but if you make them not errors, it might get through the compilation and work for you.

[rmandelb]

Once I got rid of those compiler flags, it gets a bit further:

clang++ -I. -I/Users/rmandelb/software/cfitsio/include  -O3 -fno-strict-aliasing -g3 -c Corr2.cpp -o .odir/Corr2.o
clang++ -I. -I/Users/rmandelb/software/cfitsio/include  -O3 -fno-strict-aliasing -g3 -c CorrIO.cpp -o .odir/CorrIO.o
clang++ -I. -I/Users/rmandelb/software/cfitsio/include  -O3 -fno-strict-aliasing -g3 -c CalcT.cpp -o .odir/CalcT.o
CalcT.cpp:10:8: error: declaration conflicts with target of using declaration already in scope
double norm(double x) { return x*x; }
       ^
/Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin/../lib/c++/v1/complex:953:1: note: 
      target of using declaration
norm(double __re)
^
CalcT.cpp:9:12: note: using declaration
using std::norm;
           ^
1 error generated.
make: *** [.odir/CalcT.o] Error 1

If this is not a simple fix, then I will simply work on the cluster at CMU instead of on my laptop...

[rmjarvis]

In the short term, you can replace the line in CalcT.cpp

using std::norm;

with

double norm(const std::complex<double>& x) { return std::norm(x); }

and that should fix this.

Weirdly, I just tried clang++ on my machine and didn't get this error. But I did get the struct/class difference stuff that you were seeing with -Wall, which I fixed. It looks like the norm thing is a difference in Apple's implementation of the std::complex class, which I would venture is not actually in compliance with the standard. Shame on them. But easy to work around.

[rmandelb]

That works, TY!

[rmandelb]

Hi all, following up from the comment in https://github.com/GalSim-developers/GalSim/issues/387#issuecomment-31943307 -

I tried the test with interpolation from one grid to an offset grid. In constructing the offset grid, I was careful to eliminate any rows that ended up off the edge of the original grid (since the InterpolatedImage will just set those to zero). I have some results that surprised me in two ways, and further ideas that solved one of the problems but not the other. Here's a progress update on what worked, what didn't work, and some more questions:

The first result is for the EE power spectrum, where I compare the estimated PS from the original gridded points with the PS from the interpolated points (offset grid). Here is an example of the EE power for quintic interpolant, with 10k realizations:

That bump around ell~100-200 is real, not a noise artifact, and it shows up for all interpolants. Other than that, the power spectrum seems to be nearly correct for all k up to k_max (though perhaps there's a 1-2% error). For reference, the cubic interpolant is not quite as good, and the linear interpolant really stinks. Anyway, I am quite confused about that bump at ell~100-200.

But then I saw that I have some BB power: This is not an interpolation artifact, since it's there in the direct calculation. And it's not negligible compared to the low-k EE power, either. So now we have two problems: the bump around ell~100-200 and the BB power.

I figured out the bump in the sense that I understand where it comes from but I don't understand why. Since it was there for all interpolants, I thought maybe it was an artifact of my procedure of always cutting off the last row/column of the gridded shear field (since getShear() function will set all off-grid shears to zero). As a test, I took my original set of shears, simply cut off the last row/column without doing any interpolation at all, and lo and behold, the bump at ell~100-200 appeared. What I don't understand is why does this systematically give a jump in the power? For reference, ell~100 is about 3x the minimum ell corresponding to the entire grid extent, and I can't think of anything special about that.

To mitigate this problem for now, I put an option into galsim/lensing_ps.py to interpolate on the grid assuming it is periodic, so going off-grid corresponds to wrapping around to the other side of the grid.

Then I had another thought about a useful sanity check: I realize now that my theoretical line on those plots is wrong... it is probably relevant to the case where we interpolate to non-gridded points, but not to a grid. In fact, a useful sanity check is that if we use wrapped interpolation, and we use the nearest neighbor interpolant, we should exactly recover the original EE/BB/EB power since the PS estimation tool assumes a periodic grid anyway! So I used the NN interpolant with wrapping, and confirmed that my EE/BB/EB power are now identical before/after interpolation. What is interesting but not surprising is that the exact preservation of the PS does not translate to exact preservation of the correlation function; both xi_+ and xi_- are suppressed by ~1% on all scales, presumably because once you get to the grid edge, the shear correlations suddenly jump to those for 10 degrees instead of whatever separation you actually are. This means our recommendations to users about interpolations might depend on whether they are trying to recover the real-space or Fourier-space 2-point functions.

However, I haven't figured out the B mode power that I'm seeing. I tried changing the grid to increase kmax by a factor of a few, and to decrease kmin by a factor of a few; I softened the power-spectrum cutoff at kmax so it's no longer a hard cut. None of these nothing changed the BB power that I'm seeing even before interpolation. So I'm kind of stumped and feel I must be doing something stupid somewhere.

Remaining to-dos (updated based on subsequent comments):

• [x] put wrapped interpolation into the other getFoo methods
• [x] fix the B modes [it was due to kmin_factor]
• [x] try to understand why the grid cutoff was giving a bump in the power: grid cutoff was changing the k binning and number of modes in each bin
• [x] make option to cutoff edges properly vs. going to edge;
• [x] include proper wrapped interpolation for cubic/quintic/linear in GalSim lensing engine
• [x] make command-line option for wrapping interpolation in lensing engine, or not
• [x] turn more of the variables into command-line options
• [x] redo interpolation tests on grids with / without kmin_factor, with / without edge pixels; with / without wrapped interpolation
• [x] make a routine to check interpolation to non-gridded points
• [x] write calculation for correlation function including grid effects, to interpret results of the above. Or, try Mike's suggestions.
• [x] try Lanczos in addition to others, since Gary thinks it may do better than quintic
• [x] change order of loops, so all interpolants get tried at once
• [x] turn grid-sampling factor into a variable instead of being hard-coded
• [x] give lensing engine an option to zero out power above kmax internally, rather than forcing the user to do it
• [x] make FT interpolant utility accessible to users
• [x] give lensing engine method that will return the shear correlation functions for a given PS and grid configuration
• [x] simplify interpolant testing code to use the previous modifications in lensing engine
• [x] documentation: be careful about confusion between P(k) cutoff and grid cutoff, which are completely different things for which I used the same word.
• [x] put periodic wrapping for periodic interpolation and interpolant updates into the getFoo methods besides getShear (not needed for these tests)
• [x] edit docstrings to give guidance on grid configurations for interpolation
• [x] add warning about aliasing due to input P(k) having power above kmax,
• [x] edit docstrings: accuracy of interpolation of shear correlations and power spectra
• [x] clearly explain impact of cutting off grids or even edges of non-cutoff grids, tradeoffs in getting PS vs. CF right and E<->B mode leakage; periodic stuff
• [x] edit to clarify about limitations in the CF (integration over finite k range - unrelated to interpolation).
• [ ] more general documentation cleanup
• [ ] go over all changed code to see if any cleanup is needed
• [ ] go over new script to see if any cleanup is needed
• [ ] check for possibly broken unit tests / demos due to changes
• [ ] unit tests of new functionality
• [ ] check that all the refactoring didn't break something in the earlier part of this work
• [ ] CHANGELOG, etc.
• [ ] put stuff into config?

[rmjarvis]

Is it possible that this is related to the fact that we are using an interpolant that is using points off the edge? Quintic uses 3 points to each side of the given location. So anything within 3 pixels of the edge of the InterpolatedImage will have some spurious zeros being used.

Can you try making the grid 3 (or more) pixels larger than it needs to be and see if that has any effect?

[rmandelb]

Is it possible that this

Sorry, which of the effects in my comment are you trying to explain? There were a few different things that I highlighted and I'm not sure which one you are trying to explain.

[rmjarvis]

Any of it really, but the bump or the B mode I guess. If trimming a single row/col around the edge makes a difference, then it seems indicative that the Quintic implicitly using points off the edge would also matter.

[rmandelb]

The B mode is there before interpolation, so it cannot be related to the edge effects.

The bump shows up because of a single/row col elimination even without interpolation, so it also can’t be because of an effect in the quintic interpolant using points off the edge.

I’m not arguing with the idea that the quintic/cubic using points off the edge could matter for the performance of those interpolants, I just don’t see how they could be responsible for those two issues?

[rmandelb]

BTW, in terms of how we deal with periodic interpolation off the edge, I guess I would have to modify the code to not just wrap the new positions around, but also pad the grid of shears / convergences with a few rows from around the other side before making the SBInterpolatedImage?

[rmjarvis]

I see. Right. I guess that wouldn't make sense.

And yes, that sounds like the right plan if you want to implement periodic interpolation.

[rmandelb]

I have good news and bad news, Mike: I know exactly why we have B modes (and I don't really like the answer). So it's good that I now have a deeper understanding of what's going on but, well, it's an annoying result. I should say up front that it doesn't invalidate anything we've done for GREAT3 variable shear sims; nothing in how we're evaluating results is wrong, though our reduced sensitivity is partly due to the effect that I found.

Last night it occurred to me that I had been doing everything with a default GREAT3-like grid with kmin_factor=3. When I tested the impact of kmin_factor I had tried changing it to 5 (in case we needed an even lower kmin), but nothing other than that.

So, I confirmed this morning with a few experiments that if we use kmin_factor>1 and try to produce a purely E-mode shear field, then on the internally expanded grid, we achieve the desired result of P_E=input power and P_B=0, but on the central subgrid that we later take, there is no guarantee that P_B=0, and in fact this is precisely the source of the B modes that I'm seeing!

In retrospect this actually seems kind of obvious: if you generate a shear field on some grid and then take a significantly smaller sub grid, why would you expect the E/B power to be precisely preserved?

My conclusion is that if you have a grid of some extent and spacing, and you don't want to do anything significantly fancier than what we've done here, then you have two choices:

(1) Use kmin_factor=1 and get the right input E and B power, at the expense of reduced real-space shear correlations on scales corresponding to >~1/2 the grid extent.

(2) Use kmin_factor>1 to get more realistic real-space shear correlation properties. But then the power spectra have E<->B leakage.

I don't think there's a right answer; I think what you choose to do will have to depend on your application. For GREAT3 we went with (2), without quite realizing the choice we were making. After discussing with Barney, we agreed that this choice has the properties we want... for example, for variable PSF + variable shear we had been counting on having roughly realistic tradeoffs between the variable PSF and variable shear correlations, so if we had gone with choice (1) that would have been entirely lost. However, this new finding does explain the difficulties we were having with E/B mode separation, since it's likely that we had both E mode lensing shears leaking into B mode (from the kmin_factor effect and the fact that we use a sparsely sampled grid within each field) and B mode intrinsic shears leaking into E mode. Hence it's probably partly responsible for our reduced m sensitivity. There is no disaster here, however, since we compare the submissions to the generated intrinsic shear + lensing shear fields rather than the idealized input power spectra.

I will hazard a guess that had we known about this tradeoff in advance, we would have either done the same thing that we did without knowing about it, or we would have taken the additional time to cook up a smarter scheme (but still not gone with option (1)).

For people who plan on testing their ability to recover shear power spectra rather than correlations, option (1) seems to be the best one. The docstrings need to describe this tradeoff clearly, as well as to caution in general about the impact of generating shears on a grid and then only taking part of the grid.

Finally, regarding the "bump" I was seeing: Barney pointed out something kind of embarrassing - I may have accidentally been using different k binning when I removed a single row/column, so I was comparing apples to oranges. I need to confirm, but that could explain it. I'm updating my list of things to try to reflect all of the above!

[rmandelb]

Hello everyone who is still paying attention... I have an update on this issue. I'm pleased to report that it's one that features some answers instead of more of the "what is going on?" comments like those earlier on this page!

I implemented wrapped interpolation in the lensing engine (by which I mean that if I am doing "periodic interpolation", it wraps a layer of 5 edge pixels around to make a larger internal image to be used for the interpolation, rather than treating things off the grid as zeros). Once that was done, I could do a test of interpolation to shifted grids (not completely random points yet) that compares three cases:

(1) Interpolate to a shifted grid, in which case interpolation at the edges will include some 0's. (2) Interpolate to a shifted grid, but cut off the edge pixels. (3) Interpolate to a shifted grid using periodic interpolation, wrapping the grid around. (This seems silly if we care about correlation functions, but should be great for power spectra, I would imagine!)

As a reminder, the purpose of this somewhat artificial test was to get a sense for the basic effects of the interpolants and of edge effects before moving to random points (where we cannot even do a power spectrum analysis at all).

As for how I did this test: since my findings (described in the comments above) were that kmin_factor=3 fixes the suppressed shear correlation function on large scale but leads to B mode contamination, and kmin_factor=1 is necessary to get the power spectrum right, I did the comparison of the 3 cases both with kmin_factor=1 (in which case I only asked how well the power spectrum was preserved after interpolation) and with kmin_factor=3 (in which case I only asked how well the correlation function was preserved after interpolation).

Here are results for the kmin_factor=3 case where we are testing shear correlation function recovery after interpolation to another grid that is shifted by a random sub-pixel amount.

Case 1, no special treatment of edges: For the interpolant for which I expected the best behavior of those tested, quintic, I find that both xi_+ and xi_- are suppressed on all scales by ~2%. My interpretation of this is that near the edges, the interpolated shears included some points that were set to 0, and those edge points contribute to the shear correlations on all scales at some level. The results for case 2, shown below, support this interpretation. If I change the interpolant to cubic, then the results for xi_+ are similar to those for the quintic interpolant, but xi_- has some additional suppression on scales below ~0.6 degrees:

The linear interpolant is quite horrendous in this regard, crashing and burning below 1 degree for xi_-:

but again with decent results for xi_+. Finally the nearest neighbor interpolant behaves as well as quintic. I think this makes sense if we are interpolating to an offset grid, because in that case nearest neighbor just reproduces the same grid except with one set of pixels at the edge that was interpolated from off-grid and therefore is set to zero! (I'm sure NN will not do nearly as well when we interpolate to random points.)

Case 2, do the interpolation as usual but then chop off edge pixels since they are contaminated by interpolation that used zeros:

In this case, the quintic interpolant performs very nicely and the shear correlation functions are basically unaffected except for a slight suppression of xi_- by 2% below 0.5 degree. For cubic and linear, xi_+ remains fine and the small-scale suppression in xi_- seen with the quintic interpolant becomes increasingly worse (5%, 25%). And since we've cut off the edge pixels, the NN interpolant is perfect, of course. (I actually considered this as a sanity check of the results - if it wasn't true for an offset grid with edge pixels cut off, then it must be due to a bug.)

Case 3, do periodic interpolation: as I said above, I expected this would do great things for the power spectra but it's really a silly thing to do in principle if you just care about getting the correlation function right, and the results of this test suggest I was right. The results are similar to case 2 except with some suppression of the shear correlation functions (xi_+ AND xi_-) on all scales, which I think makes sense.

--> Basic conclusions from the kmin=3 case: If you wish to reproduce the shear correlation function properly, then you should chop off pixels near the edge (where the size of the region to be chopped depends on the interpolant). If you do this, then the quintic interpolant performs very nicely, whereas the cubic and linear interpolants lead to noticable suppression of xi_- but, interestingly, not xi_+. I confess I haven't thought hard enough about why only one of those would be affected and not the other. Comments / questions about all of this are welcome; my plan is that my findings will be summarized as a set of brief recommendations for how to use interpolation of shear fields in the lensing engine docstrings - after I extend this test to randomly-distributed points.

Will post about the kmin_factor=1 case and the effect of interpolation on the shear power spectra tomorrow... I think I've babbled enough for now!

[rmandelb]

And @rmjarvis , I have a quick question about some code I need to change in SBInterpolatedImage.cpp - I want to make sure I don't break something else.

If I use kmin_factor=1 then the mean shear will be very close to zero and occasionally it's so numerically close to zero that if I try to do interpolation using an SBInterpolatedImage, it will hit the exception:

            for (int iy = img.getYMin(); iy<= img.getYMax(); ++iy, ++y) {
                int x = xStart;
                for (int ix = img.getXMin(); ix<= img.getXMax(); ++ix, ++x) {
                    double value = img(ix,iy);
                    _pimpl->vx[i]->xSet(x, y, value);
                    sum += value;
                    sumx += value*x;
                    sumy += value*y;
                    xxdbg<<"ix,iy,x,y = "<<ix<<','<<iy<<','<<x<<','<<y<<std::endl;
                    xxdbg<<"value = "<<value<<", sums = "<<sum<<','<<sumx<<','<<sumy<<std::endl;
                }
            }
            if (sum == 0.) throw std::runtime_error(
                "Trying to initialize InterpolatedImage with an empty image.");

I would like to change the test near the end for empty images. One option I thought of is

if (sum==0. && sumx==0 && sumy==0) {

Do you think this will break something? Is there a better / more efficient test?

[rmjarvis]

I would like to change the test near the end for empty images.

I've recently gotten rid of that test, since it occasionally fails when using SBInterpolatedImage for things other than real GSObjects (such as here). It looks like this change is on branch #364, so the easiest thing is probably to just merge that branch into this one.

Unless you expect to want to PR this before that one is merged in, in which case I can try to find which commit it was and you can cherry pick it.

[rmandelb]

I've recently gotten rid of that test, since it occasionally fails when using SBInterpolatedImage for things other than real GSObjects (such as here).

I see. I think that in future the way forward is to make a 2d lookup table class as we discussed on another issue, but for now that is helpful.

I'm sure #364 will be done first, so I have no issues with merging it in here, and it looks like there are no major changes in lensing_ps.py on that branch that I would have to adapt to (unlike for #364a which had significant changes there).

[barnabytprowe]

Hi Rachel, very interesting and quite encouraging I think. Mainly encouraging that the results seem convergent as you increase the polynomial order of the interpolant: I'm sure the precise details of the degree of failure for linear vs cubic vs quintic etc. depends quite a lot on the power spectrum being used.

As to why xi+ seems easier than xi-...? Nothing occurs to me instantly. Will think about it...

[rmandelb]

Yes, I suspect you’re right that the results are power-spectrum dependent. I should be sure to mention that in the documentation. These results are with a cosmological power spectrum, which seemed like the best thing to use for such tests.

[barnabytprowe]

which seemed like the best thing to use for such tests.

Totally agree.

[rmandelb]

Now, I will post results for the kmin_factor=1 case, where we are testing the reproduction of the power spectra after interpolation to a shifted grid:

Case 1, no special treatment of edges: The interpolation doesn't introduce significant B modes, though they are there at the 0.5% level on some scales. For the E mode power spectrum, the quintic interpolant leads to slight suppression below kmax (~3%), increasing to ~20% for the linear interpolant. Also, as we decrease the order of interpolation, the k value at which results deterioriate also decreases - e.g., for the linear interpolant there are few % effects all the way down to kmin:

Case 2, cut off edges: The edge cutoff doesn't introduce significant B modes (well, they are there, but <~1% of the amplitude of E modes, usually less). So I will focus on the E mode power spectrum recovery after interpolation.

When we cut off the edges, the Delta k that is used for power spectrum estimation changes, which means I can't just take ratio of power spectrum for the original grid vs. that in the interpolated grid after chopping the edges. So, we have to consider, first, the question of "what does cutting off the edges do to the power spectrum before interpolation" and then "if I cut the edges off the original grid and the interpolated one, how well does it work?" For the first question, I show the power spectrum as observed (and the binned theory) before and after cutting off the edges in this plot:

There are no ratios shown because the k values at which the power is estimated are not the same. However, there are definitely small changes in the power due to the edge cutoff - I can read the numerical values from files and report them in the docstrings.

For the second question, how does the PS look after interpolation (when I've also cut the edges from the original and interpolated grid), here are the results for quintic:

It's quite good, just a slight suppression near kmax (~2%). This is worse for cubic and linear, both in the amplitude of the suppression at kmax and in the range of k values affected (the problem starts at progressively lower k as the order of interpolation is decreased). I think none of this is surprising given the results from the correlation function tests.

Case 3, periodic interpolation: The results are basically the same as the interpolation results for case 2, except that since there's no edge cutoff, there is not even a slight hint of B modes, or of E mode power changes due to the cutoff. Of course, NN interpolation is perfect (ratio of 1 for EE power after/before interpolation) and indeed this is one thing I checked since it must be true given the way the test was set up.

--> Summary: for power spectrum recovery, the best thing to do is use periodic interpolation, with a quintic interpolant. This prevents any B modes from showing up, and gives the best accuracy of interpolation up to kmax. Not using periodic interpolation results in E-B mode mixing at some low level, and cutting off the edges changes the power spectrum by a small but non-negligible amount.

Next week I'll do the tests with random positions.

[rmandelb]

Hi all - another update here. I did the test of interpolation to random points, and then realized (as seems to be very common for me with lensing engine tests) that the test itself may not be well-posed, or at least I didn't do a good job at formulating it. I could use some feedback on what exactly is the sensible question to ask, because of the following issues:

(1) It seems like the validity of interpolation might be related to something like the ratio of "typical separation between random points to which you want to interpolated" to "grid spacing". For now, I am using the same number of grid points and random points in the same area, which means that the typical random point separation and the grid spacing are very similar. But if I doubled the number of points then, with a random distribution in the same area, the typical separation between points gets smaller and it seems like we should have a harder time getting shear correlations right for a typical pair. I think that what I'm doing makes sense: use grid spacing ~ typical random point separation, and only check for recovery of shear correlations on scales from the grid spacing and up.

(2) If you set up the same bins for the correlation function estimates for the gridded and random points, it is nonetheless the case that the effective radius within each bin will differ, particularly for the small scales where only a few possible configurations go into each bin when you use the grid. So I need a way to control for gridding effects when comparing the results for the grid and random points. For power spectrum estimates we had a way to do this automatically: we would feed the power spectrum into the estimation code and ask what the output PS would be if we binned it into our bins that we were using for PS estimation. But I don't have a way to do this for correlation functions. If the answer is that the only way to do this is to do a similar direct calculation of such effects for the correlation function, then I will do it, but since it could be painful (given that we're just using corr2 for correlation functions) I wanted to see if anyone else has a better idea.

[barnabytprowe]

Good questions. Let's talk about this on the phone...

[rmandelb]

@rmjarvis , I have a question for you.

As described in https://github.com/GalSim-developers/GalSim/issues/387#issuecomment-34681039 part (2), I want a way to account for the different pair configurations that contribute within each bin in radius. I thought of 2 ways to do this:

(a) Add a simple brute-force correlation function estimator to the code and use it to average the theoretical correlation function over the pair separations in each bin. I did this, but it's too darn slow to be very functional.

(b) Then I started thinking I could just hack up corr2 to do this for me. My impression is that it might not be too hard to do this if I modify the code in CorrEE_multi.cpp, specifically the lines that accumulate bindata.xiplus and bindata.ximinus, so that instead of using the pair's ellipticities it uses some new function that knows the actual correlation function. Am I right that it could be as simple as inserting a new function that can retrieve the actual correlation function? (e.g., if I have tabulated it, it could read in the file and do some kind of interpolation) Or is it more complicated than that?

Of course if you have a better idea for approaching this problem overall, then I would be happy to hear it. Another idea I had was to forget about a theoretical calculation, just compare the correlation function from randomly distributed points vs. that from a very very fine grid (much finer than the grid we're actually using). But I worried that the memory required could be prohibitive, and I wouldn't have a great way of assessing whether it was fine enough to be trustworthy, which was why I had opted to try a theoretical calculation.

[rmjarvis]

I'm not really sure what you are trying to do. And I haven't looked through your new code on this branch at all. But ignorance of prior art shouldn't prevent me from having an opinion, so I'll have a go. :)

I would think that the easiest way to compare a gridded and interpolated correlation function to a "true" correlation function would be to have a gridded version and one that it gridded at say 10x smaller scales. They could both come from the same underlying power spectrum calculation, but the former would just take every 10th value in each grid direction.

ps_test = galsim.PowerSpectrum(...)
ps_true = galsim.PowerSpecturm([... same thing ...])
ps_test.buildGrid(grid_spacing = dx, ngrid = ngrid, kmax_factor = 10, ...)
ps_true.buildGrid(grid_spacing = dx/10, ngrid = ngrid*10, ...)

Then internally, the two will have the same call to PowerSpectrumRealizer, so you have the same effect from whatever the fft does. But the test one will have a wider grid than the true one. So any differences must be due to interpolation effects.

The corr2 program can easily handle a million points (should take around a minute), and memory and running time scale linearly with the number of points. So you can use this when deciding how finely you want to make the gridding.

BTW, the file CorrEE_multi.cpp is a legacy file from the old version of the code. It is not used at all in the current version of corr2.

[rmandelb]

I think this example you gave corresponds to interpolating by a pretty huge factor, which I'm not sure is a sensible use case. I would assume that most people would say something like "I have 10 galaxies/arcmin^2, so I will make a regular grid that corresponds to roughly the same density of points." Do you think I am wrong and that it is helpful to test interpolation in this other mode of far denser point spacing?

Also, it means we would never do a test of interpolating shears to random points, which doesn't seem like a great idea to me. Given that we define shears on a grid using a DFT approach, it seems like having a test that goes off-grid is important, isn't it?

Since you say that you've not been following along very closely, I should point out that I did already do a test of interpolating to offset grids (i.e., make one grid with some spacing, and interpolate to another grid with the same spacing that is offset by some amount).

Also, my concern about timing and memory was for the generation of the shears via FFT, not the corr2 calculation. Once we're doing something like kmin_factor=3 then if we have a 100x100 grid (=300x300 internally) and we start subsampling by a very large factor, the FFTs get enormous.

[rmjarvis]

I was suggesting that you use random locations, but get the interpolated g1,g2 values from both ps objects. Technically, both are interpolated, but one of them is sampled much more finely than the other, so for practical purposes, we can treat those values as "true".

If you think it would be helpful to chat, you can ping me on either skype or google hangout. I'll be available for the next hour or so.

[rmandelb]

Oh, I see what you mean now. I guess as long as I use the same true PS and k range for the ps objects, that could be a good test (though expensive).

[rmjarvis]

You don't need to use 10x of course. That might be overkill. I suspect the interpolation errors scale as grid_scale**n where n might be 4 or higher. (Probably 4 for Cubic and 6 for Quintic? Just a guess.) So even a factor of 5 or 3 might be sufficient to get a good comparison.

And also, you can use a lot of sample points, so you can beat down the effect of the finite number of points pretty well. 10 million is probably reasonable.

[gbernstein]

Some delayed comments on this, which i hope are not too irrelevant given my total ignorance of the code in question.

On Feb 10, 2014, at 3:54 PM, rmandelb notifications@github.com wrote:

Hi all - another update here. I did the test of interpolation to random points, and then realized (as seems to be very common for me with lensing engine tests) that the test itself may not be well-posed, or at least I didn't do a good job at formulating it. I could use some feedback on what exactly is the sensible question to ask, because of the following issues:

(1) It seems like the validity of interpolation might be related to something like the ratio of "typical separation between random points to which you want to interpolated" to "grid spacing". For now, I am using the same number of grid points and random points in the same area, which means that the typical random point separation and the grid spacing are very similar. But if I doubled the number of points then, with a random distribution in the same area, the typical separation between points gets smaller and it seems like we should have a harder time getting shear correlations right for a typical pair. I think that what I'm doing makes sense: use grid spacing ~ typical random point separation, and only check for recovery of shear correlations on scales from the grid spacing and up.

You would not expect to ever be able to reproduce the correlation function on scales below the grid pitch (or even comparable to it) by interpolating between the grid. The reason is that the gridded realization necessarily truncates the power spectrum at the nyquist frequency (or aliases it), and only has power at the grid in k space. A perfect interpolant will produce produce a densely sampled image having power only at these k values present in the gridded realization, so its correlation function won't match the one you get from th full continuous parent power spectrum. I'd guess that random interpolation to any point density will work but must be compared to the ft of the sampled power spectrum realized on the grid.

(2) If you set up the same bins for the correlation function estimates for the gridded and random points, it is nonetheless the case that the effective radius within each bin will differ, particularly for the small scales where only a few possible configurations go into each bin when you use the grid. So I need a way to control for gridding effects when comparing the results for the grid and random points. For power spectrum estimates we had a way to do this automatically: we would feed the power spectrum into the estimation code and ask what the output PS would be if we binned it into our bins that we were using for PS estimation. But I don't have a way to do this for correlation functions. If the answer is that the only way to do this is to do a similar direct calculation of such effects for the correlation function, then I will do it, but since it could be painful (given that we're just using corr2 for correlation functions) I wanted to see if anyone else has a be! tter idea .

—

One possibility ( maybe this is what you mean already?) is to take each pair in your random points and calculate the theoretical 2pt function at this separation. Then bin the theory and the random values in the same bins, insuring you have the same effective radius in each bin. Though the caveat above would still apply, namely that the "theory" would have to be derived from a sampled version of the ps if you want to expect perfect agreement.

The above may be only marginally comprehensible, my apologies. G

Reply to this email directly or view it on GitHub.

[rmandelb]

Hi Gary - thanks for chiming in. To respond to your two points:

You would not expect to ever be able to reproduce the correlation function on scales below the grid pitch (or even comparable to it) by interpolating between the grid.

Right. That's why I'm assuming one would use the density of random points as a guide for selecting the grid pitch. But once you've made some choice, it seems like a good idea to quantify how well the different interpolants do / how close can we approach to getting the correlation function right on scales near the grid pitch (where the "true" correlation function we should compare to is one that comes from integrating the power spectrum from the original grid kmin to kmax, not from 0<k<infinity). In other words I think we must be clear in the documentation that there are two aspects to why interpolating between gridded points will modify the correlation function: the first is that using a grid necessarily means that we are truncating the power spectrum at some min/max k, and the second is that the interpolant will cause deviations in the correlation function even on scales above the grid pitch. I've been trying to design tests to quantify the 2nd point, since I assume the user can quantify the first point simply by comparing the ideal theoretical correlation function vs. the one for their chosen grid parameters.

So I think we're on the same page, right? Or is there some aspect to this that I'm missing?

Regarding your second point:

One possibility ( maybe this is what you mean already?)

Yes, that possibility you raise is precisely what I had in mind. I just couldn't come up with a computationally feasible way to do it so I was trying to weasel out of it. :) The naive way that I thought of initially is super slow (I checked) and the alternative would involve hacking Mike's corr2 code, which did not appeal to me...

Mike's suggestion may provide a way around the problem.

[rmandelb]

(The formatting of my comment just above came out funny, so I had to edit it - you may need to view on github rather than via e-mail. There are no changes in content.)

[gbernstein]

On Feb 16, 2014, at 11:23 PM, rmandelb notifications@github.com wrote:

Hi Gary - thanks for chiming in. To respond to your two points:

You would not expect to ever be able to reproduce the correlation function on scales below the grid pitch (or even comparable to it) by interpolating between the grid.

Right. That's why I'm assuming one would use the density of random points as a guide for selecting the grid pitch. But once you've made some choice, it seems like a good idea to quantify how well the different interpolants do / how close can we approach to getting the correlation function right on scales near the grid pitch (where the "true" correlation function we should compare to is one that comes from integrating the power spectrum from the original grid kmin to kmax, not from 0<k<infinity). In other words I think we must be clear in the documentation that there are two aspects to why interpolating between gridded points will modify the correlation function: the first is that using a grid necessarily means that we are truncating the power spectrum at some min/max k, and the second is that the interpolant will cause deviations in the correlation function even on scales above the grid pitch. I've been trying to design tests to quantify the 2nd point, since I a! ssume the user can quantify the first point simply by comparing the ideal theoretical correlation function vs. the one for their chosen grid parameters.

So I think we're on the same page, right? Or is there some aspect to this that I'm missing?

Yes, totally agreed. The result you expect to get is that the ps of the random points is the sampled/truncated ps times the ft of the interpolant squared. The 2pcf is of course the ft of this, which convolves the interpolant with a band limited, folded version of the parent function. The quintic we have is optimized for a different purpose and you might find a lanczos of the same size doing a slightly better job near the grid scale. Regarding your second point:

One possibility ( maybe this is what you mean already?) Yes, that possibility you raise is precisely what I had in mind. I just couldn't come up with a computationally feasible way to do it so I was trying to weasel out of it. :) The naive way that I thought of initially is super slow (I checked) and the alternative would involve hacking Mike's corr2 code, which did not appeal to me...

Mike's suggestion may provide a way around the problem.

— Reply to this email directly or view it on GitHub.

[rmandelb]

The result you expect to get is that the ps of the random points is the sampled/truncated ps times > the ft of the interpolant squared.

Yes, that's what I was expecting. Interesting point re: quintic vs. Lanczos, I can add the Lanczos with varying n to my list of interpolants to try.

[rmandelb]

I have some results that used the test that Mike proposed. As a summary, that test involves:

(1) Making a grid that has the same physical extent as the one we want, but has a 10x smaller grid spacing, or 100 times as many points. So, this grid has spacing of 0.01 degrees, with 1000x1000 points, giving a total extent of 10x10 deg^2. (2) Interpolate from that grid to random positions. This will serve as our reference case. (3) Take every 10th point in the original grid, to achieve the grid spacing we originally wanted, which was 0.1 degrees. (4) Interpolate from that grid to random positions. This will serve as our test case.

While our reference is not ideal in the sense that it includes interpolation, the scale at which the interpolant should become important is 10x smaller than the scale at which it should become important for our test case, so there is a wide range of scales where this is a good test of interpolation. The interpolants I tried were nearest, linear, cubic, quintic, lanczos3, lanczos5. I did not test how sensitive these are to the input power spectrum, just used a cosmological E-mode power spectrum (no B-mode power).

Here is a summary of my conclusions.

The base case was to do the interpolation and not worry about any edge effects in the interpolation. In that case, here's an example of what the results look like for the quintic interpolant, where the first plot is xi+ and the second is xi_-:

If I compare the interpolants, for xi_m my conclusions are that:

• nearest and linear lead to ~3% suppression from 10-30x the grid spacing, and give 10s of % suppression below 10x the grid spacing.
• cubic gives ~2% suppression from 4-30x the grid spacing, then tanks below 4x the grid spacing
• quintic gives ~1% suppression from 3-30x the grid spacing, then tanks below 3x the grid spacing
• lanczos3 gives ~2% suppression from 10-30x the grid spacing, up to a 5% enhancement at 3x the grid spacing, then tanks below 2x the grid spacing. lanczos5 is quite similar in behavior.

For xi_p, the situation is a bit weirder.

• nearest and linear interpolants gives a major enhancement below 3x the grid spacing, then a ~2% enhancement from 2x-30x the grid spacing, then a weird dip of 10% around 30x the grid spacing.
• cubic, quintic, and lanczos3,5 are quite similar to the above in all ways, except that there is a strange additional 15% dip around 2x the grid spacing.

Then I considered the case of cutting off the edges, since we expect the interpolant to bring in some zeros from off the edge of the grid, which should suppress the shears near the edges. I was wondering if that could be the cause of the extended low-level suppression we saw in the previous results. And indeed that seems to be the case - for example, for xi_-, the ~3% suppression at 10-30x the grid spacing has gone away now that we cut off the edges, though of course the larger suppression that we saw on smaller scales is still there. In all, I think that the best-behaved interpolant once we cut off the edges is lanczos5:

It still has that funny little peak around 3x the grid spacing and then drops like a rock below 2x the grid spacing, but otherwise it's quite good.

For xi_+, once we cut off the edges, some of the oddities remain though they are not huge. Below are examples for quintic, lanczos5:

The quintic is slightly better but they both have similar shape.

I think the high-level recommendations for users that comes out of this work are:

(1) If you want shears at random positions, and it's important to you that those shear correlations should be correct (i.e.: reflect the input power spectrum sampled using the DFT method from kmin->kmax) above some particular separation theta, then you should choose a grid spacing for your grid that is <~theta/3 of that scale in order to avoid a significant error due to the interpolation.

(2) Ideally the user should choose a grid configuration that extends outside the region where they want to do the interpolation, so that they can stay away from the edges of the grid.

(3) We recommend the lanczos5 interpolant, which should be set as the default. Quintic and lanczos3 are not much worse. Nearest/linear/cubic are not good.

I did have one crazy idea which I didn't try (would be happy to get feedback on this): since we think we know what the interpolant is doing, namely multiplying the input shear PS by the square of the FT of the interplant, would it make any sense for a user who knows they'll be interpolating to start by dividing their shear PS by the square of the FT of the interpolant, boosting it up in a way that will counteract the effects of the interpolant? It means the shear PS and correlation function for the gridded points will be wrong, but if the user only cares about getting it right for the random points, then could this do the trick? We could perhaps include a kwarg for this if it makes any sense, but I'm not sure it does.

Any other comments / questions / suggestions? If not, my remaining action items on this branch are some code cleanup and lots of documentation updates, no further tests to run.

[rmjarvis]

since we think we know what the interpolant is doing, namely multiplying the input shear PS by the square of the FT of the interplant, would it make any sense for a user who knows they'll be interpolating to start by dividing their shear PS by the square of the FT of the interpolant, boosting it up in a way that will counteract the effects of the interpolant?

I thought of this too. Are there any zero crossings in the interpolants' power spectra that might be problematic?

I think at least we probably ought to provide a method in the Interpolant classes that returns the power spectrum so users can apply this correction if they want to.

[rmjarvis]

Oh, and I think the conclusions and recommendations you stated above make sense.

When someone gets around to doing #465, I think that is where this default belongs. The PowerSpectrum class should use the LookupTable2D internally (or LookupImage or whatever we end up calling it), and that class should have Lanczos5 as the default interpolant.

[rmandelb]

I will check for zero crossings in the interpolants, but that’s a good point… I think there are in some cases.

As for a method to get the power spectrum, I agree: I coded one up by hand in this script, but one of my remaining action items is to turn some of these utilities that I invented for this issue into user-accessible code. (e.g., I also want to make PowerSpectrum objects capable of returning a theory prediction for xi+ and xi- for the given grid’s kmin and kmax. It still won’t be quite right because it ignores the fact that discrete vs. continuous FT issue, but it’s better than assuming it’s the xi that you get when you integrate the PS from 0 to infinity!)

[rmandelb]

I agree that the default ultimately belongs in the LookupTable2D (or whatever) class, but my plan was to have a PR for updated code on this branch that sets that as a default in lensing_ps.py so that the default does the right thing even before we create that new class.

On Feb 21, 2014, at 10:41 AM, Mike Jarvis notifications@github.com wrote:

Oh, and I think the conclusions and recommendations you stated above make sense.

When someone gets around to doing #465, I think that is where this default belongs. The PowerSpectrum class should use the LookupTable2D internally (or LookupImage or whatever we end up calling it), and that class should have Lanczos5 as the default interpolant.

— Reply to this email directly or view it on GitHub.

---

Rachel Mandelbaum rmandelb@andrew.cmu.edu http://www.andrew.cmu.edu/~rmandelb/

[rmjarvis]

my plan was to have a PR for updated code on this branch that sets that as a default in lensing_ps.py

Agreed. My suggestion was just for the eventual point when we do #465, with hash tags so that issue gets callback references to this discussion. :)

[rmandelb]

Closing #387.