Hello,

I tried running the metacal example (using the bootstrap solution) from the Wiki page and I got the following error on a lot of images :

TypeError: __init__() missing 1 required positional argument: 'value'

This error seems to come from scipy.optimize.leastsq function..

Here are the version of the different libraries used :
python = 3.6.8
scipy = 1.2.0
ngmix = v1.2 (clone from the master branch)
galsim = 2.1.4

For example, currently you can shear a shape in place using

sh=ngmix.Shape(0.2, 0.3)
sh.shear(0.05, 0.0)

instead that should do

shnew = sh.shear(0.05, 0.0)

and we could replace the old version with

sh.shear_inplace(..)

There are a lot of other places in the code where references are passed around and this is hard to reason about. For example, the bootstrappers take in observations and sometimes mutate them.

@rgruendl ran into an ngmix error when using it to measure the moments of Piff models. At least for the example he posted, this is clearly a bad Piff solution with both positive and negative lobes. However, rather than either computing some values for the moments or returning an error flag, it entered an infinite loop and died from a stack overflow. In my test it gives a RecursionError from exceeding maximum recursion depth.

import piff
import galsim
import ngmix
import numpy as np

# Test a Piff model that was giving errors when doing the ngmix fit.

# https://desar2.cosmology.illinois.edu/DESFiles/desarchive/ACT/finalcut/Y6A2_PIFF/20141228-r5286/D00392686/p01/psf/D00392686_z_c43_r5286p01_piff-model.fits
fname = 'D00392686_z_c43_r5286p01_piff-model.fits'

# Corresponding image is here, but we don't need it for this test script.
# https://desar2.cosmology.illinois.edu/DESFiles/desarchive/OPS/finalcut/Y5A1/r3567/20141228/D00392686/p02/red/immask/D00392686_z_c43_r3567p02_immasked.fits.fz

psf=piff.read(fname)

ibad = 63
ccdnum = 43
bad_star = psf.stars[ibad]
x = bad_star.image_pos.x
y = bad_star.image_pos.y
flux = bad_star.flux

stamp_size=24
b=galsim.BoundsI(int(x)-stamp_size/2, int(x)+stamp_size/2,
                 int(y)-stamp_size/2, int(y)+stamp_size/2)
model = galsim.ImageF(bounds=b, wcs=psf.wcs[ccdnum].local(bad_star.image_pos))

psf.draw(x=x, y=y, chipnum=ccdnum, flux=flux, image=model)

# The following is excerpted from 
# https://github.com/DarkEnergySurvey/despyPIFF/blob/main/python/despyPIFF/piff_qa_utils.py
# with some simplifying modifications by MJ

#
# Setting up priors
#

fwhm = 0.9  # Original version did a calculation here.  But that doesn't seem to matter.
flag = 0
dx, dy, g1, g2, flux = 0., 0., 0., 0., 0.
T_guess = (fwhm/ 2.35482)**2 * 2.
T = T_guess

cen = model.true_center - model.origin
pixel_scale = 0.263
seed = 1234  # Doesn't seem to matter
rng=np.random.RandomState(seed)

cen_prior=ngmix.priors.CenPrior(0.0,0.0,pixel_scale,pixel_scale,rng=rng)
gprior=ngmix.priors.GPriorBA(0.1,rng=rng)
Tprior=ngmix.priors.LogNormal(T,0.2,rng=rng)
Fprior=ngmix.priors.FlatPrior(-10.,1.e10,rng=rng)
prior=ngmix.joint_prior.PriorSimpleSep(cen_prior,gprior,Tprior,Fprior)

#
# Putting data in context for NGMIX
#

jac=ngmix.Jacobian(wcs=model.wcs, x=cen.x + x - int(x+0.5), y=cen.y + y -int(y+0.5))
# Original used the weight, but that also doesn't matter
obs=ngmix.Observation(image=model.array, weight=None, jacobian=jac)

psf_fitter = ngmix.fitting.Fitter(model='gauss', prior=prior)
psf_guesser = ngmix.guessers.SimplePSFGuesser(rng=rng, guess_from_moms=True)

# This fails with infinite loop.
psf_runner = ngmix.runners.PSFRunner(fitter=psf_fitter, guesser=psf_guesser, ntry=3)
res = psf_runner.go(obs=obs)

print(res)
ngmix_flag=res['flags']
assert ngmix_flag == 0

On my system, this outputs the following:

Traceback (most recent call last):
  File "ngmix_bug.py", line 69, in <module>
    res = psf_runner.go(obs=obs)
  File "/opt/anaconda3/envs/py3.8/lib/python3.8/site-packages/ngmix/runners.py", line 110, in go
    return run_psf_fitter(
  File "/opt/anaconda3/envs/py3.8/lib/python3.8/site-packages/ngmix/runners.py", line 207, in run_psf_fitter
    res = run_fitter(
  File "/opt/anaconda3/envs/py3.8/lib/python3.8/site-packages/ngmix/runners.py", line 141, in run_fitter
    guess = guesser(obs=obs)
  File "/opt/anaconda3/envs/py3.8/lib/python3.8/site-packages/ngmix/guessers.py", line 1034, in __call__
    return self._get_guess(obs=obs)
  File "/opt/anaconda3/envs/py3.8/lib/python3.8/site-packages/ngmix/guessers.py", line 1038, in _get_guess
    T, flux = self._get_T_flux(obs=obs)
  File "/opt/anaconda3/envs/py3.8/lib/python3.8/site-packages/ngmix/guessers.py", line 774, in _get_T_flux
    T, flux = self._get_T_flux_from_moms(obs=obs)
  File "/opt/anaconda3/envs/py3.8/lib/python3.8/site-packages/ngmix/guessers.py", line 818, in _get_T_flux_from_moms
    T, flux = self._get_T_flux(obs=obs)
  File "/opt/anaconda3/envs/py3.8/lib/python3.8/site-packages/ngmix/guessers.py", line 774, in _get_T_flux
    T, flux = self._get_T_flux_from_moms(obs=obs)
  File "/opt/anaconda3/envs/py3.8/lib/python3.8/site-packages/ngmix/guessers.py", line 818, in _get_T_flux_from_moms
    T, flux = self._get_T_flux(obs=obs)
  File "/opt/anaconda3/envs/py3.8/lib/python3.8/site-packages/ngmix/guessers.py", line 774, in _get_T_flux
    T, flux = self._get_T_flux_from_moms(obs=obs)
  File "/opt/anaconda3/envs/py3.8/lib/python3.8/site-packages/ngmix/guessers.py", line 818, in _get_T_flux_from_moms
    T, flux = self._get_T_flux(obs=obs)

... [snip] (Many more of these recursions)

  File "/opt/anaconda3/envs/py3.8/lib/python3.8/site-packages/ngmix/guessers.py", line 774, in _get_T_flux
    T, flux = self._get_T_flux_from_moms(obs=obs)
  File "/opt/anaconda3/envs/py3.8/lib/python3.8/site-packages/ngmix/guessers.py", line 818, in _get_T_flux_from_moms
    T, flux = self._get_T_flux(obs=obs)
  File "/opt/anaconda3/envs/py3.8/lib/python3.8/site-packages/ngmix/guessers.py", line 774, in _get_T_flux
    T, flux = self._get_T_flux_from_moms(obs=obs)
  File "/opt/anaconda3/envs/py3.8/lib/python3.8/site-packages/ngmix/guessers.py", line 806, in _get_T_flux_from_moms
    res = wt.get_weighted_moments(obs=obs, maxrad=1.0e9)
  File "/opt/anaconda3/envs/py3.8/lib/python3.8/site-packages/ngmix/gmix/gmix.py", line 626, in get_weighted_moments
    res = self.get_weighted_sums(obs, maxrad)
  File "/opt/anaconda3/envs/py3.8/lib/python3.8/site-packages/ngmix/gmix/gmix.py", line 656, in get_weighted_sums
    gmix_nb.get_weighted_sums(
  File "/opt/anaconda3/envs/py3.8/lib/python3.8/site-packages/numba/core/dispatcher.py", line 368, in _compile_for_args
    argtypes.append(self.typeof_pyval(a))
  File "/opt/anaconda3/envs/py3.8/lib/python3.8/site-packages/numba/core/dispatcher.py", line 686, in typeof_pyval
    tp = typeof(val, Purpose.argument)
  File "/opt/anaconda3/envs/py3.8/lib/python3.8/site-packages/numba/core/typing/typeof.py", line 31, in typeof
    ty = typeof_impl(val, c)
  File "/opt/anaconda3/envs/py3.8/lib/python3.8/functools.py", line 875, in wrapper
    return dispatch(args[0].__class__)(*args, **kw)
  File "/opt/anaconda3/envs/py3.8/lib/python3.8/site-packages/numba/core/typing/typeof.py", line 139, in _typeof_numpy_scalar
    return numpy_support.map_arrayscalar_type(val)
  File "/opt/anaconda3/envs/py3.8/lib/python3.8/site-packages/numba/np/numpy_support.py", line 213, in map_arrayscalar_type
    return from_dtype(dtype)
  File "/opt/anaconda3/envs/py3.8/lib/python3.8/site-packages/numba/np/numpy_support.py", line 93, in from_dtype
    return from_struct_dtype(dtype)
  File "/opt/anaconda3/envs/py3.8/lib/python3.8/site-packages/numba/np/numpy_support.py", line 520, in from_struct_dtype
    ty = from_dtype(elemdtype)
  File "/opt/anaconda3/envs/py3.8/lib/python3.8/site-packages/numba/np/numpy_support.py", line 111, in from_dtype
    return types.NestedArray(subtype, dtype.shape)
  File "/opt/anaconda3/envs/py3.8/lib/python3.8/site-packages/numba/core/types/abstract.py", line 66, in __call__
    inst = type.__call__(cls, *args, **kwargs)
  File "/opt/anaconda3/envs/py3.8/lib/python3.8/site-packages/numba/core/types/npytypes.py", line 550, in __init__
    super(NestedArray, self).__init__(dtype, ndim, 'C', name=name)
  File "/opt/anaconda3/envs/py3.8/lib/python3.8/site-packages/numba/core/types/npytypes.py", line 432, in __init__
    super(Array, self).__init__(dtype, ndim, layout, name=name)
  File "/opt/anaconda3/envs/py3.8/lib/python3.8/site-packages/numba/core/types/common.py", line 51, in __init__
    if isinstance(dtype, Buffer):
  File "/opt/anaconda3/envs/py3.8/lib/python3.8/abc.py", line 98, in __instancecheck__
    return _abc_instancecheck(cls, instance)
RecursionError: maximum recursion depth exceeded in comparison

Related to #82

lots of code duplication in _gmix.c

I think this code here generates an infinite loop for large guess widths

https://github.com/esheldon/ngmix/blob/master/ngmix/guessers.py#L308

@beckermr

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

$ python run.py 10 0.0 0.0
# of sims: 10
wcs_g1   : 0.000000
wcs_g2   : 0.000000
dudx     : 0.263000
dudy     : -0.000000
dvdx     : -0.000000
dvdy     : 0.263000
m [1e-3] : 4.628283 +/- 4.470836
c [1e-4] : -0.508443 +/- 1.581002

Now if I introduce a net g1 distortion in the WCS, I get

$ python run.py 10 0.1 0.0
# of sims: 10
wcs_g1   : 0.100000
wcs_g2   : 0.000000
dudx     : 0.237892
dudy     : -0.000000
dvdx     : -0.000000
dvdy     : 0.290757
m [1e-3] : -197.869263 +/- 5.030609
c [1e-4] : -0.769631 +/- 1.390868

Here is the test code

import sys
import numpy as np
import tqdm

import ngmix
import galsim
import numba

from ngmix import metacal
from metadetect.fitting import Moments


def run_metacal(*, n_sims, wcs_g1, wcs_g2):
    """Run metacal and measure m and c.

    The resulting m and c are printed to STDOUT.

    Parameters
    ----------
    n_sims : int
        The number of objects to simulated.
    wcs_g1 : float
        The shear on the 1-axis of the WCS Jacobian.
    wcs_g2 : float
        The shear on the 2-axis of the WCS Jacobian.
    """
    jc = galsim.ShearWCS(0.263, galsim.Shear(g1=wcs_g1, g2=wcs_g2)).jacobian()

    jacobian_dict = {
        'dudx': jc.dudx,
        'dudy': jc.dudy,
        'dvdx': jc.dvdx,
        'dvdy': jc.dvdy
    }

    swap_g1g2 = False

    res = _run_metacal(
        n_sims=n_sims,
        rng=np.random.RandomState(seed=10),
        swap_g1g2=swap_g1g2,
        **jacobian_dict)

    g1 = np.array([r['noshear']['g'][0] for r in res])
    g2 = np.array([r['noshear']['g'][1] for r in res])
    g1p = np.array([r['1p']['g'][0] for r in res])
    g1m = np.array([r['1m']['g'][0] for r in res])
    g2p = np.array([r['2p']['g'][1] for r in res])
    g2m = np.array([r['2m']['g'][1] for r in res])

    g_true = 0.02
    step = 0.01

    if swap_g1g2:
        R11 = (g1p - g1m) / 2 / step
        R22 = (g2p - g2m) / 2 / step * g_true

        m, merr, c, cerr = _jack_est(g2, R22, g1, R11)
    else:
        R11 = (g1p - g1m) / 2 / step * g_true
        R22 = (g2p - g2m) / 2 / step

        m, merr, c, cerr = _jack_est(g1, R11, g2, R22)

    print("""\
# of sims: {n_sims}
wcs_g1   : {wcs_g1:f}
wcs_g2   : {wcs_g2:f}
dudx     : {dudx:f}
dudy     : {dudy:f}
dvdx     : {dvdx:f}
dvdy     : {dvdy:f}
m [1e-3] : {m:f} +/- {msd:f}
c [1e-4] : {c:f} +/- {csd:f}""".format(
        n_sims=len(g1),
        wcs_g1=wcs_g1,
        wcs_g2=wcs_g2,
        **jacobian_dict,
        m=m/1e-3,
        msd=merr/1e-3,
        c=c/1e-4,
        csd=cerr/1e-4), flush=True)


def _run_metacal(*, n_sims, rng, swap_g1g2, dudx, dudy, dvdx, dvdy):
    """Run metacal on an image composed of stamps w/ constant noise.

    Parameters
    ----------
    n_sims : int
        The number of objects to run.
    rng : np.random.RandomState
        An RNG to use.
    swap_g1g2 : bool
        If True, set the true shear on the 2-axis to 0.02 and 1-axis to 0.0.
        Otherwise, the true shear on the 1-axis is 0.02 and on the 2-axis is
        0.0.
    dudx : float
        The du/dx Jacobian component.
    dudy : float
        The du/dy Jacobian component.
    dydx : float
        The dv/dx Jacobian component.
    dvdy : float
        The dv/dy Jacobian component.

    Returns
    -------
    result : dict
        A dictionary with each of the metacal catalogs.
    """

    stamp_size = 33
    psf_stamp_size = 33

    cen = (stamp_size - 1) / 2
    psf_cen = (psf_stamp_size - 1)/2

    noise = 1
    flux = 64000

    galsim_jac = galsim.JacobianWCS(
        dudx=dudx,
        dudy=dudy,
        dvdx=dvdx,
        dvdy=dvdy)

    if swap_g1g2:
        g1 = 0.0
        g2 = 0.02
    else:
        g1 = 0.02
        g2 = 0.0

    gal = galsim.Exponential(
        half_light_radius=0.5
    ).withFlux(
        flux
    ).shear(
        g1=g1, g2=g2)

    psf = galsim.Gaussian(fwhm=0.9).withFlux(1)

    data = []
    for ind in tqdm.trange(n_sims):
        ################################
        # make the obs

        # psf
        psf_im = psf.drawImage(
            nx=psf_stamp_size,
            ny=psf_stamp_size,
            wcs=galsim_jac).array
        psf_noise = np.sqrt(np.sum(psf_im**2)) / 500
        wgt = np.ones_like(psf_im) / psf_noise**2
        psf_jac = ngmix.Jacobian(
            x=psf_cen,
            y=psf_cen,
            dudx=dudx,
            dudy=dudy,
            dvdx=dvdx,
            dvdy=dvdy)
        psf_obs = ngmix.Observation(
            image=psf_im,
            weight=wgt,
            jacobian=psf_jac)

        # now render object
        obj = galsim.Convolve(gal, psf)
        offset = rng.uniform(low=-0.5, high=0.5, size=2)
        im = obj.drawImage(
            nx=stamp_size,
            ny=stamp_size,
            wcs=galsim_jac,
            offset=offset).array
        jac = ngmix.Jacobian(
            x=cen+offset[0],
            y=cen+offset[1],
            dudx=dudx,
            dudy=dudy,
            dvdx=dvdx,
            dvdy=dvdy)
        wgt = np.ones_like(im) / noise**2
        nse = rng.normal(size=im.shape) * noise
        im += (rng.normal(size=im.shape) * noise)
        obs = ngmix.Observation(
            image=im,
            weight=wgt,
            noise=nse,
            bmask=np.zeros_like(im, dtype=np.int32),
            ormask=np.zeros_like(im, dtype=np.int32),
            jacobian=jac,
            psf=psf_obs
        )

        # build the mbobs
        mbobs = ngmix.MultiBandObsList()
        obslist = ngmix.ObsList()
        obslist.append(obs)
        mbobs.append(obslist)

        mbobs.meta['id'] = ind+1
        # these settings do not matter that much I think
        mbobs[0].meta['Tsky'] = 1
        mbobs[0].meta['magzp_ref'] = 26.5
        mbobs[0][0].meta['orig_col'] = ind+1
        mbobs[0][0].meta['orig_row'] = ind+1

        ################################
        # run the fitters
        try:
            res = _run_metacal_fitter(mbobs, rng)
        except Exception as e:
            print(e)
            res = None

        if res is not None:
            data.append(res)

    if len(data) > 0:
        res = data
    else:
        res = None

    return res


@numba.njit
def _jack_est(g1, R11, g2, R22):
    g1bar = np.mean(g1)
    R11bar = np.mean(R11)
    g2bar = np.mean(g2)
    R22bar = np.mean(R22)
    n = g1.shape[0]
    fac = n / (n-1)
    m_samps = np.zeros_like(g1)
    c_samps = np.zeros_like(g1)

    for i in range(n):
        _g1 = fac * (g1bar - g1[i]/n)
        _R11 = fac * (R11bar - R11[i]/n)
        _g2 = fac * (g2bar - g2[i]/n)
        _R22 = fac * (R22bar - R22[i]/n)
        m_samps[i] = _g1 / _R11 - 1
        c_samps[i] = _g2 / _R22

    m = np.mean(m_samps)
    c = np.mean(c_samps)

    m_err = np.sqrt(np.sum((m - m_samps)**2) / fac)
    c_err = np.sqrt(np.sum((c - c_samps)**2) / fac)

    return m, m_err, c, c_err


def _fit_psf(psf):
    runner = ngmix.bootstrap.PSFRunner(
        psf,
        'gauss',
        1.0,
        {'maxfev': 2000, 'ftol': 1.0e-5, 'xtol': 1.0e-5}
    )
    runner.go(ntry=2)

    psf_fitter = runner.fitter
    res = psf_fitter.get_result()
    psf.update_meta_data({'fitter': psf_fitter})

    if res['flags'] == 0:
        gmix = psf_fitter.get_gmix()
        psf.set_gmix(gmix)
    else:
        from ngmix.gexceptions import BootPSFFailure
        raise BootPSFFailure("failed to fit psfs: %s" % str(res))


def _run_metacal_fitter(mbobs, rng):
    # fit the PSF
    _fit_psf(mbobs[0][0].psf)

    metacal_pars = {
        'psf': 'fitgauss',
        'types': ['noshear', '1p', '1m', '2p', '2m'],
        'use_noise_image': True,
        'step': 0.01
    }
    moments_pars = {'bmask_flags': 2**30, 'weight': {'fwhm': 1.2}}

    obs_dict = metacal.get_all_metacal(mbobs, **metacal_pars)

    # overall flags, or'ed from each moments fit
    res = {'mcal_flags': 0}
    for key in sorted(obs_dict):
        try:
            fitter = Moments(moments_pars, rng)
            fres = fitter.go([obs_dict[key]])
        except Exception as err:
            print(err)
            fres = {'flags': np.ones(1, dtype=[('flags', 'i4')])}

        res['mcal_flags'] |= fres['flags'][0]
        tres = {}
        for name in fres.dtype.names:
            no_wmom = name.replace('wmom_', '')
            tres[no_wmom] = fres[name][0]
        tres['flags'] = fres['flags'][0]  # make sure this is moved over
        res[key] = tres

    return res


if __name__ == '__main__':
    if len(sys.argv) > 2:
        wcs_g1 = float(sys.argv[2])
    else:
        wcs_g1 = 0.0

    if len(sys.argv) > 3:
        wcs_g2 = float(sys.argv[3])
    else:
        wcs_g2 = wcs_g1

    run_metacal(n_sims=int(sys.argv[1]), wcs_g1=wcs_g1, wcs_g2=wcs_g2)

It only needs the Moments fitter 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 code to compute the net e1,e2 for a mixture model isn't quite right. I believe the code for get_e1e2T should replace the computation of irr,irc,icc with:

        irr=((gm['irr'] + gm['row']**2)*gm['p']).sum()*ipsum
        irc=((gm['irc'] + gm['row']*gm['col'])*gm['p']).sum()*ipsum
        icc=((gm['icc'] + gm['col']**2)*gm['p']).sum()*ipsum

I think that should make the code compute the correct moments and then the correct e1,e2,T.

Making this change can only increase the value of T, but I think the effect on e1,e2 can be in any direction. So it is possible that this explains why ngmix PSF shapes don't seem to match up with either im3shape or HSM.

Some time after commit 5efd6b3 a bug was introduced, or some behavior changed that is affecting programs that use ngmix.

@beckermr there are some commits from you in this range as well.

numba 0.41.0 has a bug that causes ngmix to fail. Def best to avoid it.

proposal that the user can use the functions like set_cen but can't modify the underlying data without using the api. This is to guarantee consistency.

E.g. this would raise an error

data = gm.get_data()
data['p'] = 5

with gmix.writeable():
    data = gm.get_data()
    data['p'] = 5

We can easily make it work for psfs and images with different pixel scales.

1. deconvolve image by pixel
2. deconvolve psf by its different pixel
3. deconvolve image by pixel-deconvolved psf
4. do metacal shears
5. reconvolve image by original pixel and chosen pre-pixel psf

I always see a lot of print statements like this:

     singular covariance
    pars at singular: 1.47e-08  -9.88e-09  1.84e-08  4.06e-11     0.274    0.0692 
     singular covariance
    pars at singular: 1.44e-08  -3.91e-10  -3.48e-08  7.51e-08     0.274    0.0692 

These should be logging statements at DEBUG.

Several joint prior classes in joint_prior.py require the centre prior to have the method cen_prior.get_lnprob_scalar_sep(pars[0],pars[1]). None of the priors defined in priors.py have this method. This leads to crashes when executing fill_fdiff.

add a context manger that switches on writeable refs and runs update_pixels on exit

this brings up again the issue of jacobians, which are always copies in the current code

Names like LM are confusing

• LM -> Fitter
• LMCoellip -> CoellipFitter
• GMixEM -> EMFitter
• GMixEMFixCen -> EMFitterFixCen
• GMixEMFluxOnly -> EmFitterFluxOnly
• TemplateFluxFitter -> PSFFluxFitter (we can still support non-psf flux fitting)

We should move all print statements to logging calls.

import ngmix
pars_sersic = [1.5, 0.0, 0.0, 0.2, -0.1, 16.0, 100.0]
ngmix.GMixModel(pars_sersic,"sersic")

leads to the above ngmix error. I checked that "sersic" is in the dictionary of implemented models, hence the number 8 for the model. But apparently some of the actual implementation is missing.

Is this the right way to implement a model with a free Sersic index? Also, I needed to guess the order of the parameters, but different permutations of the order did not get rid of the error.

this is needed by the balrog folks

For the ksigma tests, we need to look at the mean of the flux error in the test.

cc @esheldon who pointed this out

e.g. guesses for coellipe etc.

The last run of ngmix still had the sign error in e1. I believe the only change you need to make is in this line:

I think this should be (v,u), not (u,v). i.e. u->col, v->row. What you have now does u->row, v->col, which means e1 = (Icc-irr)/(Icc+Irr) will be (Ivv-Iuu)/(Ivv+Iuu), which is the source of the sign error.

The original implementation works such that the fluxes of each gaussian are normalized relative to the total flux in the image. This is inconvenient, update it to work in the correct units

I don't think this is really a breaking change because, due to the odd units, anyone using this code would eventually reset the flux anyway.

Hi.
I'm trying to do metacal with some real data. All I have are the observed ellipticities g1, g2.

In the example metacal.py, it seems that you build a simulation of an observed image using galsim with a shear_true(is this the induced shear?) on an exponential profile.

Am I supposed to build an image with the ellipticities? Is it possible for me to do this using galsim?

I'm sorry if this is a silly question, but I don't quite understand.

I redid the normalization of the dev variance values (f in ngmix) and got slightly different results. Here is my fixed version. I double checked it and it appears to be correct.

static const double PyGMix_fvals_dev[] = {
  2.9934935706271918e-07,
  3.4651596338231207e-06,
  2.4807910570562753e-05,
  1.4307404300535354e-04,
  7.2753169298239500e-04,
  3.4582464394427260e-03,
  1.6086645440719100e-02,
  7.7006776775654429e-02,
  4.1012562102501476e-01,
  2.9812509778548648e+00
};

The minimal set of things to cover is

• Jacobians
• gmixes
• fitting.py
• bootstrap.py (boostrappers and runners)
• metacal full reproducibility test (requires most of above)

e.g. old codes that we no longer use such as the importance samplers, etc.

This will facilitate leading order corrections for psf modeling errors

Develop test cases for which we can test the expected result at the pixel level

This is tricky and may not be possible at some level due to the inherently non-analytic operations involved.

One test we could include is something like the correlated noise test that is in descwl-shear-sims

Useful for example when storing whole images.

The pixel array is only useful for the optimizations that are used in the fitting codes

The following code works correctly on v1.0.0, but fails on v1.1.0rc1.

jac = ngmix.Jacobian(dudx=0, dudy=-0.2, dvdx=-0.2, dvdy=0, x=23, y=23)
assert np.allclose(jac.get_sdet(), 0.2)

The problem is a missing abs in the calculation of sdet in Jacobian._finish_init. The bug was introduced when you changed the definition of det to no longer mean the absolute value of the determinant.

Indeed I'm not completely sure if the latter change (which seems to be an undocumented API change in v1.0 -> v1.1) was intended or not. But I'm quite sure that sdet should mean the sqrt(|det|) regardless.

(This bug bit Rebecca Chen recently when using ngmix to measure Piff shapes.)

These were badly designed from the beginning and are mostly not used any more.

I think this can be accomplished by putting in the correct pixel area factors when processing each image.

currently we are instead modifying the data so it is surface brightness (e.g. in fitvd), but might as well do this in the code.

if AM fails, try something else

This code in bootstrap.py near line 1540 (which is meant to check for a PSF fit) fails if band 0 has no observations.

    def fit_gal_psf_flux(self, normalize_psf=True):
        """
        use psf as a template, measure flux (linear)
        """


        mbo = self.mb_obs_list
        nband = len(mbo)

        if not mbo[0][0].psf.has_gmix():
            raise RuntimeError("you need to fit the psfs first")

The bootstrapper can generate empty observation lists if the PSF fits fail and no observations are kept.

I am going to port over the moments fitter from metadetect.

I think there may be a typo in the docstring for GMixBDF. It states that the constructor parameters should be the following:

Parameters
----------
pars: sequence
        [c1,c2,g1,g2,T,flux]
TdByTe: number, optional
        Optionally set TdByTe.  Defaults to 1.0

However when I tried constructing a GMixBDF object this way, it fails:

gmix.GMixBDF(pars=[0, 0, 0.25, 0.25, 1., 100.], TdByTe=2.)

---------------------------------------------------------------------------
GMixFatalError                            Traceback (most recent call last)
<ipython-input-4-fcf2c55ef9ea> in <module>()
----> 1 gmix.GMixBDF(pars=[0, 0, 0.25, 0.25, 1., 100.], TdByTe=2.)

/home/spencer/miniconda3/envs/balrog/lib/python2.7/site-packages/ngmix/gmix.pyc in __init__(self, pars, TdByTe)
    955 
    956         self._TdByTe = float(TdByTe)
--> 957         super(GMixBDF,self)._do_init(pars,'bdf')
    958 
    959     def copy(self):

/home/spencer/miniconda3/envs/balrog/lib/python2.7/site-packages/ngmix/gmix.pyc in _do_init(self, pars, model)
    848 
    849         self._set_fill_func()
--> 850         self.fill(pars)
    851 
    852     def copy(self):

/home/spencer/miniconda3/envs/balrog/lib/python2.7/site-packages/ngmix/gmix.pyc in fill(self, pars)
    309             err="model '%s' requires %s pars, got %s"
    310             err =err % (self._model_name,self._npars, npars)
--> 311             raise GMixFatalError(err)
    312 
    313         self._fill(pars)

GMixFatalError: "model 'bdf' requires 7 pars, got 6"

Presumably this is because fracdev was never assigned at construction time. Checking gmix_fill_bdf it looks like the docstring should instead read:

Parameters
----------
pars: sequence
        [c1,c2,g1,g2,T,fracdev,flux]
TdByTe: number, optional
        Optionally set TdByTe.  Defaults to 1.0

End py2 support and make a new major release

Will the older 1.* releases that support py27 in conda-forge still work?

TODO

• rework most of the bootstrappers, without backwards compatibility
• remove unused code that probably doesn't even work, e.g. importance samplers, etc.
• improve test coverage
• add pixel level tests for metacal operations

Going over the different steps in the metacal procedure, I got a bit puzzled by this line where if I understand correctly the output images are drawn:

What puzzles me is when drawn here, the shear galaxy image should have already been reconvolved with the target PSF, which already contains the pixel response (added here), so I would expect a method="no_pixel" when drawing this galaxy.

To make things clearer for me, I tried to do a pure GalSim equivalent, starting from the metacal example script here. Given some obs and obsdict obtained by metacal, to reproduce the output with GalSim I need to do:

# Grab observations
im   = galsim.InterpolatedImage(galsim.Image(obs.image, scale=0.263))
psf  = galsim.InterpolatedImage(galsim.Image(obs.psf.image, scale=0.263))
rpsf = obsdict['noshear'].psf.galsim_obj 
# Deconvolve input image by PSF+pixel
gal = galsim.Convolve([im, galsim.Deconvolve(psf)])
# Apply shear and reconvolve by target PSF !!!which should contain the pixel response as well !!!
gal_1p = galsim.Convolve([gal.shear(g1=0.01,g2=0), rpsf])
# Draw the postage stamp
rec_1p = gal_1p.drawImage(nx=42,ny=42, scale=0.263).array

which gives me the following residuals when compared to obsdict['1p'].image:

but note that I draw the galaxy without the method="no_pixel" flag despite the reconvolution PSF being already pixel convolved.

Am I missing something? or is the output metacal image convolved with the pixel response twice (which doesn't matter that much probably)?

So right now we have a lot of different APIs for fitting/bootstrapping. I propose we unify them around a single API.

Each fitter will be a class with the following init

class MyFitter(BaseFitter):
    def __init__(self, fitter_par1, fitter_par2, ..., rng=None):
        # do init here

The method for fitting will take a list of mbobs and return a list of fitter results

    def go(self, list_of_mbobs, guesses=None):
        # do fit

We will accept single obs or obslists or single mbobs and do an appropriate conversion. On return, a scalar input (a single obs, single obslist or single mbobs) will result in a scalar return. The return value will be a list (or scalar) of the fit results. These fit results should be a dictionary with any appropriate things. We always require a 'flags' entry with non-zero values indicating errors.

The go method will have an optional keyword argument that allows the user to pass their own guesses for the fit. If None then the fitter should make a reasonable guess.

No one would ever want the current default in a real world situation