Calculate the second derivatives of the interpolating function at the
tabulated points. At the endpoints, we use a natural spline
with zero second derivative on that boundary.
545540
private void sety2(double[] x, double[] y) {
        double p, qn, sig, un;
        double[] u = new double[n - 1];

        y2[0] = u[0] = 0.0;

        for (int i = 1; i < n - 1; i++) {
            sig = (x[i] - x[i - 1]) / (x[i + 1] - x[i - 1]);
            p = sig * y2[i - 1] + 2.0;
            y2[i] = (sig - 1.0) / p;
            u[i] = (y[i + 1] - y[i]) / (x[i + 1] - x[i]) - (y[i] - y[i - 1]) / (x[i] - x[i - 1]);
            u[i] = (6.0 * u[i] / (x[i + 1] - x[i - 1]) - sig * u[i - 1]) / p;
        }

        qn = un = 0.0;
        y2[n - 1] = (un - qn * u[n - 2]) / (qn * y2[n - 2] + 1.0);
        for (int k = n - 2; k >= 0; k--) {
            y2[k] = y2[k] * y2[k + 1] + u[k];
        }
    }
1158086
protected double interpolate( double boundA , double boundB )
	{
		double alphaNew;

		// interpolate minimum for rapid convergence
		if( Double.isNaN(gp) ) {
			alphaNew = SearchInterpolate.quadratic(fprev, gprev, stprev, fp, stp);
		} else {
			alphaNew = SearchInterpolate.cubic2(fprev, gprev, stprev, fp, gp, stp);
			if( Double.isNaN(alphaNew))
				alphaNew = SearchInterpolate.quadratic(fprev, gprev, stprev, fp, stp);
		}

		// order the bound
		double l,u;
		if( boundA < boundB ) {
			l=boundA;u=boundB;
		} else {
			l=boundB;u=boundA;
		}

		// enforce min/max allowed values
		if( alphaNew < l )
			alphaNew = l;
		else if( alphaNew > u )
			alphaNew = u;

		return alphaNew;
	}
234649
def second_derivative_5(var, key):
    ''''''
    global derivative_data
    import mavutil
    tnow = mavutil.mavfile_global.timestamp

    if not key in derivative_data:
        derivative_data[key] = (tnow, [var]*5)
        return 0
    (last_time, data) = derivative_data[key]
    data.pop(0)
    data.append(var)
    derivative_data[key] = (tnow, data)
    h = (tnow - last_time)
    # N=5 2nd derivative from
    # http://www.holoborodko.com/pavel/numerical-methods/numerical-derivative/smooth-low-noise-differentiators/
    ret = ((data[4] + data[0]) - 2*data[2]) / (4*h**2)
    return ret
76859
def interpolate_series(start, end, steps, how='linear'):
    
    return INTERP_DICT[how](start, end, steps)
1109270
def derivatives(self, x, y, theta_E, r_trunc, center_x=0, center_y=0):
        
        x_shift = x - center_x
        y_shift = y - center_y

        dphi_dr = self._dphi_dr(x_shift, y_shift, theta_E, r_trunc)
        dr_dx, dr_dy = self._dr_dx(x_shift, y_shift)
        f_x = dphi_dr * dr_dx
        f_y = dphi_dr * dr_dy
        return f_x, f_y
44087
def derivative(self, t, n=1):
        """"""
        angle = radians(self.theta + t*self.delta)
        phi = radians(self.rotation)
        rx = self.radius.real
        ry = self.radius.imag
        k = (self.delta*2*pi/360)**n  # ((d/dt)angle)**n

        if n % 4 == 0 and n > 0:
            return rx*cos(phi)*cos(angle) - ry*sin(phi)*sin(angle) + 1j*(
                rx*sin(phi)*cos(angle) + ry*cos(phi)*sin(angle))
        elif n % 4 == 1:
            return k*(-rx*cos(phi)*sin(angle) - ry*sin(phi)*cos(angle) + 1j*(
                -rx*sin(phi)*sin(angle) + ry*cos(phi)*cos(angle)))
        elif n % 4 == 2:
            return k*(-rx*cos(phi)*cos(angle) + ry*sin(phi)*sin(angle) + 1j*(
                -rx*sin(phi)*cos(angle) - ry*cos(phi)*sin(angle)))
        elif n % 4 == 3:
            return k*(rx*cos(phi)*sin(angle) + ry*sin(phi)*cos(angle) + 1j*(
                rx*sin(phi)*sin(angle) - ry*cos(phi)*cos(angle)))
        else:
            raise ValueError("n should be a positive integer.")
975887
def interpolate(x, y, z, interp_type='linear', hres=50000,
                minimum_neighbors=3, gamma=0.25, kappa_star=5.052,
                search_radius=None, rbf_func='linear', rbf_smooth=0,
                boundary_coords=None):
    """"""
    return interpolate_to_grid(x, y, z, interp_type=interp_type, hres=hres,
                               minimum_neighbors=minimum_neighbors, gamma=gamma,
                               kappa_star=kappa_star, search_radius=search_radius,
                               rbf_func=rbf_func, rbf_smooth=rbf_smooth,
                               boundary_coords=boundary_coords)
1036613
def derivatives(self, x, y, n_sersic, R_sersic, k_eff, center_x=0, center_y=0):
        
        x_ = x - center_x
        y_ = y - center_y
        r = np.sqrt(x_**2 + y_**2)
        if isinstance(r, int) or isinstance(r, float):
            r = max(self._s, r)
        else:
            r[r < self._s] = self._s
        alpha = -self.alpha_abs(x, y, n_sersic, R_sersic, k_eff, center_x, center_y)
        f_x = alpha * x_ / r
        f_y = alpha * y_ / r
        return f_x, f_y
509517
private static int interpolate(
        double y_minus, // = f(-1)
        double y_0, // = f(0)
        double y_plus, // = f(+1)
        double[] result // xe, ye, root1, root2
    ) {

        double a = 0.5 * (y_plus + y_minus) - y_0;
        double b = 0.5 * (y_plus - y_minus);

        double xe = -b / (2.0 * a);
        double ye = (a * xe + b) * xe + y_0;
        double dis = b * b - 4 * a * y_0;

        double root1 = Double.NaN;
        double root2 = Double.NaN;
        int count = 0;

        if (dis >= 0) {
            double dx = 0.5 * Math.sqrt(dis) / Math.abs(a);
            if (Math.abs(xe - dx) <= 1.0) {
                root1 = xe - dx;
                count++;
            }
            if (Math.abs(xe + dx) <= 1.0) {
                root2 = xe + dx;
                count++;
            }
        }

        result[0] = xe;
        result[1] = ye;
        result[2] = root1;
        result[3] = root2;

        return count;

    }
124469
def evaluate(self,       x, y, x_derivative=0, smooth=0, simple='auto'):
        

        if simple=='auto': simple = self.simple
        
        # find which values y is in between
        for n in range(0, len(self.y_values)-1):
            # if it's in between, interpolate!
            if self.y_values[n] <= y and self.y_values[n+1] >= y:
                y1 = self.y_values[n]
                y2 = self.y_values[n+1]
                z1 = self.x_splines[y1].evaluate(x, x_derivative, smooth, simple)
                z2 = self.x_splines[y2].evaluate(x, x_derivative, smooth, simple)
                return z1 + (y-y1)*(z2-z1)/(y2-y1)

        print("YARG! The y value "+str(y)+" is out of interpolation range!")            
        if y >= self.y_values[-1]: return self.x_splines[self.y_values[-1]].evaluate(x, x_derivative, smooth, simple)
        else                     : return self.x_splines[self.y_values[0]].evaluate(x, x_derivative, smooth, simple)
399652
def derivative(func, x0, dx=1.0, n=1, args=(), order=3):
    '''
    '''
    if order < n + 1:
        raise ValueError
    if order % 2 == 0:
        raise ValueError
    weights = central_diff_weights(order, n)
    tot = 0.0
    ho = order >> 1
    for k in range(order):
        tot += weights[k]*func(x0 + (k - ho)*dx, *args)
    return tot/product([dx]*n)
441477
def smooth(polylines):
    
    for c in polylines:
        if len(c) < 9:
            # smoothing wouldn't make sense here
            continue
        x = c[:, 0]
        y = c[:, 1]

        t = np.arange(x.shape[0], dtype=float)
        t /= t[-1]
        x = UnivariateSpline(t, x)(t)
        y = UnivariateSpline(t, y)(t)
        c[:, 0] = x
        c[:, 1] = y
366856
def _deriv_logaddexp2(x1, x2):
    """"""
    y1 = np.exp2(x1)
    y2 = np.exp2(x2)
    df_dx1 = y1 / (y1 + y2)
    df_dx2 = y2 / (y1 + y2)
    return np.vstack([df_dx1, df_dx2]).T
1089007
def _derY(self,x,y):
        '''
        
        '''
        x_pos, y_pos = self.findSector(x,y)
        alpha, beta = self.findCoords(x,y,x_pos,y_pos)

        # Get four corners data for each point
        xA = self.x_values[x_pos,y_pos]
        xB = self.x_values[x_pos+1,y_pos]
        xC = self.x_values[x_pos,y_pos+1]
        xD = self.x_values[x_pos+1,y_pos+1]
        yA = self.y_values[x_pos,y_pos]
        yB = self.y_values[x_pos+1,y_pos]
        yC = self.y_values[x_pos,y_pos+1]
        yD = self.y_values[x_pos+1,y_pos+1]
        fA = self.f_values[x_pos,y_pos]
        fB = self.f_values[x_pos+1,y_pos]
        fC = self.f_values[x_pos,y_pos+1]
        fD = self.f_values[x_pos+1,y_pos+1]

        # Calculate components of the alpha,beta --> x,y delta translation matrix
        alpha_x = (1-beta)*(xB-xA) + beta*(xD-xC)
        alpha_y = (1-beta)*(yB-yA) + beta*(yD-yC)
        beta_x  = (1-alpha)*(xC-xA) + alpha*(xD-xB)
        beta_y  = (1-alpha)*(yC-yA) + alpha*(yD-yB)

        # Invert the delta translation matrix into x,y --> alpha,beta
        det = alpha_x*beta_y - beta_x*alpha_y
        y_alpha = -beta_x/det
        y_beta  = alpha_x/det

        # Calculate the derivative of f w.r.t. alpha and beta
        dfda = (1-beta)*(fB-fA) + beta*(fD-fC)
        dfdb = (1-alpha)*(fC-fA) + alpha*(fD-fB)

        # Calculate the derivative with respect to x (and return it)
        dfdy = y_alpha*dfda + y_beta*dfdb
        return dfdy
109539
def endpoints_minima(self, slope_cutoff=5e-3):
        
        energies = self.energies
        try:
            sp = self.spline()
        except:
            print("Energy spline failed.")
            return None
        der = sp.derivative()
        der_energies = der(range(len(energies)))
        return {"polar": abs(der_energies[-1]) <= slope_cutoff,
                "nonpolar": abs(der_energies[0]) <= slope_cutoff}
318111
private static double cubicInterpolate(final double x0, final double y0, final double x1,
      final double y1, final double x2, final double y2, final double x3, final double y3,
      final double x) {
    final double l0Numer = (x - x1) * (x - x2) * (x - x3);
    final double l1Numer = (x - x0) * (x - x2) * (x - x3);
    final double l2Numer = (x - x0) * (x - x1) * (x - x3);
    final double l3Numer = (x - x0) * (x - x1) * (x - x2);

    final double l0Denom = (x0 - x1) * (x0 - x2) * (x0 - x3);
    final double l1Denom = (x1 - x0) * (x1 - x2) * (x1 - x3);
    final double l2Denom = (x2 - x0) * (x2 - x1) * (x2 - x3);
    final double l3Denom = (x3 - x0) * (x3 - x1) * (x3 - x2);

    final double term0 = (y0 * l0Numer) / l0Denom;
    final double term1 = (y1 * l1Numer) / l1Denom;
    final double term2 = (y2 * l2Numer) / l2Denom;
    final double term3 = (y3 * l3Numer) / l3Denom;

    return term0 + term1 + term2 + term3;
  }
736993
def derivatives(self, x, y, center_x = 0, center_y = 0, **kwargs):

        

        assert 'norm' in kwargs.keys(), "key word arguments must contain 'norm', " \
                                        "the normalization of deflection angle in units of arcsec."

        x_ = x - center_x
        y_ = y - center_y
        R = np.sqrt(x_**2 + y_**2)

        alpha = self._interp(x_, y_, **kwargs)

        cos_theta = x_ * R ** -1
        sin_theta = y_ * R ** -1

        f_x, f_y = alpha * cos_theta, alpha * sin_theta

        return f_x, f_y