Minimum and maximum angular errors
        corresponding to 1st and 2nd axes
        of PCA distribution.
320829
def angular_errors(self, degrees=True):
        
        hyp_axes = self.method(self)
        v = angular_errors(hyp_axes)
        if degrees:
            v = N.degrees(v)
        return tuple(v)
562859
def error_asymptotes(pca,**kwargs):
    
    ax = kwargs.pop("ax",current_axes())

    lon,lat = pca.plane_errors('upper', n=1000)
    ax.plot(lon,lat,'-')

    lon,lat = pca.plane_errors('lower', n=1000)
    ax.plot(lon,lat,'-')

    ax.plane(*pca.strike_dip())
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;
                }
            }
        }
    }
582889
def compute_symm_block_tridiag_covariances(H_diag, H_upper_diag):
    
    T, D, _ = H_diag.shape
    assert H_diag.ndim == 3 and H_diag.shape[2] == D
    assert H_upper_diag.shape == (T - 1, D, D)

    J_init = J_11 = J_22 = np.zeros((D, D))
    h_init = h_1 = h_2 = np.zeros((D,))

    J_21 = np.swapaxes(H_upper_diag, -1, -2)
    J_node = H_diag
    h_node = np.zeros((T, D))

    _, _, sigmas, E_xt_xtp1 = \
        info_E_step(J_init, h_init, 0,
                    J_11, J_21, J_22, h_1, h_2, np.zeros((T-1)),
                    J_node, h_node, np.zeros(T))
    return sigmas, E_xt_xtp1
986781
def _axis_in_detector(geometry):
    """"""
    du, dv = geometry.det_axes_init
    axis = geometry.axis
    c = np.array([np.vdot(axis, du), np.vdot(axis, dv)])
    cnorm = np.linalg.norm(c)

    # Check for numerical errors
    assert cnorm != 0

    return c / cnorm
585555
def _check_r(self, r):
        """"""
        if abs(np.dot(r[:, 0], r[:, 0]) - 1) > eps or \
            abs(np.dot(r[:, 0], r[:, 0]) - 1) > eps or \
            abs(np.dot(r[:, 0], r[:, 0]) - 1) > eps or \
            np.dot(r[:, 0], r[:, 1]) > eps or \
            np.dot(r[:, 1], r[:, 2]) > eps or \
            np.dot(r[:, 2], r[:, 0]) > eps:
            raise ValueError("The rotation matrix is significantly non-orthonormal.")
154473
def plot_aligned(pca, sparse=True, **kwargs):
    
    colormap = kwargs.pop('colormap',None)

    A = pca.rotated()
    # Flip the result if upside down
    if pca.normal[2] < 0:
        A[:,-1] *= -1
    minmax = [(A[:,i].min(),A[:,i].max()) for i in range(3)]
    lengths = [j-i for i,j in minmax]

    if sparse:
        i = 1
        l = len(A)
        if l > 10000:
            i = N.ceil(l/10000)
            A = A[::int(i)]

    log.info("Plotting with {} points".format(len(A)))

    w = 8
    ratio = (lengths[2]*2+lengths[1])/lengths[0]*1.5
    h = w*ratio
    if h < 3: h = 3

    r = (lengths[1]+5)/(lengths[2]+5)
    if r > 5:
        r = 5

    fig = Figure(figsize=(w,h))
    fig.canvas = FigureCanvas(fig)

    def setup_axes():
        gs = GridSpec(3,1, height_ratios=[r,1,1])
        kwargs = dict()
        axes = []
        for g in gs:
            ax = fig.add_subplot(g,**kwargs)
            kwargs['sharex'] = ax
            yield ax

    axes = list(setup_axes())

    fig.subplots_adjust(hspace=0, wspace=0.1)

    #lengths = attitude.pca.singular_values[::-1]

    titles = (
        "Plan view",
        "Long cross-section",
        "Short cross-section")

    ylabels = (
        "Meters",
        "Residuals (m)",
        "Residuals (m)")

    colors = ['cornflowerblue','red']

    hyp = sampling_axes(pca)
    vertical = vector(0,0,1)

    for title,ax,(a,b),ylabel in zip(titles,axes,
            [(0,1),(0,2),(1,2)],ylabels):

        kw = dict(linewidth=2, alpha=0.5)
        bounds = minmax[a]
        if b != 2:
            ax.plot(bounds,(0,0), c=colors[a], **kw)
            ax.plot((0,0),minmax[b], c=colors[b], **kw)
        else:
            ax.plot(bounds,(-10,-10), c=colors[a], **kw)
            v0 = N.zeros(3)
            v0[a] = 1
            axes = N.array([v0,vertical])

            ax_ = (axes@N.diag(hyp)@axes.T).T
            ax_ = N.sqrt(ax_)
            l1 = N.linalg.norm(ax_[:-1])
            l2 = N.linalg.norm(ax_[-1])
            ang_error = 2*N.degrees(N.arctan2(l2,l1))

            title += ": {:.0f} m long, angular error (95% CI): {:.2f}º".format(lengths[a], ang_error)

            bounds = minmax[0]
            x_ = N.linspace(bounds[0]*1.2,bounds[1]*1.2,100)

            err = HyperbolicErrors(hyp,x_,axes=axes)
            err.plot(ax, fc='#cccccc', alpha=0.3)

        x,y = A[:,a], A[:,b]
        kw = dict(alpha=0.5, zorder=5)

        if colormap is None:
            ax.plot(x,y,c="#555555", linestyle='None', marker='.',**kw)
        else:
            ax.scatter(x,y,c=A[:,-1], cmap=colormap, **kw)

        #ax.set_aspect("equal")

        ax.text(0.01,.99,title,
            verticalalignment='top',
            transform=ax.transAxes)
        #ax.autoscale(tight=True)
        ax.yaxis.set_ticks([])
        ax.xaxis.set_ticks_position('bottom')
        if a != 1:
            pass
            #ax.xaxis.set_ticks([])
            #ax.spines['bottom'].set_color('none')
        for spine in ax.spines.values():
            spine.set_visible(False)


    ax.text(0.99,0.99,"Max residual: {:.1f} m".format(lengths[2]),
        verticalalignment='bottom',
        ha='right',
        transform=ax.transAxes)


    ax.set_xlabel("Meters")
    return fig
231217
private static void gaussian(float a[][],
	                             int index[]) {
		int n = index.length;
		float c[] = new float[n];
		// Initialize the index
		for (int i = 0; i < n; ++i) {
			index[i] = i;
		}
		// Find the rescaling factors, one from each row
		for (int i = 0; i < n; ++i){
			float c1 = 0;
			for (int j = 0; j < n; ++j){
				float c0 = Math.abs(a[i][j]);
				if (c0 > c1) {
					c1 = c0;
				}
			}
			c[i] = c1;
		}
		// Search the pivoting element from each column
		int k = 0;
		for (int j = 0; j < n - 1; ++j){
			float pi1 = 0;
			for (int i = j; i < n; ++i){
				float pi0 = Math.abs(a[index[i]][j]);
				pi0 /= c[index[i]];
				if (pi0 > pi1) {

					pi1 = pi0;
					k = i;
				}
			}
			// Interchange rows according to the pivoting order
			int itmp = index[j];
			index[j] = index[k];
			index[k] = itmp;
			for (int i = j + 1; i < n; ++i){
				float pj = a[index[i]][j] / a[index[j]][j];
				// Record pivoting ratios below the diagonal
				a[index[i]][j] = pj;

				// Modify other elements accordingly
				for (int l = j + 1; l < n; ++l){
					a[index[i]][l] -= pj * a[index[j]][l];
				}
			}
		}
	}
1350223
def _symmetrize_correlograms(correlograms):
    """"""

    n_clusters, _, n_bins = correlograms.shape
    assert n_clusters == _

    # We symmetrize c[i, j, 0].
    # This is necessary because the algorithm in correlograms()
    # is sensitive to the order of identical spikes.
    correlograms[..., 0] = np.maximum(correlograms[..., 0],
                                      correlograms[..., 0].T)

    sym = correlograms[..., 1:][..., ::-1]
    sym = np.transpose(sym, (1, 0, 2))

    return np.dstack((sym, correlograms))
851778
def feature_reduction(data, method, n_features):
    
    if method == "PCA":
        feature_red_method = sklearn.decomposition.PCA(n_components=n_features)
        data_processed = feature_red_method.fit_transform(data)

    elif method == "agglom":
        feature_red_method = sklearn.cluster.FeatureAgglomeration(n_clusters=n_features)
        data_processed = feature_red_method.fit_transform(data)

    elif method == "ica":
        feature_red_method = sklearn.decomposition.FastICA(n_components=n_features)
        data_processed = feature_red_method.fit_transform(data)

    elif method == "kernelPCA":
        feature_red_method = sklearn.decomposition.KernelPCA(n_components=n_features, kernel='linear')
        data_processed = feature_red_method.fit_transform(data)

    elif method == "kernelPCA":
        feature_red_method = sklearn.decomposition.KernelPCA(n_components=n_features, kernel='linear')
        data_processed = feature_red_method.fit_transform(data)

    elif method == "sparsePCA":
        feature_red_method = sklearn.decomposition.SparsePCA(n_components=n_features)
        data_processed = feature_red_method.fit_transform(data)

    elif method == "incrementalPCA":
        feature_red_method = sklearn.decomposition.IncrementalPCA(n_components=n_features)
        data_processed = feature_red_method.fit_transform(data)

    elif method == "nmf":
        if np.min(data) < 0:
            data -= np.min(data)
        feature_red_method = sklearn.decomposition.NMF(n_components=n_features)
        data_processed = feature_red_method.fit_transform(data)

    else:
        feature_red_method = None
        data_processed = data.copy()

    return(data_processed)
609234
void computeRawDescriptor(double c_x, double c_y, double sigma, double orientation, TupleDesc_F64 descriptor) {
		double c = Math.cos(orientation);
		double s = Math.sin(orientation);

		float fwidthSubregion = widthSubregion;
		int sampleWidth = widthGrid*widthSubregion;
		double sampleRadius = sampleWidth/2;

		double sampleToPixels = sigma*sigmaToPixels;

		Deriv image = (Deriv)imageDerivX.getImage();

		for (int sampleY = 0; sampleY < sampleWidth; sampleY++) {
			float subY = sampleY/fwidthSubregion;
			double y = sampleToPixels*(sampleY-sampleRadius);

			for (int sampleX = 0; sampleX < sampleWidth; sampleX++) {
				// coordinate of samples in terms of sub-region.  Center of sample point, hence + 0.5f
				float subX = sampleX/fwidthSubregion;
				// recentered local pixel sample coordinate
				double x = sampleToPixels*(sampleX-sampleRadius);

				// pixel coordinate in the image that is to be sampled.  Note the rounding
				// If the pixel coordinate is -1 < x < 0 then it will round to 0 instead of -1, but the rounding
				// method below is WAY faster than Math.round() so this is a small loss.
				int pixelX = (int)(x*c - y*s + c_x + 0.5);
				int pixelY = (int)(x*s + y*c + c_y + 0.5);

				// skip pixels outside of the image
				if( image.isInBounds(pixelX,pixelY) ) {
					// spacial image derivative at this point
					float spacialDX = imageDerivX.unsafe_getF(pixelX, pixelY);
					float spacialDY = imageDerivY.unsafe_getF(pixelX, pixelY);

					double adjDX =  c*spacialDX + s*spacialDY;
					double adjDY = -s*spacialDX + c*spacialDY;

					double angle = UtilAngle.domain2PI(Math.atan2(adjDY,adjDX));

					float weightGaussian = gaussianWeight[sampleY*sampleWidth+sampleX];
					float weightGradient = (float)Math.sqrt(spacialDX*spacialDX + spacialDY*spacialDY);

					// trilinear interpolation intro descriptor
					trilinearInterpolation(weightGaussian*weightGradient,subX,subY,angle, descriptor);
				}
			}
		}
	}
416816
def PCA_components(x):
     
    # validate inputs
    try:    
        x = np.array(x)
    except:
        raise ValueError('Impossible to convert x to a numpy array.')
    # eigen values and eigen vectors of data covariance matrix
    eigen_values, eigen_vectors = np.linalg.eig(np.cov(x.T))
    # sort eigen vectors according biggest eigen value
    eigen_order = eigen_vectors.T[(-eigen_values).argsort()]
    # form output - order the eigenvalues
    return eigen_values[(-eigen_values).argsort()]
1046122
def noise2(self, x, y):
		
		# Skew input space to determine which simplex (triangle) we are in
		s = (x + y) * _F2
		i = floor(x + s)
		j = floor(y + s)
		t = (i + j) * _G2
		x0 = x - (i - t) # "Unskewed" distances from cell origin
		y0 = y - (j - t)

		if x0 > y0:
			i1 = 1; j1 = 0 # Lower triangle, XY order: (0,0)->(1,0)->(1,1)
		else:
			i1 = 0; j1 = 1 # Upper triangle, YX order: (0,0)->(0,1)->(1,1)
		
		x1 = x0 - i1 + _G2 # Offsets for middle corner in (x,y) unskewed coords
		y1 = y0 - j1 + _G2
		x2 = x0 + _G2 * 2.0 - 1.0 # Offsets for last corner in (x,y) unskewed coords
		y2 = y0 + _G2 * 2.0 - 1.0

		# Determine hashed gradient indices of the three simplex corners
		perm = self.permutation
		ii = int(i) % self.period
		jj = int(j) % self.period
		gi0 = perm[ii + perm[jj]] % 12
		gi1 = perm[ii + i1 + perm[jj + j1]] % 12
		gi2 = perm[ii + 1 + perm[jj + 1]] % 12

		# Calculate the contribution from the three corners
		tt = 0.5 - x0**2 - y0**2
		if tt > 0:
			g = _GRAD3[gi0]
			noise = tt**4 * (g[0] * x0 + g[1] * y0)
		else:
			noise = 0.0
		
		tt = 0.5 - x1**2 - y1**2
		if tt > 0:
			g = _GRAD3[gi1]
			noise += tt**4 * (g[0] * x1 + g[1] * y1)
		
		tt = 0.5 - x2**2 - y2**2
		if tt > 0:
			g = _GRAD3[gi2]
			noise += tt**4 * (g[0] * x2 + g[1] * y2)

		return noise * 70.0
54959
def susvd(x, x_obs, rho, penalties):
    

    assert type(x) == Tensor, "Input array must be a Tensor"

    while True:

        # proximal operator for the Fro. norm
        x = squared_error(x, rho, x_obs)

        # sequential singular value thresholding
        for ix, penalty in enumerate(penalties):
            x = x.unfold(ix).svt(penalty / rho).fold()

        yield x
578443
def error_perturbation_var(C, S):
    

    """
    K = S.shape[0]
    cov = tmatrix_cov(C)
    for i in range(K):
        R = S[i, :, :]
        X[i] = (R[:, :, np.newaxis] * cov * R[:, np.newaxis, :]).sum()
1185273
def spectral_centroid_and_spread(data, fs):
    """"""
    data = np.mean(data, axis=1)

    nFFT = len(data) // 2
    X = FFT(data, nFFT)

    ind = (np.arange(1, len(X) + 1)) * (fs / (2.0 * len(X)))

    Xt = X.copy()
    Xt = Xt / Xt.max()
    NUM = np.sum(ind * Xt)
    DEN = np.sum(Xt) + EPSILON

    # Centroid:
    C = (NUM / DEN)

    # Spread:
    S = np.sqrt(np.sum(((ind - C) ** 2) * Xt) / DEN)

    # Normalize:
    C = C / (fs / 2.0)
    S = S / (fs / 2.0)

    return (C, S)
114810
def pca_ellipse(data, loc=None, ax=None, **ellipse_kwargs):
  '''
  '''
  from sklearn.decomposition import PCA  # Late import; sklearn is optional
  pca = PCA(n_components=2).fit(data)
  if loc is None:
    loc = pca.mean_
  if ax is None:
    ax = plt.gca()
  cov = pca.explained_variance_ * pca.components_.T
  u,s,v = np.linalg.svd(cov)
  width,height = 2*np.sqrt(s[:2])
  angle = np.rad2deg(np.arctan2(u[1,0], u[0,0]))
  ell = Ellipse(xy=loc, width=width, height=height, angle=angle,
                **ellipse_kwargs)
  ax.add_patch(ell)
  return ell