Python misc.derivative函数代码示例

def solarAzimuth(_day, _hour = None):
    """Solar azimuth in radians"""  
    def cos_azmth(_hourAngle):
        """Definition of the cosine of the azimuth replacing the altitude by 
        it formula"""
        return (np.sin(np.arcsin(np.cos(latitude)*np.cos(declination(_day,radianes))* \
        np.cos(_hourAngle)+np.sin(latitude)*np.sin(declination(_day,radianes))))*np.sin(latitude)- \
        np.sin(declination(_day,radianes))) \
        / (np.cos(np.arcsin(np.cos(latitude)*np.cos(declination(_day,radianes))* \
        np.cos(_hourAngle)+np.sin(latitude)*np.sin(declination(_day,radianes))))*np.cos(latitude))

    if _hour is None:
        for i in range(np.size(hourAngle)):
            derivada[i] = sc.derivative(cos_azmth,hourAngle[i],dx=1e-6)
            _signo = -np.sign(derivada[i])
            azimuth[i] =_signo * np.arccos(cos_azmth(hourAngle[i]))
            azimuth_deg[i] = azimuth[i]*180/np.pi
            
    else:
         _hourAngle = (12-int(_hour))*np.pi/12
         _derivada = sc.derivative(cos_azmth,_hourAngle,dx=1e-6)
         _signo = - np.sign(_derivada)
         if cos_azmth(_hourAngle) > 1:
             _azimuth = _signo * np.arccos(1) * 180 / np.pi
         else:
             _azimuth = _signo * np.arccos(cos_azmth(_hourAngle)) * 180 / np.pi
         return _azimuth
    
    return

		def fisher(theta, i, j = None):
			"""
			Fisher information using the first order derivative

			:param theta: the theta of the density
			:param i: The ith component of the diagonal of the fisher information matrix will be returned (if j is None)
			:param j: The i,j th component of the fisher information matrix will be returned
			"""

			#Bring it in a form that we can derive
			fh = lambda ti, t0, tn, x: np.log(self.density(x, list(t0) + [ti] + list(tn)))

			# The derivative
			f_d_theta_i = lambda x: derivative(fh, theta[i], dx=1e-5, n=1, args=(theta[0:i], theta[i + 1:], x))

			if j is not None:
				f_d_theta_j = lambda x: derivative(fh, theta[j], dx=1e-5, n=1, args=(theta[0:j], theta[j + 1:], x))
				f = lambda x: np.float128(0) if fabs(self.density(x, theta)) < 1e-5 else f_d_theta_i(x) * f_d_theta_j(x) * self.density(x, theta)
			else:
				# The function to integrate
				f = lambda x: np.float128(0) if fabs(self.density(x, theta)) < 1e-5 else f_d_theta_i(x) ** 2 * self.density(x, theta)


			#First order
			result = integrate(f, self.support)
			return result

def newton(fn=0, fp=0, x=0.1, mn=0):
    t = time()
    test = 0
    if not fn and not fp:
        from scipy.misc import derivative

        func = input("equation(use x as unknown):")
        fn = lambda x: eval(func)
        fp = lambda x: derivative(fn, x)
        print("div")
    elif not fp:
        from scipy.misc import derivative

        fp = lambda x: derivative(fn, x)
    if mn:
        while fn(x, mn) != 0:
            if fp(x, mn) == 0:
                x += 1
                print("error")
                test += 1
                print(test)
                break
            x, lastx = x - (fn(x, mn) / fp(x, mn)), x
    else:
        while fn(x) != 0 or rerr > 10 ** -12:
            x, lastx = x - (fn(x) / fp(x)), x
            rerr = abs(lastx - x) / x
    return x

    def eval_function(self, individual, radius, representatives):
        distanceSum = 0.001+sum([self.share_function(self.distance(individual, r), 1, radius) for r in representatives])
        # print "distanceSum: %s" % distanceSum
        calcWithY = partial(self.function_object.calculateWithXY, individual[0])
        calcWithX = partial(self.function_object.calculateWithYX, individual[1])
        derivY = abs(misc.derivative(calcWithY, individual[1], dx=1e-6))
        derivX = abs(misc.derivative(calcWithX, individual[0], dx=1e-6))

        return self.function_object.calculate(individual)/distanceSum, distanceSum*derivX, distanceSum*derivY

def calc_curvature(func, d_point=0):
    #曲率半径、曲率を求める
    #曲率が大きいほど関数の曲がり具合が大きい

    dx1 = scmisc.derivative(func, d_point, n=1, dx=1e-6) #dxを指定しないとダメ
    dx2 = scmisc.derivative(func, d_point, n=2, dx=1e-6)
    R = (1.0 + (dx1**2.0))**(3.0/2.0) / np.fabs(dx2) #曲率半径
    curvature_of_func = 1.0 / R

    return curvature_of_func

def C_Meng(T, Tc, Pc, D, B):
    r"""Calculate the 3rd virial coefficient using the Meng-Duan-Li correlation

    Parameters
    ----------
    T : float
        Temperature [K]
    Tc : float
        Critical temperature [K]
    Pc : float
        Critical pressure
    D : float
        dipole moment [debye]
    B : list
        Second virial coefficient tuple with B, ∂B/∂T, ∂²B/∂T²

    Returns
    -------
    C : float
        Third virial coefficient [m⁶/mol²]
    C1 : float
        T(∂C/∂T) [m⁶/mol²]
    C2 : float
        T²(∂²C/∂T²) [m⁶/mol²]

    Examples
    --------
    Selected date from Table 2, pag 74 for neon

    >>> from lib.mEoS import Ne
    >>> D = Ne.momentoDipolar.Debye
    >>> B = B_Meng(273.15, Ne.Tc, Ne.Pc, Ne.f_acent, D)
    >>> "%.2f" % (C_Meng(273.15, Ne.Tc, Ne.Pc, D, B)[0]*1e-5)
    '0.24'

    References
    ----------
    [7]_ Meng, L., Duan, Y.Y. Li, L. Correations for Second and Third Virial
        Coefficients of Pure Fluids. Fluid Phase Equilibria 226 (2004) 109-120
    """
    mur = D**2*Pc/1.01325/Tc**2

    def f(T, *args):
        B = args[-1]
        Br = B*Pc/R/Tc
        Tr = T/Tc
        f0 = 1094.051 - 3334.145/Tr**0.1 + 3389.848/Tr**0.2 - 1149.58/Tr**0.3
        f1 = (2.0243-0.85902/Tr)*1e-10
        return 5.476e-3 + (Br-0.0936)**2*(f0+mur**4*f1)

    C = f(T, B[0])
    C1 = derivative(f, T, n=1, args=(T, B[1]))
    C2 = derivative(f, T, n=2, args=(T, B[2]))

    return C, C1, C2

    def gruneisen_parameter(self, temperature, volume):
        """
        Slater-gamma formulation(the default):
            gruneisen paramter = - d log(theta)/ d log(V)
                               = - ( 1/6 + 0.5 d log(B)/ d log(V) )
                               = - (1/6 + 0.5 V/B dB/dV),
                                    where dB/dV = d^2E/dV^2 + V * d^3E/dV^3

        Mie-gruneisen formulation:
            Eq(31) in doi.org/10.1016/j.comphy.2003.12.001
            Eq(7) in Blanco et. al. Joumal of Molecular Structure (Theochem)
                368 (1996) 245-255
            Also se J.P. Poirier, Introduction to the Physics of the Earth’s
                Interior, 2nd ed. (Cambridge University Press, Cambridge,
                2000) Eq(3.53)

        Args:
            temperature (float): temperature in K
            volume (float): in Ang^3

        Returns:
            float: unitless
        """
        if isinstance(self.eos, PolynomialEOS):
            p = np.poly1d(self.eos.eos_params)
            # first derivative of energy at 0K wrt volume evaluated at the
            # given volume, in eV/Ang^3
            dEdV = np.polyder(p, 1)(volume)
            # second derivative of energy at 0K wrt volume evaluated at the
            # given volume, in eV/Ang^6
            d2EdV2 = np.polyder(p, 2)(volume)
            # third derivative of energy at 0K wrt volume evaluated at the
            # given volume, in eV/Ang^9
            d3EdV3 = np.polyder(p, 3)(volume)
        else:
            func = self.ev_eos_fit.func
            dEdV = derivative(func, volume, dx=1e-3)
            d2EdV2 = derivative(func, volume, dx=1e-3, n=2, order=5)
            d3EdV3 = derivative(func, volume, dx=1e-3, n=3, order=7)

        # Mie-gruneisen formulation
        if self.use_mie_gruneisen:
            p0 = dEdV
            return (self.gpa_to_ev_ang * volume *
                    (self.pressure + p0 / self.gpa_to_ev_ang) /
                    self.vibrational_internal_energy(temperature, volume))

        # Slater-gamma formulation
        # first derivative of bulk modulus wrt volume, eV/Ang^6
        dBdV = d2EdV2 + d3EdV3 * volume
        return -(1./6. + 0.5 * volume * dBdV /
                 FloatWithUnit(self.ev_eos_fit.b0_GPa, "GPa").to("eV ang^-3"))

def radius_of_curvature(theta, func_R, *args_for_func_R):
    """Returns R_c = (R^2 + R'^2)^{3/2} / |R^2 + 2 R'^2 - R R''| 

Uses numerical differentiation.  NOT RECOMMENDED SINCE NOT ACCURATE ON
THE AXIS.  Use `axis_Rc` instead.

    """
    R = func_R(theta, *args_for_func_R)
    dR = derivative(func_R, theta,
                    dx=DX_FOR_NUMERICAL_DERIVATIVE, args=args_for_func_R)
    d2R = derivative(func_R, theta, n=2,
                     dx=DX_FOR_NUMERICAL_DERIVATIVE, args=args_for_func_R)
    return (R**2 + dR**2)**1.5 / np.abs(R**2 + 2*dR**2 - R*d2R)

def C_Orbey_Vera(T, Tc, Pc, w):
    """Calculate the third virial coefficient using the Orbey-Vera correlation

    Parameters
    ----------
    T : float
        Temperature [K]
    Tc : float
        Critical temperature [K]
    Pc : float
        Critical pressure
    w : float
        Acentric factor [-]

    Returns
    -------
    C : float
        Third virial coefficient [m⁶/mol²]
    C1 : float
        T(∂C/∂T) [m⁶/mol²]
    C2 : float
        T²(∂²C/∂T²) [m⁶/mol²]

    Examples
    --------
    Selected points from Table 2 of paper

    >>> from lib.mEoS.Benzene import Benzene as Bz
    >>> "%.1f" % (C_Orbey_Vera(0.877*Bz.Tc, Bz.Tc, Bz.Pc, Bz.f_acent)[0]*1e9)
    '41.5'
    >>> "%.1f" % (C_Orbey_Vera(1.019*Bz.Tc, Bz.Tc, Bz.Pc, Bz.f_acent)[0]*1e9)
    '35.8'

    References
    ----------
    [4]_ Orbey, H., Vera, J.H.: Correlation for the third virial coefficient
        using Tc, Pc and ω as parameters, AIChE Journal 29, 107 (1983)
    """
    def f(T):
        Tr = T/Tc
        g0 = 0.01407+0.02432/Tr**2.8-0.00313/Tr**10.5
        g1 = -0.02676+0.0177/Tr**2.8+0.04/Tr**3-0.003/Tr**6-0.00228/Tr**10.5
        g = g0+w*g1
        return g*R**2*Tc**2/Pc**2

    C = f(T)
    C1 = derivative(f, T, n=1)
    C2 = derivative(f, T, n=2)

    return C, C1, C2

def derivative_check(x,popt,left,right):
	if len(left)==0 or len(right)==0:
		return True
	if len(left)==1 or len(right)==1:
		return True
	maxlx=np.max(left[:,0])
	minlx=np.min(left[:,0])
	maxly=np.max(left[:,1])
	minly=np.min(left[:,1])
	maxrx=np.max(right[:,0])
	minrx=np.min(right[:,0])
	maxry=np.max(right[:,1])
	minry=np.min(right[:,1])
	t1=-1
	t2=-1
	pend1=[]
	pend2=[]
	pend1.append(deriv.derivative(fu,minlx,n=1,args=popt))
	pend1.append(deriv.derivative(fu,maxlx,n=1,args=popt))
	pend1.append(deriv.derivative(fu,minrx,n=1,args=popt))
	pend1.append(deriv.derivative(fu,maxrx,n=1,args=popt))
	sign=pend1[0]
	for i in pend1:
		if i<0:
			return False
	pend1=[]
	pend1.append(deriv.derivative(fu,minlx,n=2,args=popt))
	pend1.append(deriv.derivative(fu,maxlx,n=2,args=popt))
	pend1.append(deriv.derivative(fu,minrx,n=2,args=popt))
	pend1.append(deriv.derivative(fu,maxrx,n=2,args=popt))
	sign=pend1[0]
	for i in pend1:
		if (sign>0 and i<0) or (sign<0 and i>0):
			return False
	return True

    def get_corr_pred(self, eps, d_eps, sig, t_n, t_n1, alpha, q, kappa):
        #         g = lambda k: 0.8 - 0.8 * np.exp(-k)
        #         g = lambda k: 1. / (1 + np.exp(-2 * k + 6.))
        n_e, n_ip, n_s = eps.shape
        D = np.zeros((n_e, n_ip, 3, 3))
        D[:, :, 0, 0] = self.E_m
        D[:, :, 2, 2] = self.E_f
        sig_trial = sig[:, :, 1]/(1-self.g(kappa)) + self.E_b * d_eps[:,:, 1]
        xi_trial = sig_trial - q
        f_trial = abs(xi_trial) - (self.sigma_y + self.K_bar * alpha)
        elas = f_trial <= 1e-8
        plas = f_trial > 1e-8
        d_sig = np.einsum('...st,...t->...s', D, d_eps)
        sig += d_sig

        d_gamma = f_trial / (self.E_b + self.K_bar + self.H_bar) * plas
        alpha += d_gamma
        kappa += d_gamma
        q += d_gamma * self.H_bar * np.sign(xi_trial)
        w = self.g(kappa)

        sig_e = sig_trial - d_gamma * self.E_b * np.sign(xi_trial)
        sig[:, :, 1] = (1-w)*sig_e

        E_p = -self.E_b / (self.E_b + self.K_bar + self.H_bar) * derivative(self.g, kappa, dx=1e-6) * sig_e \
            + (1 - w) * self.E_b * (self.K_bar + self.H_bar) / \
            (self.E_b + self.K_bar + self.H_bar)

        D[:, :, 1, 1] = (1-w)*self.E_b*elas + E_p*plas

        return sig, D, alpha, q, kappa

  def _testCompareToExplicitDerivative(self, dtype):
    """Compare to the explicit reparameterization derivative.

    Verifies that the computed derivative satisfies
    dsample / dalpha = d igammainv(alpha, u) / dalpha,
    where u = igamma(alpha, sample).

    Args:
      dtype: TensorFlow dtype to perform the computations in.
    """
    delta = 1e-3
    np_dtype = dtype.as_numpy_dtype
    try:
      from scipy import misc  # pylint: disable=g-import-not-at-top
      from scipy import special  # pylint: disable=g-import-not-at-top

      alpha_val = np.logspace(-2, 3, dtype=np_dtype)
      alpha = constant_op.constant(alpha_val)
      sample = random_ops.random_gamma([], alpha, np_dtype(1.0), dtype=dtype)
      actual = gradients_impl.gradients(sample, alpha)[0]

      (sample_val, actual_val) = self.evaluate((sample, actual))

      u = special.gammainc(alpha_val, sample_val)
      expected_val = misc.derivative(
          lambda alpha_prime: special.gammaincinv(alpha_prime, u),
          alpha_val, dx=delta * alpha_val)

      self.assertAllClose(actual_val, expected_val, rtol=1e-3, atol=1e-3)
    except ImportError as e:
      tf_logging.warn("Cannot use special functions in a test: %s" % str(e))

def func_2nd_derivative_numerical(x):
    value = 1.
### START YOUR CODE HERE ###
    value = misc.derivative(func,x,1e-5,2);

#### END YOUR CODE HERE ####
    return value

def func_first_derivative(x):
    value = 1.
### START YOUR CODE HERE ###
    value = misc.derivative(func,x,1e-5);

#### END YOUR CODE HERE ####
    return value

    def simulateDeltaWithTime(self, timeSteps):
        opt=copy.copy(self.option)
        spotPrice=copy.copy(opt.S0)
        mechanism=self.mechanism
        
        plotDf=pd.DataFrame(np.zeros([len(timeSteps), 3]), columns=['TimeToMaturity', 'DA_Delta', 'BLS_Delta'])
        plotDf['TimeToMaturity']=timeSteps
        
        def optPrice(assetPrice):
            opt.S0=assetPrice

            orders=self.traders.getOrders(opt, self.numOrders)
            orders=mechanism.clearOrders(orders)
       
            return mechanism.getPrice(orders)
        
        for i, p in enumerate(timeSteps):
            opt.T=p            
            plotDf.ix[i]['DA_Delta']=derivative(optPrice, x0=spotPrice, dx=20)
            
            opt.S0=spotPrice
            plotDf.ix[i]['BLS_Delta']=opt.delta()
            
            print i
     
        return plotDf