def int_pot_2D_moi(self, xp, yp, x, R, h, basis_func):
     """FWD model function. Incorporates the Method of Images.
     Returns contribution of a point xp,yp, belonging to a basis source
     support centered at (0,0) to the potential measured at (x,0),
     integrated over xp,yp gives the potential generated by a
     basis source element centered at (0,0) at point (x,0)
     #Eq 20, Ness(2015)
     Parameters
     ----------
     xp, yp : floats or np.arrays
         point or set of points where function should be calculated
     x :  float
         position at which potential is being measured
     R : float
         The size of the basis function
     h : float
         thickness of slice
     basis_func : method
         Fuction of the basis source
     Returns
     -------
     pot : float
     """
     L = ((x-xp)**2 + yp**2)**(0.5)
     if L < 0.00001:
         L = 0.00001
     correction = np.arcsinh((h-(2*h*self.iters))/L) + np.arcsinh((h+(2*h*self.iters))/L)
     pot = np.arcsinh(h/L) + np.sum(self.iter_factor*correction)
     dist = np.sqrt(xp**2 + yp**2)
     pot *= basis_func(dist, R)  # Eq 20, Ness et.al.
     return pot

def plot_IQU(solution, title, col, ncol=6, coord='C'):
    # Es=solution[np.array(final_index).tolist()].reshape((4, len(final_index)/4))
    # I = Es[0] + Es[3]
    # Q = Es[0] - Es[3]
    # U = Es[1] + Es[2]
    IQUV = sol2map(solution)
    IQUV.shape = (4, IQUV.shape[0] / 4)
    I = IQUV[0]
    Q = IQUV[1]
    U = IQUV[2]
    V = IQUV[3]
    pangle = 180 * np.arctan2(Q, U) / 2 / PI
    plotcoordtmp = coord
    hpv.mollview(np.log10(I), min=0, max=4, coord=plotcoordtmp, title=title, nest=True, sub=(4, ncol, col))

    hpv.mollview((Q ** 2 + U ** 2) ** .5 / I, min=0, max=1, coord=plotcoordtmp, title=title, nest=True,
                 sub=(4, ncol, ncol + col))
    from matplotlib import cm
    cool_cmap = cm.hsv
    cool_cmap.set_under("w")  # sets background to white
    hpv.mollview(pangle, min=-90, max=90, coord=plotcoordtmp, title=title, nest=True, sub=(4, ncol, 2 * ncol + col),
                 cmap=cool_cmap)

    hpv.mollview(np.arcsinh(V) / np.log(10), min=-np.arcsinh(10. ** 4) / np.log(10),
                 max=np.arcsinh(10. ** 4) / np.log(10), coord=plotcoordtmp, title=title, nest=True,
                 sub=(4, ncol, 3 * ncol + col))
    if col == ncol:
        plt.show()

def Scale_asinh(inputArray, scale_min=None, scale_max=None, non_linear=2.0):
    """Scale_asinh(inputArray, scale_min=None, scale_max=None, non_linear=2.0)
    Performs asinh scaling of the input numpy array.
    (from Min-Su Shin, Princeton)
    
    inputArray: image data array
    scale_min (None): minimum data value
        use inputArray.min() if None
    scale_max (None): maximum data value
        use inputArray.max() if None
    non_linear (2.0): non-linearity factor
    
    >>> scaledArray = Scale_asinh(inputArray)
    """
    imageData=numpy.array(inputArray, copy=True)
    
    if scale_min == None:
        scale_min = imageData.min()
    if scale_max == None:
        scale_max = imageData.max()
    
    factor = numpy.arcsinh((scale_max - scale_min)/non_linear)
    indices0 = numpy.where(imageData < scale_min)
    indices1 = numpy.where((imageData >= scale_min) & (imageData <= scale_max))
    indices2 = numpy.where(imageData > scale_max)
    imageData[indices0] = 0.0
    imageData[indices2] = 1.0
    imageData[indices1] = numpy.arcsinh((imageData[indices1] - scale_min)/non_linear)/factor
    return imageData

 def imstretch(self):
     data = np.clip(self.data_array, self.threshold[0], self.threshold[1])
     if self.mode == "linear":
         pass
     elif self.mode == "logarithmic":
         data = np.reciprocal(1 + np.power(0.5 / data, self.factor))
     elif self.mode == "gamma":
         data = np.power(data, self.factor)
     elif self.mode == "arcsinh":
         mn = np.nanmin(data)
         mx = np.nanmax(data)
         tmp = bytescale(data, high=1.0)
         beta = np.clip(self.factor, 0.0, self.factor)
         sclbeta = (beta - mn) / (mx - mn)
         sclbeta = np.clip(sclbeta, 1.0e-12, sclbeta)
         nonlinearity = 1.0 / sclbeta
         extrema = np.arcsinh(np.array([0.0, nonlinearity]))
         data = np.clip(np.arcsinh(data * nonlinearity), extrema[0], extrema[1])
     elif self.mode == "square root":
         data = np.sqrt(np.fabs(data)) * np.sign(data)
     elif self.mode == "histogram equalization":
         imhist, bins = np.histogram(data.flatten(), 256, normed=True)
         cdf = imhist.cumsum()  # cumulative distribution function
         cdf = 255 * cdf / cdf[-1]  # normalize
         im2 = np.interp(data.flatten(), bins[:-1], cdf)
         data = im2.reshape(data.shape)
     self.scaled = bytescale(data).flatten().tolist()

def asinh(inputArray, scale_min=None, scale_max=None, non_linear=2.0):
	"""Performs asinh scaling of the input numpy array.

	@type inputArray: numpy array
	@param inputArray: image data array
	@type scale_min: float
	@param scale_min: minimum data value
	@type scale_max: float
	@param scale_max: maximum data value
	@type non_linear: float
	@param non_linear: non-linearity factor
	@rtype: numpy array
	@return: image data array
	
	"""		
    
	print "img_scale : asinh"
	imageData=numpy.array(inputArray, copy=True)
	
	if scale_min == None:
		scale_min = imageData.min()
	if scale_max == None:
		scale_max = imageData.max()
	factor = numpy.arcsinh((scale_max - scale_min)/non_linear)
	indices0 = numpy.where(imageData < scale_min)
	indices1 = numpy.where((imageData >= scale_min) & (imageData <= scale_max))
	indices2 = numpy.where(imageData > scale_max)
	imageData[indices0] = 0.0
	imageData[indices2] = 1.0
	imageData[indices1] = numpy.arcsinh((imageData[indices1] - \
	scale_min)/non_linear)/factor

	return imageData

    def load_minibatch(self, filepath, nimg, farts, gridsize, cg, num,cg_additional=1,twoclasses=False):
        """
        Load a mini batch of images and their labels. 
        Labels need to be converted to tensorflow
        format

        inputs:
        filepath -- Path where the files are located
        nimg -- Number of images in the total batch
        farts -- Fraction of artifacts
        gridsize -- Number of pixels to a side
        cg -- Coarsegraining factor 
        num -- The minibatch number 
        cg_additional -- additional coursegraining to perform on the fly
        """
            
        X = np.load('{0}/X_{1}_{2}_{3}_{4}_mb{5}.npy'.format(filepath, nimg, farts, gridsize, cg, num))
        y = np.load('{0}/y_{1}_{2}_{3}_{4}_mb{5}.npy'.format(filepath, nimg, farts, gridsize, cg, num))
        X[X==-99] = np.nan
        if cg_additional!=1:
            X = np.mean(np.mean(X.reshape([X.shape[0],gridsize//cg,gridsize//cg//cg_additional,cg_additional]),axis=3).T.reshape(gridsize//cg//cg_additional,gridsize//cg//cg_additional,cg_additional,X.shape[0]),axis=2).T.reshape([X.shape[0],(gridsize//cg//cg_additional)**2])
        X = 255*(np.arcsinh(X)-np.atleast_2d(np.arcsinh(np.nanmin(X,axis=1))).T)/np.atleast_2d((np.arcsinh(np.nanmax(X,axis=1))-np.arcsinh(np.nanmin(X,axis=1)))).T
        X[np.isnan(X)] = 0
        #X -= np.atleast_2d(np.mean(X,axis=1)).T
        #print(np.nanmean(X, axis=1))

        ey = self.convert_labels(y, twoclasses)

        return X, ey

    def test_arcsinh(self):
        import math
        from numpy import arcsinh

        for v in [float("inf"), float("-inf"), 1.0, math.e]:
            assert math.asinh(v) == arcsinh(v)
        assert math.isnan(arcsinh(float("nan")))

 def test_asinh(self):
     """Test arcsinh scaling."""
     img = scale_image(DATA, scale='asinh')
     mean, median, stddev = sigmaclip_stats(DATA, sigma=3.0)
     z = (mean + (2.0 * stddev)) / 2.
     ref = np.arcsinh(DATASCL / z) / np.arcsinh(1.0 / z)
     assert_allclose(img, ref, atol=0, rtol=1.e-5)

    def __call__(self, value):
        self.autoscale_None(value)  # set vmin, vmax if unset
        inverted = self.vmax <= self.vmin

        hi, lo = max(self.vmin, self.vmax), min(self.vmin, self.vmax)
        ra = hi - lo
        mid = lo + ra * self.bias
        mn = mid - ra * self.contrast
        mx = mid + ra * self.contrast

        if self.stretch == "linear":
            result = (value - mn) * (1.0 / (mx - mn))
            result = np.clip(result, 0, 1)
        elif self.stretch == "arcsinh":
            b = max(self.bias, 1e-5)
            c = self.contrast
            result = (value - lo) / (1.0 * (hi - lo))
            result = np.arcsinh(result / b) / np.arcsinh((b + c) / b)
            result = np.clip(result, 0, 1)
        elif self.stretch == "sqrt":
            result = (value - mn) * (1.0 / (mx - mn))
            result = np.clip(result, 0, 1)
            result = np.sqrt(result)
        else:
            raise TypeError("Invalid stretch: %s" % self.stretch)

        if inverted:
            result = 1 - result

        return result

def scale_two_arcsinh(x,up1,up2,down1,down2,m="normal"):
    
    if m != "inverse":
        if x >= 0: return up1*np.arcsinh(x*up2)
        if x < 0: return down1*np.arcsinh(x*down2)
    else:
        if x >= 0: return 1./up2*np.sinh(x/up1)
        if x < 0: return 1./down2*np.sinh(x/down1)

 def newspace(high):
     Smax = high
     K = exact
     deps = 1./size * (np.arcsinh((Smax - K)*(1/density)) - np.arcsinh(-K/density))
     eps = np.arcsinh(-K/density) + np.arange(size)*deps
     space = K + density * np.sinh(eps)
     space -= min(space)
     return space

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

    def __call__(self, values, out=None, clip=True):

        values = _prepare(values, out=out, clip=clip)

        np.true_divide(values, self.a, out=values)
        np.arcsinh(values, out=values)
        np.true_divide(values, np.arcsinh(1. / self.a), out=values)

        return values

def elec_catenary_hyperbolic_lowest_proj ( T0,  w, l, h):
    Lh0 = elec_catenary_hyperbolic_length_equal_high (T0, w, l)
    if Lh0 != 0.:
        a = l/2. - T0/w*np.arcsinh(h/Lh0) 
        b = l/2. + T0/w*np.arcsinh(h/Lh0)
    else:
        a = l/2.
        b = l/2.
    return a

def asinh(x):
    """
    Inverse hyperbolic sine
    """
    if isinstance(x, UncertainFunction):
        mcpts = np.arcsinh(x._mcpts)
        return UncertainFunction(mcpts)
    else:
        return np.arcsinh(x)

def colorImage(b,g,r,bMinusr=0.8,bMinusg=0.4,sdev=None,nonlin=5.,m=0.5,M=None):
    w = r.shape[0]/2-5
    rb = r/b
    gb = g/b
    rnorm = numpy.median(rb[w:-w,w:-w])
    gnorm = numpy.median(gb[w:-w,w:-w])
    r /= rnorm
    g /= gnorm
    r *= 10**(0.4*bMinusr)
    g *= 10**(0.4*bMinusg)

    r /= 620.
    g /= 540.
    b /= 460.

    I = (r+g+b)/3.

    if sdev is None:    
        sdev = clip(I)[1]
    m = m*sdev
    if M is None:
        M = I[w:-w,w:-w].max()
    nonlin = nonlin*sdev

    f = numpy.arcsinh((I-m)/nonlin)/numpy.arcsinh((M-m)/nonlin)
    f[I<m] = 0.
    f[I>M] = 1.
    R = r*f/I
    G = g*f/I
    B = b*f/I

    R[I<=0] = 0.
    G[I<=0] = 0.
    B[I<=0] = 0.

    R[R<=0] = 0.
    G[G<=0] = 0.
    B[B<=0] = 0.

    R[R>1] = 1.
    G[G>1] = 1.
    B[B>1] = 1.

    white = True
    if white:
        cond = (f==1)
        R[cond] = 1.
        G[cond] = 1.
        B[cond] = 1.

    arr = numpy.empty((R.shape[0],R.shape[1],3))
    arr[:,:,0] = R
    arr[:,:,1] = G
    arr[:,:,2] = B

    return arr,sdev,M,rnorm,gnorm

def ref_table(PTEN,MWPL,S,FH,AE,MBL):
    n = 754.
    max_a =(MBL-MWPL*FH)/MWPL
    INC = max_a/n
    a = np.array(np.linspace(0.0,INC*n,num=n+1))
    H = MWPL*a
    Ttop = H + MWPL*FH
    Vtop = (Ttop**2-H**2)**0.5
    ang = 90. - np.arcsin(H/Ttop)*180/np.pi
    sp_temp = -(S/2.)+(FH/2.)*(1.+(4*a**2/(S**2.-FH**2)))**0.5 
    sp_temp = [0 if i < 0 else i for i in sp_temp]
    # INITIALIZE ARRanchor_yS
    Vbot = np.array([0.]*len(a))
    Tbot = np.array([0.]*len(a))
    Tave = np.array([0.]*len(a))
    stretch = np.array([0.]*len(a))
    x = np.array([0.]*len(a))
    s = np.array([0.]*len(a))
    yp = np.array([0.]*len(a))
    sp = np.array([0.]*len(a))
    for i in range (0,len(a)):
        if sp_temp[i] == 0.: 
            Vbot[i] = 0.
            Tbot[i] = H[i]
            Tave[i] = 0.5*(Ttop[i]+Tbot[i])
            stretch[i] = 1.+Tave[i]/AE
            if i == 0:
                x[i] = S*stretch[i] -FH
            else: 
                x[i] = S*stretch[i-1] -FH*(1.+2.*a[i]/FH)**0.5+a[i]*np.arccosh(1.+FH/a[i])
            s[i] = (FH**2+2.*FH*a[i])**0.5
            sp [i] = 0.
        else: 
            s[i] = S*stretch[i-1]
            Vbot[i] = Vtop[i]-MWPL*s[i]
            Tbot[i] = (H[i]**2+Vbot[i]**2)**0.5
            Tave[i] =  0.5*(Ttop[i]+Tbot[i])
            stretch[i] = 1.+Tave[i]/AE
            sp[i] = sp_temp[i]*stretch[i-1]
            x[i] = a[i]*(np.arcsinh((S+sp[i])/a[i])-np.arcsinh(sp[i]/a[i]))*stretch[i-1]
            
        if i == 0: 
            yp[i] = -a[i]+(a[i]**2+sp[i]**2)**0.5
        else: 
            yp[i] = (-a[i]+(a[i]**2+sp[i]**2)**0.5)*stretch[i-1]
    x0 = np.interp(PTEN,Ttop,x)
    offset = x - x0
    MKH = np.array([0.]*len(a))
    MKV = np.array([0.]*len(a))
    for i in range(1,len(a)):
        MKH[i] = (H[i]-H[i-1])/(offset[i]-offset[i-1])
        MKV[i] = (Vtop[i]-Vtop[i-1])/(offset[i]-offset[i-1])
    return x0,a,x,H,sp,yp,s,Ttop,Vtop,Tbot,Vbot,Tave,stretch,ang,offset,MKH,MKV,INC

def asinhmag(flux, fluxerr,  m0 = 22.5, f0=1.0, b=0.01):
	"""
	Implements
	http://ssg.astro.washington.edu/elsst/opsim.shtml?lightcurve_mags
	"""

	mag = m0 -(2.5/np.log(10.)) * ( np.arcsinh( flux / (f0 * 2.0 * b)) + np.log(b) )
	
	magplu = m0 -(2.5/np.log(10.)) * ( np.arcsinh( (flux+fluxerr) / (f0 * 2.0 * b)) + np.log(b) )
	magmin = m0 -(2.5/np.log(10.)) * ( np.arcsinh( (flux-fluxerr) / (f0 * 2.0 * b)) + np.log(b) )
	magerr = 0.5*(magmin - magplu)
	
	return (mag, magerr)

def arcsinh(x, out=None):
    """
    Raises a ValueError if input cannot be rescaled to a dimensionless
    quantity.
    """
    if not isinstance(x, Quantity):
        return np.arcsinh(x, out)

    return Quantity(
        np.arcsinh(x.rescale(dimensionless).magnitude, out),
        dimensionless,
        copy=False
    )

def arcsinhspace(start, stop, num=50):
    """
    Return numbers spaced evenly on an arcsinh scale.

    Parameters
    ----------
    start : number
        start value
    stop : number
        stop/end value, inclusive
    num : number
        number of intervales between start and stop
    """
    return _np.sinh(_np.linspace(_np.arcsinh(start), _np.arcsinh(stop), num))