Merge pull request #121 from hvasbath/extend_solution_export

dev branch pre-release

CHANGELOG.md

@@ -4,6 +4,24 @@ 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.5]  24 Mai 2023
+Contributors: Hannes Vasyura-Bathke @hvasbath
+
+### Added
+- heart.ResultPoint: extend by event object from MAP
+- plotting.seismic,geodetic: add standardized residuals to fit plots
+- plotting.marginals: add transforms to MTQT parameters plotting marginals
+- plotting.marginals: multi-source support for correlation_hist
+- plotting.seismic: "station_variance_reductions" plot
+- added new source_type: "SFSource"
+
+### Changed
+- plotting.seismic.seismic_fits: subplots are ordered row wise distance based, channels in cols
+
+### Fixed
+- sampler.base.iter_parallel_chains: fix chunksize determination
+- minor logging messages
+
 ## [1.2.4] 14 February 2023
 Contributors: Hannes Vasyura-Bathke @hvasbath
 

beat/apps/beat.py

@@ -1918,7 +1918,7 @@ def setup(parser):
             plotting.plots_catalog[plot](problem, po)
         # except Exception as err:
         #    pass
-        except (TypeError, plotting.ModeError) as err:
+        except (TypeError, plotting.ModeError, NotImplementedError) as err:
 
             logger.warning(
                 "Could not plot %s got Error: %s \n %s"
@@ -2191,7 +2191,7 @@ def setup(parser):
 
         src_class_name = problem.config.problem_config.source_type
         for source in sources:
-            source.regularize()
+            # source.regularize()
             try:
                 sandbox.addSource(
                     talpa_source_catalog[src_class_name].fromPyrockoSource(
@@ -2316,8 +2316,12 @@ def setup(parser):
 
     var_reds = problem.get_variance_reductions(point)
 
+    events = problem.get_pyrocko_events(point)
     rpoint = heart.ResultPoint(
-        point=point, post_llk=options.post_llk, variance_reductions=var_reds
+        events=events,
+        point=point,
+        post_llk=options.post_llk,
+        variance_reductions=var_reds,
     )
     rpoint.regularize()
     rpoint.validate()

beat/config.py

@@ -59,6 +59,7 @@
 source_names = """
     ExplosionSource
     RectangularExplosionSource
+    SFSource
     DCSource
     CLVDSource
     MTSource
@@ -71,6 +72,7 @@
 source_classes = [
     gf.ExplosionSource,
     gf.RectangularExplosionSource,
+    gf.SFSource,
     gf.DCSource,
     gf.CLVDSource,
     MTSourceWithMagnitude,
@@ -110,6 +112,7 @@
 
 mt_components = ["mnn", "mee", "mdd", "mne", "mnd", "med"]
 dc_components = ["strike1", "dip1", "rake1", "strike2", "dip2", "rake2"]
+sf_components = ["fn", "fe", "fd"]
 
 kinematic_dist_vars = static_dist_vars + partial_kinematic_vars + hypo_vars
 transd_vars_dist = partial_kinematic_vars + static_dist_vars + voronoi_locations
@@ -143,6 +146,7 @@
 
 hyper_name_laplacian = "h_laplacian"
 
+sf_force = (0, 1e10)
 moffdiag = (-1.0, 1.0)
 mdiag = (-num.sqrt(2), num.sqrt(2))
 
@@ -172,6 +176,9 @@
     mne=moffdiag,
     mnd=moffdiag,
     med=moffdiag,
+    fn=sf_force,
+    fe=sf_force,
+    fd=sf_force,
     exx=(-200.0, 200.0),
     eyy=(-200.0, 200.0),
     exy=(-200.0, 200.0),
@@ -204,7 +211,7 @@
     bl_azimuth=(0, 180),
     bl_amplitude=(0.0, 0.1),
     locking_depth=(1.0, 10.0),
-    hypers=(-5.0, 8.0),
+    hypers=(-2.0, 6.0),
     ramp=(-0.005, 0.005),
     offset=(-0.05, 0.05),
     lat=(30.0, 30.5),

beat/heart.py

@@ -28,7 +28,7 @@
     Tuple,
 )
 from pyrocko.guts_array import Array
-from pyrocko.model import gnss
+from pyrocko.model import gnss, Event, get_effective_latlon
 from pyrocko.moment_tensor import to6
 from pyrocko.spit import OutOfBounds
 from scipy import linalg
@@ -436,6 +436,7 @@ class ResultPoint(Object):
     Containing point in solution space.
     """
 
+    events = List.T(Event.T(default=Event.D()), optional=True)
     post_llk = String.T(
         optional=True,
         help="describes which posterior likelihood value the point belongs to",
@@ -576,6 +577,9 @@ def results_for_export(results, datatype=None, attributes=None):
     diameter=(0.0, 100.0),
     sign=(-1.0, 1.0),
     volume_change=(-1e12, 1e12),
+    fn=(-1e20, 1e20),
+    fe=(-1e20, 1e20),
+    fd=(-1e20, 1e20),
     mnn=(-sqrt2, sqrt2),
     mee=(-sqrt2, sqrt2),
     mdd=(-sqrt2, sqrt2),

beat/models/base.py

@@ -26,7 +26,6 @@
 
 def get_hypervalue_from_point(point, observe, counter, hp_specific=False):
     hp_name = get_hyper_name(observe)
-
     if hp_name in point:
         if hp_specific:
             hp = point[hp_name][counter(hp_name)]

beat/models/geodetic.py

@@ -264,7 +264,7 @@ def save_covs(datasets, cov_mat="pred_v"):
             }
 
             outname = os.path.join(results_path, "%s_C_%s" % ("geodetic", cov_mat))
-            logger.info('"%s" to: %s' % (outname))
+            logger.info('"geodetic covariances" to: %s', outname)
             num.savez(outname, **covs)
 
         gc = self.config
@@ -627,6 +627,27 @@ def point2sources(self, point):
             # reset source time may result in store error otherwise
             source.time = 0.0
 
+    def get_pyrocko_events(self, point=None):
+        """
+        Transform sources to pyrocko events.
+
+        Returns
+        -------
+        events : list
+            of :class:`pyrocko.model.Event`
+        """
+
+        if point is not None:
+            self.point2sources(point)
+
+        events = []
+        target = self.targets[0]
+        store = self.engine.get_store(target.store_id)
+        for source in self.sources:
+            events.append(source.pyrocko_event(store=store, target=target))
+
+        return events
+
     def get_formula(self, input_rvs, fixed_rvs, hyperparams, problem_config):
         """
         Get geodetic likelihood formula for the model built. Has to be called

beat/models/polarity.py

@@ -179,6 +179,27 @@ def point2sources(self, point):
         for i, source in enumerate(self.sources):
             update_source(source, **source_points[i])
 
+    def get_pyrocko_events(self, point=None):
+        """
+        Transform sources to pyrocko events.
+
+        Returns
+        -------
+        events : list
+            of :class:`pyrocko.model.Event`
+        """
+
+        if point is not None:
+            self.point2sources(point)
+
+        events = []
+        target = self.targets[0]
+        store = self.engine.get_store(target.store_id)
+        for source in self.sources:
+            events.append(source.pyrocko_event(store=store, target=target))
+
+        return events
+
     def get_all_station_names(self):
         """
         Returns list of station names in the order of polarity maps.

beat/models/problems.py

@@ -7,6 +7,7 @@
 import theano.tensor as tt
 from pymc3 import Deterministic, Model, Potential, Uniform
 from pyrocko import util
+from pyrocko.model import get_effective_latlon
 from theano import config as tconfig
 
 from beat import config as bconfig
@@ -502,6 +503,35 @@ def point2sources(self, point):
         for composite in self.composites.values():
             self.composites[composite.name].point2sources(point)
 
+    def get_pyrocko_events(self, point=None):
+        """
+        Transform sources to pyrocko events.
+
+        Parameters
+        ----------
+        point : :func:`pymc3.Point`
+            Dictionary with model parameters, for which the events are
+            returned
+
+        Returns
+        -------
+        events : list
+            of :class:`pyrocko.model.Event`
+        """
+
+        for composite in self.composites.values():
+            if hasattr(composite, "get_pyrocko_events"):
+                events = composite.get_pyrocko_events(point)
+            else:
+                events = []
+
+            if events:
+                for event in events:
+                    elat, elon = get_effective_latlon(event)
+                    event.set_origin(elat, elon)
+
+                return events
+
     def update_weights(self, point, n_jobs=1, plot=False):
         """
         Calculate and update model prediction uncertainty covariances of

beat/models/seismic.py

@@ -714,6 +714,27 @@ def point2sources(self, point):
         for i, source in enumerate(self.sources):
             utility.update_source(source, **source_points[i])
 
+    def get_pyrocko_events(self, point=None):
+        """
+        Transform sources to pyrocko events.
+
+        Returns
+        -------
+        events : list
+            of :class:`pyrocko.model.Event`
+        """
+
+        if point is not None:
+            self.point2sources(point)
+
+        events = []
+        target = self.targets[0]
+        store = self.engine.get_store(target.store_id)
+        for source in self.sources:
+            events.append(source.pyrocko_event(store=store, target=target))
+
+        return events
+
     def get_formula(self, input_rvs, fixed_rvs, hyperparams, problem_config):
         """
         Get seismic likelihood formula for the model built. Has to be called

beat/plotting/__init__.py

@@ -20,6 +20,7 @@
     "fuzzy_mt_decomp": draw_fuzzy_mt_decomposition,
     "moment_rate": draw_moment_rate,
     "station_map": draw_station_map_gmt,
+    "station_variance_reductions": draw_station_variance_reductions,
 }
 
 
@@ -33,6 +34,7 @@
     "hudson",
     "lune",
     "fuzzy_beachball",
+    "station_variance_reductions",
 ]