Lloyd's algorithm can be though of a diagonal-only Hessian; this method
    incorporates the diagonal blocks, too.
551042
def quasi_newton_uniform_blocks(points, cells, *args, **kwargs):
    

    def get_new_points(mesh):
        # TODO need copy?
        x = mesh.node_coords.copy()
        x += update(mesh)
        # update ghosts
        x[ghosted_mesh.is_ghost_point] = ghosted_mesh.reflect_ghost(
            x[ghosted_mesh.mirrors]
        )
        return x

    ghosted_mesh = GhostedMesh(points, cells)

    runner(
        get_new_points,
        ghosted_mesh,
        *args,
        **kwargs,
        update_topology=lambda mesh: ghosted_mesh.update_topology(),
        # get_stats_mesh=lambda mesh: ghosted_mesh.get_unghosted_mesh(),
    )

    mesh = ghosted_mesh.get_unghosted_mesh()
    return mesh.node_coords, mesh.cells["nodes"]
827941
public static LinearSolverDense<DMatrixRMaj> general(int numRows , int numCols ) {
        if( numRows == numCols )
            return linear(numRows);
        else
            return leastSquares(numRows,numCols);
    }
97089
public static void hessian(GrayF32 integral, int skip , int size ,
							   GrayF32 intensity)
	{
		// todo check size with skip
//		InputSanityCheck.checkSameShape(integral,intensity);

		if( BoofConcurrency.USE_CONCURRENT) {
			ImplIntegralImageFeatureIntensity_MT.hessianBorder(integral, skip, size, intensity);
			ImplIntegralImageFeatureIntensity_MT.hessianInner(integral, skip, size, intensity);
		} else {
			ImplIntegralImageFeatureIntensity.hessianBorder(integral, skip, size, intensity);
			ImplIntegralImageFeatureIntensity.hessianInner(integral, skip, size, intensity);
		}
	}
1172687
def find_neighbors(self, grid, node, diagonal_movement=None):
        '''
        
        '''
        if not diagonal_movement:
            diagonal_movement = self.diagonal_movement
        return grid.neighbors(node, diagonal_movement=diagonal_movement)
446161
private void setupA( List<DMatrixRMaj> homographies ) {
		A.reshape(2*homographies.size(),6, false);

		DMatrixRMaj h1 = new DMatrixRMaj(3,1);
		DMatrixRMaj h2 = new DMatrixRMaj(3,1);

		DMatrixRMaj v12 = new DMatrixRMaj(1,6);
		DMatrixRMaj v11 = new DMatrixRMaj(1,6);
		DMatrixRMaj v22 = new DMatrixRMaj(1,6);

		DMatrixRMaj v11m22 = new DMatrixRMaj(1,6);

		for( int i = 0; i < homographies.size(); i++ ) {
			DMatrixRMaj H = homographies.get(i);

			CommonOps_DDRM.extract(H,0,3,0,1,h1,0,0);
			CommonOps_DDRM.extract(H,0,3,1,2,h2,0,0);

			// normalize H by the max value to reduce numerical error when computing A
			// several numbers are multiplied against each other and could become quite large/small
			double max1 = CommonOps_DDRM.elementMaxAbs(h1);
			double max2 = CommonOps_DDRM.elementMaxAbs(h2);
			double max = Math.max(max1,max2);

			CommonOps_DDRM.divide(h1,max);
			CommonOps_DDRM.divide(h2,max);

			// compute elements of A
			computeV(h1, h2, v12);
			computeV(h1, h1, v11);
			computeV(h2, h2, v22);

			CommonOps_DDRM.subtract(v11, v22, v11m22);

			CommonOps_DDRM.insert( v12    , A, i*2   , 0);
			CommonOps_DDRM.insert( v11m22 , A, i*2+1 , 0);
		}
	}
1287249
static void backsubstitution0134(DMatrixRMaj P_plus, DMatrixRMaj P , DMatrixRMaj X ,
									 double H[] ) {
		final int N = P.numRows;
		DMatrixRMaj tmp = new DMatrixRMaj(N*2, 1);

		double H6 = H[6];
		double H7 = H[7];
		double H8 = H[8];

		for (int i = 0, index = 0; i < N; i++) {
			double x = -X.data[index],y = -X.data[index+1];
			double sum = P.data[index++] * H6 + P.data[index++] * H7 + H8;
			tmp.data[i] = x * sum;
			tmp.data[i+N] = y * sum;
		}

		double h0=0,h1=0,h3=0,h4=0;
		for (int i = 0; i < N; i++) {
			double P_pls_0 = P_plus.data[i];
			double P_pls_1 = P_plus.data[i+N];

			double tmp_i = tmp.data[i];
			double tmp_j = tmp.data[i+N];

			h0 += P_pls_0*tmp_i;
			h1 += P_pls_1*tmp_i;
			h3 += P_pls_0*tmp_j;
			h4 += P_pls_1*tmp_j;
		}

		H[0] = -h0;
		H[1] = -h1;
		H[3] = -h3;
		H[4] = -h4;
	}
1210662
def convex_hull(self):
        
        bnodes, nb, na, nt = _tripack.bnodes(self.lst, self.lptr, self.lend, self.npoints)
        return self._deshuffle_simplices(bnodes[:nb] - 1)
1046776
def hessian(self, x, y, theta_E, center_x=0, center_y=0):
        
        x_shift = x - center_x
        y_shift = y - center_y
        R = (x_shift*x_shift + y_shift*y_shift)**(3./2)
        if isinstance(R, int) or isinstance(R, float):
            prefac = theta_E / max(0.000001, R)
        else:
            prefac = np.empty_like(R)
            r = R[R>0]  #in the SIS regime
            prefac[R==0] = 0.
            prefac[R>0] = theta_E / r

        f_xx = y_shift*y_shift * prefac
        f_yy = x_shift*x_shift * prefac
        f_xy = -x_shift*y_shift * prefac
        return f_xx, f_yy, f_xy
766162
private void tql2() {
    // This is derived from the Algol procedures tql2, by
    // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
    // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
    // Fortran subroutine in EISPACK.

    System.arraycopy(e, 1, e, 0, n - 1);
    e[n - 1] = 0.;

    double f = 0., tst1 = 0.;
    for(int l = 0; l < n; l++) {
      // Find small subdiagonal element
      tst1 = max(tst1, abs(d[l]) + abs(e[l]));
      int m = l;
      for(; m < n; ++m) {
        if(abs(e[m]) <= 0x1P-52 * tst1) {
          break;
        }
      }

      // If m == l, d[l] is an eigenvalue,
      // otherwise, iterate.
      if(m > l) {
        do {
          // Compute implicit shift
          f += tql2ComputeImplicitShift(l);
          // Implicit QL transformation.
          tql2ImplicitQL(l, m, d[l + 1]);
          // Check for convergence.
          // TODO: Iteration limit?
        }
        while(abs(e[l]) > 0x1P-52 * tst1);
      }
      d[l] += f;
      e[l] = 0.;
    }
  }
747818
def kappa(self, x, y, kwargs, diff=diff):
        
        f_xx, f_xy, f_yx, f_yy = self.hessian(x, y, kwargs, diff=diff)
        kappa = 1./2 * (f_xx + f_yy)
        return kappa
153893
function diagonal(s, d) {
    if (s.dy == undefined) {
        s.dy = 0;
    }
    if (d.dy == undefined) {
        d.dy = 0;
    }
    let path = `M ${s.x} ${(s.y + s.dy)}
          C ${(s.x + d.x) / 2} ${(s.y + s.dy)},
            ${(s.x + d.x) / 2} ${(d.y + d.dy)},
            ${d.x} ${(d.y + d.dy)}`;
    return path;
}
494292
public static Basic1DMatrix diagonal(int size, double diagonal) {
        double[] array = new double[size * size];

        for (int i = 0; i < size; i++) {
            array[i * size + i] = diagonal;
        }

        return new Basic1DMatrix(size, size, array);
    }
1046616
private void hqr2(Matrix H, Matrix V, Vector d, Vector e) {

        // This is derived from the Algol procedure hqr2,
        // by Martin and Wilkinson, Handbook for Auto. Comp.,
        // Vol.ii-Linear Algebra, and the corresponding
        // Fortran subroutine in EISPACK.

        // Initialize

        int nn = e.length();
        int n = nn - 1;
        int low = 0;
        int high = nn - 1;
        double eps = Math.pow(2.0, -52.0);
        double exshift = 0.0;
        double p = 0;
        double q = 0;
        double r = 0;
        double s = 0;
        double z = 0;
        double t;
        double w;
        double x;
        double y;

        // Store roots isolated by balanc and compute matrix norm

        double norm = 0.0;
        for (int i = 0; i < nn; i++) {
            if (i < low | i > high) {
                d.set(i, H.get(i, i));
                e.set(i, 0.0);
            }
            for (int j = Math.max(i - 1, 0); j < nn; j++) {
                norm = norm + Math.abs(H.get(i, j));
            }
        }

        // Outer loop over eigenvalue index

        int iter = 0;
        while (n >= low) {

            // Look for single small sub-diagonal element

            int l = n;
            while (l > low) {
                s = Math.abs(H.get(l - 1, l - 1)) 
                    + Math.abs(H.get(l, l));

                if (s == 0.0) {
                    s = norm;
                }
                if (Math.abs(H.get(l, l - 1)) < eps * s) {
                    break;
                }
                l--;
            }

            // Check for convergence
            // One root found

            if (l == n) {
                H.updateAt(n, n, Matrices.asPlusFunction(exshift));
                d.set(n, H.get(n, n));
                e.set(n, 0.0);
                n--;
                iter = 0;

                // Two roots found

            } else if (l == n - 1) {
                w = H.get(n, n - 1) * H.get(n - 1, n);
                p = (H.get(n - 1, n - 1) - H.get(n, n)) / 2.0;
                q = p * p + w;
                z = Math.sqrt(Math.abs(q));
                H.updateAt(n, n, Matrices.asPlusFunction(exshift));
                H.updateAt(n - 1, n - 1, Matrices.asPlusFunction(exshift));
                x = H.get(n, n);

                // Real pair

                if (q >= 0) {
                    if (p >= 0) {
                        z = p + z;
                    } else {
                        z = p - z;
                    }
                    d.set(n - 1, x + z);
                    d.set(n, d.get(n - 1));
                    if (z != 0.0) {
                        d.set(n, x - w / z);
                    }
                    e.set(n - 1, 0.0);
                    e.set(n, 0.0);
                    x = H.get(n, n - 1);
                    s = Math.abs(x) + Math.abs(z);
                    p = x / s;
                    q = z / s;
                    r = Math.sqrt(p * p + q * q);
                    p = p / r;
                    q = q / r;

                    // Row modification

                    for (int j = n - 1; j < nn; j++) {
                        z = H.get(n - 1, j);
                        H.set(n - 1, j, q * z + p * H.get(n, j));
                        H.set(n, j, q * H.get(n, j) - p * z);
                    }

                    // Column modification

                    for (int i = 0; i <= n; i++) {
                        z = H.get(i, n - 1);
                        H.set(i, n - 1, q * z + p * H.get(i, n));
                        H.set(i, n, q * H.get(i, n) - p * z);
                    }

                    // Accumulate transformations

                    for (int i = low; i <= high; i++) {
                        z = V.get(i, n - 1);
                        V.set(i, n - 1, q * z + p * V.get(i, n));
                        V.set(i, n, q * V.get(i, n) - p * z);
                    }

                    // Complex pair

                } else {
                    d.set(n - 1, x + p);
                    d.set(n, x + p);
                    e.set(n - 1, z);
                    e.set(n, -z);
                }
                n = n - 2;
                iter = 0;

                // No convergence yet

            } else {

                // Form shift

                x = H.get(n, n);
                y = 0.0;
                w = 0.0;
                if (l < n) {
                    y = H.get(n - 1, n - 1);
                    w = H.get(n, n - 1) * H.get(n - 1, n);
                }

                // Wilkinson's original ad hoc shift

                if (iter == 10) {
                    exshift += x;
                    for (int i = low; i <= n; i++) {
                        H.updateAt(i, i, Matrices.asMinusFunction(x));
                    }
                    s = Math.abs(H.get(n, n - 1)) 
                        + Math.abs(H.get(n - 1, n - 2));
                    x = y = 0.75 * s; // WTF ???
                    w = -0.4375 * s * s; // Are you kidding me???
                }

                // MATLAB's new ad hoc shift

                if (iter == 30) {
                    s = (y - x) / 2.0;
                    s = s * s + w;
                    if (s > 0) {
                        s = Math.sqrt(s);
                        if (y < x) {
                            s = -s;
                        }
                        s = x - w / ((y - x) / 2.0 + s);
                        for (int i = low; i <= n; i++) {
                            H.updateAt(i, i, Matrices.asMinusFunction(s));
                        }
                        exshift += s;
                        x = y = w = 0.964;
                    }
                }

                iter = iter + 1; // (Could check iteration count here.)

                // Look for two consecutive small sub-diagonal elements

                int m = n - 2;
                while (m >= l) {
                    z = H.get(m, m);
                    r = x - z;
                    s = y - z;
                    p = (r * s - w) / H.get(m + 1, m) 
                                      + H.get(m, m + 1);
                    q = H.get(m + 1, m + 1) - z - r - s;
                    r = H.get(m + 2, m + 1);
                    s = Math.abs(p) + Math.abs(q) + Math.abs(r);
                    p = p / s;
                    q = q / s;
                    r = r / s;
                    if (m == l) {
                        break;
                    }
                    if (Math.abs(H.get(m, m - 1)) * (Math.abs(q) + Math.abs(r)) < eps
                            * (Math.abs(p) * (Math.abs(H.get(m - 1, m - 1))
                                    + Math.abs(z) + Math.abs(H.get(m + 1, m + 1))))) {
                        break;
                    }
                    m--;
                }

                for (int i = m + 2; i <= n; i++) {
                    H.set(i, i - 2, 0.0);
                    if (i > m + 2) {
                        H.set(i, i - 3, 0.0);
                    }
                }

                // Double QR step involving rows l:n and columns m:n

                for (int k = m; k <= n - 1; k++) {
                    boolean notlast = (k != n - 1);
                    if (k != m) {
                        p = H.get(k, k - 1);
                        q = H.get(k + 1, k - 1);
                        r = (notlast ? H.get(k + 2, k - 1) : 0.0);
                        x = Math.abs(p) + Math.abs(q) + Math.abs(r);
                        if (x == 0.0) {
                            continue;
                        }
                        p = p / x;
                        q = q / x;
                        r = r / x;
                    }

                    s = Math.sqrt(p * p + q * q + r * r);
                    if (p < 0) {
                        s = -s;
                    }
                    if (s != 0) {
                        if (k != m) {
                            H.set(k, k - 1, -s * x);
                        } else if (l != m) {
                            H.updateAt(k, k - 1, Matrices.INV_FUNCTION);
                        }
                        p = p + s;
                        x = p / s;
                        y = q / s;
                        z = r / s;
                        q = q / p;
                        r = r / p;

                        // Row modification

                        for (int j = k; j < nn; j++) {
                            p = H.get(k, j) + q * H.get(k + 1, j);
                            if (notlast) {
                                p = p + r * H.get(k + 2, j);
                                H.updateAt(k + 2, j,
                                    Matrices.asMinusFunction(p * z));
                            }
                            H.updateAt(k, j, Matrices.asMinusFunction(p * x));
                            H.updateAt(k + 1, j, Matrices.asMinusFunction(p * y));
                        }

                        // Column modification

                        for (int i = 0; i <= Math.min(n, k + 3); i++) {
                            p = x * H.get(i, k) + y 
                                * H.get(i, k + 1);
                            if (notlast) {
                                p = p + z * H.get(i, k + 2);
                                H.updateAt(i, k + 2,
                                    Matrices.asMinusFunction(p * r));
                            }
                            H.updateAt(i, k, Matrices.asMinusFunction(p));
                            H.updateAt(i, k + 1, Matrices.asMinusFunction(p * q));
                        }

                        // Accumulate transformations

                        for (int i = low; i <= high; i++) {
                            p = x * V.get(i, k) + y 
                                * V.get(i, k + 1);
                            if (notlast) {
                                p = p + z * V.get(i, k + 2);
                                V.updateAt(i, k + 2,
                                    Matrices.asMinusFunction(p * r));
                            }
                            V.updateAt(i, k, Matrices.asMinusFunction(p));
                            V.updateAt(i, k + 1, Matrices.asMinusFunction(p * q));
                        }
                    } // (s != 0)
                } // k loop
            } // check convergence
        } // while (n >= low)

        // Backsubstitute to find vectors of upper triangular form

        if (norm == 0.0) {
            return;
        }

        for (n = nn - 1; n >= 0; n--) {
            p = d.get(n);
            q = e.get(n);

            // Real vector

            if (q == 0) {
                int l = n;
                H.set(n, n, 1.0);
                for (int i = n - 1; i >= 0; i--) {
                    w = H.get(i, i) - p;
                    r = 0.0;
                    for (int j = l; j <= n; j++) {
                        r = r + H.get(i, j) * H.get(j, n);
                    }
                    if (e.get(i) < 0.0) {
                        z = w;
                        s = r;
                    } else {
                        l = i;
                        if (e.get(i) == 0.0) {
                            if (w != 0.0) {
                                H.set(i, n, -r / w);
                            } else {
                                H.set(i, n, -r / (eps * norm));
                            }

                            // Solve real equations

                        } else {
                            x = H.get(i, i + 1);
                            y = H.get(i + 1, i);
                            q = (d.get(i) - p) * (d.get(i) - p) 
                                + e.get(i) * e.get(i);
                            t = (x * s - z * r) / q;
                            H.set(i, n, t);
                            if (Math.abs(x) > Math.abs(z)) {
                                H.set(i + 1, n, (-r - w * t) / x);
                            } else {
                                H.set(i + 1, n, (-s - y * t) / z);
                            }
                        }

                        // Overflow control

                        t = Math.abs(H.get(i, n));
                        if ((eps * t) * t > 1) {
                            for (int j = i; j <= n; j++) {
                                H.updateAt(j, n, Matrices.asDivFunction(t));
                            }
                        }
                    }
                }

                // Complex vector

            } else if (q < 0) {
                int l = n - 1;

                // Last vector component imaginary so matrix is triangular

                if (Math.abs(H.get(n, n - 1)) 
                        > Math.abs(H.get(n - 1, n))) {
                    H.set(n - 1, n - 1, q / H.get(n, n - 1));
                    H.set(n - 1, n, -(H.get(n, n) - p) 
                                            / H.get(n, n - 1));
                } else {
                    double[] cdiv = cdiv(0.0, -H.get(n - 1, n),
                         H.get(n - 1, n - 1) - p, q);
                    H.set(n - 1, n - 1, cdiv[0]);
                    H.set(n - 1, n, cdiv[1]);
                }
                H.set(n, n - 1, 0.0);
                H.set(n, n, 1.0);
                for (int i = n - 2; i >= 0; i--) {
                    double ra;
                    double sa;
                    double vr;
                    double vi;
                    ra = 0.0;
                    sa = 0.0;
                    for (int j = l; j <= n; j++) {
                        ra = ra + H.get(i, j) * H.get(j, n - 1);
                        sa = sa + H.get(i, j) * H.get(j, n);
                    }
                    w = H.get(i, i) - p;

                    if (e.get(i) < 0.0) {
                        z = w;
                        r = ra;
                        s = sa;
                    } else {
                        l = i;
                        if (e.get(i) == 0) {
                            double[] cdiv = cdiv(-ra, -sa, w, q);
                            H.set(i, n - 1, cdiv[0]);
                            H.set(i, n, cdiv[1]);
                        } else {

                            // Solve complex equations

                            x = H.get(i, i + 1);
                            y = H.get(i + 1, i);
                            vr = (d.get(i) - p) * (d.get(i) - p) 
                                 + e.get(i) * e.get(i) - q * q;
                            vi = (d.get(i) - p) * 2.0 * q;
                            if (vr == 0.0 & vi == 0.0) {
                                vr = eps
                                        * norm
                                        * (Math.abs(w) + Math.abs(q)
                                                + Math.abs(x) + Math.abs(y) + Math
                                                    .abs(z));
                            }
                            double[] cdiv = cdiv(x * r - z * ra + q * sa,
                                                 x * s - z * sa - q * ra, vr, vi);
                            H.set(i, n - 1, cdiv[0]);
                            H.set(i, n, cdiv[1]);
                            if (Math.abs(x) > (Math.abs(z) + Math.abs(q))) {
                                H.set(i + 1, n - 1, (-ra - w 
                                             * H.get(i, n - 1) + q
                                             * H.get(i, n)) / x);
                                H.set(i + 1, n, (-sa - w 
                                             * H.get(i, n) - q 
                                             * H.get(i, n - 1)) / x);
                            } else {
                                cdiv = cdiv(-r - y 
                                     * H.get(i, n - 1), -s - y 
                                     * H.get(i, n), z, q);
                                H.set(i + 1, n - 1, cdiv[0]);
                                H.set(i + 1, n, cdiv[1]);
                            }
                        }

                        // Overflow control

                        t = Math.max(Math.abs(H.get(i, n - 1)), 
                                     Math.abs(H.get(i, n)));
                        if ((eps * t) * t > 1) {
                            for (int j = i; j <= n; j++) {
                                H.updateAt(j, n - 1, Matrices.asDivFunction(t));
                                H.updateAt(j, n, Matrices.asDivFunction(t));
                            }
                        }
                    }
                }
            }
        }

        // Vectors of isolated roots

        for (int i = 0; i < nn; i++) {
            if (i < low | i > high) {
                for (int j = i; j < nn; j++) {
                    V.set(i, j, H.get(i, j));
                }
            }
        }

        // Back transformation to get eigenvectors of original matrix

        for (int j = nn - 1; j >= low; j--) {
            for (int i = low; i <= high; i++) {
                z = 0.0;
                for (int k = low; k <= Math.min(j, high); k++) {
                    z = z + V.get(i, k) * H.get(k, j);
                }
                V.set(i, j, z);
            }
        }
    }
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
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
992882
public static void multPlusTransA(int blockLength ,
                                      DSubmatrixD1 A , DSubmatrixD1 B ,
                                      DSubmatrixD1 C )
    {
        int heightA = Math.min( blockLength , A.row1 - A.row0 );

        for( int i = C.row0+blockLength; i < C.row1; i += blockLength ) {
            int heightC = Math.min( blockLength , C.row1 - i );

            int indexA = A.row0*A.original.numCols + (i-C.row0+A.col0)*heightA;

            for( int j = i; j < C.col1; j += blockLength ) {
                int widthC = Math.min( blockLength , C.col1 - j );

                int indexC = i*C.original.numCols + j*heightC;
                int indexB = B.row0*B.original.numCols + (j-C.col0+B.col0)*heightA;

                blockMultPlusTransA(A.original.data,B.original.data,C.original.data,
                            indexA,indexB,indexC,heightA,heightC,widthC);
            }
        }
    }
370350
def _to_diagonally_dominant_weighted(mat):
    """"""
    mat += np.diag(np.sum(np.abs(mat), axis=1) + 0.01)
    return mat