metacal bias for crazy WCS #72
Comments
|
@esheldon for viz |
|
you should always add noise to the psf image when you will fit using LM (which the code is doing internally for Maybe try with |
|
I am not seeing psf fit failures but I can try this. |
|
Still biased with psf noise: and biased with |
|
I'd guess something to do with the col/row <-> x,y or u,v stuff. I always find that confusing. Have you checked if actually round objects come out measured as round when the wcs has a large g1 component? That seems easier to diagnose if that also shows a problem. |
|
Another thought -- if it's just a bug in the wcs handling, it should be insensitive to the size of the galaxy. Could check with hlr=3 instead of 0.5. If that works well, but 0.5 doesn't, then it's more likely something subtle with the relative pixel sizes in the two directions. In which case, going even smaller, but much higher S/N might be instructive. |
|
So an object with hlr=3 is still biased |
|
An object with hlr=0.1 is even more biased |
|
Oh crap. This is the moments fitter. |
|
Here is a Gaussian max like fit An earlier version of my shear test above showed metacal to be biased with this as well, which makes sense given the result above. So I think there is a bug in the Jacobian or pixel handling in ngmix somewhere. :/ |
|
Scripts that run these tests are in the same spot as above. |
|
I'm going to step away for now. I will probably start in on unit tests for ngmix tomorrow. |
|
Sounds like the first unit test to write is that Gaussian fit and moments give consistent answers when the wcs is sheared. |
|
OK. I started in on some unit tests for key parts. #73 A few results:
Edit: both of these tests are without PSF or pixel effects - that is next |
|
I did the main results of sheldon & huff with admom, so I'm a bit surprised. |
|
Well, metacalibration can calibrate anything. :) |
|
Also for certain cases, like no WCS distortions, it is closer to 2e-4. |
|
oh, I thought you mean 0.001 bias with metacal |
|
Ahhhh sorry! Yes, this is just in the shape that comes out. |
|
where do we stand on this one? |
|
Likely still a bug. |
|
To expand a bit. We have tests for almost all of the relevant WCS handling code and did not find any serious bugs. So I think there is a methodological error somewhere. |
|
Here is an updated script from Anna Niemiec that shows the issue. Note changing the psf to 'fitgauss' removes the bias. The issue is somewhere in the @rmjarvis is the original author of that import numpy as np
import ngmix
import galsim
def main(seed, psf='gauss', noise=1.e-6, ntrial=100):
print("ngmix version:", ngmix.__version__)
wcs = galsim.JacobianWCS(
-0.00105142719975775,
0.16467706437987895,
0.15681099855148395,
-0.0015749298342502371
)
print("WCS:", wcs.getDecomposition())
print()
print()
shear_true = [0.02, 0.00]
rng = np.random.RandomState(seed)
# We will measure moments with a fixed gaussian weight function
weight_fwhm = 1.2
fitter = ngmix.gaussmom.GaussMom(fwhm=weight_fwhm)
psf_fitter = ngmix.gaussmom.GaussMom(fwhm=weight_fwhm)
# these "runners" run the measurement code on observations
psf_runner = ngmix.runners.PSFRunner(fitter=psf_fitter)
runner = ngmix.runners.Runner(fitter=fitter)
# this "bootstrapper" runs the metacal image shearing as well as both psf
# and object measurements
boot = ngmix.metacal.MetacalBootstrapper(
runner=runner, psf_runner=psf_runner,
rng=rng,
psf=psf,
types=['noshear', '1p', '1m', '2p', '2m'],
)
dlist = []
for i in progress(ntrial, miniters=10):
im, psf_im, obs = make_data(
rng=rng, noise=noise, shear=shear_true, wcs=wcs
)
resdict, obsdict = boot.go(obs)
for stype, sres in resdict.items():
st = make_struct(res=sres, obs=obsdict[stype], shear_type=stype)
dlist.append(st)
print()
data = np.hstack(dlist)
# selections performed separately on each shear type
w = select(data=data, shear_type='noshear')
w_1p = select(data=data, shear_type='1p')
w_1m = select(data=data, shear_type='1m')
w_2p = select(data=data, shear_type='2p')
w_2m = select(data=data, shear_type='2m')
g = data['g'][w].mean(axis=0)
gerr = data['g'][w].std(axis=0) / np.sqrt(w.size)
g1_1p = data['g'][w_1p, 0].mean()
g1_1m = data['g'][w_1m, 0].mean()
# g2_1p = data['g'][w_1p, 1].mean()
# g2_1m = data['g'][w_1m, 1].mean()
# g1_2p = data['g'][w_2p, 0].mean()
# g1_2m = data['g'][w_2m, 0].mean()
g2_2p = data['g'][w_2p, 1].mean()
g2_2m = data['g'][w_2m, 1].mean()
R11 = (g1_1p - g1_1m)/0.02
R22 = (g2_2p - g2_2m)/0.02
# R12 = (g1_2p - g1_2m)/0.02
# R21 = (g2_1p - g2_1m)/0.02
shear = np.array([g[0] / R11, g[1]/R22])
shear_err = gerr / R11
m = np.linalg.norm(shear)/np.linalg.norm(shear_true)-1
merr = np.linalg.norm(shear_err)/np.linalg.norm(shear_true)
s2n = data['s2n'][w].mean()
print('S/N: %g' % s2n)
print('R11: %g' % R11)
print('m: %g +/- %g (99.7%% conf)' % (m, merr*3))
print('c: %g +/- %g (99.7%% conf)' % (shear[1], shear_err[1]*3))
print('shear 1 = %g +/- %g' % (shear[0], shear_err[0]))
print('shear 2 = %g +/- %g' % (shear[1], shear_err[1]))
return sres, im, psf_im
def make_struct(res, obs, shear_type):
"""
make the data structure
Parameters
----------
res: dict
With keys 's2n', 'e', and 'T'
obs: ngmix.Observation
The observation for this shear type
shear_type: str
The shear type
Returns
-------
1-element array with fields
"""
dt = [
('flags', 'i4'),
('shear_type', 'U7'),
('s2n', 'f8'),
('g', 'f8', 2),
('T', 'f8'),
('Tpsf', 'f8'),
]
data = np.zeros(1, dtype=dt)
data['shear_type'] = shear_type
data['flags'] = res['flags']
if res['flags'] == 0:
data['s2n'] = res['s2n']
# for moments we are actually measureing e, the elliptity
data['g'] = res['e']
data['T'] = res['T']
else:
data['s2n'] = np.nan
data['g'] = np.nan
data['T'] = np.nan
data['Tpsf'] = np.nan
# we only have one epoch and band, so we can get the psf T from the
# observation rather than averaging over epochs/bands
data['Tpsf'] = obs.psf.meta['result']['T']
return data
def select(data, shear_type):
"""
select the data by shear type and size
Parameters
----------
data: array
The array with fields shear_type and T
shear_type: str
e.g. 'noshear', '1p', etc.
Returns
-------
array of indices
"""
# raw moments, so the T is the post-psf T. This the
# selection is > 1.2 rather than something smaller like 0.5
# for pre-psf T from one of the maximum likelihood fitters
wtype, = np.where(
(data['shear_type'] == shear_type) &
(data['flags'] == 0)
)
w, = np.where(data['T'][wtype]/data['Tpsf'][wtype] > 1.2)
print('%s kept: %d/%d' % (shear_type, w.size, wtype.size))
w = wtype[w]
return w
def make_data(rng, noise, shear, wcs):
"""
simulate an exponential object with moffat psf
the hlr of the exponential is drawn from a gaussian
with mean 0.4 arcseconds and sigma 0.2
Parameters
----------
rng: np.random.RandomState
The random number generator
noise: float
Noise for the image
shear: (g1, g2)
The shear in each component
Returns
-------
ngmix.Observation
"""
psf_noise = 1.0e-8
stamp_size = 91
# psf_stamp = 71
# scale = 0.263
psf_fwhm = 0.9
gal_hlr = 0.5
psf = galsim.Moffat(
beta=2.5, fwhm=psf_fwhm,
).shear(
g1=0.02,
g2=-0.01,
)
obj0 = galsim.Exponential(
half_light_radius=gal_hlr,
).shear(
g1=shear[0],
g2=shear[1],
)
obj = galsim.Convolve(psf, obj0)
psf_im = psf.drawImage(nx=stamp_size, ny=stamp_size, wcs=wcs).array
im = obj.drawImage(nx=stamp_size, ny=stamp_size, wcs=wcs).array
psf_im += rng.normal(scale=psf_noise, size=psf_im.shape)
im += rng.normal(scale=noise, size=im.shape)
cen = (np.array(im.shape)-1.0)/2.0
psf_cen = (np.array(psf_im.shape)-1.0)/2.0
jacobian = ngmix.Jacobian(
x=cen[1], y=cen[0], wcs=wcs.jacobian(
image_pos=galsim.PositionD(cen[1], cen[0])
),
)
psf_jacobian = ngmix.Jacobian(
x=psf_cen[1], y=psf_cen[0], wcs=wcs.jacobian(
image_pos=galsim.PositionD(psf_cen[1], psf_cen[0])
),
)
wt = im*0 + 1.0/noise**2
psf_wt = psf_im*0 + 1.0/psf_noise**2
psf_obs = ngmix.Observation(
psf_im,
weight=psf_wt,
jacobian=psf_jacobian,
)
obs = ngmix.Observation(
im,
weight=wt,
jacobian=jacobian,
psf=psf_obs,
)
return im, psf_im, obs
def progress(total, miniters=1):
last_print_n = 0
last_printed_len = 0
sl = str(len(str(total)))
mf = '%'+sl+'d/%'+sl+'d %3d%%'
for i in range(total):
yield i
num = i+1
if i == 0 or num == total or num - last_print_n >= miniters:
meter = mf % (num, total, 100*float(num) / total)
nspace = max(last_printed_len-len(meter), 0)
print('\r'+meter+' '*nspace, flush=True, end='')
last_printed_len = len(meter)
if i > 0:
last_print_n = num
print(flush=True)
if __name__ == '__main__':
main(42, psf='fitgauss', noise=1.0e-8) |
|
Here is the current version of that. I probably introduced a bug recently, the original by @rmjarvis worked for distorted psfs ngmix/ngmix/metacal/metacal.py Line 481 in 51c488b |
|
If I do not reconvolve the inferred gaussian PSF by the pixel it works. |
|
I was summoned. But I'm not sure, do you need me to look at this? |
|
Thanks Mike. I'm finding that the algorithm you developed is failing for the above case. However it works if I don't reconvolve by the pixel. This is odd because in the code as originally developed the pixel is reconvolved. |
|
(And we do expect that reconvolving by the pixel is needed) |
|
I just copy pasted your code except for usin ngmix/ngmix/metacal/metacal.py Line 845 in 51c488b |
|
I should clarify, the algorithm isn't failing. We find a bias in the recovered shear with a complex WCS. And looking over your original tests, it didn't actually have a shear recovery test, it was about noise corrections. So it may be this is not a regression, it's just that we just never tested this algorithm with a distorted WCS. |
|
And moreover, the GalSim test you got this from (testing noise propagation) never worked properly. |
|
Still, metacal should really work with a moderately distorted WCS, so this is worth tracking down. I just wouldn't trust the code you got from me necessarily. |
|
@esheldon You have not touched the code in question in a significant way in two years or more. I wonder if this has always been there... |
|
Yeah, I expect it's not a regression. I've not touched it. Its just that we didn't test it in this regime. Note the code that wraps mike's algorithm shares all its other moving parts with the fitgauss method which works fine. |
|
I propose to modify the algorithm to instead work on the pixel convolved PSF rather than the pre-pixel PSF |
|
(which works for this case) |
|
That seems fine. Or we can make the algorithm more conservative and have it expand the gaussian more. I suppose those two are ultimately the same thing. |
|
I tried expanding, that doesn't work |
|
Huh? How can that be? Running it on the PSF with a pixel should be roughly the same as expanding the PSF. Are you sure there isn't something else going on? |
|
I didn't say I understand it. But when I put the pixel back in I get a bias. If I don't put it back, even if I do no dilation at all, it works fine. |
|
Did you try different PSF sizes? |
|
Yep |
|
More evidence What I'm doing is the following. I always draw the new psf image with method="no_pixel". This is because for all algorithms other than psf='gauss' the new psf already has the pixel in it. For psf='gauss' I needed to put the pixel back in for this to work. I get the bias. However, if I instead draw the PSF without setting method='no_pixel' this works fine. So it is some issue with the following
But this is ok
|
|
It may be interesting to find out why this is happening. But there is actually no reason to be working with a pre pixel PSF in this case. We are just trying to make a round PSF that is large enough, and that is best done with the pixel in. |
I am seeing some biases in metacalibration when I use a WCS Jacobian that differs a lot from a simple pixel scale.
Here is a basic test with a simple pixel scale
Now if I introduce a net
g1distortion in the WCS, I getHere is the test code
It only needs the
Momentsfitter from the metadetect repo.I have a copy in git here: https://github.com/beckermr/misc/blob/simple-des-y3/work/sheared_wcs_wl_sims/run.py
The text was updated successfully, but these errors were encountered: