Calculates the gradient of the joint log posterior with respect to self.
        Calculation of the log posterior is restricted to the variables in calculation_set.
773813
def logp_gradient_contribution(self, calculation_set=None):
        
        # NEED some sort of check to see if the log p calculation has recently
        # failed, in which case not to continue

        return self.logp_partial_gradient(self, calculation_set) + builtins.sum(
            [child.logp_partial_gradient(self, calculation_set) for child in self.children])
330579
def calculate_gradient(self, batch_info, device, model, rollout):
        
        assert isinstance(rollout, Trajectories), "ACER algorithm requires trajectory input"

        local_epsilon = 1e-6

        evaluator = model.evaluate(rollout)

        actions = evaluator.get('rollout:actions')
        rollout_probabilities = torch.exp(evaluator.get('rollout:logprobs'))

        # We calculate the trust-region update with respect to the average model
        if self.trust_region:
            self.update_average_model(model)

        logprobs = evaluator.get('model:logprobs')
        q = evaluator.get('model:q')

        # Selected action values
        action_logprobs = select_indices(logprobs, actions)
        action_q = select_indices(q, actions)

        # We only want to propagate gradients through specific variables
        with torch.no_grad():
            model_probabilities = torch.exp(logprobs)

            # Importance sampling correction - we must find the quotient of probabilities
            rho = model_probabilities / (rollout_probabilities + local_epsilon)

            # Probability quotient only for selected actions
            actions_rho = select_indices(rho, actions)

            # Calculate policy state values
            model_state_values = (model_probabilities * q).sum(dim=1)

            trajectory_rewards = rollout.transition_tensors['rewards']
            trajectory_dones = rollout.transition_tensors['dones']

            q_retraced = self.retrace(
                trajectory_rewards,
                trajectory_dones,
                action_q.reshape(trajectory_rewards.size()),
                model_state_values.reshape(trajectory_rewards.size()),
                actions_rho.reshape(trajectory_rewards.size()),
                rollout.rollout_tensors['final_values']
            ).flatten()

            advantages = q_retraced - model_state_values
            importance_sampling_coefficient = torch.min(actions_rho, self.rho_cap * torch.ones_like(actions_rho))

            explained_variance = 1 - torch.var(q_retraced - action_q) / torch.var(q_retraced)

        # Entropy of the policy distribution
        policy_entropy = torch.mean(model.entropy(logprobs))
        policy_gradient_loss = -torch.mean(advantages * importance_sampling_coefficient * action_logprobs)

        # Policy gradient bias correction
        with torch.no_grad():
            advantages_bias_correction = q - model_state_values.view(model_probabilities.size(0), 1)
            bias_correction_coefficient = F.relu(1.0 - self.rho_cap / (rho + local_epsilon))

        # This sum is an expectation with respect to action probabilities according to model policy
        policy_gradient_bias_correction_gain = torch.sum(
            logprobs * bias_correction_coefficient * advantages_bias_correction * model_probabilities,
            dim=1
        )

        policy_gradient_bias_correction_loss = - torch.mean(policy_gradient_bias_correction_gain)

        policy_loss = policy_gradient_loss + policy_gradient_bias_correction_loss

        q_function_loss = 0.5 * F.mse_loss(action_q, q_retraced)

        if self.trust_region:
            with torch.no_grad():
                average_evaluator = self.average_model.evaluate(rollout)
                average_action_logits = average_evaluator.get('model:logprobs')

            actor_loss = policy_loss - self.entropy_coefficient * policy_entropy
            q_loss = self.q_coefficient * q_function_loss

            actor_gradient = torch.autograd.grad(-actor_loss, logprobs, retain_graph=True)[0]

            # kl_divergence = model.kl_divergence(average_action_logits, action_logits).mean()
            # kl_divergence_grad = torch.autograd.grad(kl_divergence, action_logits, retain_graph=True)

            # Analytically calculated derivative of KL divergence on logits
            # That makes it hardcoded for discrete action spaces
            kl_divergence_grad_symbolic = - torch.exp(average_action_logits) / logprobs.size(0)

            k_dot_g = (actor_gradient * kl_divergence_grad_symbolic).sum(dim=-1)
            k_dot_k = (kl_divergence_grad_symbolic ** 2).sum(dim=-1)

            adjustment = (k_dot_g - self.trust_region_delta) / k_dot_k
            adjustment_clipped = adjustment.clamp(min=0.0)

            actor_gradient_updated = actor_gradient - adjustment_clipped.view(adjustment_clipped.size(0), 1)

            # Populate gradient from the newly updated fn
            logprobs.backward(gradient=-actor_gradient_updated, retain_graph=True)
            q_loss.backward(retain_graph=True)
        else:
            # Just populate gradient from the loss
            loss = policy_loss + self.q_coefficient * q_function_loss - self.entropy_coefficient * policy_entropy

            loss.backward()

        return {
            'policy_loss': policy_loss.item(),
            'policy_gradient_loss': policy_gradient_loss.item(),
            'policy_gradient_bias_correction': policy_gradient_bias_correction_loss.item(),
            'avg_q_selected': action_q.mean().item(),
            'avg_q_retraced': q_retraced.mean().item(),
            'q_loss': q_function_loss.item(),
            'policy_entropy': policy_entropy.item(),
            'advantage_norm': torch.norm(advantages).item(),
            'explained_variance': explained_variance.item(),
            'model_prob_std': model_probabilities.std().item(),
            'rollout_prob_std': rollout_probabilities.std().item()
        }
143555
def cv_gradient(self,z):
        
        gradient = np.zeros(np.sum(self.approx_param_no))
        z_t = z.T      
        log_q = self.normal_log_q(z.T)
        log_p = self.log_p(z.T)
        grad_log_q = self.grad_log_q(z)
        gradient = grad_log_q*(log_p-log_q)

        alpha0 = alpha_recursion(np.zeros(np.sum(self.approx_param_no)), grad_log_q, gradient, np.sum(self.approx_param_no))           

        vectorized = gradient - ((alpha0/np.var(grad_log_q,axis=1))*grad_log_q.T).T

        return np.mean(vectorized,axis=1)
314955
def LOO(self, kern, X, Y, likelihood, posterior, Y_metadata=None, K=None, f_hat=None, W=None, Ki_W_i=None):
        
        Ki_f_init = np.zeros_like(Y)

        if K is None:
            K = kern.K(X)

        if f_hat is None:
            f_hat, _ = self.rasm_mode(K, Y, likelihood, Ki_f_init, Y_metadata=Y_metadata)

        if W is None:
            W = -likelihood.d2logpdf_df2(f_hat, Y, Y_metadata=Y_metadata)

        if Ki_W_i is None:
            _, _, _, Ki_W_i = self._compute_B_statistics(K, W, likelihood.log_concave)

        logpdf_dfhat = likelihood.dlogpdf_df(f_hat, Y, Y_metadata=Y_metadata)

        if W.shape[1] == 1:
            W = np.diagflat(W)

        #Eq 14, and 16
        var_site = 1./np.diag(W)[:, None]
        mu_site = f_hat + var_site*logpdf_dfhat
        prec_site = 1./var_site
        #Eq 19
        marginal_cov = Ki_W_i
        marginal_mu = marginal_cov.dot(np.diagflat(prec_site)).dot(mu_site)
        marginal_var = np.diag(marginal_cov)[:, None]
        #Eq 30 with using site parameters instead of Gaussian site parameters
        #(var_site instead of sigma^{2} )
        posterior_cav_var = 1./(1./marginal_var - 1./var_site)
        posterior_cav_mean = posterior_cav_var*((1./marginal_var)*marginal_mu - (1./var_site)*Y)

        flat_y = Y.flatten()
        flat_mu = posterior_cav_mean.flatten()
        flat_var = posterior_cav_var.flatten()

        if Y_metadata is not None:
            #Need to zip individual elements of Y_metadata aswell
            Y_metadata_flat = {}
            if Y_metadata is not None:
                for key, val in Y_metadata.items():
                    Y_metadata_flat[key] = np.atleast_1d(val).reshape(-1, 1)

            zipped_values = []

            for i in range(Y.shape[0]):
                y_m = {}
                for key, val in Y_metadata_flat.items():
                    if np.isscalar(val) or val.shape[0] == 1:
                        y_m[key] = val
                    else:
                        #Won't broadcast yet
                        y_m[key] = val[i]
                zipped_values.append((flat_y[i], flat_mu[i], flat_var[i], y_m))
        else:
            #Otherwise just pass along None's
            zipped_values = zip(flat_y, flat_mu, flat_var, [None]*Y.shape[0])

        def integral_generator(yi, mi, vi, yi_m):
            def f(fi_star):
                #More stable in the log space
                p_fi = np.exp(likelihood.logpdf(fi_star, yi, yi_m)
                              - 0.5*np.log(2*np.pi*vi)
                              - 0.5*np.square(mi-fi_star)/vi)
                return p_fi
            return f

        #Eq 30
        p_ystar, _ = zip(*[quad(integral_generator(y, m, v, yi_m), -np.inf, np.inf)
                           for y, m, v, yi_m in zipped_values])
        p_ystar = np.array(p_ystar).reshape(-1, 1)
        return np.log(p_ystar)
1129541
@Override
	public void computeGradient(S jacLeft , S jacRight ,
								DMatrixRMaj residuals, DMatrixRMaj gradient) {

		// Find the gradient using the two matrices for Jacobian
		// g = J'*r = [L,R]'*r
		x1.reshape(jacLeft.getNumCols(),1);
		x2.reshape(jacRight.getNumCols(),1);
		multTransA(jacLeft,residuals,x1);
		multTransA(jacRight,residuals,x2);

		CommonOps_DDRM.insert(x1,gradient,0,0);
		CommonOps_DDRM.insert(x2,gradient,x1.numRows,0);
	}
933194
def update_gradients_full(self, dL_dK, X, X2=None, reset=True):
        
        self.variance.gradient = np.sum(self.K(X, X2)* dL_dK)/self.variance

        #now the lengthscale gradient(s)
        dL_dr = self.dK_dr_via_X(X, X2) * dL_dK
        if self.ARD:

            tmp = dL_dr*self._inv_dist(X, X2)
            if X2 is None: X2 = X
            if use_stationary_cython:
                self.lengthscale.gradient = self._lengthscale_grads_cython(tmp, X, X2)
            else:
                self.lengthscale.gradient = self._lengthscale_grads_pure(tmp, X, X2)
        else:
            r = self._scaled_dist(X, X2)
            self.lengthscale.gradient = -np.sum(dL_dr*r)/self.lengthscale
515222
def _forward_step(self, Y, T_X, Var_X, Var_dX, rho_e, sigma2_e, weight):
        
        # We currently only implement diagonal form
        # of covariance matrix for Var_X, Var_dX and T_X, which means
        # each dimension of X is independent and their innovation noise
        # are also independent. Note that log_p_data takes this assumption.
        if Var_X.ndim == 1:
            inv_Var_X = np.diag(1 / Var_X)
            half_log_det_Var_X = np.sum(np.log(Var_X)) / 2.0
            Var_X = np.diag(Var_X)
            # the marginal variance of X
        else:
            half_log_det_Var_X = self._half_log_det(Var_X)
            inv_Var_X = np.linalg.inv(Var_X)
        if Var_dX.ndim == 1:
            inv_Var_dX = np.diag(1 / Var_dX)
            half_log_det_Var_dX = np.sum(np.log(Var_dX)) / 2.0
            Var_dX = np.diag(Var_dX)
            # the marginal variance of Delta X (the change of X from
            # previous time point)
        else:
            inv_Var_dX = np.linalg.inv(Var_dX)
            half_log_det_Var_dX = self._half_log_det(Var_dX)
        if T_X.ndim == 1:
            T_X = np.diag(T_X)
            # Transfer function of X: the expected mean of X at t+1
            # time point is T_x * X
        [n_T, n_V] = np.shape(Y)
        # numbers of time points and voxels
        mu = [None] * n_T
        # posterior mean of X, conditioned on all data up till the current
        # time point
        Gamma_inv = [None] * n_T
        # inverse of poterior Gamma.
        mu_Gamma_inv = [None] * n_T
        # mu * inv(Gamma)
        log_p_data = - np.log(np.pi * 2) * (n_T * n_V) / 2 \
            - half_log_det_Var_X - np.sum(np.log(sigma2_e)) * n_T / 2.0\
            + np.sum(np.log(1 - rho_e**2)) / 2.0 - half_log_det_Var_dX \
            * (n_T - 1)
        # This is the term to be incremented by c_n at each time step.
        # We first add all the fixed terms to it.

        # The following are a few fixed terms.
        Lambda_0 = np.dot(T_X, np.dot(inv_Var_dX, T_X.T)) \
            + np.dot(weight * rho_e**2 / sigma2_e, weight.T)
        H = np.dot(inv_Var_dX, T_X.T) + np.dot(weight * rho_e / sigma2_e,
                                               weight.T)
        Lambda_1 = inv_Var_dX + np.dot(weight / sigma2_e, weight.T)

        Gamma_inv[0] = inv_Var_X + np.dot(
            weight * (1 - rho_e**2) / sigma2_e, weight.T)
        # We might not need this and only use linalg.solve for related terms.
        mu_Gamma_inv[0] = np.dot(
            Y[0, :] * (1 - rho_e**2) / sigma2_e, weight.T)
        mu[0] = np.linalg.solve(Gamma_inv[0], mu_Gamma_inv[0])
        log_p_data -= 0.5 * np.sum(Y[0, :]**2 * (1 - rho_e**2) / sigma2_e)
        # This is the term added for the first time point.

        deltaY = Y[1:, :] - rho_e * Y[:-1, :]
        deltaY_sigma2inv_rho_weightT = np.dot(
            deltaY / sigma2_e * rho_e, weight.T)
        for t in np.arange(1, n_T):
            Gamma_tilde_inv = Lambda_0 + Gamma_inv[t - 1]
            tmp = np.linalg.solve(Gamma_tilde_inv, H.T)
            Gamma_inv[t] = Lambda_1 - np.dot(H, tmp)
            mu_Gamma_inv[t] = np.dot(deltaY[t - 1, :] / sigma2_e, weight.T) \
                + np.dot(mu_Gamma_inv[t - 1]
                         - deltaY_sigma2inv_rho_weightT[t - 1, :], tmp)
            mu[t] = np.linalg.solve(Gamma_inv[t], mu_Gamma_inv[t])
            tmp2 = mu_Gamma_inv[t - 1] - deltaY_sigma2inv_rho_weightT[t - 1, :]
            log_p_data += -self._half_log_det(Gamma_tilde_inv) \
                + np.dot(tmp2, np.linalg.solve(Gamma_tilde_inv, tmp2)) / 2.0
        log_p_data += -self._half_log_det(Gamma_inv[-1]) \
            + np.dot(mu_Gamma_inv[-1], mu[-1]) / 2.0 \
            - np.sum(deltaY**2 / sigma2_e) / 2.0
        return mu, mu_Gamma_inv, Gamma_inv, log_p_data, Lambda_0, \
            Lambda_1, H, deltaY, deltaY_sigma2inv_rho_weightT
366856
def _deriv_logaddexp2(x1, x2):
    """"""
    y1 = np.exp2(x1)
    y2 = np.exp2(x2)
    df_dx1 = y1 / (y1 + y2)
    df_dx2 = y2 / (y1 + y2)
    return np.vstack([df_dx1, df_dx2]).T
1119027
def inference(self, kern, X, likelihood, Y, mean_function=None, Y_metadata=None, K=None, variance=None, Z_tilde=None):
        

        if mean_function is None:
            m = 0
        else:
            m = mean_function.f(X)

        if variance is None:
            variance = likelihood.gaussian_variance(Y_metadata)

        YYT_factor = Y-m

        if K is None:
            K = kern.K(X)

        Ky = K.copy()
        diag.add(Ky, variance+1e-8)

        Wi, LW, LWi, W_logdet = pdinv(Ky)

        alpha, _ = dpotrs(LW, YYT_factor, lower=1)

        log_marginal =  0.5*(-Y.size * log_2_pi - Y.shape[1] * W_logdet - np.sum(alpha * YYT_factor))

        if Z_tilde is not None:
            # This is a correction term for the log marginal likelihood
            # In EP this is log Z_tilde, which is the difference between the
            # Gaussian marginal and Z_EP
            log_marginal += Z_tilde

        dL_dK = 0.5 * (tdot(alpha) - Y.shape[1] * Wi)

        dL_dthetaL = likelihood.exact_inference_gradients(np.diag(dL_dK), Y_metadata)

        return Posterior(woodbury_chol=LW, woodbury_vector=alpha, K=K), log_marginal, {'dL_dK':dL_dK, 'dL_dthetaL':dL_dthetaL, 'dL_dm':alpha}
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
88439
def logp_partial_gradient(self, variable, calculation_set=None):
        
        if (calculation_set is None) or (self in calculation_set):

            if not datatypes.is_continuous(variable):
                return zeros(shape(variable.value))

            if variable is self:
                try:
                    gradient_func = self._logp_partial_gradients['value']

                except KeyError:
                    raise NotImplementedError(
                        repr(
                            self) +
                        " has no gradient function for 'value'")

                gradient = np.reshape(
                    gradient_func.get(
                    ),
                    np.shape(
                        variable.value))
            else:
                gradient = builtins.sum(
                    [self._pgradient(variable,
                                     parameter,
                                     value) for parameter,
                     value in six.iteritems(self.parents)])

            return gradient
        else:
            return 0
518588
def _backward_step(self, deltaY, deltaY_sigma2inv_rho_weightT,
                       sigma2_e, weight, mu, mu_Gamma_inv, Gamma_inv,
                       Lambda_0, Lambda_1, H):
        
        n_T = len(Gamma_inv)
        # All the terms with hat before are parameters of posterior
        # distributions of X conditioned on data from all time points,
        # whereas the ones without hat calculated by _forward_step
        # are mean and covariance of posterior of X conditioned on
        # data up to the time point.
        Gamma_inv_hat = [None] * n_T
        mu_Gamma_inv_hat = [None] * n_T
        mu_hat = [None] * n_T

        mu_hat[-1] = mu[-1].copy()
        mu_Gamma_inv_hat[-1] = mu_Gamma_inv[-1].copy()
        Gamma_inv_hat[-1] = Gamma_inv[-1].copy()

        for t in np.arange(n_T - 2, -1, -1):
            tmp = np.linalg.solve(Gamma_inv_hat[t + 1] - Gamma_inv[t + 1]
                                  + Lambda_1, H)
            Gamma_inv_hat[t] = Gamma_inv[t] + Lambda_0 - np.dot(H.T, tmp)
            mu_Gamma_inv_hat[t] = mu_Gamma_inv[t] \
                - deltaY_sigma2inv_rho_weightT[t, :] + np.dot(
                    mu_Gamma_inv_hat[t + 1] - mu_Gamma_inv[t + 1]
                    + np.dot(deltaY[t, :] / sigma2_e, weight.T), tmp)
            mu_hat[t] = np.linalg.solve(Gamma_inv_hat[t],
                                        mu_Gamma_inv_hat[t])
        return mu_hat, mu_Gamma_inv_hat, Gamma_inv_hat
907830
def posterior(self, x):
        """"""
        return Gamma(a=self.a + 0.5 * x.size,
                     b=self.b + 0.5 * np.sum(x**2))
743193
def sample_gradient(x, y)
      d0 = sample_distance x, y - 1
      d1 = sample_distance x - 1, y
      d2 = sample_distance x + 1, y
      d3 = sample_distance x, y + 1

      [(d2 - d1) / @scale, (d3 - d0) / @scale]
    end
38561
def gradient(self):
        """"""
        functional = self
        c = self.scale

        class RosenbrockGradient(Operator):

            """The gradient operator of the Rosenbrock functional."""

            def __init__(self):
                """Initialize a new instance."""
                super(RosenbrockGradient, self).__init__(
                    functional.domain, functional.domain, linear=False)

            def _call(self, x, out):
                """Apply the gradient operator to the given point."""
                for i in range(1, self.domain.size - 1):
                    out[i] = (2 * c * (x[i] - x[i - 1]**2) -
                              4 * c * (x[i + 1] - x[i]**2) * x[i] -
                              2 * (1 - x[i]))
                out[0] = (-4 * c * (x[1] - x[0] ** 2) * x[0] +
                          2 * (x[0] - 1))
                out[-1] = 2 * c * (x[-1] - x[-2] ** 2)

            def derivative(self, x):
                """The derivative of the gradient.

                This is also known as the Hessian.
                """
                # TODO: Implement optimized version of this that does not need
                # a matrix.
                shape = (functional.domain.size, functional.domain.size)
                matrix = np.zeros(shape)

                # Straightforward computation
                for i in range(0, self.domain.size - 1):
                    matrix[i, i] = (2 * c + 2 + 12 * c * x[i] ** 2 -
                                    4 * c * x[i + 1])
                    matrix[i + 1, i] = -4 * c * x[i]
                    matrix[i, i + 1] = -4 * c * x[i]
                matrix[-1, -1] = 2 * c
                matrix[0, 0] = 2 + 12 * c * x[0] ** 2 - 4 * c * x[1]
                return MatrixOperator(matrix, self.domain, self.range)

        return RosenbrockGradient()
94993
def sde_update_gradient_full(self, gradients):
        
    
        self.variance.gradient = gradients[0]
        self.lengthscale.gradient = gradients[1]
586101
function(y) { 

      // compute and accumulate gradient wrt weights and bias of this layer
      var x = this.in_act;
      x.dw = global.zeros(x.w.length); // zero out the gradient of input Vol
      var loss = 0.0;
      if(y instanceof Array || y instanceof Float64Array) {
        for(var i=0;i<this.out_depth;i++) {
          var dy = x.w[i] - y[i];
          x.dw[i] = dy;
          loss += 0.5*dy*dy;
        }
      } else if(typeof y === 'number') {
        // lets hope that only one number is being regressed
        var dy = x.w[0] - y;
        x.dw[0] = dy;
        loss += 0.5*dy*dy;
      } else {
        // assume it is a struct with entries .dim and .val
        // and we pass gradient only along dimension dim to be equal to val
        var i = y.dim;
        var yi = y.val;
        var dy = x.w[i] - yi;
        x.dw[i] = dy;
        loss += 0.5*dy*dy;
      }
      return loss;
    }