diff --git a/pyrvt/motions.py b/pyrvt/motions.py
index cad31d4..ff9768a 100644
--- a/pyrvt/motions.py
+++ b/pyrvt/motions.py
@@ -20,10 +20,32 @@
"""
Classes and functions used to define random vibration theory (RVT) based
motions.
+
+References
+----------
+.. [AH84] Anderson, J. G., & Hough, S. E. (1984). A model for the shape of the
+ Fourier amplitude spectrum of acceleration at high frequencies. Bulletin of
+ the Seismological Society of America, 74(5), 1969-1993.
+
+.. [AB11] Atkinson, G. M., & Boore, D. M. (2011). Modifications to existing
+ ground-motion prediction equations in light of new data. Bulletin of the
+ Seismological Society of America, 101(3), 1121-1135.
+
+.. [C03] Campbell, K. W. (2003). Prediction of strong ground motion using
+ the hybrid empirical method and its use in the development of ground-motion
+ (attenuation) relations in eastern North America. Bulletin of the
+ Seismological Society of America, 93(3), 1012-1033.
+
+.. [GV76] Gasparini, D. A., & Vanmarcke, E. H. (1976). Simulated earthquake
+ motions compatible with prescribed response spectra. Massachusetts
+ Institute of Technology, Department of Civil Engineering, Constructed
+ Facilities Division.
+
"""
import numpy as np
+from scipy.stats import linregress
from scipy.interpolate import interp1d
from . import peak_calculators
@@ -31,6 +53,75 @@
DEFAULT_CALC = 'V75'
+def increasing_x(x, y=None):
+ """Check if *x* is monotonically increasing. If not, reverse it.
+
+ Parameters
+ ----------
+ x : :class:`numpy.array`
+ X values, which are checked.
+
+ y : :class:`numpy.array` or ``None``, default: ``None``
+ Y values, which are reversed if *x* is reversed
+
+ Returns
+ -------
+ If *y* is not ``None``, then returns (*x*, *y*), otherwise only returns *x*.
+ x : :class:`numpy.array`
+ X values in monotonically increasing order.
+
+ y : :class:`numpy.array`
+ Y values in same order as *x*.
+
+ Raises
+ ------
+ :class:`NotImplementedError`
+ If *x* is not monotonic.
+ """
+
+ diffs = np.diff(x)
+ if np.all(0 <= diffs):
+ # All increasing, do nothing
+ pass
+ elif np.all(diffs <= 0):
+ # All decreasing, reverse
+ x = x[::-1]
+ if y is not None:
+ y = y[::-1]
+ else:
+ raise NotImplementedError('Values are not regularly ordered.')
+
+ if y is None:
+ return x
+ else:
+ return x, y
+
+
+def log_spaced_values(min, max):
+ """Generate values with constant log-spacing.
+
+ Values are generated with 512 points per decade.
+
+ Parameters
+ ----------
+ min : float
+ Minimum value of the range.
+
+ max : float
+ Maximum value of the range.
+
+ Returns
+ -------
+ value : :class:`numpy.array`
+ Log-spaced values
+
+ """
+ lower = np.log10(min)
+ upper = np.log10(max)
+ count = np.ceil(512 * (upper - lower))
+ return np.logspace(lower, upper, count)
+
+
def compute_sdof_tf(freqs, osc_freq, osc_damping):
"""Compute the single-degree-of-freedom transfer function. When applied on
the acceleration Fourier amplitude spectrum, it provides the
@@ -61,12 +152,6 @@ def compute_sdof_tf(freqs, osc_freq, osc_damping):
def compute_stress_drop(magnitude):
"""Compute the stress drop using [AB11]_ model.
- References
- ----------
- .. [AB11] Atkinson, G. M., & Boore, D. M. (2011). Modifications to existing
- ground-motion prediction equations in light of new data. Bulletin of
- the Seismological Society of America, 101(3), 1121-1135.
-
Parameters
----------
magnitude : float
@@ -121,13 +206,13 @@ class RvtMotion(object):
Parameters
----------
- freqs : :class:`numpy.array`
+ freqs : :class:`numpy.array` or ``None``, default: ``None``
Frequency array [Hz]
- fourier_amps : :class:`numpy.array`
+ fourier_amps : :class:`numpy.array` or ``None``, default: ``None``
Absolute value of acceleration Fourier amplitudes.
- duration : float
+ duration : float or ``None``, default: ``None``
Ground motion duration [dec].
peak_calculator : str or :class:`~.peak_calculators.Calculator`, default: ``None``
@@ -147,6 +232,10 @@ def __init__(self, freqs=None, fourier_amps=None, duration=None,
self._fourier_amps = fourier_amps
self._duration = duration
+ if self._freqs is not None:
+ self._freqs, self._fourier_amps = increasing_x(
+ self._freqs, self._fourier_amps)
+
if isinstance(peak_calculator, peak_calculators.Calculator):
self.peak_calculator = peak_calculator
else:
@@ -161,7 +250,7 @@ def freqs(self):
@property
def fourier_amps(self):
"""Acceleration Fourier amplitude values."""
- return self._freqs
+ return self._fourier_amps
@property
def duration(self):
@@ -222,17 +311,67 @@ def compute_peak(self, transfer_func=None, osc_freq=None,
osc_freq=osc_freq, osc_damping=osc_damping,
full_output=False)
+ def compute_attenuation(self, min_freq, max_freq=None, full=False):
+ """Compute the site attenuation (κ) based on a log-linear fit.
+
+ This function computes the site attenuation defined by [AH84]_ as:
+
+ .. math::
+ a(f) = A_0 \exp(-\pi \kappa f) \\text( for ) f > f_E
+
+ for a single Fourier amplitude spectrum
+
+
+ Parameters
+ ----------
+ min_freq : float
+ minimum frequency of the fit
+
+ max_freq : float or ``None``, default: ``None``
+ maximum frequency of the fit. If ``None``, then the maximum
+ frequency range is used.
+
+ full : bool, default: ``False``
+ If the complete output should be returned
+
+ Returns
+ -------
+ If *full* is ``False``, then only *atten* is returned. Otherwise,
+ the tuple (*atten*, *r_value*, *freqs*, *fitted*) is return.
+
+ atten : float
+ attenuation parameter
+
+ r_sqrt : float
+ sqared correlation coefficient of the fit (R²). See
+ :function:`scipy.stats.linregress`.
+
+ freqs : :class:`numpy.array`
+ selected frequencies
+
+ fitted : :class:`numpy.array`
+ fitted values
+
+ """
+ max_freq = max_freq or self.freqs[-1]
+ mask = (min_freq <= self.freqs) & (self.freqs <= max_freq)
+
+ slope, intercept, r_value, p_value, stderr = linregress(
+ self.freqs[mask], np.log(self.fourier_amps[mask]))
+
+ atten = slope / -np.pi
+
+ if full:
+ freqs = self.freqs[mask]
+ fitted = np.exp(intercept + slope * freqs)
+ return atten, r_value ** 2, freqs, fitted
+ else:
+ return atten
+
class SourceTheoryMotion(RvtMotion):
"""Single-corner source theory model with default parameters from [C03]_.
- References
- ----------
- .. [C03] Campbell, K. W. (2003). Prediction of strong ground motion using
- the hybrid empirical method and its use in the development of
- ground-motion (attenuation) relations in eastern North America.
- Bulletin of the Seismological Society of America, 93(3), 1012-1033.
-
Parameters
----------
magnitude : float
@@ -367,17 +506,21 @@ def compute_duration(self):
return duration_source + duration_path
- def compute_fourier_amps(self, freqs):
+ def compute_fourier_amps(self, freqs=None):
"""Compute the acceleration Fourier amplitudes for a frequency range.
Parameters
----------
- freqs : numpy.array
- Frequency range
+ freqs : :class:`numpy.array` or ``None``, default: ``None``
+ Frequency range. If no frequency range is specified then
+ :function:`log_spaced_values(0.05, 200.)` is used.
"""
+ if freqs is None:
+ self._freqs = log_spaced_values(0.05, 200.)
+ else:
+ self._freqs = increasing_x(np.asarray(freqs))
- self._freqs = np.asarray(freqs)
self._duration = self.compute_duration()
# Model component
@@ -467,13 +610,8 @@ def __init__(self, osc_freqs, osc_accels_target, duration=None,
super(CompatibleRvtMotion, self).__init__(
peak_calculator=peak_calculator)
- osc_freqs = np.asarray(osc_freqs)
- osc_accels_target = np.asarray(osc_accels_target)
-
- # Order by increasing frequency
- ind = osc_freqs.argsort()
- osc_freqs = osc_freqs[ind]
- osc_accels_target = osc_accels_target[ind]
+ osc_freqs, osc_accels_target = increasing_x(
+ np.asarray(osc_freqs), np.asarray(osc_accels_target))
if duration:
self._duration = duration
@@ -485,12 +623,10 @@ def __init__(self, osc_freqs, osc_accels_target, duration=None,
osc_freqs, osc_accels_target, osc_damping)
# The frequency needs to be extended to account for the fact that the
- # osciallator transfer function has a width. The number of frequencies
+ # oscillator transfer function has a width. The number of frequencies
# depends on the range of frequencies provided.
- lower = np.log10(osc_freqs[0] / 2.)
- upper = np.log10(2 * osc_freqs[-1])
- count = int(512 * (upper - lower))
- self._freqs = np.logspace(lower, upper, count)
+ self._freqs = log_spaced_values(osc_freqs[0] / 2.,
+ 2. * osc_freqs[-1])
self._fourier_amps = np.empty_like(self._freqs)
# Indices of the first and last point with the range of the provided
@@ -501,7 +637,6 @@ def __init__(self, osc_freqs, osc_accels_target, duration=None,
# last is extend one past the usable range to allow use of first:last
# notation
last = indices[-1, 0] + 1
-
log_freqs = np.log(self._freqs)
log_osc_freqs = np.log(osc_freqs)
@@ -581,7 +716,7 @@ def _extrap(freq, freqs, fourier_amps, max_slope=None):
self.iterations += 1
def _estimate_fourier_amps(self, osc_freqs, osc_accels, osc_damping):
- """Compute an estimate of the FAS using the Vanmarcke methodology.
+ """Compute an estimate of the FAS using the [GV76]_ methodology.
Parameters
----------
@@ -601,10 +736,9 @@ def _estimate_fourier_amps(self, osc_freqs, osc_accels, osc_damping):
frequencies specifed by *osc_freqs*.
"""
- # TODO Add reference to the docstring
# Compute initial value using Vanmarcke methodology. The response is
- # first computed at the lowest frequency and then subsequently comptued
+ # first computed at the lowest frequency and then subsequently computed
# at higher frequencies.
peak_factor = 2.5
diff --git a/pyrvt/peak_calculators.py b/pyrvt/peak_calculators.py
index 0693812..d75e4fa 100644
--- a/pyrvt/peak_calculators.py
+++ b/pyrvt/peak_calculators.py
@@ -852,6 +852,7 @@ def get_peak_calculator(method, calc_kwds):
calculator : :class:`.Calculator`
"""
+ calc_kwds = calc_kwds or dict()
calculators = [
BooreJoyner1984,
diff --git a/pyrvt/tests/test_motions.py b/pyrvt/tests/test_motions.py
index ec4147d..b6a74f0 100644
--- a/pyrvt/tests/test_motions.py
+++ b/pyrvt/tests/test_motions.py
@@ -1,70 +1,96 @@
-#!/usr/bin/env python3
-# encoding: utf-8
-
-# pyRVT: Seismological random vibration theory implemented with Python
-# Copyright (C) 2013-2014 Albert R. Kottke albert.kottke@gmail.com
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program. If not, see .
-
-import numpy as np
-from numpy.testing import assert_allclose
-
-from .. import motions
-from .. import peak_calculators
-
-import matplotlib.pyplot as plt
-
-def test_compatible_rvt_motion():
- # Compute the target from the point source model.
- target = motions.SourceTheoryMotion(
- 6., 20., 'wna',
- peak_calculator=peak_calculators.DerKiureghian1985())
- target.compute_fourier_amps(np.logspace(-1.5, 2, 1024))
-
- osc_freqs = np.logspace(-1, 2, num=50)
- osc_accels_target = target.compute_osc_accels(osc_freqs, 0.05)
-
- compat = motions.CompatibleRvtMotion(
- osc_freqs, osc_accels_target,
- duration=target.duration, osc_damping=0.05,
- peak_calculator=peak_calculators.DerKiureghian1985())
-
- osc_accels_compat = compat.compute_osc_accels(osc_freqs, 0.05)
-
- # Might be off by a few percent because of difficulties with the inversion.
- assert_allclose(osc_accels_target, osc_accels_compat, rtol=0.03, atol=0.05)
-
- # fig, axes = plt.subplots(2, 1)
- #
- # ax = axes.flat[0]
- #
- # ax.plot(target.freqs, target.fourier_amps, 'b-', label='Target')
- # ax.plot(compat.freqs, compat.fourier_amps, 'r--', label='Compatible')
- #
- # ax.set_xlabel('Frequency (Hz)')
- # ax.set_xscale('log')
- # ax.set_ylabel('Fourier Ampl. (cm/s)')
- # ax.set_yscale('log')
- #
- # ax = axes.flat[1]
- #
- # ax.plot(osc_freqs, osc_resp_target, 'b-', label='Target')
- # ax.plot(osc_freqs, osc_resp_compat, 'r--', label='Compatible')
- #
- # ax.set_xlabel('Frequency (Hz)')
- # ax.set_xscale('log')
- # ax.set_ylabel('Spectral Accel. (cm/s²)')
- # ax.set_yscale('log')
- #
- # fig.savefig('compatible_fas.png', dpi=300)
+#!/usr/bin/env python3
+# encoding: utf-8
+
+# pyRVT: Seismological random vibration theory implemented with Python
+# Copyright (C) 2013-2014 Albert R. Kottke albert.kottke@gmail.com
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see .
+
+import numpy as np
+from numpy.testing import assert_allclose
+
+from .. import motions
+from .. import peak_calculators
+
+import matplotlib.pyplot as plt
+
+def test_compute_attenuation():
+ m = motions.SourceTheoryMotion(5.5, 0, 'cena', depth=1)
+ m.compute_fourier_amps()
+
+ atten = m.compute_attenuation(50)
+ assert_allclose(0.006, atten, rtol=0.01)
+
+def test_compute_attenuation_full():
+ m = motions.SourceTheoryMotion(5.5, 0, 'cena', depth=1)
+ m.compute_fourier_amps()
+
+ atten, r_value, freqs, fitted = m.compute_attenuation(50, full=True)
+
+ assert_allclose(0.006, atten, rtol=0.01)
+ assert_allclose(1.0, r_value, rtol=0.01)
+
+ # fig = plt.figure()
+ # ax = fig.add_subplot(111)
+ # ax.plot(m.freqs, m.fourier_amps, 'b-')
+ # ax.set_xlabel('Frequency (Hz)')
+ # ax.set_xscale('log')
+ # ax.set_ylabel('Amplitude')
+ # ax.set_yscale('log')
+ #
+ # fig.savefig('test')
+
+def test_compatible_rvt_motion():
+ # Compute the target from the point source model.
+ target = motions.SourceTheoryMotion(
+ 6., 20., 'wna',
+ peak_calculator=peak_calculators.DerKiureghian1985())
+ target.compute_fourier_amps(np.logspace(-1.5, 2, 1024))
+
+ osc_freqs = np.logspace(-1, 2, num=50)
+ osc_accels_target = target.compute_osc_accels(osc_freqs, 0.05)
+
+ compat = motions.CompatibleRvtMotion(
+ osc_freqs, osc_accels_target,
+ duration=target.duration, osc_damping=0.05,
+ peak_calculator=peak_calculators.DerKiureghian1985())
+
+ osc_accels_compat = compat.compute_osc_accels(osc_freqs, 0.05)
+
+ # Might be off by a few percent because of difficulties with the inversion.
+ assert_allclose(osc_accels_target, osc_accels_compat, rtol=0.03, atol=0.05)
+
+ # fig, axes = plt.subplots(2, 1)
+ #
+ # ax = axes.flat[0]
+ #
+ # ax.plot(target.freqs, target.fourier_amps, 'b-', label='Target')
+ # ax.plot(compat.freqs, compat.fourier_amps, 'r--', label='Compatible')
+ #
+ # ax.set_xlabel('Frequency (Hz)')
+ # ax.set_xscale('log')
+ # ax.set_ylabel('Fourier Ampl. (cm/s)')
+ # ax.set_yscale('log')
+ #
+ # ax = axes.flat[1]
+ #
+ # ax.plot(osc_freqs, osc_resp_target, 'b-', label='Target')
+ # ax.plot(osc_freqs, osc_resp_compat, 'r--', label='Compatible')
+ #
+ # ax.set_xlabel('Frequency (Hz)')
+ # ax.set_xscale('log')
+ # ax.set_ylabel('Spectral Accel. (cm/s²)')
+ # ax.set_yscale('log')
+ #
+ # fig.savefig('compatible_fas.png', dpi=300)
diff --git a/pyrvt/tools.py b/pyrvt/tools.py
index 893f730..a888d29 100644
--- a/pyrvt/tools.py
+++ b/pyrvt/tools.py
@@ -18,9 +18,7 @@
# along with this program. If not, see .
"""
-File: tools.py
-Author: Albert Kottke
-Description: Tools for reading/writing of files and performing operations.
+Tools for reading/writing of files and performing operations.
"""
import csv