In [2]:
import healpix_util as hu
import astropy as ap
import math, sys, time
import numpy as np
from astropy.io import fits
from astropy.table import Table
import astropy.io.ascii as ascii
from astropy.constants import c
import matplotlib.pyplot as plt
import math
import scipy.special as sp
from scipy import integrate
from astropy.stats import bayesian_blocks
from astropy.stats import histogram
from astropy.visualization import hist

%matplotlib inline

In [2]:
#Read sorted redshift, RA, DEC values of each quasar/AGN
hmqdat=ascii.read("./hmqdata_sorted_full.csv")
hmqdat


Out[2]:
<Table length=510355>
ZRADEC
float64float64float64
-0.002178.95169161.598303
-0.001143.5277939.44227
-0.001217.26922112.069508
0.02.2663936.472577
0.07.306389-28.213611
0.010.10375-23.666889
0.017.522778-31.037222
0.025.4575-31.127222
0.044.53027820.500556
0.044.99-29.338889
.........
6.2512.52777834.756389
6.29177.0137087.035639
6.305157.6130565.415278
6.43352.284722-3.033056
6.43832.555-4.939167
6.44177.06936352.863975
6.60446.3205-31.848889
6.74517.471375-30.790639
6.886357.138917-30.902778
7.085170.0061676.690083

In [6]:
#Creating output file 
findpair = open("./output/findpair.dat",'w')
# We intend to obtain average redshift of the pair and deltaZ deltatheta Zdeltatheta values first essentially finding right pairs
findpair.write("Z\t deltaZ\t deltatheta\t Zdeltatheta\n")
#Algorithm keeps removing rows after done finding quasars within deltaZ˜0.02 
start_time=time.time()
while len(hmqdat['Z'])>0:
    j=1
    while hmqdat[j]['Z']<=(hmqdat[0]['Z']+0.02):
        deltaZ=abs(hmqdat[j]['Z']-hmqdat[0]['Z'])
        Z=(hmqdat[j]['Z']+hmqdat[0]['Z'])/2.0
        RA1=hmqdat[0]['RA']*math.pi/180.0
        RA2=hmqdat[j]['RA']*math.pi/180.0
        DEC1=hmqdat[0]['DEC']*math.pi/180.0
        DEC2=hmqdat[j]['DEC']*math.pi/180.0
        deltatheta=abs(math.acos(math.sin(DEC1)*math.sin(DEC2)+math.cos(DEC1)*math.cos(DEC2)*math.cos(RA1-RA2)))
        Zdeltatheta=Z*deltatheta
        if Zdeltatheta<=0.02:
            findpair.write("%f\t %f\t %f\t %f\n" %(Z, deltaZ, deltatheta, Zdeltatheta))
        j=j+1
    hmqdat.remove_row(0)
findpair.close()
print "Total time taken to run: "
print time.time()-start_time


---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-6-ebf78828843e> in <module>()
      7 while len(hmqdat['Z'])>0:
      8     j=1
----> 9     while hmqdat[j]['Z']<=(hmqdat[0]['Z']+0.02):
     10         deltaZ=abs(hmqdat[j]['Z']-hmqdat[0]['Z'])
     11         Z=(hmqdat[j]['Z']+hmqdat[0]['Z'])/2.0

/Users/rohin/anaconda/lib/python2.7/site-packages/astropy/table/table.pyc in __getitem__(self, item)
   1196             return self.columns[item]
   1197         elif isinstance(item, (int, np.integer)):
-> 1198             return self.Row(self, item)
   1199         elif (isinstance(item, (tuple, list)) and item and
   1200               all(isinstance(x, six.string_types) for x in item)):

/Users/rohin/anaconda/lib/python2.7/site-packages/astropy/table/row.pyc in __init__(self, table, index)
     41         if index < -n or index >= n:
     42             raise IndexError('index {0} out of range for table with length {1}'
---> 43                              .format(index, len(table)))
     44 
     45     def __getitem__(self, item):

IndexError: index 1 out of range for table with length 1

In [ ]:
#Creating output file 
Ralpha = open("./output/Ralpha.dat",'w')
# We intend to obtain average redshift of the pair and deltaZ deltatheta Zdeltatheta R alpha values
Ralpha.write("Z\t deltaZ\t deltatheta\t Zdeltatheta\t R\t alpha\t \n")
#Algorithm keeps removing rows after done finding quasars within deltaZ˜0.02 
start_time=time.time()
while len(hmqdat['Z'])>0:
    j=1
    while hmqdat[j]['Z']<=(hmqdat[0]['Z']+0.02):
        deltaZ=abs(hmqdat[j]['Z']-hmqdat[0]['Z'])
        Z=(Z1+Z2)/2.0
        RA1=hmqdat[0]['RA']*math.pi/180.0
        RA2=hmqdat[j]['RA']*math.pi/180.0
        DEC1=hmqdat[0]['DEC']*math.pi/180.0
        DEC2=hmqdat[j]['DEC']*math.pi/180.0
        deltatheta=abs(math.acos(math.sin(DEC1)*math.sin(DEC2)+math.cos(DEC1)*math.cos(DEC2)*math.cos(RA1-RA2)))
        Zdeltatheta=Z*deltatheta
            if Zdeltatheta<=0.02:
                R=math.sqrt(deltaZ**2+Zdeltatheta**2)
                alpha=math.acos(deltaZ/R)
                Ralpha.write("%f\t %f\t %f\t %f\t %f\t %f\n" %(Z, deltaZ, deltatheta, Zdeltatheta, R, alpha))
        j=j+1
    hmqdat.remove_row(0)
Ralpha.close()
print "Total time taken to run: "
print time.time()-start_time

In [7]:
j


Out[7]:
1

In [8]:
hmqdat


Out[8]:
<Table length=1>
ZRADEC
float64float64float64
7.085170.0061676.690083

In [9]:
pairdat=ascii.read("./output/findpair.dat")

In [10]:
pairdat


Out[10]:
<Table length=6523551>
ZdeltaZdeltathetaZdeltatheta
float64float64float64float64
-0.00150.0010.539116-0.000809
-0.00150.0010.989785-0.001485
-0.0010.0021.952679-0.001953
-0.0010.0022.550885-0.002551
-0.0010.0022.466295-0.002466
-0.0010.0022.567833-0.002568
-0.0010.0022.530615-0.002531
-0.0010.0021.574572-0.001575
-0.0010.0022.372919-0.002373
-0.0010.0021.998327-0.001998
............
3.57550.0050.0035050.01253
3.5740.0020.005520.019727
3.57950.0110.0011660.004175
3.5820.0140.0034150.012232
3.58550.0070.0044150.01583
3.5960.0160.0020470.00736
3.59850.0170.0042210.015191
3.60.0180.003320.011952
3.5960.0060.003650.013124
3.5960.00.0029010.010434

In [19]:
pairdat.sort('Z')

In [16]:
pairdats=pairdat

In [17]:
pairdats


Out[17]:
<Table length=6523551>
ZdeltaZdeltathetaZdeltatheta
float64float64float64float64
-0.00150.0010.539116-0.000809
-0.00150.0010.989785-0.001485
-0.0010.0020.524223-0.000524
-0.0010.0020.729136-0.000729
-0.0010.0021.205034-0.001205
-0.0010.0020.938176-0.000938
-0.0010.0021.216635-0.001217
-0.0010.0020.992697-0.000993
-0.0010.0020.769144-0.000769
-0.0010.0020.686675-0.000687
............
3.57550.0050.0035050.01253
3.5780.0160.0036390.01302
3.57950.0110.0011660.004175
3.5820.0140.0034150.012232
3.58550.0070.0044150.01583
3.5960.0060.003650.013124
3.5960.0160.0020470.00736
3.5960.00.0029010.010434
3.59850.0170.0042210.015191
3.60.0180.003320.011952

In [22]:
pair_blocks=hist(pairdats['Z'],bins='blocks')



In [21]:
pairs_knuth=hist(pairdats['Z'],bins='knuth')



In [23]:
def AutoCorrelation(x):
    #x = np.asarray(x)
    y = x-x.mean()
    result = np.correlate(y, y, mode='full')
    result = result[len(result)//2:]
    result /= result[0]
    return result

In [24]:
DD=AutoCorrelation(pairs_knuth[0])

In [26]:
plt.plot(pairs_knuth[1][1:len(pairs_knuth[1])],DD)
plt.savefig("./images/DD_af_knuth.pdf")
plt.savefig("./images/DD_af_knuth.jpeg")



In [27]:
help(Table.add_column)


Help on method add_column in module astropy.table.table:

add_column(self, col, index=None, rename_duplicate=False) unbound astropy.table.table.Table method
    Add a new Column object ``col`` to the table.  If ``index``
    is supplied then insert column before ``index`` position
    in the list of columns, otherwise append column to the end
    of the list.
    
    Parameters
    ----------
    col : Column
        Column object to add.
    index : int or `None`
        Insert column before this position or at end (default)
    rename_duplicate : bool
        Uniquify column name if it already exist (default=False)
    
    Examples
    --------
    Create a table with two columns 'a' and 'b'::
    
        >>> t = Table([[1, 2, 3], [0.1, 0.2, 0.3]], names=('a', 'b'))
        >>> print(t)
         a   b
        --- ---
          1 0.1
          2 0.2
          3 0.3
    
    Create a third column 'c' and append it to the end of the table::
    
        >>> col_c = Column(name='c', data=['x', 'y', 'z'])
        >>> t.add_column(col_c)
        >>> print(t)
         a   b   c
        --- --- ---
          1 0.1   x
          2 0.2   y
          3 0.3   z
    
    Add column 'd' at position 1. Note that the column is inserted
    before the given index::
    
        >>> col_d = Column(name='d', data=['a', 'b', 'c'])
        >>> t.add_column(col_d, 1)
        >>> print(t)
         a   d   b   c
        --- --- --- ---
          1   a 0.1   x
          2   b 0.2   y
          3   c 0.3   z
    
    Add second column named 'b' with rename_duplicate::
    
        >>> t = Table([[1, 2, 3], [0.1, 0.2, 0.3]], names=('a', 'b'))
        >>> col_b = Column(name='b', data=[1.1, 1.2, 1.3])
        >>> t.add_column(col_b, rename_duplicate=True)
        >>> print(t)
         a   b  b_1
        --- --- ---
          1 0.1 1.1
          2 0.2 1.2
          3 0.3 1.3
    
    To add several columns use add_columns.


In [30]:
R=np.sqrt(pairdats['deltaZ']**2+(pairdats['Zdeltatheta'])**2)

In [31]:
R


Out[31]:
<Column name='deltaZ' dtype='float64' length=6523551>
0.00128606128869
0.00179004672537
0.0022292269767
0.0021520895157
0.002311934458
0.00208023630539
0.002285352415
0.0028272443826
0.00310334409638
0.0025454423944
0.00322552511666
0.003254806648
...
0.0198295971495
0.0141191705445
0.0134927469285
0.0206283616848
0.0117651921637
0.0185912557457
0.0173086205675
0.0176120554639
0.0144317748444
0.010431996
0.0227972339893
0.021606718955

In [34]:
Rdat=pairdats.Column(name='R',data=R)
pairdats.add_column(Rdat)

In [35]:
pairdats


Out[35]:
<Table length=6523551>
ZdeltaZdeltathetaZdeltathetaR
float64float64float64float64float64
-0.00150.0010.539116-0.0008090.00128606128869
-0.00150.0010.989785-0.0014850.00179004672537
-0.0010.0020.984608-0.0009850.0022292269767
-0.0010.0020.794663-0.0007950.0021520895157
-0.0010.0021.159759-0.001160.002311934458
-0.0010.0020.572174-0.0005720.00208023630539
-0.0010.0021.105819-0.0011060.002285352415
-0.0010.0021.998327-0.0019980.0028272443826
-0.0010.0022.372919-0.0023730.00310334409638
-0.0010.0021.574572-0.0015750.0025454423944
...............
3.57550.0050.0035050.012530.0134927469285
3.5780.0160.0036390.013020.0206283616848
3.57950.0110.0011660.0041750.0117651921637
3.5820.0140.0034150.0122320.0185912557457
3.58550.0070.0044150.015830.0173086205675
3.5960.0160.0020470.007360.0176120554639
3.5960.0060.003650.0131240.0144317748444
3.5960.00.0029010.0104340.010431996
3.59850.0170.0042210.0151910.0227972339893
3.60.0180.003320.0119520.021606718955

In [36]:
alpha=np.arccos(pairdats['deltaZ']/pairdats['R'])


/Users/rohin/anaconda/lib/python2.7/site-packages/ipykernel/__main__.py:1: RuntimeWarning: invalid value encountered in divide
  if __name__ == '__main__':

In [37]:
alpha


Out[37]:
<Column name='deltaZ' dtype='float64' length=6523551>
0.680007633186
0.978045564681
0.457471875897
0.378203831404
0.525493621222
0.278644209975
0.505074247406
0.784979738416
0.870471520917
0.666940281126
0.901984698482
0.909074940238
...
1.46976520295
1.20881436555
1.19117426594
0.683081894981
0.362646215648
0.718122811461
1.15444874925
0.431190943163
1.1420298228
1.57079632679
0.729204740345
0.586154176817

In [38]:
alphadat=pairdats.Column(name='alpha',data=alpha)
pairdats.add_column(alphadat)

In [39]:
pairdats


Out[39]:
<Table length=6523551>
ZdeltaZdeltathetaZdeltathetaRalpha
float64float64float64float64float64float64
-0.00150.0010.539116-0.0008090.001286061288690.680007633186
-0.00150.0010.989785-0.0014850.001790046725370.978045564681
-0.0010.0020.984608-0.0009850.00222922697670.457471875897
-0.0010.0020.794663-0.0007950.00215208951570.378203831404
-0.0010.0021.159759-0.001160.0023119344580.525493621222
-0.0010.0020.572174-0.0005720.002080236305390.278644209975
-0.0010.0021.105819-0.0011060.0022853524150.505074247406
-0.0010.0021.998327-0.0019980.00282724438260.784979738416
-0.0010.0022.372919-0.0023730.003103344096380.870471520917
-0.0010.0021.574572-0.0015750.00254544239440.666940281126
..................
3.57550.0050.0035050.012530.01349274692851.19117426594
3.5780.0160.0036390.013020.02062836168480.683081894981
3.57950.0110.0011660.0041750.01176519216370.362646215648
3.5820.0140.0034150.0122320.01859125574570.718122811461
3.58550.0070.0044150.015830.01730862056751.15444874925
3.5960.0160.0020470.007360.01761205546390.431190943163
3.5960.0060.003650.0131240.01443177484441.1420298228
3.5960.00.0029010.0104340.0104319961.57079632679
3.59850.0170.0042210.0151910.02279723398930.729204740345
3.60.0180.003320.0119520.0216067189550.586154176817

In [43]:
om=0.3
ol=0.7
xLCDM=1.0/np.sqrt(om*(1+pairdats['Z'])**3+ol)

In [44]:
xLCDM


Out[44]:
<Column name='Z' dtype='float64' length=6523551>
1.00067467016
1.00067467016
1.00044985352
1.00044985352
1.00044985352
1.00044985352
1.00044985352
1.00044985352
1.00044985352
1.00044985352
1.00044985352
1.00044985352
...
0.184401282435
0.184371769673
0.184312767301
0.184165396292
0.184077066095
0.183930002874
0.18372443672
0.183109985927
0.183109985927
0.183109985927
0.182964183113
0.182876792494

In [42]:
xLCDMdat=pairdats.Column(name='xLCDM',data=xLCDM)
pairdats.add_column(xLCDMdat)
pairdats


Out[42]:
<Table length=6523551>
ZdeltaZdeltathetaZdeltathetaRalphaxLCDM
float64float64float64float64float64float64float64
-0.00150.0010.539116-0.0008090.001286061288690.6800076331861.00067467016
-0.00150.0010.989785-0.0014850.001790046725370.9780455646811.00067467016
-0.0010.0020.984608-0.0009850.00222922697670.4574718758971.00044985352
-0.0010.0020.794663-0.0007950.00215208951570.3782038314041.00044985352
-0.0010.0021.159759-0.001160.0023119344580.5254936212221.00044985352
-0.0010.0020.572174-0.0005720.002080236305390.2786442099751.00044985352
-0.0010.0021.105819-0.0011060.0022853524150.5050742474061.00044985352
-0.0010.0021.998327-0.0019980.00282724438260.7849797384161.00044985352
-0.0010.0022.372919-0.0023730.003103344096380.8704715209171.00044985352
-0.0010.0021.574572-0.0015750.00254544239440.6669402811261.00044985352
.....................
3.57550.0050.0035050.012530.01349274692851.191174265940.184312767301
3.5780.0160.0036390.013020.02062836168480.6830818949810.184165396292
3.57950.0110.0011660.0041750.01176519216370.3626462156480.184077066095
3.5820.0140.0034150.0122320.01859125574570.7181228114610.183930002874
3.58550.0070.0044150.015830.01730862056751.154448749250.18372443672
3.5960.0160.0020470.007360.01761205546390.4311909431630.183109985927
3.5960.0060.003650.0131240.01443177484441.14202982280.183109985927
3.5960.00.0029010.0104340.0104319961.570796326790.183109985927
3.59850.0170.0042210.0151910.02279723398930.7292047403450.182964183113
3.60.0180.003320.0119520.0216067189550.5861541768170.182876792494

In [93]:
yLCDM0=np.sqrt(om*(1+pairdats['Z'])**3+ol)/pairdats['Z']


/Users/rohin/anaconda/lib/python2.7/site-packages/ipykernel/__main__.py:1: RuntimeWarning: divide by zero encountered in divide
  if __name__ == '__main__':

In [94]:
yLCDM0


Out[94]:
<Column name='Z' dtype='float64' length=6523551>
-666.217189807
-666.217189807
-999.550348757
-999.550348757
-999.550348757
-999.550348757
-999.550348757
-999.550348757
-999.550348757
-999.550348757
-999.550348757
-999.550348757
...
1.51733513647
1.517365741
1.51742696705
1.51758013098
1.51767209697
1.51782548603
1.51804046608
1.51868704299
1.51868704299
1.51868704299
1.51884134936
1.51893399917

In [95]:
LCDMf = lambda x: 1.0/math.sqrt(om*(1+x)**3+ol)
np.vectorize(LCDMf)

def LCDMfint(z):
    return integrate.quad(LCDMf, 0, z)

LCDMfint=np.vectorize(LCDMfint)

yLCDM=yLCDM0*LCDMfint(pairdats['Z'])[0]


/Users/rohin/anaconda/lib/python2.7/site-packages/ipykernel/__main__.py:9: RuntimeWarning: invalid value encountered in multiply

In [96]:
yLCDM


Out[96]:
<Column name='Z' dtype='float64' length=6523551>
0.999662947356
0.999662947356
0.9997751988
0.9997751988
0.9997751988
0.9997751988
0.9997751988
0.9997751988
0.9997751988
0.9997751988
0.9997751988
0.9997751988
...
2.4287612342
2.42895011303
2.42932784764
2.43027204961
2.43083847858
2.4317823732
2.43310350314
2.43706463723
2.43706463723
2.43706463723
2.43800726613
2.43857275155

In [98]:
plt.plot(np.log(pairdats['Z']),yLCDM)


/Users/rohin/anaconda/lib/python2.7/site-packages/ipykernel/__main__.py:1: RuntimeWarning: divide by zero encountered in log
  if __name__ == '__main__':
/Users/rohin/anaconda/lib/python2.7/site-packages/ipykernel/__main__.py:1: RuntimeWarning: invalid value encountered in log
  if __name__ == '__main__':
Out[98]:
[<matplotlib.lines.Line2D at 0x11aa02990>]

In [99]:
yLCDMdat=pairdats.Column(name='yLCDM',data=yLCDM)
pairdats.add_column(yLCDMdat)
pairdats


Out[99]:
<Table length=6523551>
ZdeltaZdeltathetaZdeltathetaRalphaxLCDMyLCDM
float64float64float64float64float64float64float64float64
-0.00150.0010.539116-0.0008090.001286061288690.6800076331861.000674670160.999662947356
-0.00150.0010.989785-0.0014850.001790046725370.9780455646811.000674670160.999662947356
-0.0010.0020.984608-0.0009850.00222922697670.4574718758971.000449853520.9997751988
-0.0010.0020.794663-0.0007950.00215208951570.3782038314041.000449853520.9997751988
-0.0010.0021.159759-0.001160.0023119344580.5254936212221.000449853520.9997751988
-0.0010.0020.572174-0.0005720.002080236305390.2786442099751.000449853520.9997751988
-0.0010.0021.105819-0.0011060.0022853524150.5050742474061.000449853520.9997751988
-0.0010.0021.998327-0.0019980.00282724438260.7849797384161.000449853520.9997751988
-0.0010.0022.372919-0.0023730.003103344096380.8704715209171.000449853520.9997751988
-0.0010.0021.574572-0.0015750.00254544239440.6669402811261.000449853520.9997751988
........................
3.57550.0050.0035050.012530.01349274692851.191174265940.1843127673012.42932784764
3.5780.0160.0036390.013020.02062836168480.6830818949810.1841653962922.43027204961
3.57950.0110.0011660.0041750.01176519216370.3626462156480.1840770660952.43083847858
3.5820.0140.0034150.0122320.01859125574570.7181228114610.1839300028742.4317823732
3.58550.0070.0044150.015830.01730862056751.154448749250.183724436722.43310350314
3.5960.0160.0020470.007360.01761205546390.4311909431630.1831099859272.43706463723
3.5960.0060.003650.0131240.01443177484441.14202982280.1831099859272.43706463723
3.5960.00.0029010.0104340.0104319961.570796326790.1831099859272.43706463723
3.59850.0170.0042210.0151910.02279723398930.7292047403450.1829641831132.43800726613
3.60.0180.003320.0119520.0216067189550.5861541768170.1828767924942.43857275155

In [100]:
sLCDM=pairdats['xLCDM']*np.sqrt((pairdats['deltaZ'])**2+(pairdats['yLCDM']*pairdats['Zdeltatheta'])**2)

In [101]:
sLCDMdat=pairdats.Column(name='sLCDM',data=sLCDM)
pairdats.add_column(sLCDMdat)
pairdats


Out[101]:
<Table length=6523551>
ZdeltaZdeltathetaZdeltathetaRalphaxLCDMyLCDMsLCDM
float64float64float64float64float64float64float64float64float64
-0.00150.0010.539116-0.0008090.001286061288690.6800076331861.000674670160.9996629473560.0012869625095
-0.00150.0010.989785-0.0014850.001790046725370.9780455646811.000674670160.9996629473560.00179110666594
-0.0010.0020.984608-0.0009850.00222922697670.4574718758971.000449853520.99977519880.0022303051794
-0.0010.0020.794663-0.0007950.00215208951570.3782038314041.000449853520.99977519880.00215311611837
-0.0010.0021.159759-0.001160.0023119344580.5254936212221.000449853520.99977519880.00231296456834
-0.0010.0020.572174-0.0005720.002080236305390.2786442099751.000449853520.99977519880.00208108886284
-0.0010.0021.105819-0.0011060.0022853524150.5050742474061.000449853520.99977519880.00228634775052
-0.0010.0021.998327-0.0019980.00282724438260.7849797384161.000449853520.99977519880.00282796744178
-0.0010.0022.372919-0.0023730.003103344096380.8704715209171.000449853520.99977519880.00310439404409
-0.0010.0021.574572-0.0015750.00254544239440.6669402811261.000449853520.99977519880.00254663322885
...........................
3.57550.0050.0035050.012530.01349274692851.191174265940.1843127673012.429327847640.00568556885124
3.5780.0160.0036390.013020.02062836168480.6830818949810.1841653962922.430272049610.00653002077534
3.57950.0110.0011660.0041750.01176519216370.3626462156480.1840770660952.430838478580.00275499566675
3.5820.0140.0034150.0122320.01859125574570.7181228114610.1839300028742.43178237320.00604679069142
3.58550.0070.0044150.015830.01730862056751.154448749250.183724436722.433103503140.00719225310651
3.5960.0160.0020470.007360.01761205546390.4311909431630.1831099859272.437064637230.00440122912702
3.5960.0060.003650.0131240.01443177484441.14202982280.1831099859272.437064637230.00595875619834
3.5960.00.0029010.0104340.0104319961.570796326790.1831099859272.437064637230.00465618159246
3.59850.0170.0042210.0151910.02279723398930.7292047403450.1829641831132.438007266130.00745598271651
3.60.0180.003320.0119520.0216067189550.5861541768170.1828767924942.438572751550.00626464175486

In [104]:
ascii.write(pairdats,"./output/sLCDM_correl.dat")

In [105]:
pairdats


Out[105]:
<Table length=6523551>
ZdeltaZdeltathetaZdeltathetaRalphaxLCDMyLCDMsLCDM
float64float64float64float64float64float64float64float64float64
-0.00150.0010.539116-0.0008090.001286061288690.6800076331861.000674670160.9996629473560.0012869625095
-0.00150.0010.989785-0.0014850.001790046725370.9780455646811.000674670160.9996629473560.00179110666594
-0.0010.0020.984608-0.0009850.00222922697670.4574718758971.000449853520.99977519880.0022303051794
-0.0010.0020.794663-0.0007950.00215208951570.3782038314041.000449853520.99977519880.00215311611837
-0.0010.0021.159759-0.001160.0023119344580.5254936212221.000449853520.99977519880.00231296456834
-0.0010.0020.572174-0.0005720.002080236305390.2786442099751.000449853520.99977519880.00208108886284
-0.0010.0021.105819-0.0011060.0022853524150.5050742474061.000449853520.99977519880.00228634775052
-0.0010.0021.998327-0.0019980.00282724438260.7849797384161.000449853520.99977519880.00282796744178
-0.0010.0022.372919-0.0023730.003103344096380.8704715209171.000449853520.99977519880.00310439404409
-0.0010.0021.574572-0.0015750.00254544239440.6669402811261.000449853520.99977519880.00254663322885
...........................
3.57550.0050.0035050.012530.01349274692851.191174265940.1843127673012.429327847640.00568556885124
3.5780.0160.0036390.013020.02062836168480.6830818949810.1841653962922.430272049610.00653002077534
3.57950.0110.0011660.0041750.01176519216370.3626462156480.1840770660952.430838478580.00275499566675
3.5820.0140.0034150.0122320.01859125574570.7181228114610.1839300028742.43178237320.00604679069142
3.58550.0070.0044150.015830.01730862056751.154448749250.183724436722.433103503140.00719225310651
3.5960.0160.0020470.007360.01761205546390.4311909431630.1831099859272.437064637230.00440122912702
3.5960.0060.003650.0131240.01443177484441.14202982280.1831099859272.437064637230.00595875619834
3.5960.00.0029010.0104340.0104319961.570796326790.1831099859272.437064637230.00465618159246
3.59850.0170.0042210.0151910.02279723398930.7292047403450.1829641831132.438007266130.00745598271651
3.60.0180.003320.0119520.0216067189550.5861541768170.1828767924942.438572751550.00626464175486

In [ ]:
spairs_knuth=hist(pairdats['sLCDM'][2253:len(pairdats['sLCDM']-1)],bins='knuth')
spairs_blocks=hist(pairdats['sLCDM'][2253:len(pairdats['sLCDM']-1)],bins='blocks')
DDsk=AutoCorrelation(spairs_knuth[0])
DDsb=AutoCorrelation(spairs_blocks[0])
plt.plot(spairs_knuth[1][1:len(spairs_knuth[1])],DDsk)
plt.savefig("./images/DD_sLCDM_knuth.pdf")
plt.savefig("./images/DD_sLCDM_knuth.jpeg")
plt.plot(spairs_blocks[1][1:len(spairs_blocks[1])],DDsb)
plt.savefig("./images/DD_sLCDM_blocks.pdf")
plt.savefig("./images/DD_sLCDM_blocks.jpeg")

In [129]:
spairs_1000=hist(sLCDM[2253:len(sLCDM)-1],bins=1000)
plt.savefig("./images/hist_sLCDM_1000.pdf")
plt.savefig("./images/hist_sLCDM_1000.jpeg")
DDs1000=AutoCorrelation(spairs_1000[0])
plt.figure()
plt.plot(spairs_1000[1][1:len(spairs_1000[1])],DDs1000)
plt.savefig("./images/DD_sLCDM_1000.pdf")
plt.savefig("./images/DD_sLCDM_1000.jpeg")



In [112]:
max(sLCDM)


Out[112]:
0.028187966602223606

In [114]:
np.all(np.isfinite(sLCDM))


Out[114]:
False

In [115]:
min(sLCDM)


Out[115]:
9.6817160036876236e-07

In [116]:
help(np.isfinite)


Help on ufunc object:

isfinite = class ufunc(__builtin__.object)
 |  Functions that operate element by element on whole arrays.
 |  
 |  To see the documentation for a specific ufunc, use np.info().  For
 |  example, np.info(np.sin).  Because ufuncs are written in C
 |  (for speed) and linked into Python with NumPy's ufunc facility,
 |  Python's help() function finds this page whenever help() is called
 |  on a ufunc.
 |  
 |  A detailed explanation of ufuncs can be found in the "ufuncs.rst"
 |  file in the NumPy reference guide.
 |  
 |  Unary ufuncs:
 |  =============
 |  
 |  op(X, out=None)
 |  Apply op to X elementwise
 |  
 |  Parameters
 |  ----------
 |  X : array_like
 |      Input array.
 |  out : array_like
 |      An array to store the output. Must be the same shape as `X`.
 |  
 |  Returns
 |  -------
 |  r : array_like
 |      `r` will have the same shape as `X`; if out is provided, `r`
 |      will be equal to out.
 |  
 |  Binary ufuncs:
 |  ==============
 |  
 |  op(X, Y, out=None)
 |  Apply `op` to `X` and `Y` elementwise. May "broadcast" to make
 |  the shapes of `X` and `Y` congruent.
 |  
 |  The broadcasting rules are:
 |  
 |  * Dimensions of length 1 may be prepended to either array.
 |  * Arrays may be repeated along dimensions of length 1.
 |  
 |  Parameters
 |  ----------
 |  X : array_like
 |      First input array.
 |  Y : array_like
 |      Second input array.
 |  out : array_like
 |      An array to store the output. Must be the same shape as the
 |      output would have.
 |  
 |  Returns
 |  -------
 |  r : array_like
 |      The return value; if out is provided, `r` will be equal to out.
 |  
 |  Methods defined here:
 |  
 |  __call__(...)
 |      x.__call__(...) <==> x(...)
 |  
 |  __repr__(...)
 |      x.__repr__() <==> repr(x)
 |  
 |  __str__(...)
 |      x.__str__() <==> str(x)
 |  
 |  accumulate(...)
 |      accumulate(array, axis=0, dtype=None, out=None)
 |      
 |      Accumulate the result of applying the operator to all elements.
 |      
 |      For a one-dimensional array, accumulate produces results equivalent to::
 |      
 |        r = np.empty(len(A))
 |        t = op.identity        # op = the ufunc being applied to A's  elements
 |        for i in range(len(A)):
 |            t = op(t, A[i])
 |            r[i] = t
 |        return r
 |      
 |      For example, add.accumulate() is equivalent to np.cumsum().
 |      
 |      For a multi-dimensional array, accumulate is applied along only one
 |      axis (axis zero by default; see Examples below) so repeated use is
 |      necessary if one wants to accumulate over multiple axes.
 |      
 |      Parameters
 |      ----------
 |      array : array_like
 |          The array to act on.
 |      axis : int, optional
 |          The axis along which to apply the accumulation; default is zero.
 |      dtype : data-type code, optional
 |          The data-type used to represent the intermediate results. Defaults
 |          to the data-type of the output array if such is provided, or the
 |          the data-type of the input array if no output array is provided.
 |      out : ndarray, optional
 |          A location into which the result is stored. If not provided a
 |          freshly-allocated array is returned.
 |      
 |      Returns
 |      -------
 |      r : ndarray
 |          The accumulated values. If `out` was supplied, `r` is a reference to
 |          `out`.
 |      
 |      Examples
 |      --------
 |      1-D array examples:
 |      
 |      >>> np.add.accumulate([2, 3, 5])
 |      array([ 2,  5, 10])
 |      >>> np.multiply.accumulate([2, 3, 5])
 |      array([ 2,  6, 30])
 |      
 |      2-D array examples:
 |      
 |      >>> I = np.eye(2)
 |      >>> I
 |      array([[ 1.,  0.],
 |             [ 0.,  1.]])
 |      
 |      Accumulate along axis 0 (rows), down columns:
 |      
 |      >>> np.add.accumulate(I, 0)
 |      array([[ 1.,  0.],
 |             [ 1.,  1.]])
 |      >>> np.add.accumulate(I) # no axis specified = axis zero
 |      array([[ 1.,  0.],
 |             [ 1.,  1.]])
 |      
 |      Accumulate along axis 1 (columns), through rows:
 |      
 |      >>> np.add.accumulate(I, 1)
 |      array([[ 1.,  1.],
 |             [ 0.,  1.]])
 |  
 |  at(...)
 |      at(a, indices, b=None)
 |      
 |      Performs unbuffered in place operation on operand 'a' for elements
 |      specified by 'indices'. For addition ufunc, this method is equivalent to
 |      `a[indices] += b`, except that results are accumulated for elements that
 |      are indexed more than once. For example, `a[[0,0]] += 1` will only
 |      increment the first element once because of buffering, whereas
 |      `add.at(a, [0,0], 1)` will increment the first element twice.
 |      
 |      .. versionadded:: 1.8.0
 |      
 |      Parameters
 |      ----------
 |      a : array_like
 |          The array to perform in place operation on.
 |      indices : array_like or tuple
 |          Array like index object or slice object for indexing into first
 |          operand. If first operand has multiple dimensions, indices can be a
 |          tuple of array like index objects or slice objects.
 |      b : array_like
 |          Second operand for ufuncs requiring two operands. Operand must be
 |          broadcastable over first operand after indexing or slicing.
 |      
 |      Examples
 |      --------
 |      Set items 0 and 1 to their negative values:
 |      
 |      >>> a = np.array([1, 2, 3, 4])
 |      >>> np.negative.at(a, [0, 1])
 |      >>> print(a)
 |      array([-1, -2, 3, 4])
 |      
 |      ::
 |      
 |      Increment items 0 and 1, and increment item 2 twice:
 |      
 |      >>> a = np.array([1, 2, 3, 4])
 |      >>> np.add.at(a, [0, 1, 2, 2], 1)
 |      >>> print(a)
 |      array([2, 3, 5, 4])
 |      
 |      ::
 |      
 |      Add items 0 and 1 in first array to second array,
 |      and store results in first array:
 |      
 |      >>> a = np.array([1, 2, 3, 4])
 |      >>> b = np.array([1, 2])
 |      >>> np.add.at(a, [0, 1], b)
 |      >>> print(a)
 |      array([2, 4, 3, 4])
 |  
 |  outer(...)
 |      outer(A, B)
 |      
 |      Apply the ufunc `op` to all pairs (a, b) with a in `A` and b in `B`.
 |      
 |      Let ``M = A.ndim``, ``N = B.ndim``. Then the result, `C`, of
 |      ``op.outer(A, B)`` is an array of dimension M + N such that:
 |      
 |      .. math:: C[i_0, ..., i_{M-1}, j_0, ..., j_{N-1}] =
 |         op(A[i_0, ..., i_{M-1}], B[j_0, ..., j_{N-1}])
 |      
 |      For `A` and `B` one-dimensional, this is equivalent to::
 |      
 |        r = empty(len(A),len(B))
 |        for i in range(len(A)):
 |            for j in range(len(B)):
 |                r[i,j] = op(A[i], B[j]) # op = ufunc in question
 |      
 |      Parameters
 |      ----------
 |      A : array_like
 |          First array
 |      B : array_like
 |          Second array
 |      
 |      Returns
 |      -------
 |      r : ndarray
 |          Output array
 |      
 |      See Also
 |      --------
 |      numpy.outer
 |      
 |      Examples
 |      --------
 |      >>> np.multiply.outer([1, 2, 3], [4, 5, 6])
 |      array([[ 4,  5,  6],
 |             [ 8, 10, 12],
 |             [12, 15, 18]])
 |      
 |      A multi-dimensional example:
 |      
 |      >>> A = np.array([[1, 2, 3], [4, 5, 6]])
 |      >>> A.shape
 |      (2, 3)
 |      >>> B = np.array([[1, 2, 3, 4]])
 |      >>> B.shape
 |      (1, 4)
 |      >>> C = np.multiply.outer(A, B)
 |      >>> C.shape; C
 |      (2, 3, 1, 4)
 |      array([[[[ 1,  2,  3,  4]],
 |              [[ 2,  4,  6,  8]],
 |              [[ 3,  6,  9, 12]]],
 |             [[[ 4,  8, 12, 16]],
 |              [[ 5, 10, 15, 20]],
 |              [[ 6, 12, 18, 24]]]])
 |  
 |  reduce(...)
 |      reduce(a, axis=0, dtype=None, out=None, keepdims=False)
 |      
 |      Reduces `a`'s dimension by one, by applying ufunc along one axis.
 |      
 |      Let :math:`a.shape = (N_0, ..., N_i, ..., N_{M-1})`.  Then
 |      :math:`ufunc.reduce(a, axis=i)[k_0, ..,k_{i-1}, k_{i+1}, .., k_{M-1}]` =
 |      the result of iterating `j` over :math:`range(N_i)`, cumulatively applying
 |      ufunc to each :math:`a[k_0, ..,k_{i-1}, j, k_{i+1}, .., k_{M-1}]`.
 |      For a one-dimensional array, reduce produces results equivalent to:
 |      ::
 |      
 |       r = op.identity # op = ufunc
 |       for i in range(len(A)):
 |         r = op(r, A[i])
 |       return r
 |      
 |      For example, add.reduce() is equivalent to sum().
 |      
 |      Parameters
 |      ----------
 |      a : array_like
 |          The array to act on.
 |      axis : None or int or tuple of ints, optional
 |          Axis or axes along which a reduction is performed.
 |          The default (`axis` = 0) is perform a reduction over the first
 |          dimension of the input array. `axis` may be negative, in
 |          which case it counts from the last to the first axis.
 |      
 |          .. versionadded:: 1.7.0
 |      
 |          If this is `None`, a reduction is performed over all the axes.
 |          If this is a tuple of ints, a reduction is performed on multiple
 |          axes, instead of a single axis or all the axes as before.
 |      
 |          For operations which are either not commutative or not associative,
 |          doing a reduction over multiple axes is not well-defined. The
 |          ufuncs do not currently raise an exception in this case, but will
 |          likely do so in the future.
 |      dtype : data-type code, optional
 |          The type used to represent the intermediate results. Defaults
 |          to the data-type of the output array if this is provided, or
 |          the data-type of the input array if no output array is provided.
 |      out : ndarray, optional
 |          A location into which the result is stored. If not provided, a
 |          freshly-allocated array is returned.
 |      keepdims : bool, optional
 |          If this is set to True, the axes which are reduced are left
 |          in the result as dimensions with size one. With this option,
 |          the result will broadcast correctly against the original `arr`.
 |      
 |          .. versionadded:: 1.7.0
 |      
 |      Returns
 |      -------
 |      r : ndarray
 |          The reduced array. If `out` was supplied, `r` is a reference to it.
 |      
 |      Examples
 |      --------
 |      >>> np.multiply.reduce([2,3,5])
 |      30
 |      
 |      A multi-dimensional array example:
 |      
 |      >>> X = np.arange(8).reshape((2,2,2))
 |      >>> X
 |      array([[[0, 1],
 |              [2, 3]],
 |             [[4, 5],
 |              [6, 7]]])
 |      >>> np.add.reduce(X, 0)
 |      array([[ 4,  6],
 |             [ 8, 10]])
 |      >>> np.add.reduce(X) # confirm: default axis value is 0
 |      array([[ 4,  6],
 |             [ 8, 10]])
 |      >>> np.add.reduce(X, 1)
 |      array([[ 2,  4],
 |             [10, 12]])
 |      >>> np.add.reduce(X, 2)
 |      array([[ 1,  5],
 |             [ 9, 13]])
 |  
 |  reduceat(...)
 |      reduceat(a, indices, axis=0, dtype=None, out=None)
 |      
 |      Performs a (local) reduce with specified slices over a single axis.
 |      
 |      For i in ``range(len(indices))``, `reduceat` computes
 |      ``ufunc.reduce(a[indices[i]:indices[i+1]])``, which becomes the i-th
 |      generalized "row" parallel to `axis` in the final result (i.e., in a
 |      2-D array, for example, if `axis = 0`, it becomes the i-th row, but if
 |      `axis = 1`, it becomes the i-th column).  There are three exceptions to this:
 |      
 |      * when ``i = len(indices) - 1`` (so for the last index),
 |        ``indices[i+1] = a.shape[axis]``.
 |      * if ``indices[i] >= indices[i + 1]``, the i-th generalized "row" is
 |        simply ``a[indices[i]]``.
 |      * if ``indices[i] >= len(a)`` or ``indices[i] < 0``, an error is raised.
 |      
 |      The shape of the output depends on the size of `indices`, and may be
 |      larger than `a` (this happens if ``len(indices) > a.shape[axis]``).
 |      
 |      Parameters
 |      ----------
 |      a : array_like
 |          The array to act on.
 |      indices : array_like
 |          Paired indices, comma separated (not colon), specifying slices to
 |          reduce.
 |      axis : int, optional
 |          The axis along which to apply the reduceat.
 |      dtype : data-type code, optional
 |          The type used to represent the intermediate results. Defaults
 |          to the data type of the output array if this is provided, or
 |          the data type of the input array if no output array is provided.
 |      out : ndarray, optional
 |          A location into which the result is stored. If not provided a
 |          freshly-allocated array is returned.
 |      
 |      Returns
 |      -------
 |      r : ndarray
 |          The reduced values. If `out` was supplied, `r` is a reference to
 |          `out`.
 |      
 |      Notes
 |      -----
 |      A descriptive example:
 |      
 |      If `a` is 1-D, the function `ufunc.accumulate(a)` is the same as
 |      ``ufunc.reduceat(a, indices)[::2]`` where `indices` is
 |      ``range(len(array) - 1)`` with a zero placed
 |      in every other element:
 |      ``indices = zeros(2 * len(a) - 1)``, ``indices[1::2] = range(1, len(a))``.
 |      
 |      Don't be fooled by this attribute's name: `reduceat(a)` is not
 |      necessarily smaller than `a`.
 |      
 |      Examples
 |      --------
 |      To take the running sum of four successive values:
 |      
 |      >>> np.add.reduceat(np.arange(8),[0,4, 1,5, 2,6, 3,7])[::2]
 |      array([ 6, 10, 14, 18])
 |      
 |      A 2-D example:
 |      
 |      >>> x = np.linspace(0, 15, 16).reshape(4,4)
 |      >>> x
 |      array([[  0.,   1.,   2.,   3.],
 |             [  4.,   5.,   6.,   7.],
 |             [  8.,   9.,  10.,  11.],
 |             [ 12.,  13.,  14.,  15.]])
 |      
 |      ::
 |      
 |       # reduce such that the result has the following five rows:
 |       # [row1 + row2 + row3]
 |       # [row4]
 |       # [row2]
 |       # [row3]
 |       # [row1 + row2 + row3 + row4]
 |      
 |      >>> np.add.reduceat(x, [0, 3, 1, 2, 0])
 |      array([[ 12.,  15.,  18.,  21.],
 |             [ 12.,  13.,  14.,  15.],
 |             [  4.,   5.,   6.,   7.],
 |             [  8.,   9.,  10.,  11.],
 |             [ 24.,  28.,  32.,  36.]])
 |      
 |      ::
 |      
 |       # reduce such that result has the following two columns:
 |       # [col1 * col2 * col3, col4]
 |      
 |      >>> np.multiply.reduceat(x, [0, 3], 1)
 |      array([[    0.,     3.],
 |             [  120.,     7.],
 |             [  720.,    11.],
 |             [ 2184.,    15.]])
 |  
 |  ----------------------------------------------------------------------
 |  Data descriptors defined here:
 |  
 |  identity
 |      The identity value.
 |      
 |      Data attribute containing the identity element for the ufunc, if it has one.
 |      If it does not, the attribute value is None.
 |      
 |      Examples
 |      --------
 |      >>> np.add.identity
 |      0
 |      >>> np.multiply.identity
 |      1
 |      >>> np.power.identity
 |      1
 |      >>> print(np.exp.identity)
 |      None
 |  
 |  nargs
 |      The number of arguments.
 |      
 |      Data attribute containing the number of arguments the ufunc takes, including
 |      optional ones.
 |      
 |      Notes
 |      -----
 |      Typically this value will be one more than what you might expect because all
 |      ufuncs take  the optional "out" argument.
 |      
 |      Examples
 |      --------
 |      >>> np.add.nargs
 |      3
 |      >>> np.multiply.nargs
 |      3
 |      >>> np.power.nargs
 |      3
 |      >>> np.exp.nargs
 |      2
 |  
 |  nin
 |      The number of inputs.
 |      
 |      Data attribute containing the number of arguments the ufunc treats as input.
 |      
 |      Examples
 |      --------
 |      >>> np.add.nin
 |      2
 |      >>> np.multiply.nin
 |      2
 |      >>> np.power.nin
 |      2
 |      >>> np.exp.nin
 |      1
 |  
 |  nout
 |      The number of outputs.
 |      
 |      Data attribute containing the number of arguments the ufunc treats as output.
 |      
 |      Notes
 |      -----
 |      Since all ufuncs can take output arguments, this will always be (at least) 1.
 |      
 |      Examples
 |      --------
 |      >>> np.add.nout
 |      1
 |      >>> np.multiply.nout
 |      1
 |      >>> np.power.nout
 |      1
 |      >>> np.exp.nout
 |      1
 |  
 |  ntypes
 |      The number of types.
 |      
 |      The number of numerical NumPy types - of which there are 18 total - on which
 |      the ufunc can operate.
 |      
 |      See Also
 |      --------
 |      numpy.ufunc.types
 |      
 |      Examples
 |      --------
 |      >>> np.add.ntypes
 |      18
 |      >>> np.multiply.ntypes
 |      18
 |      >>> np.power.ntypes
 |      17
 |      >>> np.exp.ntypes
 |      7
 |      >>> np.remainder.ntypes
 |      14
 |  
 |  signature
 |  
 |  types
 |      Returns a list with types grouped input->output.
 |      
 |      Data attribute listing the data-type "Domain-Range" groupings the ufunc can
 |      deliver. The data-types are given using the character codes.
 |      
 |      See Also
 |      --------
 |      numpy.ufunc.ntypes
 |      
 |      Examples
 |      --------
 |      >>> np.add.types
 |      ['??->?', 'bb->b', 'BB->B', 'hh->h', 'HH->H', 'ii->i', 'II->I', 'll->l',
 |      'LL->L', 'qq->q', 'QQ->Q', 'ff->f', 'dd->d', 'gg->g', 'FF->F', 'DD->D',
 |      'GG->G', 'OO->O']
 |      
 |      >>> np.multiply.types
 |      ['??->?', 'bb->b', 'BB->B', 'hh->h', 'HH->H', 'ii->i', 'II->I', 'll->l',
 |      'LL->L', 'qq->q', 'QQ->Q', 'ff->f', 'dd->d', 'gg->g', 'FF->F', 'DD->D',
 |      'GG->G', 'OO->O']
 |      
 |      >>> np.power.types
 |      ['bb->b', 'BB->B', 'hh->h', 'HH->H', 'ii->i', 'II->I', 'll->l', 'LL->L',
 |      'qq->q', 'QQ->Q', 'ff->f', 'dd->d', 'gg->g', 'FF->F', 'DD->D', 'GG->G',
 |      'OO->O']
 |      
 |      >>> np.exp.types
 |      ['f->f', 'd->d', 'g->g', 'F->F', 'D->D', 'G->G', 'O->O']
 |      
 |      >>> np.remainder.types
 |      ['bb->b', 'BB->B', 'hh->h', 'HH->H', 'ii->i', 'II->I', 'll->l', 'LL->L',
 |      'qq->q', 'QQ->Q', 'ff->f', 'dd->d', 'gg->g', 'OO->O']


In [117]:
np.argwhere(np.isnan(sLCDM))


Out[117]:
array([[ 206],
       [ 207],
       [ 208],
       ..., 
       [2250],
       [2251],
       [2252]])

In [119]:
sLCDM[2252]


Out[119]:
nan

In [121]:
pairdats['Z'][2252]


Out[121]:
0.0

In [ ]:
spairs1000=spairs_knuth

In [134]:
plt.figure()
plt.loglog(4.1e3*spairs_1000[1][1:len(spairs_1000[1])],DDs1000)
plt.savefig("./images/DD_sLCDMlog_1000.pdf")
plt.savefig("./images/DD_sLCDMlog_1000.jpeg")



In [135]:
plt.figure()
plt.loglog(4.1e3*spairs_1000[1][1:len(spairs_1000[1])],spairs_1000[1][1:len(spairs_1000[1])]**2*DDs1000)
plt.savefig("./images/DD_sLCDMlog_1000.pdf")
plt.savefig("./images/DD_sLCDMlog_1000.jpeg")



In [1]:
pair_blocks


---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-1-ffef691b92ba> in <module>()
----> 1 pair_blocks

NameError: name 'pair_blocks' is not defined

In [3]:
help(np.digitize)


Help on built-in function digitize in module numpy.core.multiarray:

digitize(...)
    digitize(x, bins, right=False)
    
    Return the indices of the bins to which each value in input array belongs.
    
    Each index ``i`` returned is such that ``bins[i-1] <= x < bins[i]`` if
    `bins` is monotonically increasing, or ``bins[i-1] > x >= bins[i]`` if
    `bins` is monotonically decreasing. If values in `x` are beyond the
    bounds of `bins`, 0 or ``len(bins)`` is returned as appropriate. If right
    is True, then the right bin is closed so that the index ``i`` is such
    that ``bins[i-1] < x <= bins[i]`` or bins[i-1] >= x > bins[i]`` if `bins`
    is monotonically increasing or decreasing, respectively.
    
    Parameters
    ----------
    x : array_like
        Input array to be binned. Prior to Numpy 1.10.0, this array had to
        be 1-dimensional, but can now have any shape.
    bins : array_like
        Array of bins. It has to be 1-dimensional and monotonic.
    right : bool, optional
        Indicating whether the intervals include the right or the left bin
        edge. Default behavior is (right==False) indicating that the interval
        does not include the right edge. The left bin end is open in this
        case, i.e., bins[i-1] <= x < bins[i] is the default behavior for
        monotonically increasing bins.
    
    Returns
    -------
    out : ndarray of ints
        Output array of indices, of same shape as `x`.
    
    Raises
    ------
    ValueError
        If `bins` is not monotonic.
    TypeError
        If the type of the input is complex.
    
    See Also
    --------
    bincount, histogram, unique
    
    Notes
    -----
    If values in `x` are such that they fall outside the bin range,
    attempting to index `bins` with the indices that `digitize` returns
    will result in an IndexError.
    
    .. versionadded:: 1.10.0
    
    `np.digitize` is  implemented in terms of `np.searchsorted`. This means
    that a binary search is used to bin the values, which scales much better
    for larger number of bins than the previous linear search. It also removes
    the requirement for the input array to be 1-dimensional.
    
    Examples
    --------
    >>> x = np.array([0.2, 6.4, 3.0, 1.6])
    >>> bins = np.array([0.0, 1.0, 2.5, 4.0, 10.0])
    >>> inds = np.digitize(x, bins)
    >>> inds
    array([1, 4, 3, 2])
    >>> for n in range(x.size):
    ...   print(bins[inds[n]-1], "<=", x[n], "<", bins[inds[n]])
    ...
    0.0 <= 0.2 < 1.0
    4.0 <= 6.4 < 10.0
    2.5 <= 3.0 < 4.0
    1.0 <= 1.6 < 2.5
    
    >>> x = np.array([1.2, 10.0, 12.4, 15.5, 20.])
    >>> bins = np.array([0, 5, 10, 15, 20])
    >>> np.digitize(x,bins,right=True)
    array([1, 2, 3, 4, 4])
    >>> np.digitize(x,bins,right=False)
    array([1, 3, 3, 4, 5])


In [ ]: