# The MIT License (MIT)
#
# Copyright (c) 2016-2018 Albert Kottke
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
import copy
from abc import abstractmethod
from collections.abc import Generator
import numpy as np
import numpy.typing as npt
from scipy import stats
from scipy.sparse import diags
from . import site
# Used to define the random state. A specific state can be set with:
# random_state.set_seed(42)
random_state = np.random.RandomState()
class TruncatedNorm:
"""Truncated normal random number generator.
Parameters
----------
limit : float
Standard normal limits to impose
"""
def __init__(self, limit):
self.limit = limit
@property
def limit(self):
return self._limit
@limit.setter
def limit(self, value):
self._limit = value
# Need to scale the standard deviation to achieve sample standard
# deviation based on the truncation. Given truncation of 2 standard
# deviations, the input standard deviation must be increased to
# 1.136847 to maintain a unit standard deviation for the random
# samples.
self._scale = 1 / np.sqrt(stats.truncnorm.stats(-value, value, moments="v"))
@property
def scale(self):
return self._scale
def __call__(self, size=1):
"""Random number generator that follows a truncated normal distribution.
This is the defalut random number generator used by the program. It
generates normally distributed values ranging from -2 to +2 with unit
standard deviation.
The state of the random number generator is controlled by the
``np.random.RandomState`` instance.
Parameters
----------
size : int
Number of random values to compute
Returns
-------
rvs : ndarray or scalar
Random variates of given `size`.
"""
return stats.truncnorm.rvs(
-self.limit, self.limit, scale=self._scale, size=size
)
def correlated(self, correl):
# Acceptance proportion
accept = np.diff(stats.norm.cdf([-self.limit, self.limit]))[0]
# The expected number of tries required
expected = np.ceil(1 / accept).astype(int)
while True:
# Compute the multivariate normal with a unit variance and
# specified standard deviation. Use twice the expected since
# this calculation is fast and we don't want to loop.
randvar = np.random.multivariate_normal(
[0, 0], [[1, correl], [correl, 1]], size=(2 * expected)
)
valid = np.all(np.abs(randvar) < self.limit, axis=1)
if np.any(valid):
# Return the first valid value
return randvar[valid][0]
def correlated_at_percentile(self, correl, percentile):
"""Return a deterministic correlated pair for a given percentile.
The primary variate ``z1`` is the inverse CDF of the *scaled*
truncated normal at ``percentile``. The secondary variate ``z2``
is drawn from the conditional distribution
:math:`z_2 | z_1 = \\rho z_1 + \\sqrt{1 - \\rho^2}\\, z_1`,
also clipped to ``[-limit, +limit]``.
Both variates are expressed in the same units as the samples
returned by :meth:`correlated` (unit standard deviation, truncated
at ``±limit``).
Parameters
----------
correl : float
Correlation coefficient between the two variates.
percentile : float
Quantile in ``(0, 1)`` to evaluate.
Returns
-------
pair : ndarray, shape (2,)
Deterministic correlated variate pair.
"""
# Map the percentile through the scaled truncated-normal CDF so that
# the resulting z1 sits at exactly that percentile of the marginal
# distribution used by the random sampler.
# stats.truncnorm with scale == self._scale and limits
# [-limit/scale, +limit/scale] (un-scaled) matches the marginals.
a, b = -self._limit / self._scale, self._limit / self._scale
z1 = stats.truncnorm.ppf(percentile, a, b, scale=self._scale)
# Conditional mean of z2 given z1 under bivariate normal
z2 = correl * z1 + np.sqrt(max(1.0 - correl**2, 0.0)) * z1
# Clip both to the truncation window
z1 = np.clip(z1, -self._limit, self._limit)
z2 = np.clip(z2, -self._limit, self._limit)
return np.array([z1, z2])
# Random number generator used for all random number. Limited to +/- 2,
# and the standard deviation is scaled to maintain the standard deviation
# FIXME
randnorm = TruncatedNorm(2)
class ToroThicknessVariation:
"""Toro (1995) [T95]_ thickness variation model.
The recommended values are provided as defaults to this model.
.. rubric:: References
.. [T95] Toro, G. R. (1995). Probabilistic models of site velocity
profiles for generic and site-specific ground-motion amplification
studies. Brookhaven National Laboratory Technical Report: 779574.
Parameters
----------
c_1: float, optional
:math:`c_1` model parameter.
c_2: float, optional
:math:`c_2` model parameter.
c_3: float, optional
:math:`c_3` model parameter.
"""
def __init__(self, c_1=10.86, c_2=-0.89, c_3=1.98):
self._c_1 = c_1
self._c_2 = c_2
self._c_3 = c_3
@property
def c_3(self):
return self._c_3
@property
def c_2(self):
return self._c_2
@property
def c_1(self):
return self._c_1
def iter_thickness(self, depth_total):
r"""Iterate over the varied thicknesses.
The layering is generated using a non-homogenous Poisson process. The
following routine is used to generate the layering. The rate
function, :math:`\lambda(t)`, is integrated from 0 to t to generate
cumulative rate function, :math:`\Lambda(t)`. This function is then
inverted producing :math:`\Lambda^{-1}(t)`. Random variables
are produced using the a exponential random variation with
:math:`\mu = 1` and converted to the nonhomogenous variables using
the inverted function.
Parameters
----------
depth_total: float
Total depth generated. Last thickness is truncated to achieve
this depth.
Yields
------
float
Varied thickness.
"""
total = 0
depth_prev = 0
while depth_prev < depth_total:
# Add a random exponential increment
total += np.random.exponential(1.0)
# Convert between x and depth using the inverse of \Lambda(t)
depth = (
np.power(
(self.c_2 * total) / self.c_3
+ total / self.c_3
+ np.power(self.c_1, self.c_2 + 1),
1 / (self.c_2 + 1),
)
- self.c_1
)
thickness = depth - depth_prev
if depth > depth_total:
thickness = depth_total - depth_prev
depth = depth_prev + thickness
depth_mid = (depth_prev + depth) / 2
yield thickness, depth_mid
depth_prev = depth
def __call__(self, profile):
"""Calculated a varied thickness profile.
Parameters
----------
profile : site.Profile
Profile to be varied. Not modified in place.
Returns
-------
site.Profile
Varied site profile.
"""
layers = []
for thick, depth_mid in self.iter_thickness(profile[-2].depth_base):
# Locate the proper layer and add it to the model
for layer in profile:
if layer.depth < depth_mid <= layer.depth_base:
layers.append(
site.Layer(
layer.soil_type,
thick,
layer.initial_shear_vel,
layer.damping_min,
)
)
break
else:
raise LookupError
# Add the half-space
hsl = profile[-1]
layers.append(site.Layer(hsl.soil_type, 0, hsl.initial_shear_vel))
varied = site.Profile(layers, profile.wt_depth)
return varied
class HalfSpaceDepthVariation:
def __init__(self, dist: stats.rv_continuous):
self._dist = dist
def __call__(self, profile: site.Profile) -> site.Profile:
# Update the distribution with the central value of the profile
varied_depth = self._dist.rvs()
# Find the layer
index, depth_within = profile.lookup_depth(varied_depth)
half_space = profile[-1]
if index < (len(profile) - 1):
# Variation is within the layers
layers = [layer.copy() for layer in profile[: (index + 1)]]
# Reduce the thickness of the layer above the half-space
layers[-1]._thickness = depth_within
else:
# Variation extends past the depth of the model
orig_thick = profile[-2].thickness
total_thick = orig_thick + depth_within
count = np.ceil(total_thick // orig_thick).astype(int)
thick = total_thick / count
# Don't copy half-space
layers = [layer.copy() for layer in profile[:-1]]
# Don't call the setter function as it needs a profile defined
layers[-1]._thickness = thick
parent = layers[-1]
for _ in range(count - 1):
layers.append(parent.copy())
layers.append(half_space)
return site.Profile(layers, profile.wt_depth)
class LayerThicknessVariation:
def __init__(
self,
models: list[stats.rv_continuous] | dict[int, stats.rv_continuous],
discretize_kwds: dict[str, float] | None = None,
) -> None:
self._models = models
self._discretize_kwds = discretize_kwds
class VelocityVariation:
"""Abstract model for varying the velocity."""
def __init__(self, vary_bedrock=False):
self._vary_bedrock = vary_bedrock
def __call__(self, profile: site.Profile) -> site.Profile:
"""Calculate a varied shear-wave velocity profile.
Parameters
----------
profile : site.Profile
Profile to be varied. Not modified in place.
Returns
-------
site.Profile
Varied site profile.
"""
mean = np.log(profile.initial_shear_vel)
covar = self._calc_covar_matrix(profile)
ln_vel_rand = np.random.multivariate_normal(mean, covar, check_valid="ignore")
# Limits based on the number of standard deviations
offset = randnorm.limit * np.sqrt(np.diag(covar))
ln_vel = np.clip(ln_vel_rand, mean - offset, mean + offset)
vel = np.exp(ln_vel)
varied = profile.copy()
# Update the velocities
end = None if self.vary_bedrock else -1
for i, v in enumerate(vel[:end]):
varied[i].initial_shear_vel = v
return varied
def _calc_covar_matrix(self, profile):
"""Calculate the covariance matrix.
Parameters
----------
profile : site.Profile
Input site profile
Yields
------
covar : `class`:numpy.array
Covariance matrix
"""
corr = self._calc_corr(profile)
std = self._calc_ln_std(profile)
# Modify the standard deviation by the truncated norm scale
std *= randnorm.scale
var = std**2
covar = corr * std[:-1] * std[1:]
# Main diagonal is the variance
mat = diags([covar, var, covar], [-1, 0, 1]).toarray()
return mat
@abstractmethod
def _calc_corr(self, profile):
"""Compute the adjacent-layer correlations.
Parameters
----------
profile : site.Profile
Input site profile
Yields
------
np.array
Correlation matrix
"""
raise NotImplementedError
@abstractmethod
def _calc_ln_std(self, profile):
"""Compute the standard deviation for each layer.
Parameters
----------
profile : site.Profile
Input site profile
Yields
------
np.array
Standard deviation of the shear-wave velocity
"""
raise NotImplementedError
@property
def vary_bedrock(self):
return self._vary_bedrock
[docs]
class ToroVelocityVariation(VelocityVariation):
r"""Toro (1995) [T95] velocity variation model.
Default values can be selected with :meth:`.generic_model`.
Parameters
----------
ln_std: float, optional
:math:`\sigma_{ln}` model parameter.
rho_0: float, optional
:math:`ρ_0` model parameter.
delta: float, optional
:math:`\Delta` model parameter.
rho_200: float, optional
:math:`ρ_200` model parameter.
h_0: float, optional
:math:`h_0` model parameter.
b: float, optional
:math:`b` model parameter.
vary_bedrock: bool, optional
If the velocity of the bedrock (half-space) should be varied.
"""
PARAMS = {
"Geomatrix AB": {
"ln_std": 0.46,
"rho_0": 0.96,
"delta": 13.1,
"rho_200": 0.96,
"h_0": 0.0,
"b": 0.095,
},
"Geomatrix CD": {
"ln_std": 0.38,
"rho_0": 0.99,
"delta": 8.0,
"rho_200": 1.00,
"h_0": 0.0,
"b": 0.160,
},
"USGS AB": {
"ln_std": 0.35,
"rho_0": 0.95,
"delta": 4.2,
"rho_200": 1.00,
"h_0": 0.0,
"b": 0.138,
},
"USGS CD": {
"ln_std": 0.36,
"rho_0": 0.99,
"delta": 3.9,
"rho_200": 1.00,
"h_0": 0.0,
"b": 0.293,
},
"USGS A": {
"ln_std": 0.36,
"rho_0": 0.95,
"delta": 3.4,
"rho_200": 0.42,
"h_0": 0.0,
"b": 0.063,
},
"USGS B": {
"ln_std": 0.27,
"rho_0": 0.97,
"delta": 3.8,
"rho_200": 1.00,
"h_0": 0.0,
"b": 0.293,
},
"USGS C": {
"ln_std": 0.31,
"rho_0": 0.99,
"delta": 3.9,
"rho_200": 0.98,
"h_0": 0.0,
"b": 0.344,
},
"USGS D": {
"ln_std": 0.37,
"rho_0": 0.00,
"delta": 5.0,
"rho_200": 0.50,
"h_0": 0.0,
"b": 0.744,
},
}
[docs]
def __init__(
self,
ln_std: float,
rho_0: float,
delta: float,
rho_200: float,
h_0: float,
b: float,
vary_bedrock: bool = False,
):
"""Initialize the model."""
super().__init__(vary_bedrock=vary_bedrock)
self._ln_std = ln_std
self._rho_0 = rho_0
self._delta = delta
self._rho_200 = rho_200
self._h_0 = h_0
self._b = b
def _calc_corr(self, profile: site.Profile) -> np.ndarray:
"""Compute the adjacent-layer correlations.
Parameters
----------
profile : :class:`site.Profile`
Input site profile
Yields
------
corr : :class:`numpy.array`
Adjacent-layer correlations
"""
# Toro defines the depth as the average midpoint depths of layers i and i-1.
depths_mid = np.array(profile.depth_mid)
depth = np.mean(np.c_[depths_mid[:-1], depths_mid[1:]], axis=1)
# t variable from Toro; defined as the difference of the midpoint depths
# Here the thickness is limited to 100 m to prevent underflow on the
# exponent. For this thickness, the correlation will be minor
thick = np.minimum(np.diff(depth), 100)
# Remove the depth associated with the final layer. We will set that it is
# perfectly correlated later.
depth = depth[:-1]
# Depth dependent correlation
corr_depth = self.rho_200 * np.power(
(depth + self.h_0) / (200 + self.h_0), self.b
)
corr_depth[depth > 200] = self.rho_200
# Thickness dependent correlation
corr_thick = self.rho_0 * np.exp(-thick / self.delta)
# Final correlation
# Correlation coefficient
corr = (1 - corr_depth) * corr_thick + corr_depth
# Bedrock is perfectly correlated with layer above it
corr = np.r_[corr, 1]
return corr
def _calc_ln_std(self, profile):
ln_std = self.ln_std * np.ones(len(profile))
return ln_std
@property
def ln_std(self):
return self._ln_std
@property
def rho_0(self):
return self._rho_0
@property
def delta(self):
return self._delta
@property
def rho_200(self):
return self._rho_200
@property
def h_0(self):
return self._h_0
@property
def b(self):
return self._b
@classmethod
def site_classes(cls):
return cls.PARAMS.keys()
[docs]
@classmethod
def generic_model(cls, site_class, **kwds):
"""Use generic model parameters based on site class.
Parameters
----------
site_class: str
Site classification. Possible options are:
* Geomatrix AB
* Geomatrix CD
* USGS AB
* USGS CD
* USGS A
* USGS B
* USGS C
* USGS D
See the report for definitions of the Geomatrix site
classication. USGS site classification is based on :math:`V_{s30}`:
=========== =====================
Site Class :math:`V_{s30}` (m/s)
=========== =====================
A >750 m/s
B 360 to 750 m/s
C 180 to 360 m/s
D <180 m/s
=========== =====================
Returns
-------
ToroVelocityVariation
Initialized :class:`ToroVelocityVariation` with generic parameters.
"""
p = dict(cls.PARAMS[site_class])
p.update(kwds)
return cls(**p)
@staticmethod
def usgs_site_class(vs30):
return "USGS " + np.array(list("DCBA"))[np.searchsorted([180, 360, 750], vs30)]
class DepthDependToroVelVariation(ToroVelocityVariation):
r"""Toro (1995) [T95] velocity variation model modified for a depth dependent
standard deviation that can be overridden by the soil_type name.
Default values can be selected with :meth:`.generic_model`.
Parameters
----------
depth: array_like, optional
Depths defining the standard deviation model. Default is [0, 15]
following the SPID model.
ln_std: array_like, optional
:math:`\sigma_{ln}` model parameter. Default is [0.25, 0.15]
following the SPID model.
rho_0: float, optional
:math:`ρ_0` model parameter.
delta: float, optional
:math:`\Delta` model parameter.
rho_200: float, optional
:math:`ρ_200` model parameter.
h_0: float, optional
:math:`h_0` model parameter.
b: float, optional
:math:`b` model parameter.
vary_bedrock: bool, optional
If the velocity of the bedrock (half-space) should be varied.
ln_std_map: dict[str, float], optional
Mapping between the soil_type and the defined ln_std. Default is *None*.
"""
def __init__(
self,
depth: npt.ArrayLike,
ln_std: npt.ArrayLike,
rho_0: float,
delta: float,
rho_200: float,
h_0: float,
b: float,
vary_bedrock: bool = False,
ln_std_map: dict[str, float] | None = None,
):
"""Initialize the model."""
super().__init__(
ln_std, rho_0, delta, rho_200, h_0, b, vary_bedrock=vary_bedrock
)
self.depth = depth
self.ln_std_map = ln_std_map or dict()
def _calc_ln_std(self, profile):
# Depth based values
ln_std = np.interp(
profile.depth_mid,
self.depth,
self.ln_std,
left=self.ln_std[0],
right=self.ln_std[-1],
)
# Update based on soil_type name
for name, value in self.ln_std_map.items():
for i, layer in enumerate(profile):
if name in layer.soil_type.name:
ln_std[i] = value
return ln_std
@classmethod
def generic_model(cls, site_class, /, *, ln_std_map=None, **kwds):
"""Use generic model parameters based on site class.
Parameters
----------
site_class: str
Site classification. Possible options are:
* Geomatrix AB
* Geomatrix CD
* USGS AB
* USGS CD
* USGS A
* USGS B
* USGS C
* USGS D
See the report for definitions of the Geomatrix site
classication. USGS site classification is based on :math:`V_{s30}`:
=========== =====================
Site Class :math:`V_{s30}` (m/s)
=========== =====================
A >750 m/s
B 360 to 750 m/s
C 180 to 360 m/s
D <180 m/s
=========== =====================
ln_std_map: dict[str, float], optional
Mapping between the soil_type and the defined ln_std. Default is *None*.
Returns
-------
DepthAndSoilTypeDependToroVelVariation
Initialized :class:`DepthDependToroVelVariation` with generic parameters.
"""
p = dict(cls.PARAMS[site_class])
p.update(kwds)
if "depth" not in kwds:
p["depth"] = [0, 15]
p["ln_std"] = [0.25, 0.15]
p["ln_std_map"] = ln_std_map
return cls(**p)
[docs]
class SoilTypeVariation:
"""Base class for soil-type (modulus-reduction and damping) variation.
Parameters
----------
correlation : float
Correlation coefficient between the modulus-reduction and damping
random variables.
limits_mod_reduc : list[float], optional
``[min, max]`` clipping bounds for modulus reduction.
limits_damping : list[float], optional
``[min, max]`` clipping bounds for damping.
vary_bedrock : bool, optional
Whether to include the half-space in the variation.
sample_mode : {'random', 'fixed_percentiles'}, optional
How samples are drawn when :meth:`iter_varied_profiles` iterates:
* ``'random'`` *(default)* — each call draws an independent
random realisation from the truncated-normal distribution.
* ``'fixed_percentiles'`` — each call uses a pre-specified
percentile supplied via ``sample_index`` so that the same
index always produces an identical realisation. The percentile
is selected as ``percentiles[sample_index % len(percentiles)]``,
so the list cycles when *count* is a multiple of its length.
The caller must pass the ``sample_index`` keyword argument to
:meth:`__call__`, and :func:`iter_varied_profiles` does this
automatically.
percentiles : list[float] | None, optional
Ordered sequence of quantiles in ``(0, 1)`` used in
``'fixed_percentiles'`` mode. The sequence cycles: iteration *i*
draws ``percentiles[i % len(percentiles)]``. Required (and only
used) when *sample_mode* is ``'fixed_percentiles'``.
"""
[docs]
def __init__(
self,
correlation,
limits_mod_reduc=[0.05, 1],
limits_damping=[0, 0.15],
vary_bedrock=False,
sample_mode="random",
percentiles=None,
):
_valid_modes = {"random", "fixed_percentiles"}
if sample_mode not in _valid_modes:
raise ValueError(
f"sample_mode must be one of {_valid_modes!r}, got {sample_mode!r}"
)
if sample_mode == "fixed_percentiles":
if percentiles is None or len(percentiles) == 0:
raise ValueError(
"percentiles must be a non-empty sequence when "
"sample_mode='fixed_percentiles'"
)
percentiles = list(percentiles)
if any(not (0.0 < p < 1.0) for p in percentiles):
raise ValueError(
"All percentile values must be strictly between 0 and 1"
)
self._vary_bedrock = vary_bedrock
self._correlation = correlation
self._limits_mod_reduc = list(limits_mod_reduc)
self._limits_damping = list(limits_damping)
self._sample_mode = sample_mode
self._percentiles = percentiles
[docs]
def __call__(self, soil_type, sample_index=None):
"""Return a single varied realisation of *soil_type*.
Parameters
----------
soil_type : site.SoilType
The nominal (seed) soil type to vary.
sample_index : int | None, optional
Index into :attr:`percentiles` used when
``sample_mode='fixed_percentiles'``. Ignored in
``'random'`` mode. Must be provided (and within range) in
``'fixed_percentiles'`` mode.
"""
def get_values(nlp):
try:
return nlp.values
except AttributeError:
return np.asarray(nlp).astype(float)
mod_reduc = get_values(soil_type.mod_reduc)
damping = get_values(soil_type.damping)
# A pair of correlated random variables
if self._sample_mode == "fixed_percentiles":
if sample_index is None:
raise ValueError(
"sample_index must be provided when sample_mode='fixed_percentiles'"
)
percentile = self._percentiles[sample_index]
randvar = randnorm.correlated_at_percentile(self.correlation, percentile)
else:
randvar = randnorm.correlated(self.correlation)
varied_mod_reduc, varied_damping = self._get_varied(randvar, mod_reduc, damping)
# Clip the values to the specified min/max
varied_mod_reduc = np.clip(
varied_mod_reduc, self.limits_mod_reduc[0], self.limits_mod_reduc[1]
)
varied_damping = np.clip(
varied_damping, self.limits_damping[0], self.limits_damping[1]
)
# Set the values
realization = copy.deepcopy(soil_type)
for attr_name, values in zip(
["mod_reduc", "damping"], [varied_mod_reduc, varied_damping]
):
try:
getattr(realization, attr_name).values = values
except AttributeError:
setattr(realization, attr_name, values)
return realization
[docs]
def vary_profile(self, profile: site.Profile, sample_index: int | None = None):
"""Return a profile with varied soil types.
Parameters
----------
profile : site.Profile
Input profile to vary.
sample_index : int | None, optional
Index into :attr:`percentiles` for ``sample_mode='fixed_percentiles'``.
Ignored in ``'random'`` mode.
"""
# Map of varied soil types
varied = {
str(st): self(st, sample_index=sample_index)
for st in profile.iter_soil_types()
}
# Create new layers
end = None if self.vary_bedrock else -1
layers = [
site.Layer(
varied[str(layer.soil_type)],
layer.thickness,
layer.initial_shear_vel,
layer.damping_min,
)
for layer in profile[:end]
]
# Add the unrandomized bedrock
if not self.vary_bedrock:
layers.append(profile[-1])
return site.Profile(layers, profile.wt_depth)
def _get_varied(self, randvar, mod_reduc, damping):
raise NotImplementedError
@property
def correlation(self):
return self._correlation
@property
def limits_damping(self):
return self._limits_damping
@property
def limits_mod_reduc(self):
return self._limits_mod_reduc
@property
def vary_bedrock(self):
return self._vary_bedrock
@property
def sample_mode(self):
"""Sampling mode: ``'random'`` or ``'fixed_percentiles'``."""
return self._sample_mode
@property
def percentiles(self):
"""Percentile list used in ``'fixed_percentiles'`` mode, or *None*."""
return self._percentiles
class DarendeliVariation(SoilTypeVariation):
def _get_varied(self, randvar, mod_reduc, damping):
mod_reduc_means = mod_reduc
mod_reduc_stds = self.calc_std_mod_reduc(mod_reduc_means)
varied_mod_reduc = mod_reduc_means + randvar[0] * mod_reduc_stds
damping_means = damping
damping_stds = self.calc_std_damping(damping_means)
varied_damping = damping_means + randvar[1] * damping_stds
return varied_mod_reduc, varied_damping
@staticmethod
def calc_std_mod_reduc(mod_reduc):
"""Calculate the standard deviation as a function of G/G_max.
Equation 7.29 from Darendeli (2001).
Parameters
----------
mod_reduc : array_like
Modulus reduction values.
Returns
-------
std : :class:`numpy.ndarray`
Standard deviation.
"""
mod_reduc = np.asarray(mod_reduc).astype(float)
std = np.exp(-4.23) + np.sqrt(
0.25 / np.exp(3.62) - (mod_reduc - 0.5) ** 2 / np.exp(3.62)
)
return std
@staticmethod
def calc_std_damping(damping):
"""Calculate the standard deviation as a function of damping in decimal.
Equation 7.30 from Darendeli (2001).
Parameters
----------
damping : array_like
Material damping values in decimal.
Returns
-------
std : :class:`numpy.ndarray`
Standard deviation.
"""
damping = np.asarray(damping).astype(float)
std = (np.exp(-5) + np.exp(-0.25) * np.sqrt(100 * damping)) / 100.0
return std
[docs]
class SpidVariation(SoilTypeVariation):
"""Variation defined by the EPRI SPID (2013) and documented in PNNL (2014).
EPRI SPID (2013): https://www.nrc.gov/docs/ML1233/ML12333A170.pdf
"""
[docs]
def __init__(
self,
correlation,
limits_mod_reduc=[0, 1],
limits_damping=[0, 0.15],
std_mod_reduc=0.15,
std_damping=0.30,
sample_mode="random",
percentiles=None,
):
super().__init__(
correlation,
limits_mod_reduc,
limits_damping,
sample_mode=sample_mode,
percentiles=percentiles,
)
self._std_mod_reduc = std_mod_reduc
self._std_damping = std_damping
def _get_varied(self, randvar, mod_reduc, damping):
# Vary the G/Gmax in transformed space.
# PNNL (2014) Hanford Site Wide Hazard Study is available here:
# https://www.hanford.gov/files.cfm/00_Front_Matter.pdf
# https://www.hanford.gov/files.cfm/9.0_Ground_Motion_Characterization.pdf
# Equation 9.43 of PNNL (2014)
# Here epsilon is added so that the denomiator doesn't go to zero.
f_mean = mod_reduc / (1 - mod_reduc + np.finfo(float).eps)
# Instead of constraining the standard deviation at a specific
# strain, then standard deviation is constrained at G/Gmax of 0.5.
# This is modified from Equation 9.44 of PNNL (2014).
f_std = self.std_mod_reduc * (1 / (1 - 0.5))
f_real = np.exp(randvar[0] * f_std) * f_mean
# Equation 9.45 of PNNL (2014)
varied_mod_reduc = f_real / (1 + f_real)
# Simple log distribution
varied_damping = np.exp(randvar[1] * self.std_damping) * damping
return varied_mod_reduc, varied_damping
@property
def std_damping(self):
return self._std_damping
@property
def std_mod_reduc(self):
return self._std_mod_reduc
def iter_varied_profiles(
profile: site.Profile,
count: int,
var_thickness: ToroThicknessVariation | None = None,
var_velocity: VelocityVariation | None = None,
var_soiltypes: SoilTypeVariation | None = None,
) -> Generator[site.Profile]:
"""Iterate over simulated profiles.
Parameters
----------
profile : site.Profile
seed profile
count : int
number of interations
var_thickness : ToroThicknessVariation | None
model for the thickness variation
var_velocity : VelocityVariation | None
model for the velocity variation
var_soiltypes : SoilTypeVariation | None
model for the soil type variation. When ``sample_mode='fixed_percentiles'``
*count* must be divisible by the number of configured percentiles; the
percentile list cycles so that each percentile is used
``count // len(percentiles)`` times in order.
Returns
-------
site.Profile
varied profile
"""
if var_soiltypes and var_soiltypes.sample_mode == "fixed_percentiles":
n_pct = len(var_soiltypes.percentiles)
if count % n_pct != 0:
raise ValueError(
f"count ({count}) must be divisible by the number of percentiles "
f"({n_pct}) when sample_mode='fixed_percentiles'"
)
for i in range(count):
# Copy the profile to form the realization
_profile = profile.copy()
if var_thickness:
_profile = var_thickness(_profile)
if var_velocity:
_profile = var_velocity(_profile)
if var_soiltypes:
# Determine the sample_index to forward (None in random mode)
n_pct = (
len(var_soiltypes.percentiles)
if var_soiltypes.sample_mode == "fixed_percentiles"
else 0
)
sample_index = (
i % n_pct if var_soiltypes.sample_mode == "fixed_percentiles" else None
)
_profile = var_soiltypes.vary_profile(_profile, sample_index=sample_index)
yield _profile