From acf7cda3404f3cd6b8e22c77079076e8ee1f6c36 Mon Sep 17 00:00:00 2001
From: Lucie Baumont <lucie.baumont@cea.fr>
Date: Thu, 12 Jan 2023 16:51:33 +0100
Subject: [PATCH 01/11] no longer setting weights to 1 before rescaling

---
 shapepipe/modules/ngmix_package/ngmix.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/shapepipe/modules/ngmix_package/ngmix.py b/shapepipe/modules/ngmix_package/ngmix.py
index 7cb614a42..d54b1c10a 100644
--- a/shapepipe/modules/ngmix_package/ngmix.py
+++ b/shapepipe/modules/ngmix_package/ngmix.py
@@ -813,7 +813,6 @@ def do_ngmix_metacal(
 
         weight_map = np.copy(weights[n_e])
         weight_map[np.where(flags[n_e] != 0)] = 0.
-        weight_map[weight_map != 0] = 1
 
         psf_guess = np.array([0., 0., 0., 0., psf_T, 1.])
         try:

From d5a2f07bf1b301bbd578e0c8008f5dcce61d57ff Mon Sep 17 00:00:00 2001
From: Lucie Baumont <lucie.baumont@cea.fr>
Date: Thu, 26 Jan 2023 15:32:02 +0100
Subject: [PATCH 02/11] tiny documentation change

---
 shapepipe/modules/find_exposures_package/__init__.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/shapepipe/modules/find_exposures_package/__init__.py b/shapepipe/modules/find_exposures_package/__init__.py
index 0b0debcc7..579107c26 100644
--- a/shapepipe/modules/find_exposures_package/__init__.py
+++ b/shapepipe/modules/find_exposures_package/__init__.py
@@ -14,8 +14,8 @@
 ===========
 
 Identify the exposure images that were co-added to produce the tiles
-(stacked image). The image names are listed in the tile FITS header,
-which is read by this module to extract the names.
+(stacked image). The image names are listed in the megapipe tile FITS 
+header,which is read by this module to extract the names.
 
 The output ASCII file contains the image base names (without file extension).
 

From 94552f8e59d18a9f1587bb8d769613a3a844b4b3 Mon Sep 17 00:00:00 2001
From: Lucie Baumont <lucie.baumont@cea.fr>
Date: Thu, 2 Feb 2023 15:56:29 +0100
Subject: [PATCH 03/11] force tmp directory creation

---
 shapepipe/modules/ngmix_package/#ngmix.py# | 1029 ++++++++++++++++++++
 shapepipe/pipeline/file_handler.py         |    1 +
 2 files changed, 1030 insertions(+)
 create mode 100644 shapepipe/modules/ngmix_package/#ngmix.py#

diff --git a/shapepipe/modules/ngmix_package/#ngmix.py# b/shapepipe/modules/ngmix_package/#ngmix.py#
new file mode 100644
index 000000000..fa7d15a0f
--- /dev/null
+++ b/shapepipe/modules/ngmix_package/#ngmix.py#
@@ -0,0 +1,1029 @@
+"""NGMIX.
+
+This module contains a class for ngmix shape measurement.
+
+:Author: Axel Guinot
+
+"""
+
+import re
+
+import galsim
+import ngmix
+import numpy as np
+from astropy.io import fits
+from modopt.math.stats import sigma_mad
+from ngmix.observation import MultiBandObsList, Observation, ObsList
+from numpy.random import uniform as urand
+from sqlitedict import SqliteDict
+
+from shapepipe.pipeline import file_io
+
+
+class Ngmix(object):
+    """Ngmix.
+
+    Class to handle NGMIX shapepe measurement.
+
+    Parameters
+    ----------
+    input_file_list : list
+        Input files
+    output_dir : str
+        Output directory
+    file_number_string : str
+        File numbering scheme
+    zero_point : float
+        Photometric zero point
+    pixel_scale : float
+        Pixel scale in arcsec
+    f_wcs_path : str
+        Path to merged single-exposure single-HDU headers
+    w_log : logging.Logger
+        Logging instance
+    id_obj_min : int, optional
+        First galaxy ID to process, not used if the value is set to ``-1``;
+        the default is ``-1``
+    id_obj_max : int, optional
+        Last galaxy ID to process, not used if the value is set to ``-1``;
+        the default is ``-1``
+
+    Raises
+    ------
+    IndexError
+        If the length of the input file list is incorrect
+
+    """
+
+    def __init__(
+        self,
+        input_file_list,
+        output_dir,
+        file_number_string,
+        zero_point,
+        pixel_scale,
+        f_wcs_path,
+        w_log,
+        id_obj_min=-1,
+        id_obj_max=-1
+    ):
+
+        if len(input_file_list) != 6:
+            raise IndexError(
+                f'Input file list has length {len(input_file_list)},'
+                + ' required is 6'
+            )
+
+        self._tile_cat_path = input_file_list[0]
+        self._gal_vignet_path = input_file_list[1]
+        self._bkg_vignet_path = input_file_list[2]
+        self._psf_vignet_path = input_file_list[3]
+        self._weight_vignet_path = input_file_list[4]
+        self._flag_vignet_path = input_file_list[5]
+
+        self._output_dir = output_dir
+        self._file_number_string = file_number_string
+
+        self._zero_point = zero_point
+        self._pixel_scale = pixel_scale
+
+        self._f_wcs_path = f_wcs_path
+        self._id_obj_min = id_obj_min
+        self._id_obj_max = id_obj_max
+
+        self._w_log = w_log
+
+        # Initiatlise random generator
+        seed = int(''.join(re.findall(r'\d+', self._file_number_string)))
+        np.random.seed(seed)
+        self._w_log.info(f'Random generator initialisation seed = {seed}')
+
+    @classmethod
+    def MegaCamFlip(self, vign, ccd_nb):
+        """Flip for MegaCam.
+
+        MegaPipe has CCDs that are upside down. This function flips the
+        postage stamps in these CCDs. TO DO: This will give incorrect results
+        when used with THELI ccds.  Fix this.
+
+        Parameters
+        ----------
+        vign : numpy.ndarray
+            Array containing the postage stamp to flip
+        ccd_nb : int
+            ID of the CCD containing the postage stamp
+
+        Returns
+        -------
+        numpy.ndarray
+            The flipped postage stamp
+
+        """
+        if ccd_nb < 18 or ccd_nb in [36, 37]:
+            # swap x axis so origin is on top-right
+            return np.rot90(vign, k=2)
+            #print('rotating megapipe image')
+        else:
+            # swap y axis so origin is on bottom-left
+            return vign
+
+    def get_prior(self):
+        """Get Prior.
+
+        Return prior for the different parameters.
+        Prior on ellipticity is from Bernstein and Armstrong 2014.
+        Prior on centers is a 2d gaussian.
+        Priors on size and flux are flat.
+         
+
+        Returns
+        -------
+        ngmix.priors
+            Priors for the different parameters (ellipticity,center, size, flux)
+        """
+        # Prior on ellipticity. Details do not matter, as long
+        # as it regularizes the fit. From Bernstein & Armstrong 2014
+        g_sigma = 0.4
+        rng = np.random.RandomState(22)
+        g_prior = ngmix.priors.GPriorBA(g_sigma,rng)
+
+        # 2-d Gaussian prior on the center row and column center
+        # (relative to the center of the jacobian, which
+        # would be zero) and the sigma of the Gaussians.
+        # Units same as jacobian, probably arcsec
+        row, col = 0.0, 0.0
+        row_sigma, col_sigma = self._pixel_scale, self._pixel_scale
+        cen_prior = ngmix.priors.CenPrior(row, col, row_sigma, col_sigma, rng)
+
+        # Size prior. Instead of flat, two-sided error function (TwoSidedErf)
+        # could be used, though they are not used here
+        r50minval = -10.0  # arcsec squared
+        r50maxval = 1.e6
+        r50_prior = ngmix.priors.FlatPrior(r50minval, r50maxval, rng)
+
+        # Flux prior. Bounds need to make sense for
+        # images in question
+        Fminval = -1.e4
+        Fmaxval = 1.e9
+        F_prior = ngmix.priors.FlatPrior(Fminval, Fmaxval, rng)
+
+        # Joint prior, combine all individual priors
+        prior = ngmix.joint_prior.PriorSimpleSep(
+            cen_prior,
+            g_prior,
+            r50_prior,
+            F_prior
+        )
+
+        return prior
+
+    def compile_results(self, results):
+        """Compile Results.
+
+        Prepare the results of NGMIX before saving. TODO: add snr_r and T_r
+
+        Parameters
+        ----------
+        results : dict
+            Results of NGMIX metacal
+
+        Returns
+        -------
+        dict
+            Compiled results ready to be written to a file
+            note: psfo is the original image psf from psfex or mccd
+
+        Raises
+        ------
+        KeyError
+            If SNR key not found
+
+        """
+        names = ['1m', '1p', '2m', '2p', 'noshear']
+        names2 = [
+            'id',
+            'n_epoch_model',
+            'moments_fail',
+            'ntry_fit',
+            'g1_psfo_ngmix',
+            'g2_psfo_ngmix',
+            'r50_psfo_ngmix',
+            'g1_err_psfo_ngmix',
+            'g2_err_psfo_ngmix',
+            'r50_err_psfo_ngmix',
+            'g1',
+            'g1_err',
+            'g2',
+            'g2_err',
+            'r50',
+            'r50_err',
+            'r50psf',
+            'g1_psf',
+            'g2_psf',
+            'flux',
+            'flux_err',
+            's2n',
+            'mag',
+            'mag_err',
+            'flags',
+            'mcal_flags'
+        ]
+        output_dict = {k: {kk: [] for kk in names2} for k in names}
+        for idx in range(len(results)):
+            for name in names:
+
+                mag = (
+                    -2.5 * np.log10(results[idx][name]['flux'])
+                    + self._zero_point
+                )
+                mag_err = np.abs(
+                    -2.5 * results[idx][name]['flux_err']
+                    / (results[idx][name]['flux'] * np.log(10))
+                )
+
+                output_dict[name]['id'].append(results[idx]['obj_id'])
+                output_dict[name]['n_epoch_model'].append(
+                    results[idx]['n_epoch_model']
+                )
+                output_dict[name]['moments_fail'].append(
+                    results[idx]['moments_fail']
+                )
+                output_dict[name]['ntry_fit'].append(
+                    results[idx][name]['ntry']
+                )
+                output_dict[name]['g1_psfo_ngmix'].append(
+                    results[idx]['g_PSFo'][0]
+                )
+                output_dict[name]['g2_psfo_ngmix'].append(
+                    results[idx]['g_PSFo'][1]
+                )
+                output_dict[name]['g1_err_psfo_ngmix'].append(
+                    results[idx]['g_err_PSFo'][0]
+                )
+                output_dict[name]['g2_err_psfo_ngmix'].append(
+                    results[idx]['g_err_PSFo'][1]
+                )
+                output_dict[name]['r50_psfo_ngmix'].append(
+                    results[idx]['r50_PSFo']
+                )
+                output_dict[name]['r50_err_psfo_ngmix'].append(
+                    results[idx]['r50_err_PSFo']
+                )
+                output_dict[name]['g1'].append(results[idx][name]['g'][0])
+                output_dict[name]['g1_err'].append(
+                    results[idx][name]['pars_err'][2]
+                )
+                output_dict[name]['g2'].append(results[idx][name]['g'][1])
+                output_dict[name]['g2_err'].append(
+                    results[idx][name]['pars_err'][3]
+                )
+                output_dict[name]['r50'].append(results[idx][name]['pars'][4])
+                output_dict[name]['r50_err'].append(results[idx][name]['pars_err'][4])
+                output_dict[name]['r50psf'].append(results[idx][name]['r50psf'])
+                output_dict[name]['g1_psf'].append(
+                    results[idx][name]['gpsf'][0]
+                )
+                output_dict[name]['g2_psf'].append(
+                    results[idx][name]['gpsf'][1]
+                )
+                output_dict[name]['flux'].append(results[idx][name]['flux'])
+                output_dict[name]['flux_err'].append(
+                    results[idx][name]['flux_err']
+                )
+                output_dict[name]['mag'].append(mag)
+                output_dict[name]['mag_err'].append(mag_err)
+
+                if 's2n' in results[idx][name]:
+                    output_dict[name]['s2n'].append(results[idx][name]['s2n'])
+                elif 's2n_r' in results[idx][name]:
+                    output_dict[name]['s2n'].append(
+                        results[idx][name]['s2n_r']
+                    )
+                else:
+                    raise KeyError('No SNR key (s2n, s2n_r) found in results')
+
+                output_dict[name]['flags'].append(results[idx][name]['flags'])
+                output_dict[name]['mcal_flags'].append(
+                    results[idx]['mcal_flags']
+                )
+
+        return output_dict
+
+    def save_results(self, output_dict):
+        """Save Results.
+
+        Save the results into a FITS file.
+
+        Parameters
+        ----------
+        output_dict
+            Dictionary containing the results
+
+        """
+        output_name = (
+            f'{self._output_dir}/ngmix{self._file_number_string}.fits'
+        )
+
+        f = file_io.FITSCatalogue(
+            output_name,
+            open_mode=file_io.BaseCatalogue.OpenMode.ReadWrite
+        )
+
+        for key in output_dict.keys():
+            f.save_as_fits(output_dict[key], ext_name=key.upper())
+
+    def process(self):
+        """Process.
+
+        Funcion to processs NGMIX.
+        organizes object cutouts from detection catalog in image, 
+        weight, and flag files
+        per object: 
+            gathers wcs and psf info from exposures
+            background subtracts (make this an option)
+            scales by relative zeropoints
+            runs metacal convolutions and ngmix fitting
+        Returns
+        -------
+        dict
+            Dictionary containing the NGMIX metacal results
+
+        """
+        # sextractor detection catalog for the tile
+        tile_cat = file_io.FITSCatalogue(
+            self._tile_cat_path,
+            SEx_catalogue=True,
+        )
+        tile_cat.open()
+        obj_id = np.copy(tile_cat.get_data()['NUMBER'])
+        tile_vign = np.copy(tile_cat.get_data()['VIGNET'])
+        tile_ra = np.copy(tile_cat.get_data()['XWIN_WORLD'])
+        tile_dec = np.copy(tile_cat.get_data()['YWIN_WORLD'])
+        tile_cat.close()
+
+        f_wcs_file = SqliteDict(self._f_wcs_path)
+        gal_vign_cat = SqliteDict(self._gal_vignet_path)
+        bkg_vign_cat = SqliteDict(self._bkg_vignet_path)
+        psf_vign_cat = SqliteDict(self._psf_vignet_path)
+        weight_vign_cat = SqliteDict(self._weight_vignet_path)
+        flag_vign_cat = SqliteDict(self._flag_vignet_path)
+
+        final_res = []
+        prior = self.get_prior()
+
+        count = 0
+        id_first = -1
+        id_last = -1
+
+        for i_tile, id_tmp in enumerate(obj_id):
+            if self._id_obj_min > 0 and id_tmp < self._id_obj_min:
+                continue
+            if self._id_obj_max > 0 and id_tmp > self._id_obj_max:
+                continue
+
+            if id_first == -1:
+                id_first = id_tmp
+            id_last = id_tmp
+
+            count = count + 1
+
+            gal_vign = []
+            psf_vign = []
+            sigma_psf = []
+            weight_vign = []
+            flag_vign = []
+            jacob_list = []
+            if (
+                (psf_vign_cat[str(id_tmp)] == 'empty')
+                or (gal_vign_cat[str(id_tmp)] == 'empty')
+            ):
+                continue
+            #identify exposure and ccd number from psf catalog
+            psf_expccd_name = list(psf_vign_cat[str(id_tmp)].keys())
+            for expccd_name_tmp in psf_expccd_name:
+                exp_name, ccd_n = re.split('-', expccd_name_tmp)
+
+                gal_vign_tmp = (
+                    gal_vign_cat[str(id_tmp)][expccd_name_tmp]['VIGNET']
+                )
+                if len(np.where(gal_vign_tmp.ravel() == 0)[0]) != 0:
+                    continue
+                # background subtraction
+                bkg_vign_tmp = (
+                    bkg_vign_cat[str(id_tmp)][expccd_name_tmp]['VIGNET']
+                )
+                gal_vign_sub_bkg = gal_vign_tmp - bkg_vign_tmp
+                # skip this step for THELI
+                tile_vign_tmp = (
+                    Ngmix.MegaCamFlip(np.copy(tile_vign[i_tile]), int(ccd_n))
+                )
+
+                flag_vign_tmp = (
+                    flag_vign_cat[str(id_tmp)][expccd_name_tmp]['VIGNET']
+                )
+                flag_vign_tmp[np.where(tile_vign_tmp == -1e30)] = 2**10
+                v_flag_tmp = flag_vign_tmp.ravel()
+                if len(np.where(v_flag_tmp != 0)[0]) / (51 * 51) > 1 / 3.0:
+                    continue
+
+                weight_vign_tmp = (
+                    weight_vign_cat[str(id_tmp)][expccd_name_tmp]['VIGNET']
+                )
+
+                jacob_tmp = get_jacob(
+                    f_wcs_file[exp_name][int(ccd_n)]['WCS'],
+                    tile_ra[i_tile],
+                    tile_dec[i_tile]
+                )
+
+                header_tmp = fits.Header.fromstring(
+                    f_wcs_file[exp_name][int(ccd_n)]['header']
+                )
+                #rescale images and weights by relative flux scale
+                # this should be it's own function
+                Fscale = header_tmp['FSCALE']
+
+                gal_vign_scaled = gal_vign_sub_bkg * Fscale
+                weight_vign_scaled = weight_vign_tmp * 1 / Fscale ** 2
+
+                gal_vign.append(gal_vign_scaled)
+                psf_vign.append(
+                    psf_vign_cat[str(id_tmp)][expccd_name_tmp]['VIGNET']
+                )
+                sigma_psf.append(
+                    psf_vign_cat[
+                        str(id_tmp)
+                    ][expccd_name_tmp]['SHAPES']['SIGMA_PSF_HSM']
+                )
+                weight_vign.append(weight_vign_scaled)
+                flag_vign.append(flag_vign_tmp)
+                jacob_list.append(jacob_tmp)
+                
+            #if object is observed, carry out metacal operations and run ngmix
+            if len(gal_vign) == 0:
+                continue
+            try:
+                res = do_ngmix_metacal(
+                    gal_vign,
+                    psf_vign,
+                    sigma_psf,
+                    weight_vign,
+                    flag_vign,
+                    jacob_list,
+                    prior,
+                    self._pixel_scale
+                )
+
+            except Exception as ee:
+                self._w_log.info(
+                    f'ngmix failed for object ID={id_tmp}.\nMessage: {ee}'
+                )
+                continue
+
+            res['obj_id'] = id_tmp
+            res['n_epoch_model'] = len(gal_vign)
+            final_res.append(res)
+
+        self._w_log.info(
+            f'ngmix loop over objects finished, processed {count} '
+            + f'objects, id first/last={id_first}/{id_last}'
+        )
+
+        f_wcs_file.close()
+        gal_vign_cat.close()
+        bkg_vign_cat.close()
+        flag_vign_cat.close()
+        weight_vign_cat.close()
+        psf_vign_cat.close()
+
+        # Put all results together
+        res_dict = self.compile_results(final_res)
+
+        # Save results
+        self.save_results(res_dict)
+
+
+def get_guess_shapeHSM(
+    img,
+    pixel_scale,
+    guess_flux_unit='img',
+    guess_size_type='T',
+    guess_size_unit='sky',
+    guess_centroid=True,
+    guess_centroid_unit='sky'
+):
+    r"""Get Guess using galsim hsm shapes.
+
+    Get the guess vector for the NGMIX shape measurement
+    ``[center_x, center_y, g1, g2, size_T, flux]``.
+    No guesses are given for the ellipticity ``(0, 0)``.
+
+    Parameters
+    ----------
+    img : numpy.ndarray
+        Array containing the image
+    pixel_scale : float
+        Approximation of the pixel scale
+    guess_flux_unit : str
+        If ``img`` returns the flux in pixel units, otherwise if ``sky``
+        returns the flux in :math:`{\rm arcsec}^{-2}`
+    guess_size_type : str
+        If ``T`` returns the size in quadrupole moments definition
+        :math:`2\sigma^2`, if ``sigma`` returns the moments
+        :math:`\sigma`, if ``r50`` returns the half light radius
+    guess_size_unit : str
+        If ``img`` returns the size in pixel units, otherwise if ``sky``
+        returns the size in arcsec
+    guess_centroid : bool
+        If ``True``, will return a guess on the object centroid, otherwise if
+        ``False``, will return the image centre
+    guess_centroid_unit : str
+        If ``img`` returns the centroid in pixel unit, otherwise if ``sky``
+        returns the centroid in arcsec
+
+    Returns
+    -------
+    numpy.ndarray
+        Return the guess array ``[center_x, center_y, g1, g2, size_T, flux]``
+
+    Raises
+    ------
+    GalSimHSMError
+        For an error in the computation of adaptive moments
+    ValueError
+        For invalid unit guess types
+
+    """
+    galsim_img = galsim.Image(img, scale=pixel_scale)
+
+    hsm_shape = galsim.hsm.FindAdaptiveMom(galsim_img, strict=False)
+
+    error_msg = hsm_shape.error_message
+                
+    if error_msg != '':
+        raise galsim.hsm.GalSimHSMError(
+            f'Error in adaptive moments :\n{error_msg}'
+        )
+
+    if guess_flux_unit == 'img':
+        guess_flux = hsm_shape.moments_amp
+    elif guess_flux_unit == 'sky':
+        guess_flux = hsm_shape.moments_amp / pixel_scale ** 2
+    else:
+        raise ValueError(
+            f'invalid guess_flux_unit \'{guess_flux_unit}\','
+            + ' must be one of \'img\', \'sky\''
+        )
+
+    if guess_size_unit == 'img':
+        size_unit = 1.
+    elif guess_size_unit == 'sky':
+        size_unit = pixel_scale
+    else:
+        raise ValueError(
+            'invalid guess_size_unit \'{guess_size_unit}\','
+            + 'must be one of \'img\', \'sky\''
+        )
+
+    if guess_size_type == 'sigma':
+        guess_size = hsm_shape.moments_sigma * size_unit
+    elif guess_size_type == 'T':
+        guess_size = 2 * (hsm_shape.moments_sigma * size_unit) ** 2
+    elif guess_size_type == 'r50':
+        guess_size = hsm_shape.moments_sigma * size_unit * 1.17741002252
+    else:
+        raise ValueError(
+            'invalid guess_size_type \'{guess_size_type}\','
+            + 'must be one of \'sigma\', \'T\', or \'r50\''
+        )
+
+    if guess_centroid_unit == 'img':
+        centroid_unit = 1
+    elif guess_centroid_unit == 'sky':
+        centroid_unit = pixel_scale
+    else:
+        raise ValueError(
+            f'invalid guess_centroid_unit \'{guess_centroid_unit}\','
+            + '  must be one of \'img\', \'sky\''
+        )
+
+    if guess_centroid:
+        guess_centroid = (
+            (hsm_shape.moments_centroid - galsim_img.center) * centroid_unit
+        )
+    else:
+        guess_centroid = galsim_img.center * centroid_unit
+
+    guess = np.array([
+        guess_centroid.x,
+        guess_centroid.y,
+        0.,
+        0.,
+        guess_size,
+        guess_flux
+    ])
+
+    return guess
+
+
+def make_galsimfit(obs, model, guess0, prior=None, ntry=5):
+    """Make GalSim Fit.
+
+    wrapper for ngmix image fit using simple GalSim model.
+
+    Parameters
+    ----------
+    obs : ngmix.observation.Observation
+        Image to fit
+    model : str
+        Model for fit
+    guess0 : numpy.ndarray
+        Parameters of first model guess
+    prior : ngmix.prior, optional
+        Prior for fit parameters
+    ntry : int, optional
+        Number of tries for fit, the default is ``5``
+
+    Returns
+    -------
+    dict
+        Results
+
+    Raises
+    ------
+    ngmix.BootGalFailure
+        Failure to bootstrap galaxy
+
+    """
+
+    limit = 0.1
+
+    guess = np.copy(guess0)
+    fres = {}
+    for it in range(ntry):
+        guess[0:5] += urand(low=-limit, high=limit)
+        guess[5:] *= (1 + urand(low=-limit, high=limit))
+        fres['flags'] = 1
+        #OLD_LM_PARS = {"maxfev": 1000, "ftol": 1.0e6, "xtol": 1.0e-6}
+        try:
+            fitter = ngmix.fitting.GalsimFitter(
+                model=model,
+                prior=prior
+            )
+            fres = fitter.go(obs=obs,guess=guess)
+            #print(fres.items())
+        except Exception:
+            continue
+
+        if fres['flags'] == 0:
+            break
+
+    if fres['flags'] != 0:
+        raise ngmix.gexceptions.BootGalFailure(
+            'Failed to fit galaxy image with galsimfit'
+        )
+
+    fres['ntry'] = it + 1
+    return fres
+
+
+def get_jacob(wcs, ra, dec):
+    """Get Jacobian.
+
+    Return the Jacobian of the WCS at the required position.
+
+    Parameters
+    ----------
+    wcs : astropy.wcs.WCS
+        WCS object for which we want the Jacobian
+    ra : float
+        RA position of the center of the vignet (in degrees)
+    dec : float
+        Dec position of the center of the vignet (in degress)
+
+    Returns
+    -------
+    galsim.wcs.BaseWCS.jacobian
+        Jacobian of the WCS at the required position
+
+    """
+    g_wcs = galsim.fitswcs.AstropyWCS(wcs=wcs)
+    world_pos = galsim.CelestialCoord(
+        ra=ra * galsim.angle.degrees,
+        dec=dec * galsim.angle.degrees,
+    )
+    galsim_jacob = g_wcs.jacobian(world_pos=world_pos)
+
+    return galsim_jacob
+
+
+def get_noise(gal, weight, guess, pixel_scale, thresh=1.2):
+    r"""Get Noise.
+
+    Compute the sigma of the noise from an object postage stamp.
+    Use a guess on the object size, ellipticity and flux to create a window
+    function.
+
+    Parameters
+    ----------
+    gal : numpy.ndarray
+        Galaxy image
+    weight : numpy.ndarray
+        Weight image
+    guess : list
+        Gaussian parameters fot the window function
+        ``[x0, y0, g1, g2, r50, flux]``
+    pixel_scale : float
+        Pixel scale of the galaxy image
+    thresh : float, optional
+        Threshold to cut the window function,
+        cut = ``thresh`` * :math:`\sigma_{\rm noise}`;  the default is ``1.2``
+
+    Returns
+    -------
+    float
+        Sigma of the noise on the galaxy image
+
+    """
+    img_shape = gal.shape
+    m_weight = weight != 0
+
+    sig_tmp = sigma_mad(gal[m_weight])
+
+    gauss_win = galsim.Gaussian(sigma=np.sqrt(guess[4] / 2), flux=guess[5])
+    gauss_win = gauss_win.shear(g1=guess[2], g2=guess[3])
+    gauss_win = gauss_win.drawImage(
+        nx=img_shape[0],
+        ny=img_shape[1],
+        scale=pixel_scale
+    ).array
+
+    m_weight = weight[gauss_win < thresh * sig_tmp] != 0
+
+    sig_noise = sigma_mad(gal[gauss_win < thresh * sig_tmp][m_weight])
+
+    return sig_noise
+
+
+def do_ngmix_metacal(
+    gals,
+    psfs,
+    psfs_sigma,
+    weights,
+    flags,
+    jacob_list,
+    prior,
+    pixel_scale
+):
+    """Do Ngmix Metacal.
+
+    Performs  metacalibration on a multi-epoch object and return the joint
+    shape measurement with NGMIX.
+
+    Parameters
+    ----------
+    gals : list
+        List of the galaxy vignets.  List indices run over epochs
+    psfs : list
+        List of the PSF vignets
+    psfs_sigma : list
+        List of the sigma PSFs
+    weights : list
+        List of the weight vignets
+    flags : list
+        List of the flag vignets
+    jacob_list : list
+        List of the Jacobians
+    prior : ngmix.priors
+        Priors for the fitting parameters
+    pixel_scale : float
+        pixel scale in arcsec
+
+    Returns
+    -------
+    dict
+        Dictionary containing the results of NGMIX metacal
+
+    """
+    n_epoch = len(gals)
+
+    if n_epoch == 0:
+        raise ValueError("0 epoch to process")
+
+    # Construct observation objects to pass to ngmix 
+    gal_obs_list = ObsList()
+    r50_guess_psf = []
+    psf_res_gT = {
+        'g_PSFo': np.array([0., 0.]),
+        'g_err_PSFo': np.array([0., 0.]),
+        'r50_PSFo': 0.,
+        'r50_err_PSFo': 0.
+    }
+    gal_guess = []
+    gal_guess_flag = True
+    wsum = 0
+    for n_e in range(n_epoch):
+
+        psf_jacob = ngmix.Jacobian(
+            row=(psfs[0].shape[0] - 1) / 2,
+            col=(psfs[0].shape[1] - 1) / 2,
+            wcs=jacob_list[n_e]
+        )
+        # psf observation is part of ngmix observation
+        psf_obs = Observation(psfs[n_e], jacobian=psf_jacob)
+        # convert sigma to T, which is the half-light radius
+        # NOT the usual definition T=2*sigma^2
+        psf_r50 = psfs_sigma[n_e] * 1.17741  * pixel_scale
+
+        # integrate flag info into weights
+        weight_map = np.copy(weights[n_e])
+        weight_map[np.where(flags[n_e] != 0)] = 0.
+        weight_map[weight_map != 0] = 1
+
+        # fit gaussian to psf
+        psf_guess = np.array([0., 0., 0., 0., psf_r50, 1.])
+        try:
+            psf_res = make_galsimfit(psf_obs, 'gauss', psf_guess)
+        except Exception:
+            continue
+
+        # Gal guess
+        try:
+            gal_guess_tmp = get_guess_shapeHSM(
+                gals[n_e],
+                pixel_scale,
+                guess_size_type='sigma'
+            )
+        except Exception:
+            gal_guess_flag = False
+            gal_guess_tmp = np.array([0., 0., 0., 0., 1, 100])
+
+        # Recenter jacobian if necessary
+        gal_jacob = ngmix.Jacobian(
+            row=(gals[0].shape[0] - 1) / 2 + gal_guess_tmp[0],
+            col=(gals[0].shape[1] - 1) / 2 + gal_guess_tmp[1],
+            wcs=jacob_list[n_e]
+        )
+
+        # Noise handling: this should be it's own function
+        # This code combines integrates flag information into the weights.
+        # Both THELI and Megapipe weights must be rescaled by the 
+        # inverse variance because ngmix expects inverse variance maps.  
+        if gal_guess_flag:
+            sig_noise = get_noise(
+                gals[n_e],
+                weight_map,
+                gal_guess_tmp,
+                pixel_scale,
+            )
+        else:
+            sig_noise = sigma_mad(gals[n_e])
+
+        noise_img = np.random.randn(*gals[n_e].shape) * sig_noise
+        noise_img_gal = np.random.randn(*gals[n_e].shape) * sig_noise
+
+        gal_masked = np.copy(gals[n_e])
+        if (len(np.where(weight_map == 0)[0]) != 0):
+            gal_masked[weight_map == 0] = noise_img_gal[weight_map == 0]
+
+        weight_map *= 1 / sig_noise ** 2
+
+        # gaussian fit to the original psf 
+        #(pre-counterfactual image operations)
+        w_tmp = np.sum(weight_map)
+
+        psf_res_gT['g_PSFo'] += psf_res['g'] * w_tmp
+        psf_res_gT['g_err_PSFo'] += np.array([
+            psf_res['pars_err'][2],
+            psf_res['pars_err'][3]
+        ]) * w_tmp
+        psf_res_gT['r50_PSFo'] += psf_res['pars'][4] * w_tmp
+        psf_res_gT['r50_err_PSFo'] += psf_res['pars_err'][4] * w_tmp
+        wsum += w_tmp
+
+        gal_obs = Observation(
+            gal_masked,
+            weight=weight_map,
+            jacobian=gal_jacob,
+            psf=psf_obs,
+            noise=noise_img
+        )
+
+        if gal_guess_flag:
+            gal_guess_tmp[:2] = 0
+            gal_guess.append(gal_guess_tmp)
+
+        gal_obs_list.append(gal_obs)
+        r50_guess_psf.append(psf_r50)
+        gal_guess_flag = True
+
+    if wsum == 0:
+        raise ZeroDivisionError('Sum of weights = 0, division by zero')
+
+    # Normalize PSF fit output
+    for key in psf_res_gT.keys():
+        psf_res_gT[key] /= wsum
+
+    # Gal guess handling
+    fail_get_guess = False
+    if len(gal_guess) == 0:
+        fail_get_guess = True
+        gal_pars = [0., 0., 0., 0., 1, 100]
+    else:
+        gal_pars = np.mean(gal_guess, 0)
+
+    psf_model = 'gauss'
+    gal_model = 'gauss'
+
+    # metacal specific parameters
+    metacal_pars = {
+        'types': ['noshear', '1p', '1m', '2p', '2m'],
+        'step': 0.01,
+        'psf': 'gauss',
+        'fixnoise': True,
+        'use_noise_image': True
+    }
+
+    r50guess = np.mean(r50_guess_psf)
+
+    # retry the fit twice
+    ntry = 2
+    ###### put in options for ngmix.metacal.metacal_bootstrap  obs,
+   # runner,
+   # psf_runner=None,
+   # ignore_failed_psf=True,
+   # rng=None,
+   # **metacal_kws
+#):)
+    obs_dict_mcal = ngmix.metacal.get_all_metacal(gal_obs_list, **metacal_pars)
+    res = {'mcal_flags': 0}
+
+    ntry = 5
+
+    for key in sorted(obs_dict_mcal):
+
+        fres = make_galsimfit(
+            obs_dict_mcal[key],
+            gal_model,
+            gal_pars,
+            prior=prior
+        )
+
+        res['mcal_flags'] |= fres['flags']
+        tres = {}
+
+        for name in fres.keys():
+            tres[name] = fres[name]
+        tres['flags'] = fres['flags']
+
+        wsum = 0
+        r50psf_sum = 0
+        gpsf_sum = np.zeros(2)
+        npsf = 0
+        for obs in obs_dict_mcal[key]:
+
+            if hasattr(obs, 'psf_nopix'):
+                try:
+                    psf_res = make_galsimfit(
+                        obs.psf_nopix,
+                        psf_model,
+                        np.array([0., 0., 0., 0., r50guess, 1.]),
+                        ntry=ntry
+                    )
+                except Exception:
+                    continue
+                g1, g2 = psf_res['g']
+                r50 = psf_res['pars'][4]
+            else:
+                try:
+                    psf_res = make_galsimfit(
+                        obs.psf,
+                        psf_model,
+                        np.array([0., 0., 0., 0., r50guess, 1.]),
+                    )
+                except Exception:
+                    continue
+                g1, g2 = psf_res['g']
+                r50 = psf_res['pars'][4]
+
+            # TODO we sometimes use other weights
+            twsum = obs.weight.sum()
+
+            wsum += twsum
+            gpsf_sum[0] += g1 * twsum
+            gpsf_sum[1] += g2 * twsum
+            r50psf_sum += r50 * twsum
+            npsf += 1
+
+        tres['gpsf'] = gpsf_sum / wsum
+        tres['r50psf'] = r50psf_sum / wsum
+
+        res[key] = tres
+
+    # result dictionary, keyed by the types in metacal_pars above
+    metacal_res = res
+
+    metacal_res.update(psf_res_gT)
+    metacal_res['moments_fail'] = fail_get_guess
+
+    return metacal_res
diff --git a/shapepipe/pipeline/file_handler.py b/shapepipe/pipeline/file_handler.py
index b7479f6f3..a9124fa25 100644
--- a/shapepipe/pipeline/file_handler.py
+++ b/shapepipe/pipeline/file_handler.py
@@ -986,6 +986,7 @@ def _save_num_patterns(
             )
 
         # Save file list
+        os.makedirs(os.path.dirname(output_file), exist_ok=True)
         np.save(output_file, np.array(final_file_list))
 
         del true_file_list, final_file_list

From 922d8aa2a3088557a286f615e68c5688bd09dee4 Mon Sep 17 00:00:00 2001
From: Lucie Baumont <lucie.baumont@cea.fr>
Date: Mon, 13 Feb 2023 17:56:23 +0100
Subject: [PATCH 04/11] bash script to convert theli book-keeping to shapepipe

---
 scripts/sh/find_exposures_theli.sh | 67 ++++++++++++++++++++++++++++++
 1 file changed, 67 insertions(+)
 create mode 100755 scripts/sh/find_exposures_theli.sh

diff --git a/scripts/sh/find_exposures_theli.sh b/scripts/sh/find_exposures_theli.sh
new file mode 100755
index 000000000..60fd16ce5
--- /dev/null
+++ b/scripts/sh/find_exposures_theli.sh
@@ -0,0 +1,67 @@
+#!/usr/bin/env bash
+
+# Name: find_exposures.sh
+# Description: Script to convert THELI images to shapepipe 
+#              bookkeeping conventions
+# Author: Lucie Baumont <lucie.baumont@cea.fr>
+# Date: v1.0 1/2023
+
+
+# Command line arguments
+
+## Help string
+usage="Usage: $(basename "$0") THELI_DIR OUTPUT_DIR\n
+   THELI_DIR\n
+   \thome director of THELI images, eg. /n17data/baumont/THELI_W3/W3/CFISCOLLAB_V2.0.0A\n
+   OUTPUT_DIR\n
+   \troot directory where outputs will go \n
+"
+
+## Help if no arguments
+if [ -z $1 ]; then
+        echo -ne $usage
+        exit 1
+fi
+theli_dir=$1
+output_dir=$2
+
+#make necessary directories
+mkdir -p ${output_dir}/tiles
+mkdir -p ${output_dir}/tiles/exposure_lists
+mkdir -p ${output_dir}/exposures
+mkdir -p ${output_dir}/exposures/headers
+
+for dir in ${theli_dir}/CFIS*;
+    do
+    # create link for all tiles and rename them
+    coadd_dir=$dir/r.MP9602/coadd_V2.0.0A/
+    single_dir=${coadd_dir/coadd/single}
+    echo $single_dir
+    header_dir=${coadd_dir/coadd/headers}
+    file=${coadd_dir}/*.cut.fits
+    base=$(basename ${file} _r.MP9602.V2.0.0A.swarp.cut.fits)
+    num_pattern=${base:5}
+    num=${num_pattern//[-._]/}
+    ln -s ${coadd_dir}/${base}_r.MP9602.V2.0.0A.swarp.cut.fits ${output_dir}/tiles/CFIS_$num.fits
+    ln -s ${coadd_dir}/${base}_r.MP9602.V2.0.0A.swarp.cut.weight.fits ${output_dir}/tiles/CFIS_$num.weight.fits
+    ln -s ${coadd_dir}/${base}_r.MP9602.V2.0.0A.swarp.cut.flag.fits ${output_dir}/tiles/CFIS_$num.flag.fits
+    # create exposure list, link files but exclude broken links
+     for file in $(find ${single_dir}/*.sub.fits -type l -not -xtype l); do
+        base=$(basename $file C.sub.fits)
+        weight_link=${single_dir}/${base}C.weight.fits.gz
+        if [ -L ${weight_link} ] ; then
+            if [ -e ${weight_link} ] ; then
+                echo $(basename $file C.sub.fits)
+                ln -s ${single_dir}/${base}C.sub.fits ${output_dir}/exposures
+                ln -s ${single_dir}/${base}C.weight.fits.gz ${output_dir}/exposures
+                ln -s ${header_dir}/${base}.head ${output_dir}/exposures/headers
+            fi
+         fi 
+         done >> ${output_dir}/tiles/exposure_lists/exp_numbers-$num.txt
+    #list=$(find ${single_dir}/*.sub.fits | awk -F/ '{print substr($NF,1,8)}'|uniq)
+    #echo "${list}">exposure_lists/$num.txt
+     
+    done
+
+   #uncompress step 
+   #gunzip -f ${output_dir}/exposures/*.gz

From 60690b0264c5e0d09eca60361b2a1667ba97d35d Mon Sep 17 00:00:00 2001
From: Lucie Baumont <lucie.baumont@cea.fr>
Date: Mon, 13 Feb 2023 18:55:44 +0100
Subject: [PATCH 05/11] get_images works for theli

---
 .../modules/get_images_package/get_images.py    | 17 +++++++++--------
 1 file changed, 9 insertions(+), 8 deletions(-)

diff --git a/shapepipe/modules/get_images_package/get_images.py b/shapepipe/modules/get_images_package/get_images.py
index f0a4220c6..dd993aaf6 100644
--- a/shapepipe/modules/get_images_package/get_images.py
+++ b/shapepipe/modules/get_images_package/get_images.py
@@ -62,7 +62,6 @@ def in2out_pattern(number):
 
     # remove letters in number
     number_final = re.sub('[a-zA-Z]', '', number_final)
-
     return number_final
 
 
@@ -228,16 +227,18 @@ def get_file_list(
 
             list_files_per_type = []
             for number in image_number_list:
-
+            # find all files with exposure number and extension with glob
+            # loop over this list, make sure it is robust for a single exposure
+            # allow for prefix
                 if use_output_file_pattern:
                     # Transform input to output number patterns
-
                     number_final = in2out_pattern(number)
-
-                    # Keep initial dot in extension
-                    x = in_ext[1:]
-                    x2 = re.sub(r'\.', '', x)
-                    ext_final = in_ext[0] + x2
+                    # shapepipe extensions can only be .fits or .fitsfz
+                    x = re.sub(r'.+(?=\b.fits\b)','',in_ext)
+                    ext_final=re.sub(r'\.(?!fits)','',x)
+                    #x = in_ext[1:]
+                    #x2 = re.sub(r'\.', '', x)
+                    #ext_final = in_ext[0] + x2
                     fbase = (
                         f'{self._output_file_pattern[idx]}{number_final}'
                     )

From 7ada8501da5795c98ebd0a4d54494579346e53e1 Mon Sep 17 00:00:00 2001
From: Lucie Baumont <lucie.baumont@cea.fr>
Date: Tue, 21 Feb 2023 16:23:58 +0100
Subject: [PATCH 06/11] change WCS from TAN to TPV

---
 scripts/sh/find_exposures_theli.sh               | 9 ++++++---
 shapepipe/modules/split_exp_package/split_exp.py | 8 ++++----
 2 files changed, 10 insertions(+), 7 deletions(-)

diff --git a/scripts/sh/find_exposures_theli.sh b/scripts/sh/find_exposures_theli.sh
index 60fd16ce5..c42140704 100755
--- a/scripts/sh/find_exposures_theli.sh
+++ b/scripts/sh/find_exposures_theli.sh
@@ -54,14 +54,17 @@ for dir in ${theli_dir}/CFIS*;
                 echo $(basename $file C.sub.fits)
                 ln -s ${single_dir}/${base}C.sub.fits ${output_dir}/exposures
                 ln -s ${single_dir}/${base}C.weight.fits.gz ${output_dir}/exposures
-                ln -s ${header_dir}/${base}.head ${output_dir}/exposures/headers
+                # remove header line 2 which has french accents that cause problems with astropy, change wcs keys
+                awk 'NR!~/^(2)$/ {sub(/TAN/,"TPV"); print}' ${header_dir}/${base}.head > ${output_dir}/exposures/headers/${base}.head
+                
             fi
          fi 
          done >> ${output_dir}/tiles/exposure_lists/exp_numbers-$num.txt
-    #list=$(find ${single_dir}/*.sub.fits | awk -F/ '{print substr($NF,1,8)}'|uniq)
+   
+     #list=$(find ${single_dir}/*.sub.fits | awk -F/ '{print substr($NF,1,8)}'|uniq)
     #echo "${list}">exposure_lists/$num.txt
      
     done
 
    #uncompress step 
-   #gunzip -f ${output_dir}/exposures/*.gz
+   gunzip -f ${output_dir}/exposures/*.gz
diff --git a/shapepipe/modules/split_exp_package/split_exp.py b/shapepipe/modules/split_exp_package/split_exp.py
index 2dd1a5875..06aedeba9 100644
--- a/shapepipe/modules/split_exp_package/split_exp.py
+++ b/shapepipe/modules/split_exp_package/split_exp.py
@@ -93,7 +93,7 @@ def create_hdus(
         output_suffix : str
             Suffix for the output file
         transf_coord : bool
-            Transform the WCS (``pv`` to ``sip``) if ``True``
+            Transform the WCS (``TAN`` to ``TPV``) if ``True``
         transf_int : bool
             Set data types to int if ``True``
         save_header : bool
@@ -103,11 +103,11 @@ def create_hdus(
         header_file = np.zeros(self._n_hdu, dtype='O')
 
         for idx in range(1, self._n_hdu + 1):
-
             h = fits.getheader(exp_path, idx)
             if transf_coord:
-                stp.pv_to_sip(h)
-
+                # correct WCS convention so astropy can read PVs
+                h['CTYPE1']='RA--TPV'
+                h['CTYPE2']='RA--TPV' 
             d = fits.getdata(exp_path, idx)
             if transf_int:
                 d = d.astype(np.int16)

From e679bbf6e357e96c60a089c5f911d4a4ba70e1da Mon Sep 17 00:00:00 2001
From: Lucie Baumont <lucie.baumont@cea.fr>
Date: Wed, 22 Feb 2023 13:05:23 +0100
Subject: [PATCH 07/11] untested edits to merge_header module

---
 .../modules/merge_headers_package/__init__.py |   2 +-
 .../merge_headers_package/merge_headers.py    | 119 +++++++++++++++---
 shapepipe/modules/merge_headers_runner.py     |   4 +-
 3 files changed, 105 insertions(+), 20 deletions(-)

diff --git a/shapepipe/modules/merge_headers_package/__init__.py b/shapepipe/modules/merge_headers_package/__init__.py
index 48a972c14..c8a0eb962 100644
--- a/shapepipe/modules/merge_headers_package/__init__.py
+++ b/shapepipe/modules/merge_headers_package/__init__.py
@@ -2,7 +2,7 @@
 
 This package contains the module for ``merge_headers``.
 
-:Author: Axel Guinot
+:Author: Axel Guinot, Lucie Baumont
 
 :Parent module: ``split_exp_runner``
 
diff --git a/shapepipe/modules/merge_headers_package/merge_headers.py b/shapepipe/modules/merge_headers_package/merge_headers.py
index 95d3e9acc..65fa5b623 100644
--- a/shapepipe/modules/merge_headers_package/merge_headers.py
+++ b/shapepipe/modules/merge_headers_package/merge_headers.py
@@ -1,10 +1,10 @@
 """MERGE HEADERS.
 
 This module merges the output *header* files of the ``split_exp``
-module. It creates a binnary file that contains the WCS of each CCD for each
-exposure.
+module or external scamp headers. It creates a binary file that 
+contains the WCS of each CCD for each exposure.
 
-:Author: Axel Guinot
+:Author: Axel Guinot, Lucie Baumont
 
 """
 
@@ -13,9 +13,10 @@
 
 import numpy as np
 from sqlitedict import SqliteDict
+from astropy.io.fits import Header
+from astropy.wcs import WCS
 
-
-def merge_headers(input_file_list, output_dir):
+def merge_headers(input_exposure_paths, output_dir, ext_header, nhdu):
     """Merge Headers.
 
     This function opens the files in the input file list and merges them into
@@ -23,11 +24,14 @@ def merge_headers(input_file_list, output_dir):
 
     Parameters
     ----------
-    input_file_list : list
+    input_exposure_paths : list
         List of input files
     output_dir : str
         Output path
-
+    ext_header: bool
+        Use external scamp header if ``True``
+    n_hdu: int
+        number of ccds
     Raises
     ------
     TypeError
@@ -42,24 +46,103 @@ def merge_headers(input_file_list, output_dir):
 
     # Open SqliteDict file
     final_file = SqliteDict(f'{output_dir}/log_exp_headers.sqlite')
-    # Set matching pattern
+    # define file properties for bookkeeping
+    header_dir = os.path.split(input_file_list[0][0])[0]
+    ext = os.path.splitext(input_file_list[0][0])[1]
     pattern = 'headers-'
+    keys = unique_exposure_list(input_file_list, header_dir, pattern, ext, ext_header)
+       
+    for key in keys:
+        # Load header, numpy binary or external header
+        if ext_header:
+            final_file[key] = create_joint_header(n_hdu, key, header_dir, pattern, ext)
+        else:
+            file_path = datadir+pattern+"key"+ext
+            final_file[key] = np.load(file_path, allow_pickle=True)
+    # Close file
+    final_file.commit()
+    final_file.close()
+
+def unique_exposure_list(input_file_list, header_dir, file_pattern, ext, ext_header):
+    """unique_exposure_list.
 
-    for file_path in input_file_list:
-        # Extract file path
-        file_path_scalar = file_path[0]
-        # Check file name against pattern
+    Extract unique exposure ids from file list.
+
+    Parameters
+    ----------
+    input_file_list: list
+        list of input files
+    header_dir: str
+        directory of headers
+    file_pattern: str
+        text in filename 
+    ext: str
+        file extension, eg .npy, .head
+    ext_header: bool
+        Uses external scamp header if ``True`` 
+
+    Returns
+    -------
+        list
+            List of unique exposure numbers
+    """
+    file_list = []
+    if ext_header:
+        ext = '\-\d{1,2}' + ext 
+    # extract keys
+    for file in input_file_list:
+        file_path_scalar = file[0]
         file_name = os.path.split(file_path_scalar)[1]
-        file_base_name = os.path.splitext(file_name)[0]
+        root = re.split(ext, file_name)[0]
         pattern_split = re.split(pattern, file_base_name)
         if len(pattern_split) < 2:
             raise IndexError(
                 f'Regex "{pattern}" not found in base name "{file_base_name}".'
             )
         key = pattern_split[1]
-        # Load Numpy binary file
-        final_file[key] = np.load(file_path_scalar, allow_pickle=True)
+        
+        file_list.append(key)
+    # only keep unique keys
+    unique_keys = list(set(file_list))
 
-    # Close file
-    final_file.commit()
-    final_file.close()
+    return unique_keys
+
+    def create_joint_header(n_hdu, key, header_dir, pattern, ext):
+    """create_joint_header.
+
+    Packages scamp headers for ccds into a numpy array per exposure.
+     
+    Parameters
+    ----------
+    n_hdu: int
+        total number of ccds
+    key: str
+        exposure number
+    header_dir: str
+        directory where headers are stored
+    pattern: str
+        file pattern, e.g. headers-
+    ext: str
+        header extension 
+    Returns
+    -------
+        list
+            compilation of headers corresponding to exposure number
+    """
+
+    header_file = np.zeros(n_hdu, dtype='O')
+    for idx in range(1, n_hdu + 1):
+        filepath= filedir+pattern+key+'-'+str(n_hdu)+ext
+        # check for file                                                
+        if os.path.exists(filepath):
+            h=Header.fromfile(filepath,sep='\n',endcard=False,padding=False)
+            try:
+                w = WCS(h)
+            except Exception:
+                print(f'WCS error for file {exp_path}')
+                raise
+            header_file[idx - 1] = {'WCS': w, 'header': h.tostring()}
+        else:
+            header_file[idx - 1] = 'empty'
+
+    return header_file
diff --git a/shapepipe/modules/merge_headers_runner.py b/shapepipe/modules/merge_headers_runner.py
index 9e692bba7..7dcea1c8d 100644
--- a/shapepipe/modules/merge_headers_runner.py
+++ b/shapepipe/modules/merge_headers_runner.py
@@ -20,6 +20,8 @@
 )
 def merge_headers_runner(
     input_file_list,
+    n_hdu,
+    ext_header,
     run_dirs,
     file_number_string,
     config,
@@ -36,7 +38,7 @@ def merge_headers_runner(
     w_log.info(f'output_dir = {output_dir}')
 
     # Merge header files
-    merge_headers(input_file_list, output_dir)
+    merge_headers(input_file_list, output_dir, ext_header, n_hdu)
 
     # No return objects
     return None, None

From 031a918e21f91f3d87b857191e6d773cd311d635 Mon Sep 17 00:00:00 2001
From: Lucie Baumont <lucie.baumont@cea.fr>
Date: Wed, 22 Feb 2023 17:09:42 +0100
Subject: [PATCH 08/11] convert _ to - in filenames

---
 shapepipe/modules/get_images_package/get_images.py | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/shapepipe/modules/get_images_package/get_images.py b/shapepipe/modules/get_images_package/get_images.py
index dd993aaf6..e7fa6365f 100644
--- a/shapepipe/modules/get_images_package/get_images.py
+++ b/shapepipe/modules/get_images_package/get_images.py
@@ -59,7 +59,8 @@ def in2out_pattern(number):
     # replace dots ('.') with dashes ('-') to avoid confusion
     # with file extension delimiters
     number_final = re.sub(r'\.', '-', number)
-
+    # replace underscore with dashes
+    number_final = re.sub(r'_', '-', number_final)
     # remove letters in number
     number_final = re.sub('[a-zA-Z]', '', number_final)
     return number_final
@@ -233,12 +234,11 @@ def get_file_list(
                 if use_output_file_pattern:
                     # Transform input to output number patterns
                     number_final = in2out_pattern(number)
-                    # shapepipe extensions can only be .fits or .fitsfz
+                    # remove leading suffix from suffix.fits files
                     x = re.sub(r'.+(?=\b.fits\b)','',in_ext)
-                    ext_final=re.sub(r'\.(?!fits)','',x)
-                    #x = in_ext[1:]
-                    #x2 = re.sub(r'\.', '', x)
-                    #ext_final = in_ext[0] + x2
+                    # remove all but leading .
+                    ext_final='.'+re.sub(r'\.','',x)
+                   
                     fbase = (
                         f'{self._output_file_pattern[idx]}{number_final}'
                     )

From 5daad5d1d00ead88808677d0a615a404a76b7625 Mon Sep 17 00:00:00 2001
From: Lucie Baumont <lucie.baumont@cea.fr>
Date: Wed, 22 Feb 2023 18:22:50 +0100
Subject: [PATCH 09/11] merge headers changes run so far

---
 .../merge_headers_package/merge_headers.py    | 45 ++++++++++---------
 shapepipe/modules/merge_headers_runner.py     |  8 ++--
 2 files changed, 27 insertions(+), 26 deletions(-)

diff --git a/shapepipe/modules/merge_headers_package/merge_headers.py b/shapepipe/modules/merge_headers_package/merge_headers.py
index 65fa5b623..7e9761410 100644
--- a/shapepipe/modules/merge_headers_package/merge_headers.py
+++ b/shapepipe/modules/merge_headers_package/merge_headers.py
@@ -16,7 +16,7 @@
 from astropy.io.fits import Header
 from astropy.wcs import WCS
 
-def merge_headers(input_exposure_paths, output_dir, ext_header, nhdu):
+def merge_headers(input_file_list, output_dir, pattern, ext_header, n_hdu):
     """Merge Headers.
 
     This function opens the files in the input file list and merges them into
@@ -24,13 +24,15 @@ def merge_headers(input_exposure_paths, output_dir, ext_header, nhdu):
 
     Parameters
     ----------
-    input_exposure_paths : list
+    input_file_list : list
         List of input files
     output_dir : str
         Output path
-    ext_header: bool
+    pattern : str
+        File pattern
+    ext_header : bool
         Use external scamp header if ``True``
-    n_hdu: int
+    n_hdu : int
         number of ccds
     Raises
     ------
@@ -49,7 +51,6 @@ def merge_headers(input_exposure_paths, output_dir, ext_header, nhdu):
     # define file properties for bookkeeping
     header_dir = os.path.split(input_file_list[0][0])[0]
     ext = os.path.splitext(input_file_list[0][0])[1]
-    pattern = 'headers-'
     keys = unique_exposure_list(input_file_list, header_dir, pattern, ext, ext_header)
        
     for key in keys:
@@ -70,21 +71,21 @@ def unique_exposure_list(input_file_list, header_dir, file_pattern, ext, ext_hea
 
     Parameters
     ----------
-    input_file_list: list
+    input_file_list : list
         list of input files
-    header_dir: str
+    header_dir : str
         directory of headers
-    file_pattern: str
+    file_pattern : str
         text in filename 
-    ext: str
+    ext : str
         file extension, eg .npy, .head
-    ext_header: bool
+    ext_header : bool
         Uses external scamp header if ``True`` 
 
     Returns
     -------
-        list
-            List of unique exposure numbers
+    list
+        List of unique exposure numbers
     """
     file_list = []
     if ext_header:
@@ -94,7 +95,7 @@ def unique_exposure_list(input_file_list, header_dir, file_pattern, ext, ext_hea
         file_path_scalar = file[0]
         file_name = os.path.split(file_path_scalar)[1]
         root = re.split(ext, file_name)[0]
-        pattern_split = re.split(pattern, file_base_name)
+        pattern_split = re.split(file_pattern, root)
         if len(pattern_split) < 2:
             raise IndexError(
                 f'Regex "{pattern}" not found in base name "{file_base_name}".'
@@ -107,32 +108,32 @@ def unique_exposure_list(input_file_list, header_dir, file_pattern, ext, ext_hea
 
     return unique_keys
 
-    def create_joint_header(n_hdu, key, header_dir, pattern, ext):
+def create_joint_header(n_hdu, key, header_dir, pattern, ext):
     """create_joint_header.
 
     Packages scamp headers for ccds into a numpy array per exposure.
      
     Parameters
     ----------
-    n_hdu: int
+    n_hdu : int
         total number of ccds
-    key: str
+    key : str
         exposure number
-    header_dir: str
+    header_dir : str
         directory where headers are stored
-    pattern: str
+    pattern : str
         file pattern, e.g. headers-
-    ext: str
+    ext : str
         header extension 
     Returns
     -------
-        list
-            compilation of headers corresponding to exposure number
+    list
+        compilation of headers corresponding to exposure number
     """
 
     header_file = np.zeros(n_hdu, dtype='O')
     for idx in range(1, n_hdu + 1):
-        filepath= filedir+pattern+key+'-'+str(n_hdu)+ext
+        filepath= header_dir+pattern+key+'-'+str(n_hdu)+ext
         # check for file                                                
         if os.path.exists(filepath):
             h=Header.fromfile(filepath,sep='\n',endcard=False,padding=False)
diff --git a/shapepipe/modules/merge_headers_runner.py b/shapepipe/modules/merge_headers_runner.py
index 7dcea1c8d..dc41c98a5 100644
--- a/shapepipe/modules/merge_headers_runner.py
+++ b/shapepipe/modules/merge_headers_runner.py
@@ -20,8 +20,6 @@
 )
 def merge_headers_runner(
     input_file_list,
-    n_hdu,
-    ext_header,
     run_dirs,
     file_number_string,
     config,
@@ -33,12 +31,14 @@ def merge_headers_runner(
     output_dir = run_dirs['output']
     if config.has_option(module_config_sec, 'OUTPUT_PATH'):
         output_dir = config.getexpanded(module_config_sec, 'OUTPUT_PATH')
-
+    file_pattern = config.get(module_config_sec, 'FILE_PATTERN')
+    ext_header = config.get(module_config_sec, 'EXT_HEADER')
+    n_hdu = int(config.get(module_config_sec, 'N_HDU'))
     # Log output directory
     w_log.info(f'output_dir = {output_dir}')
 
     # Merge header files
-    merge_headers(input_file_list, output_dir, ext_header, n_hdu)
+    merge_headers(input_file_list, output_dir, file_pattern, ext_header, n_hdu)
 
     # No return objects
     return None, None

From d6f4d949b34e3cb254b17d888c70908866120194 Mon Sep 17 00:00:00 2001
From: Lucie Baumont <lucie.baumont@cea.fr>
Date: Thu, 23 Feb 2023 13:41:48 +0100
Subject: [PATCH 10/11] merge_headers.py runs

---
 shapepipe/modules/merge_headers_package/merge_headers.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/shapepipe/modules/merge_headers_package/merge_headers.py b/shapepipe/modules/merge_headers_package/merge_headers.py
index 7e9761410..c3f22f3a2 100644
--- a/shapepipe/modules/merge_headers_package/merge_headers.py
+++ b/shapepipe/modules/merge_headers_package/merge_headers.py
@@ -133,7 +133,8 @@ def create_joint_header(n_hdu, key, header_dir, pattern, ext):
 
     header_file = np.zeros(n_hdu, dtype='O')
     for idx in range(1, n_hdu + 1):
-        filepath= header_dir+pattern+key+'-'+str(n_hdu)+ext
+        filepath= header_dir+'/'+pattern+key+'-'+str(idx)+ext
+        print(filepath)
         # check for file                                                
         if os.path.exists(filepath):
             h=Header.fromfile(filepath,sep='\n',endcard=False,padding=False)

From b81a17a24b4f682efa4a6138522055da3793925c Mon Sep 17 00:00:00 2001
From: Lucie Baumont <lucie.baumont@cea.fr>
Date: Thu, 23 Feb 2023 13:58:25 +0100
Subject: [PATCH 11/11] remove file that was added by accident

---
 shapepipe/modules/ngmix_package/#ngmix.py# | 1029 --------------------
 1 file changed, 1029 deletions(-)
 delete mode 100644 shapepipe/modules/ngmix_package/#ngmix.py#

diff --git a/shapepipe/modules/ngmix_package/#ngmix.py# b/shapepipe/modules/ngmix_package/#ngmix.py#
deleted file mode 100644
index fa7d15a0f..000000000
--- a/shapepipe/modules/ngmix_package/#ngmix.py#
+++ /dev/null
@@ -1,1029 +0,0 @@
-"""NGMIX.
-
-This module contains a class for ngmix shape measurement.
-
-:Author: Axel Guinot
-
-"""
-
-import re
-
-import galsim
-import ngmix
-import numpy as np
-from astropy.io import fits
-from modopt.math.stats import sigma_mad
-from ngmix.observation import MultiBandObsList, Observation, ObsList
-from numpy.random import uniform as urand
-from sqlitedict import SqliteDict
-
-from shapepipe.pipeline import file_io
-
-
-class Ngmix(object):
-    """Ngmix.
-
-    Class to handle NGMIX shapepe measurement.
-
-    Parameters
-    ----------
-    input_file_list : list
-        Input files
-    output_dir : str
-        Output directory
-    file_number_string : str
-        File numbering scheme
-    zero_point : float
-        Photometric zero point
-    pixel_scale : float
-        Pixel scale in arcsec
-    f_wcs_path : str
-        Path to merged single-exposure single-HDU headers
-    w_log : logging.Logger
-        Logging instance
-    id_obj_min : int, optional
-        First galaxy ID to process, not used if the value is set to ``-1``;
-        the default is ``-1``
-    id_obj_max : int, optional
-        Last galaxy ID to process, not used if the value is set to ``-1``;
-        the default is ``-1``
-
-    Raises
-    ------
-    IndexError
-        If the length of the input file list is incorrect
-
-    """
-
-    def __init__(
-        self,
-        input_file_list,
-        output_dir,
-        file_number_string,
-        zero_point,
-        pixel_scale,
-        f_wcs_path,
-        w_log,
-        id_obj_min=-1,
-        id_obj_max=-1
-    ):
-
-        if len(input_file_list) != 6:
-            raise IndexError(
-                f'Input file list has length {len(input_file_list)},'
-                + ' required is 6'
-            )
-
-        self._tile_cat_path = input_file_list[0]
-        self._gal_vignet_path = input_file_list[1]
-        self._bkg_vignet_path = input_file_list[2]
-        self._psf_vignet_path = input_file_list[3]
-        self._weight_vignet_path = input_file_list[4]
-        self._flag_vignet_path = input_file_list[5]
-
-        self._output_dir = output_dir
-        self._file_number_string = file_number_string
-
-        self._zero_point = zero_point
-        self._pixel_scale = pixel_scale
-
-        self._f_wcs_path = f_wcs_path
-        self._id_obj_min = id_obj_min
-        self._id_obj_max = id_obj_max
-
-        self._w_log = w_log
-
-        # Initiatlise random generator
-        seed = int(''.join(re.findall(r'\d+', self._file_number_string)))
-        np.random.seed(seed)
-        self._w_log.info(f'Random generator initialisation seed = {seed}')
-
-    @classmethod
-    def MegaCamFlip(self, vign, ccd_nb):
-        """Flip for MegaCam.
-
-        MegaPipe has CCDs that are upside down. This function flips the
-        postage stamps in these CCDs. TO DO: This will give incorrect results
-        when used with THELI ccds.  Fix this.
-
-        Parameters
-        ----------
-        vign : numpy.ndarray
-            Array containing the postage stamp to flip
-        ccd_nb : int
-            ID of the CCD containing the postage stamp
-
-        Returns
-        -------
-        numpy.ndarray
-            The flipped postage stamp
-
-        """
-        if ccd_nb < 18 or ccd_nb in [36, 37]:
-            # swap x axis so origin is on top-right
-            return np.rot90(vign, k=2)
-            #print('rotating megapipe image')
-        else:
-            # swap y axis so origin is on bottom-left
-            return vign
-
-    def get_prior(self):
-        """Get Prior.
-
-        Return prior for the different parameters.
-        Prior on ellipticity is from Bernstein and Armstrong 2014.
-        Prior on centers is a 2d gaussian.
-        Priors on size and flux are flat.
-         
-
-        Returns
-        -------
-        ngmix.priors
-            Priors for the different parameters (ellipticity,center, size, flux)
-        """
-        # Prior on ellipticity. Details do not matter, as long
-        # as it regularizes the fit. From Bernstein & Armstrong 2014
-        g_sigma = 0.4
-        rng = np.random.RandomState(22)
-        g_prior = ngmix.priors.GPriorBA(g_sigma,rng)
-
-        # 2-d Gaussian prior on the center row and column center
-        # (relative to the center of the jacobian, which
-        # would be zero) and the sigma of the Gaussians.
-        # Units same as jacobian, probably arcsec
-        row, col = 0.0, 0.0
-        row_sigma, col_sigma = self._pixel_scale, self._pixel_scale
-        cen_prior = ngmix.priors.CenPrior(row, col, row_sigma, col_sigma, rng)
-
-        # Size prior. Instead of flat, two-sided error function (TwoSidedErf)
-        # could be used, though they are not used here
-        r50minval = -10.0  # arcsec squared
-        r50maxval = 1.e6
-        r50_prior = ngmix.priors.FlatPrior(r50minval, r50maxval, rng)
-
-        # Flux prior. Bounds need to make sense for
-        # images in question
-        Fminval = -1.e4
-        Fmaxval = 1.e9
-        F_prior = ngmix.priors.FlatPrior(Fminval, Fmaxval, rng)
-
-        # Joint prior, combine all individual priors
-        prior = ngmix.joint_prior.PriorSimpleSep(
-            cen_prior,
-            g_prior,
-            r50_prior,
-            F_prior
-        )
-
-        return prior
-
-    def compile_results(self, results):
-        """Compile Results.
-
-        Prepare the results of NGMIX before saving. TODO: add snr_r and T_r
-
-        Parameters
-        ----------
-        results : dict
-            Results of NGMIX metacal
-
-        Returns
-        -------
-        dict
-            Compiled results ready to be written to a file
-            note: psfo is the original image psf from psfex or mccd
-
-        Raises
-        ------
-        KeyError
-            If SNR key not found
-
-        """
-        names = ['1m', '1p', '2m', '2p', 'noshear']
-        names2 = [
-            'id',
-            'n_epoch_model',
-            'moments_fail',
-            'ntry_fit',
-            'g1_psfo_ngmix',
-            'g2_psfo_ngmix',
-            'r50_psfo_ngmix',
-            'g1_err_psfo_ngmix',
-            'g2_err_psfo_ngmix',
-            'r50_err_psfo_ngmix',
-            'g1',
-            'g1_err',
-            'g2',
-            'g2_err',
-            'r50',
-            'r50_err',
-            'r50psf',
-            'g1_psf',
-            'g2_psf',
-            'flux',
-            'flux_err',
-            's2n',
-            'mag',
-            'mag_err',
-            'flags',
-            'mcal_flags'
-        ]
-        output_dict = {k: {kk: [] for kk in names2} for k in names}
-        for idx in range(len(results)):
-            for name in names:
-
-                mag = (
-                    -2.5 * np.log10(results[idx][name]['flux'])
-                    + self._zero_point
-                )
-                mag_err = np.abs(
-                    -2.5 * results[idx][name]['flux_err']
-                    / (results[idx][name]['flux'] * np.log(10))
-                )
-
-                output_dict[name]['id'].append(results[idx]['obj_id'])
-                output_dict[name]['n_epoch_model'].append(
-                    results[idx]['n_epoch_model']
-                )
-                output_dict[name]['moments_fail'].append(
-                    results[idx]['moments_fail']
-                )
-                output_dict[name]['ntry_fit'].append(
-                    results[idx][name]['ntry']
-                )
-                output_dict[name]['g1_psfo_ngmix'].append(
-                    results[idx]['g_PSFo'][0]
-                )
-                output_dict[name]['g2_psfo_ngmix'].append(
-                    results[idx]['g_PSFo'][1]
-                )
-                output_dict[name]['g1_err_psfo_ngmix'].append(
-                    results[idx]['g_err_PSFo'][0]
-                )
-                output_dict[name]['g2_err_psfo_ngmix'].append(
-                    results[idx]['g_err_PSFo'][1]
-                )
-                output_dict[name]['r50_psfo_ngmix'].append(
-                    results[idx]['r50_PSFo']
-                )
-                output_dict[name]['r50_err_psfo_ngmix'].append(
-                    results[idx]['r50_err_PSFo']
-                )
-                output_dict[name]['g1'].append(results[idx][name]['g'][0])
-                output_dict[name]['g1_err'].append(
-                    results[idx][name]['pars_err'][2]
-                )
-                output_dict[name]['g2'].append(results[idx][name]['g'][1])
-                output_dict[name]['g2_err'].append(
-                    results[idx][name]['pars_err'][3]
-                )
-                output_dict[name]['r50'].append(results[idx][name]['pars'][4])
-                output_dict[name]['r50_err'].append(results[idx][name]['pars_err'][4])
-                output_dict[name]['r50psf'].append(results[idx][name]['r50psf'])
-                output_dict[name]['g1_psf'].append(
-                    results[idx][name]['gpsf'][0]
-                )
-                output_dict[name]['g2_psf'].append(
-                    results[idx][name]['gpsf'][1]
-                )
-                output_dict[name]['flux'].append(results[idx][name]['flux'])
-                output_dict[name]['flux_err'].append(
-                    results[idx][name]['flux_err']
-                )
-                output_dict[name]['mag'].append(mag)
-                output_dict[name]['mag_err'].append(mag_err)
-
-                if 's2n' in results[idx][name]:
-                    output_dict[name]['s2n'].append(results[idx][name]['s2n'])
-                elif 's2n_r' in results[idx][name]:
-                    output_dict[name]['s2n'].append(
-                        results[idx][name]['s2n_r']
-                    )
-                else:
-                    raise KeyError('No SNR key (s2n, s2n_r) found in results')
-
-                output_dict[name]['flags'].append(results[idx][name]['flags'])
-                output_dict[name]['mcal_flags'].append(
-                    results[idx]['mcal_flags']
-                )
-
-        return output_dict
-
-    def save_results(self, output_dict):
-        """Save Results.
-
-        Save the results into a FITS file.
-
-        Parameters
-        ----------
-        output_dict
-            Dictionary containing the results
-
-        """
-        output_name = (
-            f'{self._output_dir}/ngmix{self._file_number_string}.fits'
-        )
-
-        f = file_io.FITSCatalogue(
-            output_name,
-            open_mode=file_io.BaseCatalogue.OpenMode.ReadWrite
-        )
-
-        for key in output_dict.keys():
-            f.save_as_fits(output_dict[key], ext_name=key.upper())
-
-    def process(self):
-        """Process.
-
-        Funcion to processs NGMIX.
-        organizes object cutouts from detection catalog in image, 
-        weight, and flag files
-        per object: 
-            gathers wcs and psf info from exposures
-            background subtracts (make this an option)
-            scales by relative zeropoints
-            runs metacal convolutions and ngmix fitting
-        Returns
-        -------
-        dict
-            Dictionary containing the NGMIX metacal results
-
-        """
-        # sextractor detection catalog for the tile
-        tile_cat = file_io.FITSCatalogue(
-            self._tile_cat_path,
-            SEx_catalogue=True,
-        )
-        tile_cat.open()
-        obj_id = np.copy(tile_cat.get_data()['NUMBER'])
-        tile_vign = np.copy(tile_cat.get_data()['VIGNET'])
-        tile_ra = np.copy(tile_cat.get_data()['XWIN_WORLD'])
-        tile_dec = np.copy(tile_cat.get_data()['YWIN_WORLD'])
-        tile_cat.close()
-
-        f_wcs_file = SqliteDict(self._f_wcs_path)
-        gal_vign_cat = SqliteDict(self._gal_vignet_path)
-        bkg_vign_cat = SqliteDict(self._bkg_vignet_path)
-        psf_vign_cat = SqliteDict(self._psf_vignet_path)
-        weight_vign_cat = SqliteDict(self._weight_vignet_path)
-        flag_vign_cat = SqliteDict(self._flag_vignet_path)
-
-        final_res = []
-        prior = self.get_prior()
-
-        count = 0
-        id_first = -1
-        id_last = -1
-
-        for i_tile, id_tmp in enumerate(obj_id):
-            if self._id_obj_min > 0 and id_tmp < self._id_obj_min:
-                continue
-            if self._id_obj_max > 0 and id_tmp > self._id_obj_max:
-                continue
-
-            if id_first == -1:
-                id_first = id_tmp
-            id_last = id_tmp
-
-            count = count + 1
-
-            gal_vign = []
-            psf_vign = []
-            sigma_psf = []
-            weight_vign = []
-            flag_vign = []
-            jacob_list = []
-            if (
-                (psf_vign_cat[str(id_tmp)] == 'empty')
-                or (gal_vign_cat[str(id_tmp)] == 'empty')
-            ):
-                continue
-            #identify exposure and ccd number from psf catalog
-            psf_expccd_name = list(psf_vign_cat[str(id_tmp)].keys())
-            for expccd_name_tmp in psf_expccd_name:
-                exp_name, ccd_n = re.split('-', expccd_name_tmp)
-
-                gal_vign_tmp = (
-                    gal_vign_cat[str(id_tmp)][expccd_name_tmp]['VIGNET']
-                )
-                if len(np.where(gal_vign_tmp.ravel() == 0)[0]) != 0:
-                    continue
-                # background subtraction
-                bkg_vign_tmp = (
-                    bkg_vign_cat[str(id_tmp)][expccd_name_tmp]['VIGNET']
-                )
-                gal_vign_sub_bkg = gal_vign_tmp - bkg_vign_tmp
-                # skip this step for THELI
-                tile_vign_tmp = (
-                    Ngmix.MegaCamFlip(np.copy(tile_vign[i_tile]), int(ccd_n))
-                )
-
-                flag_vign_tmp = (
-                    flag_vign_cat[str(id_tmp)][expccd_name_tmp]['VIGNET']
-                )
-                flag_vign_tmp[np.where(tile_vign_tmp == -1e30)] = 2**10
-                v_flag_tmp = flag_vign_tmp.ravel()
-                if len(np.where(v_flag_tmp != 0)[0]) / (51 * 51) > 1 / 3.0:
-                    continue
-
-                weight_vign_tmp = (
-                    weight_vign_cat[str(id_tmp)][expccd_name_tmp]['VIGNET']
-                )
-
-                jacob_tmp = get_jacob(
-                    f_wcs_file[exp_name][int(ccd_n)]['WCS'],
-                    tile_ra[i_tile],
-                    tile_dec[i_tile]
-                )
-
-                header_tmp = fits.Header.fromstring(
-                    f_wcs_file[exp_name][int(ccd_n)]['header']
-                )
-                #rescale images and weights by relative flux scale
-                # this should be it's own function
-                Fscale = header_tmp['FSCALE']
-
-                gal_vign_scaled = gal_vign_sub_bkg * Fscale
-                weight_vign_scaled = weight_vign_tmp * 1 / Fscale ** 2
-
-                gal_vign.append(gal_vign_scaled)
-                psf_vign.append(
-                    psf_vign_cat[str(id_tmp)][expccd_name_tmp]['VIGNET']
-                )
-                sigma_psf.append(
-                    psf_vign_cat[
-                        str(id_tmp)
-                    ][expccd_name_tmp]['SHAPES']['SIGMA_PSF_HSM']
-                )
-                weight_vign.append(weight_vign_scaled)
-                flag_vign.append(flag_vign_tmp)
-                jacob_list.append(jacob_tmp)
-                
-            #if object is observed, carry out metacal operations and run ngmix
-            if len(gal_vign) == 0:
-                continue
-            try:
-                res = do_ngmix_metacal(
-                    gal_vign,
-                    psf_vign,
-                    sigma_psf,
-                    weight_vign,
-                    flag_vign,
-                    jacob_list,
-                    prior,
-                    self._pixel_scale
-                )
-
-            except Exception as ee:
-                self._w_log.info(
-                    f'ngmix failed for object ID={id_tmp}.\nMessage: {ee}'
-                )
-                continue
-
-            res['obj_id'] = id_tmp
-            res['n_epoch_model'] = len(gal_vign)
-            final_res.append(res)
-
-        self._w_log.info(
-            f'ngmix loop over objects finished, processed {count} '
-            + f'objects, id first/last={id_first}/{id_last}'
-        )
-
-        f_wcs_file.close()
-        gal_vign_cat.close()
-        bkg_vign_cat.close()
-        flag_vign_cat.close()
-        weight_vign_cat.close()
-        psf_vign_cat.close()
-
-        # Put all results together
-        res_dict = self.compile_results(final_res)
-
-        # Save results
-        self.save_results(res_dict)
-
-
-def get_guess_shapeHSM(
-    img,
-    pixel_scale,
-    guess_flux_unit='img',
-    guess_size_type='T',
-    guess_size_unit='sky',
-    guess_centroid=True,
-    guess_centroid_unit='sky'
-):
-    r"""Get Guess using galsim hsm shapes.
-
-    Get the guess vector for the NGMIX shape measurement
-    ``[center_x, center_y, g1, g2, size_T, flux]``.
-    No guesses are given for the ellipticity ``(0, 0)``.
-
-    Parameters
-    ----------
-    img : numpy.ndarray
-        Array containing the image
-    pixel_scale : float
-        Approximation of the pixel scale
-    guess_flux_unit : str
-        If ``img`` returns the flux in pixel units, otherwise if ``sky``
-        returns the flux in :math:`{\rm arcsec}^{-2}`
-    guess_size_type : str
-        If ``T`` returns the size in quadrupole moments definition
-        :math:`2\sigma^2`, if ``sigma`` returns the moments
-        :math:`\sigma`, if ``r50`` returns the half light radius
-    guess_size_unit : str
-        If ``img`` returns the size in pixel units, otherwise if ``sky``
-        returns the size in arcsec
-    guess_centroid : bool
-        If ``True``, will return a guess on the object centroid, otherwise if
-        ``False``, will return the image centre
-    guess_centroid_unit : str
-        If ``img`` returns the centroid in pixel unit, otherwise if ``sky``
-        returns the centroid in arcsec
-
-    Returns
-    -------
-    numpy.ndarray
-        Return the guess array ``[center_x, center_y, g1, g2, size_T, flux]``
-
-    Raises
-    ------
-    GalSimHSMError
-        For an error in the computation of adaptive moments
-    ValueError
-        For invalid unit guess types
-
-    """
-    galsim_img = galsim.Image(img, scale=pixel_scale)
-
-    hsm_shape = galsim.hsm.FindAdaptiveMom(galsim_img, strict=False)
-
-    error_msg = hsm_shape.error_message
-                
-    if error_msg != '':
-        raise galsim.hsm.GalSimHSMError(
-            f'Error in adaptive moments :\n{error_msg}'
-        )
-
-    if guess_flux_unit == 'img':
-        guess_flux = hsm_shape.moments_amp
-    elif guess_flux_unit == 'sky':
-        guess_flux = hsm_shape.moments_amp / pixel_scale ** 2
-    else:
-        raise ValueError(
-            f'invalid guess_flux_unit \'{guess_flux_unit}\','
-            + ' must be one of \'img\', \'sky\''
-        )
-
-    if guess_size_unit == 'img':
-        size_unit = 1.
-    elif guess_size_unit == 'sky':
-        size_unit = pixel_scale
-    else:
-        raise ValueError(
-            'invalid guess_size_unit \'{guess_size_unit}\','
-            + 'must be one of \'img\', \'sky\''
-        )
-
-    if guess_size_type == 'sigma':
-        guess_size = hsm_shape.moments_sigma * size_unit
-    elif guess_size_type == 'T':
-        guess_size = 2 * (hsm_shape.moments_sigma * size_unit) ** 2
-    elif guess_size_type == 'r50':
-        guess_size = hsm_shape.moments_sigma * size_unit * 1.17741002252
-    else:
-        raise ValueError(
-            'invalid guess_size_type \'{guess_size_type}\','
-            + 'must be one of \'sigma\', \'T\', or \'r50\''
-        )
-
-    if guess_centroid_unit == 'img':
-        centroid_unit = 1
-    elif guess_centroid_unit == 'sky':
-        centroid_unit = pixel_scale
-    else:
-        raise ValueError(
-            f'invalid guess_centroid_unit \'{guess_centroid_unit}\','
-            + '  must be one of \'img\', \'sky\''
-        )
-
-    if guess_centroid:
-        guess_centroid = (
-            (hsm_shape.moments_centroid - galsim_img.center) * centroid_unit
-        )
-    else:
-        guess_centroid = galsim_img.center * centroid_unit
-
-    guess = np.array([
-        guess_centroid.x,
-        guess_centroid.y,
-        0.,
-        0.,
-        guess_size,
-        guess_flux
-    ])
-
-    return guess
-
-
-def make_galsimfit(obs, model, guess0, prior=None, ntry=5):
-    """Make GalSim Fit.
-
-    wrapper for ngmix image fit using simple GalSim model.
-
-    Parameters
-    ----------
-    obs : ngmix.observation.Observation
-        Image to fit
-    model : str
-        Model for fit
-    guess0 : numpy.ndarray
-        Parameters of first model guess
-    prior : ngmix.prior, optional
-        Prior for fit parameters
-    ntry : int, optional
-        Number of tries for fit, the default is ``5``
-
-    Returns
-    -------
-    dict
-        Results
-
-    Raises
-    ------
-    ngmix.BootGalFailure
-        Failure to bootstrap galaxy
-
-    """
-
-    limit = 0.1
-
-    guess = np.copy(guess0)
-    fres = {}
-    for it in range(ntry):
-        guess[0:5] += urand(low=-limit, high=limit)
-        guess[5:] *= (1 + urand(low=-limit, high=limit))
-        fres['flags'] = 1
-        #OLD_LM_PARS = {"maxfev": 1000, "ftol": 1.0e6, "xtol": 1.0e-6}
-        try:
-            fitter = ngmix.fitting.GalsimFitter(
-                model=model,
-                prior=prior
-            )
-            fres = fitter.go(obs=obs,guess=guess)
-            #print(fres.items())
-        except Exception:
-            continue
-
-        if fres['flags'] == 0:
-            break
-
-    if fres['flags'] != 0:
-        raise ngmix.gexceptions.BootGalFailure(
-            'Failed to fit galaxy image with galsimfit'
-        )
-
-    fres['ntry'] = it + 1
-    return fres
-
-
-def get_jacob(wcs, ra, dec):
-    """Get Jacobian.
-
-    Return the Jacobian of the WCS at the required position.
-
-    Parameters
-    ----------
-    wcs : astropy.wcs.WCS
-        WCS object for which we want the Jacobian
-    ra : float
-        RA position of the center of the vignet (in degrees)
-    dec : float
-        Dec position of the center of the vignet (in degress)
-
-    Returns
-    -------
-    galsim.wcs.BaseWCS.jacobian
-        Jacobian of the WCS at the required position
-
-    """
-    g_wcs = galsim.fitswcs.AstropyWCS(wcs=wcs)
-    world_pos = galsim.CelestialCoord(
-        ra=ra * galsim.angle.degrees,
-        dec=dec * galsim.angle.degrees,
-    )
-    galsim_jacob = g_wcs.jacobian(world_pos=world_pos)
-
-    return galsim_jacob
-
-
-def get_noise(gal, weight, guess, pixel_scale, thresh=1.2):
-    r"""Get Noise.
-
-    Compute the sigma of the noise from an object postage stamp.
-    Use a guess on the object size, ellipticity and flux to create a window
-    function.
-
-    Parameters
-    ----------
-    gal : numpy.ndarray
-        Galaxy image
-    weight : numpy.ndarray
-        Weight image
-    guess : list
-        Gaussian parameters fot the window function
-        ``[x0, y0, g1, g2, r50, flux]``
-    pixel_scale : float
-        Pixel scale of the galaxy image
-    thresh : float, optional
-        Threshold to cut the window function,
-        cut = ``thresh`` * :math:`\sigma_{\rm noise}`;  the default is ``1.2``
-
-    Returns
-    -------
-    float
-        Sigma of the noise on the galaxy image
-
-    """
-    img_shape = gal.shape
-    m_weight = weight != 0
-
-    sig_tmp = sigma_mad(gal[m_weight])
-
-    gauss_win = galsim.Gaussian(sigma=np.sqrt(guess[4] / 2), flux=guess[5])
-    gauss_win = gauss_win.shear(g1=guess[2], g2=guess[3])
-    gauss_win = gauss_win.drawImage(
-        nx=img_shape[0],
-        ny=img_shape[1],
-        scale=pixel_scale
-    ).array
-
-    m_weight = weight[gauss_win < thresh * sig_tmp] != 0
-
-    sig_noise = sigma_mad(gal[gauss_win < thresh * sig_tmp][m_weight])
-
-    return sig_noise
-
-
-def do_ngmix_metacal(
-    gals,
-    psfs,
-    psfs_sigma,
-    weights,
-    flags,
-    jacob_list,
-    prior,
-    pixel_scale
-):
-    """Do Ngmix Metacal.
-
-    Performs  metacalibration on a multi-epoch object and return the joint
-    shape measurement with NGMIX.
-
-    Parameters
-    ----------
-    gals : list
-        List of the galaxy vignets.  List indices run over epochs
-    psfs : list
-        List of the PSF vignets
-    psfs_sigma : list
-        List of the sigma PSFs
-    weights : list
-        List of the weight vignets
-    flags : list
-        List of the flag vignets
-    jacob_list : list
-        List of the Jacobians
-    prior : ngmix.priors
-        Priors for the fitting parameters
-    pixel_scale : float
-        pixel scale in arcsec
-
-    Returns
-    -------
-    dict
-        Dictionary containing the results of NGMIX metacal
-
-    """
-    n_epoch = len(gals)
-
-    if n_epoch == 0:
-        raise ValueError("0 epoch to process")
-
-    # Construct observation objects to pass to ngmix 
-    gal_obs_list = ObsList()
-    r50_guess_psf = []
-    psf_res_gT = {
-        'g_PSFo': np.array([0., 0.]),
-        'g_err_PSFo': np.array([0., 0.]),
-        'r50_PSFo': 0.,
-        'r50_err_PSFo': 0.
-    }
-    gal_guess = []
-    gal_guess_flag = True
-    wsum = 0
-    for n_e in range(n_epoch):
-
-        psf_jacob = ngmix.Jacobian(
-            row=(psfs[0].shape[0] - 1) / 2,
-            col=(psfs[0].shape[1] - 1) / 2,
-            wcs=jacob_list[n_e]
-        )
-        # psf observation is part of ngmix observation
-        psf_obs = Observation(psfs[n_e], jacobian=psf_jacob)
-        # convert sigma to T, which is the half-light radius
-        # NOT the usual definition T=2*sigma^2
-        psf_r50 = psfs_sigma[n_e] * 1.17741  * pixel_scale
-
-        # integrate flag info into weights
-        weight_map = np.copy(weights[n_e])
-        weight_map[np.where(flags[n_e] != 0)] = 0.
-        weight_map[weight_map != 0] = 1
-
-        # fit gaussian to psf
-        psf_guess = np.array([0., 0., 0., 0., psf_r50, 1.])
-        try:
-            psf_res = make_galsimfit(psf_obs, 'gauss', psf_guess)
-        except Exception:
-            continue
-
-        # Gal guess
-        try:
-            gal_guess_tmp = get_guess_shapeHSM(
-                gals[n_e],
-                pixel_scale,
-                guess_size_type='sigma'
-            )
-        except Exception:
-            gal_guess_flag = False
-            gal_guess_tmp = np.array([0., 0., 0., 0., 1, 100])
-
-        # Recenter jacobian if necessary
-        gal_jacob = ngmix.Jacobian(
-            row=(gals[0].shape[0] - 1) / 2 + gal_guess_tmp[0],
-            col=(gals[0].shape[1] - 1) / 2 + gal_guess_tmp[1],
-            wcs=jacob_list[n_e]
-        )
-
-        # Noise handling: this should be it's own function
-        # This code combines integrates flag information into the weights.
-        # Both THELI and Megapipe weights must be rescaled by the 
-        # inverse variance because ngmix expects inverse variance maps.  
-        if gal_guess_flag:
-            sig_noise = get_noise(
-                gals[n_e],
-                weight_map,
-                gal_guess_tmp,
-                pixel_scale,
-            )
-        else:
-            sig_noise = sigma_mad(gals[n_e])
-
-        noise_img = np.random.randn(*gals[n_e].shape) * sig_noise
-        noise_img_gal = np.random.randn(*gals[n_e].shape) * sig_noise
-
-        gal_masked = np.copy(gals[n_e])
-        if (len(np.where(weight_map == 0)[0]) != 0):
-            gal_masked[weight_map == 0] = noise_img_gal[weight_map == 0]
-
-        weight_map *= 1 / sig_noise ** 2
-
-        # gaussian fit to the original psf 
-        #(pre-counterfactual image operations)
-        w_tmp = np.sum(weight_map)
-
-        psf_res_gT['g_PSFo'] += psf_res['g'] * w_tmp
-        psf_res_gT['g_err_PSFo'] += np.array([
-            psf_res['pars_err'][2],
-            psf_res['pars_err'][3]
-        ]) * w_tmp
-        psf_res_gT['r50_PSFo'] += psf_res['pars'][4] * w_tmp
-        psf_res_gT['r50_err_PSFo'] += psf_res['pars_err'][4] * w_tmp
-        wsum += w_tmp
-
-        gal_obs = Observation(
-            gal_masked,
-            weight=weight_map,
-            jacobian=gal_jacob,
-            psf=psf_obs,
-            noise=noise_img
-        )
-
-        if gal_guess_flag:
-            gal_guess_tmp[:2] = 0
-            gal_guess.append(gal_guess_tmp)
-
-        gal_obs_list.append(gal_obs)
-        r50_guess_psf.append(psf_r50)
-        gal_guess_flag = True
-
-    if wsum == 0:
-        raise ZeroDivisionError('Sum of weights = 0, division by zero')
-
-    # Normalize PSF fit output
-    for key in psf_res_gT.keys():
-        psf_res_gT[key] /= wsum
-
-    # Gal guess handling
-    fail_get_guess = False
-    if len(gal_guess) == 0:
-        fail_get_guess = True
-        gal_pars = [0., 0., 0., 0., 1, 100]
-    else:
-        gal_pars = np.mean(gal_guess, 0)
-
-    psf_model = 'gauss'
-    gal_model = 'gauss'
-
-    # metacal specific parameters
-    metacal_pars = {
-        'types': ['noshear', '1p', '1m', '2p', '2m'],
-        'step': 0.01,
-        'psf': 'gauss',
-        'fixnoise': True,
-        'use_noise_image': True
-    }
-
-    r50guess = np.mean(r50_guess_psf)
-
-    # retry the fit twice
-    ntry = 2
-    ###### put in options for ngmix.metacal.metacal_bootstrap  obs,
-   # runner,
-   # psf_runner=None,
-   # ignore_failed_psf=True,
-   # rng=None,
-   # **metacal_kws
-#):)
-    obs_dict_mcal = ngmix.metacal.get_all_metacal(gal_obs_list, **metacal_pars)
-    res = {'mcal_flags': 0}
-
-    ntry = 5
-
-    for key in sorted(obs_dict_mcal):
-
-        fres = make_galsimfit(
-            obs_dict_mcal[key],
-            gal_model,
-            gal_pars,
-            prior=prior
-        )
-
-        res['mcal_flags'] |= fres['flags']
-        tres = {}
-
-        for name in fres.keys():
-            tres[name] = fres[name]
-        tres['flags'] = fres['flags']
-
-        wsum = 0
-        r50psf_sum = 0
-        gpsf_sum = np.zeros(2)
-        npsf = 0
-        for obs in obs_dict_mcal[key]:
-
-            if hasattr(obs, 'psf_nopix'):
-                try:
-                    psf_res = make_galsimfit(
-                        obs.psf_nopix,
-                        psf_model,
-                        np.array([0., 0., 0., 0., r50guess, 1.]),
-                        ntry=ntry
-                    )
-                except Exception:
-                    continue
-                g1, g2 = psf_res['g']
-                r50 = psf_res['pars'][4]
-            else:
-                try:
-                    psf_res = make_galsimfit(
-                        obs.psf,
-                        psf_model,
-                        np.array([0., 0., 0., 0., r50guess, 1.]),
-                    )
-                except Exception:
-                    continue
-                g1, g2 = psf_res['g']
-                r50 = psf_res['pars'][4]
-
-            # TODO we sometimes use other weights
-            twsum = obs.weight.sum()
-
-            wsum += twsum
-            gpsf_sum[0] += g1 * twsum
-            gpsf_sum[1] += g2 * twsum
-            r50psf_sum += r50 * twsum
-            npsf += 1
-
-        tres['gpsf'] = gpsf_sum / wsum
-        tres['r50psf'] = r50psf_sum / wsum
-
-        res[key] = tres
-
-    # result dictionary, keyed by the types in metacal_pars above
-    metacal_res = res
-
-    metacal_res.update(psf_res_gT)
-    metacal_res['moments_fail'] = fail_get_guess
-
-    return metacal_res