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