Rotation proposal

beat/covariance.py

@@ -1,11 +1,11 @@
-from pyrocko import gf, trace
+from pyrocko import trace
 
 import numpy as num
 from time import time
 from scipy.linalg import toeplitz
 
 import logging
-import copy
+from copy import deepcopy
 
 from beat import heart
 from beat.utility import ensure_cov_psd, running_window_rms, list2string
@@ -50,6 +50,7 @@ def exponential_data_covariance(n, dt, tzero):
         -num.abs(num.arange(n)[:, num.newaxis] -
                  num.arange(n)[num.newaxis, :]) * dt / tzero)
 
+
 def identity_data_covariance(n, dt=None, tzero=None):
     """
     Get identity covariance matrix.
@@ -277,128 +278,6 @@ def get_data_covariances(self, wmap, sample_rate, results=None):
         return cov_ds
 
 
-def model_prediction_sensitivity(engine, *args, **kwargs):
-    '''
-    Calculate the model prediction Covariance Sensitivity Kernel.
-    (numerical derivation with respect to the input source parameter(s))
-    Following Duputel et al. 2014
-
-    :Input:
-    :py:class:'engine'
-    source_parms = list of parameters with respect to which the kernel
-                   is being calculated e.g. ['strike', 'dip', 'depth']
-    !!!
-    NEEDS to have seismosizer source object parameter variable name convention
-    !!!
-    (see seismosizer.source.keys())
-
-    calculate_model_prediction_sensitivity(request, source_params, **kwargs)
-    calculate_model_prediction_sensitivity(sources,
-                                             targets, source_params, **kwargs)
-
-    Returns traces in a list[parameter][targets] for each station and channel
-    as specified in the targets. The location code of each trace is placed to
-    show the respective source parameter.
-    '''
-
-    if len(args) not in (0, 1, 2, 3):
-        raise gf.BadRequest('invalid arguments')
-
-    if len(args) == 2:
-        kwargs['request'] = args[0]
-        kwargs['source_params'] = args[1]
-
-    elif len(args) == 3:
-        kwargs.update(gf.Request.args2kwargs(args[0:1]))
-        kwargs['source_params'] = args[2]
-
-    request = kwargs.pop('request', None)
-    nprocs = kwargs.pop('nprocs', 1)
-    source_params = kwargs.pop('source_params', None)
-    h = kwargs.pop('h', None)
-
-    if request is None:
-        request = gf.Request(**kwargs)
-
-    if h is None:
-        h = num.ones(len(source_params)) * 1e-1
-
-    # create results list
-    sensitivity_param_list = []
-    sensitivity_param_trcs = []
-
-    for i in range(len(source_params)):
-        sensitivity_param_list.append([0] * len(request.targets))
-        sensitivity_param_trcs.append([0] * len(request.targets))
-
-    for ref_source in request.sources:
-        par_count = 0
-        for param in source_params:
-            print(param, 'with h = ', h[par_count])
-            calc_source_p2h = ref_source.clone()
-            calc_source_ph = ref_source.clone()
-            calc_source_mh = ref_source.clone()
-            calc_source_m2h = ref_source.clone()
-
-            setattr(calc_source_p2h, param,
-                    ref_source[param] + (2 * h[par_count]))
-            setattr(calc_source_ph, param,
-                    ref_source[param] + (h[par_count]))
-            setattr(calc_source_mh, param,
-                    ref_source[param] - (h[par_count]))
-            setattr(calc_source_m2h, param,
-                    ref_source[param] - (2 * h[par_count]))
-
-            calc_sources = [calc_source_p2h, calc_source_ph,
-                            calc_source_mh, calc_source_m2h]
-
-            response = engine.process(sources=calc_sources,
-                                      targets=request.targets,
-                                      nprocs=nprocs)
-
-            for k in range(len(request.targets)):
-                # zero padding if necessary
-                trc_lengths = num.array(
-                    [len(response.results_list[i][k].trace.data) for i in
-                     range(len(response.results_list))])
-                Id = num.where(trc_lengths != trc_lengths.max())
-
-                for l in Id[0]:
-                    response.results_list[l][k].trace.data = num.concatenate(
-                            (response.results_list[l][k].trace.data,
-                             num.zeros(trc_lengths.max() - trc_lengths[l])))
-
-                # calculate numerical partial derivative for
-                # each source and target
-                sensitivity_param_list[par_count][k] = (
-                        sensitivity_param_list[par_count][k] + (
-                            - response.results_list[0][k].trace.data +
-                            8 * response.results_list[1][k].trace.data -
-                            8 * response.results_list[2][k].trace.data +
-                            response.results_list[3][k].trace.data) /
-                        (12 * h[par_count])
-                    )
-
-            par_count = par_count + 1
-
-    # form traces from sensitivities
-    par_count = 0
-    for param in source_params:
-        for k in range(len(request.targets)):
-            sensitivity_param_trcs[par_count][k] = trace.Trace(
-                        network=request.targets[k].codes[0],
-                        station=request.targets[k].codes[1],
-                        ydata=sensitivity_param_list[par_count][k],
-                        deltat=response.results_list[0][k].trace.deltat,
-                        tmin=response.results_list[0][k].trace.tmin,
-                        channel=request.targets[k].codes[3],
-                        location=param)
-
-        par_count = par_count + 1
-
-    return sensitivity_param_trcs
-
-
 def seismic_cov_velocity_models(
         engine, sources, targets, arrival_taper, arrival_time,
         wavename, filterer, plot=False, n_jobs=1):

beat/heart.py

@@ -2831,6 +2831,121 @@ def seis_synthetics(
         raise TypeError('Outmode %s not supported!' % outmode)
 
 
+def seis_derivative(
+        engine, sources, targets, arrival_taper, arrival_times,
+        wavename, filterer, h, parameter, stencil_order=3):
+    """
+    Calculate the model prediction Covariance Sensitivity Kernel.
+    Numerical derivation with respect to the input source 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
+        of parameters with respect to which the kernel
+        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)
+    if parameter not in sources[0].keys():
+        raise AttributeError(
+            'Parameter for which the derivative was requested is not'
+            ' represented by the source.')
+
+    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
+    tmp = num.zeros((len(stencil), nsamples, ntargets), dtype='float64')
+    for i, hstep in enumerate(stencil.hsteps):
+        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, :, :], _ = seis_synthetics(
+            engine=engine,
+            sources=diff_sources,
+            targets=targets,
+            arrival_taper=arrival_taper,
+            wavename=wavename,
+            filterer=filterer,
+            arrival_times=arrival_times,
+            pre_stack_cut=True,
+            outmode='array',
+            chop_bounds=['b', 'c'])
+
+    return tmp * stencil.coefficients / stencil.denominator
+
+
+def seis_jacobian(
+        engine, sources, targets, arrival_taper, arrival_times, wavename,
+        filterer, nderivatives=30):
+    """
+    Calculate the Jacobian with respect to an input source parameter
+
+    Involves numerical calculation of the derivatives. The stable increment h
+    is automatically determined.
+
+    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`
+    parameter : str
+        of parameters with respect to which the kernel
+        is being calculated e.g. 'strike'/ 'dip'/ 'depth'
+    nderivatives : int
+        number of derivatives to test, selects the best from a stable plateau
+
+    Returns
+    -------
+    :class:`num.array` ntargets x nsamples with the first derivative
+    """
+    for h in hs:
+        jacob = seis_derivative(
+            engine, sources, targets, arrival_taper, arrival_times, wavename,
+            filterer, h, parameter, stencil_order=3)
+
+    return jacob
+
+
 def geo_synthetics(
         engine, targets, sources, outmode='stacked_array', plot=False,
         nprocs=1):

beat/sampler/base.py

@@ -4,13 +4,15 @@
 import shutil
 
 from beat import parallel
+
 from beat.backend import check_multitrace, load_multitrace, backend_catalog
-from beat.utility import list2string
+from beat.utility import list2string, ensure_cov_psd
 
 from numpy.random import seed, randint
 from numpy.random import normal, standard_cauchy, standard_exponential, \
     poisson
 import numpy as np
+from scipy import linalg
 
 from theano import function
 
@@ -158,9 +160,21 @@ def __call__(self, num_draws=None):
 
 
 class MultivariateNormalProposal(Proposal):
+    def __init__(self, scale):
+        n, m = scale.shape
+        if n != m:
+            raise ValueError("Covariance matrix is not symmetric.")
+        self.n = n
+        s = ensure_cov_psd(scale)
+        self.chol = linalg.cholesky(s, lower=True)
+
     def __call__(self, num_draws=None):
-        return np.random.multivariate_normal(
-            mean=np.zeros(self.scale.shape[0]), cov=self.scale, size=num_draws)
+        if num_draws is not None:
+            b = np.random.randn(self.n, num_draws)
+            return np.dot(self.chol, b).T
+        else:
+            b = np.random.randn(self.n)
+            return np.dot(self.chol, b)
 
 
 class MultivariateStudentTProposal(Proposal):
@@ -176,22 +190,83 @@ class MultivariateCauchyProposal(Proposal):
     of freedom equal to one.
     """
     def __call__(self, num_draws=None):
+        if not num_draws:
+            num_draws = 1
+
         return multivariate_t_rvs(
             mean=np.zeros(self.scale.shape[0]),
             cov=self.scale, df=1, size=num_draws)
 
 
+class RotationProposal(Proposal):
+    """
+    Proposal for proposing samples in the direction of the principal
+    components.
+
+    Parameters
+    ----------
+    cov: :class:`num.NdArray`
+        of sample covariance matrix
+    """
+    def __init__(self, scale):
+
+        cov = ensure_cov_psd(scale)
+        self.n = cov.shape[0]
+        vec, pc, _ = np.linalg.svd(cov)
+        self.pcsd = 0.5 * (1. / np.sqrt(np.abs(pc)))
+        self.u = vec
+        self.proposal = None
+
+    def __call__(self, num_draws=None):
+
+        if self.proposal is None:
+            raise NotImplementedError(
+                'Please use RotationProposal through its inheriting classes!')
+        else:
+            step = self.proposal(num_draws)
+            delta = 1.1 * self.pcsd * step
+            return self.u.dot(delta.T).T
+
+
+class MultivariateRotationCauchyProposal(RotationProposal):
+    """
+    Proposes steps in the direction of principal component of the given
+    sample covariance using a Cauchy distribution.
+    """
+    def __init__(self, scale):
+
+        super(MultivariateRotationCauchyProposal, self).__init__(scale)
+        self.proposal = CauchyProposal(np.ones((self.n)))
+
+
+class MultivariateRotationNormalProposal(RotationProposal):
+    """
+    Proposes steps in the direction of principal component of the given
+    sample covariance using a Normal distribution.
+    """
+    def __init__(self, scale):
+
+        super(MultivariateRotationNormalProposal, self).__init__(scale)
+        self.proposal = NormalProposal(np.ones((self.n)))
+
+
 proposal_distributions = {
     'Cauchy': CauchyProposal,
     'Poisson': PoissonProposal,
     'Normal': NormalProposal,
     'Laplace': LaplaceProposal,
     'MultivariateNormal': MultivariateNormalProposal,
     'MultivariateCauchy': MultivariateCauchyProposal,
+    'MultivariateRotationNormal': MultivariateRotationNormalProposal,
+    'MultivariateRotationCauchy': MultivariateRotationCauchyProposal,
     'DiscreteBoundedUniform': DiscreteBoundedUniformProposal}
 
 
-multivariate_proposals = ['MultivariateCauchy', 'MultivariateNormal']
+multivariate_proposals = [
+    'MultivariateCauchy',
+    'MultivariateNormal',
+    'MultivariateRotationNormal',
+    'MultivariateRotationCauchy']
 
 
 def available_proposals():