def __interpolateParameters__(self, height, latitude, tkcDeep, tkcUpper):
     # Preallocate result array and start looping through all values
     isScalar = not util.isArray(height)
     results = []
     
     if isScalar:
         results = [0]
         height = [height]
         latitude = [latitude]
     else:
         results = np.zeros(height.shape)
     
     for i in range(0, len(height)):
         # Check where the height is with respect to the interpolation limits
         if height[i] <= tkcDeep[0][-1]:
             results[i] = scp_ip.splev(height[i], tkcDeep)
         elif height[i] >= tkcUpper[0][0]:
             results[i] = scp_ip.bisplev(height[i], latitude[i], tkcUpper)
         else:
             # Interpolate between the lower and upper interpolating functions (do
             # so linearly for now)
             low = scp_ip.splev(tkcDeep[0][-1], tkcDeep)
             high = scp_ip.bisplev(tkcUpper[0][0], latitude[i], tkcUpper)
             
             results[i] = low + (high - low) * (height[i] - tkcDeep[0][-1]) / \
                 (tkcUpper[0][0] - tkcDeep[0][-1])
                 
     if isScalar:
         return results[0]
         
     return results

    def apar_intrp(x,y,z):

        dx = dx_xy[x,y]
        dy = dy_xy[x,y]
        
        
        # Interpolating down the y-axis (along the field)
        y_vals = np.array(range(ny), dtype=float)

        for j in range(len(y_vals)):
            y_vals[j] = interpolate.bisplev(x,z,tck[j])
#        print y_vals, 'sdf'
        dy_coeffs = interpolate.splrep(range(ny), y_vals, k=3)
        
        # Interpolating along the slices of data
#        intrp = interpolate.RectBivariateSpline(range(nx),range(nz),data[:,y,:],kx=3,ky=3)
        
#        tx,ty = intrp.get_knots()
#        tck = (tx,ty,intrp.get_coeffs(),3,3)
#        print y, int(y), np.shape(data[:,y,:])
        # From the cubic spline coefficients, returns derivatives
#        print 's', int(np.rint(y)), int(y)
        dervs = ( interpolate.bisplev(x,z,tck[int(np.rint(y))], dx=1, dy=0)/dx,
                  interpolate.splev(y,dy_coeffs,der=1)/dy,
                  interpolate.bisplev(x,z,tck[int(np.rint(y))], dx=0, dy=1)/dz )

        return  dervs

    def make_pixel_lut(self, dims):
        """
        Generate an x and y image which maps the array indices into
        floating point array indices (to be corrected for pixel size later)

        returns 
        FIXME - check they are the right way around
                add some sort of known splinefile testcase
        """
        # Cache the value in case of multiple calls
        if self.pixel_lut is None:
            x_im = numpy.outer(range(dims[0]), numpy.ones(dims[1]))
            y_im = numpy.outer(numpy.ones(dims[0]), range(dims[1]))
            # xcor is tck2
            x_im = numpy.add( x_im,
                              bisplev( range(dims[1]),
                                               range(dims[0]),
                                               self.tck2 ).T,
                              x_im)
            # ycor is tck1
            y_im = numpy.add( y_im,
                              bisplev( range(dims[1]),
                                               range(dims[0]),
                                               self.tck1 ).T,
                              y_im)
            self.pixel_lut = x_im, y_im
        return self.pixel_lut

def transferFactorCalculator (rho, pt):

    # Scipy bi-linear spline representation data object.
    # Cf. [https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.interpolate.bisplrep.html]
    tck = [ np.array([    1.25,   1.25,   1.25,   1.25,    6.75,    6.75,    6.75,    6.75]), # Knots, x
            np.array([  450.,   450.,   450.,   450.,   1950.,   1950.,   1950.,   1950.]),   # Knots, y
            np.array([  0.18261346,  0.2174425 ,  0.37762574,  0.10284952,  0.26108651, # Spline coefficients
                       -0.14541588,  0.05701827,  0.27361263,  0.54531852,  0.77814774,
                        0.2479843 ,  0.23509468, -0.04597834, -0.47218929,  0.01928886,
                       -0.09066243]), 
            3, # Spline degree, x
            3] # Spline degree, y

    # Limits for the transfer factor map. Return zero if outside.
    limits = { 'rho': (  1.,    7.),
               'pt':  (400., 2000.) }

            # Check limits.
    if (limits['rho'][0] <= rho <= limits['rho'][1]) and \
       (limits['pt'] [0] <= pt  <= limits['pt'] [1]):
        # Calculate and return transfer factor.
        return interpolate.bisplev(rho, pt, tck)

    # Return fallback value.
    return 0.

 def fit(self, data, poldegree, swidth, sheight, threshold):
     if int(threshold) == -1:
         threshold = (int(data.mean()) * 10) / 7
     dims = data.shape
     xList = []
     yList = []
     zList = []
     for y in xrange(0, dims[0] - 1, sheight):
         for x in xrange(0, dims[1] - 1, swidth):
             view = data[y:y + sheight, x:x + swidth]
             flatIndex = numpy.argmax(view)
             yIdx, xIdx = numpy.unravel_index(flatIndex, view.shape)
             zValue = view[yIdx, xIdx]
             if zValue <= threshold:
                 xList.append(x + xIdx)
                 yList.append(y + yIdx)
                 zList.append(zValue)
     if len(xList) < (poldegree + 1) * (poldegree + 1):
         raise ValueError("Not enough reference points.")
     tck = interpolate.bisplrep(yList, xList, zList,
                                kx=poldegree, ky=poldegree,
                                xb=0, yb=0,
                                xe=int(dims[0]), ye=int(dims[1]))
     clipmin, clipmax = data.min(), threshold
     return interpolate.bisplev(range(dims[0]), range(dims[1]),
                                tck).clip(clipmin, clipmax)

    def interpolate(self, xi, yi):
        from scipy import interpolate
        # Need to write a norm function that calculates distance from a rib...

        """
        def interp(x1,x2,x3, x1i, x2i):
            spline = interpolate.Rbf(x1, x2, x3, function='thin-plate', smooth=0)
            return spline(x1i,x2i)

        try:
            zi = interp(self.points.x, self.points.y, self.points.z, xi, yi)
        except np.linalg.linalg.LinAlgError:
            zi = interp(self.points.y, self.points.x, self.points.z, yi, xi)
        """

        # Segfaults... Problems with the way scipy is compiled?
        tck = interpolate.bisplrep(self.points.x, self.points.y, self.points.z)
        zi = interpolate.bisplev(yi, xi, tck)


        """
        spline = interpolate.Rbf(self.points.x, self.points.y, self.points.z,
                                 function='thin-plate', smooth=0)
        zi = spline(xi,yi)
        """

        return zi

 def queryDepth(self, xq, yq):
     #return 1 data point
     depth = bisplev(xq, yq, self.tck)
     if np.isnan(depth):
         print "NaN returned for depth"
     else:
         return max(self.minz,min(depth, self.maxz))

def test_scipy_approx():
    """
    Test SciPy approximation of B-Spline surface
    :return: None
    """
    terrain_data = [
        [0.0, 0.0, 0.0], [0.0, 0.5, 0.4], [0.0, 1.0, 0.0],
        [0.5, 0.0, 0.2], [0.5, 0.5, 0.8], [0.5, 1.0, 0.2],
        [1.0, 0.0, 0.0], [1.0, 0.5, 0.4], [1.0, 1.0, 0.0]]
    tX = [item[0] for item in terrain_data]
    tY = [item[1] for item in terrain_data]
    tZ = [item[2] for item in terrain_data]
    from scipy import interpolate
    print('SciPy approximation ...')
    start_time = time.time()
    tck, fp, ior, msg = interpolate.bisplrep(tX, tY, tZ, kx=2, ky=2, full_output=1)
    end_time = time.time()
    print('Computed in {0} seconds.'.format(end_time - start_time))
    occ_bspline = convert.bspline.scipy_to_occ(tck)
    # Compute difference between original terrain data and B-Spline surface
    u_num = v_num = 50
    points = [[float(i)/u_num, float(j)/u_num, 0.0] for i in range(v_num+1) for j in range(u_num+1)]
    points = [(it[0], it[1], interpolate.bisplev(it[0], it[1], tck)) for it in points]
    # points = terrain_data

    display_results(occ_bspline, points)

 def evaluatePartialDerivativeV(self, x, y):
     result = np.empty(self.coeffElems)
     
     for i in range(self.coeffElems):
         result[i] = interpolate.bisplev(x, y, self.tcks[i], dy=1)
         
     return result

 def evaluate(self, x, y):
     result = np.empty(self.coeffElems)
     
     for i in range(self.coeffElems):
         result[i] = interpolate.bisplev(x, y, self.tcks[i])
         
     return result

    def derivatives(self, alpha, Re):

        # note: direct call to bisplev will be unnecessary with latest scipy update (add derivative method)
        tck_cl = self.cl_spline.tck[:3] + self.cl_spline.degrees  # concatenate lists
        tck_cd = self.cd_spline.tck[:3] + self.cd_spline.degrees

        dcl_dalpha = bisplev(alpha, Re, tck_cl, dx=1, dy=0)
        dcd_dalpha = bisplev(alpha, Re, tck_cd, dx=1, dy=0)

        if self.one_Re:
            dcl_dRe = 0.0
            dcd_dRe = 0.0
        else:
            dcl_dRe = bisplev(alpha, Re, tck_cl, dx=0, dy=1)
            dcd_dRe = bisplev(alpha, Re, tck_cd, dx=0, dy=1)

        return dcl_dalpha, dcl_dRe, dcd_dalpha, dcd_dRe

def loggtracks(masslimit,location,fileroot, metalstr, plot=True):



    #This is an array of all the masses that are part of the filenames
    massarrstr=['120.0', '85.0', '60.0', '40.0', '25.0', '20.0', '15.0', '12.0', '9.0', '7.0', '5.0', '4.0', '3.0', '2.5', '2.0', '1.7', '1.5', '1.25', '1.0', '0.9']



    #this is for each mass, read in the file and plot the evolutionary tracks
    for imass in range(masslimit):

        filemass = massarrstr[imass]

        filename = location+fileroot+filemass

        DataIn = np.genfromtxt(filename, dtype="float", unpack=True)
    
        time = DataIn[3,:]
        timesec = time*365*24*60*60
        Teff = 10**DataIn[ 6,:]
        logTeff = [math.log10(jj) for jj in Teff]
        L = 10**DataIn[ 5,:]*Lsun
        Rsquared = (L/(4*math.pi*sigma*Teff**4))
        Mass = DataIn[4,:]*Msun
        surfaceGrav = Mass*G/(Rsquared)
        logsurfaceGrav = [math.log10(ii) for ii in surfaceGrav]


        if plot == True:

            fadeplot(Teff, logsurfaceGrav, time)
        
        q0s = np.zeros(len(Teff))
        for timestep in range(len(Teff)):
            q0s[timestep] = interpolate.bisplev(Teff[timestep],logsurfaceGrav[timestep],gridq0s)
            
        
        
        for kk in range(len(q0s)):
            if q0s[kk] < 0:
                q0s[kk] = -q0s[kk]
                q0s[kk] = np.log10(q0s[kk])


        totalQ0s = np.trapz(q0s,timesec)
        
        logtotalQ0s = np.log10(totalQ0s)
        #print "Mass of: " + massarr[imass]+" Produces "+str(totalQ0s)+" Photons"
        logq0ints[imass] = totalQ0s
        #plt.plot(timesec/timesec[-1], q0s, 'k')
        titlestr= "Z= "+metalstr
        #plt.title(titlestr)
        #plt.xlabel("Stellar Lifetime")
        #plt.ylabel("Photons / Second")
        #plt.ylim([0,2e50])
    #plt.show()    
    return;

    def aproximate_terrain(self):
        """
        Try to aproximate terrain with bspline surface
        """

        tck,fp,ior,msg = interpolate.bisplrep(self.tX, self.tY, self.tZ, kx=5, ky=5, full_output=1)
        self.tck[(self.min_x, self.min_y, self.max_x, self.max_y)] = tck
        # Compute difference between original terrain data and b-spline surface
        self.tW = [abs(it[2] - interpolate.bisplev(it[0], it[1], tck)) for it in self.terrain_data]

def interpgrid(x,y, xlist,ylist, xmap, ymap, kx=3, ky=3, s=50):
   ''' for position x,y and a 2-D mapping map(list),
       i.e., xmap[xlist,ylist],ymap[xlist,ylist] given on a grid xlist,ylist; 
       the nearest xlist, ylist positions to each x,y pair are found and 
       interpolated to yield  mapx(x,y),mapy(x,y)
         
   x,y : rank-1 arrays of data points
   xlist, ylist, xmap, ymap: rank-1 arrays of data points
   
   +
   2008-08-24 NPMK (MSSL)
   '''
   from scipy import interpolate
   # check if the input is right data type
   # ... TBD
   
   # compute the Bivariate-spline coefficients
   # kx = ky =  3 # cubic splines (smoothing)
   task = 0 # find spline for given smoothing factor
   # s = 50 # spline goes through the given points
   # eps = 1.0e-6  (0 < eps < 1)
   
   #(tck_x, ems1) 
   tck_x = interpolate.bisplrep(xlist,ylist,xmap,kx=kx,ky=ky,s=s)
   #(fp1, ier1, msg1) = ems1
   #if ier1 in [1,2,3]: 
   #   print 'an error occurred computing the bivariate spline (xmap) '
   #   print ier1, msg1
   #   # raise error
   #   return None
   tck_y = interpolate.bisplrep(xlist,ylist,ymap,kx=kx,ky=ky,s=s) 
   #(fp2, ier2, msg2) = ems2
   #if ier2 in [1,2,3]: 
   #   print 'an error occurred computing the bivariate spline (ymap) '
   #   print ier2, msg2
   #   # raise error
   #   return None
   # compute the spline    
   
   xval = interpolate.bisplev(x, y, tck_x)
   yval = interpolate.bisplev(x, y, tck_y)
   
   return xval,yval

def sample(sgridxy, im, dd=(0,0), hscaling=1. ):
    xs = sgridxy[0].flatten()
    ys = sgridxy[1].flatten()
    assert(xs.size == ys.size)
    v = np.zeros(xs.size)
    for i in np.arange(v.size):
        v[i] = (1./(hscaling**np.sum(dd)))*flt(interpolate.bisplev(xs[i],ys[i], im, dd[0], dd[1]) )

    #return v.reshape(sgridxy[0].shape)
    return v

def plt(n=25):
	x=[]
	y=[]
	qx=[]
	qy=[]
	for i in range(n):
		x.append(r())
		y.append(r())
		qx.append(sin(x[-1]))
		qy.append(cos(y[-1]))
	qxb=bisplrep(x,y,qx,s=0)
	qyb=bisplrep(x,y,qy,s=0)
	X=arange(-2,2,0.4)
	Y=arange(-2,2,0.4)
	cla()
	hold(True)
	quiver(x,y,qx,qy,pivot='tail',color='b')
	quiver2(X,Y,bisplev(X, Y,qxb),bisplev(X, Y,qyb),pivot='tail',color='r')
	hold(False)

def processing(filename,x_u,y_u,x_l,y_l,s_val,x_new_res,y_new_res,coord_opt,contour_lim):
	#load in data as a 2D matrix
	try:
   		with open(filename): pass
	except IOError:
  		return -1
	values = np.loadtxt(filename,delimiter=',')

	#Check if 95% limit will exist
	flag = False
	for row in values:
		for element in row:
			if element >= contour_lim:
				flag = True
				break
	if (flag == False):
		return -2
	
	#define data co-ordinates
	#TODO: take into account irregularly spaced data values
	if coord_opt == 'd':
		x = np.mgrid[x_l:x_u:len(values[0])*1j]
		y = np.mgrid[y_l:y_u:len(values)*1j]
	elif coord_opt == 'n':
	#request to read in co-ordinates noted in data file
		try:
   			with open(filename+"_coord"): pass
		except IOError:
  			return -3
  		else:
  			filename_coord = filename+"_coord"
  			data_coord=open(filename_coord)
  			x=((data_coord.readline()).strip()).split(',')
  			x = [float(i) for i in x ]
  			y=((data_coord.readline()).strip()).split(',')
  			y = [float(i) for i in y ]
	
	x,y = np.meshgrid(x,y)
	#interpolate using quadratic splines
	#Quadratic are used to better preserve asymptotic nature of plots
	#TODO:What value of s is optimal?
	tck = interp.bisplrep(x,y,values,kx=2,ky=2,s=s_val)
	
	#define points to interpolate over
	xnew,ynew = np.mgrid[x_l:x_u:(x_new_res*1j),y_l:y_u:(y_new_res*1j)]
	values_new = interp.bisplev(xnew[:,0],ynew[0,:],tck)

	#plot only the cls_level line
	v=np.linspace(contour_lim,contour_lim,2)
	cs = plt.contour(xnew,ynew,values_new,v)
	
	#Extract data of cls_level line
	#TODO: investigate syntax of this line
	#TODO: catch error where there is data below 95% but not enough to generate a contour
	return (cs.collections[0].get_paths()[0]).vertices

def finterp(band, t, p, param, gen, extrap=False):
   '''interpolate at time t and param p for param,gen combo.'''
   load_data(band,param,gen)
   if param == 'dm15':
      f = dm15_flux[(band,gen)]
      ef = dm15_eflux[(band,gen)]
   else:
      f = st_flux[(band,gen)]
      ef = st_eflux[(band,gen)]

   if len(num.shape(t)) == 0:
      scalar = 1
   else:
      scalar = 0
   t = num.atleast_1d(t)
   # First the evaluation mtarix:
   Z = num.atleast_2d(bisplev(t, p, f))[:,0]
   eZ = num.atleast_2d(bisplev(t, p, ef))[:,0]
   if not extrap:
      mask = num.greater_equal(t,f[0][0])*num.less_equal(t,f[0][-1])
      mask = mask*num.greater(Z, 0)
      Z = num.where(mask, Z, 1)
      eZ = num.where(mask, eZ, -1)
   else:
      t1,t2 = get_t_lim(band, param, gen)
      mask = -num.isnan(Z)
      # extrapolate lower with t^2 law
      if num.sometrue(num.less(t,t1)):
         Tp = bisplev(t1, p, f, dx=1)
         T = bisplev(t1, p, f)
         eT = bisplev(t, p, ef)
         t0 = t1 - 2*T/Tp; a = T/(t1-t0)**2
         Z = num.where(num.less(t, t1), a*num.power(t-t0,2), Z)
         eZ = num.where(num.less(t, t1), eT, eZ)
         mask = mask*num.greater(Z,0)*num.greater(t, t0)
      if num.sometrue(num.greater(t, t2)):
         # extrapolate with m = a*(t-t2)+b
         Tp = bisplev(t2, p, f, dx=1)
         T = bisplev(t2, p, f)
         eT = bisplev(t2, p, ef)
         b = -2.5*num.log10(T)
         a = -2.5/num.log(10)/T*Tp
         f = num.power(10, -0.4*(a*(t-t2)+b))
         Z = num.where(num.greater(t, t2), f, Z)
         eZ = num.where(num.greater(t, t2), eT, eZ)
      Z = num.where(mask, Z, 1)
   if scalar:
      return Z[0],eZ[0],mask[0]
   else:
      return Z,eZ,mask

def get_splined_2d_dist(fname, islog=False, num_spline_points=100):
    data = np.loadtxt(fname, skiprows=1)
    if islog:
        lnL = data[:, 2]
    else:
        lnL = np.log(data[:, 2])
    lnL = -2* (lnL - np.max(lnL))
    tck = interpolate.bisplrep(data[:, 0], data[:, 1], lnL, s=1)
    num_spline_points = complex(0, num_spline_points)
    Q_args, N_args = np.mgrid[data[0, 0]:data[-1, 0]:num_spline_points, data[0, 1]:data[-1, 1]:num_spline_points]
    lnL_splined = interpolate.bisplev(Q_args[:, 0], N_args[0, :], tck)

    return Q_args, N_args, lnL_splined

	def createDatasetForGene(self, gene_ind, plot = False):
		if gene_ind not in [3,4,5,6,7]:
			raise Exception("Wrong gene")
		'''use only wt data for now'''
		data = self.dp.normData[:,:,0,:]
		x_range = np.linspace(0, data.shape[2]-1, data.shape[2])
		t_range = np.linspace(0, data.shape[0]-1, data.shape[0])
		xv, tv = np.meshgrid(x_range, t_range)
		x = xv.flatten()
		t = tv.flatten()
		z = data[:,gene_ind,:].flatten()
		spdat = ip.bisplrep(x,t,z,s=5)
		t_der = ip.bisplev(x_range, t_range, spdat, dx=0, dy=1)
		x_der2 = ip.bisplev(x_range, t_range, spdat, dx=2, dy=0)
		input_list = []
		for g in xrange(7):
			input_list.append(data[:,g,:].flatten())
		input_list.append(x_der2.T.flatten())
		input_list = np.rollaxis(np.array(input_list), 1, 0)
		output_list, self.omax, self.omin = self.normalize(t_der.T.flatten(), -0.9, 0.9)
		
		if plot is True:
			fig = plt.figure()
			ax = fig.add_subplot(221, projection='3d')
			ax.plot_surface(xv, tv, t_der.T)
			ax = fig.add_subplot(222, projection='3d')
			ax.plot_surface(xv, tv, x_der2.T)
			ax = fig.add_subplot(223, projection='3d')
			x_range = np.linspace(0, data.shape[2]-1, 200)
			t_range = np.linspace(0, data.shape[0]-1, 200)
			xv, tv = np.meshgrid(x_range, t_range)
			plt_data = ip.bisplev(x_range, t_range, spdat)
			ax.plot_surface(xv, tv, plt_data.T)
			ax = fig.add_subplot(224)
			ax.hist(t_der.flatten(), bins=40)
			plt.show()
			exit()
		
		return input_list, output_list