Helper method to sample from the lateral surface of a cylinder.
349949
def _sampleLocationOnSide(self):
    
    z = random.uniform(-1, 1) * self.height / 2.
    sampledAngle = 2 * random.random() * pi
    x, y = self.radius * cos(sampledAngle), self.radius * sin(sampledAngle)
    return [x, y, z]
1210857
def samples(self, gp, Y_metadata=None):
        
        orig_shape = gp.shape
        gp = gp.flatten()
        # Ysim = np.random.poisson(self.gp_link.transf(gp), [samples, gp.size]).T
        # return Ysim.reshape(orig_shape+(samples,))
        Ysim = np.random.poisson(self.gp_link.transf(gp))
        return Ysim.reshape(orig_shape)
1192538
def sample_within_cc(self, cc_index, nsamples=1):
        
        polygon = self.geometries[cc_index]['polygon']
        samples = []
        while len(samples) < nsamples:
            point = PointSampler.random_point(polygon.envelope.bounds)
            if PointSampler.contains(polygon, point):
                samples.append(point)
        return samples
186505
def find_average_diameter(elements, coordinates, adjust=1, increment=0.1,
                          processes=None, **kwargs):
    """"""
    # Copy the coordinates as will perform many opertaions on them
    coordinates = deepcopy(coordinates)
    # Center of our cartesian system is always at origin
    origin = np.array([0, 0, 0])
    # Initial center of mass to reverse translation at the end
    initial_com = center_of_mass(elements, coordinates)
    # We just shift the cage to the origin.
    coordinates = shift_com(elements, coordinates)
    # We create an array of vdw radii of elements.
    elements_vdw = np.array([[atomic_vdw_radius[x.upper()]] for x in elements])
    # We calculate maximum diameter of a molecule to determine the radius
    # of a sampling sphere neccessary to enclose the whole molecule.
    shpere_radius = max_dim(elements, coordinates)[2]
    sphere_surface_area = 4 * np.pi * shpere_radius**2
    # Here we determine the number of sampling points necessary for a fine
    # sampling. Smaller molecules require more finner density of sampling
    # points on the sampling sphere's surface, whereas largen require less.
    # This formula was created so that larger molecule do not take much longer
    # to analyse, as number_sampling_points*length_of_sampling_vectors
    # results in quadratic increase of sampling time. The 250 factor was
    # specificly determined to produce close to 1 sampling point /Angstrom^2
    # for a sphere of radius ~ 24 Angstrom. We can adjust how fine is the
    # sampling by changing the adjust factor.
    number_of_points = int(np.log10(sphere_surface_area) * 250 * adjust)
    # Here I use code by Alexandre Devert for spreading points on a sphere:
    # http://blog.marmakoide.org/?p=1
    golden_angle = np.pi * (3 - np.sqrt(5))
    theta = golden_angle * np.arange(number_of_points)
    z = np.linspace(1 - 1.0 / number_of_points, 1.0 / number_of_points - 1.0,
                    number_of_points)
    radius = np.sqrt(1 - z * z)
    points = np.zeros((number_of_points, 3))
    points[:, 0] = radius * np.cos(theta) * shpere_radius
    points[:, 1] = radius * np.sin(theta) * shpere_radius
    points[:, 2] = z * shpere_radius
    # Here we analyse the vectors and retain the ones that create the molecule
    # outline.
    if processes:
        pool = Pool(processes=processes)
        parallel = [
            pool.apply_async(
                vector_analysis_reversed,
                args=(
                    point, coordinates, elements_vdw)
                ) for point in points
        ]
        results = [p.get() for p in parallel if p.get() is not None]
        pool.terminate()
    else:
        results = [
            vector_analysis_reversed(
                point, coordinates, elements_vdw)
            for point in points
        ]
    results_cleaned = [x[0] for x in results if x is not None]
    return np.mean(results_cleaned)*2
1223617
def _sampleFromVertices(self):
    
    coordinates = [
      self.dimensions[0] / 2. * random.choice([-1, 1]),
      self.dimensions[1] / 2. * random.choice([-1, 1]),
      self.dimensions[2] / 2. * random.choice([-1, 1]),
    ]
    return coordinates
382513
def bounding_cylinder(self):
        
        from . import primitives, bounds
        kwargs = bounds.minimum_cylinder(self)
        mincyl = primitives.Cylinder(mutable=False, **kwargs)
        return mincyl
143737
def _indiameter_from_fibers(self):
        
        """
        Np = self.num_pores()
        indiam = np.zeros(Np, dtype=float)
        incen = np.zeros([Np, 3], dtype=float)
        hull_pores = np.unique(self._hull_image)
        (Lx, Ly, Lz) = np.shape(self._hull_image)
        (indx, indy, indz) = np.indices([Lx, Ly, Lz])
        indx = indx.flatten()
        indy = indy.flatten()
        indz = indz.flatten()
        for i, pore in enumerate(hull_pores):
            logger.info("Processing pore: "+str(i)+" of "+str(len(hull_pores)))
            dt_pore = self._dt_image*(self._hull_image == pore)
            indiam[pore] = dt_pore.max()*2
            max_ind = np.argmax(dt_pore)
            incen[pore, 0] = indx[max_ind]
            incen[pore, 1] = indy[max_ind]
            incen[pore, 2] = indz[max_ind]
        indiam *= self.network.resolution
        incen *= self.network.resolution
        return (indiam, incen)
1302134
def granulometry_filter(image, min_radius, max_radius, mask=None):
    '''
    '''
    #
    # Do 4-connected erosion
    #
    se = np.array([[False, True, False],
                   [True, True, True],
                   [False, True, False]])
    #
    # Initialize
    #
    inverted_image = image.max() - image
    previous_opened_image = image
    eroded_image = image
    selected_granules_image = np.zeros(image.shape)
    #
    # Select granules by successive morphological openings
    #
    for i in range(max_radius+1):
        eroded_image = grey_erosion(eroded_image, mask=mask, footprint = se)
        opened_image = grey_dilation(eroded_image, inverted_image, footprint = se)
        output_image = previous_opened_image - opened_image
        if i >= min_radius:
            selected_granules_image = np.maximum(selected_granules_image, output_image)
        previous_opened_image = opened_image
    return selected_granules_image
325784
def sampleset(self):
        """"""
        x1, x4 = self.bounding_box
        dw = self.width * 0.5
        x2 = self.x_0 - dw
        x3 = self.x_0 + dw

        if self._n_models == 1:
            w = [x1, x2, x3, x4]
        else:
            w = list(zip(x1, x2, x3, x4))

        return np.asarray(w)
430190
def reconstruct_from_shape(self, shape, optimize=False):
        
        axes_lengths = list(self.elementary_axes_lengths)
        if self.ellipsis_positions != (math.inf, math.inf):
            if len(shape) < len(self.input_composite_axes) - 1:
                raise EinopsError('Expected at least {} dimensions, got {}'.format(
                    len(self.input_composite_axes) - 1, len(shape)))
        else:
            if len(shape) != len(self.input_composite_axes):
                raise EinopsError('Expected {} dimensions, got {}'.format(len(self.input_composite_axes), len(shape)))
        for input_axis, (known_axes, unknown_axes) in enumerate(self.input_composite_axes):
            before_ellipsis = input_axis
            after_ellipsis = input_axis + len(shape) - len(self.input_composite_axes)
            if input_axis == self.ellipsis_positions[0]:
                assert len(known_axes) == 0 and len(unknown_axes) == 1
                unknown_axis, = unknown_axes
                ellipsis_shape = shape[before_ellipsis:after_ellipsis + 1]
                if any(d is None for d in ellipsis_shape):
                    raise EinopsError("Couldn't infer shape for one or more axes represented by ellipsis")
                axes_lengths[unknown_axis] = _product(ellipsis_shape)
            else:
                if input_axis < self.ellipsis_positions[0]:
                    length = shape[before_ellipsis]
                else:
                    length = shape[after_ellipsis]
                known_product = 1
                for axis in known_axes:
                    known_product *= axes_lengths[axis]

                if len(unknown_axes) == 0:
                    if isinstance(length, int) and isinstance(known_product, int) and length != known_product:
                        raise EinopsError('Shape mismatch, {} != {}'.format(length, known_product))
                else:
                    if isinstance(length, int) and isinstance(known_product, int) and length % known_product != 0:
                        raise EinopsError("Shape mismatch, can't divide axis of length {} in chunks of {}".format(
                            length, known_product))

                    unknown_axis, = unknown_axes
                    axes_lengths[unknown_axis] = length // known_product

        init_shapes = axes_lengths
        reduced_axes_lengths = [dim for i, dim in enumerate(axes_lengths) if i not in self.reduced_elementary_axes]
        final_shapes = []
        for output_axis, grouping in enumerate(self.output_composite_axes):
            if output_axis == self.ellipsis_positions[1]:
                final_shapes.extend(ellipsis_shape)
            else:
                lengths = [reduced_axes_lengths[elementary_axis] for elementary_axis in grouping]
                if any(l is None for l in lengths):
                    final_shapes.append(None)
                else:
                    final_shapes.append(_product(lengths))
        reduced_axes = self.reduced_elementary_axes
        axes_reordering = self.final_axes_grouping_flat
        if optimize:
            return _optimize_transformation(init_shapes, reduced_axes, axes_reordering, final_shapes)
        else:
            return init_shapes, reduced_axes, axes_reordering, final_shapes
360890
def sample_block_tridiag(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))

    y = info_sample(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 y
296299
def _sample_aAt(self,n):
        """"""
        #Sample frequency along largest eigenvalue using ARS
        dO1s=\
            bovy_ars.bovy_ars([0.,0.],[True,False],
                              [self._meandO-numpy.sqrt(self._sortedSigOEig[2]),
                               self._meandO+numpy.sqrt(self._sortedSigOEig[2])],
                              _h_ars,_hp_ars,nsamples=n,
                              hxparams=(self._meandO,self._sortedSigOEig[2]),
                              maxn=100)
        dO1s= numpy.array(dO1s)*self._sigMeanSign
        dO2s= numpy.random.normal(size=n)*numpy.sqrt(self._sortedSigOEig[1])
        dO3s= numpy.random.normal(size=n)*numpy.sqrt(self._sortedSigOEig[0])
        #Rotate into dOs in R,phi,z coordinates
        dO= numpy.vstack((dO3s,dO2s,dO1s))
        dO= numpy.dot(self._sigomatrixEig[1][:,self._sigomatrixEigsortIndx],
                      dO)
        Om= dO+numpy.tile(self._progenitor_Omega.T,(n,1)).T
        #Also generate angles
        da= numpy.random.normal(size=(3,n))*self._sigangle
        #And a random time
        dt= self.sample_t(n)
        #Integrate the orbits relative to the progenitor
        da+= dO*numpy.tile(dt,(3,1))
        angle= da+numpy.tile(self._progenitor_angle.T,(n,1)).T
        return (Om,angle,dt)
640274
def step(self):
        """"""
        # Assign attributes to local variables for convenience.
        u = self._u
        
        if u is None:
            raise ValueError("the levelset function is not set "
                             "(use set_levelset)")
        
        data = self.data
        
        # Determine c0 and c1.
        inside = (u > 0)
        outside = (u <= 0)
        c0 = data[outside].sum() / float(outside.sum())
        c1 = data[inside].sum() / float(inside.sum())
        
        # Image attachment.
        dres = np.array(np.gradient(u))
        abs_dres = np.abs(dres).sum(0)
        #aux = abs_dres * (c0 - c1) * (c0 + c1 - 2*data)
        aux = abs_dres * (self.lambda1*(data - c1) ** 2 -
                          self.lambda2*(data - c0) ** 2)
        
        res = np.copy(u)
        res[aux < 0] = 1
        res[aux > 0] = 0
        
        # Smoothing.
        for i in range(self.smoothing):
            res = curvop(res)
        
        self._u = res
611382
def downsample_trace(mesh, tol=1.0):
    
    idx = _find_turning_points(mesh, tol)
    if mesh.depths is not None:
        return numpy.column_stack([mesh.lons[0, idx],
                                   mesh.lats[0, idx],
                                   mesh.depths[0, idx]])
    else:
        return numpy.column_stack([mesh.lons[0, idx], mesh.lats[0, idx]])
46146
function( cylinder, head, sector ) {

    if( this.headsPerTrack < 0 || this.sectorsPerTrack < 0 )
      throw new Error( 'Unspecified device geometry' )

    return ( cylinder * this.headsPerTrack + head ) *
      this.sectorsPerTrack + ( sector - 1 )

  }
1127138
def samples(self, gp, Y_metadata=None):
        
        orig_shape = gp.shape
        gp = gp.flatten()
844865
def J(self):
        '''
        
        '''
        result = self.dr_wrt(self.x, profiler=self.profiler).copy()
        if self.profiler:
            self.profiler.harvest()
        return np.atleast_2d(result) if not sp.issparse(result) else result