class Filter(FilterBase):
    """
    Filter object defining frequency range of traces after time-domain
    filtering.
    """

    lower_corner = Float.T(
        default=0.001,
        help='Lower corner frequency')
    upper_corner = Float.T(
        default=0.1,
        help='Upper corner frequency')
    order = Int.T(
        default=4,
        help='order of filter, the higher the steeper')
    stepwise = Bool.T(
        default=True,
        help='If set to true the bandpass filter is done it two'
             ' consecutive steps, first high-pass then low-pass.')

    def apply(self, trace):
        # filter traces
        # stepwise
        if self.stepwise:
            logger.debug('Stepwise HP LP filtering')
            trace.highpass(
                corner=self.lower_corner,
                order=self.order,
                demean=True)
            trace.lowpass(
                corner=self.upper_corner,
                order=self.order,
                demean=False)
        else:
            logger.debug('Single BP filtering')
            trace.bandpass(
                corner_hp=self.lower_corner,
                corner_lp=self.upper_corner,
                order=self.order) class FilterBase(Object):
    pass class ArrivalTaper(trace.Taper):
    """
    Cosine arrival Taper.
    """

    a = Float.T(default=-15.,
                help='start of fading in; [s] w.r.t. phase arrival')
    b = Float.T(default=-10.,
                help='end of fading in; [s] w.r.t. phase arrival')
    c = Float.T(default=50.,
                help='start of fading out; [s] w.r.t. phase arrival')
    d = Float.T(default=55.,
                help='end of fading out; [s] w.r.t phase arrival')

    def check_sample_rate_consistency(self, deltat):
        """
        Check if taper durations are consistent with GF sample rate.
        """
        for chop_b in (['b', 'c'], ['a', 'd']):
            duration = self.duration(chop_b)
            ratio = duration / deltat
            utility.error_not_whole(
                ratio,
                errstr='Taper duration %g of %s is inconsistent with'
                       ' sampling rate of %g! Please adjust Taper values!' % (
                           duration, utility.list2string(chop_b), deltat))

    def duration(self, chop_bounds=['b', 'c']):
        t0 = getattr(self, chop_bounds[0])
        t1 = getattr(self, chop_bounds[1])
        return t1 - t0

    def nsamples(self, sample_rate, chop_bounds=['b', 'c']):
        """
        Returns the number of samples a tapered trace would have given
        its sample rate and chop_bounds

        Parameters
        ----------
        sample_rate : float
        """
        return int(num.ceil(sample_rate * self.duration(chop_bounds)))

    @property
    def fadein(self):
        return self.b - self.a

    @property
    def fadeout(self):
        return self.d - self.c

    def get_pyrocko_taper(self, arrival_time):
        """
        Get pyrocko CosTaper object that may be applied to trace operations.

        Parameters
        ----------
        arrival_time : float
            [s] of the reference time around which the taper will be applied

        Returns
        -------
        :class:`pyrocko.trace.CosTaper`
        """
        if not self.a < self.b < self.c < self.d:
            raise ValueError('Taper values violate: a < b < c < d')

        return trace.CosTaper(
            arrival_time + self.a,
            arrival_time + self.b,
            arrival_time + self.c,
            arrival_time + self.d) class Parameter(Object):
    """
    Optimization parameter determines the bounds of the search space.
    """

    name = String.T(default='depth')
    form = String.T(default='Uniform',
                    help='Type of prior distribution to use. Options:'
                         ' "Uniform", ...')
    lower = Array.T(shape=(None,),
                    dtype=tconfig.floatX,
                    serialize_as='list',
                    default=num.array([0., 0.], dtype=tconfig.floatX))
    upper = Array.T(shape=(None,),
                    dtype=tconfig.floatX,
                    serialize_as='list',
                    default=num.array([1., 1.], dtype=tconfig.floatX))
    testvalue = Array.T(shape=(None,),
                        dtype=tconfig.floatX,
                        serialize_as='list',
                        default=num.array([0.5, 0.5], dtype=tconfig.floatX))

    def validate_bounds(self):

        supported_vars = list(physical_bounds.keys())

        if self.name not in supported_vars:
            candidate = self.name.split('_')[-1]
            if candidate in supported_vars:
                name = candidate
            elif self.name[0:2] == 'h_':
                name = 'hypers'
            elif self.name[0:11] == 'time_shifts':
                name = 'time_shift'
            else:
                raise TypeError(
                    'The parameter "%s" cannot'
                    ' be optimized for!' % self.name)
        else:
            name = self.name

        phys_b = physical_bounds[name]
        if self.lower is not None:
            for i in range(self.dimension):
                if self.upper[i] < self.lower[i]:
                    raise ValueError(
                        'The upper parameter bound for'
                        ' parameter "%s" must be higher than the lower'
                        ' bound' % self.name)

                if self.testvalue[i] > self.upper[i] or \
                        self.testvalue[i] < self.lower[i]:
                    raise ValueError(
                        'The testvalue of parameter "%s" at index "%i" has to'
                        ' be within the upper and lower bounds, please adjust!'
                        % (self.name, i))

                if self.upper[i] > phys_b[1] or \
                        self.lower[i] < phys_b[0]:
                    raise ValueError(
                        'The parameter bounds (%f, %f) for "%s" are outside of'
                        ' physically meaningful values (%f, %f)!' % (
                            self.lower[i], self.upper[i], self.name,
                            phys_b[0], phys_b[1]))
        else:
            raise ValueError(
                'Parameter bounds for "%s" have to be defined!' % self.name)

    def get_upper(self, repeat=1):
        if self.upper.size != num.sum(repeat):
            return list_repeat(self.upper, repeat)
        else:
            return self.upper

    def get_lower(self, repeat=1):
        if self.lower.size != num.sum(repeat):
            return list_repeat(self.lower, repeat)
        else:
            return self.lower

    def get_testvalue(self, repeat=1):
        if self.testvalue.size != num.sum(repeat):
            return list_repeat(self.testvalue, repeat)
        else:
            return self.testvalue

    def random(self, shape=None):
        """
        Create random samples within the parameter bounds.

        Parameters
        ----------
        shape : int or list
            of int number of draws from distribution

        Returns
        -------
        :class:`numpy.ndarray` of size (n, m)
        """
        if shape is None:
            shape = self.dimension

        lower = self.get_lower(shape)
        rands = num.random.rand(num.sum(shape))
        try:
            return (self.get_upper(shape) - lower) * rands + lower
        except ValueError:
            raise ValueError(
                'Value inconsistency shapes: {} parameter '
                'dimension {}'.format(shape, self.dimension))

    @property
    def dimension(self):
        return self.lower.size

    def bound_to_array(self):
        return num.array([self.lower, self.testval, self.upper],
                         dtype=num.float) class ReferenceLocation(gf.Location):
    """
    Reference Location for Green's Function store calculations!
    """
    station = String.T(
        default='Store_Name',
        help='This mimics the station.station attribute which determines the'
             ' store name!') class RectangularSource(gf.RectangularSource):
    """
    Source for rectangular fault that unifies the necessary different source
    objects for teleseismic and geodetic computations.
    Reference point of the depth attribute is the top-center of the fault.

    Many of the methods of the RectangularSource have been modified from
    the HalfspaceTool from GertjanVanZwieten.
    """

    @property
    def dipvector(self):
        """
        Get 3 dimensional dip-vector of the planar fault.

        Parameters
        ----------
        dip : scalar, float
            dip-angle [deg] of the fault
        strike : scalar, float
            strike-abgle [deg] of the fault

        Returns
        -------
        :class:`numpy.ndarray`
        """

        return num.array(
            [num.cos(self.dip * d2r) * num.cos(self.strike * d2r),
             -num.cos(self.dip * d2r) * num.sin(self.strike * d2r),
             num.sin(self.dip * d2r)])

    @property
    def strikevector(self):
        """
        Get 3 dimensional strike-vector of the planar fault.

        Parameters
        ----------
        strike : scalar, float
            strike-abgle [deg] of the fault

        Returns
        -------
        :class:`numpy.ndarray`
        """

        return num.array(
            [num.sin(self.strike * d2r),
             num.cos(self.strike * d2r),
             0.])

    @property
    def center(self):
        """
        Get 3d fault center coordinates. Depth attribute is top depth!

        Parameters
        ----------
        width : scalar, float
            width [m] of the fault (dip-direction)

        Returns
        -------
        :class:`numpy.ndarray` with x, y, z coordinates of the center of the
        fault
        """

        return num.array([self.east_shift, self.north_shift, self.depth]) + \
            0.5 * self.width * self.dipvector

    def center2top_depth(self, center):
        """
        Get top depth of the fault [m] given a potential center point.
        (Patches method needs input depth to
        be top_depth.)

        Parameters
        ----------
        center : vector[x, y, z], float
            coordinates [m] of the center of the fault

        Returns
        -------
        :class:`numpy.ndarray` with x, y, z coordinates of the central
        upper edge of the fault
        """

        return num.array([center[0], center[1], center[2]]) - \
            (0.5 * self.width * self.dipvector)

    @property
    def bottom_center(self):
        """
        Bottom depth of the fault [m].
        (Patches method needs input depth to be top_depth.)

        Returns
        -------
        :class:`numpy.ndarray` with x, y, z coordinates of the central
        lower edge of the fault
        """

        return num.array([self.east_shift, self.north_shift, self.depth]) + \
               (self.width * self.dipvector)

    @property
    def bottom_depth(self):
        return float(self.bottom_center[2])

    @property
    def bottom_left(self):
        return self.bottom_center - (0.5 * self.strikevector * self.length)

    @property
    def bottom_right(self):
        return self.bottom_center + (0.5 * self.strikevector * self.length)

    @property
    def top_left(self):
        return num.array([self.east_shift, self.north_shift, self.depth]) - (
            0.5 * self.strikevector * self.length)

    @property
    def top_right(self):
        return num.array([self.east_shift, self.north_shift, self.depth]) + (
            0.5 * self.strikevector * self.length)

    @property
    def corners(self):
        return num.vstack(
            [self.top_left, self.top_right,
             self.bottom_left, self.bottom_right])

    def trace_center(self):
        """
        Get trace central coordinates of the fault [m] at the surface of the
        halfspace.

        Returns
        -------
        :class:`numpy.ndarray` with x, y, z coordinates of the central
        lower edge of the fault
        """

        bc = self.bottom_center
        xtrace = bc[0] - (
                bc[2] * num.cos(d2r * self.strike) / num.tan(d2r * self.dip))
        ytrace = bc[1] + (
                bc[2] * num.sin(d2r * self.strike) / num.tan(d2r * self.dip))
        return num.array([xtrace, ytrace, 0.])

    def patches(self, nl, nw, datatype):
        """
        Cut source into n by m sub-faults and return n times m
        :class:`RectangularSource` Objects.
        Discretization starts at shallow depth going row-wise deeper.
        REQUIRES: self.depth to be TOP DEPTH!!! Returned faults also have depth
        reference at the top!

        Parameters
        ----------
        nl : int
            number of patches in length direction (strike)
        nw : int
            number of patches in width direction (dip)
        datatype : string
            'geodetic' or 'seismic' determines the source to be returned

        Returns
        -------
        :class:`pscmp.PsCmpRectangularSource` or
        :class:`pyrocko.gf.seismosizer.RectangularSource` depending on
        datatype. Depth is being updated from top_depth to center_depth.
        """

        length = self.length / float(nl)
        width = self.width / float(nw)

        patches = []
        for j in range(nw):
            for i in range(nl):
                sub_top = self.center2top_depth(self.center) + \
                    self.strikevector * ((i + 0.5 - 0.5 * nl) * length) + \
                    self.dipvector * (j * width)

                patch = RectangularSource(
                    lat=float(self.lat),
                    lon=float(self.lon),
                    east_shift=float(sub_top[0]),
                    north_shift=float(sub_top[1]),
                    depth=float(sub_top[2]),
                    strike=self.strike, dip=self.dip, rake=self.rake,
                    length=length, width=width, stf=self.stf,
                    time=self.time, slip=self.slip, anchor='top',
                    opening_fraction=self.opening_fraction)

                patches.append(patch)

        return patches

    def get_n_patches(self, patch_size=1000., dimension='length'):
        """
        Return number of patches along dimension of the fault.

        Parameters
        ----------
        patch_size : float
            patch size [m] of desired sub-patches
        dimension : str

        Returns
        -------
        int
        """
        if dimension not in ['length', 'width']:
            raise ValueError('Invalid dimension!')

        n_p = num.round(self[dimension] / patch_size, decimals=4)
        return int(num.ceil(n_p))

    def extent_source(
            self, extension_width, extension_length,
            patch_width, patch_length):
        """
        Extend fault into all directions. Rounds dimensions to have no
        half-patches.

        Parameters
        ----------
        extension_width : float
            factor to extend source in width (dip-direction)
        extension_length : float
            factor extend source in length (strike-direction)
        patch_width : float
            Width [m] of subpatch in dip-direction
        patch_length : float
            Length [m] of subpatch in strike-direction

        Returns
        -------
        dict with list of :class:`pscmp.PsCmpRectangularSource` or
            :class:`pyrocko.gf.seismosizer.RectangularSource`
        """
        s = copy.deepcopy(self)

        length = self.length
        width = self.width

        if extension_length:
            new_length = num.ceil(
                (length + (2. * length * extension_length)) / km) * km
            npl = int(num.ceil(new_length / patch_length))
            new_length = float(npl * patch_length)
        else:
            new_length = length

        if extension_width:
            new_width = num.ceil(
                (width + (2. * width * extension_width)) / km) * km
            npw = int(num.ceil(new_width / patch_width))
            new_width = float(npw * patch_width)
        else:
            new_width = width

        logger.info(
            'Fault extended to length=%f, width=%f!' % (new_length, new_width))

        orig_center = s.center

        s.update(length=new_length, width=new_width)

        top_center = s.center2top_depth(orig_center)

        if top_center[2] < 0.:
            logger.info('Fault would intersect surface!'
                        ' Setting top center to 0.!')
            trace_center = s.trace_center()
            s.update(east_shift=float(trace_center[0]),
                     north_shift=float(trace_center[1]),
                     depth=float(trace_center[2]))
        else:
            s.update(east_shift=float(top_center[0]),
                     north_shift=float(top_center[1]),
                     depth=float(top_center[2]))

        return s

    @classmethod
    def from_kite_source(cls, source, kwargs):

        d = dict(
            lat=source.lat,
            lon=source.lon,
            north_shift=source.northing,
            east_shift=source.easting,
            depth=source.depth,
            width=source.width,
            length=source.length,
            strike=source.strike,
            dip=source.dip,
            rake=source.rake,
            slip=source.slip,
            anchor='top',
            **kwargs)

        if hasattr(source, 'decimation_factor'):
            d['decimation_factor'] = source.decimation_factor

        return cls(**d) class MTSourceWithMagnitude(gf.SourceWithMagnitude):
    '''
    A moment tensor point source.
    '''

    discretized_source_class = meta.DiscretizedMTSource

    mnn = Float.T(
        default=1.,
        help='north-north component of moment tensor')

    mee = Float.T(
        default=1.,
        help='east-east component of moment tensor')

    mdd = Float.T(
        default=1.,
        help='down-down component of moment tensor')

    mne = Float.T(
        default=0.,
        help='north-east component of moment tensor')

    mnd = Float.T(
        default=0.,
        help='north-down component of moment tensor')

    med = Float.T(
        default=0.,
        help='east-down component of moment tensor')

    def __init__(self, **kwargs):
        if 'm6' in kwargs:
            for (k, v) in zip('mnn mee mdd mne mnd med'.split(),
                              kwargs.pop('m6')):
                kwargs[k] = float(v)

        Source.__init__(self, **kwargs)

    @property
    def m6(self):
        return num.array(self.m6_astuple)

    @property
    def scaled_m6(self):
        m9 = mtm.symmat6(*self.m6)
        m0_unscaled = math.sqrt(num.sum(m9.A ** 2)) / math.sqrt(2.)
        m9 /= m0_unscaled
        m6 = mtm.to6(m9)
        return m6

    @property
    def scaled_m6_dict(self):
        keys = ['mnn', 'mee', 'mdd', 'mne', 'mnd', 'med']
        return {k: m for k, m in zip(keys, self.scaled_m6.tolist())}

    @property
    def m6_astuple(self):
        return (self.mnn, self.mee, self.mdd, self.mne, self.mnd, self.med)

    @m6.setter
    def m6(self, value):
        self.mnn, self.mee, self.mdd, self.mne, self.mnd, self.med = value

    def base_key(self):
        return Source.base_key(self) + self.m6_astuple

    def discretize_basesource(self, store, target=None):
        times, amplitudes = self.effective_stf_pre().discretize_t(
            store.config.deltat, self.time)
        m0 = mtm.magnitude_to_moment(self.magnitude)
        m6s = self.scaled_m6 * m0
        return meta.DiscretizedMTSource(
            m6s=m6s[num.newaxis, :] * amplitudes[:, num.newaxis],
            **self._dparams_base_repeated(times))

    def pyrocko_moment_tensor(self):
        return mtm.MomentTensor(m=mtm.symmat6(*self.m6_astuple) * self.moment)

    def pyrocko_event(self, **kwargs):
        mt = self.pyrocko_moment_tensor()
        return Source.pyrocko_event(
            self,
            moment_tensor=self.pyrocko_moment_tensor(),
            magnitude=float(mt.moment_magnitude()),
            **kwargs)

    @classmethod
    def from_pyrocko_event(cls, ev, **kwargs):
        d = {}
        mt = ev.moment_tensor
        if mt:
            d.update(m6=list(map(float, mt.m6())))

        d.update(kwargs)
        return super(MTSourceWithMagnitude, cls).from_pyrocko_event(ev, **d)

    def get_derived_parameters(self, store=None, target=None):
        mt = mtm.MomentTensor.from_values(self.scaled_m6)
        return num.hstack(mt.both_strike_dip_rake())

    @property
    def nderived_parameters(self):
        return self.get_derived_parameters().size class MTQTSource(gf.SourceWithMagnitude):
    """
    A moment tensor point source.

    Notes
    -----
    Following Q-T parameterization after Tape & Tape 2015
    """

    discretized_source_class = meta.DiscretizedMTSource

    w = Float.T(
        default=0.,
        help='Lune latitude delta transformed to grid. '
             'Defined: -3/8pi <= w <=3/8pi. '
             'If fixed to zero the MT is deviatoric.')

    v = Float.T(
        default=0.,
        help='Lune co-longitude transformed to grid. '
             'Definded: -1/3 <= v <= 1/3. '
             'If fixed to zero together with w the MT is pure DC.')

    kappa = Float.T(
        default=0.,
        help='Strike angle equivalent of moment tensor plane.'
             'Defined: 0 <= kappa <= 2pi')

    sigma = Float.T(
        default=0.,
        help='Rake angle equivalent of moment tensor slip angle.'
             'Defined: -pi/2 <= sigma <= pi/2')

    h = Float.T(
        default=0.,
        help='Dip angle equivalent of moment tensor plane.'
             'Defined: 0 <= h <= 1')

    def __init__(self, **kwargs):
        n = 1000
        self._beta_mapping = num.linspace(0, pi, n)
        self._u_mapping = \
            (3. / 4. * self._beta_mapping) - \
            (1. / 2. * num.sin(2. * self._beta_mapping)) + \
            (1. / 16. * num.sin(4. * self._beta_mapping))

        self.lambda_factor_matrix = num.array(
            [[sqrt3, -1., sqrt2],
             [0., 2., sqrt2],
             [-sqrt3, -1., sqrt2]], dtype='float64')

        self.R = get_rotation_matrix()
        self.roty_pi4 = self.R['y'](-pi4)
        self.rotx_pi = self.R['x'](pi)

        self._lune_lambda_matrix = num.zeros((3, 3), dtype='float64')

        Source.__init__(self, **kwargs)

    @property
    def u(self):
        """
        Lunar co-latitude(beta), dependend on w
        """
        return (3. / 8.) * num.pi - self.w

    @property
    def gamma(self):
        """
        Lunar longitude, dependend on v
        """
        return v_to_gamma(self.v)

    @property
    def beta(self):
        """
        Lunar co-latitude, dependend on u
        """
        return w_to_beta(
            self.w, u_mapping=self._u_mapping, beta_mapping=self._beta_mapping)

    def delta(self):
        """
        From Tape & Tape 2012, delta measures departure of MT being DC
        Delta = Gamma = 0 yields pure DC
        """
        return (pi / 2.) - self.beta

    @property
    def rho(self):
        return mtm.magnitude_to_moment(self.magnitude) * sqrt2

    @property
    def theta(self):
        return num.arccos(self.h)

    @property
    def rot_theta(self):
        return self.R['x'](self.theta)

    @property
    def rot_kappa(self):
        return self.R['z'](-self.kappa)

    @property
    def rot_sigma(self):
        return self.R['z'](self.sigma)

    @property
    def lune_lambda(self):
        sin_beta = num.sin(self.beta)
        cos_beta = num.cos(self.beta)
        sin_gamma = num.sin(self.gamma)
        cos_gamma = num.cos(self.gamma)
        vec = num.array([sin_beta * cos_gamma, sin_beta * sin_gamma, cos_beta])
        return 1. / sqrt6 * self.lambda_factor_matrix.dot(vec) * self.rho

    @property
    def lune_lambda_matrix(self):
        num.fill_diagonal(self._lune_lambda_matrix, self.lune_lambda)
        return self._lune_lambda_matrix

    @property
    def rot_V(self):
        return self.rot_kappa.dot(self.rot_theta).dot(self.rot_sigma)

    @property
    def rot_U(self):
        return self.rot_V.dot(self.roty_pi4)

    @property
    def m9_nwu(self):
        """
        MT orientation is in NWU
        """
        return self.rot_U.dot(
            self.lune_lambda_matrix).dot(num.linalg.inv(self.rot_U))

    @property
    def m9(self):
        """
        Pyrocko MT in NED
        """
        return self.rotx_pi.dot(self.m9_nwu).dot(self.rotx_pi.T)

    @property
    def m6(self):
        return mtm.to6(self.m9)

    @property
    def m6_astuple(self):
        return tuple(self.m6.ravel().tolist())

    def base_key(self):
        return Source.base_key(self) + self.m6_astuple

    def discretize_basesource(self, store, target=None):
        times, amplitudes = self.effective_stf_pre().discretize_t(
            store.config.deltat, self.time)
        return meta.DiscretizedMTSource(
            m6s=self.m6[num.newaxis, :] * amplitudes[:, num.newaxis],
            **self._dparams_base_repeated(times))

    def pyrocko_moment_tensor(self):
        return mtm.MomentTensor(m=mtm.symmat6(*self.m6_astuple) * self.moment)

    def pyrocko_event(self, **kwargs):
        mt = self.pyrocko_moment_tensor()
        return Source.pyrocko_event(
            self,
            moment_tensor=self.pyrocko_moment_tensor(),
            magnitude=float(mt.moment_magnitude()),
            **kwargs)

    @classmethod
    def from_pyrocko_event(cls, ev, **kwargs):
        d = {}
        mt = ev.moment_tensor
        if mt:
            logger.warning(
                'From event will ignore MT components initially. '
                'Needs mapping from NED to QT space!')
            # d.update(m6=list(map(float, mt.m6())))

        d.update(kwargs)
        return super(MTQTSource, cls).from_pyrocko_event(ev, **d)

    def get_derived_parameters(self, store=None, target=None):
        """
        Returns array with mt components and dc component conversions
        """
        scaled_m6 = self.m6 / self.moment
        mt = mtm.MomentTensor.from_values(scaled_m6)
        return num.hstack((scaled_m6, num.hstack(mt.both_strike_dip_rake())))

    @property
    def nderived_parameters(self):
        return self.get_derived_parameters().size

    def __getstate__(self):
        state = self.__dict__.copy()
        state['R'] = None
        return state

    def __setstate__(self, state):
        self.__dict__.update(state)
        self.R = get_rotation_matrix() def available_noise_structures():
    return list(NoiseStructureCatalog.keys()) def argsorted(seq, reverse=False):
    def __init__(self):
    def __call__(self, string, multiplier=1):
    def __getitem__(self, key):
    def reset(self, string=None):
    def __init__(self, list_arrays, intype='numpy'):
    def __getitem__(self, key):
    def __iter__(self):
    def variables(self):
    def __init__(self, ordering, list_arrays, blacklist=[]):
    def d2l(self, dpt):
    def l2d(self, a_list):
    def l2a(self, list_arrays):
    def f3map(self, list_arrays):
    def a2l(self, array):
    def a_nd2l(self, array):
    def srmap(self, tarray):
def weed_input_rvs(input_rvs, mode, datatype):
def apply_station_blacklist(stations, blacklist):
def weed_data_traces(data_traces, stations):
def weed_targets(targets, stations, discard_targets=[]):
def downsample_trace(data_trace, deltat=None, snap=False):
def weed_stations(
        stations, event, distances=(30., 90.), remove_duplicate=False):
def transform_sources(sources, datatypes, decimation_factors=None):
def adjust_point_units(point):
def split_point(point):
def join_points(ldicts):
def check_point_keys(point, phrase):
def update_source(source, **point):
def setup_logging(project_dir, levelname, logfilename='BEAT_log.txt'):
def search_catalog(date, min_magnitude, dayrange=1.):
def RS_dipvector(source):
def strike_vector(strike, order='ENZ'):
def RS_strikevector(source):
def RS_center(source):
def adjust_fault_reference(source, input_depth='top'):
def dump_objects(outpath, outlist):
def load_objects(loadpath):
def ensure_cov_psd(cov):
def near_psd(x, epsilon=num.finfo(num.float64).eps):
def repair_covariance(x, epsilon=num.finfo(num.float64).eps):
def running_window_rms(data, window_size, mode='valid'):
def list2string(l, fill=', '):
def unique_list(l):
def join_models(global_model, crustal_model):
def split_off_list(l, off_length):
def mod_i(i, cycle):
def biggest_common_divisor(a, b):
def gather(l, key, sort=None, filter=None):
def get_fit_indexes(llk):
def check_hyper_flag(problem):
def error_not_whole(f, errstr=''):
def scalar2floatX(a, floatX=tconfig.floatX):
def scalar2int(a, floatX=tconfig.floatX):
def PsGrnArray2LayeredModel(psgrn_input_path):
def swap_columns(array, index1, index2):
def line_intersect(e1, e2, n1, n2):
def get_rotation_matrix(axes=['x', 'y', 'z']):
    def rotx(angle):
    def roty(angle):
    def rotz(angle):
def get_random_uniform(lower, upper, dimension=1):
def positions2idxs(
        positions, cell_size, min_pos=0., backend=num, dtype='int16'):
def rotate_coords_plane_normal(coords, sf):
def time_method(loop=10000):
    def timer_decorator(func):
        def wrap_func(*args, **kwargs):
def is_odd(value):
def is_even(value):
def get_data_radiant(data):
def find_elbow(data, theta=None, rotate_left=False):
class Counter(object):
class ListArrayOrdering(object):
class ListToArrayBijection(object):
    T = 1. / (num.multiply(vec, vec) * val.T)
    T = num.matrix(num.sqrt(num.diag(num.array(T).reshape((n)))))
    B = T * vec * num.diag(num.array(num.sqrt(val)).reshape((n)))
    R = {'x': rotx,
         'y': roty,
         'z': rotz}
import logging
import os
import numpy as num
from collections import OrderedDict
from pyrocko.guts import Object, List, String, Float, Int, Tuple, Bool, Dict
from pyrocko.guts import load, dump, StringChoice, ArgumentError
from pyrocko.cake import load_model
from pyrocko import trace, model, util, gf
from pyrocko.gf import RectangularSource as PyrockoRS
from pyrocko.gf.seismosizer import Cloneable, stf_classes
from beat.heart import Filter, FilterBase, ArrivalTaper, Parameter
from beat.heart import ReferenceLocation
from beat.sources import RectangularSource, MTSourceWithMagnitude, MTQTSource
from beat.covariance import available_noise_structures
from beat import utility
from theano import config as tconfig
        from beat.models.corrections import EulerPoleCorrection
        from beat.models.corrections import StrainRateCorrection
        from beat.models.corrections import RampCorrection
        from beat.inputf import load_kite_scenes
        from beat.inputf import load_and_blacklist_gnss
        from pymc3 import Uniform
_quantity_choices = ['displacement', 'velocity', 'acceleration']
_interpolation_choices = ['nearest_neighbor', 'multilinear']
_structure_choices = available_noise_structures()
_mode_choices = [geometry_mode_str, ffi_mode_str]
_regularization_choices = ['laplacian', 'none']
_correlation_function_choices = ['nearest_neighbor', 'gaussian', 'exponential']
_discretization_choices = ['uniform', 'resolution']
_initialization_choices = ['random', 'lsq']
_backend_choices = ['csv', 'bin']
_datatype_choices = ['geodetic', 'seismic']
_sampler_choices = ['PT', 'SMC', 'Metropolis']


class InconsistentParameterNaming(Exception):
    context = '"{}" Parameter name is inconsistent with its key "{}"!\n' + \
              ' Please edit the config_{}.yaml accordingly!'

    def __init__(self, keyname='', parameter='', mode=''):
        self.parameter = parameter
        self.mode = mode
        self.keyname = keyname

    def __str__(self):
        return self.context.format(self.parameter, self.keyname, self.mode)


class ConfigNeedsUpdatingError(Exception):
    context = 'Configuration file has to be updated! \n' + \
              ' Please run "beat update <project_dir --parameters=structure>'

    def __init__(self, errmess=''):
        self.errmess = errmess

    def __str__(self):
        return '\n%s\n%s' % (self.errmess, self.context)


class GFConfig(Object):
    """
    Base config for GreensFunction calculation parameters.
    """
    store_superdir = String.T(
        default='./',
        help='Absolute path to the directory where Greens Function'
             ' stores are located')
    reference_model_idx = Int.T(
        default=0,
        help='Index to velocity model to use for the optimization.'
             ' 0 - reference, 1..n - model of variations')
    n_variations = Tuple.T(
        2,
        Int.T(),
        default=(0, 1),
        help='Start and end index to vary input velocity model. '
             'Important for the calculation of the model prediction covariance'
             ' matrix with respect to uncertainties in the velocity model.')
    earth_model_name = String.T(
        default='ak135-f-continental.f',
        help='Name of the reference earthmodel, see '
             'pyrocko.cake.builtin_models() for alternatives.')
    nworkers = Int.T(
        default=1,
        help='Number of processors to use for calculating the GFs')


class NonlinearGFConfig(GFConfig):
    """
    Config for non-linear GreensFunction calculation parameters.
    Defines how the grid of Green's Functions in the respective store is
    created.
    """

    use_crust2 = Bool.T(
        default=True,
        help='Flag, for replacing the crust from the earthmodel'
             'with crust from the crust2 model.')
    replace_water = Bool.T(
        default=True,
        help='Flag, for replacing water layers in the crust2 model.')
    custom_velocity_model = gf.Earthmodel1D.T(
        default=None,
        optional=True,
        help='Custom Earthmodel, in case crust2 and standard model not'
             ' wanted. Needs to be a :py::class:cake.LayeredModel')
    source_depth_min = Float.T(
        default=0.,
        help='Minimum depth [km] for GF function grid.')
    source_depth_max = Float.T(
        default=10.,
        help='Maximum depth [km] for GF function grid.')
    source_depth_spacing = Float.T(
        default=1.,
        help='Depth spacing [km] for GF function grid.')
    source_distance_radius = Float.T(
        default=20.,
        help='Radius of distance grid [km] for GF function grid around '
             'reference event.')
    source_distance_spacing = Float.T(
        default=1.,
        help='Distance spacing [km] for GF function grid w.r.t'
             ' reference_location.')
    error_depth = Float.T(
        default=0.1,
        help='3sigma [%/100] error in velocity model layer depth, '
             'translates to interval for varying the velocity model')
    error_velocities = Float.T(
        default=0.1,
        help='3sigma [%/100] in velocity model layer wave-velocities, '
             'translates to interval for varying the velocity model')
    depth_limit_variation = Float.T(
        default=600.,
        help='Depth limit [km] for varying the velocity model. Below that '
             'depth the velocity model is not varied based on the errors '
             'defined above!')


class SeismicGFConfig(NonlinearGFConfig):
    """
    Seismic GF parameters for Layered Halfspace.
    """

    reference_location = ReferenceLocation.T(