Merge pull request #110 from hvasbath/ffi_cdf

v1.2.1 stage_posteriors improvements, multievent support

CHANGELOG.md

@@ -4,6 +4,21 @@ All notable changes to BEAT will be documented in this file.
 
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
 
+## [1.2.1] TBD 2022
+Contributors: Hannes Vasyura-Bathke @hvasbath
+
+### Added
+**FFI**
+- add calculation of coupling to derived variables at *Fault* for the summarize method
+
+**plotting**
+- *plot_projection* argument available for *stage_posteriors* plot: cdf, pdf, kde
+- update Example 8 to showcase the use of this functionality
+
+### Fixed
+- stage_posterior plot supports multipage output, figsize fixed to fractions of A4
+- multievent waveform_fits plot returns separate figures for each sub-event
+
 
 ## [1.2.0] 21 August 2022
 Contributors: Mahdi Hamidbeygi @mahdihamidbeygi, Hannes Vasyura-Bathke @hvasbath

beat/apps/beat.py

@@ -19,7 +19,7 @@
 from collections import OrderedDict
 from optparse import OptionParser
 
-from numpy import array, atleast_2d, floor, zeros
+from numpy import array, atleast_2d, floor, zeros, cumsum
 from pyrocko import model, util
 from pyrocko.gf import LocalEngine
 from pyrocko.guts import Dict, dump, load
@@ -1204,10 +1204,14 @@ def setup(parser):
                 buffer_size=sc.buffer_size,
                 progressbar=False,
             )
+
             pc = problem.config.problem_config
+            reference = pc.get_test_point()
 
             if options.calc_derived:
-                rtrace.add_derived_variables(pc.source_type, n_sources=pc.n_sources)
+                varnames, shapes = pc.get_derived_variables_shapes()
+                rtrace.add_derived_variables(varnames, shapes)
+                splitinds = cumsum([shape[0] for shape in shapes[:-1]])
 
             rtrace.setup(draws=draws, chain=-1, overwrite=True)
 
@@ -1225,6 +1229,7 @@ def setup(parser):
                 store = composite.engine.get_store(target.store_id)
             else:
                 source = composite.load_fault_geometry()
+                sources = [source]
                 engine = LocalEngine(
                     store_superdirs=[composite.config.gf_config.store_superdir]
                 )
@@ -1233,6 +1238,7 @@ def setup(parser):
             for chain in tqdm(chains):
                 for idx in idxs:
                     point = stage.mtrace.point(idx=idx, chain=chain)
+                    reference.update(point)
                     # normalize MT source, TODO put into get_derived_params
                     if isinstance(source, MTSourceWithMagnitude):
                         composite.point2sources(point)
@@ -1241,7 +1247,7 @@ def setup(parser):
                             ldicts.append(source.scaled_m6_dict)
 
                         jpoint = utility.join_points(ldicts)
-                        point.update(jpoint)
+                        reference.update(jpoint)
                         del jpoint, ldicts
 
                     derived = []
@@ -1250,12 +1256,13 @@ def setup(parser):
                         composite.point2sources(point)
                         if hasattr(source, "get_derived_parameters"):
                             for source in sources:
-                                derived.append(
-                                    source.get_derived_parameters(
-                                        store=store, target=target
-                                    )
+                                deri = source.get_derived_parameters(
+                                    point=reference,  # need to pass correction params
+                                    store=store,
+                                    target=target,
+                                    event=problem.config.event,
                                 )
-                                nderived = source.nderived_parameters
+                                derived.append(deri)
 
                         # pyrocko Rectangular source, TODO use BEAT RS ...
                         elif isinstance(source, RectangularSource):
@@ -1265,21 +1272,11 @@ def setup(parser):
                                     source.get_magnitude(store=store, target=target)
                                 )
 
-                            nderived = 1
-
-                        # FFI
-                        else:
-                            derived.append(
-                                source.get_magnitude(
-                                    point=point, store=store, target=target
-                                )
-                            )
-                            nderived = 1
-
                     lpoint = problem.model.lijection.d2l(point)
+
                     if derived:
                         lpoint.extend(
-                            map(ravel, split(vstack(derived).T, nderived, axis=0))
+                            map(ravel, split(vstack(derived).T, splitinds, axis=0))
                         )
 
                     # TODO: in PT with large buffer sizes somehow memory leak
@@ -1782,8 +1779,8 @@ def setup(parser):
             dest="plot_projection",
             # choices=['latlon', 'local', 'individual'],
             default="local",
-            help='Output projection of the plot; "latlon" or "local"'
-            'Default: "local"',
+            help='Output projection of the plot; "latlon" or "local" for maps - Default: "local";'
+            ' "pdf", "cdf" or "kde" for stage_posterior plot - Default: "pdf"',
         )
 
         parser.add_option(

beat/backend.py

@@ -44,7 +44,7 @@
 from pymc3.step_methods.arraystep import BlockedStep
 from pyrocko import util
 
-from beat.config import dc_components, mt_components, sample_p_outname, transd_vars_dist
+from beat.config import sample_p_outname, transd_vars_dist
 from beat.covariance import calc_sample_covariance
 from beat.utility import (
     ListArrayOrdering,
@@ -57,13 +57,6 @@
 logger = logging.getLogger("backend")
 
 
-derived_variables_mapping = {
-    "MTQTSource": mt_components + dc_components,
-    "MTSource": dc_components,
-    "RectangularSource": ["magnitude"],
-}
-
-
 def thin_buffer(buffer, buffer_thinning, ensure_last=True):
     """
     Reduce a list of objects by a given value.
@@ -297,19 +290,15 @@ def __len__(self):
         else:
             return self._df.shape[0] + len(self.buffer)
 
-    def add_derived_variables(self, source_type, n_sources=1):
-
-        try:
-            varnames = derived_variables_mapping[source_type]
-            logger.info(
-                "Adding derived variables %s to " "trace." % list2string(varnames)
+    def add_derived_variables(self, varnames, shapes):
+        nshapes = len(shapes)
+        nvars = len(varnames)
+        if nvars != nshapes:
+            raise ValueError(
+                "Inconsistent number of variables %i and shapes %i!" % (nvars, nshapes)
             )
-        except KeyError:
-            logger.info("No derived variables for %s" % source_type)
-            varnames = []
 
-        for varname in varnames:
-            shape = (n_sources,)
+        for varname, shape in zip(varnames, shapes):
             self.flat_names[varname] = ttab.create_flat_names(varname, shape)
             self.var_shapes[varname] = shape
             self.var_dtypes[varname] = "float64"

beat/config.py

@@ -42,6 +42,8 @@
     _domain_choices,
 )
 from beat.sources import MTQTSource, MTSourceWithMagnitude, RectangularSource
+from beat.utility import check_point_keys, list2string
+
 
 guts_prefix = "beat"
 
@@ -100,6 +102,7 @@
 ] + block_vars
 
 static_dist_vars = ["uparr", "uperp", "utens"]
+derived_dist_vars = ["coupling"]
 
 hypo_vars = ["nucleation_strike", "nucleation_dip", "time"]
 partial_kinematic_vars = ["durations", "velocities"]
@@ -110,7 +113,7 @@
 
 kinematic_dist_vars = static_dist_vars + partial_kinematic_vars + hypo_vars
 transd_vars_dist = partial_kinematic_vars + static_dist_vars + voronoi_locations
-dist_vars = static_dist_vars + partial_kinematic_vars
+dist_vars = static_dist_vars + partial_kinematic_vars + derived_dist_vars
 
 interseismic_catalog = {"geodetic": interseismic_vars}
 
@@ -130,6 +133,14 @@
     ]
 )
 
+derived_variables_mapping = {
+    "MTQTSource": mt_components + dc_components,
+    "MTSource": dc_components,
+    "RectangularSource": ["magnitude"],
+    "RectangularSourcePole": ["magnitude", "coupling"],
+}
+
+
 hyper_name_laplacian = "h_laplacian"
 
 moffdiag = (-1.0, 1.0)
@@ -1572,6 +1583,41 @@ def get_parameter_shape(self, param):
 
         return shape
 
+    def get_derived_variables_shapes(self):
+
+        source_type = self.source_type
+
+        tpoint = self.get_test_point()
+        has_pole, _ = check_point_keys(tpoint, phrase="*_pole_lat")
+
+        if has_pole:
+            source_type += "Pole"
+
+        try:
+            varnames = derived_variables_mapping[source_type]
+            shapes = []
+            for varname in varnames:
+                if self.mode == geometry_mode_str:
+                    shape = (self.n_sources,)
+                elif self.mode == ffi_mode_str:
+                    if varname == "magnitude":
+                        shape = (1,)
+                    else:
+                        shape = (self.mode_config.npatches,)
+
+                shapes.append(shape)
+
+            logger.info(
+                "Adding derived variables %s with shapes %s to "
+                "trace." % (list2string(varnames), list2string(shapes))
+            )
+        except KeyError:
+            logger.info("No derived variables for %s" % source_type)
+            varnames = []
+            shapes = []
+
+        return varnames, shapes
+
 
 class SamplerParameters(Object):
 

beat/ffi/fault.py

@@ -134,13 +134,16 @@ def set_model_resolution(self, model_resolution):
         self._model_resolution = model_resolution
 
     def get_model_resolution(self):
-        if self._model_resolution is None:
-            logger.warning(
-                "Model resolution matrix has not been calculated! "
-                "Please run beat.ffi.optimize_discretization on this fault! "
-                "Returning None!"
-            )
-        return self._model_resolution
+        if hasattr(self, "_model_resolution"):
+            if self._model_resolution is None:
+                logger.warning(
+                    "Model resolution matrix has not been calculated! "
+                    "Please run beat.ffi.optimize_discretization on this fault! "
+                    "Returning None!"
+                )
+            return self._model_resolution
+        else:
+            return None
 
     def get_subfault_key(self, index, datatype, component):
 
@@ -937,6 +940,18 @@ def needs_optimization(self):
     def is_discretized(self):
         return True if self.npatches else False
 
+    def get_derived_parameters(self, point=None, store=None, target=None, event=None):
+
+        has_pole, _ = check_point_keys(point, phrase="*_pole_lat")
+        if has_pole:
+            euler_slips = euler_pole2slips(point=point, fault=self, event=event)
+            coupling = backslip2coupling(point, euler_slips)
+        else:
+            coupling = []
+
+        magnitude = self.get_magnitude(point=point, store=store, target=target)
+        return num.hstack([magnitude, coupling])
+
 
 def write_fault_to_pscmp(
     filename, fault, point=None, force=False, export_patch_idxs=False

beat/models/geodetic.py

@@ -838,7 +838,7 @@ def point2sources(self, point):
             with random variables from solution space
         """
         tpoint = copy.deepcopy(point)
-        tpoint.update(self.fixed_rvs)
+        # tpoint.update(self.fixed_rvs)   if vars are fixed GFS are not calculated
 
         if self.nevents == 1:
             events = [self.event]  # single event