def evaluate_spline(header, row, spline, limits):
        limits = (max(min(spline.get_knots()), float(limits[0])), min(max(spline.get_knots()), float(limits[1])))

        ys = np.linspace(limits[0], limits[1], len(header) * SplineModel.samples)
        ps = np.exp(spline(ys)) * (limits[1] - limits[0]) / (len(header) * SplineModel.samples)
        ps = ps / sum(ps)
        cfs = np.cumsum(ps)

        if 'mean' in header or 'var' in header or 'skew' in header:
            mean = sum(ps * ys)
        if 'var' in header or 'skew' in header:
            var = sum(ps * np.square(ys - mean))
        
        error = 0
        for ii in range(1, len(header)):
            if isinstance(header[ii], float):
                error = error + np.abs(SplineModelConditional.find_nearest(cfs, header[ii], ys) - float(row[ii]))
            elif header[ii] == 'mean':
                error = error + np.abs(mean - float(row[ii]))
            elif header[ii] == 'mode':
                mode = ys[ps.argmax()]
                error = error + np.abs(mode - float(row[ii]))
            elif header[ii] == 'var':
                error = error + np.sqrt(np.abs(var - float(row[ii])))
            elif header[ii] == 'skew':
                skew = sum(ps * np.pow((ys - mean) / sqrt(var), 3))
                error = error + np.pow(np.abs(skew - float(row[ii])), 1.0/3)

        return error

    def beta(self,m1,d,g=1.4,i=0):
        p=-(m1*m1+2.)/m1/m1-g*np.sin(d)*np.sin(d)
        q=(2.*m1*m1+1.)/ np.pow(m1,4.)+((g+1.)*(g+1.)/4.+
                                          (g-1.)/m1/m1)*np.sin(d)*np.sin(d)
        r=-np.cos(d)*np.cos(d)/np.pow(m1,4.)

        a=(3.*q-p*p)/3.
        b=(2.*p*p*p-9.*p*q+27.*r)/27.

        test=b*b/4.+a*a*a/27.

        if (test>0.0):
            return -1.0
        elif (test==0.0):
          x1=np.sqrt(-a/3.)
          x2=x1
          x3=2.*x1
          if(b>0.0):
            x1*=-1.
            x2*=-1.
            x3*=-1.

        if(test<0.0):
          phi=np.acos(np.sqrt(-27.*b*b/4./a/a/a))
          x1=2.*np.sqrt(-a/3.)*np.cos(phi/3.)
          x2=2.*np.sqrt(-a/3.)*np.cos(phi/3.+np.pi*2./3.)
          x3=2.*np.sqrt(-a/3.)*np.cos(phi/3.+np.pi*4./3.)
          if(b>0.0):
            x1*=-1.
            x2*=-1.
            x3*=-1.

        s1=x1-p/3.
        s2=x2-p/3.
        s3=x3-p/3.

        if(s1<s2 and s1<s3):
          t1=s2
          t2=s3
        elif(s2<s1 and s2<s3):
          t1=s1
          t2=s3
        else:
          t1=s1
          t2=s2

        b1=np.asin(np.sqrt(t1))
        b2=np.asin(np.sqrt(t2))

        betas=b1
        betaw=b2
        if(b2>b1):
          betas=b2
          betaw=b1

        if(i==0):
            return betaw
        if(i==1):
            return betas

 def lowpass_prototype_order(self, wpass, wstop):
     if self.btype == "butter":
         order = (np.log(np.pow(10,0.1*self.gstop)-1.0)-(np.pow(0.1*self.gpass)-1.0))/(np.log(abs(wstop/w0 - w0/wstop)) - np.log(abs(wpass/w0 - w0/wpass))) #page 68
     elif self.ftype in ('cheby1', 'cheby2'):
         """TBI"""
     else:
         raise Exception("Unknown filter's approximation type.")
         sys.exit(1)

def gamma_trans(mat, gamma):
    gamma_mean = np.pow(mean_rgb / 255, gamma)
    tmp_mat = np.pow(mat / 255, gamma)
    gamma_mat = np.zeros(mat.shape, dtype=np.float)
    gamma_mat[:, :, 0] = tmp_mat[:, :, 2] - gamma_mean[:, :, 2]
    gamma_mat[:, :, 1] = tmp_mat[:, :, 1] - gamma_mean[:, :, 1]
    gamma_mat[:, :, 2] = tmp_mat[:, :, 0] - gamma_mean[:, :, 0]
    return gamma_mat

    def fn(self, current_neuron, best_neuron):
        """
        Calculate the value for the multi RBF function.
        :param current_neuron: The current neuron.
        :param best_neuron: The best neuron.
        :return: A percent that determines the amount of training the current
        neuron should get.  Usually 100% when it is the bestNeuron.
        """
        vector = np.zeros(len(self.displacement))
        vector_current = self.translate_coordinates(current_neuron)
        vector_best = self.translate_coordinates(best_neuron)
        for i in range(len(vector_current)):
            vector[i] = vector_current[i] - vector_best[i]

        if self.hexagon:

            if len(self.size) !=2:
                raise Exception("Hexagon lattice can only be used in two dimensions.")

            row = vector[1]
            col = vector[0]
            even_indent = 1
            odd_indent = 2.5
            indent = odd_indent if row%2==1 else even_indent

            vector[1] = int(NeighborhoodRBF.SQ75+(row * NeighborhoodRBF.SQ75))
            vector[0] = int(indent+(3*col))

        if self.type == NeighborhoodRBF.TYPE_GAUSSIAN:
            value = 0

            for i in range(len(self.size)):
                value += np.power(vector[i], 2) / (2.0 * self.width * self.width)
            return np.exp(-value)

        elif self.type == NeighborhoodRBF.TYPE_MULTIQUADRIC:
            value = 0

            for i in range(len(self.size)):
                value += np.pow(vector[i], 2) + (self.width * self.width)
            return np.sqrt(value)
        elif self.type == NeighborhoodRBF.TYPE_INVERSE_MULTIQUADRIC:
            value = 0
            for i in range(len(self.size)):
                value += np.pow(vector[i], 2) + (self.width * self.width)
            return 1 / np.sqrt(value)
        elif self.type == NeighborhoodRBF.TYPE_MEXICAN_HAT:
            # calculate the "norm", but don't take square root
            # don't square because we are just going to square it
            norm = 0
            for i in range(len(self.size)):
                norm += np.pow(vector[i], 2)

            # calculate the value
            return (1 - norm) * np.exp(-norm / 2)
        else:
            raise Exception("Invalid RBF function type: {}".format(self.type))

def ll(e_i,e_m,delt=0.1):
    """

    :param e_i: inclusion dielectric
    :param e_m: bulk dielectric
    :param delt: fraction of inclusion
    :return:
    """
    from numpy import pow
    return pow(delt*pow(e_i,1/3.)+(1-delt)*pow(e_m,1/3.),3.)

 def __call__( self, eps, a, b, c, d, ee, f, g, h, i, k ):
     '''
     Implements the response function with arrays as variables.
     first extract the variable discretizations from the orthogonal grid.
     '''
     return ( a + pow( b, 2 ) + pow( c, 3 ) + pow( d, 4 ) + pow( ee, 5 ) + pow( f, 6 ) \
             + pow( g, 7 ) + pow( h, 8 ) + pow( i, 9 ) + pow( k, 10 ) ) * eps;

def distance_geodetic(lat1,lon1,lat2,lon2):
	RADIUS = 6371 #KM

	D2R =  pi/180;

	lat1_rad = float(lat1)*D2R;
	lat2_rad = float(lat2)*D2R;
	lon1_rad = float(lon1)*D2R;
	lon2_rad = float(lon2)*D2R;

	a =    pow( sin((lat1_rad-lat2_rad)/2),2) +  cos(lat1_rad)* cos(lat2_rad)*  pow( sin((lon1_rad-lon2_rad)/2),2);

	distance =  abs(2*RADIUS* arctan2( sqrt(a),  sqrt(1-a)));
	return distance*1000; #meter

def geodetic_to_ECEF(lat,lon):
	D2R =  pi/180;
	#meters
	a = 6378137.0 
	#first eccentricity squared 
	pow_e_2 = 6.69437999014*  pow(10,-3)
	#assume 0	
	h = 0 
	lat_rad = float(lat)*D2R;
	lon_rad = float(lon)*D2R;
	N = a/ sqrt(1-pow_e_2*  pow( sin(lat_rad),2))
	x = (N+h)* cos(lat_rad)* cos(lon_rad)
	y = (N+h)* cos(lat_rad)* sin(lon_rad)
	z = (N*(1-pow_e_2)+h)* sin(lat_rad)
	return [x,y,z]

    def invTransform(self, value):
        """
        Inverse transformation function

        :param float value: Value
        :return: Modified value
        
        .. seealso::
        
            :py:meth:`transform()`
        """
        if value < 0.:
            return -np.pow(-value, self.__exponent)
        else:
            return np.pow(value, self.__exponent)

    def backward(self, bottom_data, bottom_diff, top_data, top_diff):
        padded_ratio = Array.zeros(1, bottom_data.shape[1] + self.size - 1,
                                   bottom_data.shape[2], bottom_data.shape[3])
        accum_ratio = Array.zeros(1, 1, bottom_data.shape[2],
                                  bottom_data.shape[3])
        accum_ratio_times_bottom = Array.zeros(1, 1, bottom_data.shape[2],
                                               bottom_data.shape[3])
        cache_ratio_value = 2.0 * self.apha * self.beta / self.size
        bottom_diff = np.pow(self.scale, -self.beta)
        bottom_diff *= top_diff

        inverse_pre_pad = self.size - (self.size + 1) / 2
        for n in range(bottom_data.shape[0]):
            padded_ratio[0, inverse_pre_pad] = top_diff[n] * top_data[n]
            padded_ratio[0, inverse_pre_pad] /= self.scale[n]
            accum_ratio.fill(0)
            for c in range(self.size - 1):
                accum_ratio += padded_ratio[0, c]

            for c in range(bottom_data.shape[1]):
                accum_ratio += padded_ratio[0, c + self.size - 1]
                accum_ratio_times_bottom += bottom_data[n, c] * accum_ratio
                bottom_data[n, c] += -cache_ratio_value * \
                    accum_ratio_times_bottom
                accum_ratio += -1 * padded_ratio[0, c]

    def process_chunk(self, data):

        moment_data = numpy.log(data)
        
        moments = numpy.zeros(self.mmax - self.mmin, dtype=numpy.float32)

        mean = numpy.nanmean(moment_data)

        moment_data = moment_data - mean

        if self.mmin == 1:
            temp = numpy.ones(len(moment_data), dtype=numpy.float32)
        elif self.mmin == 2:
            temp = moment_data
        else:
            temp = numpy.pow(moment_data, self.mmin-1)

        for i in range(0, self.mmax-self.mmin):
            temp = temp * moment_data
            moments[i] = numpy.nanmean(temp)

        if self.mmin == 1:
            moments[0] = mean

        return moments

    def schecterFunctionL(self,L,phistar,Lstar,alpha):
        """ Schecter function for galaxy LF
        Default to redshift ??? LBG properties
        """
        phi=phistar*numpy.pow(L/Lstar,alpha)*numpy.exp(-L/Lstar)

        return phi

    def run_till_convergence(self, tol=.05):
        scale = np.mean(np.sum(np.pow(data, 2), 1), 0)
        runs = 0
        while 1:
            runs += 1
            means_0 = array(self.means)
            self._iterate()
            divergence = np.mean(np.sum(np.pow(means_0 - means, 2), 1), 0)
            if divergence / scale < tol:
                failed = False
                break
            if runs > self.max_runs:
                failed = True
                break

        return failed

    def func(self, t, y):
        m = self.m
        k = self.k
        p = self.p

        x, v = y
        return v, (self.fext(t,x) - k*np.pow(x,p-1))/m

def squared_error(y, est):
	num_classes = y.shape[1]
	rval = 0
	for c in range(num_classes):
		est_c = est[np.where(y[:,c] == 1)[0]]
		num_class_c = np.sum(y[:,c])
		rval += np.sum(np.pow(np.subtract(np.ones([num_class_c, 1]), est_c), 2))
	return rval

def generate_shear_layer_profile(conv_u,eta,u,u_y,u_yy):
	N = len(u)-1
	for i in range(0,N+1):
		u[i] = 1.0 + conv_u * np.tanh(eta[i])
		u_y[i] = conv_u * ( 1.0 - np.pow(np.tanh(eta[i]), 2.0) )
		u_yy[i] = 0.0
	print "############# Need to calculate analytical function for u_yy ###############"
	return

def generate_wake_profile(wake_half_width,u_deficit,eta,u,u_y,u_yy):
	N = len(u)-1
	for i in range(0,N+1):
		sech_n = 1.0 / np.cosh(eta[ii] / wake_half_width)
		u[i] = 1.0 - u_deficit * pow(sech_n,2.0)
		u_y[ii] = 2.0 * u_deficit * np.pow(sech_n,2.0) * np.tanh(eta[i])
		u_yy[ii] = 0.0
	print "############# Need to calculate analytical function for u_yy ###############"
	return 

    def inverse(self, value):
        if not self.scaled():
            raise ValueError("Not invertible until scaled")
        vmin, vmax = float(self.vmin), float(self.vmax)

        if cbook.iterable(value):
            val = np.ma.asarray(value)
            return vmin * np.ma.power((vmax/vmin), val)
        else:
            return vmin * np.pow((vmax/vmin), value)

def distance(node1,node2,lat_a, lon_a):
	lat1 = float(lat_a[node1])
	lat2 = float(lat_a[node2])
	lon1 = float(lon_a[node1])
	lon2 = float(lon_a[node2])

	RADIUS = 6371 #KM

	D2R =  pi/180;

	lat1_rad = lat1*D2R;
	lat2_rad = lat2*D2R;
	lon1_rad = lon1*D2R;
	lon2_rad = lon2*D2R;

	a =    pow( sin((lat1_rad-lat2_rad)/2),2) +  cos(lat1_rad)* cos(lat2_rad)*   pow( sin((lon1_rad-lon2_rad)/2),2);

	distance =  abs(2*RADIUS* arctan2( sqrt(a),  sqrt(1-a)));
	return distance*1000; #meter