class Configuration(object):
    """
    Configuration class for redmapper.

    This class holds all the relevant configuration information, and had
    convenient methods for validating parameters as well as generating
    filenames, etc.

    """
    version = ConfigField(default=__version__, required=True)

    galfile = ConfigField(required=True)
    zredfile = ConfigField()
    halofile = ConfigField()
    randfile = ConfigField()
    catfile = ConfigField()
    specfile = ConfigField()
    specfile_train = ConfigField()
    outbase = ConfigField(required=True)
    parfile = ConfigField()
    bkgfile = ConfigField()
    bkgfile_color = ConfigField()
    zlambdafile = ConfigField()
    maskfile = ConfigField()
    depthfile = ConfigField()
    wcenfile = ConfigField()
    redgalfile = ConfigField()
    redgalmodelfile = ConfigField()
    seedfile = ConfigField()
    zmemfile = ConfigField()
    redmagicfile = ConfigField()

    calib_nproc = ConfigField(default=1, required=True)
    calib_run_nproc = ConfigField(default=1, required=True)
    calib_run_min_nside = ConfigField(default=1, required=True)

    runcat_percolation_masking = ConfigField(default=True, required=False)

    outpath = ConfigField(default='./', required=True)
    plotpath = ConfigField(default='', required=True)
    logpath = ConfigField(default='logs', required=True)

    border = ConfigField(default=0.0, required=True)
    hpix = ConfigField(default=[], required=True, isArray=True)
    nside = ConfigField(default=0, required=True)
    galfile_pixelized = ConfigField(required=True)
    use_tempfiles_to_conserve_memory = ConfigField(default=True, required=True)

    printlogging = ConfigField(default=True, required=True)

    nmag = ConfigField(required=True)
    area = ConfigField(required=True)
    limmag_catalog = ConfigField(required=True)
    limmag_ref = ConfigField(required=True)
    refmag = ConfigField(required=True)
    ref_ind = ConfigField(required=True)
    zeropoint = ConfigField(required=True)
    survey_mode = ConfigField(required=True)
    b = ConfigField(isArray=True)
    galfile_nside = ConfigField(required=True)
    bands = ConfigField(required=True)
    has_truth = ConfigField(default=False)

    area_finebin = ConfigField(default=0.001, required=True)
    area_coarsebin = ConfigField(default=0.005, required=True)
    area_nodesize = ConfigField(default=0.05, required=True)

    zrange = ConfigField(isArray=True, array_length=2, required=True)
    lval_reference = ConfigField(default=0.2, required=True)

    maskgalfile = ConfigField(default='maskgals.fit', required=True)
    mask_mode = ConfigField(default=0, required=True)
    max_maskfrac = ConfigField(default=0.2, required=True)
    covmask_nside_default = ConfigField(default=32, required=True)

    maskgal_ngals = ConfigField(default=6000, required=True)
    maskgal_nsamples = ConfigField(default=100, required=True)
    maskgal_rad_stepsize = ConfigField(default=0.1, required=True)
    maskgal_dmag_extra = ConfigField(default=0.3, required=True)
    maskgal_zred_err = ConfigField(default=0.02, required=True)

    dldr_gamma = ConfigField(default=0.6, required=True)
    rsig = ConfigField(default=0.05, required=True)
    chisq_max = ConfigField(default=20.0, required=True)
    npzbins = ConfigField(default=21, required=True)

    zred_nsamp = ConfigField(default=4, required=True)

    mstar_survey = ConfigField(default='sdss')
    mstar_band = ConfigField(default='i03')

    calib_niter = ConfigField(default=3)
    calib_zrange_cushion = ConfigField(default=0.05)

    calib_use_pcol = ConfigField(default=True)
    calib_smooth = ConfigField(default=0.003)
    calib_minlambda = ConfigField(default=5.0)
    calib_redgal_template = ConfigField()
    calib_pivotmag_nodesize = ConfigField(default=0.1)
    calib_color_nodesizes = ConfigField(isArray=True, default=np.array([0.05]))
    calib_slope_nodesizes = ConfigField(isArray=True, default=np.array([0.1]))
    calib_color_maxnodes = ConfigField(isArray=True, default=np.array([-1.0]))
    calib_covmat_maxnodes = ConfigField(isArray=True, default=np.array([-1.0]))
    calib_covmat_nodesize = ConfigField(default=0.15)
    calib_covmat_constant = ConfigField(default=0.9)
    calib_corr_nodesize = ConfigField(default=0.05)
    calib_corr_slope_nodesize = ConfigField(default=0.1)
    calib_corr_nocorrslope = ConfigField(default=True)
    calib_corr_pcut = ConfigField(default=0.9)
    calib_fit_err_ratio = ConfigField(default=False)
    calib_err_ratio_pivot = ConfigField(default=20.0)
    calib_spec_max_zerr = ConfigField(default=0.001, required=True)

    calib_color_nsig = ConfigField(default=1.5)
    calib_redspec_nsig = ConfigField(default=2.0)

    calib_colormem_r0 = ConfigField(default=0.5)
    calib_colormem_beta = ConfigField(default=0.0)
    calib_colormem_smooth = ConfigField(default=0.003)
    calib_colormem_minlambda = ConfigField(default=10.0)
    calib_colormem_zbounds = ConfigField(isArray=True, default=np.array([0.4]))
    calib_colormem_colormodes = ConfigField(isArray=True, default=np.array([1, 2]))
    calib_colormem_sigint = ConfigField(isArray=True, default=np.array([0.05, 0.03]))
    calib_pcut = ConfigField(default=0.3)
    calib_color_pcut = ConfigField(default=0.7)

    calib_zlambda_nodesize = ConfigField(default=0.04)
    calib_zlambda_slope_nodesize = ConfigField(default=0.1)
    calib_zlambda_minlambda = ConfigField(default=20.0)
    calib_zlambda_clean_nsig = ConfigField(default=5.0)
    calib_zlambda_correct_niter = ConfigField(default=3)

    calib_lumfunc_alpha = ConfigField(default=-1.0, required=True)

    zredc_binsize_fine = ConfigField(default=0.0001)
    zredc_binsize_coarse = ConfigField(default=0.005)

    bkg_chisqbinsize = ConfigField(default=0.5)
    bkg_refmagbinsize = ConfigField(default=0.2)
    bkg_zbinsize = ConfigField(default=0.02)
    bkg_zredbinsize = ConfigField(default=0.01)
    bkg_deepmode = ConfigField(default=False)
    calib_make_full_bkg = ConfigField(default=True)
    bkg_local_annuli = ConfigField(isArray=True, array_length=2,
                                   default=np.array([2.0, 3.0]))
    bkg_local_compute = ConfigField(default=False)
    bkg_local_use = ConfigField(default=False)

    zlambda_pivot = ConfigField(default=30.0, required=True)
    zlambda_binsize = ConfigField(default=0.002, required=True)
    zlambda_maxiter = ConfigField(default=20, required=True)
    zlambda_tol = ConfigField(default=0.0002, required=True)
    zlambda_topfrac = ConfigField(default=0.7, required=True)
    zlambda_epsilon = ConfigField(default=0.005, required=True)
    zlambda_parab_step = ConfigField(default=0.001, required=True)

    centerclass = ConfigField(default='CenteringBCG', required=True)
    wcen_rsoft = ConfigField(default=0.05, required=True)
    wcen_zred_chisq_max = ConfigField(default=100.0, required=True)
    wcen_minlambda = ConfigField(default=10.0, required=True)
    wcen_maxlambda = ConfigField(default=100.0, required=True)
    wcen_cal_zrange = ConfigField(isArray=True, default=np.array([0.0,1.0]))
    wcen_pivot = ConfigField(default=30.0, required=True)
    wcen_uselum = ConfigField(default=True, required=True)
    wcen_Delta0 = ConfigField(required=False, default=0.0)
    wcen_Delta1 = ConfigField(required=False, default=0.0)
    wcen_sigma_m = ConfigField(required=False, default=0.0)
    lnw_cen_sigma = ConfigField(required=False, default=-9999.0)
    lnw_cen_mean = ConfigField(required=False, default=-9999.0)
    lnw_sat_sigma = ConfigField(required=False, default=-9999.0)
    lnw_sat_mean = ConfigField(required=False, default=-9999.0)
    lnw_fg_sigma = ConfigField(required=False, default=-9999.0)
    lnw_fg_mean = ConfigField(required=False, default=-9999.0)
    phi1_mmstar_m = ConfigField(required=False, default=-9999.0)
    phi1_mmstar_slope = ConfigField(required=False, default=-9999.0)
    phi1_msig_m = ConfigField(required=False, default=-9999.0)
    phi1_msig_slope = ConfigField(required=False, default=-9999.0)

    firstpass_r0 = ConfigField(default=0.5, required=True)
    firstpass_beta = ConfigField(default=0.0, required=True)
    firstpass_niter = ConfigField(default=2, required=True)
    firstpass_minlambda = ConfigField(default=3.0, required=True)
    firstpass_centerclass = ConfigField(default='CenteringBCG', required=True)

    likelihoods_r0 = ConfigField(default=1.0, required=True)
    likelihoods_beta = ConfigField(default=0.2, required=True)
    likelihoods_use_zred = ConfigField(default=True, required=True)
    likelihoods_minlambda = ConfigField(default=3.0, required=True)

    percolation_r0 = ConfigField(default=1.0, required=True)
    percolation_beta = ConfigField(default=0.2, required=True)
    percolation_rmask_0 = ConfigField(default=1.5, required=True)
    percolation_rmask_beta = ConfigField(default=0.2, required=True)
    percolation_rmask_gamma = ConfigField(default=0.0, required=True)
    percolation_rmask_zpivot = ConfigField(default=0.3, required=True)
    percolation_lmask = ConfigField()
    percolation_niter = ConfigField(default=2, required=True)
    percolation_minlambda = ConfigField(default=3.0, required=True)
    percolation_pbcg_cut = ConfigField(default=0.5, required=True)
    percolation_maxcen = ConfigField(default=5, required=True)
    percolation_memradius = ConfigField()
    percolation_memlum = ConfigField()

    zscan_r0 = ConfigField(default=0.5, required=True)
    zscan_beta = ConfigField(default=0.0, required=True)
    zscan_zstep = ConfigField(default=0.005, required=True)
    zscan_minlambda = ConfigField(default=3.0, required=True)

    vlim_lstar = ConfigField(default=0.2, required=False)
    vlim_depthfiles = ConfigField(default=[], required=False, isList=True)
    vlim_bands = ConfigField(default=[], required=False, isList=True)
    vlim_nsigs = ConfigField(default=[], required=False, isArray=True)

    consolidate_lambda_cuts = ConfigField(default=[5.0, 20.0], required=False, isArray=True)
    consolidate_vlim_lstars = ConfigField(default=[0.2, 5.0], required=False, isList=True)
    select_scaleval = ConfigField(default=False, required=True)

    redmagic_calib_nodesize = ConfigField(default=0.05, required=True)
    redmagic_calib_corr_nodesize = ConfigField(default=0.05, required=True)
    redmagic_calib_buffer = ConfigField(default=0.05, required=True)
    redmagic_calib_zbinsize = ConfigField(default=0.02, required=True)
    redmagic_calib_chisqcut = ConfigField(default=20.0, required=True)
    redmagic_zrange = ConfigField(default=[], required=False, isArray=True)
    redmagic_calib_fractrain = ConfigField(default=0.5, required=True)
    redmagic_calib_redshift_buffer = ConfigField(default=0.05, required=True)
    redmagic_maxlum = ConfigField(default=100.0, required=True)
    redmagic_mock_truthspec = ConfigField(default=False, required=True)
    redmagic_run_afterburner = ConfigField(default=True, required=True)
    redmagic_apply_afterburner_zsamp = ConfigField(default=True, required=True)
    redmagic_n0s = ConfigField(default=[], required=True, isArray=True)
    redmagic_etas = ConfigField(default=[], required=True, isArray=True)
    redmagic_names = ConfigField(default=[], required=True, isList=True)
    redmagic_zmaxes = ConfigField(default=[], required=True, isArray=True)
    redmagic_constchis = ConfigField(default=[], required=True, isArray=True)
    redmagic_use_constchi = ConfigField(default=False, required=True)

    def __init__(self, configfile, outpath=None):
        """
        Instantiate a Configuration object

        Parameters
        ----------
        configfile: `str`
           Configuration yaml filename
        outpath: `str`, optional
           Path for output.  Default is None, which uses yaml value.

        Raises
        ------
        RuntimeError:
           When invalid configs are used
        """
        self._file_logging_started = False

        self._reset_vars()

        # First, read in the yaml file
        confdict = read_yaml(configfile)

        self.configpath = os.path.dirname(os.path.abspath(configfile))
        self.configfile = os.path.basename(configfile)

        # And now set the config variables
        self._set_vars_from_dict(confdict)

        if outpath is not None:
            self.outpath = outpath

        # Get the logger
        if len(self.hpix) == 1:
            logname = f'redmapper-{self.hpix[0]}'
        else:
            logname = 'redmapper'
        logging.basicConfig(level=logging.INFO)
        self.logger = logging.getLogger(logname)

        # validate the galfile and refmag
        type(self).__dict__['galfile'].validate('galfile')
        type(self).__dict__['refmag'].validate('refmag')

        # get galaxy file stats
        gal_stats = self._galfile_stats()
        self.galfile_area = gal_stats['area']
        if (self.area is not None):
            if self.depthfile is not None:
                self.logger.info("WARNING: You should not need to set area in the config file when you have a depth map.")
            if (np.abs(self.area - gal_stats['area']) > 1e-3):
                self.logger.info("Config area is not equal to galaxy file area.  Using config area.")
                gal_stats.pop('area')
        else:
            if self.depthfile is None and self.nside > 0 and len(self.hpix) > 0:
                raise RuntimeError("You must set a config file area if no depthfile is present and you are running a sub-region")
        self._set_vars_from_dict(gal_stats, check_none=True)

        if self.limmag_catalog is None:
            self.limmag_catalog = self.limmag_ref

        # Get wcen numbers if available
        self.set_wcen_vals()

        # Set some defaults
        if self.specfile_train is None:
            self.specfile_train = self.specfile

        # Record the cluster dtype for convenience
        self.cluster_dtype = copy.copy(cluster_dtype_base)
        self.cluster_dtype.extend([('MAG', 'f4', self.nmag),
                                   ('MAG_ERR', 'f4', self.nmag),
                                   ('PZBINS', 'f4', self.npzbins),
                                   ('PZ', 'f4', self.npzbins),
                                   ('RA_CENT', 'f8', self.percolation_maxcen),
                                   ('DEC_CENT', 'f8', self.percolation_maxcen),
                                   ('ID_CENT', 'i8', self.percolation_maxcen),
                                   ('LAMBDA_CENT', 'f4', self.percolation_maxcen),
                                   ('ZLAMBDA_CENT', 'f4', self.percolation_maxcen),
                                   ('P_CEN', 'f4', self.percolation_maxcen),
                                   ('Q_CEN', 'f4', self.percolation_maxcen),
                                   ('P_FG', 'f4', self.percolation_maxcen),
                                   ('Q_MISS', 'f4'),
                                   ('P_SAT', 'f4', self.percolation_maxcen),
                                   ('P_C', 'f4', self.percolation_maxcen)])
        self.member_dtype = copy.copy(member_dtype_base)
        self.member_dtype.extend([('MAG', 'f4', self.nmag),
                                  ('MAG_ERR', 'f4', self.nmag)])
        # also need pz stuff, etc, etc.  Will need to deal with defaults

        # Calibration size checking
        self._set_lengths(['calib_colormem_colormodes', 'calib_colormem_sigint'],
                          len(self.calib_colormem_zbounds) + 1)
        self._set_lengths(['calib_color_nodesizes', 'calib_slope_nodesizes',
                           'calib_color_maxnodes', 'calib_covmat_maxnodes'],
                           self.nmag - 1)
                           #'calib_color_order'], self.nmag - 1)

        # redmagic size checks
        self._set_lengths(['redmagic_n0s', 'redmagic_etas', 'redmagic_names',
                           'redmagic_zmaxes'], len(self.redmagic_names))

        # Vlim size checks
        self._set_lengths(['vlim_bands', 'vlim_nsigs'],
                          len(self.vlim_depthfiles))

        # will want to set defaults here...
        self.cosmo = Cosmo()

        # Redmagic stuff
        if len(self.redmagic_zrange) == 0 or self.redmagic_zrange[0] < 0.0 or self.redmagic_zrange[1] < 0.0:
            self.redmagic_zrange = self.zrange[:]

        self._set_lengths(['redmagic_zrange'], 2)

        # Finally, we can validate...
        self.validate()

        # Checks ...
        if self.maskfile is not None and self.mask_mode == 0:
            raise ValueError("A maskfile is set, but mask_mode is 0 (no mask).  Assuming this is not intended.")

        for vlim_band in self.vlim_bands:
            if vlim_band not in self.bands:
                raise ValueError("vlim_band %s not in list of bands!" % (vlim_band))

        if self.bkg_local_annuli[1] <= self.bkg_local_annuli[0]:
            raise ValueError("bkg_local_annuli[1] must be > bkg_local_annuli[0]")

        # Now set the duplicatable config parameters...
        self.d = DuplicatableConfig(self)

        # Finally, once everything is here, we can make paths
        if not os.path.exists(self.outpath):
            os.makedirs(self.outpath, exist_ok=True)
        if not os.path.exists(os.path.join(self.outpath, self.plotpath)):
            os.makedirs(os.path.join(self.outpath, self.plotpath), exist_ok=True)

    def validate(self):
        """
        Validate the configuration.

        Raises
        ------
        ValueError:
           If any config field is not legal, ValueError is raised.
        """

        for var in type(self).__dict__:
            try:
                type(self).__dict__[var].validate(var)
            except AttributeError:
                pass

    def copy(self):
        """
        Return a copy of the configuration struct
        """
        return copy.copy(self)

    def start_file_logging(self, filename=None):
        """
        Start logging to a file.

        Parameters
        ----------
        filename : `str`, optional
            Optional filename, else will be determined from outbase and hpix.
        """
        if self._file_logging_started:
            return

        if self.printlogging:
            self.logger.info("Logging is set to be to console only.")
            return

        if not os.path.exists(os.path.join(self.outpath, self.logpath)):
            os.makedirs(os.path.join(self.outpath, self.logpath), exist_ok=True)

        if filename is None:
            logfilename = os.path.join(self.outpath, self.logpath,
                                       f'redmapper_{self.d.outbase}_{self.d.hpix[0]:04}.log')
        else:
            logfilename = os.path.join(self.outpath, self.logpath, os.path.basename(filename))

        handler = logging.FileHandler(filename=logfilename)
        self.logger.addHandler(handler)

        self._file_logging_started = True

    def stop_file_logging(self):
        """Stop logging to the file.
        """
        handlers = self.logger.handlers[:]
        for handler in handlers:
            handler.close()
            self.logger.removeHandler(handler)

    def _reset_vars(self):
        """
        Internal method to reset variables to defaults
        """
        for var in type(self).__dict__:
            try:
                type(self).__dict__[var].reset()
            except AttributeError:
                pass

    def _set_vars_from_dict(self, d, check_none=False):
        """
        Internal method to set config variables from a dictionary.

        Parameters
        ----------
        d: `dict`
           Dictionary of key/value pairs to set
        check_none: `bool`, optional
           If true, don't set any variables that have a value of None in the dict.
           Default is False
        """
        for key in d:
            if check_none and d[key] is None:
                continue
            if key not in type(self).__dict__:
                raise AttributeError("Unknown config variable: %s" % (key))
            try:
                setattr(self, key, d[key])
            except TypeError:
                raise TypeError("Error with type of variable %s" % (key))

    def _set_lengths(self, l, length):
        """
        Internal method to set the validation length for a list of field names

        Parameters
        ----------
        l: `list`
           Field names to set
        length: `int`
           Validation length
        """
        for arr in l:
            if arr not in type(self).__dict__:
                raise AttributeError("Unknown config variable: %s" % (arr))
            type(self).__dict__[arr].set_length(length)


    def _galfile_stats(self):
        """
        Internal method to get statistics from the galfile

        Returns
        -------
        gal_stats: `dict`
           Dictionary with galaxy file stats

        Raises
        ------
        ValueError:
           Raise error if galaxy file is incorrect format.
        """

        hdr = fitsio.read_header(self.galfile, ext=1)
        pixelated = hdr.get("PIXELS", 0)
        fitsformat = hdr.get("FITS", 0)

        if not fitsformat:
            raise ValueError("Input galfile must describe fits files.")

        gal_stats = {}
        if not pixelated:
            # statistics are from the header
            gal_stats['galfile_pixelized'] = False

            hdrmode = hdr.get("MODE", "").rstrip()

            if hdrmode == 'SDSS':
                gal_stats['survey_mode'] = 0
            elif hdrmode == 'DES':
                gal_stats['survey_mode'] = 1
            elif hdrmode == 'LSST':
                gal_stats['survey_mode'] = 2
            else:
                raise ValueError("Input galaxy file with unknown mode: %s" % (hdrmode))

            ## FIXME: check that these are all in the header
            gal_stats['area'] = hdr.get('AREA', -100.0)
            gal_stats['limmag_ref'] = hdr.get('LIM_REF')
            gal_stats['nmag'] = hdr.get('NMAG')
            gal_stats['zeropoint'] = hdr.get('ZP')
            gal_stats['ref_ind'] = hdr.get(self.refmag.upper()+'_IND')
            gal_stats['b'] = None
            gal_stats['galfile_nside'] = 0
            gal_stats['bands'] = [None]*gal_stats['nmag']

            # figure out the bands...
            for name in hdr:
                m = re.search('(.*)_IND', name)
                if m is None:
                    continue
                if m.groups()[0] == 'REF':
                    continue
                band = m.groups()[0].lower()
                gal_stats['bands'][hdr[name]] = band

        else:
            # statistics are from the master table file
            gal_stats['galfile_pixelized'] = True

            master=fitsio.read(self.galfile, ext=1, upper=True)

            try:
                # Support for old fits reading
                mode = master['MODE'][0].rstrip().decode()
            except AttributeError:
                mode = master['MODE'][0].rstrip()
            if (mode == 'SDSS'):
                gal_stats['survey_mode'] = 0
            elif (mode == 'DES'):
                gal_stats['survey_mode'] = 1
            elif (mode == 'LSST'):
                gal_stats['survey_mode'] = 2
            else:
                raise ValueError("Input galaxy file with unknown mode: %s" % (mode))

            gal_stats['area'] = master['AREA'][0]
            gal_stats['limmag_ref'] = master['LIM_REF'][0]
            gal_stats['nmag'] = master['NMAG'][0]
            if ('B' in master.dtype.names):
                gal_stats['b'] = master['B'][0]
            else:
                gal_stats['b'] = None
            gal_stats['zeropoint'] = master['ZEROPOINT'][0]
            gal_stats['ref_ind'] = master[self.refmag.upper()+'_IND'][0]
            gal_stats['galfile_nside'] = master['NSIDE'][0]
            gal_stats['bands'] = [None]*gal_stats['nmag']

            # Figure out the bands...
            for name in master.dtype.names:
                m = re.search('(.*)_IND', name)
                if m is None:
                    continue
                if m.groups()[0] == 'REF':
                    continue
                band = m.groups()[0].lower()
                gal_stats['bands'][master[name][0]] = band

        if any(x is None for x in gal_stats['bands']):
            # Remove this, and use config values
            gal_stats.pop('bands', None)

        return gal_stats

    def set_wcen_vals(self):
        """
        Set wcen centering values.  These will be loaded from self.wcenfile if
        available.
        """

        wcen_vals = self._wcen_vals()
        if wcen_vals is not None:
            self._set_vars_from_dict(wcen_vals)

    def _wcen_vals(self):
        """
        Load in wcen values from self.wcenfile

        Returns
        -------
        vals: `dict`
           Dictionary of wcen centering parameters
        """

        if self.wcenfile is None or not os.path.isfile(self.wcenfile):
            # We don't have wcen info to load
            return None

        wcen = fitsio.read(self.wcenfile, ext=1, lower=True)

        vals = {'wcen_Delta0': wcen[0]['delta0'],
                'wcen_Delta1': wcen[0]['delta1'],
                'wcen_sigma_m': wcen[0]['sigma_m'],
                'wcen_pivot': wcen[0]['pivot'],
                'lnw_fg_mean': wcen[0]['lnw_fg_mean'],
                'lnw_fg_sigma': wcen[0]['lnw_fg_sigma'],
                'lnw_sat_mean': wcen[0]['lnw_sat_mean'],
                'lnw_sat_sigma': wcen[0]['lnw_sat_sigma'],
                'lnw_cen_mean': wcen[0]['lnw_cen_mean'],
                'lnw_cen_sigma': wcen[0]['lnw_cen_sigma']}

        # New wcen files also record the phi1 information
        if 'phi1_mmstar_m' in wcen[0].dtype.names:
            vals['phi1_mmstar_m'] = wcen[0]['phi1_mmstar_m']
            vals['phi1_mmstar_slope'] = wcen[0]['phi1_mmstar_m']
            vals['phi1_msig_m'] = wcen[0]['phi1_msig_m']
            vals['phi1_msig_slope'] = wcen[0]['phi1_msig_slope']

        return vals

    @property
    def zrange_cushioned(self):
        """Return zrange with additional cushion."""
        zrange_cushioned = self.zrange.copy()
        zrange_cushioned[0] = np.clip(zrange_cushioned[0] - self.calib_zrange_cushion, 0.05, None)
        zrange_cushioned[1] += self.calib_zrange_cushion
        return zrange_cushioned

    def redmapper_filename(self, redmapper_name, paths=None, filetype='fit',
                           withversion=False, outbase=None):
        """
        Generate a redmapper filename with all the appropriate infixes.

        Parameters
        ----------
        redmapper_name: `str`
           String describing the redmapper file type
        paths: `list` or `tuple`, optional
           List or tuple of path strings to join for filename.
           Default is None, just use self.outpath
        filetype: `str`, optional
           File extension.  Default is 'fit`
        withversion : `bool`, optional
           Add in the redmapper version string
        outbase : `str`, optional
           Override self.d.outbase

        Returns
        -------
        filename: `str`
           Properly formatted full-path redmapper filename
        """
        if outbase is None:
            outbase = self.d.outbase

        if withversion:
            outbase += '_redmapper_v%s' % (self.version)

        if paths is None:
            return os.path.join(self.outpath,
                                '%s_%s.%s' % (outbase, redmapper_name, filetype))
        else:
            if type(paths) is not list and type(paths) is not tuple:
                raise ValueError("paths must be a list or tuple")
            pars = [self.outpath]
            pars.extend(paths)
            pars.append('%s_%s.%s' % (outbase, redmapper_name, filetype))
            return os.path.join(*pars)

    def check_files(self, check_zredfile=False, check_bkgfile=False, check_bkgfile_components=False,
                    check_parfile=False, check_zlambdafile=False, check_randfile=False):
        """
        Check that all calibration files are available for a cluster finder run.

        Parameters
        ----------
        check_zredfile: `bool`, optional
           Check that the zred file is available.  Default is False.
        check_bkgfile: `bool`, optional
           Check that the background file is available.  Default is False.
        check_bkgfile_components: `bool`, optional
           Check that the background file zred/chisq components are available.
           Default is False.
        check_parfile: `bool`, optional
           Check that the red sequence parameter file is available.  Default is False.
        check_zlambdafile: `bool`, optional
           Check that the zlambda calibration file is available.  Default is False.
        check_randfile: `bool`, optional
           Check that the random file is available.  Default is False.

        Raises
        ------
        ValueError:
           Raise ValueError if any of the check files is missing.
        """
        if check_zredfile:
            if not os.path.isfile(self.zredfile):
                raise ValueError("zredfile %s not found." % (self.zredfile))

        if check_bkgfile:
            if not os.path.isfile(self.bkgfile):
                raise ValueError("bkgfile %s not found." % (self.bkgfile))

        if check_bkgfile_components:
            with fitsio.FITS(self.bkgfile) as fits:
                if 'CHISQBKG' not in [ext.get_extname() for ext in fits[1: ]]:
                    raise ValueError("bkgfile %s does not have CHISQBKG extension." % (self.bkgfile))
                if 'ZREDBKG' not in [ext.get_extname() for ext in fits[1: ]]:
                    raise ValueError("bkgfile %s does not have ZREDBKG extension." % (self.bkgfile))

        if check_parfile:
            if not os.path.isfile(self.parfile):
                raise ValueError("parfile %s not found." % (self.parfile))

        if check_zlambdafile:
            if not os.path.isfile(self.zlambdafile):
                raise ValueError("zlambdafile %s not found." % (self.zlambdafile))

        if check_randfile:
            if not os.path.isfile(self.randfile):
                raise ValueError("randfile %s not found." % (self.randfile))

    def output_yaml(self, filename):
        """
        Output config into a yaml file.

        Parameters
        ----------
        filename: `str`
           Output yaml filename
        """
        out_dict = {}
        for key in type(self).__dict__:
            if isinstance(type(self).__dict__[key], ConfigField):
                if type(self).__dict__[key]._isArray:
                    # Convert all elements to python types and make a list
                    out_dict[key] = np.ndarray.tolist(type(self).__dict__[key]._value)
                else:
                    # Try to use numpy to convert to scalar; if it doesn't work then
                    # it's not numpy and we can use the value directly
                    try:
                        out_dict[key] = type(self).__dict__[key]._value.item()
                    except (ValueError, AttributeError, TypeError):
                        out_dict[key] = type(self).__dict__[key]._value

        with open(filename, 'w') as f:
            yaml.dump(out_dict, stream=f)

    def compute_border(self):
        """
        Compute the border radii based on the largest expected cluster at the lowest
        redshift.

        Returns
        -------
        rad: `float`
           Border radius for overlapping tiles for parallel runs.
        """

        maxdist = 1.05 * self.percolation_rmask_0 * (300. / 100.)**self.percolation_rmask_beta
        radius = maxdist / (np.radians(1.) * self.cosmo.Da(0, self.zrange[0]))

        return 3.0 * radius class MedZFitter(object):
    """
    Class to fit a spline to the median value as a function of redshift.

    This can be the color or magnitude or any quantity.
    """
    def __init__(self, z_nodes, redshifts, values):
        """
        Instantiate a MedZFitter object

        Parameters
        ----------
        z_nodes: `np.array`
           Float array for redshift nodes
        redshifts: `np.array`
           Float array of input redshifts to fit
        values: `np.array`
           Float array of color values to fit
        """
        self._z_nodes = z_nodes.astype(np.float64)
        self._redshifts = redshifts.astype(np.float64)
        self._values = values.astype(np.float64)

    def fit(self, p0, min_val=-np.inf, max_val=np.inf):
        """
        Perform a spline fit to the median value as a function of redshift.

        Parameters
        ----------
        p0: `np.array`
           Initial fit parameters, with same number of elements as nodes

        Returns
        -------
        pars: `np.array`
           Spline node parameters
        """
        # add bounds if we need it...
        bounds = []
        for i in range(len(p0)):
            bounds.append([min_val, max_val])

        res = scipy.optimize.minimize(self,
                                      p0,
                                      method='L-BFGS-B',
                                      bounds=bounds,
                                      jac=False,
                                      options={'maxfun': 2000,
                                               'maxiter': 2000,
                                               'maxcor': 20,
                                               'eps': 1e-5,
                                               'gtol': 1e-8},
                                      callback=None)
        pars = res.x

        return pars

    def __call__(self, pars):
        """
        Compute the median cost function for f(pars)

        Parameters
        ----------
        pars: `np.array`
           Fit parameters, with same number of elements as nodes

        Returns
        -------
        t: `float`
           Median cost
        """
        spl = CubicSpline(self._z_nodes, pars)
        m = spl(self._redshifts)

        absdev = np.abs(self._values - m)
        t = np.sum(absdev.astype(np.float64))

        return t class RedSequenceFitter(object):
    """
    Class to fit a spline to the red sequence as a function of redshift.

    This class computes the mean value of the color, the red sequence slope,
    and the intrinsic scatter as a function of redshift.

    The node positions can be different for the mean color, the slope, and the
    intrinsic scatter.
    """
    def __init__(self, mean_nodes,
                 redshifts, colors, mag_errs, dmags=None, trunc=None,
                 slope_nodes=None, scatter_nodes=None, lupcorrs=None,
                 probs=None, bkgs=None, scatter_max=None, use_scatter_prior=False,
                 dmags_err_ratio=None):
        """
        Instantiate a RedSequenceFitter

        Parameters
        ----------
        mean_nodes: `np.array`
           Float array for mean color redshift nodes
        redshifts: `np.array`
           Float array of input redshifts for fit
        colors: `np.array`
           Float array of input colors to fit
        mag_errs: `np.array`
           Float array of input mag errors to fit [size, 2]
        dmags: `np.array`, optional
           Float array of delta-mag (refmag - pivot) (location on x axis for slope).
           Must have same length as colors if used.  Default is None (don't fit slope).
        trunc: `np.array`, optional
           Float array of delta-color used to truncate input colors.
           This outlier rejection is corrected for in the Gaussian likelihood.
           Must have same length as colors if used.  Default is None (no truncation).
        slope_nodes: `np.array`, optional
           Float array of red sequence slope node locations.
           Default is None (use mean_nodes for slope).
        scatter_nodes: `np.array`, optional
           Float array of red sequence scatter node locations.
           Default is None (use mean_nodes for scatter).
        lupcorrs: `np.array`, optional
           Float array of model corrections based on use of luptitudes.
           Must have same length as colors if used.
           Default is None (no luptitude corrections required).
        probs: `np.array`, optional
           Float array of membership probabilities.
           Must have same length as colors if used; must also supply "bkgs" if used.
           Default is None (assume all galaxies have p=1)
        bkgs: `np.array`, optional
           Float array of background likelihoods
           Must have same length as colors if used; must also supply "probs" if used.
           Default is None (assume all galaxies have bkg=0)
        scatter_max: `float`, optional
           Maximum intrinsic scatter allowed for any node.  Default is None (no max).
        use_scatter_prior: `bool`, optional
           Use Jeffry's prior on intrinsic scatter.  Default is False.
        dmags_err_ratio : `np.ndarray`, optional
            Delta-mag for error ratio computation.
        """

        self._use_scatter_prior = use_scatter_prior

        self._mean_nodes = np.atleast_1d(mean_nodes).astype(np.float64)
        if slope_nodes is None:
            self._slope_nodes = self._mean_nodes
        else:
            self._slope_nodes = np.atleast_1d(slope_nodes).astype(np.float64)
        if scatter_nodes is None:
            self._scatter_nodes = self._mean_nodes
        else:
            self._scatter_nodes = np.atleast_1d(scatter_nodes).astype(np.float64)

        self._redshifts = np.atleast_1d(redshifts).astype(np.float64)
        self._colors = np.atleast_1d(colors).astype(np.float64)
        self._mag_err2s = np.atleast_2d(mag_errs).astype(np.float64)**2.

        self._n_mean_nodes = self._mean_nodes.size
        self._n_slope_nodes = self._slope_nodes.size
        self._n_scatter_nodes = self._scatter_nodes.size

        if self._redshifts.size != self._colors.size:
            raise ValueError("Number of redshifts must be equal to colors")
        if self._redshifts.size != self._mag_err2s.shape[0]:
            raise ValueError("Number of redshifts must be equal to mag_errs.shape[1]")
        if self._mag_err2s.shape[1] != 2:
            raise ValueError("Mag_errs must by 2xNgals")

        if trunc is not None:
            self._trunc = np.atleast_1d(trunc).astype(np.float64)
            if self._redshifts.size != self._trunc.size:
                raise ValueError("Number of redshifts must be equal to truncs")
        else:
            self._trunc = None

        if dmags is not None:
            self._dmags = np.atleast_1d(dmags).astype(np.float64)
            if self._redshifts.size != self._dmags.size:
                raise ValueError("Number of redshifts must be equal to dmags")
            self._has_dmags = True
        else:
            self._dmags = np.zeros(self._redshifts.size).astype(np.float64)
            self._has_dmags = False

        if dmags_err_ratio is not None:
            self._dmags_err_ratio = np.atleast_1d(dmags_err_ratio).astype(np.float64)
            if self._redshifts.size != self._dmags_err_ratio.size:
                raise ValueError("Number of redshifts must be equal to dmags_err_ratio")
        else:
            self._dmags_err_ratio = np.zeros(self._redshifts.size)

        if lupcorrs is not None:
            self._lupcorrs = np.atleast_1d(lupcorrs).astype(np.float64)
            if self._redshifts.size != self._lupcorrs.size:
                raise ValueError("Number of redshifts must be equal to lupcorrs")
            self._has_lupcorrs = True
        else:
            self._lupcorrs = np.zeros(self._redshifts.size)
            self._has_lupcorrs = False

        if probs is not None:
            self._probs = np.atleast_1d(probs).astype(np.float64)
            if self._redshifts.size != self._probs.size:
                raise ValueError("Number of redshifts must be equal to probs")
            self._has_probs = True
        else:
            self._has_probs = False

        if bkgs is not None:
            self._bkgs = np.atleast_1d(bkgs).astype(np.float64)
            if self._redshifts.size != self._bkgs.size:
                raise ValueError("Number of redshifts must be equal to bkgs")
            self._has_bkgs = True
        else:
            self._has_bkgs = False

        if scatter_max is not None:
            self._scatter_max = np.atleast_1d(scatter_max).astype(np.float64)
            if self._scatter_max.size != self._n_scatter_nodes:
                raise ValueError("Number of scatter_max must be equal to scatter_nodes")
            self._has_scatter_max = True
        else:
            self._has_scatter_max = False

        if self._has_probs and not self._has_bkgs:
            raise ValueError("If you supply probs you must also supply bkgs")

    def fit(self, p0_mean, p0_slope, p0_scatter,
            fit_mean=False, fit_slope=False, fit_scatter=False,
            min_scatter=0.001, err_ratio_pars=None, fit_err_ratio_ind=[0, 1]):
        """
        Fit the red sequence mean, slope, and intrinsic scatter.

        Parameters
        ----------
        p0_mean: `np.array`
           Float array of initial fit parameters for mean relation.
           Must have same number of elements as mean_nodes
        p0_slope: `np.array`
           Float array of initial fit parameters for slope relation.
           Must have same number of elements as slope_nodes
        p0_scatter: `np.array`
           Float array of initial fit parameters for scatter relation.
           Must have same number of elements as scatter_nodes
        fit_mean: `bool`, optional
           Fit the mean relation? (else fix to p0_mean).  Default is False.
        fit_slope: `bool`, optional
           Fit the slope? (else fix to p0_slope).  Default is False.
        fit_scatter: `bool`, optional
           Fit the scatter? (else fix to p0_scatter).  Default is False.
        min_scatter: `float`, optional
           Minimum intrinsic scatter.  Default is 0.001.
        err_ratio_pars : `float`, optional
            Error ratio parameters.  Will be fit if fit_scatter=True and is
            not None.
        fit_err_ratio_ind : array-like, optional
            Magnitude indices to apply error ratio to fit.

        Returns
        -------
        pars_mean: `np.array`
           Mean parameters (if fit_mean is True)
        pars_slope: `np.array`
           Slope parameters (if fit_slope is True)
        pars_scatter: `np.array`
           Scatter parameters (if fit_scatter is True)
        """
        self._fit_mean = fit_mean
        self._fit_slope = fit_slope
        self._fit_scatter = fit_scatter
        self._min_scatter = min_scatter
        if err_ratio_pars is not None:
            self._has_err_ratios = True
            self._err_ratio_pars = err_ratio_pars
            self._fit_err_ratio_ind = fit_err_ratio_ind
        else:
            # Do not fit, use 1.0
            self._has_err_ratios = False
            self._err_ratio_pars = [1.0, 0.0]
            self._fit_err_ratio_ind = []

        if not self._has_dmags and self._fit_slope:
            raise ValueError("Can only do fit_slope if dmags were supplied")

        ctr = 0
        p0 = np.array([])
        bounds = []
        if self._fit_mean:
            self._mean_index = 0
            ctr += self._n_mean_nodes
            p0 = np.append(p0, p0_mean)
            for i in range(self._n_mean_nodes):
                bounds.append([-np.inf, np.inf])
        if self._fit_slope:
            self._slope_index = ctr
            ctr += self._n_slope_nodes
            p0 = np.append(p0, p0_slope)
            for i in range(self._n_slope_nodes):
                bounds.append([-np.inf, np.inf])
        if self._fit_scatter:
            self._scatter_index = ctr
            ctr += self._n_scatter_nodes
            p0 = np.append(p0, p0_scatter)
            for i in range(self._n_scatter_nodes):
                if self._has_scatter_max:
                    bounds.append([self._min_scatter, self._scatter_max[i]])
                else:
                    bounds.append([self._min_scatter, np.inf])
            if self._has_err_ratios:
                p0 = np.append(p0, self._err_ratio_pars[0])
                bounds.append([0.5, 10.0])
                p0 = np.append(p0, self._err_ratio_pars[1])
                bounds.append([-5.0, 5.0])

        if ctr == 0:
            raise ValueError("Must select at least one of fit_mean, fit_slope, fit_scatter")

        # Precompute...
        if not self._fit_mean:
            spl = CubicSpline(self._mean_nodes, p0_mean)
            self._gmean = spl(self._redshifts)
        if not self._fit_slope:
            spl = CubicSpline(self._slope_nodes, p0_slope)
            self._gslope = spl(self._redshifts)
        if not self._fit_scatter:
            spl = CubicSpline(self._scatter_nodes, p0_scatter)
            err_ratios = self._err_ratio_pars[0] + self._err_ratio_pars[1]*self._dmags_err_ratio
            if 0 in self._fit_err_ratio_ind:
                e2 = (err_ratios**2.)*self._mag_err2s[:, 0]
            else:
                e2 = self._mag_err2s[:, 0]
            if 1 in self._fit_err_ratio_ind:
                e2 += (err_ratios**2.)*self._mag_err2s[:, 1]
            else:
                e2 += self._mag_err2s[:, 1]
            self._gsig = np.sqrt(np.clip(spl(self._redshifts), self._min_scatter, None)**2. + e2)

        if not self._fit_scatter and self._trunc is not None:
            self._phi_bma = special.erf((self._trunc / self._gsig) / np.sqrt(2.))

        res = scipy.optimize.minimize(self,
                                      p0,
                                      method='L-BFGS-B',
                                      bounds=bounds,
                                      jac=False,
                                      options={'maxfun': 2000,
                                               'maxiter': 2000,
                                               'maxcor': 20,
                                               'eps': 1e-5,
                                               'gtol': 1e-8},
                                      callback=None)
        pars = res.x

        retval = []
        if self._fit_mean:
            retval.append(pars[self._mean_index: self._mean_index + self._n_mean_nodes])
        if self._fit_slope:
            retval.append(pars[self._slope_index: self._slope_index + self._n_slope_nodes])
        if self._fit_scatter:
            if self._has_err_ratios:
                retval.append(pars[self._scatter_index: self._scatter_index + self._n_scatter_nodes + 2])
            else:
                retval.append(pars[self._scatter_index: self._scatter_index + self._n_scatter_nodes])

        return retval

    def __call__(self, pars):
        """
        Compute the red sequence log-likelihood (negative for minimization).

        Parameters
        ----------
        pars: `np.array`
           Concatenated array of all the fit parameters

        Returns
        -------
        t: `float`
           Total negative log-likelihood.
        """
        if self._fit_mean:
            # We are fitting the mean
            spl = CubicSpline(self._mean_nodes, pars[self._mean_index: self._mean_index + self._n_mean_nodes])
            gmean = spl(self._redshifts)
        else:
            gmean = self._gmean

        if self._fit_slope:
            # We are fitting the slope
            spl = CubicSpline(self._slope_nodes, pars[self._slope_index: self._slope_index + self._n_slope_nodes])
            gslope = spl(self._redshifts)
        else:
            gslope = self._gslope

        if self._fit_scatter:
            # We are fitting the scatter
            spl = CubicSpline(self._scatter_nodes, pars[self._scatter_index: self._scatter_index + self._n_scatter_nodes])
            if self._has_err_ratios:
                # Always the last one
                err_ratio_pars = pars[-2: ]
            else:
                err_ratio_pars = [1.0, 0.0]
            err_ratios = err_ratio_pars[0] + err_ratio_pars[1]*self._dmags_err_ratio
            if 0 in self._fit_err_ratio_ind:
                e2 = (err_ratios**2.)*self._mag_err2s[:, 0]
            else:
                e2 = self._mag_err2s[:, 0]
            if 1 in self._fit_err_ratio_ind:
                e2 += (err_ratios**2.)*self._mag_err2s[:, 1]
            else:
                e2 += self._mag_err2s[:, 1]
            self._gsig = np.sqrt(np.clip(spl(self._redshifts), self._min_scatter, None)**2. + e2)

        if self._fit_scatter and self._trunc is not None:
            phi_bma = special.erf((self._trunc / self._gsig) / np.sqrt(2.))
        elif self._trunc is not None:
            phi_bma = self._phi_bma
        else:
            phi_bma = 1.0

        if self._has_dmags:
            model_color = gmean + gslope * self._dmags + self._lupcorrs
        else:
            model_color = gmean

        xi = (self._colors - model_color) / self._gsig

        phi = (1. / self._gsig) * (1. / np.sqrt(2. * np.pi)) * np.exp(-0.5 * xi**2.)

        with np.warnings.catch_warnings():
            np.warnings.simplefilter("ignore")
            # This might have a divide-by-zero, that's okay, we check later.
            gci = phi / phi_bma

        if self._has_probs:
            # Use probabilities and bkgs
            with np.warnings.catch_warnings():
                np.warnings.simplefilter("ignore")

                vals = np.log(self._probs * gci + (1.0 - self._probs) * self._bkgs)
        else:
            # No probabilities or bkgs
            with np.warnings.catch_warnings():
                np.warnings.simplefilter("ignore")

                vals = np.log(gci)

        bad, = np.where(~np.isfinite(vals))
        vals[bad] = -100.0

        if self._fit_scatter and self._use_scatter_prior:
            t = -(np.sum(vals) - np.sum(np.log(np.clip(pars[self._scatter_index: self._scatter_index + self._n_scatter_nodes], self._min_scatter, None))))
        else:
            t = -np.sum(vals)

        return t class RedSequenceOffDiagonalFitter(object):
    """
    Class to fit a spline to a pair of off-diagonal elements of the
    red-sequence covariance matrix.

    This can only be run after the diagonal elements of the red-sequence model
    have been constrained.
    """
    def __init__(self, nodes, redshifts, d1, d2, s1, s2, mag_errs, j, k, probs, bkgs, covmat_prior, min_eigenvalue=0.0):
        """
        Instantiate a RedSequenceOffDiagonalFitter

        Parameters
        ----------
        nodes: `np.array`
           Float array of covariance matrix redshift nodes
        redshifts: `np.array`
           Float array of input redshifts for fit
        d1: `np.array`
           Float array of color residual (color - model_color) for color 1 (j)
           Must have same number of elements as redshifts.
        d2: `np.array`
           Float array of color residual (color - model_color) for color 2 (k)
           Must have same number of elements as redshifts.
        s1: `np.array`
           Float array of intrinsic scatters for color 1 (j)
           Must have same number of elements as redshifts.
        s2: `np.array`
           Float array of intrinsic scatters for color 2 (k)
           Must have same number of elements as redshifts.
        mag_errs: `np.array`
           2d float array of magnitude errors [n_redshifts, nmag]
        j: `int`
           Index for color 1
        k: `int`
           Index for color 2
        probs: `np.array`
           Float array of membership probabilities
           Must have same number of elements as redshifts.
        bkgs: `np.array`
           Float array of background likelihood for color 1, color 2.
           Must have same number of elements as redshifts.
        covmat_prior: `float`
           Prior on covariance matrix off-diagonal elements.
        min_eigenvalue: `float`, optional
           Minimum eigenvalue of covariance matrix.  Default is 0.0.
        """
        self._nodes = np.atleast_1d(nodes).astype(np.float64)
        self._redshifts = np.atleast_1d(redshifts).astype(np.float64)
        self._d1 = np.atleast_1d(d1).astype(np.float64)
        self._d2 = np.atleast_1d(d2).astype(np.float64)
        self._s1 = np.atleast_1d(s1).astype(np.float64)
        self._s2 = np.atleast_1d(s2).astype(np.float64)
        self._probs = np.atleast_1d(probs).astype(np.float64)
        self._bkgs = np.atleast_1d(bkgs).astype(np.float64)

        if self._redshifts.size != self._d1.size:
            raise ValueError("Number of redshifts must be equal to d1")
        if self._redshifts.size != self._d2.size:
            raise ValueError("Number of redshifts must be equal to d2")
        if self._redshifts.size != self._s1.size:
            raise ValueError("Number of redshifts must be equal to s1")
        if self._redshifts.size != self._s2.size:
            raise ValueError("Number of redshifts must be equal to s2")
        if self._redshifts.size != self._probs.size:
            raise ValueError("Number of redshifts must be equal to probs")
        if self._redshifts.size != self._bkgs.size:
            raise ValueError("Number of redshifts must be equal to bkgs")

        self._j = j
        self._k = k

        if len(mag_errs.shape) != 2:
            raise ValueError("mag_errs must be 2d")
        if mag_errs.shape[0] != self._redshifts.size:
            raise ValueError("Number of redshifts must be number of mag_errs")

        self._covmat_prior = covmat_prior
        self._min_eigenvalue = min_eigenvalue

        self._c_int = np.zeros((2, 2, self._redshifts.size))
        self._c_int[0, 0, :] = self._s1**2.
        self._c_int[1, 1, :] = self._s2**2.

        self._c_noise = np.zeros_like(self._c_int)
        self._c_noise[0, 0, :] = mag_errs[:, j]**2. + mag_errs[:, j + 1]**2.
        self._c_noise[1, 1, :] = mag_errs[:, k]**2. + mag_errs[:, k + 1]**2.
        if k == (j + 1):
            self._c_noise[0, 1, :] = -mag_errs[:, k]**2.
            self._c_noise[1, 0, :] = self._c_noise[0, 1, :]

    def fit(self, p0, full_covmats=None):
        """
        Perform a spline fit to a pair of elements of the covariance matrix.

        Parameters
        ----------
        p0: `np.array`
           Initial fit parameters, with same number of elements as nodes
        full_covmats: `np.array`, optional
           Full covariance matrix node values [nmag, nmag, nnode]
           Default is None.  If not None, then minimum eigenvalue is checked during fit.

        Returns
        -------
        pars: `np.array`
           Correlation values (r) for each node.
        """

        self._full_covmats = full_covmats

        bounds = []
        for i in range(self._nodes.size):
            bounds.append((-0.9, 0.9))

        res = scipy.optimize.minimize(self,
                                      p0,
                                      method='L-BFGS-B',
                                      bounds=bounds,
                                      jac=False,
                                      options={'maxfun': 2000,
                                               'maxiter': 2000,
                                               'maxcor': 20,
                                               'eps': 1e-3,
                                               'gtol': 1e-8},
                                      callback=None)
        pars = res.x

        return pars

    def __call__(self, pars):
        """
        Compute the off-diagonal log-likelihood (negative for minimization)

        Parameters
        ----------
        pars: `np.array`
           Float array of the correlation (r) values

        Returns
        -------
        t: `float`
           Total negative log-likelihood
        """

        spl = CubicSpline(self._nodes, pars)
        r = np.clip(spl(self._redshifts), -0.9, 0.9)

        metrics = np.zeros((2, 2, self._redshifts.size))
        self._c_int[0, 1, :] = r * self._s1 * self._s2
        self._c_int[1, 0, :] = self._c_int[0, 1, :]

        if self._full_covmats is not None:
            self._full_covmats[self._j, self._k, :] = pars * np.sqrt(self._full_covmats[self._j, self._j, :]) * np.sqrt(self._full_covmats[self._k, self._k, :])
            self._full_covmats[self._k, self._j, :] = self._full_covmats[self._j, self._k, :]

        covmats = self._c_int + self._c_noise

        dets = covmats[0, 0, :] * covmats[1, 1, :] - covmats[0, 1, :] * covmats[1, 0, :]

        metrics[0, 0, :] = covmats[1, 1, :] / dets
        metrics[1, 1, :] = covmats[0, 0, :] / dets
        metrics[1, 0, :] = -covmats[0, 1, :] / dets
        metrics[0, 1, :] = -covmats[1, 0, :] / dets

        exponents = -0.5 * (metrics[0, 0, :] * self._d1 * self._d1 +
                            (metrics[0, 1, :] + metrics[1, 0, :]) * self._d1 * self._d2 +
                            metrics[1, 1, :] * self._d2 * self._d2)

        gci = (dets**(-0.5) / (2. * np.pi)) * np.exp(exponents)

        vals = np.log(self._probs * gci + (1. - self._probs ) * self._bkgs)

        bad, = np.where(~np.isfinite(vals))
        vals[bad] = -100

        t=-(np.sum(vals) - np.sum(0.5 * (pars / self._covmat_prior)**2.))

        if ~np.isfinite(t):
            t = 1e11
        else:
            wall = False

            # Check for negative eigenvalues
            if self._full_covmats is not None:
                for i in range(self._nodes.size):
                    a = self._full_covmats[:, :, i]
                    d = np.linalg.eigvalsh(a)
                    if (np.min(d) < self._min_eigenvalue):
                        wall = True
            if wall:
                t += 100000

        return t class CorrectionFitter(object):
    """
    Class for fit a spline to the zred corrections as a function of redshift.

    This class computes the zred corrections, including slope (if desired) and
    error adjustment factor.
    """
    def __init__(self, mean_nodes, redshifts, dzs, dz_errs,
                 slope_nodes=None, r_nodes=None, bkg_nodes=None,
                 probs=None, dmags=None, ws=None):
        """
        Instantiate a CorrectionFitter

        Parameters
        ----------
        mean_nodes: `np.array`
           Float array for mean correction redshift nodes
        redshifts: `np.array`
           Float array of input redshifts for fit
        dzs: `np.array`
           Float array of input delta-z (zred - ztrue) for fit
        dz_errs: `np.array`
           Float array of error on delta-z
        slope_nodes: `np.array`, optional
           Float array for slope redshift nodes.
           Default is None (use mean_nodes for slope).
        r_nodes: `np.array`, optional
           Float array for error factor redshift nodes.
           Default is None (use slope_nodes for correction factor).
        bkg_nodes: `np.array`, optional
           Float array for background/outlier redshift nodes.
           Default is None (use slope_nodes for bkg.)
        probs: `np.array`, optional
           Float array for membership probabilities.
           Default is None (no probabilities used).
        dmags: `np.array`, optional
           Float array of delta-mag (refmag - pivot) (location on x axis for slope).
           Must have same length as redshifts if used.
           Default is None (don't fit slope).
        ws: `np.array`, optional
           Float array of likelihood weighting factors.
           Default is None (no weighting factors).
        """
        self._mean_nodes = np.atleast_1d(mean_nodes).astype(np.float64)
        # Note that the slope_nodes are the default for the r, bkg as well
        if slope_nodes is None:
            self._slope_nodes = self._mean_nodes
        else:
            self._slope_nodes = np.atleast_1d(slope_nodes).astype(np.float64)
        if r_nodes is None:
            self._r_nodes = self._slope_nodes
        else:
            self._r_nodes = np.atleast_1d(r_nodes).astype(np.float64)
        if bkg_nodes is None:
            self._bkg_nodes = self._slope_nodes
        else:
            self._bkg_nodes = np.atleast_1d(bkg_nodes).astype(np.float64)

        self._n_mean_nodes = self._mean_nodes.size
        self._n_slope_nodes = self._slope_nodes.size
        self._n_r_nodes = self._r_nodes.size
        self._n_bkg_nodes = self._bkg_nodes.size

        self._redshifts = np.atleast_1d(redshifts).astype(np.float64)
        self._dzs = np.atleast_1d(dzs).astype(np.float64)
        self._dz_errs = np.atleast_1d(dz_errs).astype(np.float64)

        if self._redshifts.size != self._dzs.size:
            raise ValueError("Number of redshifts must be equal to dzs")
        if self._redshifts.size != self._dz_errs.size:
            raise ValueError("Number of redshifts must be equal to dz_errs")

        if probs is not None:
            self._probs = np.atleast_1d(probs).astype(np.float64)
            if self._redshifts.size != self._probs.size:
                raise ValueError("Number of redshifts must be equal to probs")
        else:
            self._probs = np.ones_like(self._redshifts)

        if dmags is not None:
            self._dmags = np.atleast_1d(dmags).astype(np.float64)
            if self._redshifts.size != self._dmags.size:
                raise ValueError("Number of redshifts must be equal to dmags")
        else:
            self._dmags = np.zeros_like(self._redshifts)

        if ws is not None:
            self._ws = np.atleast_1d(ws).astype(np.float64)
            if self._redshifts.size != self._ws.size:
                raise ValueError("Number of redshifts must be equal to ws")
        else:
            self._ws = np.ones_like(self._redshifts)

    def fit(self, p0_mean, p0_slope, p0_r, p0_bkg, fit_mean=False, fit_slope=False, fit_r=False, fit_bkg=False):
        """
        Fit the zred correction factor.

        Parameters
        ----------
        p0_mean: `np.array`
           Float array of initial fit parameters for mean offset relation
           Must have same number of elements as mean_nodes
        p0_slope: `np.array`
           Float array of initial fit parameters for slope relation
           Must have same number of elements as slope_nodes
        p0_r: `np.array`
           Float array of initial fit parameters for error scaling relation
           Must have same number of elements as r_nodes
        p0_bkg: `np.array`
           Float array of initial fit parameters for bkg outlier relation
           Must have same number of elements as bkg_nodes
        fit_mean: `bool`, optional
           Fit the mean relation? (else fix to p0_mean).  Default is False.
        fit_slope: `bool`, optional
           Fit the slope relation? (else fix to p0_slope).  Default is False.
        fit_r: `bool`, optional
           Fit the r relation? (else fix to p0_r).  Default is False.
        fit_bkg: `bool`, optional
           Fit the bkg relation? (else fix to p0_bkg).  Default is False.

        Returns
        -------
        pars_mean: `np.array`
           Mean parameters (if fit_mean is True)
        pars_slope: `np.array`
           Slope parameters (if fit_slope is True)
        pars_r: `np.array`
           Error scaling parameters (if fit_r is True)
        pars_bkg: `np.array`
           Background parameters (if fit_bkg is True)
        """
        self._fit_mean = fit_mean
        self._fit_slope = fit_slope
        self._fit_r = fit_r
        self._fit_bkg = fit_bkg

        ctr = 0
        p0 = np.array([])
        bounds = []
        if self._fit_mean:
            self._mean_index = 0
            ctr += self._n_mean_nodes
            p0 = np.append(p0, p0_mean)
            for i in range(self._n_mean_nodes):
                bounds.append([-np.inf, np.inf])
        if self._fit_slope:
            self._slope_index = ctr
            ctr += self._n_slope_nodes
            p0 = np.append(p0, p0_slope)
            for i in range(self._n_slope_nodes):
                bounds.append([-np.inf, np.inf])
        if self._fit_r:
            self._r_index = ctr
            ctr += self._n_r_nodes
            p0 = np.append(p0, p0_r)
            for i in range(self._n_r_nodes):
                bounds.append([0.1, 2.0])
        if self._fit_bkg:
            self._bkg_index = ctr
            ctr += self._n_bkg_nodes
            p0 = np.append(p0, p0_bkg)
            for i in range(self._n_bkg_nodes):
                bounds.append([0.0, np.inf])

        if ctr == 0:
            raise ValueError("Must select at least one of fit_mean, fit_slope")

        # Precompute
        if not self._fit_mean:
            spl = CubicSpline(self._mean_nodes, p0_mean)
            self._gmean = spl(self._redshifts)
        if not self._fit_slope:
            spl = CubicSpline(self._slope_nodes, p0_slope)
            self._gslope = spl(self._redshifts)
        if not self._fit_r:
            spl = CubicSpline(self._r_nodes, p0_r)
            self._gr = spl(self._redshifts)
        if not self._fit_bkg:
            spl = CubicSpline(self._bkg_nodes, p0_bkg)
            self._gbkg = np.clip(spl(self._redshifts), 1e-10, None)
            self._gci1 = (1. / np.sqrt(2. * np.pi * self._gbkg)) * np.exp(-self._dzs**2. / (2. * self._gbkg))

        res = scipy.optimize.minimize(self,
                                      p0,
                                      method='L-BFGS-B',
                                      bounds=bounds,
                                      jac=False,
                                      options={'maxfun': 5000,
                                               'maxiter': 5000,
                                               'maxcor': 20,
                                               'eps': 1e-5,
                                               'gtol': 1e-10},
                                      callback=None)
        pars = res.x

        retval = []
        if self._fit_mean:
            retval.append(pars[self._mean_index: self._mean_index + self._n_mean_nodes])
        if self._fit_slope:
            retval.append(pars[self._slope_index: self._slope_index + self._n_slope_nodes])
        if self._fit_r:
            retval.append(pars[self._r_index: self._r_index + self._n_r_nodes])
        if self._fit_bkg:
            retval.append(pars[self._bkg_index: self._bkg_index + self._n_bkg_nodes])

        return retval

    def __call__(self, pars):
        """
        Compute the correction log-likelihood (negative for minimization)

        Parameters
        ----------
        pars: `np.array`
           Float array of the consolidate parameters

        Returns
        -------
        t: `float`
           Total negative log-likelihood
        """

        if self._fit_mean:
            spl = CubicSpline(self._mean_nodes, pars[self._mean_index: self._mean_index + self._n_mean_nodes])
            gmean = spl(self._redshifts)
        else:
            gmean = self._gmean

        if self._fit_slope:
            spl = CubicSpline(self._slope_nodes, pars[self._slope_index: self._slope_index + self._n_slope_nodes])
            gslope = spl(self._redshifts)
        else:
            gslope = self._gslope

        if self._fit_r:
            spl = CubicSpline(self._r_nodes, pars[self._r_index: self._r_index + self._n_r_nodes])
            gr = spl(self._redshifts)
        else:
            gr = self._gr

        if self._fit_bkg:
            #if pars[self._bkg_index: self._bkg_index + self._n_bkg_nodes].min() < 0.0:
            #    return 1e11
            spl = CubicSpline(self._bkg_nodes, pars[self._bkg_index: self._bkg_index + self._n_bkg_nodes])
            gbkg = np.clip(spl(self._redshifts), 1e-10, None)
            gci1 = (1. / np.sqrt(2. * np.pi * gbkg)) * np.exp(-self._dzs**2. / (2. * gbkg))
        else:
            gbkg = self._gbkg
            gci1 = self._gci1

        var0 = (gr * self._dz_errs)**2.
        gci0 = (1. / np.sqrt(2. * np.pi * var0)) * np.exp(-(self._dzs - (gmean + gslope * self._dmags))**2. / (2. * var0))

        vals = self._ws * (self._probs * gci0 + (1. - self._probs) * gci1)

        bad, = np.where((~np.isfinite(vals)) | (vals <= 0.0))
        vals[bad] = 4e-44

        vals = np.log(vals)

        t = -np.sum(vals)

        return t class ErrorBinFitter(object):
    """Class for fitting error ratios in bins, using median statistics.

    Parameters
    ----------
    delta_col : `np.ndarray`
        Array of delta colors
    delta_mag : `np.ndarray`
        Array of delta mag (mag - pivot)
    err_0 : `np.ndarray`
        Raw or scaled error for first magnitude in color
    err_1 : np.ndarray`
        Raw or scaled error for second magnitude in color
    sigint2 : `np.ndarray`
        Intrinsic scatter squared.
    binsize : `float`
        Bin size.
    ntrial : `int`
        Number of trials for bootstrap errors.
    """
    def __init__(self, delta_col, delta_mag, err_0, err_1, sigint2, binsize=0.5, ntrial=100):
        self.delta_col = delta_col
        self.delta_mag = delta_mag
        self.err_0 = err_0
        self.err_1 = err_1
        self.sigint2 = sigint2

        h_full = esutil.stat.histogram(delta_mag, binsize=binsize, more=True)

        self.nbin = h_full['nbin']
        self.binmag = h_full['center']
        self.rev = h_full['rev']

        delta_err = np.sqrt(self.sigint2 +
                            self.err_0**2. +
                            self.err_1**2.)
        pulls = self.delta_col/delta_err

        self.mad_err = np.zeros(self.nbin)
        for i in range(self.nbin):
            bin_mads = np.zeros(ntrial)
            i1a = self.rev[self.rev[i]: self.rev[i + 1]]
            for j in range(ntrial):
                r = np.random.choice(i1a, size=i1a.size, replace=True)

                med = np.median(pulls[r])
                bin_mads[j] = 1.4826*np.median(np.abs(pulls[r] - med))

            mad = np.median(bin_mads)
            self.mad_err[i] = 1.4826*np.median(np.abs(bin_mads - mad))

        # We require a positive error estimate (if too few in a bin)
        self.use_bins = np.where(self.mad_err > 0.0)

    def fit(self, p0, scale_indices=[0]):
        """Perform a fit to the error scaling as a function of magnitude.

        Parameters
        ----------
        p0 : array-like
            Initial fit parameters, intercept and slope.
        scale_indices : `list`
            List of mag indices to scale.  Can be [0], [1], or [0, 1]

        Returns
        -------
        pars : `np.ndarray`
            Fit parameters, intercept and slope.
        """
        self._scale_indices = scale_indices

        bounds = [[0.5, 10.0],
                  [-5.0, 5.0]]

        res = scipy.optimize.minimize(self,
                                      p0,
                                      method='L-BFGS-B',
                                      bounds=bounds,
                                      jac=False,
                                      options={'maxfun': 2000,
                                               'maxiter': 2000,
                                               'maxcor': 20,
                                               'eps': 1e-5,
                                               'gtol': 1e-8},
                                      callback=None)
        return res.x

    def __call__(self, pars):
        """Compute the chi2 for the pars.

        Parameters
        ----------
        pars : `list`
            Parameters, intercept and slope.

        Returns
        -------
        chi2 : `float`
        """
        mad = np.zeros(self.nbin)

        err_ratio = pars[0] + pars[1]*self.delta_mag

        if 0 in self._scale_indices:
            scaled_err_0 = err_ratio*self.err_0
        else:
            scaled_err_0 = self.err_0
        if 1 in self._scale_indices:
            scaled_err_1 = err_ratio*self.err_1
        else:
            scaled_err_1 = self.err_1

        delta_err = np.sqrt(self.sigint2 +
                            scaled_err_0**2. +
                            scaled_err_1**2.)
        pulls = self.delta_col/delta_err

        for i in range(self.nbin):
            i1a = self.rev[self.rev[i]: self.rev[i + 1]]
            med = np.median(pulls[i1a])
            mad[i] = 1.4826*np.median(np.abs(pulls[i1a] - med))

        chi2 = np.sum((mad[self.use_bins] - 1.0)**2./self.mad_err[self.use_bins]**2., dtype=np.float64)
        return chi2 class RedSequenceColorPar(object):
    """
    Class which describes the color-based red-sequence model.

    This is the fundamental basis of the redmapper red sequence model.
    """

    def __init__(self, filename, zbinsize=None, minsig=0.01, fine=False, zrange=None, config=None, limmag=None):
        """
        Instantiate a RedSequenceColorPar object.

        The primary way to create a RedSequenceColorPar object is by specifying
        a filename to read the red-sequence parameters from.  The default reads
        in with coarse interpolation, which is faster, uses less memory, and is
        used when computing zreds.  The 'fine=True' option reads fine-scale
        binning which is slower, uses more memory, and is used when computing
        richness and cluster z_lambda.

        If a red-sequence file is not available, a placeholder
        RedSequenceColorPar can be created by setting 'filename=None' and
        specifying a Configuration.  Some features of the RedSequenceColorPar
        (including mstar calculations) are available in this mode.

        Parameters
        ----------
        filename: `str`
           Filename of the fits file to load parameters from
        zbinsize: `float`, optional
           Redshift binning to interpolate model.
           Default is None, which uses ZBINCOAR (coarse binning)
           from the fits header, or ZBINFINE (fine binning) if fine=True
        minsig: `float`, optional
           Minimum intrinsic scatter.  Default is 0.01 mag.
        fine: `bool`, optional
           Use fine binning for interpolation.  Default is False.
        zrange: `np.array`, optional
           Redshift range to do interpolation.
           Default is None, which means use the config (if specified)
           or ZRANGE0, ZRANGE1 from the filename header
        config: `redmapper.Configuration`, optional
           Only required if filename=None, and then used for placeholder
           RedSequenceColorPar
        limmag: `float`, optional
           Maximum magnitude to do red-sequence interpolation.
        """

        if filename is None:
            if config is None:
                raise ValueError("Must have either filename or config")
            if limmag is None:
                limmag = config.limmag_catalog
            else:
                limmag = limmag
            if zrange is None:
                zrange = config.zrange
            alpha = config.calib_lumfunc_alpha
            mstar_survey = config.mstar_survey
            mstar_band = config.mstar_band
            ncol = config.nmag - 1
            templatefile = config.calib_redgal_template
            bands = config.bands

            # A compromize zbinsize
            zbinsize = 0.001

            has_file = False
        else:
            pars,hdr=fitsio.read(filename, ext=1, header=True, upper=True)
            try:
                if limmag is None:
                    limmag = hdr['LIMMAG']
                else:
                    limmag = limmag
                if (zrange is None):
                    zrange = np.array([hdr['ZRANGE0'],hdr['ZRANGE1']])
                alpha = hdr['ALPHA']
                mstar_survey = hdr['MSTARSUR']
                mstar_band = hdr['MSTARBAN']
                ncol = hdr['NCOL']
                if 'TEMPLATE' in hdr:
                    templatefile = hdr['TEMPLATE']
                    if ',' in hdr['BANDS']:
                        bands = hdr['BANDS'].rstrip().split(',')
                    else:
                        bands = list(hdr['BANDS'].rstrip())
                else:
                    templatefile = None
            except:
                raise ValueError("Missing field from parameter header.")
            has_file = True

        if len(zrange) != 2:
            raise ValueError("zrange must have 2 elements");

        if zbinsize is None:
            try:
                if fine:
                    zbinsize=hdr['ZBINFINE']
                else:
                    zbinsize=hdr['ZBINCOAR']
            except:
                raise ValueError("Missing field from parameter header.")

        try:
            lowzmode=hdr['LOWZMODE']
        except:
            lowzmode = 0

        nmag = ncol + 1
        self.nmag = nmag

        if has_file:
            bvalues=np.zeros(nmag)
            try:
                for i in range(nmag):
                    bvalues[i] = hdr['BVALUE%1d' % (i+1)]
            except:
                bvalues[:] = 0.0

            do_lupcorr = False
            if bvalues.min() > 0.0:
                # Only do luptitude corrections if we have non-zero b values.
                do_lupcorr = True

            try:
                ref_ind = hdr['REF_IND']
            except:
                try:
                    ref_ind = hdr['I_IND']
                except:
                    raise ValueError("Need REF_IND or I_IND")

        if templatefile is not None:
            rg = RedGalInitialColors(templatefile)

        nz = np.round((zrange[1]-zrange[0])/zbinsize).astype('i4') #nr of bins
        self.z = zbinsize*np.arange(nz) + zrange[0]                #z bins

        # append a high-end overflow bin
        # for computation this will be the same as the top bin, and at the end
        # we set it to some absurdly large number.
        self.z = np.append(self.z,self.z[self.z.size-1])
        nz=nz+1

        self.zbinsize=zbinsize
        self.zbinscale=int(1./zbinsize)

        ms=MStar(mstar_survey,mstar_band)

        refmagbinsize=0.01
        if (lowzmode):
            refmagrange=np.array([10.0,limmag],dtype='f4')
            # FIXME
            lumrefmagrange=np.array([10.0,ms(zrange[1])-2.5*np.log10(0.1)])
        else:
            refmagrange=np.array([12.0,limmag],dtype='f4')
            lumrefmagrange=np.array([12.0,ms(zrange[1])-2.5*np.log10(0.1)])
        self.refmagbins = np.arange(refmagrange[0], refmagrange[1], refmagbinsize, dtype='f8')
        self.lumrefmagbins=np.arange(lumrefmagrange[0],lumrefmagrange[1],refmagbinsize,dtype='f8')

        # and for fast look-ups...
        self.refmagbins = np.append(self.refmagbins,self.refmagbins[self.refmagbins.size-1])
        self.lumrefmagbins = np.append(self.lumrefmagbins,self.lumrefmagbins[self.lumrefmagbins.size-1])

        self.refmagbinsize = refmagbinsize
        self.refmagbinscale = int(1./refmagbinsize)
        self.refmaginteger = (self.refmagbins*self.refmagbinscale).astype(np.int64)
        self.lumrefmaginteger = (self.lumrefmagbins*self.refmagbinscale).astype(np.int64)

        if has_file:
            # is this an old or new structure?
            if 'PIVOTMAG_Z' in pars.dtype.names:
                refmag_name = 'REFMAG'
                pivotmag_name = 'PIVOTMAG'
            else :
                refmag_name = 'IMAG'
                pivotmag_name = 'REFMAG'

            # mark the extrapolated values
            self.extrapolated = np.zeros(nz,dtype=bool)
            loz,=np.where(self.z < np.min(pars[pivotmag_name+'_Z']))
            hiz,=np.where(self.z > np.max(pars[pivotmag_name+'_Z']))
            if (loz.size > 0) : self.extrapolated[loz] = True
            if (hiz.size > 0) : self.extrapolated[hiz] = True

            # set the pivotmag
            self.pivotmag = np.zeros(self.z.size, dtype=np.float64)
            spl=CubicSpline(pars[0][pivotmag_name+'_Z'],pars[0][pivotmag_name])
            self.pivotmag[:] = spl(self.z)

            # and the max/min refmag
            spl=CubicSpline(pars[0][pivotmag_name+'_Z'], pars[0]['MAX'+refmag_name])
            self.maxrefmag = spl(self.z)
            spl=CubicSpline(pars[0][pivotmag_name+'_Z'], pars[0]['MIN'+refmag_name])
            self.minrefmag = spl(self.z)

            # c/slope
            self.c = np.zeros((nz,ncol),dtype=np.float64)
            self.slope = np.zeros((nz,ncol),dtype=np.float64)
            for j in range(ncol):
                jstring='%02d' % (j)
                spl=CubicSpline(pars[0]['Z'+jstring],pars[0]['C'+jstring])
                self.c[:,j] = spl(self.z)
                spl=CubicSpline(pars[0]['ZS'+jstring],pars[0]['SLOPE'+jstring])
                self.slope[:,j] = spl(self.z)

            # sigma/covmat
            self.sigma = np.zeros((ncol,ncol,nz),dtype=np.float64)
            self.covmat = np.zeros((ncol,ncol,nz),dtype=np.float64)

            self.mag_err_ratio_intercept = np.ones(nmag, dtype=np.float64)
            self.mag_err_ratio_slope = np.zeros(nmag, dtype=np.float64)
            self.mag_err_ratio_pivot = 20.0  # arbitrary, with zero slope.
            if 'MAG_ERR_RATIO_INT' in pars.dtype.names:
                self.mag_err_ratio_intercept[:] = pars[0]['MAG_ERR_RATIO_INT'][:]
                self.mag_err_ratio_slope[:] = pars[0]['MAG_ERR_RATIO_SLOPE'][:]
                self.mag_err_ratio_pivot = pars[0]['MAG_ERR_RATIO_PIVOT']

            # diagonals
            for j in range(ncol):
                spl=CubicSpline(pars[0]['COVMAT_Z'],pars[0]['SIGMA'][j,j,:])
                self.sigma[j,j,:] = np.clip(spl(self.z), minsig, None)

                self.covmat[j,j,:] = self.sigma[j,j,:]*self.sigma[j,j,:]

            # off-diagonals
            for j in range(ncol):
                for k in range(j+1,ncol):
                    spl=CubicSpline(pars[0]['COVMAT_Z'],pars[0]['SIGMA'][j,k,:])
                    self.sigma[j,k,:] = spl(self.z)

                    too_high,=np.where(self.sigma[j,k,:] > 0.99)
                    if (too_high.size > 0):
                        self.sigma[j,k,too_high] = 0.99
                    too_low,=np.where(self.sigma[j,k,:] < -0.99)
                    if (too_low.size > 0):
                        self.sigma[j,k,too_low] = -0.99

                    self.sigma[k,j,:] = self.sigma[j,k,:]

                    self.covmat[j,k,:] = self.sigma[k,j,:] * self.sigma[j,j,:] * self.sigma[k,k,:]
                    self.covmat[k,j,:] = self.covmat[j,k,:]

            # volume factor
            spl=CubicSpline(pars[0]['VOLUME_FACTOR_Z'],pars[0]['VOLUME_FACTOR'])
            self.volume_factor = spl(self.z)

            # corrections
            spl=CubicSpline(pars[0]['CORR_Z'],pars[0]['CORR'])
            self.corr = spl(self.z)
            spl=CubicSpline(pars[0]['CORR_SLOPE_Z'],pars[0]['CORR_SLOPE'])
            self.corr_slope = spl(self.z)

            spl=CubicSpline(pars[0]['CORR_Z'],pars[0]['CORR2'])
            self.corr2 = spl(self.z)
            spl=CubicSpline(pars[0]['CORR_SLOPE_Z'],pars[0]['CORR2_SLOPE'])
            self.corr2_slope = spl(self.z)

            if 'CORR_R' in pars.dtype.names:
                # protect against stupidity
                if (pars[0]['CORR_R'][0] <= 0.0) :
                    self.corr_r = np.ones(nz)
                else:
                    spl=CubicSpline(pars[0]['CORR_SLOPE_Z'],pars[0]['CORR_R'])
                    self.corr_r = spl(self.z)

                test,=np.where(self.corr_r < 0.5)
                if (test.size > 0) : self.corr_r[test] = 0.5

                if (pars[0]['CORR2_R'][0] <= 0.0):
                    self.corr2_r = np.ones(nz)
                else:
                    spl=CubicSpline(pars[0]['CORR_SLOPE_Z'],pars[0]['CORR2_R'])
                    self.corr2_r = spl(self.z)

                test,=np.where(self.corr2_r < 0.5)
                if (test.size > 0) : self.corr2_r[test] = 0.5

            else :
                self.corr_r = np.ones(nz)
                self.corr2_r = np.ones(nz)

            # Deal with extrapolation in the case where we have a template
            if templatefile is not None:
                not_extrap, = np.where(~self.extrapolated)

                # Pin the pivotmag.

                self.pivotmag[loz] = self.pivotmag[not_extrap[0]]
                self.pivotmag[hiz] = self.pivotmag[not_extrap[-1]]

                # Leave the max/min refmag extrapolated

                # Pin the slope.
                for j in range(ncol):
                    self.slope[loz, j] = self.slope[not_extrap[0], j]
                    self.slope[hiz, j] = self.slope[not_extrap[-1], j]

                # Leave the covariance matrix extrapolated?

                # Replace the intercept with the template values (offset)
                for j in range(ncol):
                    try:
                        delta = (rg(bands[j], bands[j + 1], self.z[not_extrap[0]]) -
                                 self.c[not_extrap[0], j])
                        self.c[loz, j] = rg(bands[j], bands[j + 1], self.z[loz]) - delta

                        delta = (rg(bands[j], bands[j + 1], self.z[not_extrap[-1]]) -
                                 self.c[not_extrap[-1], j])
                        self.c[hiz, j] = rg(bands[j], bands[j + 1], self.z[hiz]) - delta
                    except ValueError:
                        # This is fine, we won't use the template extrapolation
                        pass

        # mstar
        # create LUT
        self._mstar = ms(self.z)

        # luminosity function integrations
        self.lumnorm = np.zeros((self.lumrefmagbins.size,nz))
        self.alpha = alpha
        for i in range(nz):
            f = schechter_pdf(self.lumrefmagbins, alpha=self.alpha, mstar=self._mstar[i])
            self.lumnorm[:, i] = refmagbinsize*np.cumsum(f, dtype=np.float64)

        if has_file:
            # lupcorr (annoying!)
            self.lupcorr = np.zeros((self.refmagbins.size,nz,ncol),dtype='f8')
            if (do_lupcorr):
                bnmgy = bvalues*1e9

                for i in range(nz):
                    mags = np.zeros((self.refmagbins.size,nmag))
                    lups = np.zeros((self.refmagbins.size,nmag))

                    mags[:,ref_ind] = self.refmagbins

                    # go redward
                    for j in range(ref_ind+1,nmag):
                        mags[:,j] = mags[:,j-1] - (self.c[i,j-1]+self.slope[i,j-1]*(mags[:,ref_ind]-self.pivotmag[i]))
                    # blueward
                    for j in range(ref_ind-1,-1,-1):
                        mags[:,j] = mags[:,j+1] + (self.c[i,j]+self.slope[i,j]*(mags[:,ref_ind]-self.pivotmag[i]))

                    # and the luptitude conversion
                    for j in range(nmag):
                        flux = 10.**((mags[:,j]-22.5)/(-2.5))
                        lups[:,j] = 2.5*np.log10(1.0/bvalues[j]) - np.arcsinh(0.5*flux/bnmgy[j])/(0.4*np.log(10.0))

                    magcol = mags[:,0:ncol] - mags[:,1:ncol+1]
                    lupcol = lups[:,0:ncol] - lups[:,1:ncol+1]

                    self.lupcorr[:,i,:] = lupcol - magcol

        # set top overflow bins to very large number
        self.z[self.z.size-1] = 1000.0
        self.zinteger=np.round(self.z*self.zbinscale).astype(np.int64)
        self.refmagbins[self.refmagbins.size-1] = 1000.0
        self.refmaginteger = np.round(self.refmagbins*self.refmagbinscale).astype(np.int64)
        self.lumrefmagbins[self.lumrefmagbins.size-1] = 1000.0
        self.lumrefmaginteger = np.round(self.lumrefmagbins*self.refmagbinscale).astype(np.int64)
        self.ncol = ncol
        self.alpha = alpha
        self.mstar_survey = mstar_survey
        self.mstar_band = mstar_band
        self.limmag = limmag

    def mstar(self,z):
        """
        Look up mstar at a set of redshifts

        Parameters
        ----------
        z: `np.array`
           Float array of redshifts

        Returns
        -------
        mstar: `np.array`
           Float array of mstar values.
        """
        # lookup and return mstar.
        zind = self.zindex(z)
        return self._mstar[zind]

    def zindex(self,z):
        """
        Look up the redshift index for the RedSequenceColorPar interpolated
        arrays.

        Parameters
        ----------
        z: `np.array`
           Float array of redshifts

        Returns
        -------
        zindex: `np.array`
           Integer array of redshift indices
        """
        # return the z index/indices with rounding.

        zind = np.searchsorted(self.zinteger,np.round(np.atleast_1d(z)*self.zbinscale).astype(np.int64))
        if (zind.size == 1):
            return np.ndarray.item(zind)
        else:
            return zind

    def refmagindex(self,refmag):
        """
        Look up the reference magnitude index for the RedSequenceColorPar
        interpolated arrays.

        Parameters
        ----------
        refmag: `np.array`
           Float array of refmag values

        Returns
        -------
        indices: `np.array`
           Integer array of refmag indices
        """
        # return the refmag index/indices with rounding

        refmagind = np.searchsorted(self.refmaginteger,np.round(np.atleast_1d(refmag)*self.refmagbinscale).astype(np.int64))
        if (refmagind.size == 1):
            return np.ndarray.item(refmagind)
        else:
            return refmagind

    def lumrefmagindex(self,lumrefmag):
        """
        Look up the luminosity table reference magnitude index.

        Parameters
        ----------
        lumrefmag: `np.array`
           Float array of refmags for luminosity table.

        Returns
        -------
        indices: `np.array`
           Integer array of lumrefmag indices
        """
        lumrefmagind = np.searchsorted(self.lumrefmaginteger,np.round(np.atleast_1d(lumrefmag)*self.refmagbinscale).astype(np.int64))
        if (lumrefmagind.size == 1):
            return np.ndarray.item(lumrefmagind)
        else:
            return lumrefmagind

    def calculate_chisq_redshifts(self, galaxy, zs, calc_lkhd=False, z_is_index=False):
        """
        Compute chisq for a single galaxy at a series of redshifts.  Optimized for the
        redshift case.

        Parameters
        ----------
        galaxy: `redmapper.Galaxy`
           The galaxy to compute chisq values
        zs: `np.array`
           Float array of redshifts or integer array of bins
        z_is_index: `bool`, optional
           The zs are indices and not redshifts.  Default is False.
        calc_lkhd: `bool`, optional
           Calculate likelihood rather than chisq.  Default is False.

        Returns
        -------
        chisqs: `np.array`
           Float array of chisq values.
        """

        calc_chisq = not calc_lkhd

        if z_is_index:
            zinds = zs
        else:
            zinds = self.zindex(zs)

        magind = self.refmagindex(galaxy.refmag)
        galcolor = galaxy.galcol

        mag_err = galaxy.mag_err.copy()
        dmag = galaxy.refmag - self.mag_err_ratio_pivot
        mag_err *= (self.mag_err_ratio_intercept[:] + self.mag_err_ratio_slope[:]*dmag)

        return compute_chisq(self.covmat[:,:,zinds], self.c[zinds,:],
                             self.slope[zinds,:], self.pivotmag[zinds],
                             np.array(galaxy.refmag), mag_err,
                             galcolor,
                             refmagerr=np.array(galaxy.refmag_err),
                             lupcorr=self.lupcorr[magind,zinds,:],
                             calc_chisq=calc_chisq, calc_lkhd=calc_lkhd)


    def calculate_chisq(self, galaxies, z, calc_lkhd=False, z_is_index=False):
        """
        Compute chisq for (a) a set of galaxies at redshift z or (b) a galaxy
        at an array of redshifts or (c) many galaxies at many redshifts.

        Parameters
        ----------
        galaxies: `redmapper.GalaxyCatalog`
           Catalog of galaxies to compute chisq values.
        z: `np.array`
           Float array of redshifts or integer array of redshift indices
        z_is_index: `bool`, optional
           The zs are indices and not redshifts.  Default is False.
        calc_lkhd: `bool`, optional
           Calculate likelihood rather than chisq.  Default is False.

        Returns
        -------
        chisqs: `np.array`
           Float array of chisq values.
        """

        if calc_lkhd:
            calc_chisq = False
        else:
            calc_chisq = True

        if z_is_index:
            zind = z
        else:
            zind = self.zindex(z)
        magind = self.refmagindex(galaxies.refmag)
        galcolor = galaxies.galcol

        # Need to check if this is a single galaxy AND redshift
        if np.atleast_1d(zind).size == 1 and len(galaxies) == 1:
            zind = np.atleast_1d(zind)
            magind = np.atleast_1d(magind)

        mag_err = galaxies.mag_err.copy()
        dmags = galaxies.refmag - self.mag_err_ratio_pivot
        for j in range(self.nmag):
            mag_err[:, j] *= self.mag_err_ratio_intercept[j] + self.mag_err_ratio_slope[j]*dmags

        return compute_chisq(self.covmat[:,:,zind], self.c[zind,:],
                             self.slope[zind,:], self.pivotmag[zind],
                             galaxies.refmag, mag_err,
                             galcolor,
                             refmagerr=galaxies.refmag_err,
                             lupcorr=self.lupcorr[magind,zind,:],
                             calc_chisq=calc_chisq, calc_lkhd=calc_lkhd)


    def plot_redsequence_diag(self, fig, ind, bands):
        """
        Plot the diagonal elements of the red-sequence model

        Parameters
        ----------
        fig: `matplotlib.Figure`
           Figure to add subplots to plot red-sequence model
        ind: `int`
           Color index to plot
        bands: `list`
           List of string names of bands for labeling
        """

        not_extrap, = np.where(~self.extrapolated)

        ax = fig.add_subplot(221)
        ax.plot(self.z[: -1], self.c[: -1, ind], 'r--')
        ax.plot(self.z[not_extrap], self.c[not_extrap, ind], 'r-')
        ax.set_xlabel('Redshift')
        ax.set_ylabel('<%s - %s>' % (bands[ind], bands[ind + 1]))

        ax = fig.add_subplot(222)
        ax.plot(self.z[: -1], self.slope[: -1, ind], 'r--')
        ax.plot(self.z[not_extrap], self.slope[not_extrap, ind], 'r-')
        ax.set_xlabel('Redshift')
        ax.set_ylabel('(%s - %s) slope' % (bands[ind], bands[ind + 1]))

        ax = fig.add_subplot(223)
        ax.plot(self.z[: -1], self.sigma[ind, ind, : -1], 'r--')
        ax.plot(self.z[not_extrap], self.sigma[ind, ind, not_extrap], 'r-')
        ax.set_ylim(bottom=0.0)
        ax.set_xlabel('Redshift')
        ax.set_ylabel('(%s - %s) sigma' % (bands[ind], bands[ind + 1]))

        fig.tight_layout()

    def plot_redsequence_offdiags(self, fig, bands):
        """
        Plot the off-diagonal elements of the red-sequence model

        Parameters
        ----------
        fig: `matplotlib.Figure`
           Figure to add subplots to plot red-sequence model
        bands: `list`
           List of string names of bands for labeling
        """

        noff = (self.ncol * self.ncol - self.ncol) // 2

        nrow = (noff + 1) // 2

        not_extrap, = np.where(~self.extrapolated)

        ctr = 1
        for j in range(self.ncol):
            for k in range(j + 1, self.ncol):
                ax = fig.add_subplot(nrow, 2, ctr)
                ax.plot(self.z[: -1], self.sigma[j, k, : -1], 'r--')
                ax.plot(self.z[not_extrap], self.sigma[j, k, not_extrap], 'r-')
                ax.set_ylim(-1.0, 1.0)
                ax.set_xlabel('Redshift')
                ax.set_ylabel('Corr %s-%s / %s-%s' % (bands[j], bands[j + 1], bands[k], bands[k + 1]))
                ctr += 1

        fig.tight_layout()

    def __repr__(self):
        return "Representation here." class ColorBackground(object):
    """
    Color-only galaxy background class.

    This class describes the binned, interpolateable background term b(x), where x describes the color(s) and reference magnitude of the galaxy.

    This contains all individual colors and color/color cross terms in a
    matrix.
    """

    def __init__(self, filename, usehdrarea=False):
        """
        Instantiate a color background

        Parameters
        ----------
        filename: `string`
           Color background filename
        usehdrarea: `bool`, optional
           Use area from the header instead of as configured.  Default is False.
        """

        refmagbinsize = 0.01
        colbinsize = 0.01
        area = 1.0

        fits = fitsio.FITS(filename)

        self.bkgs = {}

        started = False

        for ext in fits[1: ]:
            extname = ext.get_extname()

            parts = extname.split('_')
            iind = int(parts[0])
            jind = int(parts[1])

            key = iind * 100 + jind

            obkg = Entry.from_fits_ext(ext)

            # Create the refmag bins
            refmagbins = np.arange(obkg.refmagrange[0], obkg.refmagrange[1], refmagbinsize)
            nrefmagbins = refmagbins.size

            # Deal with area if necessary
            if not started:
                if usehdrarea:
                    if 'areas' in obkg.dtype.names:
                        areas = interpol(obkg.areas, obkg.refmagbins, refmagbins)
                    else:
                        hdr = ext.read_header()
                        areas = np.zeros(nrefmagbins) + hdr['AREA']
                else:
                    areas = np.zeros(nrefmagbins) + area

                started = True

            if iind == jind:
                # We are on a diagonal

                # First do the refmag
                ncolbins = obkg.colbins.size
                bc_new = np.zeros((nrefmagbins, ncolbins))
                for i in range(ncolbins):
                    bc_new[:, i] = np.clip(interpol(obkg.bc[:, i], obkg.refmagbins, refmagbins), 0.0, None)

                bc = bc_new.copy()

                # Now do the color
                colbins = np.arange(obkg.colrange[0], obkg.colrange[1], colbinsize)
                ncolbins = colbins.size

                bc_new = np.zeros((nrefmagbins, ncolbins))
                for j in range(nrefmagbins):
                    bc_new[j, :] = np.clip(interpol(bc[j, :], obkg.colbins, colbins), 0.0, None)

                n = np.sum(bc, axis=1, dtype=np.float64) * (colbinsize / obkg.colbinsize)

                sigma_g = bc_new.copy()
                for j in range(ncolbins):
                    sigma_g[:, j] = bc_new[:, j] / areas

                self.bkgs[key] = {'col1': iind,
                                  'col2': jind,
                                  'refmagindex': 1,
                                  'colbins': colbins,
                                  'colrange': obkg.colrange,
                                  'colbinsize': colbinsize,
                                  'refmagbins': refmagbins,
                                  'refmagrange': obkg.refmagrange,
                                  'refmagbinsize': refmagbinsize,
                                  'bc': bc_new,
                                  'n': n,
                                  'sigma_g': sigma_g}
            else:
                # We are on an off-diagonal

                # start with the refmag
                ncol1bins = obkg.col1bins.size
                ncol2bins = obkg.col2bins.size
                bc_new = np.zeros((nrefmagbins, ncol2bins, ncol1bins))
                for i in range(ncol1bins):
                    for j in range(ncol2bins):
                        bc_new[:, j, i] = np.clip(interpol(obkg.bc[:, j, i], obkg.refmagbins, refmagbins), 0.0, None)

                bc = bc_new.copy()

                # color1
                col1bins = np.arange(obkg.col1range[0], obkg.col1range[1], colbinsize)
                ncol1bins = col1bins.size

                bc_new = np.zeros((nrefmagbins, ncol2bins, ncol1bins))
                for j in range(ncol2bins):
                    for k in range(nrefmagbins):
                        bc_new[k, j, :] = np.clip(interpol(bc[k, j, :], obkg.col1bins, col1bins), 0.0, None)

                bc = bc_new.copy()

                col2bins = np.arange(obkg.col2range[0], obkg.col2range[1], colbinsize)
                ncol2bins = col2bins.size

                bc_new = np.zeros((nrefmagbins, ncol2bins, ncol1bins))
                for i in range(ncol1bins):
                    for k in range(nrefmagbins):
                        bc_new[k, :, i] = np.clip(interpol(bc[k, :, i], obkg.col2bins, col2bins), 0.0, None)

                temp = np.sum(bc_new, axis=1, dtype=np.float64) * colbinsize
                n = np.sum(temp, axis=1, dtype=np.float64) * colbinsize

                sigma_g = bc_new.copy()

                for j in range(ncol1bins):
                    for k in range(ncol2bins):
                        sigma_g[:, k, j] = bc_new[:, k, j] / areas

                self.bkgs[key] = {'col1': iind,
                                  'col2': jind,
                                  'refmag_index': 2,
                                  'col1bins': col1bins,
                                  'col1range': obkg.col1range,
                                  'col1binsize': colbinsize,
                                  'col2bins': col2bins,
                                  'col2range': obkg.col2range,
                                  'col2binsize': colbinsize,
                                  'refmagbins': refmagbins,
                                  'refmagrange': obkg.refmagrange,
                                  'refmagbinsize': refmagbinsize,
                                  'bc': bc_new,
                                  'n': n,
                                  'sigma_g': sigma_g}

    def sigma_g_diagonal(self, bkg_index, colors, refmags):
        """
        Compute sigma_g(color, refmag) for a diagonal (single-color) background.

        Parameters
        ----------
        bkg_index: `int`
           Color index along the diagonal
        colors: `np.array`
           Float array of colors
        refmags: `np.array`
           Float array of reference magnitudes

        Returns
        -------
        sigma_g: `np.array`
           Sigma_g(x) for input values
        """
        bkg = self.bkgs[bkg_index * 100 + bkg_index]

        colindex = np.searchsorted(bkg['colbins'], colors - bkg['colbinsize'])
        refmagindex = np.searchsorted(bkg['refmagbins'], refmags - bkg['refmagbinsize'])
        # check for overruns
        badcol, = np.where((colindex < 0) | (colindex >= bkg['colbins'].size))
        colindex[badcol] = 0
        badrefmag, = np.where((refmagindex < 0) | (refmagindex >= bkg['refmagbins'].size))
        refmagindex[badrefmag] = 0

        sigma_g = bkg['sigma_g'][refmagindex, colindex]
        sigma_g[badcol] = np.inf
        sigma_g[badrefmag] = np.inf
        badcombo, = np.where(sigma_g == 0.0)
        sigma_g[badcombo] = np.inf

        return sigma_g

    def lookup_diagonal(self, bkg_index, colors, refmags, doRaise=True):
        """
        Look up the normalized background value b(color, refmag) for a diagonal
        (single-color) background.

        Parameters
        ----------
        bkg_index: `int`
           Color index along the diagonal
        colors: `np.array`
           Float array of colors
        refmags: `np.array`
           Float array of reference magnitudes
        doRaise: `bool`, optional
           Raise exception if color background does not exist.  Default is True.

        Returns
        -------
        b: `np.array`
           Background(x) for input values
        """
        try:
            bkg = self.bkgs[bkg_index * 100 + bkg_index]
        except KeyError:
            if doRaise:
                raise KeyError("Could not find a color background for %d" % (bkg_index))
            else:
                return np.zeros(colors.size)

        refmagindex = np.clip(np.searchsorted(bkg['refmagbins'], refmags - bkg['refmagbinsize']), 0, bkg['refmagbins'].size - 1)
        col_index = np.clip(np.searchsorted(bkg['colbins'], colors - bkg['colbinsize']), 0, bkg['colbins'].size - 1)

        return bkg['bc'][refmagindex, col_index] / bkg['n'][refmagindex]

    def get_colrange(self, bkg_index):
        """
        Get the background color range for a given background index.

        Parameters
        ----------
        bkg_index: `int`
           Color index along the diagonal

        Returns
        -------
        range: `np.array`
           2-element float array with color min, max
        """
        bkg = self.bkgs[bkg_index * 100 + bkg_index]
        return bkg['colrange']

    def lookup_offdiag(self, bkg_index1, bkg_index2, colors1, colors2, refmags, doRaise=True):
        """
        Look up the normalized background value b(color1, color2, refmag) for
        an off-diagonal (two-color) background.

        Parameters
        ----------
        bkg_index1: `int`
           Color index for color 1
        bkg_index2: `int`
           Color index for color 2
        colors1: `np.array`
           Float array of colors (1)
        colors2: `np.array`
           Float array of colors (2)
        refmags: `np.array`
           Float array of reference mags
        doRaise: `bool`, optional
           Raise exception if color background does not exist.  Default is True.
        """
        key = bkg_index1 * 100 + bkg_index2
        if key not in self.bkgs:
            key = bkg_index2 * 100 + bkg_index1

        try:
            bkg = self.bkgs[key]
        except KeyError:
            if doRaise:
                raise KeyError("Could not find a color background for %d, %d" % (bkg_index1, bkg_index1))
            else:
                return np.zeros(colors1.size)

        refmagindex = np.clip(np.searchsorted(bkg['refmagbins'], refmags - bkg['refmagbinsize']), 0, bkg['refmagbins'].size - 1)
        col_index1 = np.clip(np.searchsorted(bkg['col1bins'], colors1 - bkg['col1binsize']), 0, bkg['col1bins'].size - 1)
        col_index2 = np.clip(np.searchsorted(bkg['col2bins'], colors2 - bkg['col2binsize']), 0, bkg['col2bins'].size - 1)

        return bkg['bc'][refmagindex, col_index2, col_index1] / bkg['n'][refmagindex] class GalaxyCatalog(Catalog):
    """
    Class to describe a redmapper galaxy Catalog.

    This contains ways of reading galaxy catalogs, getting colors, and matching
    galaxies to clusters.
    """
    entry_class = Galaxy

    def __init__(self, *arrays, **kwargs):
        """
        Instantiate a GalaxyCatalog

        Parameters
        ----------
        *arrays: `np.ndarray`
           Parameters for arrays to build galaxy catalog (e.g. galaxy and zred
           information)
        depth: `int`, optional
           HTM matcher depth, default is 10.
        """
        super(GalaxyCatalog, self).__init__(*arrays)
        self._htm_matcher = None
        self.depth = 10 if 'depth' not in kwargs else kwargs['depth']

    @classmethod
    def from_galfile(cls, filename, zredfile=None, nside=0, hpix=[], border=0.0, truth=False,
                     use_tempfile=False, refmag_range=[-1000.0, 1000.0], chisq_max=1e100):
        """
        Generate a GalaxyCatalog from a redmapper "galfile."

        Parameters
        ----------
        filename: `str`
           Filename of the redmapper "galfile" galaxy file.
           This file may be a straight fits file or (recommended) a
           galaxy table "master_table.fit" summary file.
        zredfile: `str`, optional
           Filename of the redmapper zred "zreds_master_table.fit" summary file.
        nside: `int`, optional
           Nside of healpix sub-region to read in.  Default is 0 (full catalog).
        hpix: `list`, optional
           Healpix numbers (ring format) of sub-region to read in.
           Default is [] (full catalog).
        border: `float`, optional
           Border around hpix (in degrees) to read in.  Default is 0.0.
        truth: `bool`, optional
           Read in truth information if available (e.g. mocks)?  Default is False.
        use_tempfile : `bool`, optional
            Use a tempfile to store galaxies to conserve memory.
        refmag_range : `list`, optional
            Refmag min, max to cut if using tempfile.
        chisq_max : `float`, optional
            Maximum chisq to cut if reading zreds and if using tempfile.
        """
        if zredfile is not None:
            use_zred = True
        else:
            use_zred = False

        if not isinstance(hpix, Iterable):
            _hpix = [hpix]
        else:
            _hpix = hpix

        # do we have appropriate keywords
        if len(_hpix) > 0 and nside is None:
            raise ValueError("If hpix is specified, must also specify nside")
        if border < 0.0:
            raise ValueError("Border must be >= 0.0.")
        # ensure that nside is valid, and hpix is within range (if necessary)
        if nside > 0:
            if not hp.isnsideok(nside):
                raise ValueError("Nside not valid")
            if len(_hpix) > 0:
                for _hp in _hpix:
                    if _hp < 0 or _hp >= hp.nside2npix(nside):
                        raise ValueError("hpix %d is out of range." % (hp))

        if border > 0.0 and len(_hpix) > 0:
            if len(_hpix) != 1:
                raise NotImplementedError("Cannot read a boundary around a pixel list.")

        # check that the file is there and the right format
        # this will raise an exception if it's not there.
        hdr = fitsio.read_header(filename, ext=1)
        pixelated = hdr.get("PIXELS", 0)
        fitsformat = hdr.get("FITS", 0)

        # check zredfile
        if use_zred:
            zhdr = fitsio.read_header(zredfile, ext=1)
            zpixelated = zhdr.get("PIXELS", 0)

        if not pixelated:
            cat = fitsio.read(filename, ext=1, upper=True)
            if use_zred:
                zcat = fitsio.read(zredfile, ext=1, upper=True)
                if zcat.size != cat.size:
                    raise ValueError("zredfile is a different length (%d) than catfile (%d)" % (zcat.size, cat.size))
                return cls(cat, zcat)
            else:
                return cls(cat)
        else:
            if use_zred:
                if not zpixelated:
                    raise ValueError("galfile is pixelated but zredfile is not")

        # this is to keep us from trying to use old IDL galfiles
        if not fitsformat:
            raise ValueError("Input galfile must describe fits files.")

        # now we can read in the galaxy table summary file...
        tab = Entry.from_fits_file(filename, ext=1)
        nside_tab = tab.nside
        if nside > nside_tab:
            raise ValueError("""Requested nside (%d) must not be larger than
                                    table nside (%d).""" % (nside, nside_tab))

        if use_zred:
            ztab = Entry.from_fits_file(zredfile, ext=1)
            zpath = os.path.dirname(zredfile)

        # which files do we want to read?
        path = os.path.dirname(os.path.abspath(filename))

        indices = get_subpixel_indices(tab, hpix=_hpix, border=border, nside=nside)

        # Make sure all the zred files are there
        if use_zred:
            # The default mode, copied from the IDL code, is that we just don't
            # read in any galaxy pixels that don't have an associated zred.
            # I don't know if this is what we want going forward, but I'll leave
            # it like this at the moment.
            # Also, we are assuming that the files actually match up in terms of length, etc.
            mark = np.zeros(indices.size, dtype=bool)
            for i, f in enumerate(ztab.filenames[indices]):
                try:
                    fname = os.path.join(zpath, f.decode())
                except AttributeError:
                    fname = os.path.join(zpath, f)
                if os.path.isfile(fname):
                    mark[i] = True

            bad, = np.where(~mark)
            if bad.size == indices.size:
                raise ValueError("There are no zred files associated with the galaxy pixels.")

            indices = np.delete(indices, bad)

        # Prepare border cuts if necessary
        trim_border = False
        if len(_hpix) == 1 and nside > 0 and border > 0.0:
            trim_border = True
            nside_cutref = 512
            boundaries = hp.boundaries(nside, _hpix[0], step=nside_cutref//nside)

            # Need to get the sub-pixels, use nest and bit-shifting.
            bit_shift = 2*int(np.round(np.log2(nside_cutref/nside)))
            inhpix_nest = np.arange(2**bit_shift, dtype=np.int32) + np.left_shift(hp.ring2nest(nside, _hpix[0]), bit_shift)
            inhpix = hp.nest2ring(nside_cutref, inhpix_nest)

            for i in range(boundaries.shape[1]):
                pixint = hp.query_disc(nside_cutref, boundaries[:, i], np.radians(border), inclusive=True, fact=8)
                inhpix = np.append(inhpix, pixint)
            inhpix = np.unique(inhpix)

        # create the catalog array to read into

        # and need to be able to cut?
        try:
            first_fname = os.path.join(path, tab.filenames[indices[0]].decode())
        except AttributeError:
            first_fname = os.path.join(path, tab.filenames[indices[0]])
        elt = fitsio.read(first_fname, ext=1, rows=0, lower=True)
        dtype_in = elt.dtype.descr
        if not truth:
            mark = []
            for dt in dtype_in:
                if (dt[0] != 'ztrue' and dt[0] != 'm200' and dt[0] != 'central' and
                    dt[0] != 'halo_id'):
                    mark.append(True)
                else:
                    mark.append(False)

            dtype = [dt for i, dt in enumerate(dtype_in) if mark[i]]
            columns = [dt[0] for dt in dtype]
        else:
            dtype = dtype_in
            columns = None

        cat_fields = [dt[0] for dt in dtype]

        # New more memory efficient reading...
        ngals = np.sum(tab.ngals[indices])

        if use_zred:
            try:
                fname = os.path.join(zpath, ztab.filenames[indices[0]].decode())
            except AttributeError:
                fname = os.path.join(zpath, ztab.filenames[indices[0]])

            zelt = fitsio.read(fname, ext=1, rows=0, lower=True)
            zcat_fields = [dt[0] for dt in zelt.dtype.descr]

            dtype.extend(zelt.dtype.descr)

        if use_tempfile:
            # Start the temporary file
            fd, tempFile = tempfile.mkstemp(suffix='.fits')
            os.close(fd)

            tempfits = fitsio.FITS(tempFile, mode='rw', clobber=True)
            tempfits.create_table_hdu(dtype=dtype)
        else:
            # Make the regular catalog
            cat = np.zeros(np.sum(tab.ngals[indices]), dtype=dtype)

        # read the files
        ctr = 0
        for index in indices:
            try:
                fname = os.path.join(path, tab.filenames[index].decode())
            except AttributeError:
                fname = os.path.join(path, tab.filenames[index])

            if use_tempfile:
                # Read into temporary buffer
                tempcat = np.zeros(tab.ngals[index], dtype=dtype)
                tempcat[cat_fields][:] = fitsio.read(fname, ext=1, lower=True, columns=columns)
            else:
                cat[cat_fields][ctr: ctr + tab.ngals[index]] = fitsio.read(fname, ext=1, lower=True, columns=columns)

            if use_zred:
                # Note that this effectively checks that the numbers of rows in each file match properly (though the exception will be cryptic...)
                try:
                    fname = os.path.join(zpath, ztab.filenames[index].decode())
                except AttributeError:
                    fname = os.path.join(zpath, ztab.filenames[index])

                if use_tempfile:
                    # Read into temporary buffer
                    tempcat[zcat_fields][:] = fitsio.read(fname, ext=1, lower=True)
                else:
                    cat[zcat_fields][ctr: ctr + tab.ngals[index]] = fitsio.read(fname, ext=1, lower=True)

            if use_tempfile:
                # Cut down the tempcat based on refmag low/high and chisq_max
                guse = ((tempcat['refmag'] > refmag_range[0]) & (tempcat['refmag'] < refmag_range[1]))
                if use_zred:
                    guse &= (tempcat['chisq'] < chisq_max)

                if guse.sum() > 0:
                    # cut down the tempcat if necessary and spool to tempfile
                    if trim_border:
                        ipring = hp.ang2pix(nside_cutref, tempcat['ra'], tempcat['dec'], lonlat=True)
                        _, indices = esutil.numpy_util.match(inhpix, ipring[guse])

                        tempfits[1].append(tempcat[guse][indices])
                    else:
                        tempfits[1].append(tempcat[guse])

            ctr += tab.ngals[index]

        if use_tempfile:
            # close out the tempfile
            tempfits.close()

            # read in the tempfile
            cat = fitsio.read(tempFile, ext=1)

            # delete the tempfile
            os.remove(tempFile)

            return(cls(cat))
        else:
            if trim_border:
                ipring = hp.ang2pix(nside_cutref, cat['ra'], cat['dec'], lonlat=True)
                _, indices = esutil.numpy_util.match(inhpix, ipring)

                return cls(cat[indices])
            else:
                # No cuts
                return cls(cat)

    @property
    def galcol(self):
        """
        Get the array of galaxy colors.

        Returns
        -------
        galcol: `np.array`
           Float array of colors [ngal, nmag - 1]
        """
        galcol = self.mag[:,:-1] - self.mag[:,1:]
        return galcol

    @property
    def galcol_err(self):
        """
        Get the array of galaxy color errors.

        Returns
        -------
        galcol_err: `np.array`
           Float array of errors [ngal, nmag - 1]
        """
        galcol_err = np.sqrt(self.mag_err[:, :-1]**2. + self.mag_err[:, 1:]**2.)
        return galcol_err

    def add_zred_fields(self, nsamp):
        """
        Add default columns for storing zreds.

        Note that this will do nothing if the columns are already there.

        Modifies GalaxyCatalog in place.

        Parameters
        ----------
        nsamp: `int`
           Number of samples of zred to record
        """
        dtype_augment = [dt for dt in zred_extra_dtype(nsamp) if dt[0].lower() not in self.dtype.names]
        if len(dtype_augment) > 0:
            self.add_fields(dtype_augment)

    def match_one(self, ra, dec, radius):
        """
        Match one ra/dec position to the galaxy catalog.

        This is typically used when finding all the neighbors around a cluster
        over a relatively large area.

        Parameters
        ----------
        ra: `float`
           Right ascension to match to.
        dec: `float`
           Declination to match to.
        radius: `float`
           Match radius (degrees)

        Returns
        -------
        indices: `np.array`
           Integer array of GalaxyCatalog indices that match to input ra/dec
        dists: `np.array`
           Float array of distance (degrees) from each galaxy in indices
        """
        if self._htm_matcher is None:
            self._htm_matcher = Matcher(self.depth, self.ra, self.dec)

        _, indices, dists = self._htm_matcher.match(ra, dec, radius, maxmatch=0)

        return indices, dists

    def match_many(self, ras, decs, radius, maxmatch=0):
        """
        Match many ra/dec positions to the galaxy catalog.

        This is typically used when matching a galaxy catalog to another
        catalog (e.g. a spectroscopic catalog).  Running this with a large
        number of positions with a large radii is possible to run out of
        memory!

        Parameters
        ----------
        ras: `np.array`
           Float arrays of right ascensions to match to.
        decs: `np.array`
           Float arrays of declinations to match to.
        radius: `np.array` or `float`
           Float array or float match radius in degrees.
        maxmatch: `int`
           Maximum number of galaxy matches to each ra/dec.
           Default is 0 (no maximum)

        Returns
        -------
        i0: `np.array`
           Integer array of indices matched to input ra/dec
        i1: `np.array`
           Integer array of indices matched to galaxy ra/dec
        dists: `np.array`
           Float array of match distances for each i0/i1 pair (degrees).
        """
        if self._htm_matcher is None:
            self._htm_matcher = Matcher(self.depth, self.ra, self.dec)

        return self._htm_matcher.match(ras, decs, radius, maxmatch=maxmatch) class Catalog(DataObject):
    """
    A Catalog is an extension of a DataObject.  It can be decomposed into
    individual Entry objects.  This class is used to describe catalogs of all
    sorts, including galaxy catalogs, cluster catalogs, and other ndarrays.
    """

    entry_class = Entry

    @property
    def size(self):
        """
        Return the size of the Catalog.
        """
        return self._ndarray.size

    def append(self, append_cat):
        """
        Append a number of rows to the catalog, in-place.

        Parameters
        ----------
        append_cat: `redmapper.Catalog` or `np.ndarray`
           Catalog to append
        """
        if isinstance(append_cat, Catalog):
            self._ndarray = np.append(self._ndarray, append_cat._ndarray)
        else:
            self._ndarray = np.append(self._ndarray, append_cat)

    def extend(self, n_new):
        """
        Extend catalog with zero-filled rows, in-place.

        Parameters
        ----------
        n_new: `int`
           Number of new rows to append
        """
        temp = np.zeros(n_new, dtype=self._ndarray.dtype)
        self._ndarray = np.append(self._ndarray, temp)

    def __len__(self): return self.size

    def __getitem__(self, key):
        if isinstance(key, int):
            return self.entry_class(self._ndarray.__getitem__(key))
        return type(self)(self._ndarray.__getitem__(key))

    def __setitem__(self, key, val):
        self._ndarray.__setitem__(key, val) class Entry(DataObject):
    """
    An Entry is an extension of a DataObject.  It is intended to be used as a
    single Entry in a Catalog object.
    """
    def __init__(self, array):
        """
        Instantiate an Entry.

        Parameters
        ----------
        array: `np.ndarray` of length 1
        """
        if array.size != 1:
            raise ValueError("Input array must have length one.")
        # If this is an array of length 1, we want it to be a scalar-ish
        if len(array.shape) == 0:
            super(Entry, self).__init__(array)
        else:
            super(Entry, self).__init__(array[0])

    # @classmethod
    # def from_dict(cls, dict): pass

    def add_fields(self, newdtype):
        """
        Add fields to an Entry.

        Parameters
        ----------
        newdtype: data-type
           `np.dtype` description
        """
        array = np.zeros(self._ndarray.size, newdtype)
        self._lower_array(array)
        self._ndarray = self._merge_arrays([self._ndarray, array])[0]

    def __getattr__(self, attr):
        try:
            return self._ndarray[attr.lower()]
        except:
            return object.__getattribute__(self, attr) class ZredColor(object):
    """
    Class to compute zred galaxy photometric redshift using a color-based
    red-sequence parametrization.
    """

    def __init__(self, zredstr, sigint=0.001, do_correction=True,
                 use_photoerr=True, zrange=None):
        """
        Instantiate a ZredColor object.

        Parameters
        ----------
        zredstr: `redmapper.RedSequenceColorPar`
           Red sequence parametrization
        sigint: `float`, optional
           Intrinsic scatter floor.  Default is 0.001.
        do_correction: `bool`, optional
           Apply zred correction terms?  Default is True.
        use_photoerr: `bool`, optional
           Use photometric errors in computing zred?  Default is True.
        zrange: `list`, optional
           Redshift range.  Useful for testing.  Default is None (use
           zredstr redshift range).
        """
        self.zredstr = zredstr

        self.sigint = sigint
        self.do_correction = do_correction
        self.use_photoerr = use_photoerr
        self.zrange = zrange

        self.nz = self.zredstr.z.size - 1
        self.notextrap, = np.where(~self.zredstr.extrapolated)

        if self.zrange is None:
            self.zbinstart = 0
            self.zbinstop = self.nz - 1
        else:
            u, = np.where((self.zredstr.z > self.zrange[0]) &
                          (self.zredstr.z < self.zrange[1]))
            self.zbinstart = u[0]
            self.zbinstop = u[-1]

    def compute_zreds(self, galaxies):
        """
        Compute zreds for a catalog of galaxies.

        Will set galaxies.zred, galaxies.zred_e, etc.

        Parameters
        ----------
        galaxies: `redmapper.GalaxyCatalog`
           Catalog of galaxies to compute zred.
        """

        for galaxy in galaxies:
            self.compute_zred(galaxy, no_corrections=True)

        if self.do_correction:
            # Bulk processing
            olddzs = np.zeros(galaxies.size)
            dzs = np.zeros_like(olddzs)
            iteration = 0

            pivotmags = interpol(self.zredstr.pivotmag, self.zredstr.z, galaxies.zred_uncorr)

            while (iteration < 5):
                olddzs[:] = dzs
                dzs[:] = interpol(self.zredstr.corr, self.zredstr.z, galaxies.zred_uncorr + olddzs) + (galaxies.refmag - pivotmags) * interpol(self.zredstr.corr_slope, self.zredstr.z, galaxies.zred_uncorr + olddzs)
                iteration += 1

            galaxies.zred = galaxies.zred_uncorr + dzs
            galaxies.zred_e = galaxies.zred_uncorr_e * interpol(self.zredstr.corr_r, self.zredstr.z, galaxies.zred)

            dz2s = interpol(self.zredstr.corr2, self.zredstr.z, galaxies.zred_uncorr)
            r2s = interpol(self.zredstr.corr2_r, self.zredstr.z, galaxies.zred_uncorr)

            galaxies.zred2 = galaxies.zred_uncorr + dz2s
            galaxies.zred2_e = galaxies.zred_uncorr_e * r2s


    def compute_zred(self, galaxy, no_corrections=False):
        """
        Compute zred for a single galaxy.

        Will set galaxy.zred, galaxy.zred_e, etc.

        Parameters
        ----------
        galaxy: `redmapper.Galaxy`
           Galaxy to compute zred
        no_corrections: `bool`, optional
           Do not apply redshift corrections.  Default is False.
        """

        lndist = np.zeros(self.nz) - 1e12
        chisq = np.zeros(self.nz) + 1e12

        zbins_limited, = np.where((galaxy.refmag < self.zredstr.maxrefmag) &
                                  (galaxy.refmag > self.zredstr.minrefmag) &
                                  (self.zredstr.z < 100.0))

        if zbins_limited.size < 2:
            self._reset_bad_values(galaxy)
            return

        neighbors = 10
        zbins = np.arange(*np.clip([zbins_limited[0] - neighbors,
                                    zbins_limited[-1] + neighbors],
                                   0, self.nz))

        lndist[zbins], chisq[zbins] = self._calculate_lndist(galaxy, zbins)

        # move from log space to regular space
        maxlndist = np.max(lndist[zbins])
        dist = np.zeros_like(lndist)
        with np.errstate(invalid='ignore', over='ignore'):
            dist[zbins] = np.exp(lndist[zbins] - maxlndist)

        # fix infinities and NaNs
        bad, = np.where(~np.isfinite(dist))
        dist[bad] = 0.0

        # take the maximum where not extrapolated
        ind_temp = np.argmax(dist[self.notextrap])
        ind = self.notextrap[ind_temp]

        calcinds, = np.where(dist > 1e-5)

        if calcinds.size >= 3:
            tdist = scipy.integrate.trapz(dist[calcinds], self.zredstr.z[calcinds])
            zred_temp = scipy.integrate.trapz(dist[calcinds] * self.zredstr.z[calcinds],
                                              self.zredstr.z[calcinds]) / tdist
            zred_e = scipy.integrate.trapz(dist[calcinds] * self.zredstr.z[calcinds]**2.,
                                           self.zredstr.z[calcinds]) / tdist - zred_temp**2.
        else:
            tdist = np.sum(dist[calcinds])
            zred_temp = np.sum(dist[calcinds] * self.zredstr.z[calcinds]) / tdist
            zred_e = np.sum(dist[calcinds] * self.zredstr.z[calcinds]**2.) / tdist - zred_temp**2.

        if zred_e < 0.0:
            zred_e = 1.0
        else:
            zred_e = np.sqrt(zred_e)

        zred_e = zred_e if zred_e > 0.005 else 0.005

        # Now fit a parabola to get the perfect zred
        ind_temp = np.argmax(dist[self.notextrap])
        ind = self.notextrap[ind_temp]

        zred = zred_temp.copy()

        neighbors = 2
        use, = np.where(lndist > -1e10)
        if use.size >= neighbors * 2 + 1:
            # If it hits a wall, then move in the other direction to ensure we have at least neighbors*2+1 points
            minuse = use.min()
            maxuse = use.max()

            minindex = minuse if minuse > ind - neighbors else ind - neighbors
            maxindex = maxuse if maxuse < ind + neighbors else ind + neighbors
            if minindex == minuse:
                maxindex = np.clip(minuse + 2 * neighbors, None, maxuse)
            elif maxindex == maxuse:
                minindex = np.clip(maxuse - 1 - 2*neighbors, minuse, None)

            if ((maxindex - minindex + 1) >= 5):
                X = np.zeros((maxindex - minindex + 1, 3))
                X[:, 1] = self.zredstr.z[minindex:maxindex + 1]
                X[:, 0] = X[:, 1] * X[:, 1]
                X[:, 2] = 1
                y = lndist[minindex: maxindex + 1]

                fit = np.matmul(np.matmul(np.linalg.inv(np.matmul(X.T, X)), X.T), y)

                if fit[0] < 0.0:
                    ztry = -fit[1] / (2.0 * fit[0])
                    # Don't let it move too far, or it's a bad fit
                    if (np.abs(ztry - zred) < 2.0*zred_e):
                        zred = ztry

        # And compute values at the real zred peak
        x = (self.zredstr.z - zred) / zred_e
        newdist = np.exp(-0.5 * x * x)

        bad, = np.where((lndist < -1e10) | (~np.isfinite(lndist)))
        newdist[bad] = 0.0
        lndist[bad] = -1e11

        if calcinds.size >= 3:
            # Note there maybe should be a distcorr here, but this is not
            #  actually computed in the IDL code (bug?)
            lkhd = scipy.integrate.trapz(newdist[calcinds] * (lndist[calcinds]), self.zredstr.z[calcinds]) / scipy.integrate.trapz(newdist[calcinds], self.zredstr.z[calcinds])
        else:
            lkhd = np.sum(newdist[calcinds] * lndist[calcinds]) / np.sum(newdist[calcinds])

        # Get chisq at the closest bin position
        zbin = np.argmin(np.abs(zred - self.zredstr.z))
        chisq = chisq[zbin]

        if not np.isfinite(lkhd):
            self._reset_bad_values(galaxy)
            return

        # And sample the uncorrected p(z)
        gdzbins, = np.where((dist > 1e-10) & (np.isfinite(dist)))

        if gdzbins.size < 3:
            # We cannot do a proper p(z)
            zred_samp = np.zeros(galaxy.zred_samp.size) + zred
        else:
            pz = dist.copy()
            n = scipy.integrate.simps(pz[gdzbins], self.zredstr.z[gdzbins])
            pz /= n

            pdf = scipy.interpolate.interp1d(self.zredstr.z[gdzbins], pz[gdzbins], kind='quadratic',
                                             bounds_error=False, fill_value=0.0)
            zred_samp = sample_from_pdf(pdf,
                                        [self.zredstr.z[gdzbins[0]], self.zredstr.z[gdzbins[-1]]],
                                        0.0001,
                                        galaxy.zred_samp.size)

        # And apply the corrections
        zred2 = np.zeros(1) + zred
        zred2_e = np.zeros(1) + zred_e
        zred_uncorr = np.zeros(1) + zred
        zred_uncorr_e = np.zeros(1) + zred_e

        if self.do_correction and not no_corrections:
            olddz = -1.0
            dz = 0.0
            iteration = 0

            pivotmag = interpol(self.zredstr.pivotmag, self.zredstr.z, zred)

            while np.abs(olddz - dz) > 1e-3 and iteration < 10:
                olddz = copy.copy(dz)
                dz = interpol(self.zredstr.corr, self.zredstr.z, zred + olddz) + (galaxy.refmag - pivotmag) * interpol(self.zredstr.corr_slope, self.zredstr.z, zred + olddz)
                iteration += 1

            zred = zred + dz

            # evaluate error correction at "z_true"
            zred_e *= interpol(self.zredstr.corr_r, self.zredstr.z, zred)

            # And the zred2 correction
            dz = interpol(self.zredstr.corr2, self.zredstr.z, zred2) + (galaxy.refmag - pivotmag) * (interpol(self.zredstr.corr2_slope, self.zredstr.z, zred2))
            # this is evaluated at zred0
            r2 = interpol(self.zredstr.corr2_r, self.zredstr.z, zred2)

            zred2 += dz
            zred2_e *= r2

        # Finally store the values

        galaxy.zred = zred
        galaxy.zred_e = zred_e
        galaxy.zred2 = zred2
        galaxy.zred2_e = zred2_e
        galaxy.zred_uncorr = zred_uncorr
        galaxy.zred_uncorr_e = zred_uncorr_e
        galaxy.zred_samp = zred_samp
        galaxy.chisq = chisq
        galaxy.lkhd = lkhd

        # and we're done

    def _calculate_lndist(self, galaxy, zbins):
        """
        Calculate the log-likelihood for a list of redshift bins.

        Parameters
        ----------
        galaxy: `redmapper.Galaxy`
           Galaxy to compute likelihoods
        zbins: `np.array`
           Integer array of redshift bins

        Returns
        -------
        lndist: `np.array`
           Log-likelihood for the redshift bins
        chisq: `np.array`
           Fit chi-squared for the redshift bins
        """

        # Note we need to deal with photoerr...
        if zbins.size > 1:
            # we have many bins...
            chisq = self.zredstr.calculate_chisq_redshifts(galaxy, zbins, z_is_index=True, calc_lkhd=False)
        else:
            # we have a single bin... hack this
            chisq = self.zredstr.calculate_chisq(galaxy, np.array([zbins[0], zbins[0]]), z_is_index=True, calc_lkhd=False)[0]

        lndist = -0.5 * chisq

        #with np.errstate(invalid='ignore'):
        lndistcorr = np.log((10.**(0.4 * (self.zredstr.alpha + 1.0) *
                                   (self.zredstr._mstar[zbins] - galaxy.refmag)) *
                             np.exp(-10.**(0.4 * (self.zredstr._mstar[zbins] - galaxy.refmag)))) *
                            self.zredstr.volume_factor[zbins])

        lndist += lndistcorr

        bad, = np.where(~np.isfinite(lndist))
        lndist[bad] = -1e11

        return (lndist, chisq)

    def _reset_bad_values(self, galaxy):
        """
        Reset all galaxy zred values to "bad" (-1s)

        Parameters
        ----------
        galaxy: `redmapper.Galaxy`
           Galaxy to set zred values.
        """

        galaxy.lkhd = -1000.0
        galaxy.zred = -1.0
        galaxy.zred_e = -1.0
        galaxy.zred2 = -1.0
        galaxy.zred2_e = -1.0
        galaxy.zred_uncorr = -1.0
        galaxy.zred_uncorr_e = -1.0
        galaxy.chisq = -1.0 def make_nodes(zrange, nodesize, maxnode=None):
    """
    Make spline nodes over a given range, mostly but not entirely uniform.

    This will return nodes that fill the space desired, and will make the
    highest redshift node spacing larger to fill the space.  In addition, the
    maxnode maximum useful node can be specified, and this will create one
    final node at the high redshift end.  This is useful when you have a color
    (e.g. u-g) that is only fittable over a restricted redshift range.

    Parameters
    ----------
    zrange: `np.array` or `list`
       2-element float array with redshift range
    nodesize: `float`
       Default (and minimum) node spacing
    maxnode: `float`, optional
       Maximum useful node.  Default is None (do full redshift range)

    Returns
    -------
    nodes: `np.array`
       Float array of redshift nodes.
    """

    if maxnode is None or maxnode < 0.0:
        _maxnode = zrange[1]
    else:
        _maxnode = np.clip(maxnode, 0.0, zrange[1])

    # Start with a simple arange
    nodes = np.arange(zrange[0], _maxnode, nodesize)

    # Should we extend the last bin?
    if ((_maxnode - nodes.max()) > (nodesize / 2. + 0.01)):
        # one more node!
        nodes = np.append(nodes, _maxnode)
    elif ~np.allclose(np.max(nodes), _maxnode):
        nodes[-1] = _maxnode

    # and finally check if maxnode was lower
    if _maxnode < zrange[1]:
        nodes = np.append(nodes, zrange[1])

    return nodes class CubicSpline(object):
    """
    CubicSpline interpolation class.
    """
    def __init__(self, x, y, yp=None, fixextrap=False):
        """
        Instantiate a CubicSpline object.

        Parameters
        ----------
        x: `np.array`
           Float array of node positions
        y: `np.array`
           Float array of node values
        yp: `str`
           Type of spline.  Default is None, which is "natural"
        fixextrap: `bool`, optional
           Fix the extrapolation at the end of the node positions.
           Default is False.
        """
        npts = len(x)
        mat = np.zeros((3, npts))
        # enforce continuity of 1st derivatives
        mat[1,1:-1] = (x[2:  ]-x[0:-2])/3.
        mat[2,0:-2] = (x[1:-1]-x[0:-2])/6.
        mat[0,2:  ] = (x[2:  ]-x[1:-1])/6.
        bb = np.zeros(npts)
        bb[1:-1] = ((y[2:  ]-y[1:-1])/(x[2:  ]-x[1:-1]) -
                    (y[1:-1]-y[0:-2])/(x[1:-1]-x[0:-2]))
        if yp is None: # natural cubic spline
            mat[1,0] = 1.
            mat[1,-1] = 1.
            bb[0] = 0.
            bb[-1] = 0.
        elif yp == '3d=0':
            mat[1, 0] = -1./(x[1]-x[0])
            mat[0, 1] =  1./(x[1]-x[0])
            mat[1,-1] =  1./(x[-2]-x[-1])
            mat[2,-2] = -1./(x[-2]-x[-1])
            bb[ 0] = 0.
            bb[-1] = 0.
        else:
            mat[1, 0] = -1./3.*(x[1]-x[0])
            mat[0, 1] = -1./6.*(x[1]-x[0])
            mat[2,-2] =  1./6.*(x[-1]-x[-2])
            mat[1,-1] =  1./3.*(x[-1]-x[-2])
            bb[ 0] = yp[0]-1.*(y[ 1]-y[ 0])/(x[ 1]-x[ 0])
            bb[-1] = yp[1]-1.*(y[-1]-y[-2])/(x[-1]-x[-2])
        y2 = solve_banded((1,1), mat, bb)
        self.x, self.y, self.y2 = (x, y, y2)

        self.fixextrap = fixextrap

    def splint(self,x):
        """
        Compute spline interpolation.

        Parameters
        ----------
        x: `np.array`
           Float array of x values to compute interpolation

        Returns
        -------
        y: `np.array`
           Spline interpolated values at x
        """
        npts = len(self.x)
        lo = np.searchsorted(self.x, x)-1
        lo = np.clip(lo, 0, npts-2)
        hi = lo + 1
        dx = self.x[hi] - self.x[lo]
        a = (self.x[hi] - x)/dx
        b = (x-self.x[lo])/dx
        y = (a*self.y[lo]+b*self.y[hi]+
             ((a**3-a)*self.y2[lo]+(b**3-b)*self.y2[hi])*dx**2./6.)
        return y

    def __call__(self, x):
        """
        Compute spline interpolation.

        Parameters
        ----------
        x: `np.array`
           Float array of x values to compute interpolation

        Returns
        -------
        y: `np.array`
           Spline interpolated values at x
        """

        if not self.fixextrap:
            return self.splint(x)
        else:
            vals = self.splint(x)

            lo, = np.where(x < self.x[0])
            vals[lo] = self.y[0]

            hi, = np.where(x > self.x[-1])
            vals[hi] = self.y[-1]

            return vals def interpol(v, x, xout):
    """
    Port of IDL interpol.py.  Does fast and simple linear interpolation.

    Parameters
    ----------
    v: `np.array`
       Float array of y (dependent) values to interpolate between
    x: `np.array`
       Float array of x (independent) values to interpolate between
    xout: `np.array`
       Float array of x values to compute interpolated values

    Returns
    -------
    yout: `np.array`
       Float array of y output values associated with xout
    """

    m = v.size
    nOut = m

    s = np.clip(np.searchsorted(x, xout) - 1, 0, m - 2)

    diff = v[s + 1] - v[s]

    return (xout - x[s]) * diff / (x[s + 1] - x[s]) + v[s] class RedGalInitialColors(object):
    """
    Class to describe the RedGalInitialColors(z, color) look-up table.
    """

    def __init__(self, redgal_template):
        """
        Instantiate a RedGalInitialColors object.

        Parameters
        ----------
        redgal_template: `str`
           Name of redgal_template file to get look-up table.
        """

        from . import Catalog
        from pkg_resources import resource_filename
        from pkg_resources import resource_exists

        self._template_file = None

        module = __name__.split('.')[0]
        if resource_exists(module, 'data/initcolors/%s' % (redgal_template)):
            self._template_file = resource_filename(module, 'data/initcolors/%s' % (redgal_template))
        elif os.path.isfile(redgal_template):
            self._template_file = os.path.abspath(redgal_template)
        else:
            raise IOError("Could not find redgal_template file %s in resource or path." % (redgal_template))

        self._template = Catalog.from_fits_file(self._template_file, ext=1)
        template_hdr = fitsio.read_header(self._template_file, ext=1)

        if ',' in template_hdr['BANDS']:
            self._template_bands = template_hdr['BANDS'].rstrip().split(',')
        else:
            self._template_bands = list(template_hdr['BANDS'].rstrip())

        self._color_names = []
        for i in range(len(self._template_bands) - 1):
            self._color_names.append(self._color_name(self._template_bands[i],
                                                      self._template_bands[i + 1]))

    @property
    def zmax(self):
        return self._template.z.max()

    def __call__(self, band1, band2, z):
        """
        Return the initial color guess (band1 - band2) at redshift z.

        Parameters
        ----------
        band1: `str`
           The first band in the color
        band2: `str`
           The second band in the color
        z: `np.array`
           Float array of redshifts

        Returns
        -------
        color: `np.array`
           Initial red galaxy color at redshifts z
        """

        name = self._color_name(band1, band2)
        if name not in self._color_names:
            raise ValueError("Color %s-%s not in initial template file %s." % (band1, band2,
                                                                               self._template_file))

        index = self._color_names.index(name)
        return interpol(self._template.color[:, index], self._template.z, z)

    def _color_name(self, band1, band2):
        """
        Return the color name associated with two bands.

        Parameters
        ----------
        band1: `str`
           The first band in the color
        band2: `str`
           The second band in the color

        Returns
        -------
        colorname: `str`
           Name of the color to use as a key.
        """

        name = '%s-%s' % (band1.upper(), band2.upper())
        return name
import os
import numpy as np
import fitsio
import time
import esutil
        import matplotlib.pyplot as plt
from scipy.optimize import least_squares
from ..configuration import Configuration
from ..fitters import MedZFitter, RedSequenceFitter, RedSequenceOffDiagonalFitter, CorrectionFitter, ErrorBinFitter
from ..redsequence import RedSequenceColorPar
from ..color_background import ColorBackground
from ..galaxy import GalaxyCatalog
from ..catalog import Catalog, Entry
from ..zred_color import ZredColor
from ..utilities import make_nodes, CubicSpline, interpol, RedGalInitialColors
            self.config = conf

        self._galfile = galfile

    def run(self, doRaise=True):
        """
        Run the red-sequence calibration.

        Parameters
        ----------
        doRaise: `bool`, optional
           Raise an error if background cannot be computed for any galaxies
           Default is True. Can be set to False for certain testing.
        """

        gals = GalaxyCatalog.from_galfile(self._galfile)

        if self.config.calib_use_pcol:
            use, = np.where((gals.z > self.config.zrange[0]) &
                            (gals.z < self.config.zrange[1]) &
                            (gals.pcol > self.config.calib_pcut))
        else:
            use, = np.where((gals.z > self.config.zrange[0]) &
                            (gals.z < self.config.zrange[1]) &
                            (gals.p > self.config.calib_pcut))

        if use.size == 0:
            raise RuntimeError("No good galaxies in %s!" % (self._galfile))

        gals = gals[use]

        nmag = self.config.nmag
        ncol = nmag - 1

        # Reference mag nodes for pivot
        pivotnodes = make_nodes(self.config.zrange, self.config.calib_pivotmag_nodesize)

        # Covmat nodes
        covmatnodes = make_nodes(self.config.zrange, self.config.calib_covmat_nodesize)

        # correction nodes
        corrnodes = make_nodes(self.config.zrange, self.config.calib_corr_nodesize)

        # correction slope nodes
        corrslopenodes = make_nodes(self.config.zrange, self.config.calib_corr_slope_nodesize)

        # volume factor (hard coded)
        volnodes = make_nodes(self.config.zrange, 0.01)

        # Start building the par dtype
        dtype = [('pivotmag_z', 'f4', pivotnodes.size),
                 ('pivotmag', 'f4', pivotnodes.size),
                 ('minrefmag', 'f4', pivotnodes.size),
                 ('maxrefmag', 'f4', pivotnodes.size),
                 ('medcol', 'f4', (pivotnodes.size, ncol)),
                 ('medcol_width', 'f4', (pivotnodes.size, ncol)),
                 ('medcol_err_ratio', 'f4', (ncol, )),
                 ('covmat_z', 'f4', covmatnodes.size),
                 ('sigma', 'f4', (ncol, ncol, covmatnodes.size)),
                 ('covmat_amp', 'f4', (ncol, ncol, covmatnodes.size)),
                 ('covmat_slope', 'f4', (ncol, ncol, covmatnodes.size)),
                 ('mag_err_ratio_int', 'f4', (nmag, )),
                 ('mag_err_ratio_slope', 'f4', (nmag, )),
                 ('mag_err_ratio_pivot', 'f4'),
                 ('corr_z', 'f4', corrnodes.size),
                 ('corr', 'f4', corrnodes.size),
                 ('corr_slope_z', 'f4', corrslopenodes.size),
                 ('corr_slope', 'f4', corrslopenodes.size),
                 ('corr_r', 'f4', corrslopenodes.size),
                 ('corr2', 'f4', corrnodes.size),
                 ('corr2_slope', 'f4', corrslopenodes.size),
                 ('corr2_r', 'f4', corrslopenodes.size),
                 ('volume_factor_z', 'f4', volnodes.size),
                 ('volume_factor', 'f4', volnodes.size)]

        # And for each color, make the nodes
        node_dict = {}
        self.ztag = [None] * ncol
        self.ctag = [None] * ncol
        self.zstag = [None] * ncol
        self.stag = [None] * ncol
        for j in range(ncol):
            self.ztag[j] = 'z%02d' % (j)
            self.ctag[j] = 'c%02d' % (j)
            self.zstag[j] = 'zs%02d' % (j)
            self.stag[j] = 'slope%02d' % (j)

            node_dict[self.ztag[j]] = make_nodes(self.config.zrange, self.config.calib_color_nodesizes[j],
                                                 maxnode=self.config.calib_color_maxnodes[j])
            node_dict[self.zstag[j]] = make_nodes(self.config.zrange, self.config.calib_slope_nodesizes[j],
                                                  maxnode=self.config.calib_color_maxnodes[j])

            dtype.extend([(self.ztag[j], 'f4', node_dict[self.ztag[j]].size),
                          (self.ctag[j], 'f4', node_dict[self.ztag[j]].size),
                          (self.zstag[j], 'f4', node_dict[self.zstag[j]].size),
                          (self.stag[j], 'f4', node_dict[self.zstag[j]].size)])

        # Make the pars ... and fill them with the defaults
        self.pars = Entry(np.zeros(1, dtype=dtype))

        self.pars.pivotmag_z = pivotnodes
        self.pars.covmat_z = covmatnodes
        self.pars.corr_z = corrnodes
        self.pars.corr_slope_z = corrslopenodes
        self.pars.volume_factor_z = volnodes
        self.pars.mag_err_ratio_pivot = self.config.calib_err_ratio_pivot

        for j in range(ncol):
            self.pars._ndarray[self.ztag[j]] = node_dict[self.ztag[j]]
            self.pars._ndarray[self.zstag[j]] = node_dict[self.zstag[j]]

        # And a special subset of color galaxies
        if self.config.calib_use_pcol:
            coluse, = np.where(gals.pcol > self.config.calib_color_pcut)
        else:
            coluse, = np.where(gals.p > self.config.calib_color_pcut)

        colgals = gals[coluse]

        # And a placeholder zredstr which allows us to do stuff

        self.zredstr = RedSequenceColorPar(None, config=self.config)