Returns the first derivative of the function with respect to Z at each
        value in (x,y,z).  Only called internally by HARKinterpolator3D._derZ.
326598
def _derZ(self,x,y,z):
        '''
        
        '''
        m = len(x)
        temp = np.zeros((m,self.funcCount))
        for j in range(self.funcCount):
            temp[:,j] = self.functions[j](x,y,z)
        temp[np.isnan(temp)] = np.inf
        i = np.argmin(temp,axis=1)
        y = temp[np.arange(m),i]
        dfdz = np.zeros_like(x)
        for j in range(self.funcCount):
            c = i == j
            dfdz[c] = self.functions[j].derivativeZ(x[c],y[c],z[c])
        return dfdz
378482
public Vector3d hermite(Vector3dc t0, Vector3dc v1, Vector3dc t1, double t, Vector3d dest) {
        double t2 = t * t;
        double t3 = t2 * t;
        dest.x = (x + x - v1.x() - v1.x() + t1.x() + t0.x()) * t3 + (3.0 * v1.x() - 3.0 * x - t0.x() - t0.x() - t1.x()) * t2 + x * t + x;
        dest.y = (y + y - v1.y() - v1.y() + t1.y() + t0.y()) * t3 + (3.0 * v1.y() - 3.0 * y - t0.y() - t0.y() - t1.y()) * t2 + y * t + y;
        dest.z = (z + z - v1.z() - v1.z() + t1.z() + t0.z()) * t3 + (3.0 * v1.z() - 3.0 * z - t0.z() - t0.z() - t1.z()) * t2 + z * t + z;
        return dest;
    }
863604
def derivative(self, x, der=1):
        
        from scipy.interpolate import splev
        return splev(x, self._sp, der=der)
637458
def derivative(self, z, x, y, fase):
        """"""
        return deriv_H(self, z, x, y, fase)
897043
def _2ndderiv_xyz(self,x,y,z,i,j):
        """"""
        return -np.pi*self._rhoc_M*self.a**3*self._b*self._c *\
            _2ndDerivInt(x,y,z,self._a2,self._b2*self._a2,self._c2*self._a2,self.n,i,j)
1165289
def differentiate(self, n=1):
        
        ak = self.coefficients()
        a_, b_ = self.domain()
        for _ in range(n):
            ak = self.differentiator(ak)
        return self.from_coeff((2./(b_-a_))**n*ak, domain=self.domain())
922909
def zmeshgrid(d):
    
    lin = xfun(2, d)
    one = ones(2, d)

    xx = zkronv(lin, one)
    yy = zkronv(one, lin)

    return xx, yy
726891
def hessian(self, x, y, center_x = 0, center_y = 0, **kwargs):
        

        diff = 1e-6
        alpha_ra, alpha_dec = self.derivatives(x, y, center_x = center_x, center_y = center_y,
                                               **kwargs)

        alpha_ra_dx, alpha_dec_dx = self.derivatives(x + diff, y, center_x = center_x, center_y = center_y,
                                               **kwargs)
        alpha_ra_dy, alpha_dec_dy = self.derivatives(x, y + diff, center_x = center_x, center_y = center_y,
                                               **kwargs)

        dalpha_rara = (alpha_ra_dx - alpha_ra) / diff
        dalpha_radec = (alpha_ra_dy - alpha_ra) / diff
        dalpha_decdec = (alpha_dec_dy - alpha_dec) / diff

        f_xx = dalpha_rara
        f_yy = dalpha_decdec
        f_xy = dalpha_radec

        return f_xx, f_yy, f_xy
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
1369553
def forward(self, x):
        
        fwd = np.zeros((self.num_matrices, self.N, self.N), settings.float_type)
        indices = np.tril_indices(self.N, 0)
        z = np.zeros(len(indices[0])).astype(int)
        for i in range(self.num_matrices):
            fwd[(z + i,) + indices] = x[i, :]
        return fwd.squeeze(axis=0) if self.squeeze else fwd
975208
def differentiator(self, A):
        
        m = len(A)
        SA = (A.T* 2*np.arange(m)).T
        DA = np.zeros_like(A)
        if m == 1: # constant
            return np.zeros_like(A[0:1])
        if m == 2: # linear
            return A[1:2,]
        DA[m-3:m-1,] = SA[m-2:m,]
        for j in range(m//2 - 1):
            k = m-3-2*j
            DA[k] = SA[k+1] + DA[k+2]
            DA[k-1] = SA[k] + DA[k+1]
        DA[0] = (SA[1] + DA[2])*0.5
        return DA
1242338
def _2ndderiv_xyz(self,x,y,z,i,j):
        """"""
        return 4.*numpy.pi*self._b*self._c\
            *_2ndDerivInt(x,y,z,
                          lambda m: self._mdens(m),
                          lambda m: self._mdens_deriv(m),
                          self._b2,self._c2,i,j,glx=self._glx,glw=self._glw)
292810
def F_(self, X):
        
        if self._interpol:
            if not hasattr(self, '_F_interp'):

                if self._lookup:
                    x = self._x_lookup
                    F_x = self._f_lookup
                else:
                    x = np.linspace(0, self._max_interp_X, self._num_interp_X)
                    F_x = self._F(x)
                self._F_interp = interp.interp1d(x, F_x, kind='linear', axis=-1, copy=False, bounds_error=False,
                                                 fill_value=0, assume_sorted=True)
            return self._F_interp(X)
        else:
            return self._F(X)
938812
def f_inv(self, z, max_iterations=250, y=None):
        

        z = z.copy()
        y = np.ones_like(z)

        it = 0
        update = np.inf
        while np.abs(update).sum() > 1e-10 and it < max_iterations:
            fy = self.f(y)
            fgrady = self.fgrad_y(y)
            update = (fy - z) / fgrady
            y -= self.rate * update
            it += 1
        #if it == max_iterations:
        #    print("WARNING!!! Maximum number of iterations reached in f_inv ")
        #    print("Sum of roots: %.4f" % np.sum(fy - z))
        return y
975779
def _der(self,*args):
        '''
        
        '''
        if len(args) > 0:
            if _isscalar(args[0]):
                return 0.0
            else:
                shape = args[0].shape
                return np.zeros(shape)
        else:
            return 0.0
853706
def derivative(self,*args):
        '''
        
        '''
        if self.i_dim == 0:
            return np.ones_like(*args[0])
        else:
            return np.zeros_like(*args[0])
85829
def function(self, x, y, amp, R_sersic, n_sersic, center_x=0, center_y=0):
        
        #if n_sersic < 0.2:
        #    n_sersic = 0.2
        #if R_sersic < 10.**(-6):
        #    R_sersic = 10.**(-6)
        R_sersic = np.maximum(0, R_sersic)
        x_shift = x - center_x
        y_shift = y - center_y
        R = np.sqrt(x_shift*x_shift + y_shift*y_shift)
        if isinstance(R, int) or isinstance(R, float):
            R = max(self._smoothing, R)
        else:
            R[R < self._smoothing] = self._smoothing
        _, bn = self.k_bn(n_sersic, R_sersic)
        R_frac = R/R_sersic
        #R_frac = R_frac.astype(np.float32)
        if isinstance(R, int) or isinstance(R, float):
            if R_frac > 100:
                result = 0
            else:
                exponent = -bn*(R_frac**(1./n_sersic)-1.)
                result = amp * np.exp(exponent)
        else:
            R_frac_real = R_frac[R_frac <= 100]
            exponent = -bn*(R_frac_real**(1./n_sersic)-1.)
            result = np.zeros_like(R)
            result[R_frac <= 100] = amp * np.exp(exponent)
        return np.nan_to_num(result)