Adding polarity package_V1

hvasbath · Jul 31, 2021 · 6c13773 · 6c13773
1 parent 51a64ee
commit 6c13773
Show file tree

Hide file tree

Showing 10 changed files with 109,024 additions and 5 deletions.
diff --git a/beat/config.py b/beat/config.py
@@ -99,6 +99,10 @@
     'geodetic': interseismic_vars}
 
 geometry_catalog = {
+##Mahdi
+    'polarity': source_catalog,
+##
+
     'geodetic': source_catalog,
     'seismic': source_catalog}
 
@@ -231,7 +235,7 @@ def multi_event_seismic_data_name(nevent=0):
 _discretization_choices = ['uniform', 'resolution']
 _initialization_choices = ['random', 'lsq']
 _backend_choices = ['csv', 'bin']
-_datatype_choices = ['geodetic', 'seismic']
+_datatype_choices = ['geodetic', 'seismic', 'polarity']         ##Mahdi
 _sampler_choices = ['PT', 'SMC', 'Metropolis']
 
 
@@ -696,6 +700,25 @@ def init_waveforms(self, wavenames=['any_P']):
         for wavename in wavenames:
             self.waveforms.append(WaveformFitConfig(name=wavename))
 
+##Mahdi
+class PolarityConfig(Object):
+    st_polarity = List.T(default=[])
+    vmdir = String.T(default='./', 
+                     optional=True,
+                     help="If not given, velocity model from Green's function will be extract")
+    # def __init__(self, st_polarity, vmdir=[]):
+        # if not vmdir:
+
+    def get_hypernames(self):
+        return []
+    def obspol(self):
+        obs_ampl = num.zeros(len(self.st_polarity))
+        station = [None] * obs_ampl.size
+        for i, st_pol in enumerate(self.st_polarity):
+            obs_ampl[i] = float(st_pol.split(" ")[-1])
+            station[i] = st_pol.split(" ")[0]
+        return num.array(station), obs_ampl
+##
 
 class CorrectionConfig(Object):
 
@@ -1592,6 +1615,8 @@ def _mapid(self):
 
 
 datatype_catalog = {
+    ##Mahdi
+    'polarity': PolarityConfig,
     'geodetic': GeodeticConfig,
     'seismic': SeismicConfig}
 
@@ -1618,6 +1643,9 @@ class BEATconfig(Object, Cloneable):
         default=None, optional=True)
     seismic_config = SeismicConfig.T(
         default=None, optional=True)
+    ##Mahdi
+    polarity_config = PolarityConfig.T(default=None, optional=True)
+    ##
     sampler_config = SamplerConfig.T(default=SamplerConfig.D())
     hyper_sampler_config = SamplerConfig.T(
         default=SamplerConfig.D(),
@@ -1849,7 +1877,13 @@ def init_config(name, date=None, min_magnitude=6.0, main_path='./',
                 c.seismic_config.gf_config.replace_water = False
         else:
             c.seismic_config = None
+        ##Mahdi
+        if 'polarity' in datatypes:
+            c.polarity_config = PolarityConfig()
 
+        else:
+            c.polarity_config = None
+        ##
     elif mode == ffi_mode_str:
 
         if source_type != 'RectangularSource':

diff --git a/beat/heart.py b/beat/heart.py
@@ -1235,6 +1235,53 @@ class GeodeticResult(Object):
     processed_res = Array.T(shape=(None,), dtype=num.float, optional=True)
     llk = Float.T(default=0., optional=True)
 
+##Mahdi
+class PolarityTarget(gf.meta.Receiver):
+    '''
+    A seismogram computation request for a single component, including
+    its post-processing parmeters.
+    '''
+
+    codes = Tuple.T(
+        4, String.T(), default=('', 'STA', '', 'Z'),
+        help='network, station, location and channel codes to be set on '
+             'the response trace.')
+
+    elevation = Float.T(
+        default=0.0,
+        help='station surface elevation in [m]')
+
+    vmodel = gf.meta.StringID.T(
+        optional=True,
+        help='ID of Green\'s function store to use for the computation. '
+             'If not given, the processor may use a system default.')
+
+    azimuth = Float.T(
+        optional=True,
+        help='azimuth of sensor component in [deg], clockwise from north. '
+             'If not given, it is guessed from the channel code.')
+    distance = Float.T(
+        optional=True,
+        help='azimuth of sensor component in [deg], clockwise from north. '
+             'If not given, it is guessed from the channel code.')
+    takeoff_angle = Float.T(
+        optional=True,
+        help='azimuth of sensor component in [deg], clockwise from north. '
+             'If not given, it is guessed from the channel code.')
+
+def init_polarity_targets(stations, polarityconfig):
+    D2R = num.pi / 180.0
+    crust_ind = 0
+    obs_stations = polarityconfig.obspol()[0]
+    targets = []
+    for i, station in enumerate(stations):
+        if station.get_channel('Z') and station.station in obs_stations:
+            targets.append(PolarityTarget(
+                            codes=(station.network, station.station,'%i' % crust_ind, 'Z'),  # n, s, l, c
+                            lat=station.lat, lon=station.lon, azimuth=station.azimuth*D2R, distance=station.dist_deg,
+                            elevation=float(station.elevation), vmodel=polarityconfig.vmdir))
+    return targets
+##
 
 def init_seismic_targets(
         stations, earth_model_name='ak135-f-average.m', channels=['T', 'Z'],
@@ -2254,7 +2301,32 @@ def get_phase_arrival_time(engine, source, target, wavename=None, snap=True):
         deltat = 1. / store.config.sample_rate
         atime = trace.t2ind(atime, deltat, snap=round) * deltat
     return atime
-
+##Mahdi
+def get_takeoff_angle(source, target):
+
+    D2R = num.pi / 180.0
+    # try:
+    #     mod = target.vmodel.earthmodel_1d
+    # except:
+    mod = cake.load_model(target.vmodel)
+
+    if mod.discontinuity('moho').name:
+        definition = 'p,P,p\\,P\\,Pv_(cmb)p'
+    else:
+        definition = 'p,P,p\\,P\\'
+
+    dist = target.distance_to(source)
+
+    phases = gf.TPDef(id='any_P', definition=definition)
+
+    # Calculate arrivals
+    rays = mod.arrivals(
+        phases=phases.phases,
+        distances=[dist*cake.m2d],
+        zstart=source.depth,
+        zstop=target.depth)
+    return rays[0].takeoff_angle() * D2R
+##
 
 def get_phase_taperer(
         engine, source, wavename, target, arrival_taper, arrival_time=num.nan):
@@ -2746,6 +2818,31 @@ def get_waveform_mapping(
             channels=channels,
             deltat=self._deltat)
 
+##Mahdi
+class DataPolarityCollection(object):
+
+    def __init__(self, observed_stations, observed_amps):
+        self.stations = num.array(observed_stations)
+        self.amplitudes = observed_amps
+        self._targets = OrderedDict()
+        self._target2index = None
+        self._station2index = None
+
+    @property
+    def n_t(self):
+        return len(self._targets.keys())
+
+    def add_targets(self, targets, replace=False, force=False):
+        DataWaveformCollection.add_targets(self, targets, replace=False, force=False)
+    def sorting(self):
+        amplitudes_inorder = num.zeros(self.n_t)
+        for i, target in enumerate(self._targets):
+            amplitudes_inorder[i] = self.amplitudes[target[1] == self.stations]
+        self.amplitudes = amplitudes_inorder
+    def prepare_data(self, source):
+        for i, target in enumerate(self._targets):
+            self._targets[target].takeoff_angle = get_takeoff_angle(source, self._targets[target])   
+##
 
 def concatenate_datasets(datasets):
     """
@@ -2780,6 +2877,16 @@ def concatenate_datasets(datasets):
     los_vectors = Bij.f3map(_lv_list).astype(tconfig.floatX)
     return datasets, los_vectors, odws, Bij
 
+##Mahdi
+def pol_datahandler(polarity_config, seismic_data_path):
+    stations = utility.load_objects(seismic_data_path)[0]
+    targets = init_polarity_targets(stations, polarity_config)
+    obs_sts, obs_amps = polarity_config.obspol()
+    datahandler = DataPolarityCollection(obs_sts, obs_amps)
+    datahandler.add_targets(targets)
+    datahandler.sorting()
+    return datahandler
+##
 
 def init_datahandler(
         seismic_config, seismic_data_path='./', responses_path=None):
@@ -3145,6 +3252,43 @@ def seis_synthetics(
     else:
         raise TypeError('Outmode %s not supported!' % outmode)
 
+##Mahdi
+
+def GAMMA(takeoffangles, azimuths):
+    return num.array([num.sin(takeoffangles)*num.cos(azimuths),
+                        num.sin(takeoffangles)*num.sin(azimuths),
+                        num.cos(takeoffangles)])
+# def p_amplitude(moment, GAMMAS):
+#     GAMMAS = GAMMAS.reshape(GAMMAS.size,1)
+#     temp = num.diag(GAMMAS.T.dot(moment).dot(GAMMAS))
+#     return temp
+
+def syn_polarity(sources, targets):  
+    toas = num.array([tr.takeoff_angle for tr in targets])
+    azis = num.array([tr.azimuth for tr in targets])
+    gamma_co = GAMMA(toas, azis)
+    pols = gamma_co.T.dot(sources[0].m9).dot(gamma_co).copy()
+    pols[pols>=0] = 1
+    pols[pols<0] = -1
+
+
+    # from theano import tensor as tt
+    # pred_pol = tt.zeros(len(sources)*len(targets))
+    # for source in sources:
+    #     for i, tr in enumerate(targets):
+    #         gamma_co = GAMMA(tr.takeoff_angle, tr.azimuth)
+    #         tmp1 = gamma_co.T.dot(sources[0].m9).dot(gamma_co)
+    #         if tmp1 == 0:                
+    #             tmp2 =  1
+    # pols = []
+    # for source in sources:
+    #     pols.append(p_amplitude(source.m9, gamma_co).copy())
+    # pols = pred_pol.copy()
+    # pols[pols>0] = 1
+    # pols[pols<0] = -1
+    return pols
+
+##
 
 def geo_synthetics(
         engine, targets, sources, outmode='stacked_array', plot=False,

diff --git a/beat/models/distributions.py b/beat/models/distributions.py
@@ -136,7 +136,20 @@ def multivariate_normal_chol(
                 (tt.dot(tmp, tmp))))
 
     return logpts
-
+##Mahdi
+
+def cumulative_normal(x, s=num.sqrt(2)):
+    # import scipy.special as ss
+    # Cumulative distribution function for the standard normal distribution
+    return 0.5 + 0.5 * tt.erf(x/s)
+
+def polarity_llk(observed_amplitude, pred_polarity, gamma, sigma):
+    n_t = len(observed_amplitude)
+    llks = tt.zeros((n_t), tconfig.floatX)
+    PI = gamma + (1-2*gamma)*cumulative_normal(observed_amplitude/sigma)
+    llks = tt.add(tt.mul(tt.log(PI),tt.true_div(tt.add(1,pred_polarity),2)), tt.mul(tt.log(tt.sub(1,PI)),tt.true_div(tt.sub(1,pred_polarity),2)))
+    return llks.sum()
+###
 
 def hyper_normal(datasets, hyperparams, llks, hp_specific=False):
     """

diff --git a/beat/models/problems.py b/beat/models/problems.py
@@ -49,6 +49,7 @@ def __str__(self):
 
 
 geometry_composite_catalog = {
+    'polarity': seismic.PolarityComposite,      ##Mahdi
     'seismic': seismic.SeismicGeometryComposite,
     'geodetic': geodetic.GeodeticGeometryComposite}
 

diff --git a/beat/models/seismic.py b/beat/models/seismic.py
@@ -19,7 +19,7 @@
 from beat import heart, covariance as cov
 from beat.models.base import \
     ConfigInconsistentError, Composite, FaultGeometryNotFoundError
-from beat.models.distributions import multivariate_normal_chol, get_hyper_name
+from beat.models.distributions import multivariate_normal_chol, get_hyper_name, polarity_llk
 
 from pymc3 import Uniform, Deterministic
 from collections import OrderedDict
@@ -936,7 +936,55 @@ def update_weights(
                 dataset.covariance.update_slog_pdet()
                 wmap.weights[i].set_value(choli)
 
+##Mahdi
+class SeismicPolarityComposite(Composite):
 
+    def __init__(self, polc, project_dir, sources, events, hypers=False):
+        super(SeismicPolarityComposite, self).__init__(events)
+
+        self.name = 'polarity'
+        self._like_name = 'pol_like'
+        self.synthesizers = {}
+        self.sources = sources    
+        self.config = polc
+        self.sigma = 0.16
+        self.gamma = 0.0
+        self.targets = []
+        self.data_handlers = []
+        for i in range(self.nevents):
+            seismic_data_path = os.path.join(project_dir, 
+                                             bconfig.multi_event_seismic_data_name(i))
+            self.data_handlers.append(heart.pol_datahandler(polarity_config=polc,
+                                                            seismic_data_path=seismic_data_path))
+
+    def point2sources(self, point):
+        return None
+    def export():
+        return None
+    def assemble_results():
+        return None
+    def get_formula(self, input_rvs, fixed_rvs, hyperparams, problem_config):
+
+        self.input_rvs = input_rvs
+        self.fixed_rvs = fixed_rvs
+        self.input_rvs.update(fixed_rvs)
+
+        plogpts = []
+        for i in range(len(self.data_handlers)):
+            data = self.data_handlers[i]
+            data.prepare_data(self.sources[i])
+            sources = [self.sources[i]]
+            self.synthesizers[i] = theanof.PolSynthesizer(sources, data._targets.values())
+            llk = polarity_llk(data.amplitudes, self.synthesizers[i](self.input_rvs), self.gamma, self.sigma)           ## Gamma, sigma
+            # synthetics[i] = heart.syn_polarity(self.sources[i], data._targets.values())
+            # plogpts.append(heart.polarity_llk(data.amplitudes, synthetics[i], 0.0 , 0.16))
+            plogpts.append(llk)
+        llks = Deterministic(self._like_name, tt.concatenate((plogpts)))
+        return llks.sum()
+    def get_synthetics(self, source, targets):
+
+        return None
+##
 class SeismicDistributerComposite(SeismicComposite):
     """
     Comprises how to solve the seismic (kinematic) linear forward model.

diff --git a/beat/theanof.py b/beat/theanof.py
@@ -409,7 +409,41 @@ def infer_shape(self, node, input_shapes):
         ncol = int(num.ceil(
             store.config.sample_rate * self.arrival_taper.duration()))
         return [(nrow, ncol), (nrow,)]
+##Mahdi
+class PolSynthesizer(theano.Op):
+
+    __props__ = ('sources', 'targets')
 
+    def __init__(self, sources, targets):
+        self.sources = tuple(sources)
+        self.targets = tuple(targets)
+        self.nobs = len(targets) * len(sources)
+
+    def __getstate__(self):
+        self.engine.close_cashed_stores()
+        return self.__dict__
+
+    def __setstate__(self, state):
+        self.__dict__.update(state)
+
+    def make_node(self, inputs):
+        inlist = []
+        for i in inputs.values():
+            inlist.append(tt.as_tensor_variable(i))
+
+        out = tt.as_tensor_variable(num.zeros((2,2)))
+        outlist = [out.type()]
+        return theano.Apply(self, inlist, outlist)
+
+    def perform(self, node, inputs, output):
+        synths = output[0]
+
+        synths[0] = heart.syn_polarity(self.sources, self.targets)
+
+    def infer_shape(self, node, input_shapes):
+        return [(self.nobs,2)]
+
+###
 
 class SeisDataChopper(theano.Op):
     """