returns bootstrap means  for Directional data
     Parameters
     _________________
     DIs : nested list of Dec,Inc pairs
     nb : number of bootstrap pseudosamples

     Returns
    -------
     BDIs:   nested list of bootstrapped mean Dec,Inc pairs
770792
def di_boot(DIs, nb=5000):
    
#
# now do bootstrap to collect BDIs  bootstrap means
#
    BDIs = []  # number of bootstraps, list of bootstrap directions
#

    for k in range(nb):  # repeat nb times
        #        if k%50==0:print k,' out of ',nb
        pDIs = pseudo(DIs)  # get a pseudosample
        bfpars = fisher_mean(pDIs)  # get bootstrap mean bootstrap sample
        BDIs.append([bfpars['dec'], bfpars['inc']])
    return BDIs
399810
def get_mean_and_stddevs(self, sctx, rctx, dctx, imt, stddev_types):
        
        return self.kwargs[str(imt)].get_mean_and_stddevs(
            sctx, rctx, dctx, imt, stddev_types)
1025464
def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
        
        assert all(stddev_type in self.DEFINED_FOR_STANDARD_DEVIATION_TYPES
                   for stddev_type in stddev_types)

        mean = np.zeros_like(sites.vs30)
        stddevs = [np.zeros_like(sites.vs30) for _ in stddev_types]

        idx_rock = sites.vs30 >= self.ROCK_VS30
        idx_soil = sites.vs30 < self.ROCK_VS30

        if idx_rock.any():
            C = self.COEFFS_ROCK[imt]
            self._compute_mean(C, rup.mag, dists.rhypo, rup.hypo_depth, mean,
                               idx_rock)
            self._compute_std(C, stddevs, idx_rock)

        if idx_soil.any():
            C = self.COEFFS_SOIL[imt]
            self._compute_mean(C, rup.mag, dists.rhypo, rup.hypo_depth, mean,
                               idx_soil)
            self._compute_std(C, stddevs, idx_soil)

        return mean, stddevs
372682
def _get_mean(self, sites, rup, dists, imt, stddev_types):
        
        # Zhao et al. 2006 - Vs30 + Rrup
        mean_zh06, stds1 = super().get_mean_and_stddevs(sites, rup, dists, imt,
                                                        stddev_types)
        #
        # Atkinson and Macias (2009) - Rrup
        gmpe = AtkinsonMacias2009()
        mean_am09, stds2 = gmpe.get_mean_and_stddevs(sites, rup, dists, imt,
                                                     stddev_types)
        #
        # Abrahamson et al. (2015) - Rrup + vs30 + backarc
        gmpe = AbrahamsonEtAl2015SInter()
        mean_ab15, stds3 = gmpe.get_mean_and_stddevs(sites, rup, dists, imt,
                                                     stddev_types)
        #
        # Ghofrani and Atkinson (2014) - Rrup + vs30
        gmpe = GhofraniAtkinson2014()
        mean_ga14, stds4 = gmpe.get_mean_and_stddevs(sites, rup, dists, imt,
                                                     stddev_types)
        # Computing adjusted mean and stds
        cff = self.SITE_COEFFS[imt]
        mean_adj = (np.log(np.exp(mean_zh06)*cff['mf'])*0.1 +
                    mean_am09*0.5 + mean_ab15*0.2 +
                    np.log(np.exp(mean_ga14)*cff['mf'])*0.2)
        return mean_adj
549272
def get_mean_and_stddevs(self, sctx, rctx, dctx, imt, stddev_types):
        
        # Get original mean and standard deviations
        mean, stddevs = super().get_mean_and_stddevs(
            sctx, rctx, dctx, imt, stddev_types)
        cff = SInterCan15Mid.SITE_COEFFS[imt]
        mean += np.log(cff['mf'])
        return mean, stddevs
817908
def bootstrap_indexes(data, n_samples=10000):
    
    for _ in xrange(n_samples):
        yield randint(data.shape[0], size=(data.shape[0],))
981445
def insert_size_stats(dists):
    
    med = numpy.median(dists)
    filter_dists = list(filter(lambda x: x < med + 10 * med, dists))
    median = numpy.median(filter_dists)
    return {"mean": float(numpy.mean(filter_dists)), "std": float(numpy.std(filter_dists)),
            "median": float(median),
            "mad": float(numpy.median([abs(x - median) for x in filter_dists]))}
252611
def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
        
        assert all(stddev_type in self.DEFINED_FOR_STANDARD_DEVIATION_TYPES
                   for stddev_type in stddev_types)

        # GMPE differentiates strike-slip, reverse and normal ruptures,
        # but combines normal and strike-slip into one category. See page 180.
        is_reverse = (45 <= rup.rake <= 135)

        stddevs = [numpy.zeros_like(sites.vs30) for _ in stddev_types]
        means = numpy.zeros_like(sites.vs30)

        [rocks_i] = (sites.vs30 > self.ROCK_VS30).nonzero()
        if len(rocks_i):
            rrup = dists.rrup.take(rocks_i)
            mean_rock = self._get_mean_rock(rup.mag, rup.rake, rrup,
                                            is_reverse, imt)
            means.put(rocks_i, mean_rock)
            for stddev_arr in stddevs:
                stddev_rock = self._get_stddev_rock(rup.mag, imt)
                stddev_arr.put(rocks_i, stddev_rock)

        [soils_i] = (sites.vs30 <= self.ROCK_VS30).nonzero()
        if len(soils_i):
            rrup = dists.rrup.take(soils_i)
            mean_soil = self._get_mean_deep_soil(rup.mag, rup.rake, rrup,
                                                 is_reverse, imt)
            means.put(soils_i, mean_soil)
            for stddev_arr in stddevs:
                stddev_soil = self._get_stddev_deep_soil(rup.mag, imt)
                stddev_arr.put(soils_i, stddev_soil)

        return means, stddevs
1246578
def _get_stddevs(self, C, imt, rup, sites, stddev_types, sa1180, dists):
        
        std_intra = self._get_intra_event_std(C, rup.mag, sa1180, sites.vs30,
                                              sites.vs30measured, dists.rrup)
        std_inter = self._get_inter_event_std(C, rup.mag, sa1180, sites.vs30)
        stddevs = []
        for stddev_type in stddev_types:
            assert stddev_type in self.DEFINED_FOR_STANDARD_DEVIATION_TYPES
            if stddev_type == const.StdDev.TOTAL:
                stddevs.append(np.sqrt(std_intra ** 2 +
                                       std_inter ** 2))
            elif stddev_type == const.StdDev.INTRA_EVENT:
                stddevs.append(std_intra)
            elif stddev_type == const.StdDev.INTER_EVENT:
                stddevs.append(std_inter)
        return stddevs
149274
def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
        
        # Check that the GSIM supports the standard deviations requested
        assert all(stddev_type in self.DEFINED_FOR_STANDARD_DEVIATION_TYPES
                   for stddev_type in stddev_types)
        coe = self.COEFFS[imt]
        mean = (self._get_magnitude_scaling(coe, rup.mag) +
                self._get_distance_scaling(coe, dists.rhypo) +
                self._get_azimuth_correction(coe, dists.azimuth))
        # Convert to g
        if imt.name in "SA PGA":
            mean = np.log(np.exp(mean) / (100.0 * g))
        # Compute std
        stddevs = self._compute_std(coe, stddev_types, dists.azimuth.shape)
        return mean, stddevs
1183718
def get_stddevs(self, mean, stddev_types):
        
        mean = np.exp(mean)
        sigma = 0.39 + np.zeros(mean.shape)
        sigma[mean < 0.068] = 0.55
        idx = np.logical_and(mean >= 0.068, mean <= 0.21)
        sigma[idx] = 0.173- 0.140 * np.log(mean[idx])
        stddevs = []
        for stddev in stddev_types:
            if stddev == const.StdDev.TOTAL:
                stddevs.append(sigma)
        return stddevs
1093004
def _friends_bootstrap_radius(args):
    """"""

    # Unzipping.
    points, ftype = 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]

    # Construct KDTree to enable quick nearest-neighbor lookup for
    # our resampled objects.
    kdtree = spatial.KDTree(points_in)

    if ftype == 'balls':
        # Compute distances from our "missing" points its closest neighbor
        # among the resampled points using the Euclidean norm
        # (i.e. "radius" of n-sphere).
        dists, ids = kdtree.query(points_out, k=1, eps=0, p=2)
    elif ftype == 'cubes':
        # Compute distances from our "missing" points its closest neighbor
        # among the resampled points using the Euclidean norm
        # (i.e. "half-side-length" of n-cube).
        dists, ids = kdtree.query(points_out, k=1, eps=0, p=np.inf)

    # Conservative upper-bound on radius.
    dist = max(dists)

    return dist
492095
def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
        
        # extract dictionaries of coefficients specific to required
        # intensity measure type and for PGA
        C = self.COEFFS[imt]
        C_PGA = self.COEFFS[PGA()]

        # Get mean and standard deviation of PGA on rock (Vs30 1100 m/s^2)
        pga1100 = np.exp(self.get_mean_values(C_PGA, sites, rup, dists, None))
        # Get mean and standard deviations for IMT
        mean = self.get_mean_values(C, sites, rup, dists, pga1100)
        if imt.name == "SA" and imt.period <= 0.25:
            # According to Campbell & Bozorgnia (2013) [NGA West 2 Report]
            # If Sa (T) < PGA for T < 0.25 then set mean Sa(T) to mean PGA
            # Get PGA on soil
            pga = self.get_mean_values(C_PGA, sites, rup, dists, pga1100)
            idx = mean <= pga
            mean[idx] = pga[idx]
        # Get standard deviations
        stddevs = self._get_stddevs(C,
                                    C_PGA,
                                    rup,
                                    sites,
                                    pga1100,
                                    stddev_types)
        return mean, stddevs
960558
def _get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
        
        # distances
        distsl = copy.copy(dists)
        distsl.rjb, distsl.rrup = \
            utils.get_equivalent_distances_east(rup.mag, dists.repi)
        #
        # Pezeshk et al. 2011 - Rrup
        mean1, stds1 = super().get_mean_and_stddevs(sites, rup, distsl, imt,
                                                    stddev_types)
        mean1 = self.apply_correction_to_BC(mean1, imt, distsl)
        #
        # Atkinson 2008 - Rjb
        gmpe = Atkinson2008prime()
        mean2, stds2 = gmpe.get_mean_and_stddevs(sites, rup, distsl, imt,
                                                 stddev_types)
        #
        # Silva et al. 2002 - Rjb
        gmpe = SilvaEtAl2002SingleCornerSaturation()
        mean4, stds4 = gmpe.get_mean_and_stddevs(sites, rup, distsl, imt,
                                                 stddev_types)
        mean4 = self.apply_correction_to_BC(mean4, imt, distsl)
        #
        # Silva et al. 2002 - Rjb
        gmpe = SilvaEtAl2002DoubleCornerSaturation()
        mean5, stds5 = gmpe.get_mean_and_stddevs(sites, rup, distsl, imt,
                                                 stddev_types)
        mean5 = self.apply_correction_to_BC(mean5, imt, distsl)
        #
        # distances
        distsl.rjb, distsl.rrup = \
            utils.get_equivalent_distances_east(rup.mag, dists.repi, ab06=True)
        #
        # Atkinson and Boore 2006 - Rrup
        gmpe = AtkinsonBoore2006Modified2011()
        mean3, stds3 = gmpe.get_mean_and_stddevs(sites, rup, distsl, imt,
                                                 stddev_types)
        # Computing adjusted mean and stds
        mean_adj = mean1*0.2 + mean2*0.2 + mean3*0.2 + mean4*0.2 + mean5*0.2

        # Note that in this case we do not apply a triangular smoothing on
        # distance as explained at page 996 of Atkinson and Adams (2013)
        # for the calculation of the standard deviation
        stds_adj = np.log(np.exp(stds1)*0.2 + np.exp(stds2)*0.2 +
                          np.exp(stds3)*0.2 + np.exp(stds4)*0.2 +
                          np.exp(stds5)*0.2)
        #
        return mean_adj, stds_adj
86290
def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
        # pylint: disable=too-many-arguments
        

        # extract dictionary of coefficients specific to required
        # intensity measure type
        coeffs = self.COEFFS[imt]
        coeffs.update(self.CONSTS)

        # equation (1) is in terms of common logarithm
        log_mean = (self._compute_magnitude(rup, coeffs) +
                    self._compute_distance(dists, coeffs) +
                    self._get_site_amplification(sites, coeffs) +
                    self._get_mechanism(rup, coeffs))
        # so convert to g and thence to the natural logarithm
        mean = log_mean*np.log(10.0) - np.log(g)

        # convert standard deviations from common to natural logarithm
        log_stddevs = self._get_stddevs(coeffs, stddev_types, len(sites.vs30))
        stddevs = log_stddevs*np.log(10.0)

        return mean, stddevs
688655
def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
        
        # extracting dictionary of coefficients specific to required
        # intensity measure type.
        C = self.COEFFS[imt]

        mean = (self._get_magnitude_scaling_term(C, rup.mag) +
                self._get_distance_scaling_term(C, rup.mag, dists.rrup) +
                self._get_style_of_faulting_term(C, rup.rake) +
                self._get_site_scaling_term(C, sites.vs30))

        stddevs = self._get_stddevs(imt,
                                    rup.mag,
                                    len(dists.rrup),
                                    stddev_types)
        return mean, stddevs
960310
def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
        
        # extracting dictionary of coefficients specific to required
        # intensity measure type.

        C = self.COEFFS[imt]
        imean = (self._compute_magnitude(rup, C) +
                 self._compute_distance(rup, dists, C) +
                 self._get_site_amplification(sites, C) +
                 self._compute_forearc_backarc_term(C, sites, dists, rup))

        istddevs = self._get_stddevs(C,
                                     stddev_types,
                                     num_sites=len(sites.vs30))

        # Convert units to g,
        # but only for PGA and SA (not PGV):
        if imt.name in "SA PGA":
            mean = np.log((10.0 ** (imean - 2.0)) / g)
        else:
            # PGV:
            mean = np.log(10.0 ** imean)
        # Return stddevs in terms of natural log scaling
        stddevs = np.log(10.0 ** np.array(istddevs))
        # mean_LogNaturale = np.log((10 ** mean) * 1e-2 / g)
        return mean, stddevs