Evaluates the cost Hessian.
174072
def _d2f(self, x):
        
        d2f_dPg2 = lil_matrix((self._ng, 1)) # w.r.t p.u. Pg
        d2f_dQg2 = lil_matrix((self._ng, 1)) # w.r.t p.u. Qg]

        for i in self._ipol:
            p_cost = list(self._gn[i].p_cost)
            d2f_dPg2[i, 0] = polyval(polyder(p_cost, 2),
                self._Pg.v0[i] * self._base_mva) * self._base_mva**2
#            for i in ipol:
#                d2f_dQg2[i] = polyval(polyder(list(gn[i].p_cost), 2),
#                                      Qg.v0[i] * base_mva) * base_mva**2

        i = r_[range(self._Pg.i1, self._Pg.iN + 1),
               range(self._Qg.i1, self._Qg.iN + 1)]

        d2f = csr_matrix((vstack([d2f_dPg2, d2f_dQg2]).toarray().flatten(),
                          (i, i)), shape=(self._nxyz, self._nxyz))
        return d2f
1023880
def build_train_op_from_cost(self, cost):
    """"""
    self.lrn_rate = tf.constant(self.hps.lrn_rate, tf.float32,
                                name='learning_rate')
    self.momentum = tf.constant(self.hps.momentum, tf.float32,
                                name='momentum')

    trainable_variables = tf.trainable_variables()
    grads = tf.gradients(cost, trainable_variables)
    devs = {v.device for v in trainable_variables}
    assert len(devs) == 1, ('There should be no trainable variables'
                            ' on any device other than the last GPU.')

    optimizer = tf.train.MomentumOptimizer(self.lrn_rate, self.momentum)

    gv_pairs = zip(grads, trainable_variables)
    gv_pairs = [gv for gv in gv_pairs if gv[0] is not None]
    devs = {gv[1].device for gv in gv_pairs}
    assert len(devs) == 1, ('There should be no gradients wrt'
                            ' vars on other GPUs.')

    apply_op = optimizer.apply_gradients(
        gv_pairs,
        global_step=self.global_step, name='train_step')

    train_ops = [apply_op]
    train_op = tf.group(*train_ops)
    return train_op
445598
public static <T extends ImageGray<T>>
	void hessian( T integral, int skip , int size ,
				  GrayF32 intensity) {

		if( integral instanceof GrayF32) {
			IntegralImageFeatureIntensity.hessian((GrayF32)integral,skip,size,intensity);
		} else if( integral instanceof GrayS32) {
			IntegralImageFeatureIntensity.hessian((GrayS32)integral,skip,size,intensity);
		} else {
			throw new IllegalArgumentException("Unsupported input type");
		}
	}
517990
def _cost_method(self, *args, **kwargs):
        

        cost_val = np.sum(np.abs(self.weights * self._linear.op(args[0])))

        if 'verbose' in kwargs and kwargs['verbose']:
            print(' - L1 NORM (X):', cost_val)

        return cost_val
208523
private double optimize(int k0, int kmax, double sumx, double sumx2, double xp, double yp, double sumxy, double sumy) {
    int k = kmax - k0 + 1;
    return (sumxy - xp * sumy - yp * sumx + k * xp * yp) / (sumx2 - 2 * sumx * xp + k * xp * xp);
    // return (-xp * yp * + yp * sumx - sumxy + xp * sumy) / (-xp * xp *
    // kmax - sumx2 + 2 * xp * sumx);
  }
1058449
public static void hess(Matrix A, ExecutorService threadpool)
    {
        if(!A.isSquare())
            throw new ArithmeticException("Only square matrices can be converted to Upper Hessenberg form");
        int m = A.rows();
        /**
         * Space used to store the vector for updating the columns of A
         */
        DenseVector columnUpdateTmp = new DenseVector(m);
        double[] vk = new double[m];
        /**
         * Space used for updating the sub matrix at step i
         */
        double[] subMatrixUpdateTmp = new double[m];
        double tmp;//Used for temp values
        for(int i = 0; i < m-2; i++)
        {
            //Holds the norm, sqrt{a_i^2 + ... + a_m^2}
            double s = 0.0;
            //First step of the loop done outside to do extra bit
            double sigh = A.get(i+1, i);//Holds the multiplication factor
            vk[i+1] = sigh;
            s += sigh*sigh;
            sigh = sigh > 0 ? 1 : -1;//Sign dosnt change the squaring, so we do it first
            
            
            
            for(int j = i+2; j < m; j++)
            {
                tmp = A.get(j, i);
                vk[j] = tmp;
                s += tmp*tmp;
            }
            
            double s1 = -sigh*Math.sqrt(s);
            //Now re use s to quickly get the norm of vk, since it will be almost the same vector
            s -= vk[i+1]*vk[i+1];
            vk[i+1] -= s1;
            s += vk[i+1]*vk[i+1];
            double s1Inv = 1.0/Math.sqrt(s);//Re use to store the norm of vk. Do the inverse to multiply quickly instead of divide
            for(int j = i+1; j < m; j++)
                vk[j] *= s1Inv;
            
            //Update sub sub matrix A[i+1:m, i:m]
            //NOTE: The first column that will be altered can be done ourslves, since we know the value set (s1) and that all below it will ber zero
            Matrix subA = new SubMatrix(A, i+1, i, m, m);
            DenseVector vVec = new DenseVector(vk, i+1, m);
            Vec tmpV = new DenseVector(subMatrixUpdateTmp, i, m);
            tmpV.zeroOut();
            vVec.multiply(subA, tmpV);
            
            if(threadpool == null)
                OuterProductUpdate(subA, vVec, tmpV, -2.0);
            else
                OuterProductUpdate(subA, vVec, tmpV, -2.0, threadpool);
            //Zero out ourselves after.
            //TODO implement so we dont compute the first row
            A.set(i+1, i, s1);
            for(int j = i+2; j < m; j++)
                A.set(j, i, 0.0);
            
            
            //Update the columns of A[0:m, i+1:m] 
            subA = new SubMatrix(A, 0, i+1, m, m);
            columnUpdateTmp.zeroOut();
            subA.multiply(vVec, 1.0, columnUpdateTmp);
            if(threadpool == null)
                OuterProductUpdate(subA, columnUpdateTmp, vVec, -2.0);
            else
                OuterProductUpdate(subA, columnUpdateTmp, vVec, -2.0, threadpool);
        }
    }
434713
def optimize(self, x0, target):
        """"""
        x = x0
        for i in range(self.maxiter):
            g = self.g(x, target)
            h = self.h(x, target)
            if i == 0:
                alpha = 0
                m = g
            else:
                alpha = - np.dot(m, np.dot(h, g)) / np.dot(m, np.dot(h, m))
                m = g + np.dot(alpha, m)
            t = - np.dot(m, g) / np.dot(m, np.dot(h, m))
            delta = np.dot(t, m)
            x = x + delta
            if np.linalg.norm(delta) < self.tol:
                break
        return x
69986
def build_cost(self, labels, logits):
    
    op = logits.op
    if "softmax" in str(op).lower():
      logits, = op.inputs

    with tf.variable_scope('costs'):
      xent = tf.nn.softmax_cross_entropy_with_logits(
          logits=logits, labels=labels)
      cost = tf.reduce_mean(xent, name='xent')
      cost += self._decay()
      cost = cost

    return cost
542790
def ellipticity(self):
        '''
        
        '''
        eig = np.linalg.eig(self.field_hessian)
        eig.sort()
        return eig[0]/eig[1] - 1
769600
def hessian_component(self, index1, index2):
        """"""
        result = np.zeros((3, 3), float)
        if index1 == index2:
            for index3 in range(self.numc):
                if self.scaling[index1, index3] > 0:
                    d_1 = 1/self.distances[index1, index3]
                    for (se, ve), (sg, vg), (sh, vh) in zip(
                        self.yield_pair_energies(index1, index3),
                        self.yield_pair_gradients(index1, index3),
                        self.yield_pair_hessians(index1, index3)
                    ):
                        result += (
                            +sh*self.dirouters[index1, index3]*ve
                            +sg*(np.identity(3, float) - self.dirouters[index1, index3])*ve*d_1
                            +sg*np.outer(self.directions[index1, index3],  vg)
                            +sg*np.outer(vg, self.directions[index1, index3])
                            +se*vh
                        )*self.scaling[index1, index3]
        elif self.scaling[index1, index2] > 0:
            d_1 = 1/self.distances[index1, index2]
            for (se, ve), (sg, vg), (sh, vh) in zip(
                self.yield_pair_energies(index1, index2),
                self.yield_pair_gradients(index1, index2),
                self.yield_pair_hessians(index1, index2)
            ):
                result -= (
                    +sh*self.dirouters[index1, index2]*ve
                    +sg*(np.identity(3, float) - self.dirouters[index1, index2])*ve*d_1
                    +sg*np.outer(self.directions[index1, index2],  vg)
                    +sg*np.outer(vg, self.directions[index1, index2])
                    +se*vh
                )*self.scaling[index1, index2]
        return result
686642
def update(self, gradient, step):
        """"""
        self.old_gradient = self.gradient
        self.gradient = gradient
        N = len(self.gradient)
        if self.inv_hessian is None:
            # update the direction
            self.direction = -self.gradient
            self.status = "SD"
            # new guess of the inverse hessian
            self.inv_hessian = np.identity(N, float)
        else:
            # update the direction
            self.direction = -np.dot(self.inv_hessian, self.gradient)
            self.status = "QN"
            # new guess of the inverse hessian (BFGS)
            y = self.gradient - self.old_gradient
            s = step
            sy = abs(np.dot(s, y))+1e-5
            A = np.outer(-y/sy, s)
            A.ravel()[::N+1] += 1
            self.inv_hessian = (
                np.dot(np.dot(A.transpose(), self.inv_hessian), A) +
                np.outer(s/sy, s)
            )
585702
def hessian(self, x, y, coeffs, beta, center_x=0, center_y=0):
        
        shapelets = self._createShapelet(coeffs)
        r, phi = param_util.cart2polar(x, y, center=np.array([center_x, center_y]))
        kappa_shapelets=self._kappaShapelets(shapelets, beta)
        gamma1_shapelets, gamma2_shapelets=self._gammaShapelets(shapelets, beta)
        kappa_value=self._shapeletOutput(r, phi, beta, kappa_shapelets)
        gamma1_value=self._shapeletOutput(r, phi, beta, gamma1_shapelets)
        gamma2_value=self._shapeletOutput(r, phi, beta, gamma2_shapelets)
        f_xx = kappa_value + gamma1_value
        f_xy = gamma2_value
        f_yy = kappa_value - gamma1_value
        return f_xx, f_yy, f_xy
549109
def _cost_method(self, *args, **kwargs):
        

        cost_val = self.thresh * nuclear_norm(cube2matrix(args[0]))

        if 'verbose' in kwargs and kwargs['verbose']:
            print(' - NUCLEAR NORM (X):', cost_val)

        return cost_val
867246
def get_hessian(self):
        """"""
        force_const = self.fields.get("Cartesian Force Constants")
        if force_const is None:
            return None
        N = len(self.molecule.numbers)
        result = np.zeros((3*N, 3*N), float)
        counter = 0
        for row in range(3*N):
            result[row, :row+1] = force_const[counter:counter+row+1]
            result[:row+1, row] = force_const[counter:counter+row+1]
            counter += row + 1
        return result
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
959730
def hessian(self, x, y, kwargs, diff=diff):
        
        alpha_ra, alpha_dec = self.alpha(x, y, kwargs)

        alpha_ra_dx, alpha_dec_dx = self.alpha(x + diff, y, kwargs)
        alpha_ra_dy, alpha_dec_dy = self.alpha(x, y + diff, kwargs)

        dalpha_rara = (alpha_ra_dx - alpha_ra)/diff
        dalpha_radec = (alpha_ra_dy - alpha_ra)/diff
        dalpha_decra = (alpha_dec_dx - alpha_dec)/diff
        dalpha_decdec = (alpha_dec_dy - alpha_dec)/diff

        f_xx = dalpha_rara
        f_yy = dalpha_decdec
        f_xy = dalpha_radec
        f_yx = dalpha_decra
        return f_xx, f_xy, f_yx, f_yy
376163
def yield_pair_hessians(self, index1, index2):
        """"""
        d_1 = 1/self.distances[index1, index2]
        d_3 = d_1**3
        if self.charges is not None:
            c1 = self.charges[index1]
            c2 = self.charges[index2]
            yield 2*c1*c2*d_3, np.zeros((3, 3))
        if self.dipoles is not None:
            d_5 = d_1**5
            d_7 = d_1**7
            p1 = self.dipoles[index1]
            p2 = self.dipoles[index2]
            yield 12*d_5*np.dot(p1, p2), np.zeros((3, 3))
            yield -90*d_7, np.outer(p1, p2) + np.outer(p2, p1)
            if self.charges is not None:
                yield 12*c1*d_5, np.zeros((3, 3))
                yield 12*c2*d_5, np.zeros((3, 3))