Internal method used to compute the expansion factor(s) for a collection
    of bounding ellipsoids using bootstrapping.
338688
def _ellipsoids_bootstrap_expand(args):
    """"""

    # Unzipping.
    points, pointvol, vol_dec, vol_check = args
    rstate = np.random

    # Resampling.
    npoints, ndim = points.shape
    idxs = rstate.randint(npoints, size=npoints)  # resample
    idx_in = np.unique(idxs)  # selected objects
    sel = np.ones(npoints, dtype='bool')
    sel[idx_in] = False
    idx_out = np.where(sel)[0]  # "missing" objects
    if len(idx_out) < 2:  # edge case
        idx_out = np.append(idx_out, [0, 1])
    points_in, points_out = points[idx_in], points[idx_out]

    # Compute bounding ellipsoids.
    ell = bounding_ellipsoid(points_in, pointvol=pointvol)
    ells = _bounding_ellipsoids(points_in, ell, pointvol=pointvol,
                                vol_dec=vol_dec, vol_check=vol_check)

    # Compute normalized distances to missing points.
    dists = [min([el.distance(p) for el in ells]) for p in points_out]

    # Compute expansion factor.
    expand = max(1., max(dists))

    return expand
1325181
private void calcInterfacesCrystal(StructureInterfaceList set, double cutoff) {


		// initialising debugging vars
		long start = -1;
		long end = -1;
		int trialCount = 0;
		int skippedRedundant = 0;
		int skippedAUsNoOverlap = 0;
		int skippedChainsNoOverlap = 0;
		int skippedSelfEquivalent = 0;

		// The bounding boxes of all AUs of the unit cell
		UnitCellBoundingBox bbGrid = new UnitCellBoundingBox(numOperatorsSg, numPolyChainsAu);;
		// we calculate all the bounds of each of the asym units, those will then be reused and translated
		bbGrid.setBbs(structure, ops, INCLUDE_HETATOMS);


		// if not crystallographic there's no search to do in other cells, only chains within "AU" will be checked
		if (!searchBeyondAU) numCells = 0;

		boolean verbose = logger.isDebugEnabled();

		if (verbose) {
			trialCount = 0;
			start= System.currentTimeMillis();
			int neighbors = (2*numCells+1)*(2*numCells+1)*(2*numCells+1)-1;
			int auTrials = (numPolyChainsAu*(numPolyChainsAu-1))/2;
			int trials = numPolyChainsAu*numOperatorsSg*numPolyChainsAu*neighbors;
			logger.debug("Chain clash trials within original AU: "+auTrials);
			logger.debug(
					"Chain clash trials between the original AU and the neighbouring "+neighbors+
					" whole unit cells ("+numCells+" neighbours)" +
					"(2x"+numPolyChainsAu+"chains x "+numOperatorsSg+"AUs x "+neighbors+"cells) : "+trials);
			logger.debug("Total trials: "+(auTrials+trials));
		}

		List<Chain> polyChains = structure.getPolyChains();

		for (int a=-numCells;a<=numCells;a++) {
			for (int b=-numCells;b<=numCells;b++) {
				for (int c=-numCells;c<=numCells;c++) {

					Point3i trans = new Point3i(a,b,c);
					Vector3d transOrth = new Vector3d(a,b,c);
					if (a!=0 || b!=0 || c!=0) {
						// we avoid doing the transformation for 0,0,0 (in case it's not crystallographic)
						this.crystallographicInfo.getCrystalCell().transfToOrthonormal(transOrth);
					}

					UnitCellBoundingBox bbGridTrans = bbGrid.getTranslatedBbs(transOrth);

					for (int n=0;n<numOperatorsSg;n++) {

						// short-cut strategies
						// 1) we skip first of all if the bounding boxes of the AUs don't overlap
						if (!bbGrid.getAuBoundingBox(0).overlaps(bbGridTrans.getAuBoundingBox(n), cutoff)) {
							skippedAUsNoOverlap++;
							continue;
						}

						// 2) we check if we didn't already see its equivalent symmetry operator partner
						CrystalTransform tt = new CrystalTransform(this.crystallographicInfo.getSpaceGroup(), n);
						tt.translate(trans);
						if (isRedundantTransform(tt)) {
							skippedRedundant++;
							continue;
						}
						addVisitedTransform(tt);


						boolean selfEquivalent = false;

						// 3) an operator can be "self redundant" if it is the inverse of itself (involutory, e.g. all pure 2-folds with no translation)
						if (tt.isEquivalent(tt)) {
							logger.debug("Transform {} is equivalent to itself, will skip half of i-chains to j-chains comparisons", tt.toString());
							// in this case we can't skip the operator, but we can skip half of the matrix comparisons e.g. j>i
							// we set a flag and do that within the loop below
							selfEquivalent = true;
						}

						StringBuilder builder = null;
						if (verbose) builder = new StringBuilder(String.valueOf(tt)).append(" ");

						// Now that we know that boxes overlap and operator is not redundant, we have to go to the details
						int contactsFound = 0;

						for (int j=0;j<numPolyChainsAu;j++) {

							for (int i=0;i<numPolyChainsAu;i++) { // we only have to compare the original asymmetric unit to every full cell around

								if(selfEquivalent && (j>i)) {
									// in case of self equivalency of the operator we can safely skip half of the matrix
									skippedSelfEquivalent++;
									continue;
								}
								// special case of original AU, we don't compare a chain to itself
								if (n==0 && a==0 && b==0 && c==0 && i==j) continue;

								// before calculating the AtomContactSet we check for overlap, then we save putting atoms into the grid
								if (!bbGrid.getChainBoundingBox(0,i).overlaps(bbGridTrans.getChainBoundingBox(n,j),cutoff)) {
									skippedChainsNoOverlap++;
									if (verbose) {
										builder.append(".");
									}
									continue;
								}

								trialCount++;

								// finally we've gone through all short-cuts and the 2 chains seem to be close enough:
								// we do the calculation of contacts
								Chain chaini = polyChains.get(i);
								Chain chainj = polyChains.get(j);

								if (n!=0 || a!=0 || b!=0 || c!=0) {
									Matrix4d mJCryst = new Matrix4d(ops[n]);
									translate(mJCryst, transOrth);
									chainj = (Chain)chainj.clone();
									Calc.transform(chainj,mJCryst);
								}

								StructureInterface interf = calcContacts(chaini, chainj, cutoff, tt, builder);
								if (interf == null) {
									continue;
								}

								contactsFound++;
								if(this.hasNcsOps()) {
									StructureInterface interfNcsRef = findNcsRef(interf);
									set.addNcsEquivalent(interf,interfNcsRef);
								} else {
									set.add(interf);
								}
							}
						}

						if( verbose ) {
							if (a==0 && b==0 && c==0 && n==0)
								builder.append(" "+contactsFound+"("+(numPolyChainsAu*(numPolyChainsAu-1))/2+")");
							else if (selfEquivalent)
								builder.append(" "+contactsFound+"("+(numPolyChainsAu*(numPolyChainsAu+1))/2+")");
							else
								builder.append(" "+contactsFound+"("+numPolyChainsAu*numPolyChainsAu+")");

							logger.debug(builder.toString());
						}
					}
				}
			}
		}

		end = System.currentTimeMillis();
		logger.debug("\n"+trialCount+" chain-chain clash trials done. Time "+(end-start)/1000+"s");
		logger.debug("  skipped (not overlapping AUs)       : "+skippedAUsNoOverlap);
		logger.debug("  skipped (not overlapping chains)    : "+skippedChainsNoOverlap);
		logger.debug("  skipped (sym redundant op pairs)    : "+skippedRedundant);
		logger.debug("  skipped (sym redundant self op)     : "+skippedSelfEquivalent);
		logger.debug("Found "+set.size()+" interfaces.");
	}
1130196
private static int[] computeOrigins(int[] sizes, int offset) {
        int count = sizes.length;
        int[] origins = new int[count + 1];
        origins[0] = offset;
        for (int i = 1; i <= count; i++) {
            origins[i] = origins[i - 1] + sizes[i - 1];
        }
        return origins;
    }
156222
def bootstrap_noise(data, func, n=10000, std=1, symmetric=True):
    
    boot_dist = []
    arr = N.zeros(data.shape)
    for i in range(n):
        if symmetric:
            # Noise on all three axes
            arr = N.random.randn(*data.shape)*std
        else:
            # Only z-coordinate noise
            arr[:,-1] = N.random.randn(data.shape[0])*std
        boot_dist.append(func(data+arr))
    return N.array(boot_dist)
1105900
void computeNodeInfo(List<EllipseRotated_F64> ellipses , List<Node> cluster ) {

		// create an info object for each member inside of the cluster
		listInfo.reset();
		for (int i = 0; i < cluster.size(); i++) {
			Node n = cluster.get(i);
			EllipseRotated_F64 t = ellipses.get( n.which );

			NodeInfo info = listInfo.grow();
			info.reset();
			info.ellipse = t;
		}

		addEdgesToInfo(cluster);
		pruneNearlyIdenticalAngles();
		findLargestAnglesForAllNodes();
	}
626521
def denormalize_bbox(bbox, rows, cols):
    
    if rows == 0:
        raise ValueError('Argument rows cannot be zero')
    if cols == 0:
        raise ValueError('Argument cols cannot be zero')

    x_min, y_min, x_max, y_max = bbox[:4]
    denormalized_bbox = [x_min * cols, y_min * rows, x_max * cols, y_max * rows]
    return denormalized_bbox + list(bbox[4:])
1109137
def geom_symm_match(g, atwts, ax, theta, do_refl):
    

    # Imports
    import numpy as np
    from scipy import linalg as spla

    # Convert g and atwts to n-D vectors
    g = make_nd_vec(g, nd=None, t=np.float64, norm=False)
    atwts = make_nd_vec(atwts, nd=None, t=np.float64, norm=False)

    # Ensure proper dimensionality
    if not g.shape[0] == 3 * atwts.shape[0]:
        raise ValueError("Size of 'g' is not 3*size of 'atwts'")
    ## end if

    # Calculate transformed geometry
    gx = symm_op(g, ax, theta, do_refl)

    # Push g to a column vector
    g = g.reshape((g.shape[0],1))

    # Augment g and gx with imaginary atomic weights
    ex_wts = atwts.repeat(3,axis=0).T.reshape((atwts.shape[0]*3,1)) * 1.j
    g = np.add(g, ex_wts)
    gx = np.add(gx, ex_wts)

##        # Define calc as the outer product of the augmented vectors
##        calc = np.dot(g.reshape((g.shape[0],1)), \
##                            np.reciprocal(gx.reshape((1,gx.shape[0]))))
##
##        # Calculate the complex magnitude of each element and take log10,
##        #  then abs again
##        calc = np.abs(np.log10(np.abs(calc)))

    # Expand g and gx as column vectors of coordinates
    calc_g = g.reshape((g.shape[0] // 3, 3))
    calc_gx = gx.reshape((gx.shape[0] // 3, 3))
##
##        # Expand each into a square matrix of identical column vectors
##        calc_g = calc_g.repeat(g.shape[0], axis=1)
##        calc_gx = gx.repeat(gx.shape[0], axis=1)

    # Calc is the absolute distance between the calc-ed values,
    #  scaled by the maximum of the individual atom distances or unity.

    # Calculate the unscaled distances
    calc = [[spla.norm(np.subtract(calc_g[i,:], calc_gx[j,:])) \
                            for j in range(calc_gx.shape[0])] \
                            for i in range(calc_g.shape[0])]

    # Calculate the scale factors
    scale_g = np.array([spla.norm(calc_g[i,:]) for i in \
                range(calc_g.shape[0])]).reshape((calc_g.shape[0],1)) \
                .repeat(calc_g.shape[0], axis=1)
    scale_gx = np.array([spla.norm(calc_gx[j,:]) for j in \
                range(calc_g.shape[0])]).reshape((1,calc_gx.shape[0])) \
                .repeat(calc_gx.shape[0], axis=0)
    scale = np.maximum(np.maximum(scale_g, scale_gx),
                np.ones_like(scale_g, dtype=np.float64))

    # Scale calc
    calc = np.divide(calc, scale)

    # Take the minimum of each row
    mins = np.min(calc, axis=1)

    # Take the maximum of the minima for the final factor
    fac = np.max(mins)

    # Using the atomic weights for checking matching can result in 'fac'
    #  being greater than unity. Return the minimum of fac and unity.
    fac = min(fac, 1.0)
    return fac
1199606
def cylinders_from_ellipses(ellipses):
    """"""
    ellipses = np.asarray(ellipses)
    ellipsoids = np.zeros((ellipses.shape[0], 10))
    ellipsoids[:, [0, 1, 2, 4, 5, 7]] = ellipses
    ellipsoids[:, 3] = 100000.0

    return ellipsoids
132311
def expand_grid_coords(dataset, dim):
    
    arrays = [dataset.interface.coords(dataset, d.name, True)
              for d in dataset.kdims]
    idx = dataset.get_dimension_index(dim)
    return cartesian_product(arrays, flat=False)[idx].T
850666
def calculate_bounding_box(data):
    
    mins = data.min(0)
    maxes = data.max(0)
    return mins, maxes
933873
private void calcInnerBounds() {
        final java.awt.Insets INSETS = getInsets();
        INNER_BOUNDS.setBounds(INSETS.left, INSETS.top, (getWidth() - INSETS.left - INSETS.right), (getHeight() - INSETS.top - INSETS.bottom));
    }
411558
def _fit(
    indexed_params,
    penalization,
    lam,
    lam_perturb,
    lam_scale_,
    estimator,
    penalty_name,
    subsample,
    bootstrap,
    prng,
    X=None,
):
    
    index = indexed_params

    if isinstance(X, np.ndarray):
        local_X = X
    else:
        local_X = X.value

    n_samples, n_features = local_X.shape

    prec_is_real = False
    while not prec_is_real:
        boot_lam = None
        if penalization == "subsampling":
            pass
        elif penalization == "random":
            boot_lam = _fix_weights(_random_weights, n_features, lam, lam_perturb, prng)
        elif penalization == "fully-random":
            boot_lam = _fix_weights(_fully_random_weights, n_features, lam_scale_, prng)
        else:
            raise NotImplementedError(
                (
                    "Only penalization = 'subsampling', "
                    "'random', and 'fully-random' have "
                    "been implemented. Found {}.".format(penalization)
                )
            )

        # new instance of estimator
        new_estimator = clone(estimator)
        if boot_lam is not None:
            new_estimator.set_params(**{penalty_name: boot_lam})

        # fit estimator
        num_subsamples = int(subsample * n_samples)
        rp = bootstrap(n_samples, num_subsamples, prng)
        new_estimator.fit(local_X[rp, :])

        # check that new_estimator.precision_ is real
        # if not, skip this boot_lam and try again
        if isinstance(new_estimator.precision_, list):
            prec_real_bools = []
            for prec in new_estimator.precision_:
                prec_real_bools.append(np.all(np.isreal(prec)))

            prec_is_real = np.all(np.array(prec_real_bools) is True)

        elif isinstance(new_estimator.precision_, np.ndarray):
            prec_is_real = np.all(np.isreal(new_estimator.precision_))

        else:
            raise ValueError("Estimator returned invalid precision_.")

    return index, (boot_lam, rp, new_estimator)
1012078
def _modified_shepp_logan_ellipsoids(ellipsoids):
    
    intensities = [1.0, -0.8, -0.2, -0.2, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]

    # Add minimal numbers to ensure that the result is nowhere negative.
    # This is needed due to numerical issues.
    intensities[2] += 5e-17
    intensities[3] += 5e-17

    assert len(ellipsoids) == len(intensities)

    for ellipsoid, intensity in zip(ellipsoids, intensities):
        ellipsoid[0] = intensity
674639
public void computeBasis( int numComponents ) {
        if( numComponents > A.getNumCols() )
            throw new IllegalArgumentException("More components requested that the data's length.");
        if( sampleIndex != A.getNumRows() )
            throw new IllegalArgumentException("Not all the data has been added");
        if( numComponents > sampleIndex )
            throw new IllegalArgumentException("More data needed to compute the desired number of components");

        this.numComponents = numComponents;

        // compute the mean of all the samples
        for( int i = 0; i < A.getNumRows(); i++ ) {
            for( int j = 0; j < mean.length; j++ ) {
                mean[j] += A.get(i,j);
            }
        }
        for( int j = 0; j < mean.length; j++ ) {
            mean[j] /= A.getNumRows();
        }

        // subtract the mean from the original data
        for( int i = 0; i < A.getNumRows(); i++ ) {
            for( int j = 0; j < mean.length; j++ ) {
                A.set(i,j,A.get(i,j)-mean[j]);
            }
        }

        // Compute SVD and save time by not computing U
        SingularValueDecomposition<DMatrixRMaj> svd =
                DecompositionFactory_DDRM.svd(A.numRows, A.numCols, false, true, false);
        if( !svd.decompose(A) )
            throw new RuntimeException("SVD failed");

        V_t = svd.getV(null,true);
        DMatrixRMaj W = svd.getW(null);

        // Singular values are in an arbitrary order initially
        SingularOps_DDRM.descendingOrder(null,false,W,V_t,true);

        // strip off unneeded components and find the basis
        V_t.reshape(numComponents,mean.length,true);
    }
460452
def denormalize_bboxes(bboxes, rows, cols):
    """"""
    return [denormalize_bbox(bbox, rows, cols) for bbox in bboxes]
1300054
def _compute_cell_extents_grid(bounding_rect=(0.03, 0.03, 0.97, 0.97),
                                   num_rows=2, num_cols=6,
                                   axis_pad=0.01):
        

        left, bottom, width, height = bounding_rect
        height_padding = axis_pad * (num_rows + 1)
        width_padding = axis_pad * (num_cols + 1)
        cell_height = float((height - height_padding) / num_rows)
        cell_width = float((width - width_padding) / num_cols)

        cell_height_padded = cell_height + axis_pad
        cell_width_padded = cell_width + axis_pad

        extents = list()
        for row in range(num_rows - 1, -1, -1):
            for col in range(num_cols):
                extents.append((left + col * cell_width_padded,
                                bottom + row * cell_height_padded,
                                cell_width, cell_height))

        return extents
515122
def _fit_ellipsoid_partial(locations, cy):
    """"""
    a = np.vstack([locations[:, 0]**2,
                   locations[:, 1]**2 - 2 * locations[:, 1] * cy,
                   locations[:, 2]**2,
                   locations[:, 0]*2,
                   locations[:, 2]*2]).T
    x = sp.linalg.lstsq(a, np.ones(locations.shape[0]))[0]
    c = [-x[3] / x[0], cy, -x[4] / x[2]]
    gam = 1 + x[3]**2 / x[0] + x[4]**2 / x[2]
    r = np.sqrt([gam / x[0], gam / x[1], gam / x[2]])
    return c, r