def ln_diff_erfs(x1, x2, return_sign=False):
x1 = np.require(x1).real<EOL>x2 = np.require(x2).real<EOL>if x1.size==<NUM_LIT:1>:<EOL><INDENT>x1 = np.reshape(x1, (<NUM_LIT:1>, <NUM_LIT:1>))<EOL><DEDENT>if x2.size==<NUM_LIT:1>:<EOL><INDENT>x2 = np.reshape(x2, (<NUM_LIT:1>, <NUM_LIT:1>))<EOL><DEDENT>if x1.shape==x2.shape:<EOL><INDENT>v = np.zeros_like(x1)<EOL><DEDENT>else:<EOL><INDENT>if x1.size==<NUM_LIT:1>:<EOL><INDENT>v = np.zeros(x2.shape)<EOL><DEDENT>elif x2.size==<NUM_LIT:1>:<EOL><INDENT>v = np.zeros(x1.shape)<EOL><DEDENT>else:<EOL><INDENT>raise ValueError("<STR_LIT>")<EOL><DEDENT><DEDENT>if x1.size == <NUM_LIT:1>:<EOL><INDENT>x1 = np.tile(x1, x2.shape)<EOL><DEDENT>if x2.size == <NUM_LIT:1>:<EOL><INDENT>x2 = np.tile(x2, x1.shape)<EOL><DEDENT>sign = np.sign(x1 - x2)<EOL>if x1.size == <NUM_LIT:1>:<EOL><INDENT>if sign== -<NUM_LIT:1>:<EOL><INDENT>swap = x1<EOL>x1 = x2<EOL>x2 = swap<EOL><DEDENT><DEDENT>else:<EOL><INDENT>I = sign == -<NUM_LIT:1><EOL>swap = x1[I]<EOL>x1[I] = x2[I]<EOL>x2[I] = swap<EOL><DEDENT>with np.errstate(divide='<STR_LIT:ignore>'):<EOL><INDENT>I0 = np.logical_or(np.logical_and(x1><NUM_LIT:0>, x2<<NUM_LIT:0>), np.logical_and(x2><NUM_LIT:0>, x1<<NUM_LIT:0>)) <EOL>I1 = (x1 == x2)<EOL>I2 = np.logical_and(x1 > <NUM_LIT:0>, np.logical_and(np.logical_not(I0),<EOL>np.logical_not(I1)))<EOL>I3 = np.logical_and(np.logical_and(np.logical_not(I0),<EOL>np.logical_not(I1)),<EOL>np.logical_not(I2))<EOL>_x2 = x2.flatten()<EOL>_x1 = x1.flatten()<EOL>for group, flags in zip((<NUM_LIT:0>, <NUM_LIT:1>, <NUM_LIT:2>, <NUM_LIT:3>), (I0, I1, I2, I3)):<EOL><INDENT>if np.any(flags):<EOL><INDENT>if not x1.size==<NUM_LIT:1>:<EOL><INDENT>_x1 = x1[flags]<EOL><DEDENT>if not x2.size==<NUM_LIT:1>:<EOL><INDENT>_x2 = x2[flags]<EOL><DEDENT>if group==<NUM_LIT:0>:<EOL><INDENT>v[flags] = np.log( erf(_x1) - erf(_x2) )<EOL><DEDENT>elif group==<NUM_LIT:1>:<EOL><INDENT>v[flags] = -np.inf<EOL><DEDENT>elif group==<NUM_LIT:2>:<EOL><INDENT>v[flags] = np.log(erfcx(_x2)<EOL>-erfcx(_x1)*np.exp(_x2**<NUM_LIT:2><EOL>-_x1**<NUM_LIT:2>)) - _x2**<NUM_LIT:2><EOL><DEDENT>elif group==<NUM_LIT:3>:<EOL><INDENT>v[flags] = np.log(erfcx(-_x1)<EOL>-erfcx(-_x2)*np.exp(_x1**<NUM_LIT:2><EOL>-_x2**<NUM_LIT:2>))-_x1**<NUM_LIT:2><EOL><DEDENT><DEDENT><DEDENT><DEDENT>if return_sign:<EOL><INDENT>return v, sign<EOL><DEDENT>else:<EOL><INDENT>if v.size==<NUM_LIT:1>:<EOL><INDENT>if sign==-<NUM_LIT:1>:<EOL><INDENT>v = v.view('<STR_LIT>')<EOL>v += np.pi*<NUM_LIT><EOL><DEDENT><DEDENT>else:<EOL><INDENT>v = v.view('<STR_LIT>')<EOL>v[I] += np.pi*<NUM_LIT><EOL><DEDENT><DEDENT>return v<EOL>
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.