Get a function to compute the multiple step sizes for a single element of the Hessian.
979474
def _get_numdiff_hessian_steps_func(objective_func, nmr_steps, step_ratio):
    """"""
    return SimpleCLFunction.from_string('''
        /**
         * Compute one element of the Hessian for a number of steps.
         * 
         * This uses the initial steps in the data structure, indexed by the parameters to change (px, py).
         *
         * Args:
         *  data: the data container
         *  x_tmp: the array with the input parameters, needs to be writable, although it will return
         *         the same values.
         *  f_x_input: the objective function value at the original set of parameters  
         *  px: the index of the first parameter to perturbate
         *  py: the index of the second parameter to perturbate
         *  steps: storage location for the output steps
         *  initial_step: the initial steps, array of same length as x_temp
         */
        void _numdiff_hessian_steps(void* data, local mot_float_type* x_tmp, 
                                    mot_float_type f_x_input,
                                    uint px, uint py, 
                                    local double* steps, 
                                    global float* initial_step){

            double step_x;
            double step_y;
            double tmp;
            bool is_first_workitem = get_local_id(0) == 0;

            if(px == py){
                for(uint step_ind = 0; step_ind < ''' + str(nmr_steps) + '''; step_ind++){
                    step_x = initial_step[px] / pown(''' + str(float(step_ratio)) + ''', step_ind);

                    tmp = (
                          _numdiff_hessian_eval_step_mono(data, x_tmp, px, 2 * step_x)
                        + _numdiff_hessian_eval_step_mono(data, x_tmp, px, -2 * step_x)
                        - 2 * f_x_input
                    ) / (4 * step_x * step_x);    

                    if(is_first_workitem){
                        steps[step_ind] = tmp;
                    }
                }
            }
            else{
                for(uint step_ind = 0; step_ind < ''' + str(nmr_steps) + '''; step_ind++){
                    step_x = initial_step[px] / pown(''' + str(float(step_ratio)) + ''', step_ind);
                    step_y = initial_step[py] / pown(''' + str(float(step_ratio)) + ''', step_ind);

                    tmp = (
                          _numdiff_hessian_eval_step_bi(data, x_tmp, px, step_x, py, step_y)
                        - _numdiff_hessian_eval_step_bi(data, x_tmp, px, step_x, py, -step_y)
                        - _numdiff_hessian_eval_step_bi(data, x_tmp, px, -step_x, py, step_y)
                        + _numdiff_hessian_eval_step_bi(data, x_tmp, px, -step_x, py, -step_y)
                    ) / (4 * step_x * step_y);

                    if(is_first_workitem){
                        steps[step_ind] = tmp;
                    }                       
                }
            }
        }
    ''', dependencies=[SimpleCLFunction.from_string('''
        /**
         * Evaluate the model with a perturbation in one dimensions.
         *
         * Args:
         *  data: the data container
         *  x_tmp: the array with the input parameters, needs to be writable, although it will return
         *         the same values.
         *  perturb_dim0: the index (into the x_tmp parameters) of the parameter to perturbate
         *  perturb_0: the added perturbation of the index corresponding to ``perturb_dim_0``
         *
         * Returns:
         *  the function evaluated at the parameters plus their perturbation.
         */
        double _numdiff_hessian_eval_step_mono(
                void* data, local mot_float_type* x_tmp, 
                uint perturb_dim_0, mot_float_type perturb_0){

            mot_float_type old_0;
            double return_val;

            if(get_local_id(0) == 0){
                old_0 = x_tmp[perturb_dim_0];
                x_tmp[perturb_dim_0] += perturb_0;
            }
            barrier(CLK_LOCAL_MEM_FENCE);

            return_val = ''' + objective_func.get_cl_function_name() + '''(x_tmp, data);
            barrier(CLK_LOCAL_MEM_FENCE);

            if(get_local_id(0) == 0){
                x_tmp[perturb_dim_0] = old_0;
            }
            barrier(CLK_LOCAL_MEM_FENCE);

            return return_val;
        }
    '''), SimpleCLFunction.from_string('''
        /**
         * Evaluate the model with a perturbation in two dimensions.
         *
         * Args:
         *  data: the data container
         *  x_tmp: the array with the input parameters, needs to be writable, although it will return
         *         the same values.
         *  perturb_dim_0: the index (into the x_tmp parameters) of the first parameter to perturbate
         *  perturb_0: the added perturbation of the index corresponding to ``perturb_dim_0``
         *  perturb_dim_1: the index (into the x_tmp parameters) of the second parameter to perturbate
         *  perturb_1: the added perturbation of the index corresponding to ``perturb_dim_1``
         *
         * Returns:
         *  the function evaluated at the parameters plus their perturbation.
         */
        double _numdiff_hessian_eval_step_bi(
                void* data, local mot_float_type* x_tmp, 
                uint perturb_dim_0, mot_float_type perturb_0,
                uint perturb_dim_1, mot_float_type perturb_1){

            mot_float_type old_0;
            mot_float_type old_1;
            double return_val;

            if(get_local_id(0) == 0){
                old_0 = x_tmp[perturb_dim_0];
                old_1 = x_tmp[perturb_dim_1];

                x_tmp[perturb_dim_0] += perturb_0;
                x_tmp[perturb_dim_1] += perturb_1;
            }
            barrier(CLK_LOCAL_MEM_FENCE);

            return_val = ''' + objective_func.get_cl_function_name() + '''(x_tmp, data);
            barrier(CLK_LOCAL_MEM_FENCE);

            if(get_local_id(0) == 0){
                x_tmp[perturb_dim_0] = old_0;
                x_tmp[perturb_dim_1] = old_1;
            }
            barrier(CLK_LOCAL_MEM_FENCE);

            return return_val;
        }
    ''')])
662532
def smallStepsPerLargeStep(self, smallStepsPerLargeStep):
        
        self._smallStepsPerLargeStep = smallStepsPerLargeStep
        self._smallStepFactor = np.power(self.largeStepFactor, 1.0 / smallStepsPerLargeStep)
417071
def hessian(self, x, y, amp, sigma_x, sigma_y, center_x = 0, center_y = 0):
        
        f_ = self.function(x, y, amp, sigma_x, sigma_y, center_x, center_y)
        f_xx = f_ * ( (-1./sigma_x**2) + (center_x-x)**2/sigma_x**4 )
        f_yy = f_ * ( (-1./sigma_y**2) + (center_y-y)**2/sigma_y**4 )
        f_xy = f_ * (center_x-x)/sigma_x**2 * (center_y-y)/sigma_y**2
        return f_xx, f_yy, f_xy
1330118
protected DoubleMatrix2D getHessF0(DoubleMatrix1D X) {
		double[][] hess = request.getF0().hessian(X.toArray());
		if(hess == FunctionsUtils.ZEROES_2D_ARRAY_PLACEHOLDER){
			return F2.make(X.size(), X.size());
		}else{
			return F2.make(hess);
		}
	}
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]
32165
def xstep(self):
        """

        self.YU[:] = self.Y - self.U
        YUf = sl.rfftn(self.YU, None, self.cri.axisN)
        YUf0 = self.block_sep0(YUf)
        YUf1 = self.block_sep1(YUf)

        b = self.rho * np.sum(np.conj(self.GDf) * YUf1, axis=-1)
        if self.cri.Cd > 1:
            b = np.sum(b, axis=self.cri.axisC, keepdims=True)
        b += self.DSf + self.rho*YUf0

        # Concatenate multiple GDf components on axisC. For
        # single-channel signals, and multi-channel signals with a
        # single-channel dictionary, we end up with sl.solvemdbi_ism
        # solving a linear system of rank dimN+1 (corresponding to the
        # dictionary and a gradient operator per spatial dimension) plus
        # an identity. For multi-channel signals with a multi-channel
        # dictionary, we end up with sl.solvemdbi_ism solving a linear
        # system of rank C.d (dimN+1) (corresponding to the dictionary
        # and a gradient operator per spatial dimension for each
        # channel) plus an identity.

        # The structure of the linear system to be solved depends on the
        # number of channels in the signal and dictionary. Both branches are
        # the same in the single-channel signal case (the choice of handling
        # it via the 'else' branch is somewhat arbitrary).
        if self.cri.C > 1 and self.cri.Cd == 1:
            # Concatenate multiple GDf components on the final axis
            # of GDf (that indexes the number of gradient operators). For
            # multi-channel signals with a single-channel dictionary,
            # sl.solvemdbi_ism has to solve a linear system of rank dimN+1
            # (corresponding to the dictionary and a gradient operator per
            # spatial dimension)
            DfGDf = np.concatenate(
                [self.Df[..., np.newaxis],] +
                [np.sqrt(self.rho)*self.GDf[..., k, np.newaxis] for k
                 in range(self.GDf.shape[-1])], axis=-1)
            self.Xf[:] = sl.solvemdbi_ism(DfGDf, self.rho, b[..., np.newaxis],
                                          self.cri.axisM, -1)[..., 0]
        else:
            # Concatenate multiple GDf components on axisC. For multi-channel
            # signals with a multi-channel dictionary, sl.solvemdbi_ism has
            # to solve a linear system of rank C.d (dimN+1) (corresponding to
            # the dictionary and a gradient operator per spatial dimension
            # for each channel) plus an identity.
            DfGDf = np.concatenate(
                [self.Df,] + [np.sqrt(self.rho)*self.GDf[..., k] for k
                              in range(self.GDf.shape[-1])],
                axis=self.cri.axisC)
            self.Xf[:] = sl.solvemdbi_ism(DfGDf, self.rho, b, self.cri.axisM,
                                          self.cri.axisC)

        self.X = sl.irfftn(self.Xf, self.cri.Nv, self.cri.axisN)

        if self.opt['LinSolveCheck']:
            if self.cri.C > 1 and self.cri.Cd == 1:
                Dop = lambda x: sl.inner(DfGDf, x[..., np.newaxis],
                                         axis=self.cri.axisM)
                DHop = lambda x: sl.inner(np.conj(DfGDf), x, axis=-1)
                ax = DHop(Dop(self.Xf))[..., 0] + self.rho*self.Xf
            else:
                Dop = lambda x: sl.inner(DfGDf, x, axis=self.cri.axisM)
                DHop = lambda x: sl.inner(np.conj(DfGDf), x,
                                          axis=self.cri.axisC)
                ax = DHop(Dop(self.Xf)) + self.rho*self.Xf
            self.xrrs = sl.rrs(ax, b)
        else:
            self.xrrs = None
146261
def checkgrad_block(self, analytic_hess, numeric_hess, verbose=False, step=1e-6, tolerance=1e-3, block_indices=None, plot=False):
        
        if analytic_hess.dtype is np.dtype('object'):
            #Make numeric hessian also into a block matrix
            real_size = get_block_shapes(analytic_hess)
            num_elements = np.sum(real_size)
            if (num_elements, num_elements) == numeric_hess.shape:
                #If the sizes are the same we assume they are the same
                #(we have not fixed any values so the numeric is the whole hessian)
                numeric_hess = get_blocks(numeric_hess, real_size)
            else:
                #Make a fake empty matrix and fill out the correct block
                tmp_numeric_hess = get_blocks(np.zeros((num_elements, num_elements)), real_size)
                tmp_numeric_hess[block_indices] = numeric_hess.copy()
                numeric_hess = tmp_numeric_hess

        if block_indices is not None:
            #Extract the right block
            analytic_hess = analytic_hess[block_indices]
            numeric_hess = numeric_hess[block_indices]
        else:
            #Unblock them if they are in blocks and you aren't checking a single block (checking whole hessian)
            if analytic_hess.dtype is np.dtype('object'):
                analytic_hess = unblock(analytic_hess)
                numeric_hess = unblock(numeric_hess)

        ratio = numeric_hess / (numpy.where(analytic_hess==0, 1e-10, analytic_hess))
        difference = numpy.abs(analytic_hess - numeric_hess)

        check_passed = numpy.all((numpy.abs(1 - ratio)) < tolerance) or numpy.allclose(numeric_hess, analytic_hess, atol = tolerance)

        if verbose:
            if block_indices:
                print("\nBlock {}".format(block_indices))
            else:
                print("\nAll blocks")

            header = ['Checked', 'Max-Ratio', 'Min-Ratio', 'Min-Difference', 'Max-Difference']
            header_string = map(lambda x: ' | '.join(header), [header])
            separator = '-' * len(header_string[0])
            print('\n'.join([header_string[0], separator]))
            min_r = '%.6f' % float(numpy.min(ratio))
            max_r = '%.6f' % float(numpy.max(ratio))
            max_d = '%.6f' % float(numpy.max(difference))
            min_d = '%.6f' % float(numpy.min(difference))
            cols = [max_r, min_r, min_d, max_d]

            if check_passed:
                checked = "\033[92m  True \033[0m"
            else:
                checked = "\033[91m  False \033[0m"

            grad_string = "{} | {}  | {} |    {}    |   {} ".format(checked, cols[0], cols[1], cols[2], cols[3])
            print(grad_string)

            if plot:
                from matplotlib import pyplot as pb
                fig, axes = pb.subplots(2, 2)
                max_lim = numpy.max(numpy.vstack((analytic_hess, numeric_hess)))
                min_lim = numpy.min(numpy.vstack((analytic_hess, numeric_hess)))
                msa = axes[0,0].matshow(analytic_hess, vmin=min_lim, vmax=max_lim)
                axes[0,0].set_title('Analytic hessian')
                axes[0,0].xaxis.set_ticklabels([None])
                axes[0,0].yaxis.set_ticklabels([None])
                axes[0,0].xaxis.set_ticks([None])
                axes[0,0].yaxis.set_ticks([None])
                msn = axes[0,1].matshow(numeric_hess, vmin=min_lim, vmax=max_lim)
                pb.colorbar(msn, ax=axes[0,1])
                axes[0,1].set_title('Numeric hessian')
                axes[0,1].xaxis.set_ticklabels([None])
                axes[0,1].yaxis.set_ticklabels([None])
                axes[0,1].xaxis.set_ticks([None])
                axes[0,1].yaxis.set_ticks([None])
                msr = axes[1,0].matshow(ratio)
                pb.colorbar(msr, ax=axes[1,0])
                axes[1,0].set_title('Ratio')
                axes[1,0].xaxis.set_ticklabels([None])
                axes[1,0].yaxis.set_ticklabels([None])
                axes[1,0].xaxis.set_ticks([None])
                axes[1,0].yaxis.set_ticks([None])
                msd = axes[1,1].matshow(difference)
                pb.colorbar(msd, ax=axes[1,1])
                axes[1,1].set_title('difference')
                axes[1,1].xaxis.set_ticklabels([None])
                axes[1,1].yaxis.set_ticklabels([None])
                axes[1,1].xaxis.set_ticks([None])
                axes[1,1].yaxis.set_ticks([None])
                if block_indices:
                    fig.suptitle("Block: {}".format(block_indices))
                pb.show()

        return check_passed
26943
def nominal_step(x=None):
    """"""
    if x is None:
        return 1.0
    return np.log1p(np.abs(x)).clip(min=1.0)
804933
public static void hessianNaive(GrayF32 integral, int skip , int size ,
									GrayF32 intensity)
	{
		final int w = intensity.width;
		final int h = intensity.height;

		// get convolution kernels for the second order derivatives
		IntegralKernel kerXX = DerivativeIntegralImage.kernelDerivXX(size,null);
		IntegralKernel kerYY = DerivativeIntegralImage.kernelDerivYY(size,null);
		IntegralKernel kerXY = DerivativeIntegralImage.kernelDerivXY(size,null);

		float norm = 1.0f/(size*size);

		for( int y = 0; y < h; y++ ) {
			for( int x = 0; x < w; x++ ) {

				int xx = x*skip;
				int yy = y*skip;

				computeHessian(integral, intensity, kerXX, kerYY, kerXY, norm, y, yy, x, xx);
			}
		}
	}
123902
def hessian(self, x, y, coeffs, beta, center_x=0, center_y=0):
        
        shapelets = self._createShapelet(coeffs)
        n_order = self._get_num_n(len(coeffs))
        dxx_shapelets = self._dxx_shapelets(shapelets, beta)
        dyy_shapelets = self._dyy_shapelets(shapelets, beta)
        dxy_shapelets = self._dxy_shapelets(shapelets, beta)
        n = len(np.atleast_1d(x))
        if n <= 1:
            f_xx = self._shapeletOutput(x, y, beta, dxx_shapelets, precalc=False)
            f_yy = self._shapeletOutput(x, y, beta, dyy_shapelets, precalc=False)
            f_xy = self._shapeletOutput(x, y, beta, dxy_shapelets, precalc=False)
        else:
            H_x, H_y = self.pre_calc(x, y, beta, n_order+2, center_x, center_y)
            f_xx = self._shapeletOutput(H_x, H_y, beta, dxx_shapelets)
            f_yy = self._shapeletOutput(H_x, H_y, beta, dyy_shapelets)
            f_xy = self._shapeletOutput(H_x, H_y, beta, dxy_shapelets)
        return f_xx, f_yy, f_xy
836155
def eval_hessian(self, *args, **kwargs):
        
        # Evaluate the hessian model and use the resulting Ans namedtuple as a
        # dict. From this, take the relevant components.
        eval_hess_dict = self.hessian_model(*args, **kwargs)._asdict()
        hess = [[[np.broadcast_to(eval_hess_dict.get(D(var, p1, p2), 0),
                                  eval_hess_dict[var].shape)
                    for p2 in self.params]
                for p1 in self.params]
            for var in self
        ]
        # Use numpy to broadcast these arrays together and then stack them along
        # the parameter dimension. We do not include the component direction in
        # this, because the components can have independent shapes.
        for idx, comp in enumerate(hess):
            hess[idx] = np.stack(np.broadcast_arrays(*comp))

        Ans = variabletuple('Ans', self.keys())
        return Ans(*hess)
659844
def hinit(func, x, t, pos_neg, f0, iord, hmax, rtol, atol, args):
    
    sk = atol + rtol * np.fabs(x)
    dnf = np.sum(np.square(f0 / sk), axis=0)
    dny = np.sum(np.square(x / sk), axis=0)

    h = np.sqrt(dny / dnf) * 0.01

    h = np.min([h, np.fabs(hmax)])
    h = custom_sign(h, pos_neg)

    # perform an explicit Euler step
    xx1 = x + h * f0
    f1 = np.array(func(xx1, t[0] + h, *args))

    # estimate the second derivative of the solution
    der2 = np.sum(np.square((f1 - f0) / sk), axis=0)
    der2 = np.sqrt(der2) / h

    # step size is computed such that h ** iord * max_d(norm(f0), norm(der2)) = 0.01
    der12 = np.max([np.fabs(der2), np.sqrt(dnf)])
    h1 = np.power(0.01 / der12, 1.0 / iord)

    h = np.min([100.0 * np.fabs(h), np.min([np.fabs(h1), np.fabs(hmax)])])

    return custom_sign(h, pos_neg), f0, f1, xx1
349564
function getGridDimension(diff, step) {

        // return step if diff = 0
        if (!diff) return step;

        var absDiff = Math.abs(diff);
        var numSteps = Math.round(absDiff / step);

        // return absDiff if less than one step apart
        if (!numSteps) return absDiff;

        // otherwise, return corrected step
        var roundedDiff = numSteps * step;
        var remainder = absDiff - roundedDiff;
        var stepCorrection = remainder / numSteps;

        return step + stepCorrection;
    }
1350945
def _getHessian(self):
        
        assert self.init,        'GP not initialised'
        assert self.fast==False, 'Not supported for fast implementation'

        if self.cache['Hessian'] is None:
            ParamMask=self.gp.getParamMask()['covar']
            std=sp.zeros(ParamMask.sum())
            H=self.gp.LMLhess_covar()
            It= (ParamMask[:,0]==1)
            self.cache['Hessian']=H[It,:][:,It]

        return self.cache['Hessian']
791235
def hessian(self, x, y, amp, sigma, e1, e2, center_x=0, center_y=0):
        
        alpha_ra, alpha_dec = self.derivatives(x, y, amp, sigma, e1, e2, center_x, center_y)
        diff = self._diff
        alpha_ra_dx, alpha_dec_dx = self.derivatives(x + diff, y, amp, sigma, e1, e2, center_x, center_y)
        alpha_ra_dy, alpha_dec_dy = self.derivatives(x, y + diff, amp, sigma, e1, e2, center_x, center_y)

        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
657239
def xstep(self):
        
        """

        self.YU[:] = self.Y - self.U
        self.block_sep0(self.YU)[:] += self.S
        Zf = sl.rfftn(self.YU, None, self.cri.axisN)
        Z0f = self.block_sep0(Zf)
        Z1f = self.block_sep1(Zf)

        DZ0f = np.conj(self.Df) * Z0f
        DZ0fBQ = sl.dot(self.B.dot(self.Q).T, DZ0f, axis=self.cri.axisC)
        Z1fQ = sl.dot(self.Q.T, Z1f, axis=self.cri.axisC)
        b = DZ0fBQ + Z1fQ

        Xh = sl.solvedbd_sm(self.gDf, (self.mu / self.rho) * self.GHGf + 1.0,
                            b, self.c, axis=self.cri.axisM)
        self.Xf[:] = sl.dot(self.Q, Xh, axis=self.cri.axisC)
        self.X = sl.irfftn(self.Xf, self.cri.Nv, self.cri.axisN)

        if self.opt['LinSolveCheck']:
            DDXf = np.conj(self.Df) *  sl.inner(self.Df, self.Xf,
                                                axis=self.cri.axisM)
            DDXfBB = sl.dot(self.B.T.dot(self.B), DDXf, axis=self.cri.axisC)
            ax = self.rho * (DDXfBB + self.Xf) + \
                 self.mu * self.GHGf * self.Xf
            b = self.rho * (sl.dot(self.B.T, DZ0f, axis=self.cri.axisC)
                            + Z1f)
            self.xrrs = sl.rrs(ax, b)
        else:
            self.xrrs = None
1209908
private void computeDir(double[] dir, double[] fg, QNInfo qn)
      throws SurpriseConvergence {
    System.arraycopy(fg, 0, dir, 0, fg.length);

    int mmm = qn.size();
    double[] as = new double[mmm];

    for (int i = mmm - 1; i >= 0; i--) {
      as[i] = qn.getRho(i) * ArrayMath.innerProduct(qn.getS(i), dir);
      plusAndConstMult(dir, qn.getY(i), -as[i], dir);
    }

    // multiply by hessian approximation
    qn.applyInitialHessian(dir);

    for (int i = 0; i < mmm; i++) {
      double b = qn.getRho(i) * ArrayMath.innerProduct(qn.getY(i), dir);
      plusAndConstMult(dir, qn.getS(i), as[i] - b, dir);
    }

    ArrayMath.multiplyInPlace(dir, -1);
  }