# Copyright (c) 2013, GPy authors (see AUTHORS.txt).
# Licensed under the BSD 3-clause license (see LICENSE.txt)
#Only works for scipy 0.12+
from scipy.special import erfcx, erf
import numpy as np
[docs]def ln_diff_erfs(x1, x2, return_sign=False):
"""Function for stably computing the log of difference of two erfs in a numerically stable manner.
:param x1 : argument of the positive erf
:type x1: ndarray
:param x2 : argument of the negative erf
:type x2: ndarray
:return: tuple containing (log(abs(erf(x1) - erf(x2))), sign(erf(x1) - erf(x2)))
Based on MATLAB code that was written by Antti Honkela and modified by David Luengo and originally derived from code by Neil Lawrence.
"""
x1 = np.require(x1).real
x2 = np.require(x2).real
if x1.size==1:
x1 = np.reshape(x1, (1, 1))
if x2.size==1:
x2 = np.reshape(x2, (1, 1))
if x1.shape==x2.shape:
v = np.zeros_like(x1)
else:
if x1.size==1:
v = np.zeros(x2.shape)
elif x2.size==1:
v = np.zeros(x1.shape)
else:
raise ValueError("This function does not broadcast unless provided with a scalar.")
if x1.size == 1:
x1 = np.tile(x1, x2.shape)
if x2.size == 1:
x2 = np.tile(x2, x1.shape)
sign = np.sign(x1 - x2)
if x1.size == 1:
if sign== -1:
swap = x1
x1 = x2
x2 = swap
else:
I = sign == -1
swap = x1[I]
x1[I] = x2[I]
x2[I] = swap
with np.errstate(divide='ignore'):
# switch off log of zero warnings.
# Case 0: arguments of different sign, no problems with loss of accuracy
I0 = np.logical_or(np.logical_and(x1>0, x2<0), np.logical_and(x2>0, x1<0)) # I1=(x1*x2)<0
# Case 1: x1 = x2 so we have log of zero.
I1 = (x1 == x2)
# Case 2: Both arguments are non-negative
I2 = np.logical_and(x1 > 0, np.logical_and(np.logical_not(I0),
np.logical_not(I1)))
# Case 3: Both arguments are non-positive
I3 = np.logical_and(np.logical_and(np.logical_not(I0),
np.logical_not(I1)),
np.logical_not(I2))
_x2 = x2.flatten()
_x1 = x1.flatten()
for group, flags in zip((0, 1, 2, 3), (I0, I1, I2, I3)):
if np.any(flags):
if not x1.size==1:
_x1 = x1[flags]
if not x2.size==1:
_x2 = x2[flags]
if group==0:
v[flags] = np.log( erf(_x1) - erf(_x2) )
elif group==1:
v[flags] = -np.inf
elif group==2:
v[flags] = np.log(erfcx(_x2)
-erfcx(_x1)*np.exp(_x2**2
-_x1**2)) - _x2**2
elif group==3:
v[flags] = np.log(erfcx(-_x1)
-erfcx(-_x2)*np.exp(_x1**2
-_x2**2))-_x1**2
# TODO: switch back on log of zero warnings.
if return_sign:
return v, sign
else:
if v.size==1:
if sign==-1:
v = v.view('complex64')
v += np.pi*1j
else:
# Need to add in a complex part because argument is negative.
v = v.view('complex64')
v[I] += np.pi*1j
return v