1# Copyright (c) 2012, 2013 Ricardo Andrade
2# Copyright (c) 2015 James Hensman
3# Licensed under the BSD 3-clause license (see LICENSE.txt)
4
5import numpy as np
6from scipy.special import ndtr as std_norm_cdf
7
8#define a standard normal pdf
9_sqrt_2pi = np.sqrt(2*np.pi)
10def std_norm_pdf(x):
11    x = np.clip(x,-1e150,1e150)
12    return np.exp(-np.square(x)/2)/_sqrt_2pi
13
14def inv_std_norm_cdf(x):
15    """
16    Inverse cumulative standard Gaussian distribution
17    Based on Winitzki, S. (2008)
18    """
19    z = 2*x -1
20    ln1z2 = np.log(1-z**2)
21    a = 8*(np.pi -3)/(3*np.pi*(4-np.pi))
22    b = 2/(np.pi * a) + ln1z2/2
23    inv_erf = np.sign(z) * np.sqrt( np.sqrt(b**2 - ln1z2/a) - b )
24    return np.sqrt(2) * inv_erf
25
26def logPdfNormal(z):
27    """
28    Robust implementations of log pdf of a standard normal.
29
30     @see [[https://github.com/mseeger/apbsint/blob/master/src/eptools/potentials/SpecfunServices.h original implementation]]
31     in C from Matthias Seeger.
32    """
33    return -0.5 * (M_LN2PI + z * z)
34
35def cdfNormal(z):
36    """
37    Robust implementations of cdf of a standard normal.
38
39     @see [[https://github.com/mseeger/apbsint/blob/master/src/eptools/potentials/SpecfunServices.h original implementation]]
40     in C from Matthias Seeger.
41    */
42    """
43    if (abs(z) < ERF_CODY_LIMIT1):
44        # Phi(z) approx (1+y R_3(y^2))/2, y=z/sqrt(2)
45        return 0.5 * (1.0 + (z / M_SQRT2) * _erfRationalHelperR3(0.5 * z * z))
46    elif (z < 0.0):
47        # Phi(z) approx N(z)Q(-z)/(-z), z<0
48        return np.exp(logPdfNormal(z)) * _erfRationalHelper(-z) / (-z)
49    else:
50        return 1.0 - np.exp(logPdfNormal(z)) * _erfRationalHelper(z) / z
51
52
53
54def logCdfNormal(z):
55    """
56    Robust implementations of log cdf of a standard normal.
57
58     @see [[https://github.com/mseeger/apbsint/blob/master/src/eptools/potentials/SpecfunServices.h original implementation]]
59     in C from Matthias Seeger.
60    """
61    if (abs(z) < ERF_CODY_LIMIT1):
62        # Phi(z) approx  (1+y R_3(y^2))/2, y=z/sqrt(2)
63        return np.log1p((z / M_SQRT2) * _erfRationalHelperR3(0.5 * z * z)) - M_LN2
64    elif (z < 0.0):
65        # Phi(z) approx N(z)Q(-z)/(-z), z<0
66        return logPdfNormal(z) - np.log(-z) + np.log(_erfRationalHelper(-z))
67    else:
68        return np.log1p(-(np.exp(logPdfNormal(z))) * _erfRationalHelper(z) / z)
69
70
71
72def derivLogCdfNormal(z):
73    """
74     Robust implementations of derivative of the log cdf of a standard normal.
75
76     @see [[https://github.com/mseeger/apbsint/blob/master/src/eptools/potentials/SpecfunServices.h original implementation]]
77     in C from Matthias Seeger.
78    """
79    if (abs(z) < ERF_CODY_LIMIT1):
80        # Phi(z) approx (1 + y R_3(y^2))/2, y = z/sqrt(2)
81        return 2.0 * np.exp(logPdfNormal(z)) / (1.0 + (z / M_SQRT2) * _erfRationalHelperR3(0.5 * z * z))
82    elif (z < 0.0):
83        # Phi(z) approx N(z) Q(-z)/(-z), z<0
84        return -z / _erfRationalHelper(-z)
85    else:
86        t = np.exp(logPdfNormal(z))
87        return t / (1.0 - t * _erfRationalHelper(z) / z)
88
89
90def _erfRationalHelper(x):
91    assert x > 0.0, "Arg of erfRationalHelper should be >0.0; was {}".format(x)
92
93    if (x >= ERF_CODY_LIMIT2):
94        """
95         x/sqrt(2) >= 4
96
97         Q(x)   = 1 + sqrt(pi) y R_1(y),
98         R_1(y) = poly(p_j,y) / poly(q_j,y),  where  y = 2/(x*x)
99
100         Ordering of arrays: 4,3,2,1,0,5 (only for numerator p_j; q_5=1)
101         ATTENTION: The p_j are negative of the entries here
102         p (see P1_ERF)
103         q (see Q1_ERF)
104        """
105        y = 2.0 / (x * x)
106
107        res = y * P1_ERF[5]
108        den = y
109        i = 0
110
111        while (i <= 3):
112            res = (res + P1_ERF[i]) * y
113            den = (den + Q1_ERF[i]) * y
114            i += 1
115
116        # Minus, because p(j) values have to be negated
117        return 1.0 - M_SQRTPI * y * (res + P1_ERF[4]) / (den + Q1_ERF[4])
118    else:
119        """
120         x/sqrt(2) < 4, x/sqrt(2) >= 0.469
121
122         Q(x)   = sqrt(pi) y R_2(y),
123         R_2(y) = poly(p_j,y) / poly(q_j,y),   y = x/sqrt(2)
124
125         Ordering of arrays: 7,6,5,4,3,2,1,0,8 (only p_8; q_8=1)
126         p (see P2_ERF)
127         q (see Q2_ERF
128        """
129        y = x / M_SQRT2
130        res = y * P2_ERF[8]
131        den = y
132        i = 0
133
134        while (i <= 6):
135            res = (res + P2_ERF[i]) * y
136            den = (den + Q2_ERF[i]) * y
137            i += 1
138
139        return M_SQRTPI * y * (res + P2_ERF[7]) / (den + Q2_ERF[7])
140
141def _erfRationalHelperR3(y):
142    assert y >= 0.0, "Arg of erfRationalHelperR3 should be >=0.0; was {}".format(y)
143
144    nom = y * P3_ERF[4]
145    den = y
146    i = 0
147    while (i <= 2):
148        nom = (nom + P3_ERF[i]) * y
149        den = (den + Q3_ERF[i]) * y
150        i += 1
151    return (nom + P3_ERF[3]) / (den + Q3_ERF[3])
152
153ERF_CODY_LIMIT1 = 0.6629
154ERF_CODY_LIMIT2 = 5.6569
155M_LN2PI         = 1.83787706640934533908193770913
156M_LN2           = 0.69314718055994530941723212146
157M_SQRTPI        = 1.77245385090551602729816748334
158M_SQRT2         = 1.41421356237309504880168872421
159
160#weights for the erfHelpers (defined here to avoid redefinitions at every call)
161P1_ERF = [
1623.05326634961232344e-1, 3.60344899949804439e-1,
1631.25781726111229246e-1, 1.60837851487422766e-2,
1646.58749161529837803e-4, 1.63153871373020978e-2]
165Q1_ERF = [
1662.56852019228982242e+0, 1.87295284992346047e+0,
1675.27905102951428412e-1, 6.05183413124413191e-2,
1682.33520497626869185e-3]
169P2_ERF = [
1705.64188496988670089e-1, 8.88314979438837594e+0,
1716.61191906371416295e+1, 2.98635138197400131e+2,
1728.81952221241769090e+2, 1.71204761263407058e+3,
1732.05107837782607147e+3, 1.23033935479799725e+3,
1742.15311535474403846e-8]
175Q2_ERF = [
1761.57449261107098347e+1, 1.17693950891312499e+2,
1775.37181101862009858e+2, 1.62138957456669019e+3,
1783.29079923573345963e+3, 4.36261909014324716e+3,
1793.43936767414372164e+3, 1.23033935480374942e+3]
180P3_ERF = [
1813.16112374387056560e+0, 1.13864154151050156e+2,
1823.77485237685302021e+2, 3.20937758913846947e+3,
1831.85777706184603153e-1]
184Q3_ERF = [
1852.36012909523441209e+1, 2.44024637934444173e+2,
1861.28261652607737228e+3, 2.84423683343917062e+3]
187