Compute gradient and Hessian matrix with respect to `w`.
432214
def update(self, w, offset=0):
        """"""
        time = self.time
        x = self.x
        exp_xw = numpy.exp(offset + numpy.dot(x, w))
        n_samples, n_features = x.shape

        gradient = numpy.zeros((1, n_features), dtype=float)
        hessian = numpy.zeros((n_features, n_features), dtype=float)

        inv_n_samples = 1. / n_samples
        risk_set = 0
        risk_set_x = 0
        risk_set_xx = 0
        k = 0
        # iterate time in descending order
        for i in range(n_samples):
            ti = time[i]
            while k < n_samples and ti == time[k]:
                risk_set += exp_xw[k]

                # preserve 2D shape of row vector
                xk = x[k:k + 1]
                risk_set_x += exp_xw[k] * xk

                # outer product
                xx = numpy.dot(xk.T, xk)
                risk_set_xx += exp_xw[k] * xx

                k += 1

            if self.event[i]:
                gradient -= (x[i:i + 1] - risk_set_x / risk_set) * inv_n_samples

                a = risk_set_xx / risk_set
                z = risk_set_x / risk_set
                # outer product
                b = numpy.dot(z.T, z)

                hessian += (a - b) * inv_n_samples

        if self.alpha > 0:
            gradient += self.alpha * inv_n_samples * w

            diag_idx = numpy.diag_indices(n_features)
            hessian[diag_idx] += self.alpha * inv_n_samples

        self.gradient = gradient.ravel()
        self.hessian = hessian
989122
def G(w, Xs):
    ''''''
    n = len(Xs)
    P = projection_matrix(w)
    Ys = [np.dot(P, X) for X in Xs]
    A = calc_A(Ys)
    A_hat = calc_A_hat(A, skew_matrix(w))

    
    u = sum(np.dot(Y, Y) for Y in Ys) / n
    v = np.dot(A_hat, sum(np.dot(Y, Y) * Y for Y in Ys)) / np.trace(np.dot(A_hat, A))

    return sum((np.dot(Y, Y) - u - 2 * np.dot(Y, v)) ** 2 for Y in Ys)
282436
def gamma(self, x, y, kwargs, diff=diff):
        
        f_xx, f_xy, f_yx, f_yy = self.hessian(x, y, kwargs, diff=diff)
        gamma1 = 1./2 * (f_xx - f_yy)
        gamma2 = f_xy
        return gamma1, gamma2
232860
function regularizeNone(weights, gradientCount) {
  const [gradient, count] = gradientCount;

  for (let i = 0; i < gradient.length; i++) {
    let grad = (gradient[i] || 0) / count;
    weights[i] = (weights[i] || 0) - grad;
  }
}
675628
def gradients_X_X2(self, dL_dK, X, X2):
        """"""
        return self._comp_grads(dL_dK, X, X2)[3:]
1087566
def update_gradients_diag(self, dL_dKdiag, X):
        """"""
        self.variance.gradient = np.sum(dL_dKdiag)
        self.period.gradient = 0
        self.lengthscale.gradient = 0
731555
def make_gradients(dims=DEFAULT_DIMS):
    
    return np.meshgrid(
        np.linspace(0.0, 1.0, dims[0]),
        np.linspace(0.0, 1.0, dims[1])
    )
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]
1195563
def fgrad_y_psi(self, y, return_covar_chain=False):
        
        mpsi = self.psi

        w, s, r, d = self.fgrad_y(y, return_precalc=True)
        gradients = np.zeros((y.shape[0], y.shape[1], len(mpsi), 4))
        for i in range(len(mpsi)):
            a,b,c  = mpsi[i]
            gradients[:, :, i, 0] = (b * (1.0/np.cosh(s[i])) ** 2).T
            gradients[:, :, i, 1] = a * (d[i] - 2.0 * s[i] * r[i] * (1.0/np.cosh(s[i])) ** 2).T
            gradients[:, :, i, 2] = (-2.0 * a * (b ** 2) * r[i] * ((1.0 / np.cosh(s[i])) ** 2)).T
        gradients[:, :, 0, 3] = 1.0

        if return_covar_chain:
            covar_grad_chain = np.zeros((y.shape[0], y.shape[1], len(mpsi), 4))
            for i in range(len(mpsi)):
                a,b,c = mpsi[i]
                covar_grad_chain[:, :, i, 0] = (r[i]).T
                covar_grad_chain[:, :, i, 1] = (a * (y + c) * ((1.0 / np.cosh(s[i])) ** 2).T)
                covar_grad_chain[:, :, i, 2] = a * b * ((1.0 / np.cosh(s[i])) ** 2).T
            covar_grad_chain[:, :, 0, 3] = y
            return gradients, covar_grad_chain

        return gradients
978014
private void computeGradients( float kernelRadius, int kernelWidth ) {

        // generate the gaussian convolution masks
        float kernel[] = new float[kernelWidth];
        float diffKernel[] = new float[kernelWidth];
        int kwidth;
        for( kwidth = 0; kwidth < kernelWidth; kwidth++ ) {
            float g1 = gaussian(kwidth, kernelRadius);
            if (g1 <= GAUSSIAN_CUT_OFF && kwidth >= 2)
                break;
            float g2 = gaussian(kwidth - 0.5f, kernelRadius);
            float g3 = gaussian(kwidth + 0.5f, kernelRadius);
            kernel[kwidth] = (g1 + g2 + g3) / 3f
                    / (2f * (float) Math.PI * kernelRadius * kernelRadius);
            diffKernel[kwidth] = g3 - g2;
        }

        int initX = kwidth - 1;
        int maxX = width - (kwidth - 1);
        int initY = width * (kwidth - 1);
        int maxY = width * (height - (kwidth - 1));

        // perform convolution in x and y directions
        for( int x = initX; x < maxX; x++ ) {
            for( int y = initY; y < maxY; y += width ) {
                int index = x + y;
                float sumX = data[index] * kernel[0];
                float sumY = sumX;
                int xOffset = 1;
                int yOffset = width;
                for( ; xOffset < kwidth; ) {
                    sumY += kernel[xOffset] * (data[index - yOffset] + data[index + yOffset]);
                    sumX += kernel[xOffset] * (data[index - xOffset] + data[index + xOffset]);
                    yOffset += width;
                    xOffset++;
                }

                yConv[index] = sumY;
                xConv[index] = sumX;
            }

        }

        for( int x = initX; x < maxX; x++ ) {
            for( int y = initY; y < maxY; y += width ) {
                float sum = 0f;
                int index = x + y;
                for( int i = 1; i < kwidth; i++ )
                    sum += diffKernel[i] * (yConv[index - i] - yConv[index + i]);

                xGradient[index] = sum;
            }

        }

        for( int x = kwidth; x < width - kwidth; x++ ) {
            for( int y = initY; y < maxY; y += width ) {
                float sum = 0.0f;
                int index = x + y;
                int yOffset = width;
                for( int i = 1; i < kwidth; i++ ) {
                    sum += diffKernel[i] * (xConv[index - yOffset] - xConv[index + yOffset]);
                    yOffset += width;
                }

                yGradient[index] = sum;
            }

        }

        initX = kwidth;
        maxX = width - kwidth;
        initY = width * kwidth;
        maxY = width * (height - kwidth);
        for( int x = initX; x < maxX; x++ ) {
            for( int y = initY; y < maxY; y += width ) {
                int index = x + y;
                int indexN = index - width;
                int indexS = index + width;
                int indexW = index - 1;
                int indexE = index + 1;
                int indexNW = indexN - 1;
                int indexNE = indexN + 1;
                int indexSW = indexS - 1;
                int indexSE = indexS + 1;

                float xGrad = xGradient[index];
                float yGrad = yGradient[index];
                float gradMag = hypot(xGrad, yGrad);

                // perform non-maximal supression
                float nMag = hypot(xGradient[indexN], yGradient[indexN]);
                float sMag = hypot(xGradient[indexS], yGradient[indexS]);
                float wMag = hypot(xGradient[indexW], yGradient[indexW]);
                float eMag = hypot(xGradient[indexE], yGradient[indexE]);
                float neMag = hypot(xGradient[indexNE], yGradient[indexNE]);
                float seMag = hypot(xGradient[indexSE], yGradient[indexSE]);
                float swMag = hypot(xGradient[indexSW], yGradient[indexSW]);
                float nwMag = hypot(xGradient[indexNW], yGradient[indexNW]);
                float tmp;
                /*
                 * An explanation of what's happening here, for those who want
                 * to understand the source: This performs the "non-maximal
                 * supression" phase of the Canny edge detection in which we
                 * need to compare the gradient magnitude to that in the
                 * direction of the gradient; only if the value is a local
                 * maximum do we consider the point as an edge candidate.
                 * 
                 * We need to break the comparison into a number of different
                 * cases depending on the gradient direction so that the
                 * appropriate values can be used. To avoid computing the
                 * gradient direction, we use two simple comparisons: first we
                 * check that the partial derivatives have the same sign (1)
                 * and then we check which is larger (2). As a consequence, we
                 * have reduced the problem to one of four identical cases that
                 * each test the central gradient magnitude against the values at
                 * two points with 'identical support'; what this means is that
                 * the geometry required to accurately interpolate the magnitude
                 * of gradient function at those points has an identical
                 * geometry (upto right-angled-rotation/reflection).
                 * 
                 * When comparing the central gradient to the two interpolated
                 * values, we avoid performing any divisions by multiplying both
                 * sides of each inequality by the greater of the two partial
                 * derivatives. The common comparand is stored in a temporary
                 * variable (3) and reused in the mirror case (4).
                 * 
                 */
                if (xGrad * yGrad <= (float) 0 /*(1)*/
                ? Math.abs(xGrad) >= Math.abs(yGrad) /*(2)*/
                ? (tmp = Math.abs(xGrad * gradMag)) >= Math.abs(yGrad * neMag - (xGrad + yGrad)
                        * eMag) /*(3)*/
                        && tmp > Math.abs(yGrad * swMag - (xGrad + yGrad) * wMag) /*(4)*/
                : (tmp = Math.abs(yGrad * gradMag)) >= Math.abs(xGrad * neMag - (yGrad + xGrad)
                        * nMag) /*(3)*/
                        && tmp > Math.abs(xGrad * swMag - (yGrad + xGrad) * sMag) /*(4)*/
                : Math.abs(xGrad) >= Math.abs(yGrad) /*(2)*/
                ? (tmp = Math.abs(xGrad * gradMag)) >= Math.abs(yGrad * seMag + (xGrad - yGrad)
                        * eMag) /*(3)*/
                        && tmp > Math.abs(yGrad * nwMag + (xGrad - yGrad) * wMag) /*(4)*/
                : (tmp = Math.abs(yGrad * gradMag)) >= Math.abs(xGrad * seMag + (yGrad - xGrad)
                        * sMag) /*(3)*/
                        && tmp > Math.abs(xGrad * nwMag + (yGrad - xGrad) * nMag) /*(4)*/
                ) {
                    magnitude[index] = gradMag >= MAGNITUDE_LIMIT
                            ? MAGNITUDE_MAX
                            : (int) (MAGNITUDE_SCALE * gradMag);
                    // NOTE: The orientation of the edge is not employed by this
                    // implementation. It is a simple matter to compute it at
                    // this point as: Math.atan2(yGrad, xGrad);
                } else {
                    magnitude[index] = 0;
                }
            }
        }
    }
675501
private void computeGradientAndHessian(DMatrixRMaj param  )
    {
        // residuals = f(x) - y
        function.compute(param, residuals);

        computeNumericalJacobian(param,jacobian);

        CommonOps_DDRM.multTransA(jacobian, residuals, g);
        CommonOps_DDRM.multTransA(jacobian, jacobian,  H);

        CommonOps_DDRM.extractDiag(H,Hdiag);
    }
463897
def gradF_t(X, Y, S, M_E, E, m0, rho):
    ''' 
    '''
    n, r = X.shape
    m, r = Y.shape

    XS = np.dot(X, S)
    YS = np.dot(Y, S.T)
    XSY = np.dot(XS, Y.T)

    Qx = np.dot(np.dot(X.T, ((M_E - XSY) * E)), YS) / n
    Qy = np.dot(np.dot(Y.T, ((M_E - XSY) * E).T), XS) / m

    W = np.dot((XSY - M_E) * E, YS) + np.dot(X, Qx) + rho * Gp(X, m0, r)
    Z = np.dot(((XSY - M_E) * E).T, XS) + np.dot(Y, Qy) + rho * Gp(Y, m0, r)

    return W, Z
443947
def C(w, Xs):
    '''
    '''
    n = len(Xs)
    P = projection_matrix(w)
    Ys = [np.dot(P, X) for X in Xs]
    A = calc_A(Ys)
    A_hat = calc_A_hat(A, skew_matrix(w))

    return np.dot(A_hat, sum(np.dot(Y, Y) * Y for Y in Ys)) / np.trace(np.dot(A_hat, A))
179950
def hessfunc(self, p):
        
        self._set_stochastics(p)
        for i in xrange(self.len):

            di = self.diff(i)
            self.hess[i, i] = self.diff(i, 2)

            if i < self.len - 1:

                for j in xrange(i + 1, self.len):
                    dij = self.diff2(i, j)

                    self.hess[i, j] = dij
                    self.hess[j, i] = dij

        return -1. * self.hess
1118344
def _get_numdiff_hessian_element_func(objective_func, nmr_steps, step_ratio):
    """"""
    return SimpleCLFunction.from_string('''
        /**
         * Compute the Hessian using (possibly) multiple steps with various interpolations. 
         */ 
        void _numdiff_hessian_element(
                void* data, local mot_float_type* x_tmp, mot_float_type f_x_input,
                uint px, uint py, global float* initial_step, global double* derivative, 
                global double* error, local double* scratch){

            const uint nmr_steps = ''' + str(nmr_steps) + ''';
            uint nmr_steps_remaining = nmr_steps;

            local double* scratch_ind = scratch;
            local double* steps = scratch_ind;      scratch_ind += nmr_steps; 
            local double* errors = scratch_ind;     scratch_ind += nmr_steps - 1;
            local double* steps_tmp = scratch_ind;  scratch_ind += nmr_steps;
            
            if(get_local_id(0) == 0){
                for(int i = 0; i < nmr_steps - 1; i++){
                    errors[i] = 0;
                }
            }
            barrier(CLK_LOCAL_MEM_FENCE);
            
            _numdiff_hessian_steps(data, x_tmp, f_x_input, px, py, steps, initial_step);

            if(nmr_steps_remaining > 1){
                nmr_steps_remaining = _numdiff_hessian_richardson_extrapolation(steps); 
                barrier(CLK_LOCAL_MEM_FENCE);
            }
            
            if(nmr_steps_remaining >= 3){
                nmr_steps_remaining = _numdiff_wynn_extrapolation(steps, errors, nmr_steps_remaining);
                barrier(CLK_LOCAL_MEM_FENCE);                
            }
            
            if(nmr_steps_remaining > 1){
                _numdiff_find_best_step(steps, errors, steps_tmp, nmr_steps_remaining);
                barrier(CLK_LOCAL_MEM_FENCE);
            }
            
            if(get_local_id(0) == 0){
                *derivative = steps[0];
                *error = errors[0];
            }
        }
    ''', dependencies=[
        _get_numdiff_hessian_steps_func(objective_func, nmr_steps, step_ratio),
        _get_numdiff_hessian_richardson_extrapolation_func(nmr_steps, step_ratio),
        _get_numdiff_wynn_extrapolation_func(),
        _get_numdiff_find_best_step_func()
    ])
319001
def fit(self, n_iter, eta=0.01):
        """"""
        for _ in range(n_iter):
            gradient = self.compute_gradient()
            self.gradient_step(gradient, eta)
1263071
def _make_linear_gradient(width, height, angle=90.):
    """"""
    import numpy as np
    X, Y = np.meshgrid(np.linspace(0, 1, width), np.linspace(0, 1, height))
    theta = np.radians(angle % 360)
    c, s = np.cos(theta), np.sin(theta)
    if 0 <= theta and theta < 0.5 * np.pi:
        Z = np.abs(c * X + s * Y)
    elif 0.5 * np.pi <= theta and theta < np.pi:
        Z = np.abs(c * (X - width) + s * Y)
    elif np.pi <= theta and theta < 1.5 * np.pi:
        Z = np.abs(c * (X - width) + s * (Y - height))
    elif 1.5 * np.pi <= theta and theta < 2.0 * np.pi:
        Z = np.abs(c * X + s * (Y - height))
    return (Z - Z.min()) / (Z.max() - Z.min())