5
\$\begingroup\$

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))

def load_data_file(name, skip_header=None):
    fname = os.path.join(os.path.dirname(__file__), 'data', name)
    return np.recfromcsv(fname, skip_header=skip_header).view(np.recarray)

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

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

    # Load the coefficients for the model
    COEFF = model.load_data_file('chiou_youngs_2013.csv', 2)

    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(
        bc=model.load_data_file('atkinson_boore_2006-bc.csv', 2),
        rock=model.load_data_file('atkinson_boore_2006-rock.csv', 2),
    )

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

    COEFF_SITE = model.load_data_file('atkinson_boore_2006-site.csv', 2)
    COEFF_SF = model.load_data_file('atkinson_boore_2006-sf.csv', 2)

    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
\$\endgroup\$
4
\$\begingroup\$

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! ^_^
\$\endgroup\$
  • \$\begingroup\$ 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. \$\endgroup\$ – arkottke Nov 29 '13 at 13:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.