# Path: redmapper/utilities.py
# def gaussFunction(x, *p):
#     """
#     Compute a normalizes Gaussian G(x) for a given set of parameters
# 
#     Parameters
#     ----------
#     x: `np.array`
#        Float array of x values
#     A: `float`
#        Normalization of the Gaussian
#     mu: `float`
#        Mean value (mu) for the Gaussian
#     sigma: `float`
#        Gaussian sigma
# 
#     Returns
#     -------
#     pdf: `np.array`
#        Float array of Gaussian pdf values
#     """
#     A, mu, sigma = p
#     return A*np.exp(-(x-mu)**2./(2.*sigma**2))
# 
# Path: redmapper/utilities.py
# 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]


import fitsio
import esutil
import numpy as np
from .utilities import gaussFunction
from .utilities import interpol
        pdis = np.sqrt(pdis**2. + self.config.wcen_rsoft**2.)

        lum = 10.**((self.cluster.mstar - self.cluster.neighbors.refmag) / (2.5))

        # Put a floor on w when we have a strange candidate at the edge that doesn't
        # match any good galaxies
        w = np.zeros(use.size) + 1e-3
        for i in range(use.size):
            # need to filter on r_lambda...
            subgal, = np.where(i1 == i)
            if subgal.size > 0:
                inside, = np.where(pdis[subgal] < self.cluster.r_lambda)
                if inside.size > 0:
                    indices = u[i2[subgal[inside]]]
                    if self.config.wcen_uselum:
                        w[i] = np.log(np.sum(self.cluster.neighbors.p[indices] * lum[indices] /
                                             pdis[subgal[inside]]) /
                                      ((1. / self.cluster.r_lambda) *
                                       np.sum(self.cluster.neighbors.p[indices] * lum[indices])))
                    else:
                        w[i] = np.log(np.sum(self.cluster.neighbors.p[indices] /
                                             pdis[subgal[inside]]) /
                                      ((1. / self.cluster.r_lambda) *
                                       np.sum(self.cluster.neighbors.p[indices])))

        sigscale = np.sqrt((np.clip(self.cluster.Lambda, None, self.config.wcen_maxlambda) / self.cluster.scaleval) / self.config.wcen_pivot)

        # scale with richness for Poisson errors
        sig = self.config.lnw_cen_sigma / sigscale


# Path: redmapper/utilities.py
# def gaussFunction(x, *p):
#     """
#     Compute a normalizes Gaussian G(x) for a given set of parameters
# 
#     Parameters
#     ----------
#     x: `np.array`
#        Float array of x values
#     A: `float`
#        Normalization of the Gaussian
#     mu: `float`
#        Mean value (mu) for the Gaussian
#     sigma: `float`
#        Gaussian sigma
# 
#     Returns
#     -------
#     pdf: `np.array`
#        Float array of Gaussian pdf values
#     """
#     A, mu, sigma = p
#     return A*np.exp(-(x-mu)**2./(2.*sigma**2))
# 
# Path: redmapper/utilities.py
# 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]


# Path: redmapper/centering.py
import fitsio
import esutil
import numpy as np
from .utilities import gaussFunction
from .utilities import interpol

        pdis = np.sqrt(pdis**2. + self.config.wcen_rsoft**2.)

        lum = 10.**((self.cluster.mstar - self.cluster.neighbors.refmag) / (2.5))

        # Put a floor on w when we have a strange candidate at the edge that doesn't
        # match any good galaxies
        w = np.zeros(use.size) + 1e-3
        for i in range(use.size):
            # need to filter on r_lambda...
            subgal, = np.where(i1 == i)
            if subgal.size > 0:
                inside, = np.where(pdis[subgal] < self.cluster.r_lambda)
                if inside.size > 0:
                    indices = u[i2[subgal[inside]]]
                    if self.config.wcen_uselum:
                        w[i] = np.log(np.sum(self.cluster.neighbors.p[indices] * lum[indices] /
                                             pdis[subgal[inside]]) /
                                      ((1. / self.cluster.r_lambda) *
                                       np.sum(self.cluster.neighbors.p[indices] * lum[indices])))
                    else:
                        w[i] = np.log(np.sum(self.cluster.neighbors.p[indices] /
                                             pdis[subgal[inside]]) /
                                      ((1. / self.cluster.r_lambda) *
                                       np.sum(self.cluster.neighbors.p[indices])))

        sigscale = np.sqrt((np.clip(self.cluster.Lambda, None, self.config.wcen_maxlambda) / self.cluster.scaleval) / self.config.wcen_pivot)

        # scale with richness for Poisson errors
        sig = self.config.lnw_cen_sigma / sigscale


        fw = gaussFunction(np.log(w),