def fit(hm, t=None):
    X, Y = np.meshgrid(np.arange(0.0, hm.shape[1], 1.0), np.arange(0.0, hm.shape[0], 1.0))

    if t is not None:
        return interpolate.bisplrep(X, Y, hm, kx=3, ky=3, task=-1, tx=t[0], ty=t[1])
    else:
        return interpolate.bisplrep(X, Y, hm, kx=3, ky=3)

    def __init__(self, tables, smooth=0.05, **params):
        logging.info("Loading atmospheric flux table %s" %tables)
        
        #Load the data table
        table = np.loadtxt(open_resource(tables)).T

        #columns in Honda files are in the same order 
        cols = ['energy']+primaries 

        flux_dict = dict(zip(cols, table))
        for key in flux_dict.iterkeys():
            
            #There are 20 lines per zenith range
            flux_dict[key] = np.array(np.split(flux_dict[key], 20))
            if not key=='energy':
                flux_dict[key] = flux_dict[key].T
    
        #Set the zenith and energy range 
        flux_dict['energy'] = flux_dict['energy'][0]
        flux_dict['coszen'] = np.linspace(0.95, -0.95, 20)
    
        #Now get a spline representation of the flux table.
        logging.debug('Make spline representation of flux')
        # do this in log of energy and log of flux (more stable)
        logE, C = np.meshgrid(np.log10(flux_dict['energy']), flux_dict['coszen'])

        self.spline_dict = {}
        for nutype in primaries:
            #Get the logarithmic flux
            log_flux = np.log10(flux_dict[nutype]).T
            #Get a spline representation
            spline =  bisplrep(logE, C, log_flux, s=smooth)
            #and store
            self.spline_dict[nutype] = spline

    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 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 fitspline(xy, z, smoothing=None):
    """
    Return a tuple (t, c, k) describing the spline as described in the
    Numpy documentation.
    """
    return bisplrep(xy[...,0].ravel(), xy[...,1].ravel(), z.ravel(),
                    s=smoothing)

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 get_tck(x,y,z,kk_default=3,logger=None):
        logger = logger or logging.getLogger(__name__)

        # estimate max number of allowed spline order
        mm = len(x)
        kk_max = np.int(np.floor(np.sqrt(mm/2.0)-1))
        kk = np.min([kk_default,kk_max])

        result = bisplrep(x=x,y=y,z=z,kx=kk,ky=kk,full_output=1)
        if result[2]>0:
            logger.info("Interpolation problem:%s" % result[-1])
            logger.info("Now, try to adjust s")
            result = bisplrep(x=x,y=y,z=z,kx=kk,ky=kk,s=result[1],full_output=1)
            if result[2]>0:
                raise ValueError("Interpolation problem:%s" % result[-1])
        return result[0]

    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 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 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 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 init_spline(self,dhalo,psi,z):
        """Compute knots and coefficients of an interpolating spline
        given a grid of points in halo distance (dhalo) and offset
        angle (psi) at which the LoS integral has been computed.
        """

        kx = 2
        ky = 2
        self._psi_min = psi.min()
        self._tck = bisplrep(dhalo,psi,np.log10(z),s=0.0,kx=kx,ky=ky,
                             nxest=int(kx+np.sqrt(len(z.flat))),
                             nyest=int(ky+np.sqrt(len(z.flat))))

def process_functions(functions, UV_data, nodes):
    """
    Processes coefficient functions of the DE/problem to create directly
    callable functions from Python.
    """

    default_lambda = "lambda x,y:"
    global_vars = None
    for name in functions:
        if functions[name] == '?':
            functions[name] = '0'
        elif functions[name] == "x" or functions[name] == "y":
            #If it is indicated that the provided U & V values to be used
            if not global_vars:
                x, y = [0] * nodes.__len__(), [0] * nodes.__len__()
                for i, node in enumerate(nodes):
                    x[i] = node[0]
                    y[i] = node[1]

                from scipy.interpolate import bisplrep, bisplev
                # Fit a bivariate B-spline to U and V values to calculate
                # values that are not on the nodes  This "global_vars"
                # dictionary is provided to eval for the lambda's to work
                # properly
                global_vars = {
                    "x_tck": bisplrep(x, y, UV_data[0]),
                    "y_tck": bisplrep(x, y, UV_data[1]),
                    "bisplev": bisplev
                }

            functions[name] = eval(
                "lambda x,y: bisplev(x, y, {0}_tck)".format(functions[name]),
                global_vars
            )
            continue

        functions[name] = default_lambda + functions[name]
        functions[name] = eval(functions[name])
    return functions

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 get_optical_path_map(self,size=(20, 20),  mask=None):
     """Return the optical path of the rays hitting the detector.
     
     This method uses the optical path of the rays hitting the surface, to 
     create a optical path map. The returned value is an interpolation of the
     values obtained by the rays.
     
     Warning: 
         If the rays hitting the surface are produced by more than one 
         optical source, the returned map migth not be valid.  
     
     *Atributes*
     
     *size*
         Tuple (nx,ny) containing the number of samples of the returned map.
         The map size will be the same as the CCD
     
     *mask*
         Shape instance containig the mask of the apperture. If not given, 
         the mask will be automatically calculated.
     
     *Return value*
     
     A masked array as defined in the numpy.ma module, containig the optical paths
     """
     
     X,Y,Z=self.get_optical_path_data()    
 
     rv=bisplrep(X,Y,Z)
     nx, ny=size
     xs, ys=self.size
     xi=-xs/2.
     xf=-xi
     yi=-ys/2.
     yf=-yi
     
     xd=linspace(xi, xf,nx)
     yd=linspace(yi, yf,ny)
     data=bisplev(xd,yd,rv)
     
     if mask!=None:
         assert(isinstance(mask, Shape))
         X, Y=meshgrid(xd, yd)
         m= ~mask.hit((X, Y, 0))
         retval= ma.array(data, mask=m)
     else:
         retval=data
     return retval

def contourpk(x,y,f, levels=None,xb=None,xe=None,yb=None,ye=None,s=60,kx=1,ky=1,dolabels=True, **kwargs):
   '''Make a contour plot with 1-D array inputs for X, Y and F. This is a  
   wrapper to convert lists of points (X,Y,Z) in 2-D arrays, then calls contour()
   
   Parameters
   ----------
   X, Y: ndarrays[:], 1D on a 2D plane
     coordinates X, Y
   Z : ndarray[:], 1D function on X,Y
   
   kwargs : dict
   -------------
    - **xb,xe,yb,ye** : float
      limits x,y for bispline interpolation valid region
    - **s** : float
      smoothing parameter for bisplrep
    - **kx, ky** : int
      order for the interpolation 
    - **dolabels** : bool
      labels on the contours if true
    - **levels** : list
      contour levels 
             
   Note
   ----
   warning: X, Y axis may have been interchanged
   ''' 
   import numpy
   from scipy import interpolate
   from pylab import contour, plt
   x1, x2, y1, y2 = min(x), max(x), min(y), max(y)
   xx = numpy.linspace(x1, x2)
   yy = numpy.linspace(y1, y2)
   X, Y = numpy.meshgrid(xx, yy)
   shp = X.shape
   task = 0
   tck = interpolate.bisplrep(x,y,f,kx=kx,ky=ky,s=s,xb=xb,xe=xe,yb=yb,ye=ye)     
   Z = interpolate.bisplev(xx, yy, tck)
   if levels == None:
      C = contour(Y, X, Z,**kwargs)
   else:
      C = contour(Y, X, Z, levels=levels,**kwargs)  
   if dolabels:     
      plt.clabel(C, inline=1,fontsize=10)    
   return Y,X,Z,tck, C

 def _data_from_ndvar(self, ndvar):
     v = ndvar.get_data(('sensor',))
     locs = ndvar.sensor.get_locs_2d(self._proj)
     if self._interpolation == 'spline':
         tck = interpolate.bisplrep(locs[:, 1], locs[:, 0], v, kx=5, ky=5)
         return interpolate.bisplev(self._grid, self._grid, tck)
     else:
         isnan = np.isnan(v)
         if np.any(isnan):
             nanmap = interpolate.griddata(locs, isnan, self._mgrid,
                                           method=self._interpolation)
             mask = nanmap > 0.5
             v = np.where(isnan, 0, v)
             vmap = interpolate.griddata(locs, v, self._mgrid,
                                         method=self._interpolation)
             np.place(vmap, mask, np.NaN)
             return vmap
         return interpolate.griddata(locs, v, self._mgrid,
                                     method=self._interpolation)

def plot_splined_contours(fname, label=None, colors=None, islog=False, linestyle='-', 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=0)
    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)
    
    my_levels = np.array([0.1, 2.3, 6.17, 11.8])
    
    a = plt.contour(Q_args, N_args, lnL_splined, my_levels, colors=colors, linestyles=linestyle)

    if label is not None:
        plt.plot(Q_args[int(abs(num_spline_points))/2, 0], N_args[0, int(abs(num_spline_points))/2], linestyle, color=colors, label=label)

    return a

	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

def contour_logL(xvals, yvals, logL, resample=False):
    """
    draw 1, 2, and 3-sigma contours from a gridded log-likelihood
    """
    #resample logL
    if resample:
        x,y = np.meshgrid(xvals, yvals)
        x.resize(x.size)
        y.resize(y.size)
        logL = logL.copy()
        logL.resize(logL.size)
        
        tck = interpolate.bisplrep(x, y, logL)
        
        xvals = np.linspace(xvals[0], xvals[-1], 20)
        yvals = np.linspace(yvals[0], yvals[-1], 20)
        
        logL = interpolate.bisplev(xvals, yvals, tck)
    
    #normalize logL: start by making the max logL=0
    L = np.exp(logL - logL.max())
    L /= L.sum()

    #assume a well-behaved peak.  Find 1-, 2-, and 3-sigma values
    Llist = L.copy().reshape(L.size)
    Llist.sort()
    levels = Llist.cumsum()

    i1 = levels.searchsorted(1 - 0.63)
    i2 = levels.searchsorted(1 - 0.95)
    i3 = levels.searchsorted(1 - 0.997)

    v1 = Llist[i1]
    v2 = Llist[i2]
    v3 = Llist[i3]

    pylab.contour(xvals, yvals, L, [v1,v2,v3])