An iterative version of `pca.plane_errors`,
    which computes an error surface for a plane.
977202
def iterative_plane_errors(axes,covariance_matrix, **kwargs):
    
    sheet = kwargs.pop('sheet','upper')
    level = kwargs.pop('level',1)
    n = kwargs.pop('n',100)

    cov = N.sqrt(N.diagonal(covariance_matrix))
    u = N.linspace(0, 2*N.pi, n)

    scales = dict(upper=1,lower=-1,nominal=0)
    c1 = scales[sheet]
    c1 *= -1 # We assume upper hemisphere
    if axes[2,2] < 0:
        c1 *= -1

    def sdot(a,b):
        return sum([i*j for i,j in zip(a,b)])

    def step_func(a):
        e = [
            N.cos(a)*cov[0],
            N.sin(a)*cov[1],
            c1*cov[2]]
        d = [sdot(e,i)
            for i in axes.T]
        x,y,z = d[2],d[0],d[1]
        r = N.sqrt(x**2 + y**2 + z**2)
        lat = N.arcsin(z/r)
        lon = N.arctan2(y, x)
        return lon,lat

    # Get a bundle of vectors defining
    # a full rotation around the unit circle
    return N.array([step_func(i)
        for i in u])
430819
def _xval_model_error(self, fit_spectra, reject_outliers, fit_lb, fit_ub,
                           fitter, func):
        
        set1 = fit_spectra[::2]
        set2 = fit_spectra[1::2]
        errs = []
        signal_select = [] 
        models = []
        signals = []
        for this_set in [set1, set2]:
            choose_transients, model, signal, params, this_idx =\
                self._fit_helper(this_set, reject_outliers,
                                 fit_lb, fit_ub, fitter)
            models.append(np.nanmean(model[choose_transients], 0))
            signals.append(np.nanmean(signal[choose_transients], 0))

            signal_select.append(signal[choose_transients])
                
        #Cross-validation error estimation:
        model_err = np.mean([ut.rmse(models[0], signals[1]),
                              ut.rmse(models[1], signals[0])])
        # Also for the signal:
        signal_err = ut.rmse(np.nanmean(signal_select[0], 0),
                             np.nanmean(signal_select[1], 0))
        # Based on the errors, choose a function. Also report errors:
        return model_err, signal_err
1164780
def sparse_surface(self):
        
        if self._method == 'ray':
            func = voxelize_ray
        elif self._method == 'subdivide':
            func = voxelize_subdivide
        else:
            raise ValueError('voxelization method incorrect')

        voxels, origin = func(
            mesh=self._data['mesh'],
            pitch=self._data['pitch'],
            max_iter=self._data['max_iter'][0])
        self._cache['origin'] = origin

        return voxels
267874
def pca(adata, **kwargs) -> Union[Axes, List[Axes], None]:
    
    return plot_scatter(adata, 'pca', **kwargs)
710041
def fitPlane(points, c="g", bc="darkgreen"):
    
    data = np.array(points)
    datamean = data.mean(axis=0)
    res = np.linalg.svd(data - datamean)
    dd, vv = res[1], res[2]
    xyz_min = points.min(axis=0)
    xyz_max = points.max(axis=0)
    s = np.linalg.norm(xyz_max - xyz_min)
    n = np.cross(vv[0], vv[1])
    pla = vs.Plane(datamean, n, s, s, c, bc)
    pla.info["normal"] = n
    pla.info["center"] = datamean
    pla.info["variance"] = dd[2]
    return pla
596760
void constructLinearSystem(Planar<I> input , Se3_F32 g ) {
		int numBands = imageType.getNumBands();

		// first precompute everything that does not depend on pixel values
		inboundsPixels = 0;
		for (int i = 0; i < keypixels.size(); i++) {
			Pixel p = keypixels.data[i];

			// Apply the known warp
			SePointOps_F32.transform(g, p.p3, S);

			if( S.z <= 0 ) {
				p.valid = false;
				continue;
			}

			// Compute projected warped pixel coordinate on image I_1
			p.proj.x = (S.x / S.z) * fx + cx;
			p.proj.y = (S.y / S.z) * fy + cy;

			// make sure it's in the bounds
			if (p.proj.x < 0 || p.proj.x > input.width - 1 || p.proj.y < 0 || p.proj.y > input.height - 1) {
				p.valid = false;
				continue;
			} else {
				p.valid = true;
			}
			inboundsPixels++;

			// pi matrix derivative relative to t at S
			float ZZ = S.z * S.z;

			p.dP11 = fx / S.z;
			p.dP13 = -S.x * fx / ZZ;
			p.dP22 = fy / S.z;
			p.dP23 = -S.y * fy / ZZ;
		}

		// how compute the components which require
		errorOptical = 0;
		int row = 0;
		for (int band = 0; band < numBands; band++) {
			interpDX.setImage(derivX.getBand(band));
			interpDY.setImage(derivY.getBand(band));
			interpI.setImage(input.getBand(band));

			for (int i = 0; i < keypixels.size(); i++) {
				Pixel p = keypixels.data[i];

				if( !p.valid )
					continue;

				// Apply the known warp
				// TODO precompute?
				SePointOps_F32.transform(g, p.p3, S);

				// sample pixel values at warped location in I_1
				// NOTE: This could be highly optimized.  Compute and save interpolation weights once per input
				//       instead of for each band and image (current,dx,dy)
				// TODO create a special bilinear class for this?
				float current = interpI.get( p.proj.x, p.proj.y);
				float dx      = interpDX.get(p.proj.x, p.proj.y);
				float dy      = interpDY.get(p.proj.x, p.proj.y);

				// B = grad^T * dPI/dt = shape(1,3)
				float b1 = dx*p.dP11;
				float b2 = dy*p.dP22;
				float b3 = dx*p.dP13 + dy*p.dP23;

				// C * A(S'(x)) = shape(1,6)
				int indexA = row * 6;
				A.data[indexA++] = -b2*S.z + b3*S.y;
				A.data[indexA++] =  b1*S.z - b3*S.x;
				A.data[indexA++] = -b1*S.y + b2*S.x;
				A.data[indexA++] = b1;
				A.data[indexA++] = b2;
				A.data[indexA]   = b3;

				float error = -(current - p.bands[band]);
				y.data[row] = error;

				errorOptical += Math.abs(error);
				row += 1;
			}
		}
		errorOptical /= row;

		A.numRows = row;
		y.numRows = row;
	}
91265
def normal_errors(axes, covariance_matrix, **kwargs):
    
    level = kwargs.pop('level',1)
    traditional_layout = kwargs.pop('traditional_layout',True)
    d = N.diagonal(covariance_matrix)
    ell = ellipse(**kwargs)

    if axes[2,2] < 0:
        axes *= -1

    # Not sure where this factor comes from but it
    # seems to make things work better
    c1 = 2
    axis_lengths = d[:2]
    f = N.linalg.norm(
        ell*axis_lengths,axis=1)

    e0 = -ell.T*d[2]*c1
    e = N.vstack((e0,f))

    _ = dot(e.T,axes).T

    if traditional_layout:
        lon,lat = stereonet_math.cart2sph(_[2],_[0],-_[1])
    else:
        lon,lat = stereonet_math.cart2sph(-_[1],_[0],_[2])
    return list(zip(lon,lat))
1207906
function() {
            var EPS = 1e-5;
            var errors = [];
            if (this.isSelfIntersecting(true)) {
                errors.push("Self intersects");
            }
            var pointcount = {};
            this.sides.map(function(side) {
                function mappoint(p) {
                    var tag = p.x + " " + p.y;
                    if (!(tag in pointcount)) pointcount[tag] = 0;
                    pointcount[tag] ++;
                }
                mappoint(side.vertex0.pos);
                mappoint(side.vertex1.pos);
            });
            for (var tag in pointcount) {
                var count = pointcount[tag];
                if (count & 1) {
                    errors.push("Uneven number of sides (" + count + ") for point " + tag);
                }
            }
            var area = this.area();
            if (area < EPS*EPS) {
                errors.push("Area is " + area);
            }
            if (errors.length > 0) {
                var ertxt = "";
                errors.map(function(err) {
                    ertxt += err + "\n";
                });
                throw new Error(ertxt);
            }
        }
396902
def _embedding_kernel_pca(matrix, dimensions=3, affinity_matrix=None,
                          sigma=1):
    
    if affinity_matrix is None:
        aff = rbf(matrix, sigma)
    else:
        aff = affinity_matrix
    kpca = sklearn.decomposition.KernelPCA(kernel='precomputed',
                                           n_components=dimensions)
    return kpca.fit_transform(aff)
812633
def array_scanlines_interlace(self, pixels):
        

        # http://www.w3.org/TR/PNG/#8InterlaceMethods
        # Array type.
        fmt = 'BH'[self.bitdepth > 8]
        # Value per row
        vpr = self.width * self.planes

        # Each iteration generates a scanline starting at (x, y)
        # and consisting of every xstep pixels.
        for lines in adam7_generate(self.width, self.height):
            for x, y, xstep in lines:
                # Pixels per row (of reduced image)
                ppr = int(math.ceil((self.width - x) / float(xstep)))
                # Values per row (of reduced image)
                reduced_row_len = ppr * self.planes
                if xstep == 1:
                    # Easy case: line is a simple slice.
                    offset = y * vpr
                    yield pixels[offset: offset + vpr]
                    continue
                # We have to step by xstep,
                # which we can do one plane at a time
                # using the step in Python slices.
                row = array(fmt)
                # There's no easier way to set the length of an array
                row.extend(pixels[0:reduced_row_len])
                offset = y * vpr + x * self.planes
                end_offset = (y + 1) * vpr
                skip = self.planes * xstep
                for i in range(self.planes):
                    row[i::self.planes] = \
                        pixels[offset + i: end_offset: skip]
                yield row
1102150
double computeError( List<AssociatedPair> points , Se3_F64 worldToCamera ) {

		double error = 0;

		for (int i = 0; i < points.size(); i++) {
			AssociatedPair pair = points.get(i);
			tmpP.set(pair.p1.x,pair.p1.y,0);
			SePointOps_F64.transform(worldToCamera,tmpP,tmpP);

			error += pair.p2.distance2(tmpP.x/tmpP.z,tmpP.y/tmpP.z);
		}

		return Math.sqrt(error/points.size());
	}
1107051
private void triangulateFeatures(List<AssociatedTriple> inliers,
									 DMatrixRMaj P1, DMatrixRMaj P2, DMatrixRMaj P3) {
		List<DMatrixRMaj> cameraMatrices = new ArrayList<>();
		cameraMatrices.add(P1);
		cameraMatrices.add(P2);
		cameraMatrices.add(P3);

		// need elements to be non-empty so that it can use set().  probably over optimization
		List<Point2D_F64> triangObs = new ArrayList<>();
		triangObs.add(null);
		triangObs.add(null);
		triangObs.add(null);

		Point4D_F64 X = new Point4D_F64();
		for (int i = 0; i < inliers.size(); i++) {
			AssociatedTriple t = inliers.get(i);

			triangObs.set(0,t.p1);
			triangObs.set(1,t.p2);
			triangObs.set(2,t.p3);

			// triangulation can fail if all 3 views have the same pixel value. This has been observed in
			// simulated 3D scenes
			if( triangulator.triangulate(triangObs,cameraMatrices,X)) {
				structure.points[i].set(X.x,X.y,X.z,X.w);
			} else {
				throw new RuntimeException("Failed to triangulate a point in the inlier set?! Handle if this is common");
			}
		}
	}
331303
def compute_absolute_error(self, predicted_data, record, dataframe_record):
        ''''''
        absolute_error = abs(record['DDG'] - predicted_data[self.ddg_analysis_type])
        dataframe_record['AbsoluteError'] = absolute_error
366825
def _photometricErrors(self, n_per_bin=100, plot=False):
        

        self.catalog.spatialBin(self.roi)

        if len(self.catalog.mag_1) < n_per_bin:
            logger.warning("Catalog contains fewer objects than requested to calculate errors.")
            n_per_bin = int(len(self.catalog.mag_1) / 3)

        # Band 1
        mag_1_thresh = self.mask.mask_1.mask_roi_sparse[self.catalog.pixel_roi_index] - self.catalog.mag_1
        sorting_indices = np.argsort(mag_1_thresh)
        mag_1_thresh_sort = mag_1_thresh[sorting_indices]
        mag_err_1_sort = self.catalog.mag_err_1[sorting_indices]

        # ADW: Can't this be done with np.median(axis=?)
        mag_1_thresh_medians = []
        mag_err_1_medians = []
        for i in range(0, int(len(mag_1_thresh) / float(n_per_bin))):
            mag_1_thresh_medians.append(np.median(mag_1_thresh_sort[n_per_bin * i: n_per_bin * (i + 1)]))
            mag_err_1_medians.append(np.median(mag_err_1_sort[n_per_bin * i: n_per_bin * (i + 1)]))
        
        if mag_1_thresh_medians[0] > 0.:
            mag_1_thresh_medians = np.insert(mag_1_thresh_medians, 0, -99.)
            mag_err_1_medians = np.insert(mag_err_1_medians, 0, mag_err_1_medians[0])
        
        self.photo_err_1 = scipy.interpolate.interp1d(mag_1_thresh_medians, mag_err_1_medians,
                                                      bounds_error=False, fill_value=mag_err_1_medians[-1])

        # Band 2
        mag_2_thresh = self.mask.mask_2.mask_roi_sparse[self.catalog.pixel_roi_index] - self.catalog.mag_2
        sorting_indices = np.argsort(mag_2_thresh)
        mag_2_thresh_sort = mag_2_thresh[sorting_indices]
        mag_err_2_sort = self.catalog.mag_err_2[sorting_indices]

        mag_2_thresh_medians = []
        mag_err_2_medians = []
        for i in range(0, int(len(mag_2_thresh) / float(n_per_bin))):
            mag_2_thresh_medians.append(np.median(mag_2_thresh_sort[n_per_bin * i: n_per_bin * (i + 1)]))
            mag_err_2_medians.append(np.median(mag_err_2_sort[n_per_bin * i: n_per_bin * (i + 1)]))

        if mag_2_thresh_medians[0] > 0.:
            mag_2_thresh_medians = np.insert(mag_2_thresh_medians, 0, -99.)
            mag_err_2_medians = np.insert(mag_err_2_medians, 0, mag_err_2_medians[0])
        
        self.photo_err_2 = scipy.interpolate.interp1d(mag_2_thresh_medians, mag_err_2_medians,
                                                      bounds_error=False, fill_value=mag_err_2_medians[-1])
695968
def PCA(Y, components):
	

	N,D = Y.shape
	sv = linalg.svd(Y, full_matrices=0);
	[s0, w0] = [sv[0][:, 0:components], np.dot(np.diag(sv[1]), sv[2]).T[:, 0:components]]
	v = s0.std(axis=0)
	s0 /= v;
	w0 *= v;
	return [s0, w0]

	if N>D:
		sv = linalg.svd(Y, full_matrices=0);
		[s0, w0] = [sv[0][:, 0:components], np.dot(np.diag(sv[1]), sv[2]).T[:, 0:components]]
		v = s0.std(axis=0)
		s0 /= v;
		w0 *= v;
		return [s0, w0]
	else:
		K=np.cov(Y)
		sv = linalg.eigh(K)
		std_var = np.sqrt(sv[0])
		pc = sv[1]*std_var[np.newaxis(),0]
		#import ipdb
		#ipdb.set_trace()
		return [pc,std_var]
1352048
def getTraitCovarStdErrors(self,term_i):
        
        assert self.init,        'GP not initialised'
        assert self.fast==False, 'Not supported for fast implementation'

        if self.P==1:
            out = (2*self.getScales()[term_i])**2*self._getLaplaceCovar()[term_i,term_i]
        else:
            C = self.vd.getTerm(term_i).getTraitCovar()
            n_params = C.getNumberParams()
            par_index = 0
            for term in range(term_i-1):
                par_index += self.vd.getTerm(term_i).getNumberScales()
            Sigma1 = self._getLaplaceCovar()[par_index:(par_index+n_params),:][:,par_index:(par_index+n_params)]
            out = sp.zeros((self.P,self.P))
            for param_i in range(n_params):
                out += C.Kgrad_param(param_i)**2*Sigma1[param_i,param_i]
                for param_j in range(param_i):
                    out += 2*abs(C.Kgrad_param(param_i)*C.Kgrad_param(param_j))*Sigma1[param_i,param_j]
        out = sp.sqrt(out)
        return out
1097575
public Se3_F64 estimateNoOutliers( List<Point2D3D> observations ) {

		// Compute a single solution using EPNP
		// 10 iterations is what JavaDoc recommends, but might need to be tuned.
		// 0 test points.  This parameters is actually ignored because EPNP only returns a single solution
		Estimate1ofPnP pnp = FactoryMultiView.pnp_1(EnumPNP.EPNP, 10, 0);

		Se3_F64 worldToCamera = new Se3_F64();
		pnp.process(observations,worldToCamera);

		// For some applications the EPNP solution might be good enough, but let's refine it
		RefinePnP refine = FactoryMultiView.pnpRefine(1e-8,200);

		Se3_F64 refinedWorldToCamera = new Se3_F64();

		if( !refine.fitModel(observations,worldToCamera,refinedWorldToCamera) )
			throw new RuntimeException("Refined failed! Input probably bad...");

		return refinedWorldToCamera;
	}