function (k, m, P) return Jacobian of measurement function, it is
            passed into p_h.
            k (iteration number), starts at 0
            m: point where Jacobian is evaluated
            P: parameter for Jacobian, usually covariance matrix.
163316
def Ak(self, k, m_pred, P_pred):  # returns state iteration matrix
        

        return self.A[:, :, self.index[self.A_time_var_index, k]]
110902
def create_function_f_J(self):
        """"""
        return ca.Function(
            'J',
            [self.t, self.x, self.y, self.m, self.p, self.c, self.ng, self.nu],
            [ca.jacobian(self.f_x_rhs, self.x)],
            ['t', 'x', 'y', 'm', 'p', 'c', 'ng', 'nu'], ['J'], self.func_opt)
1292017
def partial_derivative(self, X, y=0):
        
        self.check_fit()

        U, V = self.split_matrix(X)

        if self.theta == 0:
            return V

        else:
            A = np.power(V, -self.theta - 1)
            B = np.power(V, -self.theta) + np.power(U, -self.theta) - 1
            h = np.power(B, (-1 - self.theta) / self.theta)
            return np.multiply(A, h) - y
287606
def _pca(x, k=2):
    ""
    x = x-torch.mean(x,0)
    U,S,V = torch.svd(x.t())
    return torch.mm(x,U[:,:k])
759036
def damping_kraus_map(p=0.10):
    
    damping_op = np.sqrt(p) * np.array([[0, 1],
                                        [0, 0]])

    residual_kraus = np.diag([1, np.sqrt(1 - p)])
    return [residual_kraus, damping_op]
220012
def create_rl_bin_method(op_name, klass, jacobians={}):
    
    # Make left and right versions.
    for prefix in ['r', '']:
        # This function will became the methods.
        op_modules = [operator, builtins]
        op_names = [op_name, op_name + '_']

        op_function_base = find_element(
            op_names,
            op_modules,
            error_on_fail=True)
        # many such functions do not take keyword arguments, so we need to wrap
        # them

        def op_function(a, b):
            return op_function_base(a, b)

        def new_method(self, other, prefix=prefix):
            if not check_special_methods():
                raise NotImplementedError(
                    'Special method %s called on %s, but special methods have been disabled. Set pymc.special_methods_available to True to enable them.' %
                    (op_name, str(self)))
            # This code will create one of two Deterministic objects.
            if prefix == 'r':
                parents = {'a': other, 'b': self}

            else:
                parents = {'a': self, 'b': other}
            jacobian_formats = {'a': 'broadcast_operation',
                                'b': 'broadcast_operation'}
            return pm.Deterministic(op_function,
                                    'A Deterministic returning the value of %s(%s,%s)' % (
                                        prefix + op_name, self.__name__, str(
                                            other)),
                                    '(' + '_'.join(
                                        [self.__name__,
                                         prefix + op_name,
                                         str(other)]) + ')',
                                    parents,
                                    trace=False,
                                    plot=False,
                                    jacobians=jacobians,
                                    jacobian_formats=jacobian_formats)
        # Convert the functions into methods for klass.
        new_method.__name__ = '__' + prefix + op_name + '__'
        setattr(
            klass,
            new_method.__name__,
            UnboundMethodType(
                new_method,
                None,
                klass))
164521
def step(self, y, u, t, h):
		
		k1 = h * self.func(t, y, u)
		k2 = h * self.func(t + .5*h, y + .5*h*k1, u)
		k3 = h * self.func(t + .5*h, y + .5*h*k2, u)
		k4 = h * self.func(t + h, y + h*k3, u)
		return y + (k1 + 2*k2 + 2*k3 + k4) / 6.0
447213
function householderTransform (x) {
  // Based on "Matrix Computations" by Golub, Van Loan
  const n = x.dimension;
  const v = x.clone();
  v.x = 1;
  let sigma = 0;
  for (let i = 1; i < n; ++i) {
    const vi = v.getComponent(i);
    sigma += vi * vi;
  }
  let beta = 0;
  if (sigma !== 0) {
    const x1 = x.x;
    const mu = _Math.sqrt(x1 * x1 + sigma);
    const v1 = (x1 <= 0 ? x1 - mu : -sigma / (x1 + mu));
    v.x = v1;
    beta = 2 * v1 * v1 / (sigma + v1 * v1);
    v.multiplyScalar(1.0 / v1);
  }
  return {
    v,
    beta
  };
}
383666
def inference(self, kern, X, likelihood, Y, mean_function=None, Y_metadata=None):
        
        assert mean_function is None, "inference with a mean function not implemented"

        # Compute K
        K = kern.K(X)

        #Find mode
        if self.bad_fhat or self.first_run:
            Ki_f_init = np.zeros_like(Y)
            self.first_run = False
        else:
            Ki_f_init = self._previous_Ki_fhat

        Ki_f_init = np.zeros_like(Y)# FIXME: take this out

        f_hat, Ki_fhat = self.rasm_mode(K, Y, likelihood, Ki_f_init, Y_metadata=Y_metadata)

        #Compute hessian and other variables at mode
        log_marginal, woodbury_inv, dL_dK, dL_dthetaL = self.mode_computations(f_hat, Ki_fhat, K, Y, likelihood, kern, Y_metadata)

        self._previous_Ki_fhat = Ki_fhat.copy()
        return Posterior(woodbury_vector=Ki_fhat, woodbury_inv=woodbury_inv, K=K), log_marginal, {'dL_dK':dL_dK, 'dL_dthetaL':dL_dthetaL}
847403
def _lmder1_freudenstein_roth():
    """"""

    def func(params, vec):
        vec[0] = -13 + params[0] + ((5 - params[1]) * params[1] - 2) * params[1]
        vec[1] = -29 + params[0] + ((1 + params[1]) * params[1] - 14) * params[1]

    def jac(params, jac):
        jac[0] = 1
        jac[1,0] = params[1] * (10 - 3 * params[1]) - 2
        jac[1,1] = params[1] * (2 + 3 * params[1]) - 14

    guess = np.asfarray([0.5, -2])

    _lmder1_driver(2, func, jac, guess,
                   0.200124960962e+02, 0.699887517585e+01,
                   [0.114124844655e+02, -0.896827913732e+00])
    _lmder1_driver(2, func, jac, guess * 10,
                   0.124328339489e+05, 0.699887517449e+01,
                   [0.114130046615e+02, -0.896796038686e+00])
    _lmder1_driver(2, func, jac, guess * 100,
                   0.11426454595762e+08, 0.699887517243e+01,
                   [0.114127817858e+02, -0.896805107492e+00])
474478
function(f, x0, options) {
    options = options || {}

    this.f = f
    this.n = x0.length
    // cache the init guess
    this.xInitGuess = x0
    this.x = MathUtils.clone(x0)

    this.df = options.df || MathUtils.getApproximateGradientFunction(this.f)
    this.maxIterator = options.maxIterator || 200
    this.err = options.err || 1E-6

    // the gradient of the function evaluated at x[k]: g[k]  (x[0] = x0)
    this.g = this.df(this.x)

    // the inverse of approximate Hessian matrix: B[k]  (B[0] = I)
    this.B = MathUtils.getIdentityMatrix(this.n)

    // direction: p[k]
    this.p = []

    this.stepsize = 0
    this.convergence = Infinity
    this.isConverges = false

    this.iterator = -1
}
1137687
def K(self, X, X2, target):
        """"""
        self._K_computations(X, X2)
        target += self.variance*self._K_dvar
478768
def K(self, X, X2=None):
        """"""
        if X2 is None:
            X2 = X

        base = np.pi * (X[:, None, :] - X2[None, :, :]) / self.period
        exp_dist = np.exp( -0.5* np.sum( np.square(  np.sin( base ) / self.lengthscale ), axis = -1 ) )

        return self.variance * exp_dist
246146
def hessian(self, x, y, amp, sigma, center_x=0, center_y=0):
        
        x_ = x - center_x
        y_ = y - center_y
        r = np.sqrt(x_**2 + y_**2)
        sigma_x, sigma_y = sigma, sigma
        if isinstance(r, int) or isinstance(r, float):
            r = max(r, self.ds)
        else:
            r[r <= self.ds] = self.ds
        d_alpha_dr = -self.d_alpha_dr(r, amp, sigma_x, sigma_y)
        alpha = self.alpha_abs(r, amp, sigma)

        f_xx = -(d_alpha_dr/r + alpha/r**2) * x_**2/r + alpha/r
        f_yy = -(d_alpha_dr/r + alpha/r**2) * y_**2/r + alpha/r
        f_xy = -(d_alpha_dr/r + alpha/r**2) * x_*y_/r
        return f_xx, f_yy, f_xy
1129582
def _lmder1_rosenbrock():
    """"""

    def func(params, vec):
        vec[0] = 10 * (params[1] - params[0]**2)
        vec[1] = 1 - params[0]

    def jac(params, jac):
        jac[0,0] = -20 * params[0]
        jac[0,1] = -1
        jac[1,0] = 10
        jac[1,1] = 0

    guess = np.asfarray([-1.2, 1])
    norm1s = [0.491934955050e+01, 0.134006305822e+04, 0.1430000511923e+06]

    for i in range(3):
        _lmder1_driver(2, func, jac, guess * 10**i,
                       norm1s[i], 0, [1, 1])
344534
def mul_at(self, other, k):
        ''''''
        if k > (self.degree + other.degree) or k > self.degree: return 0 # optimization: if the required coefficient is above the maximum coefficient of the resulting polynomial, we can already predict that and just return 0

        term = 0

        for i in _range(min(len(self), len(other))):
            coef1 = self.coefficients[-(k-i+1)]
            coef2 = other.coefficients[-(i+1)]
            if coef1 == 0 or coef2 == 0: continue # log(0) is undefined, skip (and in addition it's a nice optimization)
            term += coef1 * coef2
        return term
106260
def _mm(n_items, data, initial_params, alpha, max_iter, tol, mm_fun):
    
    if initial_params is None:
        params = np.zeros(n_items)
    else:
        params = initial_params
    converged = NormOfDifferenceTest(tol=tol, order=1)
    for _ in range(max_iter):
        nums, denoms = mm_fun(n_items, data, params)
        params = log_transform((nums + alpha) / (denoms + alpha))
        if converged(params):
            return params
    raise RuntimeError("Did not converge after {} iterations".format(max_iter))