Gnss seis joint

.github/workflows/build-wheels.yaml

@@ -8,7 +8,7 @@ jobs:
     runs-on: ${{ matrix.os }}
     strategy:
       matrix:
-        os: [ubuntu-20.04, ubuntu-22.04, ubuntu-24.04, macOS-13, macOS-14]
+        os: [ubuntu-22.04, ubuntu-24.04, macOS-13, macOS-14]
         # Exclude windows-2019
 
     steps:

beat/apps/beat.py

@@ -1061,7 +1061,6 @@ def setup(parser):
     problem = load_model(project_dir, options.mode, options.hypers)
 
     step = problem.init_sampler(hypers=options.hypers)
-
     if options.hypers:
         estimate_hypers(step, problem)
     else:

beat/covariance.py

@@ -465,7 +465,7 @@ def model_prediction_sensitivity(engine, *args, **kwargs):
         kwargs["source_params"] = args[2]
 
     request = kwargs.pop("request", None)
-    nprocs = kwargs.pop("nprocs", 1)
+    nthreads = kwargs.pop("nthreads", 1)
     source_params = kwargs.pop("source_params", None)
     h = kwargs.pop("h", None)
 
@@ -505,7 +505,7 @@ def model_prediction_sensitivity(engine, *args, **kwargs):
             ]
 
             response = engine.process(
-                sources=calc_sources, targets=request.targets, nprocs=nprocs
+                sources=calc_sources, targets=request.targets, nthreads=nthreads
             )
 
             for i_k in range(len(request.targets)):

beat/defaults.py

@@ -242,7 +242,7 @@ def __getitem__(self, k):
         physical_bounds=(-10.0, 10.0), default_bounds=(-2.0, 6.0), unit=u_hyp
     ),
     "ramp": Bounds(
-        physical_bounds=(-0.005, 0.005), default_bounds=(-0.005, 0.005), unit=u_rad
+        physical_bounds=(-0.1, 0.1), default_bounds=(-0.005, 0.005), unit=u_rad
     ),
     "offset": Bounds(
         physical_bounds=(-0.05, 0.05), default_bounds=(-0.05, 0.05), unit=u_m

beat/ffi/base.py

@@ -451,11 +451,12 @@ def put(self, entries, targetidx, patchidx, durations, starttimes):
         if len(entries.shape) < 2:
             raise ValueError("Entries have to be 2d arrays!")
 
-        if entries.shape[1] != self.nsamples:
+        synthetics_nsamples = entries.shape[1]
+        if synthetics_nsamples != self.nsamples:
             raise GFLibraryError(
                 "Trace length of entries is not consistent with the library"
                 " to be filled! Entries length: %i Library: %i."
-                % (entries.shape[0], self.nsamples)
+                % (synthetics_nsamples, self.nsamples)
             )
 
         self._check_setup()

beat/heart.py

@@ -23,6 +23,7 @@
     Object,
     String,
     StringChoice,
+    StringUnion,
     Tuple,
 )
 from pyrocko.guts_array import Array
@@ -425,6 +426,60 @@ def apply(self, trace):
         trace.set_ydata(new_trace.ydata)
 
 
+class DynamicTarget(gf.Target):
+    response = trace.PoleZeroResponse.T(default=None, optional=True)
+    domain = StringChoice.T(
+        default="time",
+        choices=_domain_choices,
+        help="Domain for signal processing and likelihood calculation.",
+    )
+
+    @property
+    def nslcd_id(self):
+        return tuple(list(self.codes) + [self.domain])
+
+    @property
+    def nslcd_id_str(self):
+        return utility.list2string(self.nslcd_id)
+
+    def update_response(self, magnification, damping, period):
+        z, p, k = proto2zpk(magnification, damping, period, quantity="displacement")
+        # b, a = zpk2tf(z, p, k)
+
+        if self.response:
+            self.response.zeros = z
+            self.response.poles = p
+            self.response.constant = k
+        else:
+            logger.debug("Initializing new response!")
+            self.response = trace.PoleZeroResponse(zeros=z, poles=p, constant=k)
+
+    def update_target_times(self, sources=None, taperer=None):
+        """
+        Update the target attributes tmin and tmax to do the stacking
+        only in this interval. Adds twice taper fade in time to each taper
+        side.
+
+        Parameters
+        ----------
+        source : list
+            containing :class:`pyrocko.gf.seismosizer.Source` Objects
+        taperer : :class:`pyrocko.trace.CosTaper`
+        """
+
+        if sources is None or taperer is None:
+            self.tmin = None
+            self.tmax = None
+        else:
+            tolerance = 2 * (taperer.b - taperer.a)
+            self.tmin = taperer.a - tolerance
+            self.tmax = taperer.d + tolerance
+
+
+class TargetUnion(StringUnion):
+    members = [DynamicTarget.T(), String.T()]
+
+
 class ResultPoint(Object):
     """
     Containing point in solution space.
@@ -442,7 +497,7 @@ class ResultPoint(Object):
         help="Point in Solution space for which result is produced.",
     )
     variance_reductions = Dict.T(
-        String.T(),
+        TargetUnion.T(),
         Float.T(),
         default={},
         optional=True,
@@ -1009,56 +1064,6 @@ def get_xdata(self):
         return num.arange(len(self.ydata), dtype=num.float64) * self.deltaf + self.fmin
 
 
-class DynamicTarget(gf.Target):
-    response = trace.PoleZeroResponse.T(default=None, optional=True)
-    domain = StringChoice.T(
-        default="time",
-        choices=_domain_choices,
-        help="Domain for signal processing and likelihood calculation.",
-    )
-
-    @property
-    def nslcd_id(self):
-        return tuple(list(self.codes) + [self.domain])
-
-    @property
-    def nslcd_id_str(self):
-        return utility.list2string(self.nslcd_id)
-
-    def update_response(self, magnification, damping, period):
-        z, p, k = proto2zpk(magnification, damping, period, quantity="displacement")
-        # b, a = zpk2tf(z, p, k)
-
-        if self.response:
-            self.response.zeros = z
-            self.response.poles = p
-            self.response.constant = k
-        else:
-            logger.debug("Initializing new response!")
-            self.response = trace.PoleZeroResponse(zeros=z, poles=p, constant=k)
-
-    def update_target_times(self, sources=None, taperer=None):
-        """
-        Update the target attributes tmin and tmax to do the stacking
-        only in this interval. Adds twice taper fade in time to each taper
-        side.
-
-        Parameters
-        ----------
-        source : list
-            containing :class:`pyrocko.gf.seismosizer.Source` Objects
-        taperer : :class:`pyrocko.trace.CosTaper`
-        """
-
-        if sources is None or taperer is None:
-            self.tmin = None
-            self.tmax = None
-        else:
-            tolerance = 2 * (taperer.b - taperer.a)
-            self.tmin = taperer.a - tolerance
-            self.tmax = taperer.d + tolerance
-
-
 class GeodeticDataset(gf.meta.MultiLocation):
     """
     Overall geodetic data set class
@@ -3565,7 +3570,7 @@ def seis_synthetics(
     filterer=None,
     reference_taperer=None,
     plot=False,
-    nprocs=1,
+    nthreads=1,
     outmode="array",
     pre_stack_cut=False,
     taper_tolerance_factor=0.0,
@@ -3592,9 +3597,8 @@ def seis_synthetics(
     filterer : :class:`Filterer`
     plot : boolean
         flag for looking at traces
-    nprocs : int
-        number of processors to use for synthetics calculation
-        --> currently no effect !!!
+    nthreads : int
+        number of threads to use for synthetics calculation
     outmode : string
         output format of synthetics can be 'array', 'stacked_traces',
         'data' returns traces unstacked including post-processing,
@@ -3651,7 +3655,7 @@ def seis_synthetics(
 
     t_2 = time()
     try:
-        response = engine.process(sources=sources, targets=targets, nprocs=nprocs)
+        response = engine.process(sources=sources, targets=targets, nthreads=nthreads)
         t_1 = time()
     except IndexError:
         for source in sources:
@@ -4152,7 +4156,7 @@ def fft_transforms(
 
 
 def geo_synthetics(
-    engine, targets, sources, outmode="stacked_array", plot=False, nprocs=1
+    engine, targets, sources, outmode="stacked_array", plot=False, nthreads=1
 ):
     """
     Calculate synthetic displacements for a given static fomosto Greens
@@ -4168,9 +4172,8 @@ def geo_synthetics(
         containing :class:`pyrocko.gf.seismosizer.Target` Objects
     plot : boolean
         flag for looking at synthetics - not implemented yet
-    nprocs : int
-        number of processors to use for synthetics calculation
-        --> currently no effect !!!
+    nthreads : int
+        number of threads to use for synthetics calculation
     outmode : string
         output format of synthetics can be: 'array', 'arrays',
         'stacked_array','stacked_arrays'
@@ -4194,7 +4197,7 @@ def geo_synthetics(
         for target in targets:
             print(target)
 
-    response = engine.process(sources, targets)
+    response = engine.process(sources, targets, nthreads=nthreads)
 
     ns = len(sources)
     nt = len(targets)

beat/inputf.py

@@ -427,7 +427,7 @@ def rotate_traces_and_stations(datatraces, stations, event):
             traces = station2traces[station.station]
         except KeyError:
             logger.warning(
-                'Did not find data traces for station "%s"' % stations.station
+                'Did not find data traces for station "%s"' % station.station
             )
             continue
 

beat/models/geodetic.py

@@ -133,9 +133,12 @@ def n_t(self):
 
     def get_all_dataset_ids(self, hp_name):
         """
-        Return unique GNSS stations and radar acquisitions.
+        Return unique GNSS stations and radar acquisitions for a hyperparameter.
         """
-        return [dataset.id for dataset in self.datasets]
+        hp_dataset_typ = hp_name.split("_")[1]
+        return [
+            dataset.id for dataset in self.datasets if dataset.typ == hp_dataset_typ
+        ]
 
     def analyse_noise(self, tpoint=None):
         """
@@ -716,6 +719,7 @@ def get_synthetics(self, point):
             targets=self.targets,
             sources=self.sources,
             outmode="stacked_arrays",
+            nthreads=4,
         )
 
         synths = []

beat/models/seismic.py

@@ -558,7 +558,7 @@ def get_standardized_residuals(
             )
             ydata = result.processed_res.get_ydata()
             choli = num.linalg.inv(dataset.covariance.chol(num.exp(hp * 2.0)))
-            stdz_residuals[target.nslcd_id_str] = choli.dot(ydata)
+            stdz_residuals[target] = choli.dot(ydata)
         return stdz_residuals
 
     def get_variance_reductions(
@@ -603,7 +603,7 @@ def get_variance_reductions(
         hp_specific = self.config.dataset_specific_residual_noise_estimation
 
         var_reds = OrderedDict()
-        for result, tr in zip(results, self.datasets):
+        for result, tr, target in zip(results, self.datasets, self.targets):
             nslcd_id_str = result.processed_obs.nslcd_id_str
 
             hp = get_hypervalue_from_point(point, tr, counter, hp_specific=hp_specific)
@@ -617,11 +617,10 @@ def get_variance_reductions(
 
             logger.debug("nom %f, denom %f" % (float(nom), float(denom)))
 
-            var_reds[nslcd_id_str] = float(1 - (nom / denom))
+            var_reds[target] = float(1 - (nom / denom))
 
             logger.debug(
-                "Variance reduction for %s is %f"
-                % (nslcd_id_str, var_reds[nslcd_id_str])
+                "Variance reduction for %s is %f" % (nslcd_id_str, var_reds[target])
             )
 
             if 0:
@@ -810,7 +809,8 @@ def get_formula(self, input_rvs, fixed_rvs, hyperparams, problem_config):
                 engine=self.engine,
                 sources=sources,
                 targets=wmap.targets,
-                event=self.events[wc.event_idx],
+                events=self.events,
+                event_idx=wc.event_idx,
                 arrival_taper=wc.arrival_taper,
                 arrival_times=wmap._arrival_times,
                 wavename=wmap.name,
@@ -858,7 +858,7 @@ def get_synthetics(self, point, **kwargs):
         outmode = kwargs.pop("outmode", "stacked_traces")
         chop_bounds = kwargs.pop("chop_bounds", ["a", "d"])
         order = kwargs.pop("order", "list")
-        nprocs = kwargs.pop("nprocs", 4)
+        nthreads = kwargs.pop("nthreads", 4)
         force = kwargs.pop("force", False)
 
         self.point2sources(point)
@@ -868,6 +868,7 @@ def get_synthetics(self, point, **kwargs):
         obs = []
         for i, wmap in enumerate(self.wavemaps):
             wc = wmap.config
+            # print(wc)
             if not wmap.is_prepared or force:
                 wmap.prepare_data(
                     source=self.events[wc.event_idx],
@@ -912,7 +913,7 @@ def get_synthetics(self, point, **kwargs):
                 arrival_times=arrival_times,
                 outmode=outmode,
                 chop_bounds=chop_bounds,
-                nprocs=nprocs,
+                nthreads=nthreads,
                 # plot=True,
                 **kwargs,
             )

beat/plotting/marginals.py

@@ -469,6 +469,17 @@ def remove_var(varnames, varname):
                             except KeyError:
                                 pass
 
+                            except TypeError:
+                                # likelihood returns a single scalar
+                                line = lines[v]
+                                ax.axvline(x=line, color="white", lw=1.0)
+                                ax.axvline(
+                                    x=line,
+                                    color="black",
+                                    linestyle="dashed",
+                                    lw=1.0,
+                                )
+
                         if posterior != "None":
                             if posterior == "all":
                                 for k, idx in posterior_idxs.items():