Differentiate a set of Chebyshev polynomial expansion
           coefficients
           Originally from http://www.scientificpython.net/pyblog/chebyshev-differentiation
            + (lots of) bug fixing + pythonisation
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
521587
def derivatives(self, x, y, grid_interp_x=None, grid_interp_y=None, f_=None, f_x=None, f_y=None, f_xx=None, f_yy=None, f_xy=None):
        
        #self._check_interp(grid_interp_x, grid_interp_y, f_, f_x, f_y, f_xx, f_yy, f_xy)
        n = len(np.atleast_1d(x))
        if n <= 1 and np.shape(x) == ():
        #if type(x) == float or type(x) == int or type(x) == type(np.float64(1)) or len(x) <= 1:
            f_x_out = self.f_x_interp(x, y, grid_interp_x, grid_interp_y, f_x)
            f_y_out = self.f_y_interp(x, y, grid_interp_x, grid_interp_y, f_y)
            return f_x_out[0][0], f_y_out[0][0]
        else:
            if self._grid and n >= self._min_grid_number:
                x_, y_ = util.get_axes(x, y)
                f_x_out = self.f_x_interp(x_, y_, grid_interp_x, grid_interp_y, f_x)
                f_y_out = self.f_y_interp(x_, y_, grid_interp_x, grid_interp_y, f_y)
                f_x_out = util.image2array(f_x_out)
                f_y_out = util.image2array(f_y_out)
            else:
                #n = len(x)
                f_x_out, f_y_out = np.zeros(n), np.zeros(n)
                for i in range(n):
                    f_x_out[i] = self.f_x_interp(x[i], y[i], grid_interp_x, grid_interp_y, f_x)
                    f_y_out[i] = self.f_y_interp(x[i], y[i], grid_interp_x, grid_interp_y, f_y)
        return f_x_out, f_y_out
963117
def integrate(self):
        
        coeffs = self.coefficients()
        a,b = self.domain()
        int_coeffs = 0.5*(b-a)*poly.chebyshev.chebint(coeffs)
        antiderivative = self.from_coeff(int_coeffs, domain=self.domain())
        return antiderivative - antiderivative(a)
81410
function(x, y) {

    var r = {};

    var i = degree(x);
    var j = degree(y);
    var trace = [];

    while (i >= j) {

      var tmp = r[i - j] = FIELD['div'](x[i] || 0, y[j] || 0);

      for (var k in y) {
        x[+k + i - j] = FIELD['sub'](x[+k + i - j] || 0, FIELD['mul'](y[k] || 0, tmp));
      }

      if (Polynomial['trace'] !== null) {

        var tr = {};
        for (var k in y) {
          tr[+k + i - j] = FIELD['mul'](y[k] || 0, tmp);
        }
        trace.push(new Polynomial(tr));
      }

      i = degree(x);
    }

    // Add rest
    if (Polynomial['trace'] !== null) {
      trace.push(new Polynomial(x));
      Polynomial['trace'] = trace;
    }
    return r;
  }
34442
function multiplyLikeTerms(node, polynomialOnly=false) {
  if (!Node.Type.isOperator(node)) {
    return Node.Status.noChange(node);
  }
  let status;

  if (!polynomialOnly) {
    status = arithmeticSearch(node);
    if (status.hasChanged()) {
      status.changeType = ChangeTypes.MULTIPLY_COEFFICIENTS;
      return status;
    }

    status = multiplyFractionsSearch(node);
    if (status.hasChanged()) {
      status.changeType = ChangeTypes.MULTIPLY_COEFFICIENTS;
      return status;
    }
  }

  status = multiplyPolynomialTerms(node);
  if (status.hasChanged()) {
    status.changeType = ChangeTypes.MULTIPLY_COEFFICIENTS;
    return status;
  }

  return Node.Status.noChange(node);
}
152239
def get_support_variables(polynomial):
    
    support = []
    if is_number_type(polynomial):
        return support
    for monomial in polynomial.expand().as_coefficients_dict():
        mon, _ = __separate_scalar_factor(monomial)
        symbolic_support = flatten(split_commutative_parts(mon))
        for s in symbolic_support:
            if isinstance(s, Pow):
                base = s.base
                if is_adjoint(base):
                    base = base.adjoint()
                support.append(base)
            elif is_adjoint(s):
                support.append(s.adjoint())
            elif isinstance(s, Operator):
                support.append(s)
    return support
414411
def _partial_diff(var, *params):
    
    if isinstance(var, sympy.Derivative):
        return sympy.Derivative(var.expr, *(var.variables + params))
    else:
        return D(var, *params)
1282771
def xyz2lab(x__, y__, z__):
    

    def f__(arr):
        """f__ function
        """

        return np.where(arr > 216.0 / 24389.0,
                        arr ** (1.0/3.0),
                        (1.0 / 3.0) * (29.0 / 6.0) * (29.0 / 6.0) * arr
                        + 4.0 / 29.0)
    fy_ = f__(y__ / 100.0)
    l__ = 116 * fy_ - 16
    a__ = 500.0 * (f__(x__ / 95.047) - fy_)
    b__ = 200.0 * (fy_ - f__(z__ / 108.883))
    return l__, a__, b__
337604
def derivatives(self, x, y, Rs, theta_Rs, e1, e2, center_x=0, center_y=0):
        
        phi_G, q = param_util.ellipticity2phi_q(e1, e2)
        x_shift = x - center_x
        y_shift = y - center_y
        cos_phi = np.cos(phi_G)
        sin_phi = np.sin(phi_G)
        e = min(abs(1. - q), 0.99)
        xt1 = (cos_phi*x_shift+sin_phi*y_shift)*np.sqrt(1 - e)
        xt2 = (-sin_phi*x_shift+cos_phi*y_shift)*np.sqrt(1 + e)
        R_ = np.sqrt(xt1**2 + xt2**2)
        rho0_input = self.nfw._alpha2rho0(theta_Rs=theta_Rs, Rs=Rs)
        if Rs < 0.0000001:
            Rs = 0.0000001
        f_x_prim, f_y_prim = self.nfw.nfwAlpha(R_, Rs, rho0_input, xt1, xt2)
        f_x_prim *= np.sqrt(1 - e)
        f_y_prim *= np.sqrt(1 + e)
        f_x = cos_phi*f_x_prim-sin_phi*f_y_prim
        f_y = sin_phi*f_x_prim+cos_phi*f_y_prim
        return f_x, f_y
1178318
def derivativeX(self,*args):
        '''
        
        '''
        if self.n_dims >= 4:
            j = 1
        else:
            j = 0
        if self.i_dim == j:
            return np.ones_like(*args[0])
        else:
            return np.zeros_like(*args[0])
367770
def derivatives(self, x, y, theta_E, s, q):
        
        if q >= 1:
            q = 0.999999
        psi = self._psi(x, y, q, s)
        f_x = theta_E / np.sqrt(1. - q ** 2) * np.arctan(np.sqrt(1. - q ** 2) * x / (psi+s))
        f_y = theta_E / np.sqrt(1. - q ** 2) * np.arctanh(np.sqrt(1. - q ** 2) * y / (psi + q**2*s))
        return f_x, f_y
695386
def _set_params(self):
        
        deriv0 = np.poly1d(self.eos_params)
        deriv1 = np.polyder(deriv0, 1)
        deriv2 = np.polyder(deriv1, 1)
        deriv3 = np.polyder(deriv2, 1)

        for x in np.roots(deriv1):
            if x > 0 and deriv2(x) > 0:
                v0 = x**(-3./2.)
                break
        else:
            raise EOSError("No minimum could be found")

        derivV2 = 4./9. * x**5. * deriv2(x)
        derivV3 = (-20./9. * x**(13./2.) * deriv2(x) - 8./27. *
                   x**(15./2.) * deriv3(x))
        b0 = derivV2 / x**(3./2.)
        b1 = -1 - x**(-3./2.) * derivV3 / derivV2

        # e0, b0, b1, v0
        self._params = [deriv0(v0**(-2./3.)), b0, b1, v0]
305165
def get_support(variables, polynomial):
    
    support = []
    if is_number_type(polynomial):
        support.append([0] * len(variables))
        return support
    for monomial in polynomial.expand().as_coefficients_dict():
        tmp_support = [0] * len(variables)
        mon, _ = __separate_scalar_factor(monomial)
        symbolic_support = flatten(split_commutative_parts(mon))
        for s in symbolic_support:
            if isinstance(s, Pow):
                base = s.base
                if is_adjoint(base):
                    base = base.adjoint()
                tmp_support[variables.index(base)] = s.exp
            elif is_adjoint(s):
                tmp_support[variables.index(s.adjoint())] = 1
            elif isinstance(s, (Operator, Symbol)):
                tmp_support[variables.index(s)] = 1
        support.append(tmp_support)
    return support
997844
def hessian(self):
        
        return [[[sympy.diff(partial_dv, param) for param in self.params]
                 for partial_dv in comp] for comp in self.jacobian]
611147
def _get_derivative(self, C, sa1180, vs30):
        
        derAmp = np.zeros_like(vs30)
        n = self.CONSTS['n']
        c = C['c']
        b = C['b']
        idx = vs30 < C['vlin']
        derAmp[idx] = (b * sa1180[idx] * (-1./(sa1180[idx]+c) +
                       1./(sa1180[idx] + c*(vs30[idx]/C['vlin'])**n)))
        return derAmp
333196
def derivatives(self, x, y, coeffs, beta, center_x=0, center_y=0):
        
        shapelets = self._createShapelet(coeffs)
        n_order = self._get_num_n(len(coeffs))
        dx_shapelets = self._dx_shapelets(shapelets, beta)
        dy_shapelets = self._dy_shapelets(shapelets, beta)
        n = len(np.atleast_1d(x))
        if n <= 1:
            f_x = self._shapeletOutput(x, y, beta, dx_shapelets, precalc=False)
            f_y = self._shapeletOutput(x, y, beta, dy_shapelets, precalc=False)
        else:
            H_x, H_y = self.pre_calc(x, y, beta, n_order+1, center_x, center_y)
            f_x = self._shapeletOutput(H_x, H_y, beta, dx_shapelets)
            f_y = self._shapeletOutput(H_x, H_y, beta, dy_shapelets)
        return f_x, f_y
1364948
def hessian(self, x, y, theta_E, s, q):
        
        alpha_ra, alpha_dec = self.derivatives(x, y, theta_E, s, q)
        diff = self._diff
        alpha_ra_dx, alpha_dec_dx = self.derivatives(x + diff, y, theta_E, s, q)
        alpha_ra_dy, alpha_dec_dy = self.derivatives(x, y + diff, theta_E, s, q)

        f_xx = (alpha_ra_dx - alpha_ra) / diff
        f_xy = (alpha_ra_dy - alpha_ra) / diff
        # f_yx = (alpha_dec_dx - alpha_dec)/diff
        f_yy = (alpha_dec_dy - alpha_dec) / diff
        return f_xx, f_yy, f_xy