From 5170a0b1c6a9fa7fd9b8d5a9ead80edfcf733dca Mon Sep 17 00:00:00 2001
From: Sebastian Wehling-Benatelli <sebastian.wehling@rub.de>
Date: Mon, 1 Feb 2016 06:13:17 +0100
Subject: [PATCH] [new] implementation of a probability density function
 representation of the pick (untested)

---
 pylot/core/pick/utils.py | 111 +++++++++++++++++++++++++++++++++++++++
 1 file changed, 111 insertions(+)

diff --git a/pylot/core/pick/utils.py b/pylot/core/pick/utils.py
index d7584bff..5fa19b4c 100644
--- a/pylot/core/pick/utils.py
+++ b/pylot/core/pick/utils.py
@@ -133,6 +133,117 @@ def earllatepicker(X, nfac, TSNR, Pick1, iplot=None, stealthMode = False):
     return EPick, LPick, PickError
 
 
+def gauss_parameter(te, tm, tl, eta):
+    '''
+    takes three onset times and returns the parameters sig1, sig2, a1 and a2
+    to represent the pick as a probability density funtion (PDF) with two
+    Gauss branches
+    :param te:
+    :param tm:
+    :param tl:
+    :param eta:
+    :return:
+    '''
+
+    sig1 = (tm - te) / np.sqrt(2 * np.log(1 / eta))
+    sig2 = (tl - tm) / np.sqrt(2 * np.log(1 / eta))
+
+    a1 = 2 / (1 + sig2 / sig1)
+    a2 = 2 / (1 + sig1 / sig2)
+
+    return sig1, sig2, a1, a2
+
+
+def exp_parameter(te, tm, tl, eta):
+    '''
+    takes three onset times te, tm and tl and returns the parameters sig1,
+    sig2 and a to represent the pick as a probability density function (PDF)
+    with two exponential decay branches
+    :param te:
+    :param tm:
+    :param tl:
+    :param eta:
+    :return:
+    '''
+
+    sig1 = np.log(eta) / (te - tm)
+    sig2 = np.log(eta) / (tm - tl)
+    a = 1 / (1 / sig1 + 1 / sig2)
+
+    return sig1, sig2, a
+
+
+def gauss_branches(x, mu, sig1, sig2, a1, a2):
+    '''
+    function gauss_branches takes an axes x, a center value mu, two sigma
+    values sig1 and sig2 and two scaling factors a1 and a2 and return a
+    list containing the values of a probability density function (PDF)
+    consisting of gauss branches
+    :param x:
+    :type x:
+    :param mu:
+    :type mu:
+    :param sig1:
+    :type sig1:
+    :param sig2:
+    :type sig2:
+    :param a1:
+    :type a1:
+    :param a2:
+    :returns fun_vals: list with function values along axes x
+    '''
+    fun_vals = []
+    for k in x:
+        if k < mu:
+            fun_vals.append(a1 * 1 / (np.sqrt(2 * np.pi) * sig1) * np.exp(-((k - mu) / sig1)**2 / 2 ))
+        else:
+            fun_vals.append(a2 * 1 / (np.sqrt(2 * np.pi) * sig2) * np.exp(-((k - mu) / sig2)**2 / 2))
+    return fun_vals
+
+
+def exp_branches(x, mu, sig1, sig2, a):
+    '''
+    function exp_branches takes an axes x, a center value mu, two sigma
+    values sig1 and sig2 and a scaling factor a and return a
+    list containing the values of a probability density function (PDF)
+    consisting of exponential decay branches
+    :param x:
+    :param mu:
+    :param sig1:
+    :param sig2:
+    :param a:
+    :returns fun_vals: list with function values along axes x:
+    '''
+    fun_vals = []
+    for k in x:
+        if k < mu:
+            fun_vals.append(a * np.exp(sig1 * (k - mu)))
+        else:
+            fun_vals.append(a * np.exp(-sig2 * (k - mu)))
+    return fun_vals
+
+
+def pick_pdf(t, te, tm, tl, type='gauss', eta=0.01):
+    '''
+    
+    :param t:
+    :param te:
+    :param tm:
+    :param tl:
+    :param type:
+    :param eta:
+    :param args:
+    :return:
+    '''
+
+    parameter = dict(gauss=gauss_parameter, exp=exp_parameter)
+    branches = dict(gauss=gauss_branches, exp=exp_branches)
+
+    params = parameter[type](te, tm, tl, eta)
+
+    return branches[type](t, tm, *params)
+
+
 def fmpicker(Xraw, Xfilt, pickwin, Pick, iplot=None):
     '''
     Function to derive first motion (polarity) of given phase onset Pick.