Python numdifftools.Jacobian方法代码示例

# 需要导入模块: import numdifftools [as 别名]
# 或者: from numdifftools import Jacobian [as 别名]
def score(self, params):
        """
        Notes
        -----
        Need to generalize for AR(p) and for a constant.
        Not correct yet.  Returns numerical gradient.  Depends on package
        numdifftools.
        """
        y = self.endog
        ylag = self.exog
        nobs = self.nobs
        diffsumsq = sumofsq(y-np.dot(ylag,params))
        dsdr = 1/nobs * -2 *np.sum(ylag*(y-np.dot(ylag,params))[:,None])+\
                2*params*ylag[0]**2
        sigma2 = 1/nobs*(diffsumsq-ylag[0]**2*(1-params**2))
        gradient = -nobs/(2*sigma2)*dsdr + params/(1-params**2) + \
                1/sigma2*np.sum(ylag*(y-np.dot(ylag, params))[:,None])+\
                .5*sigma2**-2*diffsumsq*dsdr+\
                ylag[0]**2*params/sigma2 +\
                ylag[0]**2*(1-params**2)/(2*sigma2**2)*dsdr
        if self.penalty:
            pass
        j = Jacobian(self.loglike)
        return j(params)
#        return gradient 

# 需要导入模块: import numdifftools [as 别名]
# 或者: from numdifftools import Jacobian [as 别名]
def get_jacobian(function, point, minima, maxima):

    wrapper, scaled_deltas, scaled_point, orders_of_magnitude, n_dim = _get_wrapper(
        function, point, minima, maxima
    )

    # Compute the Jacobian matrix at best_fit_values
    jacobian_vector = nd.Jacobian(wrapper, scaled_deltas, method="central")(
        scaled_point
    )

    # Transform it to numpy matrix

    jacobian_vector = np.array(jacobian_vector)

    # Now correct back the Jacobian for the scales

    jacobian_vector /= orders_of_magnitude

    return jacobian_vector[0] 

# 需要导入模块: import numdifftools [as 别名]
# 或者: from numdifftools import Jacobian [as 别名]
def run_gradient_and_hessian_benchmarks(gradient_funs, hessian_funs,
                                        problem_sizes=(4, 8, 16, 32, 64, 96)):

    symbols = ('-kx', ':k>', ':k<', '--k^', '--kv', '-kp', '-ks',
               'b', '--b', '-b+', 'r', '--r', '-r+')

    results_gradients = compute_gradients(gradient_funs, problem_sizes)
    results_hessians = compute_hessians(hessian_funs, problem_sizes)

    print(results_gradients.shape)

    for i, txt in enumerate(['run times', 'errors']):
        objects = [('Jacobian ' + txt, gradient_funs, results_gradients[..., i].T),
                   ('Hessian ' + txt, hessian_funs, results_hessians[..., i].T)]
        if i == 0:
            plot_runtimes(objects, problem_sizes, symbols)
        else:
            plot_errors(objects, problem_sizes, symbols) 

# 需要导入模块: import numdifftools [as 别名]
# 或者: from numdifftools import Jacobian [as 别名]
def marginal_likelihood(self, t, T = None):
        """
        Calculate log marginal likelihood of the process under given data, using Laplace's approximation.
        This method uses a Gaussian approximation around the currently fit parameters (i.e. expects MAP
        parameters to already have been fit, e.g. through :meth:`fit`).

        :param numpy.array[float] t: Observation timestamps of the process up to time T. 1-d array of timestamps.
            must be sorted (asc)
        :param T: (optional) maximum time
        :type T: float or None

        :return: the log marginal likelihood
        :rtype: float
        """
        t, T = self._prep_t_T(t, T)

        f = self._log_posterior(t, T)
        g = self._log_posterior_grad(t, T)

        xopt = np.array(self.get_params())
        # xopt = self._fit_grad_desc(t, T).x
        H = nd.Jacobian(g)(xopt)

        return f(xopt) + 1.5 * np.log(2 * np.pi) - .5 * np.linalg.slogdet(H)[1] 

# 需要导入模块: import numdifftools [as 别名]
# 或者: from numdifftools import Jacobian [as 别名]
def score(self, params):
        """
        Score vector for Arma model
        """
        #return None
        #print params
        jac = ndt.Jacobian(self.loglike, stepMax=1e-4)
        return jac(params)[-1] 

# 需要导入模块: import numdifftools [as 别名]
# 或者: from numdifftools import Jacobian [as 别名]
def hessian(self, params):
        """
        Hessian of arma model.  Currently uses numdifftools
        """
        #return None
        Hfun = ndt.Jacobian(self.score, stepMax=1e-4)
        return Hfun(params)[-1] 

# 需要导入模块: import numdifftools [as 别名]
# 或者: from numdifftools import Jacobian [as 别名]
def example_function_jacobian_fixtures():
    def f(x):
        """f:R^3 -> R^2"""
        x1, x2, x3 = x[0], x[1], x[2]
        y1, y2 = np.sin(x1) + np.cos(x2), np.exp(x3)
        return np.array([y1, y2])

    def fprime(x):
        """Jacobian(f)(x):R^3 -> R^(2x3)"""
        x1, x2, x3 = x[0], x[1], x[2]
        jac = np.array([[np.cos(x1), -np.sin(x2), 0], [0, 0, np.exp(x3)]])
        return jac

    return {"func": f, "func_prime": fprime} 

# 需要导入模块: import numdifftools [as 别名]
# 或者: from numdifftools import Jacobian [as 别名]
def test_first_derivative_jacobian_richardson(example_function_jacobian_fixtures):
    f = example_function_jacobian_fixtures["func"]
    fprime = example_function_jacobian_fixtures["func_prime"]

    true_grad = fprime(np.ones(3))
    numdifftools_grad = Jacobian(f, order=2, n=3, method="central")(np.ones(3))
    grad = first_derivative(f, np.ones(3), n_steps=3, method="central")

    aaae(numdifftools_grad, grad)
    aaae(true_grad, grad) 

# 需要导入模块: import numdifftools [as 别名]
# 或者: from numdifftools import Jacobian [as 别名]
def V_jacobina(F, X):
    import numdifftools as nda

    V_jacobina = nda.Jacobian(F)

    return V_jacobina(X) 

# 需要导入模块: import numdifftools [as 别名]
# 或者: from numdifftools import Jacobian [as 别名]
def curl(f, x):
    jac = nd.Jacobian(f)(x)
    return sp.array([jac[1, 0] - jac[0, 1]])  # 2D curl 

# 需要导入模块: import numdifftools [as 别名]
# 或者: from numdifftools import Jacobian [as 别名]
def get_fjac(f, input_vector_convention='row'):
    '''
        Get the numerical Jacobian of the vector field function.
        If the input_vector_convention is 'row', it means that fjac takes row vectors
        as input, otherwise the input should be an array of column vectors. Note that
        the returned Jacobian would behave exactly the same if the input is an 1d array.

        The column vector convention is slightly faster than the row vector convention.
        So the matrix of row vector convention is converted into column vector convention
        under the hood.

        No matter the input vector convention, the returned Jacobian is of the following
        format:
                df_1/dx_1   df_1/dx_2   df_1/dx_3   ...
                df_2/dx_1   df_2/dx_2   df_2/dx_3   ...
                df_3/dx_1   df_3/dx_2   df_3/dx_3   ...
                ...         ...         ...         ...
    '''
    fjac = nd.Jacobian(lambda x: f(x.T).T)
    if input_vector_convention == 'row' or input_vector_convention == 0:
        def f_aux(x):
            x = x.T
            return fjac(x)

        return f_aux
    else:
        return fjac 

# 需要导入模块: import numdifftools [as 别名]
# 或者: from numdifftools import Jacobian [as 别名]
def elementwise_jacobian_transformation(fjac, X, qi, qj):
    """Inverse transform low dimension Jacobian matrix (:math:`\partial F_i / \partial x_j`) back to original space.
    The formula used to inverse transform Jacobian matrix calculated from low dimension (PCs) is:
                                            :math:`Jac = Q J Q^T`,
    where `Q, J, Jac` are the PCA loading matrix, low dimensional Jacobian matrix and the inverse transformed high
    dimensional Jacobian matrix. This function takes only one row from Q to form qi or qj.

    Parameters
    ----------
        fjac: `function`:
            The function for calculating numerical Jacobian matrix.
        X: `np.ndarray`:
            The samples coordinates with dimension n_obs x n_PCs, from which Jacobian will be calculated.
        Qi: `np.ndarray`:
            One sampled gene's PCs loading matrix with dimension n' x n_PCs, from which local dimension Jacobian matrix
            (k x k) will be inverse transformed back to high dimension.
        Qj: `np.ndarray`
            Another gene's (can be the same as those in Qi or different) PCs loading matrix with dimension  n' x n_PCs,
            from which local dimension Jacobian matrix (k x k) will be inverse transformed back to high dimension.

    Returns
    -------
        ret `np.ndarray`
            The calculated vector of Jacobian matrix (:math:`\partial F_i / \partial x_j`) for each cell.
    """

    Js = fjac(X)
    ret = np.zeros(len(X))
    for i in tqdm(range(len(X)), "calculating Jacobian for each cell"):
        J = Js[:, :, i]
        ret[i] = qi @ J @ qj
    return ret 

# 需要导入模块: import numdifftools [as 别名]
# 或者: from numdifftools import Jacobian [as 别名]
def _divergence(f, x):
    """Divergence of the reconstructed vector field function f evaluated at x"""
    jac = nd.Jacobian(f)(x)
    return np.trace(jac) 

# 需要导入模块: import numdifftools [as 别名]
# 或者: from numdifftools import Jacobian [as 别名]
def compute_divergence(f_jac, X, vectorize=True):
    """calculate divergence for many samples by taking the trace of a Jacobian matrix"""
    if vectorize:
        J = f_jac(X)
        div = np.trace(J)
    else:
        div = np.zeros(len(X))
        for i in tqdm(range(len(X)), desc="Calculating divergence"):
            J = f_jac(X[i])
            div[i] = np.trace(J)

    return div 

# 需要导入模块: import numdifftools [as 别名]
# 或者: from numdifftools import Jacobian [as 别名]
def _curl(f, x):
    """Curl of the reconstructed vector field f evaluated at x in 3D"""
    jac = nd.Jacobian(f)(x)
    return np.array([jac[2, 1] - jac[1, 2], jac[0, 2] - jac[2, 0], jac[1, 0] - jac[0, 1]])