# Path: redmapper/utilities.py
# 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]


import numpy as np
import scipy.optimize
import esutil
import warnings
from scipy import special
from .utilities import CubicSpline, interpol
            if self._has_err_ratios:
                retval.append(pars[self._scatter_index: self._scatter_index + self._n_scatter_nodes + 2])
            else:
                retval.append(pars[self._scatter_index: self._scatter_index + self._n_scatter_nodes])

        return retval

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

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

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

        if self._fit_slope:
            # We are fitting the slope

# Path: redmapper/utilities.py
# 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]


# Path: redmapper/fitters.py
import numpy as np
import scipy.optimize
import esutil
import warnings
from scipy import special
from .utilities import CubicSpline, interpol

            if self._has_err_ratios:
                retval.append(pars[self._scatter_index: self._scatter_index + self._n_scatter_nodes + 2])
            else:
                retval.append(pars[self._scatter_index: self._scatter_index + self._n_scatter_nodes])

        return retval

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

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

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

        if self._fit_slope:
            # We are fitting the slope

            spl = CubicSpline(self._slope_nodes, pars[self._slope_index: self._slope_index + self._n_slope_nodes])