def _surfdens(self,R,z,phi=0.,t=0.):
     """
     NAME:
        _surfdens
     PURPOSE:
        evaluate the surface density for this potential
     INPUT:
        R - Galactocentric cylindrical radius
        z - vertical height
        phi - azimuth
        t - time
     OUTPUT:
        the surface density
     HISTORY:
        2018-08-19 - Written - Bovy (UofT)
     """
     r= numpy.sqrt(R**2.+z**2.)
     x= r/self.a
     Rpa= numpy.sqrt(R**2.+self.a**2.)
     Rma= numpy.sqrt(R**2.-self.a**2.+0j)
     if Rma == 0:
         za= z/self.a
         return self.a**2./2.*((2.-2.*numpy.sqrt(za**2.+1)
                                +numpy.sqrt(2.)*za\
                                    *numpy.arctan(za/numpy.sqrt(2.)))/z
                               +numpy.sqrt(2*za**2.+2.)\
                                *numpy.arctanh(za/numpy.sqrt(2.*(za**2.+1)))
                               /numpy.sqrt(self.a**2.+z**2.))
     else:
         return self.a**2.*(numpy.arctan(z/x/Rma)/Rma
                            +numpy.arctanh(z/x/Rpa)/Rpa
                            -numpy.arctan(z/Rma)/Rma
                            +numpy.arctan(z/Rpa)/Rpa).real

    def deflections(self, xin, yin):
        from numpy import arctanh, arctan, arctan2, log, sin, cos

        #x,y = self.align_coords(xin,yin)
        x = xin - self.x
        y = yin - self.y
        b, rs = self.b, self.rs
        X = b / rs
        if X < 1.:
            amp = X**2 / (8 * arctanh(
                ((1 - X) / (1 + X))**0.5) / (1 - X**2)**0.5 + 4 * log(X / 2.))
        elif X == 1:
            amp = 0.25 / (1. + log(0.5))
        else:
            amp = X**2 / (8 * arctan(
                ((X - 1) / (1 + X))**0.5) / (X**2 - 1)**0.5 + 4 * log(X / 2.))

        r2 = (x**2 + y**2) / rs**2
        r = r2**0.5
        F = r * 0.
        F[r < 1.] = arctanh((1 - r2[r < 1.])**0.5) / (1 - r2[r < 1.])**0.5
        F[r == 1.] = 1.
        F[r > 1.] = arctan((r2[r > 1.] - 1.)**0.5) / (r2[r > 1.] - 1)**0.5

        dr = 4 * amp * rs * (log(r / 2) + F) / r
        A = arctan2(y, x)
        return dr * cos(A), dr * sin(A)

def slRSA_m_1Ss(ds, model, omit, partial_dsm = None, radius=3, cmetric='pearson'):
    '''one subject

    Executes slRSA on single subjects and returns tuple of arrays of 1-p's [0], and fisher Z transformed r's [1]

    ds: pymvpa dsets for 1 subj
    model: model DSM to be correlated with neural DSMs per searchlight center
    partial_dsm: model DSM to be partialled out of model-neural DSM correlation
    omit: list of targets omitted from pymvpa datasets
    radius: sl radius, default 3
    cmetric: default pearson, other optin 'spearman'
    '''        

    if __debug__:
        debug.active += ["SLC"]

    for om in omit:
        ds = ds[ds.sa.targets != om] # cut out omits
        print('Target |%s| omitted from analysis' % (om))
    ds = mean_group_sample(['targets'])(ds) #make UT ds
    print('Mean group sample computed at size:',ds.shape,'...with UT:',ds.UT)

    print('Beginning slRSA analysis...')
    if partial_dsm == None: tdcm = rsa.TargetDissimilarityCorrelationMeasure(squareform(model), comparison_metric=cmetric)
    elif partial_dsm != None: tdcm = rsa.TargetDissimilarityCorrelationMeasure(squareform(model), comparison_metric=cmetric, partial_dsm = squareform(partial_dsm))
    sl = sphere_searchlight(tdcm,radius=radius)
    slmap = sl(ds)
    if partial_dsm == None:
        print('slRSA complete with map of shape:',slmap.shape,'...p max/min:',slmap.samples[0].max(),slmap.samples[0].min(),'...r max/min',slmap.samples[1].max(),slmap.samples[1].min())
        return 1-slmap.samples[1],np.arctanh(slmap.samples[0])
    else:
        print('slRSA complete with map of shape:',slmap.shape,'...r max/min:',slmap.samples[0].max(),slmap.samples[0].min())
        return 1-slmap.samples[1],np.arctanh(slmap.samples[0])

def getAlphaNIEXS(lensParam,X):
	einRad = lensParam[0]
	f = lensParam[1]
	fp = np.sqrt(1 - f*f)
	bc = lensParam[2]
	t = lensParam[3]
	xshear = (lensParam[4] + 1j*lensParam[5])
	
	#### The rotated coordinates	
	Xr = X * np.exp(-1j*t)
	x1 = Xr.real
	x2 = Xr.imag

	#### The impact parameter b
	b = x1 + 1j*f*x2
	#### The b-squared term
	bsq = np.absolute(b)
	bsq = np.multiply(bsq,bsq)
	#### The differentiation of bsq with respect to x
	bsqx = x1 + 1j*f*f*x2

	#### The deflection angle
	alpha = einRad*(np.sqrt(f)/fp)*( np.arctanh( fp * np.sqrt(bsq+(bc*bc)) / bsqx) - np.arctanh(fp*bc/(f*Xr)) )
	alpha = np.conj(alpha)*np.exp(1j*t) - xshear*np.conj(X)
	
	#### This part can return messages like
	#### "divide by zero encountered in divide"
	#### "invalid value encountered in divide"
	#### these messages can be ignored
	return Xr, alpha

    def shear(self, r=None, zs=None):
        rs = self.rs()
        k = self.rhoC()*self.deltaC()
        sigmaC = self.sigmaC(zs)
        x = r/rs

        if isinstance(r,np.ndarray):
            f = np.piecewise(x, [x>1., x<1., x==1.], [lambda x: (rs*k)*((8*np.arctan(((x-1)/(1+x))**0.5)/(x**2*(x**2-1)**0.5)) \
                    + (4*np.log(x/2.)/x**2) \
                    - (2./(x**2-1)) \
                    + (4*np.arctan(((x-1)/(1+x))**0.5)/((x**2-1)**1.5))), \
                    lambda x: (rs*k)*((8*np.arctanh(((1-x)/(1+x))**0.5)/(x**2*(1-x**2)**0.5)) \
                    + (4*np.log(x/2.)/x**2) \
                    - (2./(x**2-1)) \
                    + (4*np.arctanh(((1-x)/(1+x))**0.5)/((x**2-1)*(1-x**2)**0.5))), \
                    lambda x: (rs*k)*(10./3+4.*np.log(0.5))])
        else:
            if (x<1.):
                f = (rs*k)*((8*np.arctanh(((1-x)/(1+x))**0.5)/(x**2*(1-x**2)**0.5)) \
                    + (4*np.log(x/2.)/x**2) \
                    - (2./(x**2-1)) \
                    + (4*np.arctanh(((1-x)/(1+x))**0.5)/((x**2-1)*(1-x**2)**0.5)))
            elif (x>1.):
                f = (rs*k)*((8*np.arctan(((x-1)/(1+x))**0.5)/(x**2*(x**2-1)**0.5)) \
                    + (4*np.log(x/2.)/x**2) \
                    - (2./(x**2-1)) \
                    + (4*np.arctan(((x-1)/(1+x))**0.5)/((x**2-1)**1.5)))
            else:
                f = (rs*k)*(10./3+4.*np.log(0.5))

        #print f/sigmaC
        #print isinstance(f,np.ndarray)
        return f/sigmaC

def transfrom(N, E):
    N = np.float(N)
    xi = N/(A1*k0)
    eta = (E-E0)/(A1*k0)
    print xi, eta
    xip1 = h1*np.sin(2*xi)*np.cosh(2*eta)
    xip2 = h2*np.sin(4*xi)*np.cosh(4*eta)
    xip3 = h3*np.sin(6*xi)*np.cosh(6*eta)
    xip4 = h4*np.sin(8*xi)*np.cosh(8*eta)
    print xip1, xip2, xip3, xip4
    etap1 = h1*np.cos(2*xi)*np.sinh(2*eta)
    etap2 = h2*np.cos(4*xi)*np.sinh(4*eta)
    etap3 = h3*np.cos(6*xi)*np.sinh(6*eta)
    etap4 = h4*np.cos(8*xi)*np.sinh(8*eta)
    print etap1, etap2, etap4, etap4
    xip  = xi-xip1-xip2-xip3-xip4
    etap = eta-etap1-etap2-etap3-etap4
    print xip, etap
    beta = np.arcsin(sech(etap)*np.sin(xip))
    l = np.arcsin(np.tanh(etap)/np.cos(beta))
    print beta, l
    Q = np.arcsinh(np.tan(beta))
    print Q
    #e = np.e # ??
    e = np.sqrt(0.006694380023)
    Qp = Q + e*np.arctanh(e*np.tanh(Q))

    for i in xrange(3):
        Qp = Q + e*np.arctanh(e*np.tanh(Qp))

    print Qp

    rhoo = np.arctan(np.sinh(Qp))
    lamb = lamb0 + l
    return rhoo*180.0/np.pi, lamb*180.0/np.pi

def generateHeatMaps(GroupDF,goodsubj,GroupTrain=[]):

    numberOfICs=10
    columnNames=[]
    for rsnNumber in range(numberOfICs):
            columnNames.append('RSN%d' % rsnNumber)


    heatmapDF=GroupDF[GroupDF.Subject_ID.isin(goodsubj)].groupby(['Subject_ID','FB','TR']).mean()
    hmDiff=np.zeros((10,10,len(unique(GroupDF[GroupDF.Subject_ID.isin(goodsubj)]['Subject_ID']))))
    hmFB=hmDiff.copy()
    hmNFB=hmDiff.copy()
    if len(GroupTrain)>0:
        heatmapTrainDF=GroupTrain[GroupTrain.Subject_ID.isin(goodsubj)].groupby(['Subject_ID','TR']).mean()
        hmTrain=hmDiff.copy()
        hmFB_Train=hmDiff.copy()
        hmNFB_Train=hmDiff.copy()

    for indx,subj in enumerate(unique(GroupDF[GroupDF.Subject_ID.isin(goodsubj)]['Subject_ID'])):
        hmFB[:,:,indx]=heatmapDF.loc[subj,'FEEDBACK'][columnNames].corr()
        hmNFB[:,:,indx]=heatmapDF.loc[subj,'NOFEEDBACK'][columnNames].corr()
        hmDiff[:,:,indx]=np.arctanh(heatmapDF.loc[subj,'FEEDBACK'][columnNames].corr())*np.sqrt(405)-np.arctanh(heatmapDF.loc[subj,'NOFEEDBACK'][columnNames].corr())*np.sqrt(405)
        if len(GroupTrain)>0:
            hmTrain[:,:,indx]=heatmapTrainDF.loc[subj][columnNames].corr()
            hmFB_Train[:,:,indx]=np.arctanh(heatmapDF.loc[subj,'FEEDBACK'][columnNames].corr())*np.sqrt(405)-np.arctanh(heatmapTrainDF.loc[subj][columnNames].corr())*np.sqrt(175)
            hmNFB_Train[:,:,indx]=np.arctanh(heatmapDF.loc[subj,'NOFEEDBACK'][columnNames].corr())*np.sqrt(405)-np.arctanh(heatmapTrainDF.loc[subj][columnNames].corr())*np.sqrt(175)

    if len(GroupTrain)>0:
        return hmFB,hmNFB,hmDiff,hmTrain,hmFB_Train,hmNFB_Train
    else:
        return hmFB,hmNFB,hmDiff

    def kappa(self, r):
        from numpy import arctanh, arctan, arctan2, log, sin, cos, pi, logspace
        x = self.b / self.rs

        if x < 1.:
            norm = x**2 / (4 * arctanh(
                ((1 - x) / (1 + x))**0.5) / (1 - x**2)**0.5 + 2 * log(x / 2))
        elif x == 1.:
            norm = 1. / (2. + 2 * log(0.5))
        else:
            norm = x**2 / (4 * arctan(
                ((x - 1) / (x + 1))**0.5) / (x**2 - 1)**0.5 + 2 * log(x / 2))

        x = r / self.rs
        A = x * 0.
        C = x < 1.
        X = x[C].copy()
        A[C] = (1. - 2 * arctanh(
            ((1. - X) / (1. + X))**0.5) / (1 - X**2)**0.5) / (X**2 - 1.)
        C = x == 1.
        A[C] = 1. / 3
        C = x > 1.
        X = x[C].copy()
        A[C] = (1. - 2 * arctan(
            ((X - 1.) / (1. + X))**0.5) / (X**2 - 1.)**0.5) / (X**2 - 1.)
        return norm * A

def coherr(C,J1,J2,p=0.05,Nsp1=None,Nsp2=None):
    """
    Function to compute lower and upper confidence intervals on
    coherency (absolute value of coherence).

    C:            coherence (real or complex)
    J1,J2:        tapered fourier transforms
    p:            the target P value (default 0.05)
    Nsp1:         number of spikes in J1, used for finite size correction.
    Nsp2:         number of spikes in J2, used for finite size correction.
                  Default is None, for no correction

    Outputs:
    CI:           confidence interval for C, N x 2 array, (lower, upper)
    phi_std:      stanard deviation of phi, N array
    """
    from numpy import iscomplexobj, absolute, fix, zeros, setdiff1d, real, sqrt,\
         arctanh, tanh
    from scipy.stats import t

    J1 = _combine_trials(J1)
    J2 = _combine_trials(J2)
    N,K = J1.shape
    assert J1.shape==J2.shape, "J1 and J2 must have the same dimensions."
    assert N == C.size, "S and J lengths don't match"
    if iscomplexobj(C): C = absolute(C)

    pp = 1 - p/2
    dof = 2*K
    dof1 = dof if Nsp1 is None else fix(2.*Nsp1*dof/(2.*Nsp1+dof))
    dof2 = dof if Nsp2 is None else fix(2.*Nsp2*dof/(2.*Nsp2+dof))
    dof = min(dof1,dof2)

    Cerr = zeros((N,2))
    tcrit = t(dof-1).ppf(pp).tolist()
    atanhCxyk = zeros((N,K))
    phasefactorxyk = zeros((N,K),dtype='complex128')

    for k in xrange(K):
        indxk = setdiff1d(range(K),[k])
        J1k = J1[:,indxk]
        J2k = J2[:,indxk]
        eJ1k = real(J1k * J1k.conj()).sum(1)
        eJ2k = real(J2k * J2k.conj()).sum(1)
        eJ12k = (J1k.conj() * J2k).sum(1)
        Cxyk = eJ12k/sqrt(eJ1k*eJ2k)
        absCxyk = absolute(Cxyk)
        atanhCxyk[:,k] = sqrt(2*K-2)*arctanh(absCxyk)
        phasefactorxyk[:,k] = Cxyk / absCxyk

    atanhC = sqrt(2*K-2)*arctanh(C);
    sigma12 = sqrt(K-1)* atanhCxyk.std(1)

    Cu = atanhC + tcrit * sigma12
    Cl = atanhC - tcrit * sigma12
    Cerr[:,0] = tanh(Cl / sqrt(2*K-2))
    Cerr[:,1] = tanh(Cu / sqrt(2*K-2))
    phistd = (2*K-2) * (1 - absolute(phasefactorxyk.mean(1)))
    return Cerr, phistd

def get_ttest_map(map_name):
    subjs = map(lambda n: 'sub%03d' % n, SUBJECTS)
    vols = load_maps(map_name, subjs)
    # se = 1.0/np.sqrt(1199-6)
    popmean = np.arctanh(0.03)
    ttest = stats.ttest_1samp(np.arctanh(vols), popmean, axis=0)[0]
    ttest[ttest < 4] = np.NaN
    return get_volume(ttest, 4, 15)

 def function_A(vega, a, b, t):
     """Helper function A(t) defined in thesis."""
     kappa = 2 * b ** (-2) / vega
     val = (1 - sqrt( 1 + kappa * (1 - exp( -vega * t )))) \
         / kappa + 1 / sqrt( 1 + kappa) \
         * (arctanh( sqrt( 1 + kappa * ( 1 - exp( -vega * t)))/ \
         sqrt(1 + kappa)) - arctanh(1 / sqrt( 1 + kappa))) 
     return val

def test_arctanh():
    a = afnumpy.random.random((2,3))
    b = numpy.array(a)
    fassert(afnumpy.arctanh(a), numpy.arctanh(b))
    c = afnumpy.random.random((2,3))
    d = numpy.array(a)
    fassert(afnumpy.arctanh(a, out=c), numpy.arctanh(b, out=d))
    fassert(c, d)

def e1e2_to_g1g2(e1, e2):
    """
    convert e1,e2 to reduced shear style ellipticity

    parameters
    ----------
    e1,e2: tuple of scalars
        shapes in (ixx-iyy)/(ixx+iyy) style space

    outputs
    -------
    g1,g2: scalars
        Reduced shear space shapes

    """

    e = numpy.sqrt(e1*e1 + e2*e2)
    if isinstance(e1, numpy.ndarray):
        w,=numpy.where(e >= 1.0)
        if w.size != 0:
            raise GMixRangeError("some e were out of bounds")

        eta=numpy.arctanh(e)
        g = numpy.tanh(0.5*eta)

        numpy.clip(g, 0.0, 0.99999999, g)

        g1=numpy.zeros(g.size)
        g2=numpy.zeros(g.size)
        w,=numpy.where(e != 0.0)
        if w.size > 0:
            fac = g[w]/e[w]

            g1[w] = fac*e1[w]
            g2[w] = fac*e2[w]

    else:
        if e >= 1.:
            raise GMixRangeError("e out of bounds: %s" % e)
        if e == 0.0:
            g1,g2=0.0,0.0

        else:

            eta=numpy.arctanh(e)
            g = numpy.tanh(0.5*eta)

            if g >= 1.:
                # round off?
                g = 0.99999999


            fac = g/e

            g1 = fac*e1
            g2 = fac*e2

    return g1,g2

    def make_from_examples(cls, X, low, high, directed = True):
        #for every pair of examples (i,j)
        #make a feature that takes on the value low at i
        #and value high at j
        #if directed, order of (i,j) matters, otherwise it does not


        m,n =  X.shape

        if directed:
            h = m **2 - m
        else:
            h = m * (m-1)/2
        W = N.zeros((n,h))
        b = N.zeros(h)
        idx = 0

        inv_low = N.arctanh(low)
        if N.abs(N.tanh(inv_low)-low) > 1e-6:
            assert False
        #

        inv_high = N.arctanh(high)

        for i in xrange(X.shape[0]):
            if directed:
                r = xrange(m)
            else:
                r = xrange(i+1,m)

            for j in r:
                if i == j:
                    continue

                diff = X[j,:] - X[i,:]
                direction = diff / N.sqrt(N.square(diff).sum())
                pi = N.dot(X[i,:],direction)
                pj = N.dot(X[j,:],direction)

                wmag =  (inv_high - inv_low) / (pj - pi)

                b[idx] = (pj*inv_low - pi*inv_high) / (pj - pi)
                W[:,idx] = wmag * direction

                #check it
                ival = N.tanh(N.dot(W[:,idx],X[i,:])+b[idx])
                jval = N.tanh(N.dot(W[:,idx],X[j,:])+b[idx])


                assert abs(ival-low) < 1e-6
                assert abs(jval-high) < 1e-6

                idx += 1

        assert idx == h


        return TanhFeatureExtractor(W,b)

def g1g2_to_e1e2(g1, g2):
    """
    convert reduced shear g1,g2 to standard ellipticity
    parameters e1,e2

    uses eta representation but could also use
        e1 = 2*g1/(1 + g1**2 + g2**2)
        e2 = 2*g2/(1 + g1**2 + g2**2)

    parameters
    ----------
    g1,g2: scalars
        Reduced shear space shapes

    outputs
    -------
    e1,e2: tuple of scalars
        shapes in (ixx-iyy)/(ixx+iyy) style space
    """
    g=numpy.sqrt(g1*g1 + g2*g2)

    if isinstance(g1, numpy.ndarray):
        w,=numpy.where(g >= 1.0)
        if w.size != 0:
            raise GMixRangeError("some g were out of bounds")

        eta = 2*numpy.arctanh(g)
        e = numpy.tanh(eta)

        numpy.clip(e, 0.0, 0.99999999, e)

        e1=numpy.zeros(g.size)
        e2=numpy.zeros(g.size)
        w,=numpy.where(g != 0.0)
        if w.size > 0:
            fac = e[w]/g[w]

            e1[w] = fac*g1[w]
            e2[w] = fac*g2[w]

    else:
        if g >= 1.:
            raise GMixRangeError("g out of bounds: %s" % g)
        if g == 0.0:
            return (0.0, 0.0)

        eta = 2*numpy.arctanh(g)
        e = numpy.tanh(eta)
        if e >= 1.:
            # round off?
            e = 0.99999999

        fac = e/g

        e1 = fac*g1
        e2 = fac*g2

    return e1,e2

def diffcorrcoeftest(rvalue1, rvalue2, N1, N2):
    ''' tests for difference between 2 r values by performing fishers transformation and performing t-test on the z values. returns z and two tailed p '''
    r_z1=np.arctanh(rvalue1) #equivalent to 0.5 * np.log((1 + rvalue1)/(1 - rvalue1))
    r_z2=np.arctanh(rvalue2)
    se_diff_r = np.sqrt(1.0/(N1 - 3.0) + 1.0/(N2 - 3.0))
    diff = r_z1 - r_z2
    z = abs(diff / se_diff_r)
    p = (1 - scipy.stats.norm.cdf(z))*2
    return round(z,3), round(p,3)

def atanh(x):
    """
    Inverse hyperbolic tangent
    """
    if isinstance(x, UncertainFunction):
        mcpts = np.arctanh(x._mcpts)
        return UncertainFunction(mcpts)
    else:
        return np.arctanh(x)

def fisherTransformation(r1, r2, n1, n2):
    z1 = numpy.arctanh(r1)
    z2 = numpy.arctanh(r2)
    
    se = numpy.sqrt(1/(float(n1)-3) + 1/(float(n2)-3))
    
    z = (z1-z2)/se
    
    return z

    def test_arctanh(self):
        import math
        from numpy import arctanh

        for v in [0.99, 0.5, 0, -0.5, -0.99]:
            assert math.atanh(v) == arctanh(v)
        for v in [2.0, -2.0]:
            assert math.isnan(arctanh(v))
        for v in [1.0, -1.0]:
            assert arctanh(v) == math.copysign(float("inf"), v)

    def compute_S(x):
        ''' compute the nonlinearity of the belief prop system
            S(x) is obtained by applying arctanh to x at i in S_ '''

        global S_, q

        y = np.ones(q,1)
        for i in (i,q):
            for j in nonzero(S_[i,:]):
                y[i] *= np.arctanh(x[j]/2)
        y = 2 * np.arctanh(y)