def execute(self):

        c12axis = self.calculate_c12axis()

        if self.interp_from_htc:
            cni = self.c12axis_init.shape[0]
            if self.c12axis_init[-1, 2] > 1.:
                self.c12axis_init /= self.blade_length

            # interpolate blade_ae distribution onto c12 distribution
            self.blade_ae.c12axis = np.zeros((cni, 4))
            for i in range(4):
                tck = pchip(c12axis[:, 2], c12axis[:,i])
                self.blade_ae.c12axis[:, i] = tck(self.c12axis_init[:, 2])
        else:
            ds_root = 1. / self.blade_ni_span
            ds_tip = 1. / self.blade_ni_span / 3.
            dist = np.array([[0., ds_root, 1],
                             [1., ds_tip, self.blade_ni_span]])
            x = distfunc(dist)
            self.blade_ae.c12axis = np.zeros((x.shape[0], 4))
            for i in range(4):
                tck = pchip(c12axis[:, 2], c12axis[:,i])
                self.blade_ae.c12axis[:, i] = tck(x)

        # scale main axis according to radius
        self.blade_ae.c12axis[:,:3] *= self.blade_length

        self.blade_ae.radius = self.blade_length + self.hub_radius
        self.blade_ae.s = self.bladegeom.smax * self.bladegeom.s * self.blade_length
        self.blade_ae.rthick = self.bladegeom.rthick * 100.
        self.blade_ae.chord = self.bladegeom.chord * self.blade_length
        self.blade_ae.aeset = np.ones(len(self.blade_ae.s))

 def test_endslopes(self):
     # this is a smoke test for gh-3453: PCHIP interpolator should not
     # set edge slopes to zero if the data do not suggest zero edge derivatives
     x = np.array([0.0, 0.1, 0.25, 0.35])
     y1 = np.array([279.35, 0.5e3, 1.0e3, 2.5e3])
     y2 = np.array([279.35, 2.5e3, 1.50e3, 1.0e3])
     for pp in (pchip(x, y1), pchip(x, y2)):
         for t in (x[0], x[-1]):
             assert_(pp(t, 1) != 0)

    def test_cast(self):
        # regression test for integer input data, see gh-3453
        data = np.array([[0, 4, 12, 27, 47, 60, 79, 87, 99, 100], [-33, -33, -19, -2, 12, 26, 38, 45, 53, 55]])
        xx = np.arange(100)
        curve = pchip(data[0], data[1])(xx)

        data1 = data * 1.0
        curve1 = pchip(data1[0], data1[1])(xx)

        assert_allclose(curve, curve1, atol=1e-14, rtol=1e-14)

def bdrate(file1, file2, anchorfile):
    if anchorfile:
        anchor = flipud(loadtxt(anchorfile));
    a = flipud(loadtxt(file1));
    b = flipud(loadtxt(file2));
    rates = [0.06,0.2];
    qa = a[:,0]
    qb = b[:,0]
    ra = a[:,2]*8./a[:,1]
    rb = b[:,2]*8./b[:,1]
    bdr = zeros((4,4))
    ret = {}
    for m in range(0,len(met_index)):
        try:
            ya = a[:,3+m];
            yb = b[:,3+m];
            if anchorfile:
                yr = anchor[:,3+m];
            #p0 = interp1d(ra, ya, interp_type)(rates[0]);
            #p1 = interp1d(ra, ya, interp_type)(rates[1]);
            if anchorfile:
                p0 = yr[0]
                p1 = yr[-1]
            else:
                minq = 20
                maxq = 55
                try:
                    minqa_index = qa.tolist().index(minq)
                    maxqa_index = qa.tolist().index(maxq)
                    minqb_index = qb.tolist().index(minq)
                    maxqb_index = qb.tolist().index(maxq)
                except ValueError:
                    q_not_found = True
                    minqa_index = -1
                    maxqa_index = 0
                    minqb_index = -1
                    maxqb_index = 0
                p0 = max(ya[maxqa_index],yb[maxqb_index])
                p1 = min(ya[minqa_index],yb[minqb_index])
            a_rate = pchip(ya, log(ra))(arange(p0,p1,abs(p1-p0)/5000.0));
            b_rate = pchip(yb, log(rb))(arange(p0,p1,abs(p1-p0)/5000.0));
            if not len(a_rate) or not len(b_rate):
                bdr = NaN;
            else:
                bdr=100 * (exp(mean(b_rate-a_rate))-1);
        except ValueError:
            bdr = NaN
        except linalg.linalg.LinAlgError:
            bdr = NaN
        except IndexError:
            bdr = NaN
        if abs(bdr) > 1000:
            bdr = NaN
        ret[m] = bdr
    return ret

    def interpolate(self):
        pitch = np.zeros((self.nframes))
        pitch[:] = self.samp_values
        pitch2 = medfilt(self.samp_values, self.SMOOTH_FACTOR)

        # This part in the original code is kind of confused and caused
        # some problems with the extrapolated points before the first
        # voiced frame and after the last voiced frame. So, I made some
        # small modifications in order to make it work better.
        edges = self.edges_finder(pitch)
        first_sample = pitch[0]
        last_sample = pitch[-1]

        if len(np.nonzero(pitch2)[0]) < 2:
            pitch[pitch == 0] = self.PTCH_TYP
        else:
            nz_pitch = pitch2[pitch2 > 0]
            pitch2 = scipy_interp.pchip(np.nonzero(pitch2)[0],
                                        nz_pitch)(range(self.nframes))
            pitch[pitch == 0] = pitch2[pitch == 0]
        if self.SMOOTH > 0:
            pitch = medfilt(pitch, self.SMOOTH_FACTOR)
        try:
            if first_sample == 0:
                pitch[:edges[0]-1] = pitch[edges[0]]
            if last_sample == 0:
                pitch[edges[-1]+1:] = pitch[edges[-1]]
        except:
            pass
        self.samp_interp = pitch

def pchipPlot(X, Y, x_new, **args):
    title = 'pchip_plot'
    if(args.has_key('wn')):
        if(type(args['wn']) == int):
	    title = title+str(args['wn'])
            plt.figure(title)
    else:
        plt.figure(title)
    pc = pchip(X,Y)
    plt.plot(X,Y, linestyle = '-')
    plt.plot(x_new, pc(x_new), linestyle = '--')
    plt.legend(['original', 'Spline pchip '], loc='best')
    plt.ylabel(r'epsilon', size = 12)
    a=plt.gca()
    a.set_yscale('log')
    a.set_xscale('log')
    #Here Im going to save the file on disk if the
    #user want and show the plot on the screen
    if(args.has_key('save')):
        if(type(args['save']) == bool and args['save'] ==  True):
            plt.savefig(title.strip()+'.png')
            plt.draw()
    if(args.has_key('view')):
        if(type(args['view']) == bool and args['view'] == True):
            plt.show()

    def solve_nonlinear(self, params, unknowns, resids):

        # we need to dig into the _ByObjWrapper val to get the array
        # values out
        # pf_in = {name: val['val'].val for name, val in params.iteritems()}
        pf_in = {}
        pf_in['s'] = params['s']
        pf_in['x'] = params['x']
        pf_in['y'] = params['y']
        pf_in['z'] = params['z']
        pf_in['rot_x'] = params['rot_x']
        pf_in['rot_y'] = params['rot_y']
        pf_in['rot_z'] = params['rot_z']
        pf_in['chord'] = params['chord']
        pf_in['rthick'] = params['rthick']
        pf_in['p_le'] = params['p_le']

        if _PGL_installed:
            pf = redistribute_planform(pf_in, s=self.s_new, spline_type=self.spline_type)
        else:
            pf = {}
            for k, v in pf_in.iteritems():
                spl = pchip(pf_in['s'], v)
                pf[k] = spl(self.s_new)

        for k, v in pf.iteritems():
            unknowns[k+self._suffix] = v
        unknowns['athick'+self._suffix] = pf['chord'] * pf['rthick']

    def solve_nonlinear(self, params, unknowns, resids):

        # we need to dig into the _ByObjWrapper val to get the array
        # values out
        # pf_in = {name: val['val'].val for name, val in params.iteritems()}
        pf_in = {}
        pf_in["s"] = params["s"]
        pf_in["x"] = params["x"]
        pf_in["y"] = params["y"]
        pf_in["z"] = params["z"]
        pf_in["rot_x"] = params["rot_x"]
        pf_in["rot_y"] = params["rot_y"]
        pf_in["rot_z"] = params["rot_z"]
        pf_in["chord"] = params["chord"]
        pf_in["rthick"] = params["rthick"]
        pf_in["p_le"] = params["p_le"]

        if _PGL_installed:
            pf = redistribute_planform(pf_in, s=self.s_new, spline_type=self.spline_type)
        else:
            pf = {}
            for k, v in pf_in.iteritems():
                spl = pchip(pf_in["s"], v)
                pf[k] = spl(self.s_new)

        for k, v in pf.iteritems():
            unknowns[k + self._suffix] = v
        unknowns["athick" + self._suffix] = pf["chord"] * pf["rthick"]

    def execute(self):

        if self.interp_from_htc:
            c12axis = self.c12axis_init.copy()
        else:
            c12axis = self.calculate_c12axis()

        ds_root = 1. / self.blade_ni_span
        ds_tip = 1. / self.blade_ni_span / 3.
        dist = np.array([[0., ds_root, 1], [1., ds_tip, self.blade_ni_span]])

        x = distfunc(dist)
        self.c12axis = np.zeros((x.shape[0], 4))
        for i in range(4):
            tck = pchip(c12axis[:, 2], c12axis[:, i])
            self.c12axis[:, i] = tck(x)

        # scale main axis according to radius
        self.c12axis[:, :3] *= self.blade_length

        l = ((self.c12axis[1:, 0]-self.c12axis[:-1, 0])**2 +
             (self.c12axis[1:, 1]-self.c12axis[:-1, 1])**2 +
             (self.c12axis[1:, 2]-self.c12axis[:-1, 2])**2)**.5

        self.blade_ae.s = self.bladegeom.s * sum(l) / self.bladegeom.s[-1]
        self.blade_ae.rthick = self.bladegeom.rthick * 100.
        self.blade_ae.chord = self.bladegeom.chord * self.blade_length
        self.blade_ae.aeset = np.ones(len(self.blade_ae.s))

    def test_nag(self):
        # Example from NAG C implementation,
        # http://nag.com/numeric/cl/nagdoc_cl25/html/e01/e01bec.html
        # suggested in gh-5326 as a smoke test for the way the derivatives
        # are computed (see also gh-3453)
        from scipy._lib.six import StringIO
        dataStr = '''
          7.99   0.00000E+0
          8.09   0.27643E-4
          8.19   0.43750E-1
          8.70   0.16918E+0
          9.20   0.46943E+0
         10.00   0.94374E+0
         12.00   0.99864E+0
         15.00   0.99992E+0
         20.00   0.99999E+0
        '''
        data = np.loadtxt(StringIO(dataStr))
        pch = pchip(data[:,0], data[:,1])

        resultStr = '''
           7.9900       0.0000
           9.1910       0.4640
          10.3920       0.9645
          11.5930       0.9965
          12.7940       0.9992
          13.9950       0.9998
          15.1960       0.9999
          16.3970       1.0000
          17.5980       1.0000
          18.7990       1.0000
          20.0000       1.0000
        '''
        result = np.loadtxt(StringIO(resultStr))
        assert_allclose(result[:,1], pch(result[:,0]), rtol=0., atol=5e-5)

 def test_two_points(self):
     # regression test for gh-6222: pchip([0, 1], [0, 1]) fails because
     # it tries to use a three-point scheme to estimate edge derivatives,
     # while there are only two points available.
     # Instead, it should construct a linear interpolator.
     x = np.linspace(0, 1, 11)
     p = pchip([0, 1], [0, 2])
     assert_allclose(p(x), 2*x, atol=1e-15)

def splinePlot(p, X, Y, cl, l):
  p.plot(X, Y, 'o', mfc=cl, c=cl)
  #f = interp1d(X, Y, kind='pchip')
  #f = splrep(X, Y, s=0)
  f = pchip(X, Y)
  nx = np.linspace(X[0], X[-1], 100)
  ny = [f(x) for x in nx]
  p.plot(nx, ny, '-', color=cl, label=l)

 def __init__(self, dz, w1, w2, w3, a, b):
     # INFO 1kcal = 4184e20 kg A^2 s^-2
     # TODO check bounds of a, b
     self.dz = dz
     self.x = np.array([-dz, -dz+1, -a, -b, 0, b, a, dz-1, dz])
     self.y = np.array([0, 0, w1, w2, w3, w2, w1, 0, 0])
     self.pmf = pchip(self.x, self.y)
     self.derivative = self.pmf.derivative(1)

 def __call__(self, z):
     """Parameters: z is a number, sequence or array.
     This method makes an instance f of LinInterp callable,
     so f(z) returns the interpolation value(s) at z.
     """
     if self.kind == 'pchip':
         return pchip(self.X, self.Y)(z)
     else:
         return interp1d(self.X, self.Y, kind=self.kind,
                         bounds_error=False)(z)

    def execute(self):

        self.pfOut.s = self.x.copy()
        self.pfOut.blade_length = self.pfIn.blade_length
        self.pfIn._compute_s()
        for name in self.pfIn.list_vars():
            var = getattr(self.pfIn, name)
            if not isinstance(var, np.ndarray): continue
            tck = pchip(self.pfIn.s, var)
            newvar = tck(self.x) 
            setattr(self.pfOut, name, newvar)

    def test_all_zeros(self):
        x = np.arange(10)
        y = np.zeros_like(x)

        # this should work and not generate any warnings
        with warnings.catch_warnings():
            warnings.filterwarnings('error')
            pch = pchip(x, y)

        xx = np.linspace(0, 9, 101)
        assert_equal(pch(xx), 0.)

def itpl(PTS,numPTS,degree,sID,eID):
    t=range(len(PTS))
    ipl_t=np.linspace(sID,eID,numPTS*(eID-sID)/(len(PTS)-1)+1)
    newX=pchip(t,[row[0] for row in PTS])(ipl_t)
    newY=pchip(t,[row[1] for row in PTS])(ipl_t)
    newZ=pchip(t,[row[2] for row in PTS])(ipl_t)

    if degree==1:
        newX=interp1d(t,[row[0] for row in PTS],kind='slinear')(ipl_t)
        newY=interp1d(t,[row[1] for row in PTS],kind='slinear')(ipl_t)
        newZ=interp1d(t,[row[2] for row in PTS],kind='slinear')(ipl_t)
    elif degree==2:
        newX=interp1d(t,[row[0] for row in PTS],kind='quadratic')(ipl_t)
        newY=interp1d(t,[row[1] for row in PTS],kind='quadratic')(ipl_t)
        newZ=interp1d(t,[row[2] for row in PTS],kind='quadratic')(ipl_t)
    elif degree==3:
        newX=interp1d(t,[row[0] for row in PTS],kind='cubic')(ipl_t)
        newY=interp1d(t,[row[1] for row in PTS],kind='cubic')(ipl_t)
        newZ=interp1d(t,[row[2] for row in PTS],kind='cubic')(ipl_t)
    elif degree==4:
        newX=pchip(t,[row[0] for row in PTS])(ipl_t)
        newY=pchip(t,[row[1] for row in PTS])(ipl_t)
        newZ=pchip(t,[row[2] for row in PTS])(ipl_t)
    else:
        newX=interp1d(t,[row[0] for row in PTS],kind='quadratic')(ipl_t)
        newY=interp1d(t,[row[1] for row in PTS],kind='quadratic')(ipl_t)
        newZ=interp1d(t,[row[2] for row in PTS],kind='quadratic')(ipl_t)

    outPTS=[]
    for i in range(0,len(newX)):
        outPTS.append([newX[i],newY[i],newZ[i]])

    return outPTS

    def __call__(self, x, Cx, C):
        """
        params:
        ----------
        x: array
            array with new x-distribution
        xp: array
            array with x-coordinates of spline control points
        yp: array
            array with y-coordinates of spline control points

        returns
        ---------
        ynew: array
            resampled points
        """
        spl = pchip(Cx, C)
        return spl(x)

    def execute(self):
        self.dict = HAWC2InputDict()
        self.dict.read(self.htc_master_file)
        self.htc = self.dict.htc
        self.vartrees.body_order = self.dict.body_order

        for section in self.htc:
            if section.name == 'simulation':
                self.add_simulation(section)
            elif section.name == 'wind':
                self.add_wind(section)
            elif section.name == 'aero':
                self.add_aero(section)
            elif section.name == 'aerodrag':
                self.add_aerodrag(section)
            elif section.name == 'new_htc_structure':
                self.add_structure(section)
            elif section.name == 'output':
                self.add_output(section)
            elif section.name == 'dll':
                self.add_dlls(section)
            elif section.name == 'hawcstab2':
                self.add_hawcstab2(section)

        # count number of blades
        for iblade in range(1, 10):
            if 'blade'+str(iblade) not in self.vartrees.body_order:
                self.vartrees.rotor.nblades = iblade-1
                break

        # copy blade twist from c2_def to blade_ae vartree
        if hasattr(self.vartrees.main_bodies, 'blade1'):
            from scipy.interpolate import pchip
            c12 = self.vartrees.main_bodies.blade1.c12axis
            tck = pchip(c12[:, 2], c12[:, 3])
            twist = tck(self.vartrees.blade_ae.s)
            self.vartrees.blade_ae.twist = twist
            self.vartrees.blade_ae.c12axis = c12.copy()
            self.vartrees.blade_structure = self.vartrees.main_bodies.blade1.beam_structure

# List of Earth masses converted to M_Jup
x_new = toJupMass(MEarth)

# np.linspace(min,max,n) creates a list of n values between min and max
jupX = np.linspace(min(MJup),max(MJup),500)
earthX = np.linspace(min(x_new),max(x_new),500)

# convert radius and masses to a density
IceD = toDensity(MEarth, IceR)
RockD = toDensity(MEarth, RockR)
IronD = toDensity(MEarth, IronR)

lineColor = 'LightGray'

# Interpolate a line between the data points
interp = pchip(np.array(x_new),np.array(IceD))
ax.loglog(earthX, interp(earthX), c=lineColor, lw=4.0)
interp = pchip(np.array(x_new),np.array(RockD))
ax.loglog(earthX, interp(earthX), c=lineColor, lw=4.0)
interp = pchip(np.array(x_new),np.array(IronD))
ax.loglog(earthX, interp(earthX), c=lineColor, lw=4.0)

# Convert Jupiter radii and masses to density
H_HE_D = []
for idx in range(0,len(MJup)):
   g_mass = MJup[idx] * 1.898e30
   cm_rad = H_HE[idx] * 6.9911e9
   density = g_mass/( (4.0/3.0) * constants.pi * (cm_rad**3))
   H_HE_D.append(density)

interp = pchip(np.array(MJup),np.array(H_HE_D))