Compute the GMFs for the given realization and
        yields arrays of the dtype (sid, eid, imti, gmv), one for rupture
345887
def gen_gmfs(self):
        
        self.sig_eps = []
        for computer in self.computers:
            rup = computer.rupture
            sids = computer.sids
            eids_by_rlz = rup.get_eids_by_rlz(self.rlzs_by_gsim)
            data = []
            for gs, rlzs in self.rlzs_by_gsim.items():
                num_events = sum(len(eids_by_rlz[rlzi]) for rlzi in rlzs)
                if num_events == 0:
                    continue
                # NB: the trick for performance is to keep the call to
                # compute.compute outside of the loop over the realizations
                # it is better to have few calls producing big arrays
                array, sig, eps = computer.compute(gs, num_events)
                array = array.transpose(1, 0, 2)  # from M, N, E to N, M, E
                for i, miniml in enumerate(self.min_iml):  # gmv < minimum
                    arr = array[:, i, :]
                    arr[arr < miniml] = 0
                n = 0
                for rlzi in rlzs:
                    eids = eids_by_rlz[rlzi]
                    e = len(eids)
                    if not e:
                        continue
                    for ei, eid in enumerate(eids):
                        gmf = array[:, :, n + ei]  # shape (N, M)
                        tot = gmf.sum(axis=0)  # shape (M,)
                        if not tot.sum():
                            continue
                        sigmas = sig[:, n + ei]
                        self.sig_eps.append((eid, sigmas, eps[:, n + ei]))
                        for sid, gmv in zip(sids, gmf):
                            if gmv.sum():
                                data.append((rlzi, sid, eid, gmv))
                    n += e
            yield numpy.array(data, self.gmv_dt)
875481
def mtakfac(m):
  
  u, s, v = msvd(np.asarray(m, dtype=complex))
  f = np.sqrt(u.conj().T @ v.conj())
  w = v @ f
  return w, s
1354018
def _jmomentsurfaceIntegrand(vz,vR,vT,R,z,df,sigmaR1,gamma,sigmaz1,n,m,o): #pragma: no cover because this is too slow; a warning is shown
    """"""
    return df(R,vR*sigmaR1,vT*sigmaR1*gamma,z,vz*sigmaz1,use_physical=False,
              func= (lambda x,y,z: x**n*y**m*z**o))
1151159
def downsample_with_averaging(array, factor):
    
    factor = tuple(factor)
    output_shape = tuple(int(math.ceil(s / f)) for s, f in zip(array.shape, factor))
    temp = np.zeros(output_shape, dtype=np.float32)
    counts = np.zeros(output_shape, np.int)
    for offset in np.ndindex(factor):
        part = array[tuple(np.s_[o::f] for o, f in zip(offset, factor))]
        indexing_expr = tuple(np.s_[:s] for s in part.shape)
        temp[indexing_expr] += part
        counts[indexing_expr] += 1
    return np.cast[array.dtype](temp / counts)
833120
def movMF(
    X,
    n_clusters,
    posterior_type="soft",
    force_weights=None,
    n_init=10,
    n_jobs=1,
    max_iter=300,
    verbose=False,
    init="random-class",
    random_state=None,
    tol=1e-6,
    copy_x=True,
):
    
    if n_init <= 0:
        raise ValueError(
            "Invalid number of initializations."
            " n_init=%d must be bigger than zero." % n_init
        )
    random_state = check_random_state(random_state)

    if max_iter <= 0:
        raise ValueError(
            "Number of iterations should be a positive number,"
            " got %d instead" % max_iter
        )

    best_inertia = np.infty
    X = as_float_array(X, copy=copy_x)
    tol = _tolerance(X, tol)

    if hasattr(init, "__array__"):
        init = check_array(init, dtype=X.dtype.type, copy=True)
        _validate_center_shape(X, n_clusters, init)

        if n_init != 1:
            warnings.warn(
                "Explicit initial center position passed: "
                "performing only one init in k-means instead of n_init=%d" % n_init,
                RuntimeWarning,
                stacklevel=2,
            )
            n_init = 1

    # defaults
    best_centers = None
    best_labels = None
    best_weights = None
    best_concentrations = None
    best_posterior = None
    best_inertia = None

    if n_jobs == 1:
        # For a single thread, less memory is needed if we just store one set
        # of the best results (as opposed to one set per run per thread).
        for it in range(n_init):
            # cluster on the sphere
            (centers, weights, concentrations, posterior, labels, inertia) = _movMF(
                X,
                n_clusters,
                posterior_type=posterior_type,
                force_weights=force_weights,
                max_iter=max_iter,
                verbose=verbose,
                init=init,
                random_state=random_state,
                tol=tol,
            )

            # determine if these results are the best so far
            if best_inertia is None or inertia < best_inertia:
                best_centers = centers.copy()
                best_labels = labels.copy()
                best_weights = weights.copy()
                best_concentrations = concentrations.copy()
                best_posterior = posterior.copy()
                best_inertia = inertia
    else:
        # parallelisation of movMF runs
        seeds = random_state.randint(np.iinfo(np.int32).max, size=n_init)
        results = Parallel(n_jobs=n_jobs, verbose=0)(
            delayed(_movMF)(
                X,
                n_clusters,
                posterior_type=posterior_type,
                force_weights=force_weights,
                max_iter=max_iter,
                verbose=verbose,
                init=init,
                random_state=random_state,
                tol=tol,
            )
            for seed in seeds
        )

        # Get results with the lowest inertia
        centers, weights, concentrations, posteriors, labels, inertia = zip(*results)
        best = np.argmin(inertia)
        best_labels = labels[best]
        best_inertia = inertia[best]
        best_centers = centers[best]
        best_concentrations = concentrations[best]
        best_posterior = posteriors[best]
        best_weights = weights[best]

    return (
        best_centers,
        best_labels,
        best_inertia,
        best_weights,
        best_concentrations,
        best_posterior,
    )
422598
def infer_dm(self, m, s, ds):
        
        OptimizedInverseModel.infer_dm(self, ds)
        if len(self.fmodel.dataset) == 0:
            return [[0.0]*self.dim_out]
        else:
            _, index = self.fmodel.dataset.nn_dims(m, np.hstack((s, ds)), range(len(m)), range(self.dim_x, self.dim_x + self.dim_y), k=1)
            guesses = [self.fmodel.dataset.get_dims(index[0], dims=range(len(m), self.dim_x))]
            result = []
                    
            for g in guesses:
                res = scipy.optimize.minimize(lambda dm:self._error_dm(m, dm, s), g,
                                              args        = (),
                                              method      = self.algo,
                                              options     = self.conf
                                             )
    
    
    
                d = self._error_dm(m, res.x, s)
                result.append((d, res.x))
            return [xi for fi, xi in sorted(result)][0]
398007
def gen_gmfs(gmf_set):
    
    for gmf in gmf_set:
        gmf_node = Node('gmf')
        gmf_node['IMT'] = gmf.imt
        if gmf.imt == 'SA':
            gmf_node['saPeriod'] = str(gmf.sa_period)
            gmf_node['saDamping'] = str(gmf.sa_damping)
        gmf_node['ruptureId'] = gmf.event_id
        sorted_nodes = sorted(gmf)
        gmf_node.nodes = (
            Node('node', dict(gmv=n.gmv, lon=n.location.x, lat=n.location.y))
            for n in sorted_nodes)
        yield gmf_node
587831
def gramian(self):
        
        if self.mode == 'spark':
            rdd = self.values.tordd()

            from pyspark.accumulators import AccumulatorParam

            class MatrixAccumulator(AccumulatorParam):
                def zero(self, value):
                    return zeros(shape(value))

                def addInPlace(self, val1, val2):
                    val1 += val2
                    return val1

            global mat
            init = zeros((self.shape[1], self.shape[1]))
            mat = rdd.context.accumulator(init, MatrixAccumulator())

            def outer_sum(x):
                global mat
                mat += outer(x, x)

            rdd.values().foreach(outer_sum)
            return self._constructor(mat.value, index=self.index)

        if self.mode == 'local':
            return self._constructor(dot(self.values.T, self.values), index=self.index)
1372193
def build_riskinputs(self, kind):
        
        logging.info('Building risk inputs from %d realization(s)', self.R)
        imtls = self.oqparam.imtls
        if not set(self.oqparam.risk_imtls) & set(imtls):
            rsk = ', '.join(self.oqparam.risk_imtls)
            haz = ', '.join(imtls)
            raise ValueError('The IMTs in the risk models (%s) are disjoint '
                             "from the IMTs in the hazard (%s)" % (rsk, haz))
        self.riskmodel.taxonomy = self.assetcol.tagcol.taxonomy
        with self.monitor('building riskinputs', autoflush=True):
            riskinputs = list(self._gen_riskinputs(kind))
        assert riskinputs
        logging.info('Built %d risk inputs', len(riskinputs))
        return riskinputs
370753
def compute_output(self, o, output_shape=None):
    
    if self.combine_dims:
      o = mtf.transpose(o, o.shape - self.o_dims + self.o_dims)
      o = mtf.replace_dimensions(o, self.o_dims, self.wo.shape.dims[0])
      reduced_dims = [self.wo.shape.dims[0]]
    else:
      reduced_dims = self.o_dims
    return mtf.einsum(
        [o, self.wo], output_shape=output_shape, reduced_dims=reduced_dims)
1254718
def gen_rupture_getters(dstore, slc=slice(None),
                        concurrent_tasks=1, hdf5cache=None):
    
    if dstore.parent:
        dstore = dstore.parent
    csm_info = dstore['csm_info']
    trt_by_grp = csm_info.grp_by("trt")
    samples = csm_info.get_samples_by_grp()
    rlzs_by_gsim = csm_info.get_rlzs_by_gsim_grp()
    rup_array = dstore['ruptures'][slc]
    maxweight = numpy.ceil(len(rup_array) / (concurrent_tasks or 1))
    nr, ne = 0, 0
    for grp_id, arr in general.group_array(rup_array, 'grp_id').items():
        if not rlzs_by_gsim[grp_id]:
            # this may happen if a source model has no sources, like
            # in event_based_risk/case_3
            continue
        for block in general.block_splitter(arr, maxweight):
            rgetter = RuptureGetter(
                hdf5cache or dstore.filename, numpy.array(block), grp_id,
                trt_by_grp[grp_id], samples[grp_id], rlzs_by_gsim[grp_id])
            rgetter.weight = getattr(block, 'weight', len(block))
            yield rgetter
            nr += len(block)
            ne += rgetter.num_events
    logging.info('Read %d ruptures and %d events', nr, ne)
1362846
def by_grp(self):
        
        dic = {}  # grp -> [(gsim_idx, rlzis), ...]
        for sm in self.csm_info.source_models:
            for sg in sm.src_groups:
                if not sg.eff_ruptures:
                    continue
                rlzs_by_gsim = self.get_rlzs_by_gsim(sg.trt, sm.ordinal)
                if not rlzs_by_gsim:
                    continue
                dic['grp-%02d' % sg.id] = numpy.array(
                    list(rlzs_by_gsim.values()))
        return dic
1097008
def _ttm_me_compute(self, V, edims, sdims, transp):
        
        shapeY = np.copy(self.shape)

        # Determine size of Y
        for n in np.union1d(edims, sdims):
            shapeY[n] = V[n].shape[1] if transp else V[n].shape[0]

        # Allocate Y (final result) and v (vectors for elementwise computations)
        Y = zeros(shapeY)
        shapeY = array(shapeY)
        v = [None for _ in range(len(edims))]

        for i in range(np.prod(shapeY[edims])):
            rsubs = unravel_index(shapeY[edims], i)
845455
def mdaOnes(shap, dtype=numpy.float, mask=None):
    
    res = MaskedDistArray(shap, dtype)
    res[:] = 1
    res.mask = mask
    return res
192344
def init(self):
        
        if hasattr(self, 'computers'):  # init already called
            return
        with hdf5.File(self.rupgetter.filename, 'r') as parent:
            self.weights = parent['weights'].value
        self.computers = []
        for ebr in self.rupgetter.get_ruptures(self.srcfilter):
            sitecol = self.sitecol.filtered(ebr.sids)
            try:
                computer = calc.gmf.GmfComputer(
                    ebr, sitecol, self.oqparam.imtls, self.cmaker,
                    self.oqparam.truncation_level, self.correl_model)
            except FarAwayRupture:
                # due to numeric errors, ruptures within the maximum_distance
                # when written, can be outside when read; I found a case with
                # a distance of 99.9996936 km over a maximum distance of 100 km
                continue
            self.computers.append(computer)
1126258
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(C, rup.mag) +
                self.get_distance_term(C, rup, dists.rjb, imt) +
                self.get_site_amplification(C, sites))
        # GMPE originally in cm/s/s - convert to g
        if imt.name in "PGA SA":
            mean -= np.log(100.0 * g)
        stddevs = self.get_stddevs(C, dists.rjb.shape, stddev_types, sites)
        if self.sigma_mu_epsilon:
            # Apply the epistemic uncertainty factor (sigma_mu) multiplied by
            # the number of standard deviations
            sigma_mu = self.get_sigma_mu_adjustment(C, imt, rup, dists)
            # Cap sigma_mu at 0.5 ln units
            sigma_mu[sigma_mu > 0.5] = 0.5
            # Sigma mu should not be less than the standard deviation of the
            # fault-to-fault variability
            sigma_mu[sigma_mu < C["tau_fault"]] = C["tau_fault"]
            mean += (self.sigma_mu_epsilon * sigma_mu)
        return mean, stddevs
743924
def _fix_tbit_dtype(self, array, colnums):
        
        descr = array.dtype.descr
        for i, colnum in enumerate(colnums):
            npy_type, isvar, istbit = self._get_tbl_numpy_dtype(colnum)
            if (istbit):
                coldescr = list(descr[i])
                coldescr[1] = '?'
                descr[i] = tuple(coldescr)

        return array.view(descr)