def pearson4cdf(X, m, nu, a, _lambda, mu, sigma):
    # pearson4pdf
    #   p = pearson4pdf(X,m,nu,a,lambda)
    #
    #   Returns the pearson type IV probability density function with
    #   parameters m, nu, a and lambda at the values of X.
    #
    #   Example
    #
    #   See also
    #       pearson4pdf betapdf normpdf
    #       pearspdf pearsrnd
    #
    Xx = (X - _lambda) / a
    
    if Xx < -sqrt(3):    
        p1 = fx(X, m, nu, a, _lambda) * a / (2 * m - 1) * (1j - Xx) * hyp2f1(1, m + nu / 2 * 1j, 2 * m, 2 / (1 - 1j * Xx))    
        p = float(p1.real)   
    elif Xx > sqrt(3):    
        p1 = 1 - fx(-X, m, -nu, a, -_lambda) * a / (2 * m - 1) * (1j + Xx) * hyp2f1(1, m - nu / 2 * 1j, 2 * m, 2 / (1 + 1j * Xx))  
        p = float(p1.real)    
    elif Xx < 0 and Xx > -sqrt(3) and abs(nu) < (4 - 2 * sqrt(3)) * m:    
        p1 = norm.cdf(X, mu, sigma)    
        p = float(p1.real)   
    elif Xx < 0 and Xx > -sqrt(3) and abs(nu) > (4 - 2 * sqrt(3)) * m:    
        p1 = (1 - exp(-(nu + 1j * 2 * m) * pi)) ** (-1) - (1j * a * fx(X, m, nu, a, _lambda)) / (1j * nu - 2 * m + 2) * (1 + Xx ** 2) * hyp2f1(1, 2 - 2 * m, 2 - m + 1j * nu / 2, (1 + 1j * Xx) / 2)    
        p = float(p1.real)   
    elif Xx > 0 and Xx < sqrt(3) and abs(nu) < (4 - 2 * sqrt(3)) * m:    
        p1 = norm.cdf(X, mu, sigma)    
        p = float(p1.real)  
    else:     
        p1 = 1 - (1 - exp(-(-nu + 1j * 2 * m) * pi)) ** (-1) + (1j * a * fx(-X, m, -nu, a, -_lambda)) / (1j * (-nu) - 2 * m + 2) * (1 + (-Xx) ** 2) * hyp2f1(1,2-2*m,2-m-1j*nu/2,(1-1j*Xx)/2)    
        p = float(p1.real)
    
    return p

def e_ratio(a,b,e,x):
  # Get S
  bt2 = mp.beta(a,b-1.0)             # Beta function
  bix = mp.betainc(a,b+1.0,0.0,e)    # Incomplete Beta function
  hf = mp.hyp2f1(1.0,a,a+b-1.0,-1.0) # 2F1, Gauss' hypergeometric function
  hfre = mp.re(hf)
  Sval = bix - x*bt2*hfre
  # Get U
  c1 = mp.mpc(1.0 + a)
  c2 = mp.mpc(-b)
  c3 = mp.mpc(1.0)
  c4 = mp.mpc(2.0 + a)
  Uval = mp.appellf1(c1,c2,c3,c4,e,-e)
  Ure = mp.re(Uval)
  # Get P & Q
  Pval = mp.hyp2f1(a+1.0,1.0-b,a+2.0,e) # 2F1, Gauss' hypergeometric function
  Pre = mp.re(Pval)
  Qval = mp.hyp2f1(a+1.0,2.0-b,a+2.0,e) # 2F1, Gauss' hypergeometric function
  Qre = mp.re(Qval)
  # Get T
  Tval = ( (e**(1.0+a)) / (1.0+a) )*( 3.0*Pre + 2.0*Qre - Ure )
  Tval = Tval + 4.0*Sval
  # Get Rval (ratio)
  Rval = 0.25*(1.0-e*e)*( (1.0-e)**(1.0-b) )*( e**(1.0-a) )*Tval
  return Rval

def viavel(beta, c, alpha, m, n):
    if c == 0.:
        return 2*alpha*beta/(2*n*m+1)
    if c == 1.:
        return -1.
    if m == 0.:
        return 1
    k = 2*n*m-1
    q = 2*alpha*beta
    return (1-c)*float(hyp2f1(1,1,k+2,c))-q*(-k-2+(k+1)*float(hyp2f1(1,1,k+3,c)))/(k+2) - 1 

def test_hyp2f1_strange_points():
    pts = [(2, -1, -1, 0.7), (2, -2, -2, 0.7)]
    kw = dict(eliminate=True)
    dataset = [p + (float(mpmath.hyp2f1(*p, **kw)),) for p in pts]
    dataset = np.array(dataset, dtype=np.float_)

    FuncData(sc.hyp2f1, dataset, (0, 1, 2, 3), 4, rtol=1e-10).check()

def test_hyp2f1_real_some_points():
    pts = [
        (1,2,3,0),
        (1./3, 2./3, 5./6, 27./32),
        (1./4, 1./2, 3./4, 80./81),
        (2,-2,-3,3),
        (2,-3,-2,3),
        (2,-1.5,-1.5,3),
        (1,2,3,0),
        (0.7235, -1, -5, 0.3),
        (0.25, 1./3, 2, 0.999),
        (0.25, 1./3, 2, -1),
        (2,3,5,0.99),
        (3./2,-0.5,3,0.99),
        (2,2.5,-3.25,0.999),
        (-8, 18.016500331508873, 10.805295997850628, 0.90875647507000001),
        (-10,900,-10.5,0.99),
        (-10,900,10.5,0.99),
        (-1,2,1,1.0),
        (-1,2,1,-1.0),
        (-3,13,5,1.0),
        (-3,13,5,-1.0),
    ]
    dataset = [p + (float(mpmath.hyp2f1(*p)),) for p in pts]
    dataset = np.array(dataset, dtype=np.float_)

    olderr = np.seterr(invalid='ignore')
    try:
        FuncData(sc.hyp2f1, dataset, (0,1,2,3), 4, rtol=1e-10).check()
    finally:
        np.seterr(**olderr)

def _student_t_cdf(df, t, dps=None):
    if dps is None:
        dps = mpmath.mp.dps
    with mpmath.workdps(dps):
        df, t = mpmath.mpf(df), mpmath.mpf(t)
        fac = mpmath.hyp2f1(0.5, 0.5*(df + 1), 1.5, -t**2/df)
        fac *= t*mpmath.gamma(0.5*(df + 1))
        fac /= mpmath.sqrt(mpmath.pi*df)*mpmath.gamma(0.5*df)
        return 0.5 + fac

def zk(z, k):
    """
    modified dispersion function for Kappa distribution.
    (Mace and Hellberg, 1995)
    
    """
    i = mp.mpc(0, 1)
    coeff = i * (k + 0.5) * (k-0.5) / (mp.sqrt(k**3) * (k+1))
    return coeff * hyp2f1(1, 2*k+2, k+2, (1-z/(i * mp.sqrt(k)))/2)

def test_hyp2f1_real_some():
    dataset = []
    for a in [-10, -5, -1.8, 1.8, 5, 10]:
        for b in [-2.5, -1, 1, 7.4]:
            for c in [-9, -1.8, 5, 20.4]:
                for z in [-10, -1.01, -0.99, 0, 0.6, 0.95, 1.5, 10]:
                    try:
                        v = float(mpmath.hyp2f1(a, b, c, z))
                    except:
                        continue
                    dataset.append((a, b, c, z, v))
    dataset = np.array(dataset, dtype=np.float_)
    FuncData(sc.hyp2f1, dataset, (0, 1, 2, 3), 4, rtol=1e-9).check()

def test_hyp2f1_some_points_2():
    # Taken from mpmath unit tests -- this point failed for mpmath 0.13 but
    # was fixed in their SVN since then
    pts = [(112, (51, 10), (-9, 10), -0.99999), (10, -900, 10.5, 0.99), (10, -900, -10.5, 0.99)]

    def fev(x):
        if isinstance(x, tuple):
            return float(x[0]) / x[1]
        else:
            return x

    dataset = [tuple(map(fev, p)) + (float(mpmath.hyp2f1(*p)),) for p in pts]
    dataset = np.array(dataset, dtype=np.float_)

    FuncData(sc.hyp2f1, dataset, (0, 1, 2, 3), 4, rtol=1e-10).check()

def test_hyp2f1_real_some():
    dataset = []
    for a in [-10, -5, -1.8, 1.8, 5, 10]:
        for b in [-2.5, -1, 1, 7.4]:
            for c in [-9, -1.8, 5, 20.4]:
                for z in [-10, -1.01, -0.99, 0, 0.6, 0.95, 1.5, 10]:
                    try:
                        v = float(mpmath.hyp2f1(a, b, c, z))
                    except:
                        continue
                    dataset.append((a, b, c, z, v))
    dataset = np.array(dataset, dtype=np.float_)

    olderr = np.seterr(invalid='ignore')
    try:
        FuncData(sc.hyp2f1, dataset, (0,1,2,3), 4, rtol=1e-9,
                 ignore_inf_sign=True).check()
    finally:
        np.seterr(**olderr)

def test_hyp2f1_real_some_points():
    pts = [
        (1, 2, 3, 0),
        (1.0 / 3, 2.0 / 3, 5.0 / 6, 27.0 / 32),
        (1.0 / 4, 1.0 / 2, 3.0 / 4, 80.0 / 81),
        (2, -2, -3, 3),
        (2, -3, -2, 3),
        (2, -1.5, -1.5, 3),
        (1, 2, 3, 0),
        (0.7235, -1, -5, 0.3),
        (0.25, 1.0 / 3, 2, 0.999),
        (0.25, 1.0 / 3, 2, -1),
        (2, 3, 5, 0.99),
        (3.0 / 2, -0.5, 3, 0.99),
        (2, 2.5, -3.25, 0.999),
        (-8, 18.016500331508873, 10.805295997850628, 0.90875647507000001),
        (-10, 900, -10.5, 0.99),
        (-10, 900, 10.5, 0.99),
    ]
    dataset = [p + (float(mpmath.hyp2f1(*p)),) for p in pts]
    dataset = np.array(dataset, dtype=np.float_)

    FuncData(sc.hyp2f1, dataset, (0, 1, 2, 3), 4, rtol=1e-10).check()

def test_hyp2f1_real_random():
    dataset = []

    npoints = 500
    dataset = np.zeros((npoints, 5), np.float_)

    np.random.seed(1234)
    dataset[:,0] = np.random.pareto(1.5, npoints)
    dataset[:,1] = np.random.pareto(1.5, npoints)
    dataset[:,2] = np.random.pareto(1.5, npoints)
    dataset[:,3] = 2*np.random.rand(npoints) - 1

    dataset[:,0] *= (-1)**np.random.randint(2, npoints)
    dataset[:,1] *= (-1)**np.random.randint(2, npoints)
    dataset[:,2] *= (-1)**np.random.randint(2, npoints)

    for ds in dataset:
        if mpmath.__version__ < '0.14':
            # mpmath < 0.14 fails for c too much smaller than a, b
            if abs(ds[:2]).max() > abs(ds[2]):
                ds[2] = abs(ds[:2]).max()
        ds[4] = float(mpmath.hyp2f1(*tuple(ds[:4])))

    FuncData(sc.hyp2f1, dataset, (0,1,2,3), 4, rtol=1e-9).check()

def z(ZZ):
    mpmath.mp.dps = 20
    return N.exp(-N.log(1728*ZZ) / 5) * mpmath.hyp2f1(31.0/60, 11.0/60, 6.0/5, 1.0/ZZ) /\
                                        mpmath.hyp2f1(19.0/60, -1.0/60, 4.0/5, 1.0/ZZ)

as yet, since it tries to import the libraries from *this* folder, rather than
installation (which doesn't work because the fortran code isn't installed.)
"""
import inspect
import os

#LOCATION = "/".join(os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe()))).split("/")[:-1])
import sys
#sys.path.insert(0, LOCATION)
import numpy as np
from halomod import ProjectedCF
from halomod.integrate_corr import projected_corr_gal
from astropy.units import Mpc
from mpmath import gamma,hyp2f1

hyp2f1A = np.frompyfunc(lambda a,b,c,z: float(hyp2f1(a,b,c,z)), 4, 1)

def wprp_pl(rp,r0,g):
    return rp*(rp/r0)**-g * gamma(0.5)*gamma((g-1)/2.0)/gamma(g/2.0)

def wprp_pl_lim(rp,r0,g,rmax):
    return (1/gamma(g/2)) * (h.rp*h.rlim/r0)**-g * gamma((g-1)/2) * \
           (gamma(0.5) * h.rp * h.rlim**g - h.rp**g *h.rlim * gamma(g/2)*\
           hyp2f1A(0.5,(g-1)/2,(g+1)/2,h.rp.value**2/h.rlim.value**2))

class TestProjCorr():
    def __init__(self):
        self.rp = np.logspace(-2,1.2,50)
        self.gamma = 1.85 # Values from S1 sample of Beutler+2011
        self.r0 = 5.14

 def test_hyp2f1(self):
     assert_mpmath_equal(sc.hyp2f1,
                         _exception_to_nan(lambda a, b, c, x: mpmath.hyp2f1(a, b, c, x, **HYPERKW)),
                         [Arg(), Arg(), Arg(), Arg()])

 def test_hyp2f1_complex(self):
     # Scipy's hyp2f1 seems to have performance and accuracy problems
     assert_mpmath_equal(lambda a, b, c, x: sc.hyp2f1(a.real, b.real, c.real, x),
                         _exception_to_nan(lambda a, b, c, x: mpmath.hyp2f1(a, b, c, x, **HYPERKW)),
                         [Arg(-1e2, 1e2), Arg(-1e2, 1e2), Arg(-1e2, 1e2), ComplexArg()],
                         n=10)