In [1]:
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 [ ]: