Returns the mean, standard deviation, mean gradient and standard deviation gradient at X for all the MCMC samples.
194842
def predict_withGradients(self, X):
        
        if X.ndim==1: X = X[None,:]
        ps = self.model.param_array.copy()
        means = []
        stds = []
        dmdxs = []
        dsdxs = []
        for s in self.hmc_samples:
            if self.model._fixes_ is None:
                self.model[:] = s
            else:
                self.model[self.model._fixes_] = s
            self.model._trigger_params_changed()
            m, v = self.model.predict(X)
            std = np.sqrt(np.clip(v, 1e-10, np.inf))
            dmdx, dvdx = self.model.predictive_gradients(X)
            dmdx = dmdx[:,:,0]
            dsdx = dvdx / (2*std)
            means.append(m)
            stds.append(std)
            dmdxs.append(dmdx)
            dsdxs.append(dsdx)
        self.model.param_array[:] = ps
        self.model._trigger_params_changed()
        return means, stds, dmdxs, dsdxs
981517
def sample_normal(x, y)
      gradient_x, gradient_y = sample_gradient x, y
      length = Gosu::distance 0, 0, gradient_x, gradient_y
      if length == 0
        [0, 0] # This could be NaN in edge cases.
      else
        [gradient_x / length, gradient_y / length]
      end
    end
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
623554
def update_gradients_full(self, dL_dK, X, X2=None):
        """"""
        if X2 is None:
            X2 = X

        base = np.pi * (X[:, None, :] - X2[None, :, :]) / self.period

        sin_base = np.sin( base )
        exp_dist = np.exp( -0.5* np.sum( np.square(  sin_base / self.lengthscale ), axis = -1 ) )

        dwl = self.variance * (1.0/np.square(self.lengthscale)) * sin_base*np.cos(base) * (base / self.period)

        dl = self.variance * np.square( sin_base) / np.power( self.lengthscale, 3)

        self.variance.gradient = np.sum(exp_dist * dL_dK)
        #target[0] += np.sum( exp_dist * dL_dK)

        if self.ARD1: # different periods
            self.period.gradient = (dwl * exp_dist[:,:,None] * dL_dK[:, :, None]).sum(0).sum(0)
        else:  # same period
            self.period.gradient = np.sum(dwl.sum(-1) * exp_dist * dL_dK)

        if self.ARD2: # different lengthscales
            self.lengthscale.gradient = (dl * exp_dist[:,:,None] * dL_dK[:, :, None]).sum(0).sum(0)
        else: # same lengthscales
            self.lengthscale.gradient = np.sum(dl.sum(-1) * exp_dist * dL_dK)
1074001
def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
        

        # get mean and std using the superclass
        mean, stddevs = super().get_mean_and_stddevs(
            sites, rup, dists, imt, stddev_types)

        # correction factor (see Atkinson and Boore, 2011; equation 5 at
        # page 1126 and nga08_gm_tmr.for line 508
        corr_fact = 10.0**(np.max([0, 3.888 - 0.674 * rup.mag]) -
                           (np.max([0, 2.933 - 0.510 * rup.mag]) *
                            np.log10(dists.rjb + 10.)))

        return np.log(np.exp(mean)*corr_fact), stddevs
675628
def gradients_X_X2(self, dL_dK, X, X2):
        """"""
        return self._comp_grads(dL_dK, X, X2)[3:]
685535
def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
        

        sites.vs30 = 700 * np.ones(len(sites.vs30))

        mean, stddevs = super().get_mean_and_stddevs(
            sites, rup, dists, imt, stddev_types)

        C = CauzziFaccioli2008SWISS01.COEFFS

        tau_ss = 'tau'
        log_phi_ss = np.log(10)
        mean, stddevs = _apply_adjustments(
            C, self.COEFFS_FS_ROCK[imt], tau_ss,
            mean, stddevs, sites, rup, dists.rhypo, imt, stddev_types,
            log_phi_ss)

        return mean, np.log(10 ** np.array(stddevs))
1224930
def get_deviations(self):
        

        mean_vec = self.mean()

        df = self.loc[:,:].copy()
        for col in df.columns:
            df.loc[:,col] -= mean_vec[col]
        return type(self).from_dataframe(pst=self.pst,df=df)
933194
def update_gradients_full(self, dL_dK, X, X2=None, reset=True):
        
        self.variance.gradient = np.sum(self.K(X, X2)* dL_dK)/self.variance

        #now the lengthscale gradient(s)
        dL_dr = self.dK_dr_via_X(X, X2) * dL_dK
        if self.ARD:

            tmp = dL_dr*self._inv_dist(X, X2)
            if X2 is None: X2 = X
            if use_stationary_cython:
                self.lengthscale.gradient = self._lengthscale_grads_cython(tmp, X, X2)
            else:
                self.lengthscale.gradient = self._lengthscale_grads_pure(tmp, X, X2)
        else:
            r = self._scaled_dist(X, X2)
            self.lengthscale.gradient = -np.sum(dL_dr*r)/self.lengthscale
670889
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]
        # Implements mean model (equation 12)
        mean = (self._compute_magnitude(rup, C) +
                self._compute_distance(dists, C) +
                self._get_site_amplification(sites, rup, C) +
                self._get_mechanism(rup, C))

        stddevs = self._get_stddevs(rup, np.exp(mean), stddev_types, sites)

        return mean, stddevs
110424
def _get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types, 
        mag_conversion_sigma=0.0):
        
        # extract dictionaries of coefficients specific to required
        # intensity measure type
        C = self.COEFFS[imt]

        # clip distance at 4 km, minimum distance for which the equation is
        # valid (see section 2.2.4, page 201). This also avoids singularity
        # in the equation
        rhypo = dists.rhypo
        rhypo[rhypo < 4.] = 4.

        mean = C['a'] * rup.mag + C['b'] * rhypo - np.log10(rhypo)

        mean[sites.vs30 >= 800] += C['c1']
        mean[sites.vs30 < 800] += C['c2']

        # convert from log10 to ln, and from cm/s2 to g
        mean = mean * np.log(10) - 2 * np.log(10) - np.log(g)

        stddevs = self._get_stddevs(C, stddev_types, rhypo.shape[0], 
            mag_conversion_sigma=mag_conversion_sigma)

        return mean, stddevs
381182
def gradients_X(self, dL_dK, X, X2, target):
        """"""
        # First account for gradients arising from presence of X in exponent.
        self._K_computations(X, X2)
        if X2 is None:
            _K_dist = 2*(X[:, None, :] - X[None, :, :])
        else:
            _K_dist = X[:, None, :] - X2[None, :, :] # don't cache this in _K_co
        gradients_X = (-2.*self.variance)*np.transpose((self._K_dvar/self._w2)[:, :, None]*_K_dist, (1, 0, 2))
        target += np.sum(gradients_X*dL_dK.T[:, :, None], 0)
        # Now account for gradients arising from presence of X in lengthscale.
        self._dK_computations(dL_dK)
        if X2 is None:
            target += 2.*self.mapping.df_dX(self._dL_dl[:, None], X)
        else:
            target += self.mapping.df_dX(self._dL_dl[:, None], X)
552411
def examples(self):
        
        # TODO: adapt for tensor-valued functions

        # Get the points and calculate some statistics on them
        mins = self.domain.min()
        maxs = self.domain.max()
        means = (maxs + mins) / 2.0
        stds = (maxs - mins) / 4.0
        ndim = getattr(self.domain, 'ndim', 0)

        # Zero and One
        yield ('Zero', self.zero())
        yield ('One', self.one())

        # Indicator function in first dimension
        def step_fun(x):
            return (x[0] > means[0])

        yield ('Step', self.element(step_fun))

        # Indicator function on hypercube
        def cube_fun(x):
            result = True
            for points, mean, std in zip(x, means, stds):
                result = np.logical_and(result, points < mean + std)
                result = np.logical_and(result, points > mean - std)
            return result

        yield ('Cube', self.element(cube_fun))

        # Indicator function on a ball
        if ndim > 1:  # Only if ndim > 1, don't duplicate cube
            def ball_fun(x):
                r = sum((xi - mean) ** 2 / std ** 2
                        for xi, mean, std in zip(x, means, stds))
                return r < 1.0

            yield ('Ball', self.element(ball_fun))

        # Gaussian function
        def gaussian_fun(x):
            r2 = sum((xi - mean) ** 2 / (2 * std ** 2)
                     for xi, mean, std in zip(x, means, stds))
            return np.exp(-r2)

        yield ('Gaussian', self.element(gaussian_fun))

        # Gradient in each dimensions
        for axis in range(ndim):
            def gradient_fun(x):
                return (x[axis] - mins[axis]) / (maxs[axis] - mins[axis])

            yield ('Grad {}'.format(axis), self.element(gradient_fun))

        # Gradient in all dimensions
        if ndim > 1:  # Only if ndim > 1, don't duplicate grad 0
            def all_gradient_fun(x):
                return sum((xi - xmin) / (xmax - xmin)
                           for xi, xmin, xmax in zip(x, mins, maxs))

            yield ('Grad all', self.element(all_gradient_fun))
485759
def _get_stddevs(self, C, stddev_types, num_sites, mag, c1_rrup,
                     log_phi_ss, mean_phi_ss):
        
        phi_ss = _compute_phi_ss(C, mag, c1_rrup, log_phi_ss, mean_phi_ss)

        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(
                               C['tau'] * C['tau'] +
                               phi_ss * phi_ss) +
                               np.zeros(num_sites))

            elif stddev_type == const.StdDev.INTRA_EVENT:
                stddevs.append(phi_ss + np.zeros(num_sites))

            elif stddev_type == const.StdDev.INTER_EVENT:
                stddevs.append(C['tau'] + np.zeros(num_sites))
        return stddevs
138161
def gradients_XX_diag(self, dL_dK_diag, X):
        
        dL_dK_diag = dL_dK_diag.copy().reshape(-1, 1, 1)
        assert (dL_dK_diag.size == X.shape[0]) or (dL_dK_diag.size == 1), "dL_dK_diag has to be given as row [N] or column vector [Nx1]"

        l4 =  np.ones(X.shape[1])*self.lengthscale**2
        return dL_dK_diag * (np.eye(X.shape[1]) * -self.dK2_drdr_diag()/(l4))[None, :,:]
505948
def get_optimization_gradients(self):
        """"""
        grad_array = self.fields.get("Opt point       1 Gradient at each geome")
        if grad_array is None:
            return []
        else:
            return np.reshape(grad_array, (-1, len(self.molecule.numbers), 3))
547341
def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
        
        mean, stddevs = super().get_mean_and_stddevs(
            sites, rup, dists, imt, stddev_types)

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

        mean[idx_rock] += 0.275
        mean[idx_soil] += 0.31

        return mean, stddevs