Function calculates sigmoid function (g) for every sample. With g it calculates part of Hessian matrix and gradient, aggregates and output them. It also calculates J function which is needed for checking the convergence of parameters theta.
975472
def map_fit(interface, state, label, inp):
    
    import numpy as np
    out = interface.output(0)

    H, J, grad = 0, 0, 0
    for row in inp:
        row = row.strip().split(state["delimiter"])  # split row
        if len(row) > 1:  # check if row is empty
            # add intercept term to every sample
            x = np.array([1] + [(0 if v in state["missing_vals"] else float(v)) for i, v in enumerate(row) if
                                i in state["X_indices"]])
            # map label value to 0 or 1. If label does not match set error
            y = 0 if state["y_map"][0] == row[state["y_index"]] else 1 if state["y_map"][1] == row[
                state["y_index"]] else "Error"

            g = 1. / (1 + np.exp(-np.dot(x, state["thetas"])))  # sigmoid function
            grad += x * (g - y)  # gradient
            H += np.multiply(np.outer(x, x), g * (1 - g))  # Hessian matrix
            J -= np.log(g) if y == 1 else np.log(1 - g)  # J cost function
    out.add("grad", grad)
    out.add("J", J)
    for i, row in enumerate(H):
        out.add(i, row)
1204850
def probability_density(self, X):
        
        self.check_fit()

        U, V = self.split_matrix(X)

        if self.theta == 0:
            return np.multiply(U, V)

        else:
            num = np.multiply(np.multiply(-self.theta, self._g(1)), 1 + self._g(np.add(U, V)))
            aux = np.multiply(self._g(U), self._g(V)) + self._g(1)
            den = np.power(aux, 2)
            return num / den
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
221584
def _g_a (self, x, a, b, s):
        
        return 2*b/(1.0+np.exp(s*(x-a)))
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;
        }
    ''')])
405450
def get_G(self, k, add_noise=True):
        
        G = np.ones((self.N - 1, self.N - 1))
        if (add_noise):
            noise = pi * 0.1 * np.random.rand(
                (self.N - 1) * (self.N - 1)).reshape((self.N - 1, self.N - 1))
        other_indices = np.delete(range(self.N), k)
        for idx, i in enumerate(other_indices):
            for jdx, j in enumerate(other_indices):
                if (add_noise and
                        i != j):  # do not add noise on diagonal elements.
                    thetak_ij = self.get_inner_angle(k,
                                                     (i, j)) + noise[idx, jdx]
                else:
                    thetak_ij = self.get_inner_angle(k, (i, j))
                G[idx, jdx] = cos(thetak_ij)
                G[jdx, idx] = cos(thetak_ij)
        return G
509111
private boolean computeSearchDirection() {
		// Compute the function's gradient
		function.computeGradient(temp0_Nx1.data);

		// compute the change in gradient
		for( int i = 0; i < N; i++ ) {
			y.data[i] = temp0_Nx1.data[i] - g.data[i];
			g.data[i] = temp0_Nx1.data[i];
		}

		// Update the inverse Hessian matrix
		if( iterations != 0 ) {
			EquationsBFGS.inverseUpdate(B, s, y, temp0_Nx1, temp1_Nx1);
		}

		// compute the search direction
		CommonOps_DDRM.mult(-1,B,g, searchVector);

		// use the line search to find the next x
		if( !setupLineSearch(fx, x.data, g.data, searchVector.data) ) {
			// the search direction has a positive derivative, meaning the B matrix is
			// no longer SPD.  Attempt to fix the situation by resetting the matrix
			resetMatrixB();
			// do the search again, it can't fail this time
			CommonOps_DDRM.mult(-1,B,g, searchVector);
			setupLineSearch(fx, x.data, g.data, searchVector.data);
		} else if(Math.abs(derivAtZero) <= gtol ) {
			if( verbose != null ) {
				verbose.printf("finished select direction, gtest=%e\n",Math.abs(derivAtZero));
			}

			// the input might have been modified by the function.  So copy it
			System.arraycopy(function.getCurrentState(),0,x.data,0,N);
			return terminateSearch(true);
		}

		mode = 1;
		iterations++;
		return false;
	}
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
200451
public double computeAverage (Function<Triple<GeneralDataset<L, F>,GeneralDataset<L, F>,SavedState>,Double> function) 
  {
    double sum = 0;
    Iterator<Triple<GeneralDataset<L, F>,GeneralDataset<L, F>,SavedState>> foldIt = iterator();
    while (foldIt.hasNext()) {
      sum += function.apply(foldIt.next());
    }
    return sum / kfold;
  }
707982
def F_t(X, Y, S, M_E, E, m0, rho):
    '''  '''
    r = X.shape[1]
    out1 = (((np.dot(np.dot(X, S), Y.T) - M_E) * E)**2).sum() / 2
    out2 = rho * G(Y, m0, r)
    out3 = rho * G(X, m0, r)

    return out1 + out2 + out3
442663
def _g(self, h, xp, s):
        """"""

        nphi = sum(self.phi)

        return (nphi / 2.0) * log(2 * pi) + nphi * \
            log(s) + 0.5 * sum((h - xp) ** 2) / (s ** 2)
1059980
def hessian(self, x, y, kappa_ext, ra_0=0, dec_0=0):
        
        gamma1 = 0
        gamma2 = 0
        kappa = kappa_ext
        f_xx = kappa + gamma1
        f_yy = kappa - gamma1
        f_xy = gamma2
        return f_xx, f_yy, f_xy
110284
private double computeOutputError(double output, double[] gradient) {

        double error = 0.0;
        double out = outputLayer.output[0];
        double g = output - out;

        error += (0.5*g * g);


        gradient[0] = g;

        return error;
    }
641132
def hessian(self, x, y, theta_E, r_trunc, center_x=0, center_y=0):
        
        x_shift = x - center_x
        y_shift = y - center_y
        dphi_dr = self._dphi_dr(x_shift, y_shift, theta_E, r_trunc)
        d2phi_dr2 = self._d2phi_dr2(x_shift, y_shift, theta_E, r_trunc)
        dr_dx, dr_dy = self._dr_dx(x, y)
        d2r_dx2, d2r_dy2, d2r_dxy = self._d2r_dx2(x_shift, y_shift)
        f_xx = d2r_dx2*dphi_dr + dr_dx**2*d2phi_dr2
        f_yy = d2r_dy2*dphi_dr + dr_dy**2*d2phi_dr2
        f_xy = d2r_dxy*dphi_dr + dr_dx*dr_dy*d2phi_dr2
        return f_xx, f_yy, f_xy
1094609
public static double sigmoid(double x) {
    if (x<0) {
      double num = Math.exp(x);
      return num / (1.0 + num);
    }
    else {
      double den = 1.0 + Math.exp(-x);
      return 1.0 / den;
    }
  }
741166
function bgrat(a, b, x) {
      var T, fTs, u, logH, logM, lnx2s, lnx2sn, ps, cs, j, ser;

      T = a + 0.5 * (b - 1);
      fTs = 1 / (4 * T * T);
      u = -T * Math.log(x);
      logH = Math.log(b) + log1p(gam1(b)) + b * Math.log(u) - u;
      logM = logH - (algdiv(b, a) + b * Math.log(T));
      cs = [1];   // cs[n] = 1 / (2n + 1)!
      ps = [1];   // ps[n] = p_n in (9.3)
      lnx2s = Math.log(x) * 0.5;   // (ln(x) / 2) ^ 2
      lnx2s *= lnx2s;
      lnx2sn = 1;   // lnx2s ^ n
      j = gratR(b, u, logH);
      ser =
      Math.log(
         series(function(n) {
            var mbmn; // holds m*b - n

            mbmn = -n;
            if (n > 0) {
               j = (b + 2 * n - 2) * (b + 2 * n - 2 + 1) * j +
                     (u + b + 2 * n - 2 + 1) * lnx2sn;
               j *= fTs;
               cs[n] = cs[n - 1] / (2 * n) / (2 * n + 1);
               ps[n] = (b - 1) * cs[n] + 1 / n * series(function(m) {
                  if (m === 0) { return 0; }
                  mbmn += b;
                  return mbmn * cs[m] * ps[n - m];
               }, n);
               lnx2sn *= lnx2s;
            }
            return ps[n] * j;
         })
      );
      return logM + ser;
   }
935006
def checkgrad(f, fprime, x, *args,**kw_args):
    
    LG.debug("Checking gradient ...")
    import numpy as np

    # using machine precision to choose h
    eps = np.finfo(float).eps
    step = np.sqrt(eps)*(x.min())

    # shake things up a bit by taking random steps for each x dimension
    h = step*np.sign(np.random.uniform(-1, 1, x.size))

    f_ph = f(x+h, *args, **kw_args)
    f_mh = f(x-h, *args, **kw_args)
    numerical_gradient = (f_ph - f_mh)/(2*h)
    analytical_gradient = fprime(x, *args, **kw_args)
    ratio = (f_ph - f_mh)/(2*np.dot(h, analytical_gradient))

    h = np.zeros_like(x)
    for i in range(len(x)):
        pdb.set_trace()
    h[i] = step
    f_ph = f(x+h, *args, **kw_args)
    f_mh = f(x-h, *args, **kw_args)
    numerical_gradient = (f_ph - f_mh)/(2*step)
    analytical_gradient = fprime(x, *args, **kw_args)[i]
    ratio = (f_ph - f_mh)/(2*step*analytical_gradient)
    h[i] = 0
    LG.debug("[%d] numerical: %f, analytical: %f, ratio: %f" % (i, numerical_gradient,analytical_gradient,ratio))