# Organization of complex empirical equations

I am working on coding up a set of models that are based on empirical models. These models tend to be quite complicated. In some cases, the logic used by the equations makes it impossible to compute the results using numpy fast arrays. As a result, I have decided to use list comprehension to do the math. However, this results in the functions being indent 3 levels. Is this bad? Is there a better way to organize this code? Any other comments on the style or organization would be most welcome.

Currently, the code is organized with classes, which I think makes the most sense. Any opinions on this would be welcome too.

Here are a few examples of the code:

import logging
import math

import numpy as np

class Model(object):
'''Abstract class for ground motion prediction models.'''

INDICES_PSA = []
INDEX_PGA = None
INDEX_PGV = None
INDEX_PGD = None

LIMITS = dict()

def __init__(self, name, abbrev, **kwds):
'''Compute the response predicted the model. No default implementation.

Input names:
mag: float
moment magnitude of the event
dist_closest: float
closest distance to the rupture (km)
v_s30: float
time-averaged shear-wave velocity over the top 30 m of the site (m/s)
'''
super(Model, self).__init__()

self.name = name
self.abbrev = abbrev
self.parameters = kwds

self._ln_resp = None
self._ln_std = None

self._check_inputs(**self.parameters)

def interp_spec_accels(self, periods):
return np.interp(periods,
self.periods, self._ln_resp[self.INDICES_PSA])

@property
def periods(self):
return self.PERIODS[self.INDICES_PSA]

@property
def spec_accels(self):
return self._resp(self.INDICES_PSA)

@property
def pga(self):
'''Return the PGA [g]'''
return self._resp(self.INDEX_PGA)

@property
def pgv(self):
'''Return the PGV [cm/s]'''
return self._resp(self.INDEX_PGV, FROM_GRAVITY)

@property
def pgd(self):
'''Return the PGD [cm]'''
return self._resp(self.INDEX_PGD, FROM_GRAVITY)

def _resp(self, index, scale=1):
if not index is None:
return np.exp(self._ln_resp[index]) * scale

def _check_inputs(self, **kwds):
for key, (minimum, maximum) in self.LIMITS.items():
if not key in kwds:
continue

value = kwds[key]
if not minimum is None and value < minimum:
logging.warning('{} ({:g}) is less than the recommended limit ({:g}).'.format(
key, minimum, value))
elif not maximum is None and maximum < value:
logging.warning('{} ({:g}) is greater than the recommended limit ({:g}).'.format(
key, maximum, value))

fname = os.path.join(os.path.dirname(__file__), 'data', name)

class ChiouYoungs2013(model.Model):
'''Chiou and Youngs (2013) model.'''

# Reference velocity (m/s)
V_REF = 1130.

# Load the coefficients for the model

PERIODS = COEFF['period']

# Period independent coefficients
COEFF_C_2 = 1.06
COEFF_C_4 = -2.1
COEFF_C_4a = -0.5
COEFF_C_RB = 50
COEFF_C_8 = 0.2153
COEFF_C_8a = 0.2695

INDICES_PSA = np.arange(23)

LIMITS = dict(
dist_closest=(0, 300),
mag=(3.5, 8.5),
v_s30=(180, 1500),
z_tor=(0, 20),
)

def __init__(self, **kwds):
'''Compute the response predicted the Chiou and Youngs (2013) ground
motion model.

Input:
depth_1_0_centered: float
depth (m) to a shear-wave velocity of 1,000 m/sec centered to
the California average
depth_tor: float
Depth to the top of the rupture (km)
dist_jb: float
Joyner-Boore distance to the rupture plane (km)
dist_closest: float
closest distance to the rupture (km)
dist_x: float
site coordinate (km) measured perpendicular to the fault strike
from the fault line with the down-dip direction being positive
(see Figure 3.12 in Chiou and Youngs (2008).
dpp_centered: float
direct point parameter for directivity effect (see Chiou and
Spudich (2013)) centered on the earthquake-specific average DPP
for California. A value of 0 provides the average directivity.
dip_angle: float
Fault dip angle (deg)
flag_hw: int
Hanging-wall flag.
1: R_x >= 0
0: R_x < 0
flag_nm: int
Normal faulting flag.
1: -120° <= rake angle <= -60° (excludes normal-oblique)
0: otherwise
flag_rv: int
Reverse-faulting flag.
1: 30° <= rake angle <= 150° (combined reverse and
reverse-oblique)
0: otherwise
mag: float
moment magnitude of the event
v_s30: float
time-averaged shear-wave velocity over the top 30 m of the site
(m/s)
'''

dist_jb = kwds['dist_jb']
dist_closest = kwds['dist_closest']
dist_x = kwds['dist_x']
dip_angle = kwds['dip_angle']
flag_hw = kwds['flag_hw']
flag_nm = kwds['flag_nm']
flag_rv = kwds['flag_rv']
mag = kwds['mag']
v_s30 = kwds['v_s30']

dpp_centered = kwds.get('dpp_centered', 0.)
depth_tor = kwds.get('depth_tor',
self.calc_depth_tor(mag, flag_rv))
depth_1_0_centered = kwds.get('depth_1_0_centered',
self.calc_depth_1_0(v_s30, 'california'))

def calc_ln_resp_ref(period, c_1, c_1a, c_1b, c_1c, c_1d, c_n, c_m,
c_3, c_5, c_hm, c_6, c_7, c_7b, c_8b, c_9, c_9a,
c_9b, c_11, c_11b, c_gamma1, c_gamma2, c_gamma3,
phi_1, phi_2, phi_3, phi_4, phi_5, phi_6, tau_1,
tau_2, sigma_1, sigma_2, sigma_3):
cosh_mag = math.cosh(2 * max(mag - 4.5, 0))
return (
c_1
+ (c_1a + c_1c / cosh_mag) * flag_rv
+ (c_1b + c_1d / cosh_mag) * flag_nm
+ (c_7 + c_7b / cosh_mag) * depth_tor
+ (c_11 + c_11b / cosh_mag) * math.cos(dip_angle) ** 2
+ self.COEFF_C_2 * (mag - 6)
+ (self.COEFF_C_2 - c_3) / c_n * math.log(1 + math.exp(c_n * (c_m - mag)))
+ self.COEFF_C_4 * math.log(dist_closest + c_5 * math.cosh(c_6 * max(mag - c_hm, 0)))
+ (self.COEFF_C_4a - self.COEFF_C_4) * math.log(math.sqrt(dist_closest ** 2 + self.COEFF_C_RB ** 2))
+ (c_gamma1 + c_gamma2 / math.cosh(max(mag - c_gamma3, 0))) * dist_closest
+ (self.COEFF_C_8 * max(1 - max(dist_closest - 40, 0) / 30, 0) * min(max(mag - 5.5, 0) / 0.8, 1)
* math.exp(-self.COEFF_C_8a * (mag - c_8b) ** 2) * dpp_centered)
+ (c_9 * flag_hw * math.cos(dip_angle) * (c_9a + (1 - c_9a) * math.tanh(dist_x / c_9b))
* (1 - math.sqrt(dist_jb ** 2 + depth_tor ** 2) / (dist_closest + 1)))
)

def calc_ln_resp(ln_resp_ref, period, c_1, c_1a, c_1b, c_1c, c_1d, c_n,
c_m, c_3, c_5, c_hm, c_6, c_7, c_7b, c_8b, c_9, c_9a,
c_9b, c_11, c_11b, c_gamma1, c_gamma2, c_gamma3,
phi_1, phi_2, phi_3, phi_4, phi_5, phi_6, tau_1,
tau_2, sigma_1, sigma_2, sigma_3):
return (ln_resp_ref
+ phi_1 * min(math.log(v_s30 / 1130.), 0)
+ (phi_2 * (math.exp(phi_3 * (min(v_s30, 1130.) - 360.)) -
math.exp(phi_3 * (1130. - 360.)))
* math.log((math.exp(ln_resp_ref) + phi_4) / phi_4))
+ phi_5 * (1 - math.exp(-depth_1_0_centered / phi_6)))

def calc_ln_std(ln_resp_ref, period, c_1, c_1a, c_1b, c_1c, c_1d, c_n,
c_m, c_3, c_5, c_hm, c_6, c_7, c_7b, c_8b, c_9, c_9a,
c_9b, c_11, c_11b, c_gamma1, c_gamma2, c_gamma3, phi_1,
phi_2, phi_3, phi_4, phi_5, phi_6, tau_1, tau_2,
sigma_1, sigma_2, sigma_3):
between = tau_1 + (tau_2 - tau_1) / 2.25 * (min(max(mag, 5.), 7.25) - 5.)

nl_0 = (phi_2 * (math.exp(phi_3 * (min(v_s30, 1130.) - 360.))
- math.exp(phi_3 * (1130. - 360.)))
* (ln_resp_ref / (ln_resp_ref + phi_4)))

# FIXME
logging.warning('CY do not specify what f_inf and f_meas should be in their equation')
f_inf = 1
f_meas = 1

within_nl = (sigma_1
+ (sigma_2 - sigma_1) / 2.25 * (min(max(mag, 5.), 7.25) - 5.) * math.sqrt(
sigma_3 * f_inf + 0.7 * f_meas + (1 + nl_0) ** 2))

return math.sqrt((1 - nl_0) ** 2 * between ** 2 + within_nl ** 2)

ln_resp_ref = [calc_ln_resp_ref(*c) for c in self.COEFF]

self._ln_resp = np.array([calc_ln_resp(lrr, *c) for (lrr, c) in zip(ln_resp_ref, self.COEFF)])
self._ln_std = np.array([calc_ln_std(lrr, *c) for (lrr, c) in zip(ln_resp_ref, self.COEFF)])

def _check_inputs(self, **kwds):
super(ChiouYoungs2013, self)._check_inputs(**kwds)

assert flag_nm in [0, 1]
assert flag_rv in [0, 1]

if flag_nm or flag_rv:
_min, _max = 3.5, 8.0
if not (_min <= kwds['mag'] <= _max):
logging.warning(
'Magnitude ({}) exceeds recommened bounds ({} to {})'\
' for a reverse or normal earthquake!'.format(
kwds['mag'], _min, _max))
else:
_min, _max = 3.5, 8.5
if not (_min <= kwds['mag'] <= _max):
logging.warning(
'Magnitude ({}) exceeds recommened bounds ({} to {})'\
' for a strike-slip earthquake!'.format(
kwds['mag'], _min, _max))

def calc_depth_1_0(self, v_s30, region='california'):
'''Calculate an estimate of the depth to 1 km/sec based on Vs30 and region.

Input:
v_s30: float
time-averaged shear-wave velocity over the top 30 m of the site
(m/s)
region: str
region for the V_s30-Z_1.0 correlation. Potential regions:
california (default)
japan
Returns:
depth_1_0_centered: float
depth (m) to a shear-wave velocity of 1,000 m/sec.

'''
assert region in ['california', 'japan']

if region == 'california':
power = 4
bar = 571
foo = -7.15 / power
elif region == 'japan':
power = 2
bar = 412
foo = -5.23 / power

return foo * math.log((v_s30 ** power + bar ** power)
/ ((1360. ** power + bar ** power)))

def calc_depth_tor(self, mag, flag_rv):
'''Calculate an estimate of the depth to top of rupture (km).

Input:
mag: float
moment magnitude of the event
flag_rv: int
Reverse-faulting flag.
1: 30° <= rake angle <= 150° (combined reverse and
reverse-oblique)
0: otherwise

Returns:
depth_tor: float
Estimated depth to top of rupture (km)
'''

if flag_rv:
# Reverse and reverse-oblique faulting
foo = 2.704 - 1.226 * max(mag - 5.849, 0)
else:
# Combined strike-slip and normal faulting
foo = 2.673 - 1.136 * max(mag - 4.970, 0)

return max(foo, 0) ** 2

class AtkinsonBoore2006(model.Model):
'''Atkinson and Boore (2006) ground motion prediction model.

Developed for the Eastern North America with a reference velocity of 760
or 2000 m/s.'''

# Load the coefficients for the model
COEFF = dict(
)

PERIODS = COEFF['bc']['period']

INDEX_PGD = 0
INDEX_PGV = 1
INDEX_PGA = 2
INDICES_PSA = np.arange(3, 27)

def __init__(self, **kwds):
'''Compute the response predicted the Atkinson and Boore (2006) ground
motion model.

Input:
mag: float
moment magnitude of the event
dist_closest: float
closest distance to the rupture (km)
v_s30: float
time-averaged shear-wave velocity over the top 30 m of the site (m/s)
'''
super(AtkinsonBoore2006, self).__init__(
'Atkinson and Boore (2006)', 'AB06', **kwds)

mag = kwds['mag']
dist_closest = kwds['dist_closest']
v_s30 = kwds['v_s30']

def compute_log10_resp(period, c_1, c_2, c_3, c_4, c_5, c_6, c_7, c_8, c_9, c_10):
R0 = 10.0
R1 = 70.0
R2 = 140.0

f0 = max(np.log10(R0 / dist_closest), 0)
f1 = min(np.log10(dist_closest), np.log10(R1))
f2 = max(np.log10(dist_closest / R2), 0)

log10_resp = (c_1 + c_2 * mag + c_3 * mag ** 2 + (c_4 + c_5 * mag)
* f1 + (c_6 + c_7 * mag) * f2 + (c_8 + c_9 * mag) * f0 +
c_10 * dist_closest)

return log10_resp

def compute_log10_site(pga_bc, period, b_lin, b_1, b_2):
VS_1 = 180.
VS_2 = 300.
VS_REF = 760.

if v_s30 <= VS_1:
b_nl = b_1
elif VS_1 < v_s30 <= VS_2:
b_nl = ((b_1 - b_2) * np.log(v_s30 / VS_2) / np.log(VS_1 / VS_2))
elif VS_2 < v_s30 <= VS_REF:
b_nl = (b_2 * np.log(v_s30 / VS_REF) / np.log(VS_2 / VS_REF))
else:
# Vs30 > VS_REF
b_nl = 0

pga_bc = max(pga_bc, 60.)

log10_site = np.log10(np.exp(b_lin * np.log(v_s30 / VS_REF) + b_nl *
np.log(pga_bc / 100.)))

return log10_site

def compute_log10_stress_factor(stress_drop, period, delta, m1, mh):
foo = delta + 0.05
bar = 0.05 + delta * max(mag - m1, 0) / (mh - m1)

log10_stress_factor = min(2., stress_drop / 140.) * min(foo, bar)

return log10_stress_factor

COEFF = self.COEFF['bc'] if v_s30 else self.COEFF['rock']

# Compute the log10 PSA in units of cm/sec/sec
log10_resp = np.array([compute_log10_resp(*c) for c in COEFF])

# Apply the stress correction factor as recommended by Atkinson and
# Boore (2011)
if mag >= 5:
stress_drop = 10. ** (3.45 - 0.2 * mag)

log10_stress_factor = [compute_log10_stress_factor(
stress_drop, *c) for c in self.COEFF_SF]

log10_resp += np.interp(COEFF.period,
self.COEFF_SF.period, log10_stress_factor)

if v_s30:
# Compute the site amplification
pga_bc = (10 ** log10_resp[self.INDEX_PGA])

log10_site = [compute_log10_site(pga_bc, *c)
for c in self.COEFF_SITE]

# Amplification is specified at periods the differ from the ground
# motion model so we have to interpolate to a common period
# spacing before adding the influence of the site
log10_site = np.interp(COEFF.period,
self.COEFF_SITE.period, log10_site)
log10_resp += log10_site

# Convert from cm/sec/sec to gravity
log10_resp += np.log10(model.TO_GRAVITY)

# Convert to log-base 10
self._ln_resp = np.log(10 ** log10_resp)
self._ln_std = 0.30


These models tend to be quite complicated.

Maybe that's standard in your field, but they definitely are complicated, and converting those equations into code will never be pretty (I'm looking at you, calc_ln_resp_ref!).

However, this results in the functions being indent 3 levels. Is this bad?

No, it's not bad per se. You did an excellent job to keep the nesting under control with those helper functions, and the code is actually easy to read and follow.

Is there a better way to organize this code? Any other comments on the style or organization would be most welcome.

The code actually looks good. A few suggestions:

1. Feel free to organize the different models into different files if your current file is going to continue growing.
2. Add a docstring to the module (your existing docstrings in classes and methods are great).
3. Get more confident and keep up the good job! ^_^
• Thanks for the feedback. I have one model per file because there will be probably 15+ models with more being added as they are published. – arkottke Nov 29 '13 at 13:52