a horrible expression for the integral to compute potential

        :param x: R/Rs
        :param tau: t/Rs
        :type x: float >0
171276
def _h(self, X, tau):

        

        def cos_func(y):
            if isinstance(y, float) or isinstance(y, int):
                if y > 1:
                    return np.arccosh(y)
                else:
                    return np.arccos(y)
            else:
                values = np.ones_like(y)
                inds1 = np.where(y < 1)
                inds2 = np.where(y > 1)
                values[inds1] = np.arccos(y[inds1])
                # values[inds2] = np.arccosh(y[inds2])
                values[inds2] = np.arccos((1 - y[inds2]) ** .5)
                return values

        def nfw_func(x):
            if isinstance(x,float) or isinstance(x,int):
                if x<1:
                    return np.log(x / 2.) ** 2 - np.arccosh(1. / x) ** 2
                else:
                    return np.log(x / 2.)**2 - np.arccos(1. /x)**2
            else:
                inds1 = np.where(x < 1)
                inds2 = np.where(x >= 1)
                y = np.zeros_like(x)
                y[inds1] = np.log(x[inds1] / 2.) ** 2 - np.arccosh(1. / x[inds1]) ** 2
                y[inds2] = np.log(x[inds2] / 2.) ** 2 - np.arccosh(1. / x[inds2]) ** 2

                return y

        def tnfw_func(x,tau):

            u = x**2
            t2 = tau**2
            Lx = self.L(u**.5, tau)

            return (t2 + 1) ** -2 * (
                2 * t2 * np.pi * (tau - (t2 + u) ** .5 + tau * np.log(tau + (t2 + u) ** .5))
                +
                2 * (t2 - 1) * tau * (t2 + u) ** .5 * Lx
                +
                t2 * (t2 - 1) * Lx ** 2
                +
                4 * t2 * (u - 1) * self.F(u**.5)
                +
                t2 * (t2 - 1) * (cos_func(u ** -0.5)) ** 2
                +
                t2 * ((t2 - 1) * np.log(tau) - t2 - 1) * np.log(u)
                -
                t2 * (
                (t2 - 1) * np.log(tau) * np.log(4 * tau) + 2 * np.log(0.5 * tau) - 2 * tau * (tau - np.pi) * np.log(
                    tau * 2)))

        c = 1e-9
        rescale = 1

        if np.any(X-c<1):

            warnings.warn('Truncated NFW potential not yet implemented for x<1. Using the expression for the NFW '
                          'potential in this regime isntead.')
            rescale = tnfw_func(1.00001,tau)*nfw_func(1)**-1

        if isinstance(X,float) or isinstance(X,int):
            if X<1:
                X = min(X,0.001)
                return nfw_func(X)*rescale
            else:
                return tnfw_func(X,tau)
        else:

            X[np.where(X<0.001)] = 0.001

            values = np.zeros_like(X)
            inds1,inds2 = np.where(X<1),np.where(X>=1)
            values[inds1] = nfw_func(X[inds1])*rescale
            values[inds2] = tnfw_func(X[inds2],tau)

            return values
31101
def cumulative_distribution(self, X, U=0):
        
        self.check_fit()

        low_bounds = self.model.dataset.mean() - (5 * self.model.dataset.std())
        return self.model.integrate_box_1d(low_bounds, X) - U
1181830
def prox_hard(X, step, thresh=0):
    
    thresh_ = _step_gamma(step, thresh)
    below = np.abs(X) < thresh_
    X[below] = 0
    return X
1285746
def sum(x, weights=None):
    '''
    
    '''
    x = to_potential(x)
    if is_const_potential(x): return PotentialConstant(np.sum(x.c))
    else:                     return PotentialSum(x, weights=weights)
1087219
def rlogistic(mu, tau, size=None):
    

    u = np.random.random(size)
    return mu + np.log(u / (1 - u)) / tau
250766
def rhalf_normal(tau, size=None):
    

    return abs(np.random.normal(0, np.sqrt(1 / tau), size))
1310726
def Exponential(x, a, tau, y0):
    
    return np.exp(x / tau) * a + y0
983748
def rmse_2pi(x, xhat):
    ''' '''
    real_diff = from_0_to_pi(x - xhat)
    np.square(real_diff, out=real_diff)
    sum_ = np.sum(real_diff)
    return sqrt(sum_ / len(x))
510007
function(x, mean, std) {
        return 0.5 * (1 + erf((x - mean) / Math.sqrt(2 * std * std)));
    }
212602
def xstep(self):
        
        """

        self.YU[:] = self.Y - self.U

        b = self.DSf + self.rho*sl.rfftn(self.YU, None, self.cri.axisN)
        if self.cri.Cd == 1:
            self.Xf[:] = sl.solvedbi_sm(self.Df, self.mu + self.rho,
                                        b, self.c, self.cri.axisM)
        else:
            self.Xf[:] = sl.solvemdbi_ism(self.Df, self.mu + self.rho, b,
                                          self.cri.axisM, self.cri.axisC)

        self.X = sl.irfftn(self.Xf, self.cri.Nv, self.cri.axisN)

        if self.opt['LinSolveCheck']:
            Dop = lambda x: sl.inner(self.Df, x, axis=self.cri.axisM)
            if self.cri.Cd == 1:
                DHop = lambda x: np.conj(self.Df) * x
            else:
                DHop = lambda x: sl.inner(np.conj(self.Df), x,
                                          axis=self.cri.axisC)
            ax = DHop(Dop(self.Xf)) + (self.mu + self.rho)*self.Xf
            self.xrrs = sl.rrs(ax, b)
        else:
            self.xrrs = None
109953
def _at_percentile(tau, var_tau, percentile):
    
    assert (percentile >= 0.0) and (percentile <= 1.0)
    c_val = _scaling(tau, var_tau)
    k_val = _dof(tau, var_tau)
    return np.sqrt(c_val * chi2.ppf(percentile, df=k_val))
1224812
def fprop(self, x):
    
    logits = tf.py_func(self.fprop_np, [x], tf.float32)
    return {self.O_LOGITS: logits}
555525
def prox_xline(x, step):
    """"""
    if not np.isscalar(x):
        x= x[0]
    if x > 0.5:
        return np.array([0.5])
    else:
        return np.array([x])
386186
def gelu(x):
  
  cdf = 0.5 * (1.0 + tf.tanh(
      (np.sqrt(2 / np.pi) * (x + 0.044715 * tf.pow(x, 3)))))
  return x * cdf
740835
def Itax(mt, x, t):
    
    return (Sx(mt, x + 1) - Sx(mt, x + t + 1)) / Dx(mt, x)
319652
def _std(self,x):
        
        return np.nanstd(x.values,ddof=self._ddof)
659427
def atanh(x):
    
    if isinstance(x, UncertainFunction):
        mcpts = np.arctanh(x._mcpts)
        return UncertainFunction(mcpts)
    else:
        return np.arctanh(x)