Merge pull request #112 from hvasbath/seis_derivative2

v1.2.2 bug fixes, GF backend versioning, traceplot improvements

CHANGELOG.md

@@ -4,7 +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.1] TBD 2022
+
+## [1.2.2] TBD 2022
+Contributors: Hannes Vasyura-Bathke @hvasbath
+
+### Added
+- plotting.marginals.traceplot: source_idxs argument can be slice e.g. 1:10 to take mean of patch attributes
+- heart.seis_derivative: function to numerically calculate derivatives wrt source attributes for waveforms
+- utility: functions for slice to string conversion and vice versa
+- config; NonlinearGFConfig added version attribute to specify backend version to use for GF calculation
+
+### Changed
+- plotting.marginals.traceplot: CDF plotting of multiple distributions in a single subplot marks quantiles only once
+
+### Fixed
+- plotting.marginals.traceplot: multipage output repeating variables
+
+
+## [1.2.1] 14 September 2022
 Contributors: Hannes Vasyura-Bathke @hvasbath
 
 ### Added

beat/apps/beat.py

@@ -35,7 +35,7 @@
 from beat.sampler.pt import SamplingHistory
 from beat.sampler.smc import sample_factor_final_stage
 from beat.sources import MTQTSource, MTSourceWithMagnitude
-from beat.utility import list2string
+from beat.utility import list2string, string2slice
 
 logger = logging.getLogger("beat")
 
@@ -1862,7 +1862,14 @@ def setup(parser):
     if len(options.source_idxs) < 1:
         source_idxs = None
     else:
-        source_idxs = [int(idx) for idx in options.source_idxs]
+        source_idxs = []
+        for str_idx in options.source_idxs:
+            try:
+                idx = int(str_idx)
+            except ValueError:
+                idx = string2slice(str_idx)
+
+            source_idxs.append(idx)
 
     po = plotting.PlotOptions(
         plot_projection=options.plot_projection,

beat/config.py

@@ -380,6 +380,10 @@ class NonlinearGFConfig(GFConfig):
         "depth the velocity model is not varied based on the errors "
         "defined above!",
     )
+    version = String.T(
+        default="",
+        help="Version number of the backend codes. If not defined, default versions will be used.",
+    )
 
 
 class SeismicGFConfig(NonlinearGFConfig):

beat/heart.py

@@ -2160,7 +2160,12 @@ def get_fomosto_baseconfig(gfconfig, event, station, waveforms, crust_ind):
 
 
 def choose_backend(
-    fomosto_config, code, source_model, receiver_model, gf_directory="qseis2d_green"
+    fomosto_config,
+    code,
+    source_model,
+    receiver_model,
+    gf_directory="qseis2d_green",
+    version=None,
 ):
     """
     Get backend related config.
@@ -2178,7 +2183,8 @@ def choose_backend(
         raise TypeError('For near-field phases the "qseis" backend has to be used!')
 
     if code == "qseis":
-        version = "2006a"
+        default_version = "2006a"
+        code_version = version or default_version
         if "slowest" in waveforms or distances.min() < 10:
             logger.info(
                 "Receiver and source"
@@ -2204,16 +2210,17 @@ def choose_backend(
             wavelet_duration_samples=0.001,
             sw_flat_earth_transform=sw_flat_earth_transform,
             sw_algorithm=sw_algorithm,
-            qseis_version=version,
+            qseis_version=code_version,
         )
 
     elif code == "qssp":
         source_model = copy.deepcopy(receiver_model)
         receiver_model = None
-        version = "2010"
+        default_version = "2010"
+        code_version = version or default_version
 
         conf = qssp.QSSPConfig(
-            qssp_version=version,
+            qssp_version=code_version,
             slowness_max=float(num.max(slowness_taper)),
             toroidal_modes=True,
             spheroidal_modes=True,
@@ -2223,6 +2230,8 @@ def choose_backend(
     else:
         raise NotImplementedError("Backend not supported: %s" % code)
 
+    logger.info("Using modelling code: %s with version: %s", code, code_version)
+
     # fill remaining fomosto params
     fc.earthmodel_1d = source_model
     fc.earthmodel_receiver_1d = receiver_model
@@ -2322,7 +2331,12 @@ def seis_construct_gf(
             gf_directory = os.path.join(sf.store_superdir, "base_gfs_%i" % crust_ind)
 
             conf = choose_backend(
-                fomosto_config, sf.code, source_model, receiver_model, gf_directory
+                fomosto_config,
+                sf.code,
+                source_model,
+                receiver_model,
+                gf_directory,
+                version=sf.version,
             )
 
             fomosto_config.validate()
@@ -2463,7 +2477,8 @@ def geo_construct_gf(event, geodetic_config, crust_ind=0, execute=True, force=Fa
     """
     from pyrocko.fomosto import psgrn_pscmp as ppp
 
-    version = "2008a"
+    default_version = "2008a"
+
     gfc = geodetic_config.gf_config
 
     # extract source crustal profile and check for water layer
@@ -2477,9 +2492,9 @@ def geo_construct_gf(event, geodetic_config, crust_ind=0, execute=True, force=Fa
 
     c = ppp.PsGrnPsCmpConfig()
 
-    c.pscmp_config.version = version
+    c.pscmp_config.version = gfc.version or default_version
 
-    c.psgrn_config.version = version
+    c.psgrn_config.version = gfc.version or default_version
     c.psgrn_config.sampling_interval = gfc.sampling_interval
     c.psgrn_config.gf_depth_spacing = gfc.medium_depth_spacing
     c.psgrn_config.gf_distance_spacing = gfc.medium_distance_spacing
@@ -3797,6 +3812,131 @@ def seis_synthetics(
         raise TypeError("Outmode %s not supported!" % outmode)
 
 
+spatial_derivative_parameters = {"dn": "north_shift", "de": "east_shift", "dd": "depth"}
+
+
+def seis_derivative(
+    engine,
+    sources,
+    targets,
+    arrival_taper,
+    arrival_times,
+    wavename,
+    filterer,
+    h,
+    parameter,
+    stencil_order=3,
+    outmode="tapered_data",
+):
+    """
+    Calculate numerical derivative with respect to source or spatial parameter
+    Parameters
+    ----------
+    engine : :class:`pyrocko.gf.seismosizer.LocalEngine`
+    sources : list
+        containing :class:`pyrocko.gf.seismosizer.Source` Objects
+        reference source is the first in the list!!!
+    targets : list
+        containing :class:`pyrocko.gf.seismosizer.Target` Objects
+    arrival_taper : :class:`ArrivalTaper`
+    arrival_times : list or:class:`numpy.ndarray`
+        containing the start times [s] since 1st.January 1970 to start
+        tapering
+    wavename : string
+        of the tabulated phase that determines the phase arrival
+    filterer : :class:`Filterer`
+    h : float
+        distance for derivative calculation
+    parameter : str
+        parameter with respect to which the derivative
+        is being calculated e.g. 'strike', 'dip', 'depth'
+    stencil_order : int
+        order N of numerical stencil differentiation, available; 3 or 5
+    Returns
+    -------
+    :class:`num.array` ntargets x nsamples with the first derivative
+    """
+
+    ntargets = len(targets)
+
+    available_params = list(sources[0].keys()) + list(
+        spatial_derivative_parameters.keys()
+    )
+    if parameter not in available_params:
+        raise AttributeError(
+            "Parameter for which the derivative was requested is neither"
+            " represented by the source nor the target. Supported parameters:"
+            " %s" % utility.list2string(available_params)
+        )
+
+    calc_sources = copy.deepcopy(sources)
+    store = engine.get_store(targets[0].store_id)
+    nsamples = int(num.ceil(store.config.sample_rate * arrival_taper.duration()))
+
+    stencil = utility.StencilOperator(h=h, order=stencil_order)
+
+    # loop over stencil steps
+    n_stencil_steps = len(stencil)
+    tmp = num.zeros((n_stencil_steps, ntargets, nsamples), dtype="float64")
+    for i, hstep in enumerate(stencil.hsteps):
+
+        if parameter in spatial_derivative_parameters:
+            target_param_name = spatial_derivative_parameters[parameter]
+            diff_targets = []
+            diff_sources = sources
+            for target in targets:
+                target_diff = copy.deepcopy(target)
+                target_param = getattr(target, target_param_name)
+                setattr(target_diff, target_param_name, target_param + hstep)
+                diff_targets.append(target_diff)
+
+            arrival_times = num.repeat(arrival_times, n_stencil_steps)
+        else:
+            diff_targets = targets
+            diff_sources = []
+            for source in calc_sources:
+                source_diff = source.clone()
+                source_param = source[parameter]
+                setattr(source_diff, parameter, source_param + hstep)
+                diff_sources.append(source_diff)
+
+        tmp[i, :, :], tmins = seis_synthetics(
+            engine=engine,
+            sources=diff_sources,
+            targets=diff_targets,
+            arrival_taper=arrival_taper,
+            wavename=wavename,
+            filterer=filterer,
+            arrival_times=arrival_times,
+            pre_stack_cut=True,
+            outmode="array",
+            chop_bounds=["b", "c"],
+        )
+
+    diff_array = (tmp * stencil.coefficients).sum(axis=0) / stencil.denominator
+
+    if outmode == "array":
+        return diff_array
+    elif outmode == "tapered_data":
+        out_traces = []
+        for i, target in enumerate(targets):
+            store = engine.get_store(target.store_id)
+            network, station, location, channel = target.codes
+            strain_trace = trace.Trace(
+                tmin=tmins[i],
+                ydata=diff_array[i, :],
+                network=network,
+                station=station,
+                location="{}{}".format(parameter, location),
+                channel=channel,
+                deltat=store.config.deltat,
+            )
+            out_traces.append(strain_trace)
+        return out_traces
+    else:
+        raise IOError("Outmode %s not supported!" % outmode)
+
+
 def radiation_weights_p(takeoff_angles, azimuths):
     """
     Station dependent propagation coefficients for P waves

beat/plotting/common.py

@@ -310,6 +310,8 @@ def histplot_op(
     """
 
     cumulative = kwargs.pop("cumulative", False)
+    nsources = kwargs.pop("nsources", False)
+    isource = kwargs.pop("isource", 0)
 
     if color is not None and cmap is not None:
         logger.debug("Using color for histogram edgecolor ...")
@@ -396,30 +398,39 @@ def histplot_op(
 
             for quantile, value in sigma_quants.items():
                 quantile /= 100.0
-                x = [leftb, value, value]
-                y = [quantile, quantile, 0.0]
+                if nsources == 1:
+                    x = [leftb, value, value]
+                    y = [quantile, quantile, 0.0]
+                else:
+                    x = [leftb, rightb]
+                    y = [quantile, quantile]
 
-                ax.plot(x, y, "--k", linewidth=0.5)
                 fontsize = 6
-                xval = (value - num.abs(leftb)) / 2 + leftb
-
-                ax.text(
-                    xval,
-                    quantile,
-                    "{}%".format(int(quantile * 100)),
-                    fontsize=fontsize,
-                    horizontalalignment="center",
-                    verticalalignment="bottom",
-                )
-                # if quantile > 0.3:
-                ax.text(
-                    value,
-                    quantile / 2,
-                    "%.3f" % value,
-                    fontsize=fontsize,
-                    horizontalalignment="left",
-                    verticalalignment="bottom",
-                )
+
+                if isource + 1 == nsources:
+                    # plot for last hist in axis
+                    ax.plot(x, y, "--k", linewidth=0.5)
+
+                    xval = (value - leftb) / 2 + leftb
+
+                    ax.text(
+                        xval,
+                        quantile,
+                        "{}%".format(int(quantile * 100)),
+                        fontsize=fontsize,
+                        horizontalalignment="center",
+                        verticalalignment="bottom",
+                    )
+
+                if nsources == 1:
+                    ax.text(
+                        value,
+                        quantile / 2,
+                        "%.3f" % value,
+                        fontsize=fontsize,
+                        horizontalalignment="left",
+                        verticalalignment="bottom",
+                    )
 
 
 def kde2plot_op(ax, x, y, grid=200, **kwargs):