Python base.RuptureContext类代码示例

 def test_get_mean_and_stddevs_good(self):
     """
     Tests the full execution of the GMPE tables for valid data
     """
     gsim = GMPETable(gmpe_table=self.TABLE_FILE)
     rctx = RuptureContext()
     rctx.mag = 6.0
     rctx.rake = 90.0
     dctx = DistancesContext()
     # Test values at the given distances and those outside range
     dctx.rjb = np.array([0.5, 1.0, 10.0, 100.0, 500.0])
     sctx = SitesContext()
     stddevs = [const.StdDev.TOTAL]
     expected_mean = np.array([20.0, 20.0, 10.0, 5.0, 1.0E-19])
     # PGA
     mean, sigma = gsim.get_mean_and_stddevs(sctx, rctx, dctx,
                                             imt_module.PGA(),
                                             stddevs)
     np.testing.assert_array_almost_equal(np.exp(mean), expected_mean, 5)
     np.testing.assert_array_almost_equal(sigma[0], 0.25 * np.ones(5), 5)
     # SA
     mean, sigma = gsim.get_mean_and_stddevs(sctx, rctx, dctx,
                                             imt_module.SA(1.0),
                                             stddevs)
     np.testing.assert_array_almost_equal(np.exp(mean), expected_mean, 5)
     np.testing.assert_array_almost_equal(sigma[0], 0.4 * np.ones(5), 5)

 def test_get_mean_and_stddevs(self):
     """
     Tests mean and standard deviations without amplification
     """
     gsim = GMPETable(gmpe_table=self.TABLE_FILE)
     rctx = RuptureContext()
     rctx.mag = 6.0
     dctx = DistancesContext()
     # Test values at the given distances and those outside range
     dctx.rjb = np.array([0.5, 1.0, 10.0, 100.0, 500.0])
     sctx = SitesContext()
     stddevs = [const.StdDev.TOTAL]
     expected_mean = np.array([2.0, 2.0, 1.0, 0.5, 1.0E-20])
     # PGA
     mean, sigma = gsim.get_mean_and_stddevs(sctx, rctx, dctx,
                                             imt_module.PGA(),
                                             stddevs)
     np.testing.assert_array_almost_equal(np.exp(mean), expected_mean, 5)
     np.testing.assert_array_almost_equal(sigma[0], 0.5 * np.ones(5), 5)
     # SA
     mean, sigma = gsim.get_mean_and_stddevs(sctx, rctx, dctx,
                                             imt_module.SA(1.0),
                                             stddevs)
     np.testing.assert_array_almost_equal(np.exp(mean), expected_mean, 5)
     np.testing.assert_array_almost_equal(sigma[0], 0.8 * np.ones(5), 5)
     # PGV
     mean, sigma = gsim.get_mean_and_stddevs(sctx, rctx, dctx,
                                             imt_module.PGV(),
                                             stddevs)
     np.testing.assert_array_almost_equal(np.exp(mean),
                                          10. * expected_mean,
                                          5)
     np.testing.assert_array_almost_equal(sigma[0], 0.5 * np.ones(5), 5)

 def test_get_mean_table(self, idx=0):
     """
     Test the retrieval of the mean amplification tables for a given
     magnitude and IMT
     """
     rctx = RuptureContext()
     rctx.mag = 6.0
     # PGA
     expected_table = np.ones([10, 2])
     expected_table[:, self.IDX] *= 1.5
     np.testing.assert_array_almost_equal(
         self.amp_table.get_mean_table(imt_module.PGA(), rctx),
         expected_table)
     # SA
     expected_table[:, self.IDX] = 2.0 * np.ones(10)
     np.testing.assert_array_almost_equal(
         self.amp_table.get_mean_table(imt_module.SA(0.5), rctx),
         expected_table)
     # SA (period interpolation)
     interpolator = interp1d(np.log10(self.amp_table.periods),
                             np.log10(np.array([1.5, 2.0, 0.5])))
     period = 0.3
     expected_table[:, self.IDX] = (
         10.0 ** interpolator(np.log10(period))) * np.ones(10.)
     np.testing.assert_array_almost_equal(
         self.amp_table.get_mean_table(imt_module.SA(period), rctx),
         expected_table)

 def _get_event_context(self, idx, nodal_plane_index=1):
     """
     Returns the event contexts for a specific event
     """
     idx = idx[0]
     rctx = RuptureContext()
     rup = self.records[idx]
     setattr(rctx, 'mag', rup.event.magnitude.value)
     if nodal_plane_index == 2:
         setattr(rctx, 'dip',
             rup.event.mechanism.nodal_planes.nodal_plane_2['dip'])
         setattr(rctx, 'rake', 
             rup.event.mechanism.nodal_planes.nodal_plane_2['rake'])
     else:
         setattr(rctx, 'dip',
             rup.event.mechanism.nodal_planes.nodal_plane_1['dip'])
         setattr(rctx, 'rake',
             rup.event.mechanism.nodal_planes.nodal_plane_1['rake'])
     if not rctx.rake:
         rctx.rake = rup.event.mechanism.get_rake_from_mechanism_type()
     if rup.event.rupture:
         setattr(rctx, 'ztor', rup.event.rupture.depth)
         setattr(rctx, 'width', rup.event.rupture.width)
     setattr(rctx, 'hypo_depth', rup.event.depth)
     return rctx

 def _get_event_context(self, idx, nodal_plane_index=1):
     """
     Returns the event contexts for a specific event
     """
     idx = idx[0]
     rctx = RuptureContext()
     rup = self.records[idx]
     setattr(rctx, 'mag', rup.event.magnitude.value)
     if nodal_plane_index == 2:
         setattr(rctx, 'strike',
             rup.event.mechanism.nodal_planes.nodal_plane_2['strike'])
         setattr(rctx, 'dip',
             rup.event.mechanism.nodal_planes.nodal_plane_2['dip'])
         setattr(rctx, 'rake',
             rup.event.mechanism.nodal_planes.nodal_plane_2['rake'])
     else:
         setattr(rctx, 'strike', 0.0)
         setattr(rctx, 'dip', 90.0)
         rctx.rake = rup.event.mechanism.get_rake_from_mechanism_type()
     if rup.event.rupture.surface:
         setattr(rctx, 'ztor', rup.event.rupture.surface.get_top_edge_depth())
         setattr(rctx, 'width', rup.event.rupture.surface.width)
         setattr(rctx, 'hypo_loc', rup.event.rupture.surface.get_hypo_location(1000))
     else:
         setattr(rctx, 'ztor', rup.event.depth)
         # Use the PeerMSR to define the area and assuming an aspect ratio
         # of 1 get the width
         setattr(rctx, 'width',
                 np.sqrt(DEFAULT_MSR.get_median_area(rctx.mag, 0)))
         # Default hypocentre location to the middle of the rupture
         setattr(rctx, 'hypo_loc', (0.5, 0.5))
     setattr(rctx, 'hypo_depth', rup.event.depth)
     setattr(rctx, 'hypo_lat', rup.event.latitude)
     setattr(rctx, 'hypo_lon', rup.event.longitude)
     return rctx

    def test_dist_not_in_increasing_order(self):
        sctx = SitesContext()
        rctx = RuptureContext()
        dctx = DistancesContext()

        rctx.mag = 5.
        dctx.rhypo = numpy.array([150, 100])
        mean_150_100, _ = self.GSIM_CLASS().get_mean_and_stddevs(
            sctx, rctx, dctx, SA(0.1, 5), [StdDev.TOTAL]
        )

        dctx.rhypo = numpy.array([100, 150])
        mean_100_150, _ = self.GSIM_CLASS().get_mean_and_stddevs(
            sctx, rctx, dctx, SA(0.1, 5), [StdDev.TOTAL]
        )
        self.assertAlmostEqual(mean_150_100[1], mean_100_150[0])
        self.assertAlmostEqual(mean_150_100[0], mean_100_150[1])

 def test_get_amplification_factors(self):
     """
     Tests the amplification tables
     """
     rctx = RuptureContext()
     rctx.mag = 6.0
     dctx = DistancesContext()
     # Takes distances at the values found in the table (not checking
     # distance interpolation)
     dctx.rjb = np.copy(self.amp_table.distances[:, 0, 0])
     # Test Vs30 is 700.0 m/s midpoint between the 400 m/s and 1000 m/s
     # specified in the table
     sctx = SitesContext()
     sctx.vs30 = 700.0 * np.ones_like(dctx.rjb)
     stddevs = [const.StdDev.TOTAL]
     expected_mean = np.ones_like(dctx.rjb)
     expected_sigma = np.ones_like(dctx.rjb)
     # Check PGA and PGV
     mean_amp, sigma_amp = self.amp_table.get_amplification_factors(
         imt_module.PGA(), sctx, rctx, dctx.rjb, stddevs)
     np.testing.assert_array_almost_equal(
         mean_amp,
         midpoint(1.0, 1.5) * expected_mean)
     np.testing.assert_array_almost_equal(
         sigma_amp[0],
         0.9 * expected_mean)
     mean_amp, sigma_amp = self.amp_table.get_amplification_factors(
         imt_module.PGV(), sctx, rctx, dctx.rjb, stddevs)
     np.testing.assert_array_almost_equal(
         mean_amp,
         midpoint(1.0, 0.5) * expected_mean)
     np.testing.assert_array_almost_equal(
         sigma_amp[0],
         0.9 * expected_mean)
     # Sa (0.5)
     mean_amp, sigma_amp = self.amp_table.get_amplification_factors(
         imt_module.SA(0.5), sctx, rctx, dctx.rjb, stddevs)
     np.testing.assert_array_almost_equal(
         mean_amp,
         midpoint(1.0, 2.0) * expected_mean)
     np.testing.assert_array_almost_equal(
         sigma_amp[0],
         0.9 * expected_mean)

 def test_get_mean_stddevs_unsupported_stddev(self):
     """
     Tests the execution of the GMPE with an unsupported standard deviation
     type
     """
     gsim = GMPETable(gmpe_table=self.TABLE_FILE)
     rctx = RuptureContext()
     rctx.mag = 6.0
     dctx = DistancesContext()
     # Test values at the given distances and those outside range
     dctx.rjb = np.array([0.5, 1.0, 10.0, 100.0, 500.0])
     sctx = SitesContext()
     sctx.vs30 = 1000. * np.ones(5)
     stddevs = [const.StdDev.TOTAL, const.StdDev.INTER_EVENT]
     with self.assertRaises(ValueError) as ve:
         gsim.get_mean_and_stddevs(sctx, rctx, dctx, imt_module.PGA(),
                                   stddevs)
     self.assertEqual(str(ve.exception),
                      "Standard Deviation type Inter event not supported")

    def check_gmpe_adjustments(self, adj_gmpe_set, original_gmpe):
        """
        Takes a set of three adjusted GMPEs representing the "low", "middle"
        and "high" stress drop adjustments for Germany and compares them
        against the original "target" GMPE for a variety of magnitudes
        and styles of fauling.
        """
        low_gsim, mid_gsim, high_gsim = adj_gmpe_set
        tot_std = [const.StdDev.TOTAL]
        for imt in self.imts:
            for mag in self.mags:
                for rake in self.rakes:
                    rctx = RuptureContext()
                    rctx.mag = mag
                    rctx.rake = rake
                    rctx.hypo_depth = 10.
                    # Get "original" values
                    mean = original_gmpe.get_mean_and_stddevs(self.sctx, rctx,
                                                              self.dctx, imt,
                                                              tot_std)[0]
                    mean = np.exp(mean)
                    # Get "low" adjustments (0.75 times the original)
                    low_mean = low_gsim.get_mean_and_stddevs(self.sctx, rctx,
                                                             self.dctx, imt,
                                                             tot_std)[0]
                    np.testing.assert_array_almost_equal(
                        np.exp(low_mean) / mean, 0.75 * np.ones_like(low_mean))

                    # Get "middle" adjustments (1.25 times the original)
                    mid_mean = mid_gsim.get_mean_and_stddevs(self.sctx, rctx,
                                                             self.dctx, imt,
                                                             tot_std)[0]
                    np.testing.assert_array_almost_equal(
                        np.exp(mid_mean) / mean, 1.25 * np.ones_like(mid_mean))

                    # Get "high" adjustments (1.5 times the original)
                    high_mean = high_gsim.get_mean_and_stddevs(self.sctx, rctx,
                                                               self.dctx, imt,
                                                               tot_std)[0]
                    np.testing.assert_array_almost_equal(
                        np.exp(high_mean) / mean,
                        1.5 * np.ones_like(high_mean))

 def test_get_sigma_table(self):
     """
     Test the retrieval of the standard deviation modification tables
     for a given magnitude and IMT
     """
     rctx = RuptureContext()
     rctx.mag = 6.0
     # PGA
     expected_table = np.ones([10, 2])
     expected_table[:, self.IDX] *= 0.8
     stddevs = ["Total"]
     pga_table = self.amp_table.get_sigma_tables(imt_module.PGA(),
                                                 rctx,
                                                 stddevs)[0]
     np.testing.assert_array_almost_equal(pga_table, expected_table)
     # SA (for coverage)
     sa_table = self.amp_table.get_sigma_tables(imt_module.SA(0.3),
                                                rctx,
                                                stddevs)[0]
     np.testing.assert_array_almost_equal(sa_table, expected_table)

 def test_rhypo_smaller_than_15(self):
     # test the calculation in case of rhypo distances less than 15 km
     # (for rhypo=0 the distance term has a singularity). In this case the
     # method should return values equal to the ones obtained by clipping
     # distances at 15 km.
     sctx = SitesContext()
     sctx.vs30 = numpy.array([800.0, 800.0, 800.0])
     rctx = RuptureContext()
     rctx.mag = 5.0
     rctx.rake = 0
     dctx = DistancesContext()
     dctx.rhypo = numpy.array([0.0, 10.0, 16.0])
     dctx.rhypo.flags.writeable = False
     mean_0, stds_0 = self.GSIM_CLASS().get_mean_and_stddevs(
         sctx, rctx, dctx, PGA(), [StdDev.TOTAL])
     setattr(dctx, 'rhypo', numpy.array([15.0, 15.0, 16.0]))
     mean_15, stds_15 = self.GSIM_CLASS().get_mean_and_stddevs(
         sctx, rctx, dctx, PGA(), [StdDev.TOTAL])
     numpy.testing.assert_array_equal(mean_0, mean_15)
     numpy.testing.assert_array_equal(stds_0, stds_15)

 def get_response_spectrum(self, magnitude, distance, periods, rake=90, vs30=800, damping=0.05):
     """
     """
     responses = np.zeros((len(periods),))
     p_damping = damping * 100
     rup = RuptureContext()
     rup.mag = magnitude
     rup.rake = rake
     dists = DistancesContext()
     dists.rjb = np.array([distance])
     sites = SitesContext()
     sites.vs30 = np.array([vs30])
     stddev_types = [StdDev.TOTAL]
     for i, period in enumerate(periods):
         if period == 0:
             imt = _PGA()
         else:
             imt = _SA(period, p_damping)
         responses[i] = np.exp(self._gmpe.get_mean_and_stddevs(sites, rup, dists, imt, stddev_types)[0][0])
     return ResponseSpectrum(periods, responses, unit='g', damping=damping)

    def test_equality(self):
        sctx1 = SitesContext()
        sctx1.vs30 = numpy.array([500., 600., 700.])
        sctx1.vs30measured = True
        sctx1.z1pt0 = numpy.array([40., 50., 60.])
        sctx1.z2pt5 = numpy.array([1, 2, 3])

        sctx2 = SitesContext()
        sctx2.vs30 = numpy.array([500., 600., 700.])
        sctx2.vs30measured = True
        sctx2.z1pt0 = numpy.array([40., 50., 60.])
        sctx2.z2pt5 = numpy.array([1, 2, 3])

        self.assertTrue(sctx1 == sctx2)

        sctx2 = SitesContext()
        sctx2.vs30 = numpy.array([500., 600.])
        sctx2.vs30measured = True
        sctx2.z1pt0 = numpy.array([40., 50., 60.])
        sctx2.z2pt5 = numpy.array([1, 2, 3])

        self.assertTrue(sctx1 != sctx2)

        sctx2 = SitesContext()
        sctx2.vs30 = numpy.array([500., 600., 700.])
        sctx2.vs30measured = False
        sctx2.z1pt0 = numpy.array([40., 50., 60.])
        sctx2.z2pt5 = numpy.array([1, 2, 3])

        self.assertTrue(sctx1 != sctx2)

        sctx2 = SitesContext()
        sctx2.vs30 = numpy.array([500., 600., 700.])
        sctx2.vs30measured = True
        sctx2.z1pt0 = numpy.array([40., 50., 60.])

        self.assertTrue(sctx1 != sctx2)

        rctx = RuptureContext()
        rctx.mag = 5.
        self.assertTrue(sctx1 != rctx)

    def test_mag_greater_8pt5(self):
        gmpe = SadighEtAl1997()

        sctx = SitesContext()
        rctx = RuptureContext()
        dctx = DistancesContext()

        rctx.rake =  0.0
        dctx.rrup = numpy.array([0., 1.])
        sctx.vs30 = numpy.array([800., 800.])

        rctx.mag = 9.0
        mean_rock_9, _ = gmpe.get_mean_and_stddevs(
            sctx, rctx, dctx, PGA(), [StdDev.TOTAL]
        )
        rctx.mag = 8.5
        mean_rock_8pt5, _ = gmpe.get_mean_and_stddevs(
            sctx, rctx, dctx, PGA(), [StdDev.TOTAL]
        )
        numpy.testing.assert_allclose(mean_rock_9, mean_rock_8pt5)

        sctx.vs30 = numpy.array([300., 300.])
        rctx.mag = 9.0
        mean_soil_9, _ = gmpe.get_mean_and_stddevs(
            sctx, rctx, dctx, PGA(), [StdDev.TOTAL]
        )
        rctx.mag = 8.5
        mean_soil_8pt5, _ = gmpe.get_mean_and_stddevs(
            sctx, rctx, dctx, PGA(), [StdDev.TOTAL]
        )
        numpy.testing.assert_allclose(mean_soil_9, mean_soil_8pt5)

    def test_mag_dist_outside_range(self):
        sctx = SitesContext()
        rctx = RuptureContext()
        dctx = DistancesContext()

        # rupture with Mw = 3 (Mblg=2.9434938048208452) at rhypo = 1 must give
        # same mean as rupture with Mw = 4.4 (Mblg=4.8927897867183798) at
        # rhypo = 10
        rctx.mag = 2.9434938048208452
        dctx.rhypo = numpy.array([1])
        mean_mw3_d1, _ = self.GSIM_CLASS().get_mean_and_stddevs(
            sctx, rctx, dctx, SA(0.1, 5), [StdDev.TOTAL]
        )

        rctx.mag = 4.8927897867183798
        dctx.rhypo = numpy.array([10])
        mean_mw4pt4_d10, _ = self.GSIM_CLASS().get_mean_and_stddevs(
            sctx, rctx, dctx, SA(0.1, 5), [StdDev.TOTAL]
        )

        self.assertAlmostEqual(float(mean_mw3_d1), float(mean_mw4pt4_d10))

        # rupture with Mw = 9 (Mblg = 8.2093636421088814) at rhypo = 1500 km
        # must give same mean as rupture with Mw = 8.2
        # (Mblg = 7.752253535347597) at rhypo = 1000
        rctx.mag = 8.2093636421088814
        dctx.rhypo = numpy.array([1500.])
        mean_mw9_d1500, _ = self.GSIM_CLASS().get_mean_and_stddevs(
            sctx, rctx, dctx, SA(0.1, 5), [StdDev.TOTAL]
        )

        rctx.mag = 7.752253535347597
        dctx.rhypo = numpy.array([1000.])
        mean_mw8pt2_d1000, _ = self.GSIM_CLASS().get_mean_and_stddevs(
            sctx, rctx, dctx, SA(0.1, 5), [StdDev.TOTAL]
        )

        self.assertAlmostEqual(mean_mw9_d1500, mean_mw8pt2_d1000)