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 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 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) 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] def read_members(catfile):
    """
    Read members associated to a catalog file.

    Parameters
    ----------
    catfile: `str`
       Filename of a cluster catalog.

    Returns
    -------
    members: `redmapper.GalaxyCatalog`
       Member catalog as a GalaxyCatalog
    """
    import os

    from .galaxy import GalaxyCatalog

    memfile = '%s_members.fit' % (catfile.rstrip('.fit'))

    if not os.path.isfile(memfile):
        raise RuntimeError("Could not find member file %s" % (memfile))

    return GalaxyCatalog.from_fits_file(memfile) class SpecPlot(object):
    """
    Class to make plots comparing the cluster spec-z with cluster z_lambda (photo-z).

    In the case of real data, the most likely central is used as a proxy for
    the cluster redshift (so miscenters look like wrong redshifts).  In the
    case of mock data, the mean redshift of the members can be used instead.
    """

    def __init__(self, conf, binsize=0.02, nsig=4.0):
        """
        Instantiate a SpecPlot.

        Parameters
        ----------
        conf: `redmapper.Configuration` or `str`
           Configuration object or filename
        binsize: `float`, optional
           Redshift smoothing bin size.  Default is 0.02.
        nsig: `float`, optional
           Number of sigma to be considered a redshift outlier.  Default is 4.0.
        """
        if not isinstance(conf, Configuration):
            self.config = Configuration(conf)
        else:
            self.config = conf

        self.binsize = binsize
        self.nsig = nsig
        self._withversion = False

    def plot_cluster_catalog(self, cat, title=None, figure_return=False, withversion=False):
        """
        Plot the spectroscopic comparison for a cluster catalog, using the
        default catalog values.

        This plot will compare the catalog cg_spec_z (most likely central
        galaxy redshift) with the cluster z_lambda.

        Parameters
        ----------
        cat: `redmapper.ClusterCatalog`
           Cluster catalog to plot.
        title: `str`, optional
           Title string for plot.  Default is None.
        figure_return: `bool`, optional
           Return the figure instead of saving a png.  Default is False.
        withversion: `bool`, optional
           Plots should be saved with the version string.

        Returns
        -------
        fig: `matplotlib.Figure` or `None`
           Figure to show, if figure_return is True.
        """

        return self.plot_values(cat.cg_spec_z, cat.z_lambda, cat.z_lambda_e, title=title,
                                figure_return=figure_return, withversion=withversion)

    def plot_cluster_catalog_from_members(self, cat, mem, title=None, figure_return=False,
                                          withversion=False):
        """
        Plot the spectroscopic comparison for a cluster catalog, using the
        average member redshift.

        This plot will compare the catalog median redshift with the cluster
        z_lambda.  Generally only useful for simulated catalogs where you have
        complete coverage.

        Parameters
        ----------
        cat: `redmapper.ClusterCatalog`
           Cluster catalog to plot.
        mem: `redmapper.Catalog`
           Member catalog associated with cat.
        title: `str`, optional
           Title string for plot.  Default is None.
        figure_return: `bool`, optional
           Return the figure instead of saving a png.  Default is False.
        withversion: `bool`, optional
           Plots should be saved with the version string.

        Returns
        -------
        fig: `matplotlib.Figure` or `None`
           Figure to show, if figure_return is True.
        """

        # Not sure why this is necessary on the mocks, but there are some -1s...
        ok, = np.where(mem.zspec > 0.0)
        mem = mem[ok]

        a, b = esutil.numpy_util.match(cat.mem_match_id, mem.mem_match_id)

        h, rev = esutil.stat.histogram(a, rev=True)
        mem_zmed = np.zeros(cat.size)
        for i in range(h.size):
            if h[i] == 0:
                continue
            i1a = rev[rev[i]: rev[i + 1]]
            mem_zmed[i] = np.median(mem.zspec[i1a])

        return self.plot_values(mem_zmed, cat.z_lambda, cat.z_lambda_e, title=title,
                                figure_return=figure_return, withversion=withversion)

    def plot_values(self, z_spec, z_phot, z_phot_e, name=r'z_\lambda', specname=r'z_{\mathrm{spec}}',
                    title=None, figure_return=False, calib_zrange=None, withversion=False):
        """
        Make a pretty spectrscopic plot from an arbitrary list of values.

        Parameters
        ----------
        z_spec: `np.array`
           Float array of spectroscopic redshifts
        z_phot: `np.array`
           Float array of photometric redshifts
        z_phot_e: `np.array`
           Float array of photometric redshift errors
        name: `str`, optional
           Name of photo-z field for label.  Default is 'z_lambda'.
        title: `str`, optional
           Title string.  Default is None
        figure_return: `bool`, optional
           Return the figure instead of saving a png.  Default is False.
        calib_zrange: `np.array` or `list`
           2-element array with calibration redshift range to mark.  Default is None.
        withversion: `bool`, optional
           Plots should be saved with the version string.

        Returns
        -------
        fig: `matplotlib.Figure` or `None`
           Figure to show, if figure_return is True.
        """

        # We will need to remove any underscores in name for some usage...
        name_clean = name.replace('_', '')
        name_clean = name_clean.replace('^', '')

        use, = np.where(z_spec > 0.0)

        import matplotlib.pyplot as plt
        from matplotlib.ticker import MultipleLocator

        plot_xrange = np.array([0.0, self.config.zrange[1] + 0.1])

        fig = plt.figure(figsize=(8, 6))
        fig.clf()

        ax = fig.add_subplot(211)

        # Need to do the map, then need to get the levels and make a contour
        z_bins, z_map = self._make_photoz_map(z_spec[use], z_phot[use])

        nlevs = 5
        levbinsize = (1.0 - 0.1) / nlevs
        levs = np.arange(nlevs) * levbinsize + 0.1
        levs = np.append(levs, 10)
        levs = -levs[::-1]

        colors = ['#000000', '#333333', '#666666', '#9A9A9A', '#CECECE', '#FFFFFF']

        ax.contourf(z_bins, z_bins, -z_map.T, levs, colors=colors)
        ax.plot(plot_xrange, plot_xrange, 'b--', linewidth=3)
        ax.set_xlim(plot_xrange)
        ax.set_ylim(plot_xrange)
        ax.tick_params(axis='y', which='major', labelsize=14, length=5, left=True, right=True, direction='in')
        ax.tick_params(axis='y', which='minor', left=True, right=True, direction='in')
        ax.tick_params(axis='x', labelbottom=False, which='major', length=5, direction='in', bottom=True, top=True)
        ax.tick_params(axis='x', which='minor', bottom=True, top=True, direction='in')
        minorLocator = MultipleLocator(0.05)
        ax.yaxis.set_minor_locator(minorLocator)
        minorLocator = MultipleLocator(0.02)
        ax.xaxis.set_minor_locator(minorLocator)
        ax.set_ylabel(r'$%s$' % (specname), fontsize=16)

        if calib_zrange is not None:
            if len(calib_zrange) == 2:
                ylim = ax.get_ylim()
                ax.plot([calib_zrange[0], calib_zrange[0]], ylim, 'k:')
                ax.plot([calib_zrange[1], calib_zrange[1]], ylim, 'k:')

        # Plot the outliers
        bad, = np.where(np.abs(z_phot[use] - z_spec[use]) / z_phot_e[use] > self.nsig)
        if bad.size > 0:
            ax.plot(z_phot[use[bad]], z_spec[use[bad]], 'r*')

        fracout = float(bad.size) / float(use.size)

        fout_label = r'$f_\mathrm{out} = %7.4f$' % (fracout)
        ax.annotate(fout_label, (plot_xrange[0] + 0.1, self.config.zrange[1]),
                    xycoords='data', ha='left', va='top', fontsize=16)

        ax.set_title(title)

        # Compute the bias / scatter
        h, rev = esutil.stat.histogram(z_phot[use], min=0.0, max=self.config.zrange[1]-0.001, rev=True, binsize=self.binsize)
        bins = np.arange(h.size) * self.binsize + 0.0 + self.binsize/2.

        bias = np.zeros(h.size)
        scatter = np.zeros_like(bias)
        errs = np.zeros_like(bias)

        gd, = np.where(h >= 3)
        for ind in gd:
            i1a = rev[rev[ind]: rev[ind + 1]]
            bias[ind] = np.median(z_spec[use[i1a]] - z_phot[use[i1a]])
            scatter[ind] = 1.4862 * np.median(np.abs((z_spec[use[i1a]] - z_phot[use[i1a]]) - bias[ind]))
            errs[ind] = np.median(z_phot_e[use[i1a]])

        # Make the bottom plot
        ax2 = fig.add_subplot(212, sharex=ax)

        ax2.plot(plot_xrange, [0.0, 0.0], 'b--', linewidth=2)
        ax2.plot(bins, bias, 'm-.', label=r'$\mathrm{Bias}$', linewidth=3)
        ax2.plot(bins, scatter/(1. + bins), 'r:', label=r'$\sigma_z / (1 + z)$', linewidth=3)
        ax2.plot(bins, errs/(1. + bins), 'c--', label=r'$\sigma_{%s} / (1 + z)$' % (name_clean), linewidth=3)

        ax2.legend(loc=4, fontsize=14)

        ax2.set_xlim(plot_xrange)
        ax2.set_ylim([-0.03, 0.024])

        if calib_zrange is not None:
            if len(calib_zrange) == 2:
                ylim = ax.get_ylim()
                ax2.plot([calib_zrange[0], calib_zrange[0]], [-0.03, 0.024], 'k:')
                ax2.plot([calib_zrange[1], calib_zrange[1]], [-0.03, 0.024], 'k:')

        ax2.tick_params(axis='both', which='major', labelsize=14, length=5, left=True, right=True, top=True, bottom=True, direction='in')
        ax2.tick_params(axis='y', which='minor', left=True, right=True, direction='in')
        ax2.tick_params(axis='x', which='minor', bottom=True, top=True, direction='in')
        ax2.set_xlabel(r'$%s$' % (name), fontsize=16)
        ax2.set_ylabel(r'$%s - %s$' % (specname, name), fontsize=16)
        minorLocator = MultipleLocator(0.02)
        ax2.xaxis.set_minor_locator(minorLocator)
        minorLocator = MultipleLocator(0.002)
        ax2.yaxis.set_minor_locator(minorLocator)

        fig.subplots_adjust(hspace=0.0)

        fig.tight_layout()

        if not figure_return:
            self._withversion = withversion
            fig.savefig(self.filename)
            plt.close(fig)
        else:
            return fig

    @property
    def filename(self):
        """
        Get the filename that should be used to save the figure.

        Returns
        -------
        filename: `str`
           Formatted filename to save figure.
        """
        return self.config.redmapper_filename('zspec', paths=(self.config.plotpath,),
                                              withversion=self._withversion, filetype='png')

    def _make_photoz_map(self, z_spec, z_photo,
                         dzmax=0.1, nbins_coarse=20, nbins_fine=200, zrange=[0.0, 1.2]):
        """
        Internal method to make a photo-z map for contour plotting.

        Parameters
        ----------
        z_spec: `np.array`
           Float array of z_spec values.
        z_photo: `np.array`
           Float array of photometric redshift values.
        dzmax: `float`, optional
           Maximum delta-z when fitting in bin.  Default is 0.1
        nbins_coarse: `int`, optional
           Number of bins for initial coarse binning.  Default is 20
        nbins_fine: `int`, optional
           Number of bins for smooth fine binning for plot.  Default is 200
        zrange: `list`, optional
           Redshift range to build map.  Default is [0.0, 1.2]

        Returns
        -------
        z_bins: `np.array`
           Float array of redshift bins for map
        z_map_smooth: `np.array`
           Float array of image of smoothed map.
        """

        zmin = z_photo.min()
        zmax = z_photo.max()

        zbinsize = (zmax - zmin) / nbins_coarse
        zbins = np.arange(nbins_coarse) * zbinsize + zmin
        dzbinsize = (2. * dzmax)/(nbins_coarse)
        dzbins = np.arange(nbins_coarse) * dzbinsize - dzmax

        z = zbins + 0.5 * zbinsize
        dz = dzbins + 0.5 * dzbinsize

        spl = np.zeros((nbins_coarse, nbins_coarse))

        for i in range(nbins_coarse):
            use, = np.where((z_photo >= zbins[i]) & (z_photo < (zbins[i] + zbinsize)))

            if use.size == 0:
                continue

            dzvec = z_spec[use] - z_photo[use]

            for j in range(nbins_coarse):
                use2, = np.where((dzvec >= dzbins[j]) & (dzvec < (dzbins[j] + dzbinsize)))
                spl[i, j] = use2.size

            p0 = [use.size, np.median(dzvec), np.std(dzvec)]
            try:
                with np.warnings.catch_warnings():
                    np.warnings.simplefilter('error')
                    coeff, varMatrix = scipy.optimize.curve_fit(gaussFunction, dz, spl[i, :], p0=p0)
            except (RuntimeError, RuntimeWarning, scipy.optimize.OptimizeWarning):
                # set to the default?
                coeff = p0
            spl[i, :] = spl[i, :] / coeff[0]

        maxdz = 0.5 * (z[1] - z[0])

        zbinsize = (zrange[1] - zrange[0]) / nbins_fine
        z_bins = np.arange(nbins_fine) * zbinsize + zrange[0]

        z_map_values = np.zeros((nbins_fine, nbins_fine))

        for i in range(nbins_fine):
            delta = np.abs(z_bins[i] - z)
            index = np.argmin(delta)
            val = delta[index]

            if (val < maxdz):
                dz1 = z_bins - z_bins[i]
                ss = CubicSpline(dz, spl[index, :])
                y = ss(dz1)
                bad, = np.where((np.abs(dz1) > 0.1) | (y < 0.0))
                y[bad] = 0.0
                if y.max() == 0:
                    z_map_values[i, :] = 0.0
                else:
                    z_map_values[i, :] = y / y.max()

        bad = np.where(z_map_values < 1e-3)
        z_map_values[bad] = 0.0

        z_map_smooth = scipy.ndimage.uniform_filter(z_map_values, size=5, mode='nearest')

        return z_bins, z_map_smooth class NzPlot(object):
    """
    Class to make a plot with cluster catalog n(z) with comoving coordinates.
    """

    def __init__(self, conf, binsize=0.02):
        """
        Instantiate a NzPlot.

        Parameters
        ----------
        conf: `redmapper.Configuration` or `str`
           Configuration object or filename
        binsize: `float`, optional
           Redshift smoothing bin size.  Default is 0.02.
        """
        if not isinstance(conf, Configuration):
            self.config = Configuration(conf)
        else:
            self.config = conf

        self.binsize = binsize
        self._redmapper_name = 'nz'
        self._withversion = False

    def plot_cluster_catalog(self, cat, areastr, nosamp=False, withversion=False):
        """
        Plot the n(z) for a cluster catalog, using the default catalog values.

        This plot will sample redshifts from z_lambda +/- z_lambda_e.

        Parameters
        ----------
        cat: `redmapper.ClusterCatalog`
           Cluster catalog to plot.
        areastr: `redmapper.Catalog`
           Area structure, with .z and .area
        nosamp: `bool`, optional
           Do not sample from z_lambda.  Default is False.
        withversion: `bool`, optional
           Plots should be saved with the version string.
        """

        if nosamp:
            zsamp = cat.z_lambda
        else:
            zsamp = np.zeros(cat.size)

            nsampbin = 100
            npzbins = self.config.npzbins

            for i, cluster in enumerate(cat):
                xvals = ((np.arange(nsampbin, dtype=np.float64) / nsampbin) *
                         (cluster.pzbins[-1] - cluster.pzbins[0]) + cluster.pzbins[0])
                yvals = interpol(cluster.pz, cluster.pzbins, xvals)
                pdf = yvals / np.sum(yvals)
                cdf = np.cumsum(pdf, dtype=np.float64)
                cdfi = (cdf * xvals.size).astype(np.int32)
                rand = (np.random.uniform(size=1)*nsampbin).astype(np.int32)
                test, = np.where(cdfi >= rand)
                zsamp[i] = xvals[test[0]]

        self.plot_nz(zsamp, areastr, self.config.zrange,
                     xlabel=r'$z_{\lambda}$',
                     ylabel=r'$n\,(1e4\,\mathrm{clusters} / \mathrm{Mpc}^{3})$',
                     redmapper_name='nz', withversion=withversion)

    def plot_nz(self, z, areastr, zrange, xlabel=None, ylabel=None,
                title=None, redmapper_name='nz', calib_zrange=None,
                withversion=False):
        """
        Plot the n(z) for an arbitrary list of objects

        Parameters
        ----------
        z: `np.array`
           Float array of redshifts
        areastr: `redmapper.Catalog`
           Structure describing area as a function of redshift
        zrange: `np.array` or `list`
           Redshift range to plot
        xlabel: `str`, optional
           Plot x label.  Default is None.
        ylabel: `str`, optional
           Plot y label.  Default is None.
        title: `str`, optional
           Plot title.  Default is None.
        redmapper_name: `str`, optional
           Name to put into filename.  Default is 'nz'
        calib_zrange: `np.array` or `list`
           Calibration redshift range to overplot
        withversion: `bool`, optional
           Plots should be saved with the version string.
        """

        import matplotlib.pyplot as plt

        hist = esutil.stat.histogram(z, min=zrange[0], max=zrange[1]-0.0001, binsize=self.binsize, more=True)
        h = hist['hist']
        zbins = hist['center']

        indices = np.clip(np.searchsorted(areastr.z, zbins), 0, areastr.size - 1)

        vol = np.zeros(zbins.size)
        for i in range(zbins.size):
            vol[i] = (self.config.cosmo.V(np.clip(zbins[i] - self.binsize/2., zrange[0], None),
                                          np.clip(zbins[i] + self.binsize/2., None, zrange[1])) *
                      (areastr.area[indices[i]] / 41252.961))

        dens = h.astype(np.float32) / vol
        err = np.sqrt(dens * vol) / vol

        u, = np.where((dens > 0.0) & (err > 0.0))
        u2, = np.where((dens[u] / err[u]) > 5.0)

        fig = plt.figure(1, figsize=(8, 6))
        fig.clf()

        ax = fig.add_subplot(111)

        ax.errorbar(zbins[u[u2]], dens[u[u2]]*1e4, yerr=err[u[u2]]*1e4, fmt='r.', markersize=8)
        if xlabel is not None:
            ax.set_xlabel(xlabel, fontsize=16)
        if ylabel is not None:
            ax.set_ylabel(ylabel, fontsize=16)
        if title is not None:
            ax.set_title(title, fontsize=16)
        ax.tick_params(axis='both',which='major',labelsize=14)

        if calib_zrange is not None:
            if len(calib_zrange) == 2:
                ylim = ax.get_ylim()
                ax.plot([calib_zrange[0], calib_zrange[0]], ylim, 'k:')
                ax.plot([calib_zrange[1], calib_zrange[1]], ylim, 'k:')

        fig.tight_layout()

        self._redmapper_name = redmapper_name
        self._withversion = withversion
        fig.savefig(self.filename)
        plt.close(fig)

    def plot_redmagic_catalog(self, cat, name, eta, n0, areastr, sample=True,
                              zrange=None, calib_zrange=None, extraname=None,
                              withversion=False):
        """
        Plot the n(z) for a redmagic catalog.

        Parameters
        ----------
        cat: `redmapper.Catalog`
           Galaxy catalog of redmagic galaxies
        name: `str`
           Name (mode) of redmagic catalog
        eta: `float`
           Luminosity cut of catalog (for labeling)
        n0: `float`
           Target n0 of catalog (for labeling)
        areastr: `redmapper.Catalog`
           Structure describing area as a function of redshift
        sample: `bool`, optional
           Sample the p(z)s? Default is True.
        zrange: `np.array` or `list`, optional
           Redshift range to plot.  Default is None (full range)
        calib_zrange: `np.array` or `list`, optional
           Redshift range used for calibration to overplot.  Default is None.
        extraname: `str`, optional
           Extra name to insert (E.g. 'calib').  Default is None
        withversion: `bool`, optional
           Plots should be saved with the version string.
        """

        # Take the first sample
        zsamp = cat.zredmagic_samp[:, 0]

        if zrange is None:
            zrange = self.config.redmagic_zrange

        if extraname is None:
            redmapper_name = 'redmagic_%s_%3.1f-%02d_nz' % (name, eta, int(n0))
        else:
            redmapper_name = 'redmagic_%s_%s_%3.1f-%02d_nz' % (extraname, name, eta, int(n0))

        self.plot_nz(zsamp, areastr, zrange,
                     xlabel=r'$z_{\mathrm{redmagic}}$',
                     ylabel=r'$n\,(1e4\,\mathrm{galaxies} / \mathrm{Mpc}^{3})$',
                     title='%s: %3.1f-%02d' % (name, eta, int(n0)),
                     redmapper_name=redmapper_name,
                     calib_zrange=calib_zrange, withversion=withversion)

    @property
    def filename(self):
        """
        Get the filename that should be used to save the figure.

        Returns
        -------
        filename: `str`
           Formatted filename to save figure.
        """
        return self.config.redmapper_filename(self._redmapper_name, paths=(self.config.plotpath,),
                                              withversion=self._withversion, filetype='png') class VolumeLimitMask(object):
    """
    A class to describe a volume limit mask.

    This is based on combining the red sequence model with depth maps in
    different bands, to determine the highest redshift that a typical
    red-sequence galaxy can be observed with the specified luminosity.

    """

    def __init__(self, config, vlim_lstar, vlimfile=None, use_geometry=False, withversion=True):
        """
        Instantiate a VolumeLimitMask

        If the mask described by maskfile already exists, it will be read in
        directly.  If it does not exist, it will be generated from the depth
        files described in the config parameters, and then stored in vlimfile.

        Parameters
        ----------
        config: `redmapper.Configuration`
           Configuration object
        vlim_lstar: `float`
           Luminosity cutoff (units of L*) for red-sequence volume limit
           computation
        vlimfile: `str`, optional
           Filename to store volume limit mask.  Default is None, which
           means generate the filename of the 'vlim_zmask' type.
        use_geometry: `bool`, optional
           Use the geometric mask info only.  Only use if necessary.
           Default is False.
        withversion: `bool`, optional
           Output filename with redmapper version string.
        """
        self.config = config
        self.vlim_lstar = vlim_lstar

        if vlimfile is None:
            self.vlimfile = self.config.redmapper_filename('vl%02d_vlim_zmask' %
                                                            (int(self.vlim_lstar*10)),
                                                           withversion=withversion)
        else:
            self.vlimfile = vlimfile

        if os.path.isfile(self.vlimfile):
            self._read_mask()
        else:
            if use_geometry:
                self._build_geometry_mask()
            else:
                self._build_mask()
            self._read_mask()

    def _read_mask(self):
        """
        Read an existing volume-limit mask into the VolumeLimitMask structure.
        """

        hdr = fitsio.read_header(self.vlimfile, ext=1)
        if 'PIXTYPE' not in hdr or hdr['PIXTYPE'] != 'HEALSPARSE':
            raise RuntimeError("Need to specify vlimfile in healsparse format.")

        cov_hdr = fitsio.read_header(self.vlimfile, ext='COV')
        nside_coverage = cov_hdr['NSIDE']

        if len(self.config.d.hpix) > 0:
            covpixels = get_healsparse_subpix_indices(self.config.d.nside, self.config.d.hpix,
                                                      self.config.border, nside_coverage)
        else:
            covpixels = None

        self.sparse_vlimmap = healsparse.HealSparseMap.read(self.vlimfile, pixels=covpixels)

        self.nside = self.sparse_vlimmap.nside_sparse
        self.subpix_nside = self.config.d.hpix
        self.subpix_hpix = self.config.d.nside
        self.subpix_border = self.config.border

    def _build_mask(self):
        """
        Build a VolumeLimitMask from the parameters in the config file, and
        store the mask in self.vlimfile
        """

        # Make some checks to make sure we can build a volume limit mask
        if self.config.depthfile is None or not os.path.isfile(self.config.depthfile):
            raise RuntimeError("Cannot create a volume limit mask without a depth file")
        for fname in self.config.vlim_depthfiles:
            if not os.path.isfile(fname):
                raise RuntimeError("Could not find specified vlim_depthfile %s" % (fname))

        # Read in the red-sequence parameters
        zredstr = RedSequenceColorPar(self.config.parfile, fine=True)

        # create the redshift bins
        zbinsize = 0.001 # arbitrary fine bin
        nzbins = int(np.ceil((self.config.zrange[1] - self.config.zrange[0]) / zbinsize))
        # Note that we want to start one step above the low redshift range
        zbins = np.arange(nzbins) * zbinsize + self.config.zrange[0] + zbinsize

        # magnitude limits
        limmags = zredstr.mstar(zbins) - 2.5 * np.log10(self.vlim_lstar)

        # get the reference index
        ref_ind = self.config.bands.index(self.config.refmag)

        # Read in the primary depth structure
        sparse_depthmap = healsparse.HealSparseMap.read(self.config.depthfile)

        dtype_vlimmap = [('fracgood', 'f4'),
                         ('zmax', 'f4')]

        sparse_vlimmap = healsparse.HealSparseMap.make_empty(sparse_depthmap.nside_coverage,
                                                             sparse_depthmap.nside_sparse,
                                                             dtype=dtype_vlimmap,
                                                             primary='fracgood')

        validPixels = sparse_depthmap.valid_pixels
        depthValues = sparse_depthmap.get_values_pix(validPixels)
        vlimmap = np.zeros(validPixels.size, dtype=dtype_vlimmap)
        vlimmap['fracgood'] = depthValues['fracgood']

        lo, = np.where(depthValues['m50'] <= limmags.min())
        vlimmap['zmax'][lo] = zbins.min()
        hi, = np.where(depthValues['m50'] >= limmags.max())
        vlimmap['zmax'][hi] = zbins.max()
        mid, = np.where((depthValues['m50'] > limmags.min()) & (depthValues['m50'] < limmags.max()))
        if mid.size > 0:
            l = np.searchsorted(limmags, depthValues['m50'][mid], side='right')
            vlimmap['zmax'][mid] = zbins[l]

        # Read in any additional depth maps
        for i, depthfile in enumerate(self.config.vlim_depthfiles):
            sparse_depthmap2, hdr2 = healsparse.HealSparseMap.read(depthfile, header=True)

            validPixels2 = sparse_depthmap2.valid_pixels
            depthValues2 = sparse_depthmap2.get_values_pix(validPixels2)

            nsig = hdr2['NSIG']
            zp = hdr2['ZP']

            # find mag name thing...
            # Note this is validated in the config read
            map_ind = self.config.bands.index(self.config.vlim_bands[i])

            # match pixels
            a, b = esutil.numpy_util.match(validPixels, validPixels2)

            n2 = self.config.vlim_nsigs[i]**2.
            flim_in = 10.**((depthValues2['limmag'][b] - zp) / (-2.5))
            fn = np.clip((flim_in**2. * depthValues2['exptime'][b]) / (nsig**2.) - flim_in, 0.001, None)
            flim_mask = (n2 + np.sqrt(n2**2. + 4.*depthValues2['exptime'][b] * n2 * fn)) / (2.*depthValues2['exptime'][b])
            lim_mask = np.zeros(vlimmap.size)
            lim_mask[a] = zp - 2.5*np.log10(flim_mask)

            zinds = np.searchsorted(zredstr.z, zbins, side='right')

            limmags_temp = zredstr.mstar(zbins) - 2.5*np.log10(self.vlim_lstar)
            refmag_lim = limmags_temp.copy()

            if (map_ind == ref_ind):
                self.config.logger.info('Warning: vlim_band %s is the same as the reference band!  Skipping...' % (self.config.vlim_bands[i]))
            else:
                if map_ind < ref_ind:
                    # Need to go blueward
                    for jj in range(ref_ind - 1, map_ind - 1, -1):
                        limmags_temp += (zredstr.c[zinds, jj] + zredstr.slope[zinds, jj] * (refmag_lim - zredstr.pivotmag[zinds]))
                else:
                    # Need to go redward
                    for jj in range(ref_ind, map_ind):
                        limmags_temp -= (zredstr.c[zinds, jj] + zredstr.slope[zinds, jj] * (refmag_lim - zredstr.pivotmag[zinds]))

            # adjust zmax with zmax_temp
            zmax_temp = np.zeros(vlimmap.size)

            lo, = np.where(lim_mask <= limmags_temp.min())
            zmax_temp[lo] = zbins.min()
            hi, = np.where(lim_mask >= limmags_temp.max())
            zmax_temp[hi] = zbins.max()
            mid, = np.where((lim_mask > limmags_temp.min()) & (lim_mask < limmags_temp.max()))
            if mid.size > 0:
                l = np.clip(np.searchsorted(limmags_temp, lim_mask[mid], side='right'), 0, zbins.size - 1)
                zmax_temp[mid] = zbins[l]

            limited, = np.where(zmax_temp < vlimmap['zmax'])
            vlimmap['zmax'][limited] = zmax_temp[limited]


        gd, = np.where(vlimmap['zmax'] > zbins[0])

        sparse_vlimmap.update_values_pix(validPixels[gd], vlimmap[gd])

        sparse_vlimmap.write(self.vlimfile)

    def _build_geometry_mask(self):
        """
        Build a VolumeLimitMask from the geometric mask from the config file,
        and store the mask in self.vlimfile
        """

        if self.config.maskfile is None or not os.path.isfile(self.config.maskfile):
            raise RuntimeError("Cannot create a geometry volume limit mask without a mask file")

        sparse_mask = healsparse.HealSparseMap.read(self.config.maskfile)

        dtype_vlimmap = [('fracgood', 'f4'),
                         ('zmax', 'f4')]

        sparse_vlimmap = healsparse.HealSparseMap.make_empty(sparse_mask.nside_coverage,
                                                             sparse_mask.nside_sparse,
                                                             dtype=dtype_vlimmap,
                                                             primary='fracgood')

        validPixels = sparse_mask.valid_pixels
        maskValues = sparse_mask.get_values_pix(validPixels)
        vlimmap = np.zeros(validPixels.size, dtype=dtype_vlimmap)
        vlimmap['fracgood'] = maskValues
        vlimmap['zmax'] = self.config.zrange[1]

        sparse_vlimmap.update_values_pix(validPixels, vlimmap)

        sparse_vlimmap.write(self.vlimfile)

    def calc_zmax(self, ras, decs, get_fracgood=False):
        """
        Calculate the maximum redshifts associated with a set of ra/decs.

        Parameters
        ----------
        ras: `np.array` or `float`
           Float array of right ascensions
        decs: `np.array` or `float`
           Float array of declinations
        get_fracgood: `bool`, optional
           Also retrieve the fracgood pixel coverage.  Default is False.

        Returns
        -------
        zmax: `np.array`
           Float array of maximum redshifts
        fracgood: `np.array`
           Float array of fracgood, if get_fracgood=True
        """

        if (len(ras) != len(decs)):
            raise ValueError("ras, decs must be same length")

        values = self.sparse_vlimmap.get_values_pos(ras, decs, lonlat=True)

        bad, = np.where(np.abs(decs) > 90.0)
        values['zmax'][bad] = hp.UNSEEN
        values['fracgood'][bad] = 0.0

        if not get_fracgood:
            return np.clip(values['zmax'], 0.0, None)
        else:
            return (np.clip(values['zmax'], 0.0, None), values['fracgood'])

    def get_areas(self):
        """
        Retrieve the area structure (area as a function of redshift) associated
        with the volume-limit mask.

        Returns
        -------
        astr: `redmapper.Catalog`
           Area structure catalog, with .z and .area
        """

        zbinsize = self.config.area_finebin
        zbins = np.arange(self.config.zrange[0], self.config.zrange[1], zbinsize)

        astr = Catalog(np.zeros(zbins.size, dtype=[('z', 'f4'),
                                                   ('area', 'f4')]))
        astr.z = zbins

        pixsize = hp.nside2pixarea(self.nside, degrees=True)

        validPixels = self.sparse_vlimmap.valid_pixels
        zmax = self.sparse_vlimmap.get_values_pix(validPixels)['zmax']
        st = np.argsort(zmax)

        fracgoods = self.sparse_vlimmap.get_values_pix(validPixels)['fracgood'][st]

        inds = np.searchsorted(zmax[st], zbins, side='right')

        lo = (inds <= 0)
        astr.area[lo] = np.sum(fracgoods.astype(np.float64)) * pixsize

        if np.sum(~lo) > 0:
            carea = pixsize * np.cumsum(fracgoods, dtype=np.float64)
            astr.area[~lo] = carea[carea.size - inds[~lo]]

        return astr class VolumeLimitMaskFixed(object):
    """
    A class to describe a volume limit mask with a fixed redshift maximum.

    This class is used as a placeholder when there is no depth information to
    construct a true VolumeLimitMask.
    """

    def __init__(self, config):
        """
        Instantiate a VolumeLimitMaskFixed

        The maximum redshift is set by config.zrange[1] (the max redshift in
        the config file).

        Parameters
        ----------
        config: `redmapper.Configuration`
           Configuration object
        """
        self.z_max = config.zrange[1]
        self.zrange = config.zrange
        self.zbinsize = config.area_finebin
        self.area = config.area

    def calc_zmax(self, ra, dec):
        """
        Calculate the maximum redshifts associated with a set of ra/decs.

        Parameters
        ----------
        ra: `np.array` or `float`
           Float array of right ascensions
        dec: `np.array` or `float`
           Float array of right ascensions

        Returns
        -------
        zmax: `float`
           Maximum redshift from the config file.
        """

        return self.z_max

    def get_areas(self):
        """
        Retrieve the area structure (area as a function of redshift) associated
        with the fixed-redshift mask.

        Returns
        -------
        astr: `redmapper.Catalog`
           Area structure catalog, with .z and .area
        """

        zbins = np.arange(self.zrange[0], self.zrange[1], self.zbinsize)

        astr = Catalog(np.zeros(zbins.size, dtype=[('z', 'f4'),
                                                   ('area', 'f4')]))
        astr.z = zbins

        astr.area = self.area

        return astr class RedmagicSelector(object):
    """
    Class to select redMaGiC galaxies.
    """
    def __init__(self, conf, vlim_masks=None):
        """
        Instantiate a RedmagicSelector

        Parameters
        ----------
        conf: `redmapper.Configuration` or `str
           Configuration object or config filename
        vlim_masks: `OrderedDict`, optional
           Dictionary of vlim_masks.  Will read in if not set.
        """

        if not isinstance(conf, Configuration):
            self.config = Configuration(conf)
        else:
            self.config = conf

        redmagicfilepath = os.path.dirname(self.config.redmagicfile)

        self.calib_data = OrderedDict()
        with fitsio.FITS(self.config.redmagicfile) as fits:
            # Number of modes is number of binary extentions
            self.n_modes = len(fits) - 1

            for ext in range(self.n_modes):
                data = Entry(fits[ext + 1].read())

                try:
                    name = data.name.decode().rstrip()
                except AttributeError:
                    name = data.name.rstrip()

                self.calib_data[name] = data

        self.modes = self.calib_data.keys()

        self.zredstr = RedSequenceColorPar(self.config.parfile, fine=True)

        if vlim_masks is None:
            self.vlim_masks = OrderedDict()

            for mode in self.modes:
                try:
                    vmaskfile = self.calib_data[mode].vmaskfile.decode().rstrip()
                except AttributeError:
                    vmaskfile = self.calib_data[mode].vmaskfile.rstrip()

                if vmaskfile == '':
                    # There is no vmaskfile, we need to do a fixed area one
                    self.vlim_masks[mode] = VolumeLimitMaskFixed(self.config)
                else:
                    vmaskfile = os.path.join(redmagicfilepath,
                                             os.path.basename(vmaskfile))
                    if not os.path.isfile(vmaskfile):
                        raise RuntimeError("Could not find vmaskfile %s.  Must be in same path as redmagic calibration file %s." % (vmaskfile, os.path.abspath(self.config.redmagicfile)))

                    #if os.path.isfile(vmaskfile):
                    #    vmaskfile = vmaskfile
                    #elif os.path.isfile(os.path.join(self.config.configpath, vmaskfile)):
                    #    vmaskfile = os.path.join(self.config.configpath, vmaskfile)
                    #else:
                    #    raise RuntimeError("Could not find vmaskfile %s" % (vmaskfile))
                    self.vlim_masks[mode] = VolumeLimitMask(self.config,
                                                            self.calib_data[mode].etamin,
                                                            vlimfile=vmaskfile)
        else:
            # Check that it's an OrderedDict?  Must it be ordered?
            self.vlim_masks = vlim_masks

        self.spec = None

    def select_redmagic_galaxies(self, gals, mode, return_indices=False):
        """
        Select redMaGiC galaxies from a galaxy catalog, according to the mode.

        Parameters
        ----------
        gals: `redmapper.GalaxyCatalog`
           Catalog of galaxies for redMaPPer
        mode: `str`
           redMaGiC mode to select
        return_indices: `bool`, optional
           Return the indices of the galaxies selected.  Default is False.

        Returns
        -------
        redmagic_catalog: `redmapper.GalaxyCatalog`
           Catalog of redMaGiC galaxies
        indices: `np.ndarray`
           Integer array of selection (if return_indices is True)
        """

        # Check if we have to decode mode (py2/py3)
        if hasattr(mode, 'decode'):
            _mode = mode.decode()
        else:
            _mode = mode

        if _mode not in self.modes:
            raise RuntimeError("Requested redMaGiC mode %s not available." % (_mode))

        calstr = self.calib_data[_mode]

        # Takes in galaxies...
        # Which are the possibly red galaxies?
        lstar_cushion = calstr.lstar_cushion
        z_cushion = calstr.z_cushion
        z_buffer = calstr.buffer

        mstar_init = self.zredstr.mstar(gals.zred_uncorr)

        cut_zrange = [calstr.cost_zrange[0] - z_cushion - z_buffer,
                      calstr.cost_zrange[1] + z_cushion + z_buffer]
        minlstar = np.clip(np.min(calstr.etamin) - lstar_cushion, 0.1, None)

        red_poss_mask = ((gals.zred_uncorr > cut_zrange[0]) &
                         (gals.zred_uncorr < cut_zrange[1]) &
                         (gals.chisq < calstr.maxchi) &
                         (gals.refmag < (mstar_init - 2.5*np.log10(minlstar))))

        # Creates a new catalog...
        zredmagic = np.copy(gals.zred_uncorr)
        zredmagic_e = np.copy(gals.zred_uncorr_e)
        try:
            zredmagic_samp = np.copy(gals.zred_samp)
        except (ValueError, AttributeError) as e:
            # Sample from zred + zred_e (not optimal, for old catalogs)
            zredmagic_samp = np.zeros((zredmagic.size, 1))
            zredmagic_samp[:, 0] = np.random.normal(loc=zredmagic,
                                                    scale=zredmagic_e,
                                                    size=zredmagic.size)

        spl = CubicSpline(calstr.nodes, calstr.cmax, fixextrap=True)
        chi2max = np.clip(spl(gals.zred_uncorr), 0.1, calstr.maxchi)

        if calstr.run_afterburner:
            spl = CubicSpline(calstr.corrnodes, calstr.bias, fixextrap=True)
            offset = spl(gals.zred_uncorr)
            zredmagic -= offset

            if calstr.apply_afterburner:
                for i in range(zredmagic_samp.shape[1]):
                    zredmagic_samp[:, i] -= offset

            spl = CubicSpline(calstr.corrnodes, calstr.eratio, fixextrap=True)
            zredmagic_e *= spl(gals.zred_uncorr)

        # Compute mstar
        mstar = self.zredstr.mstar(zredmagic)

        # Compute the maximum redshift
        vmask = self.vlim_masks[_mode]
        zmax = vmask.calc_zmax(gals.ra, gals.dec)

        # Do the redmagic selection
        gd, = np.where((gals.chisq < chi2max) &
                       (gals.refmag < (mstar - 2.5 * np.log10(calstr.etamin))) &
                       (zredmagic < zmax) &
                       (red_poss_mask))

        redmagic_catalog = GalaxyCatalog(np.zeros(gd.size, dtype=[('id', 'i8'),
                                                                  ('ra', 'f8'),
                                                                  ('dec', 'f8'),
                                                                  ('refmag', 'f4'),
                                                                  ('refmag_err', 'f4'),
                                                                  ('mag', 'f4', self.config.nmag),
                                                                  ('mag_err', 'f4', self.config.nmag),
                                                                  ('lum', 'f4'),
                                                                  ('zredmagic', 'f4'),
                                                                  ('zredmagic_e', 'f4'),
                                                                  ('zredmagic_samp', 'f4', self.config.zred_nsamp),
                                                                  ('chisq', 'f4'),
                                                                  ('zspec', 'f4')]))

        if gd.size == 0:
            if return_indices:
                return redmagic_catalog, gd
            else:
                return redmagic_catalog

        redmagic_catalog.id = gals.id[gd]
        redmagic_catalog.ra = gals.ra[gd]
        redmagic_catalog.dec = gals.dec[gd]
        redmagic_catalog.refmag = gals.refmag[gd]
        redmagic_catalog.refmag_err = gals.refmag_err[gd]
        redmagic_catalog.mag[:, :] = gals.mag[gd, :]
        redmagic_catalog.mag_err[:, :] = gals.mag_err[gd, :]
        redmagic_catalog.zredmagic = zredmagic[gd]
        redmagic_catalog.zredmagic_e = zredmagic_e[gd]
        redmagic_catalog.zredmagic_samp = zredmagic_samp[gd, :]
        redmagic_catalog.chisq = gals.chisq[gd]

        # Compute the luminosity
        redmagic_catalog.lum = 10.**((mstar[gd] - redmagic_catalog.refmag) / 2.5)

        # In the future, add absolute magnitude calculations, but that will
        # require some k-corrections.

        # Compute the zspec (check this)
        if 'ztrue' in gals.dtype.names:
            # We have truth zspec
            redmagic_catalog.zspec = gals.ztrue[gd]
        elif 'zspec' in gals.dtype.names:
            # We have already done a zspec match
            redmagic_catalog.zspec = gals.zspec[gd]
        else:
            # We need to do a zspec match here
            if self.spec is None:
                self.config.logger.info("Reading in spectroscopic information...")

                self.spec = GalaxyCatalog.from_fits_file(self.config.specfile)
                use, = np.where(self.spec.z_err < 0.001)
                self.spec = self.spec[use]
                self.config.logger.info("Done reading in spectroscopic information.")

            redmagic_catalog.zspec[:] = -1.0
            i0, i1, dists = self.spec.match_many(redmagic_catalog.ra, redmagic_catalog.dec, 3./3600., maxmatch=1)
            redmagic_catalog.zspec[i0] = self.spec.z[i1]

        if return_indices:
            return redmagic_catalog, gd
        else:
            return redmagic_catalog class RunRedmagicTask(object):
    """
    Class to run redmagic on a full catalog.
    """

    def __init__(self, configfile, path=None):
        """
        Instantiate a RunRedmagicTask

        Parameters
        ----------
        configfile: `str`
           Configuration yaml filename
        path: `str`, optional
l           Output path.  Default is None, use same absolute
           path as configfile.
        """
        if path is None:
            outpath = os.path.dirname(os.path.abspath(configfile))
        else:
            outpath = path

        self.config = Configuration(configfile, outpath=path)

    def run(self, modes=None, clobber=False, do_plots=True, n_randoms=None, rng=None):
        """
        Run redMaGiC selection over a full catalog.

        The modes are optional, if not specified all the modes
        will be run.

        Parameters
        ----------
        modes: `list`, optional
           List of strings of modes to run
        clobber: `bool`, optional
           Overwrite any existing files.  Default is False.
        do_plots: `bool`, optional
           Make the output plots
        n_randoms: `int`, optional
           If None, then 10x the number of redmagic galaxies are generated.
           If 0, then no randoms are generated.
           If >0, then that many randoms are generated.
        rng : `np.random.RandomState`, optional
           Pre-set random number generator.  Default is None.
        """
        self.config.start_file_logging()

        if rng is None:
            rng = np.random.RandomState()

        if not self.config.galfile_pixelized:
            raise ValueError("Code only runs with pixelized galfile.")

        # Prepare the redMaGiC selector
        selector = RedmagicSelector(self.config)

        if modes is None:
            modes = selector.modes

        n_modes = len(modes)

        # Check if files exist, clobber is False
        filenames = [''] * n_modes
        for i, mode in enumerate(modes):
            filenames[i] = self.config.redmapper_filename('redmagic_%s' % (mode), withversion=True)

            if os.path.isfile(filenames[i]) and not clobber:
                raise RuntimeError("redMaGiC file %s already exists, and clobber is False" % (filenames[i]))

        # Loop over all pixels in the galaxy table
        tab = Entry.from_fits_file(self.config.galfile)

        started = [False] * n_modes

        self.config.logger.info("Making redMaGiC selection for %d modes and %d pixels" % (n_modes, tab.hpix.size))
        if self.config.has_truth:
            self.config.logger.info("Using truth information for zspec")

        for i, pix in enumerate(tab.hpix):
            gals = GalaxyCatalog.from_galfile(self.config.galfile,
                                              zredfile=self.config.zredfile,
                                              nside=tab.nside,
                                              hpix=pix,
                                              border=0.0,
                                              truth=self.config.has_truth)

            # Loop over all modes
            for j, mode in enumerate(modes):
                # Select the red galaxies
                red_gals = selector.select_redmagic_galaxies(gals, mode)

                # Spool out redMaGiC galaxies
                if not started[j]:
                    # write a new file (and overwrite if necessary, since we
                    # already did the clobber check)
                    red_gals.to_fits_file(filenames[j], clobber=True)
                    started[j] = True
                else:
                    with fitsio.FITS(filenames[j], mode='rw') as fits:
                        fits[1].append(red_gals._ndarray)

        # Load in catalogs and make plots!
        if do_plots:
            import matplotlib.pyplot as plt

            for j, mode in enumerate(modes):
                gals = GalaxyCatalog.from_fits_file(filenames[j])

                nzplot = NzPlot(self.config, binsize=self.config.redmagic_calib_zbinsize)
                nzplot.plot_redmagic_catalog(gals, mode, selector.calib_data[mode].etamin,
                                             selector.calib_data[mode].n0,
                                             selector.vlim_masks[mode].get_areas(),
                                             zrange=selector.calib_data[mode].zrange,
                                             withversion=True)

                okspec, = np.where(gals.zspec > 0.0)
                if okspec.size > 0:
                    specplot = SpecPlot(self.config)
                    fig = specplot.plot_values(gals.zspec[okspec], gals.zredmagic[okspec],
                                               gals.zredmagic_e[okspec],
                                               name=r'z_{\mathrm{redmagic}}',
                                               title='%s: %3.1f-%02d' %
                                               (mode, selector.calib_data[mode].etamin,
                                                int(selector.calib_data[mode].n0)),
                                               figure_return=True)
                    fig.savefig(self.config.redmapper_filename('redmagic_zspec_%s_%3.1f-%02d' %
                                                               (mode, selector.calib_data[mode].etamin,
                                                                int(selector.calib_data[mode].n0)),
                                                               paths=(self.config.plotpath,),
                                                               withversion=True,
                                                               filetype='png'))
                    plt.close(fig)

        # Need a check on when to kick out
        if (n_randoms is not None and n_randoms == 0):
            # We are not generating randoms
            return

        # Random generation
        for j, mode in enumerate(modes):
            if isinstance(selector.vlim_masks[mode], VolumeLimitMaskFixed):
                self.config.logger.info("Cannot construct randoms for %s, because we don't have geometry/depth info." % (mode))
                continue

            gals = GalaxyCatalog.from_fits_file(filenames[j])

            rand_generator = RedmagicGenerateRandoms(self.config, selector.vlim_masks[mode], gals)
            if n_randoms is None:
                _n_randoms = gals.size * 10
            else:
                _n_randoms = n_randoms

            randfile = self.config.redmapper_filename('redmagic_%s_randoms' % (mode), withversion=True)

            rand_generator.generate_randoms(_n_randoms, randfile, clobber=clobber, rng=rng)

        self.config.stop_file_logging()
from collections import OrderedDict
from ..configuration import Configuration
from ..fitters import MedZFitter
from ..redsequence import RedSequenceColorPar
from ..galaxy import GalaxyCatalog
from ..catalog import Catalog, Entry
from ..utilities import make_nodes, CubicSpline, interpol, read_members
from ..plotting import SpecPlot, NzPlot
from ..volumelimit import VolumeLimitMask, VolumeLimitMaskFixed
from .redmagic_selector import RedmagicSelector
from .redmagictask import RunRedmagicTask
import os
import numpy as np
import fitsio
import time
import scipy.optimize
import esutil
import healpy as hp
import copy
        import matplotlib.pyplot as plt
            r = np.random.random(gd.size)
            st = np.argsort(r)
            afterburner_use = gd[st[0: ntrain]]

            # Compute the volume
            vmask = vlim_masks[self.config.redmagic_names[i]]
            astr = vlim_areas[self.config.redmagic_names[i]]

            aind = np.searchsorted(astr.z, zbins)
            z_areas = astr.area[aind]

            volume = np.zeros(zbins.size)
            for j in range(zbins.size):
                volume[j] = (self.config.cosmo.V(zbins[j] - self.config.redmagic_calib_zbinsize/2.,
                                                 zbins[j] + self.config.redmagic_calib_zbinsize/2.) *
                                  (z_areas[j] / 41252.961))

            zmax = vmask.calc_zmax(gals.ra[red_poss], gals.dec[red_poss])

            # Compute the density based on the histogram above and the volume
            dens = h.astype(np.float64) / volume

            if not self.config.redmagic_use_constchi:
                bad, = np.where(dens < self.config.redmagic_n0s[i] * 1e-4)
                if bad.size > 0:
                    self.config.logger.info("Warning: not enough galaxies at z=%s" % (zbins[bad].__str__()))

            # get starting values
            cmaxvals = np.zeros(nodes.size)

            if self.config.redmagic_use_constchi:
                cmaxvals[:] = self.config.redmagic_constchis[i]
                _constchi = self.config.redmagic_constchis[i]
            else:
                _constchi = None
                aind = np.searchsorted(astr.z, nodes)
                test_areas = astr.area[aind]

                test_vol = np.zeros(nodes.size)
                for j in range(nodes.size):
                    test_vol[j] = (self.config.cosmo.V(nodes[j] - self.config.redmagic_calib_zbinsize/2.,
                                                       nodes[j] + self.config.redmagic_calib_zbinsize/2.) *
                                   (test_areas[j] / 41252.961))

                for j in range(nodes.size):
                    zrange = [nodes[j] - self.config.redmagic_calib_zbinsize/2.,
                              nodes[j] + self.config.redmagic_calib_zbinsize/2.]
                    if j == 0:
                        zrange = [nodes[j],
                                  nodes[j] + self.config.redmagic_calib_zbinsize]
                    elif j == nodes.size - 1:
                        zrange = [nodes[j] - self.config.redmagic_calib_zbinsize,
                                  nodes[j]]

                    u, = np.where((gals.zuse[red_poss] > zrange[0]) &
                                  (gals.zuse[red_poss] < zrange[1]))

                    st = np.argsort(gals.chisq[red_poss[u]])
                    test_den = np.arange(u.size) / test_vol[j]
                    ind = np.clip(np.searchsorted(test_den, self.config.redmagic_n0s[i] * 1e-4), 0, test_den.size - 1)
                    cmaxvals[j] = gals.chisq[red_poss[u[st[ind]]]]

            # minimize chisquared parameters
            rmfitter = RedmagicParameterFitter(nodes, corrnodes,
                                               gals.zuse[red_poss], gals.zuse_e[red_poss],
                                               gals.chisq[red_poss], mstar_init[red_poss],
                                               gals.zcal[red_poss], gals.zcal_e[red_poss],
                                               gals.refmag[red_poss], zsamp[red_poss],
                                               zmax,
                                               self.config.redmagic_etas[i],
                                               self.config.redmagic_n0s[i],
                                               volume, cost_zrange,
                                               self.config.redmagic_calib_zbinsize,
                                               zredstr,
                                               maxchi=self.config.redmagic_calib_chisqcut,
                                               ab_use=afterburner_use,
                                               ab_apply=self.config.redmagic_apply_afterburner_zsamp,
                                               constchi=_constchi)

            if not self.config.redmagic_use_constchi:
                self.config.logger.info("Fitting first pass...")
                cmaxvals = rmfitter.fit(cmaxvals, afterburner=False)

            # default is no bias or eratio
            biasvals = np.zeros(corrnodes.size)
            eratiovals = np.ones(corrnodes.size)

            # run with afterburner
            if self.config.redmagic_run_afterburner:
                self.config.logger.info("Fitting with afterburner...")

                if not self.config.redmagic_use_constchi:
                    # Let's look at 5 iterations here...
                    for k in range(5):
                        self.config.logger.info("Afterburner iteration %d" % (k))
                        # Fit the bias and eratio...
                        biasvals, eratiovals = rmfitter.fit_bias_eratio(cmaxvals, biasvals, eratiovals)
                        cmaxvals = rmfitter.fit(cmaxvals, biaspars=biasvals, eratiopars=eratiovals, afterburner=True)

                # And a last fit of bias/eratio
                biasvals, eratiovals = rmfitter.fit_bias_eratio(cmaxvals, biasvals, eratiovals)

            rmfitter = None

            # Record the calibrations
            calstr.cmax[:] = cmaxvals
            calstr.bias[:] = biasvals
            calstr.eratio[:] = eratiovals

            calstr.to_fits_file(self.config.redmagicfile, clobber=False, extname=self.config.redmagic_names[i])

            # Do the redmagic selection on our calibration galaxies
            selector = RedmagicSelector(self.config, vlim_masks=vlim_masks)
            redgals, gd = selector.select_redmagic_galaxies(gals, self.config.redmagic_names[i],
                                                            return_indices=True)
            gals.zredmagic[gd] = redgals.zredmagic
            gals.zredmagic_e[gd] = redgals.zredmagic_e
            gals.zredmagic_samp[gd, :] = redgals.zredmagic_samp

            # make pretty plots

            nzplot = NzPlot(self.config, binsize=self.config.redmagic_calib_zbinsize)