Graph Analyses of the SDSS Skeleton Galaxies¶
Problem Description:¶
Measuring graph statistics of the SDSS Skeleton galaxies using python-igraph and pyspark's parallel computation capability.
Key Steps :¶
- generating cKDTree using 3-d positions of SDSS galaxies
- generating graph edges using cKDTree for a certain linking length, such 1.0 Mpc/h
- parallelize graph calculations for various
linklenusing PySpark/UDF. For this, we will broadcast two dataframes, one the galaxy cataloggaldfand the other kdtreegalkdtree
In [1]:
# generate edges
import sys
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.spatial import cKDTree
# plot settings
plt.rc('font', family='serif')
plt.rc('font', serif='Times New Roman')
plt.rcParams.update({'font.size': 16})
plt.rcParams['mathtext.fontset'] = 'stix'
In [2]:
# galaxy dataframe
galdf = pd.read_csv("./D2018Mar_SDSSvolskbfHR4.txt",delimiter=r"\s+")
galdf.describe()
Out[2]:
In [3]:
galdf.head()
Out[3]:
In [4]:
# how to "select" data using multiple criteria
# Key Point : use `(condition1) & (condition2)`
galdf[(galdf["pz"] < 20.0) & (galdf["pz"] > -20.0)].pz.values
Out[4]:
In [5]:
#
# 3d point distribution on a sphere
#
fig = plt.figure(figsize=(12,12))
#plt.axis([-1,1,-1,1])
plt.axes().set_aspect('equal', 'datalim')
plt.title(r'SDSS Skeleton Galaxies')
plt.xlabel(r'$x$')
plt.ylabel(r'$y$')
plt.scatter(galdf[np.abs(galdf.pz.values) < 20.].px.values,\
galdf[np.abs(galdf.pz.values) < 20.].py.values,s=1.5,marker=".", color='grey')
plt.show()
In [6]:
#
# using earth-map api_s for utilizing various projections
#
"""
from mpl_toolkits.basemap import Basemap
fig = plt.figure(figsize=(13,12))
m = Basemap(projection='ortho',lat_0=27,lon_0=260)
x, y = m(galdf.ra.values,galdf.dec.values)
m.scatter(x,y,s=2,marker='o',color='grey',alpha=0.5)
m.drawmeridians(np.arange(0,360,30))
m.drawparallels(np.arange(-90,90,30))
plt.title(r'SDSS Skeleton Galaxies')
#fig.savefig("visSDSS.png")
plt.show()
"""
Out[6]:
2. Building a KDtree and Generating FOF networks for various linking lengths¶
- Defining a function measuring angular separation between two polar coordinates for 2-d data
In [7]:
# define a spherical distance function
# angular distance between two polar coordinates : (phi,theta)
def polarangdist (p1,t1,p2,t2,debug=False) :
innerproduct = np.sin(np.radians(t1))*np.cos(np.radians(p1))*np.sin(np.radians(t2))*np.cos(np.radians(p2))\
+ np.sin(np.radians(t1))*np.sin(np.radians(p1))*np.sin(np.radians(t2))*np.sin(np.radians(p2))\
+ np.cos(np.radians(t1)) * np.cos(np.radians(t2))
if np.abs(innerproduct) > 1.:
if np.abs(innerproduct) < 1.001:
innerproduct = 1.0
if debug:
print "Minor roundoff issue in polarangdist : set innerproduct = 1.0"
else:
print "Funny things happen in polarangdist"
sys.exit(1)
return np.degrees(np.arccos(innerproduct))
- Making a cKDTree using 3d coordinates
In [8]:
# now generate a kdtree for finding nearest neighbor
galkdtree = cKDTree(galdf[['px','py','pz']])
sys.getsizeof(galkdtree)
Out[8]:
In [9]:
sys.getsizeof(galdf)
Out[9]:
In [10]:
# let's test how scipy/kdtree works
galkdtree.query(galdf.loc[4,['px','py','pz']],k=[2])
Out[10]:
Generating a FOF network and testing this Python/igraph with C/C++ igraph as a sanity check¶
In [11]:
#
# Generate a FOF network using the KDtree above
#
#
from igraph import *
#angcut = 0.1 # in degree
lengcut = 5.0 # in [1/h Mpc]
inearest = 1
#numtotal = len(galdf.index)
numtotal = 5
thisangdist = 0.
edgelists = []
for ivertex in range(numtotal):
inearest=2 # this is the first nearest neighbor in kdtree
dummy, inow = galkdtree.query(galdf.loc[ivertex,['px','py','pz']],k=[inearest])
"""
thisangdist= polarangdist(hotdf.loc[ivertex,['phi']].values[0],hotdf.loc[ivertex,['theta']].values[0],\
hotdf.loc[inow[0],['phi']].values[0],hotdf.loc[inow[0],['theta']].values[0],\
debug=True)"""
#while (thisangdist <= angcut):
while (dummy[0] <= lengcut):
#print ivertex, inow[0], " : ", dummy[0], inearest, " : sphdist = ", thisangdist
#print ivertex, inow[0], " : ", dummy[0], inearest
#print galdf.loc[ivertex,['ra']].values[0]
#print galdf.loc[ivertex,['ra','dec','px','py','pz']]
#print galdf.loc[inow[0],['ra']].values[0]
#print galdf.loc[inow[0],['ra','dec','px','py','pz']]
edgelists.append([ivertex,inow[0]])
inearest = inearest + 1 # the next nearest one
dummy, inow = galkdtree.query(galdf.loc[ivertex,['px','py','pz']],k=[inearest])
"""
thisangdist = polarangdist(hotdf.loc[ivertex,['phi']].values[0],hotdf.loc[ivertex,['theta']].values[0],\
hotdf.loc[inow[0],['phi']].values[0],hotdf.loc[inow[0],['theta']].values[0],\
debug=True)"""
print edgelists[:20]
It is time to parallelize the graph calculations using PySpark.¶
In [12]:
import gc
gc.collect()
Out[12]:
In [13]:
from pyspark import SparkContext
from pyspark.sql import SQLContext
#sc = SparkContext(master='local[3]', appName='calgraph')
sqlsc = SQLContext(sc)
#sc.setCheckpointDir("./checkpoints")
#sc.setCheckpointDir("hdfs://localhost:8020/myhdfs/spark/checkpoints")
sc.setCheckpointDir("hdfs://master:54310/tmp/spark/checkpoints")
In [14]:
## Define basic index and linking length dataframe
#ilink = np.arange(1,20)
#linklen = np.arange(1,20)*0.2 + 0.2
#ilink = np.arange(1,4)
#linklen = np.arange(1,4)*0.5
ilink = np.arange(1,60)
linklen = np.arange(1,60)*0.1 + 0.1
print ilink
print linklen
numdata = len(linklen)
linkdf = pd.DataFrame(list(zip(ilink, linklen)),\
columns=['ilink','linklen'])
Broadcast the kdtree to workers / Load pandasDF to sparkDF¶
In [15]:
## Now, setting for pyspark
bcastgalkdtree = sc.broadcast(galkdtree)
bcastgaldf =sc.broadcast(galdf)
sparkdf = sqlsc.createDataFrame(linkdf).persist()
sparkdf.printSchema()
In [16]:
sparkdf.show(5,truncate=True)
- Key Feature : broadcasted
objectneeds this prefix(?).valuefor accessing to its own attributes.
In [17]:
## Define an UDF for applying this to the `linklen` column
def testPassingBcastVariable(linklen, curgaldf, curtree):
dummy, inow = curtree.value.query(curgaldf.value.loc[4,['px','py','pz']],k=[2])
#dummy, inow = curtree.query([104.,179.,-136.],k=[2])
return dummy[0].item()
def temp(linklen, curgaldf, curtree):
dummy, inow = curtree.query(curgaldf.loc[4,['px','py','pz']],k=[2])
#dummy, inow = curtree.query([104.,179.,-136.],k=[2])
return dummy[0].item()
In [18]:
import pyspark.sql.functions as F
import pyspark.sql.types as T
from pyspark import Row
from functools import partial
testudf = F.udf(partial(testPassingBcastVariable, curtree = bcastgalkdtree, curgaldf = bcastgaldf),\
T.DoubleType())
#graphudf = F.udf(getGraphStatistics, T.DoubleType())
In [19]:
tmpudf = F.udf(lambda x: x*x, T.DoubleType())
In [20]:
#getGraphStatistics(0.5,galdf,galkdtree)
type(temp(0.5,galdf,galkdtree))
Out[20]:
In [21]:
galkdtree.query([104.,179.,-136.],k=[2])
Out[21]:
In [22]:
bcastgalkdtree.value.query([104.,179.,-136.],k=[2])
Out[22]:
In [23]:
bcastgaldf.value.loc[4,['px','py','pz']]
Out[23]:
In [24]:
sparkdf.withColumn('tmp',tmpudf('linklen')).show(5)
In [25]:
testdf = sparkdf.withColumn('gstat',testudf('linklen'))
testdf.persist().printSchema()
In [26]:
print len(bcastgaldf.value.index)
Now, define the right UDF and let calculations done!¶
- Key Point :
gStatsgets inputs of one column variablelinklenfrom dataframe and two broadcasted dataframescurgaldfandcurtree. Hence, got a column and apply this UDF with two passed broadcasted dataframe arguments.
- Key Note : The
numpytypes are not compatible withSpark's UDF. Hence,returnwithnumpyObject.item()by converting them to generic python types.
In [27]:
# curgaldf and curtree are Broadcasted Variables. Hence, they need ".value"
def gStats(linklen, curgaldf, curtree):
lengcut = linklen
inearest = 1
numtotal = len(curgaldf.value.index)
thisangdist = 0.
edgelists = []
g = Graph()
for ivertex in range(numtotal):
inearest=2 # this is the first nearest neighbor in kdtree
dummy, inow = curtree.value.query(curgaldf.value.loc[ivertex,['px','py','pz']],k=[inearest])
"""thisangdist= polarangdist(hotdf.loc[ivertex,['phi']].values[0],hotdf.loc[ivertex,['theta']].values[0],\
hotdf.loc[inow[0],['phi']].values[0],hotdf.loc[inow[0],['theta']].values[0],\
debug=False)"""
#while (thisangdist <= angcut):
while (dummy[0] <= lengcut):
edgelists.append([ivertex,inow[0]])
inearest = inearest + 1 # the next nearest one
dummy, inow = curtree.value.query(curgaldf.value.loc[ivertex,['px','py','pz']],k=[inearest])
"""thisangdist = polarangdist(hotdf.loc[ivertex,['phi']].values[0],hotdf.loc[ivertex,['theta']].values[0],\
hotdf.loc[inow[0],['phi']].values[0],hotdf.loc[inow[0],['theta']].values[0],\
debug=False)"""
# make a graph and do the jobs
g = Graph(edgelists)
transitivity = g.transitivity_undirected()
diameter = np.double(g.diameter()).item()
localcc = g.transitivity_avglocal_undirected()
numclique = np.double(g.clique_number()).item()
subcomp = g.components()
gcompfrac = np.double(np.max(subcomp.sizes()))/np.double(g.vcount())
numsubcomptwo = len(list(filter(lambda size: size >= 2, subcomp.sizes())))
maxdg = np.max(g.degree())
dgcent = np.double(np.sum(maxdg-g.degree()))/((g.vcount() - 2) * (g.vcount() - 1))
return [diameter,transitivity,localcc,numclique,gcompfrac.item(),\
np.double(numsubcomptwo).item(),dgcent.item()]
In [28]:
"""
%%time
tmpout = gStats(1.0,bcastgaldf,bcastgalkdtree)
"""
Out[28]:
In [29]:
# check the return types, which should be compatible for Spark/UDF's return types
#type(tmpout[6])
In [30]:
#tmpout
In [31]:
gstatudf = F.udf(partial(gStats, curtree = bcastgalkdtree, curgaldf = bcastgaldf),\
T.ArrayType(T.DoubleType()))
In [32]:
sparkdf.show()
In [33]:
# Let's do the jobs on workers
%time
gstatdf = sparkdf.withColumn('gstat',gstatudf('linklen')).cache()
In [34]:
%time gstatdf.collect()
Out[34]:
In [35]:
#show each DoubleType in the `gstat` ArrayType
%time gstatdf.select(F.col("gstat").getItem(1)).show()
3. Results¶
- Now, collect the results
gstatdf(spark/dataframe) togstatpd(pandas/dataframe)
In [40]:
gstatpd = gstatdf.select("ilink","linklen").toPandas()
In [41]:
gstatpd.head()
Out[41]:
In [42]:
tmppd = gstatdf.select(F.col("gstat").getItem(1).alias("transitivity")).toPandas()
In [43]:
re = pd.merge(gstatpd,\
gstatdf.select(F.col("gstat").getItem(0).alias("diameter")).toPandas(),\
left_index=True,right_index=True)
re = pd.merge(re,\
gstatdf.select(F.col("gstat").getItem(1).alias("transitivity")).toPandas(),\
left_index=True,right_index=True)
re = pd.merge(re,\
gstatdf.select(F.col("gstat").getItem(2).alias("localcc")).toPandas(),\
left_index=True,right_index=True)
re = pd.merge(re,\
gstatdf.select(F.col("gstat").getItem(3).alias("numclique")).toPandas(),\
left_index=True,right_index=True)
re = pd.merge(re,\
gstatdf.select(F.col("gstat").getItem(4).alias("gcomp")).toPandas(),\
left_index=True,right_index=True)
re = pd.merge(re,\
gstatdf.select(F.col("gstat").getItem(5).alias("subcomp")).toPandas(),\
left_index=True,right_index=True)
re = pd.merge(re,\
gstatdf.select(F.col("gstat").getItem(6).alias("dgcent")).toPandas(),\
left_index=True,right_index=True)
re.head()
Out[43]:
In [44]:
re.describe()
Out[44]:
In [45]:
import pickle
with open('SDSSGraphSpark.pickle','wb') as f:
pickle.dump(re,f)
f.close() #keep dumping the current results to overwrite the pickle.
- Now plot the results
In [46]:
import matplotlib.pyplot as plt
# plot settings
plt.rc('font', family='serif')
plt.rc('font', serif='Times New Roman')
plt.rcParams.update({'font.size': 16})
plt.rcParams['mathtext.fontset'] = 'stix'
lenmax = 6.0
fig = plt.figure(figsize=(12,11))
plt.subplot(321)
plt.axis([0,lenmax,0,1.0])
#plt.title("Diameter")
#plt.yscale('log')
#plt.xscale('log')
plt.xlabel(r'Linking Length [Degree]')
plt.ylabel(r'Giant Component Fraction')
plt.scatter(re.linklen.values,re.gcomp.values,s=10)
plt.plot(re.linklen.values,re.gcomp.values)
#plt.text(0.2,0.4,'CMB COLD Spots',color='b')
plt.subplot(322)
plt.axis([0,lenmax,0,500])
#plt.title("Diameter")
#plt.yscale('log')
#plt.xscale('log')
plt.xlabel(r'Linking Length [Degree]')
plt.ylabel(r'Diameter')
plt.scatter(re.linklen.values,re.diameter.values,s=10)
plt.plot(re.linklen.values,re.diameter.values)
plt.subplot(323)
plt.axis([0,lenmax,0.5,0.7])
#plt.title("Diameter")
#plt.yscale('log')
#plt.xscale('log')
plt.xlabel(r'Linking Length [Degree]')
plt.ylabel(r'Triangular Configurations')
plt.plot(re.linklen.values,re.transitivity.values,label='transitivity')
plt.scatter(re.linklen.values,re.transitivity.values,s=10)
plt.plot(re.linklen.values,re.localcc.values,color='green',label='average lcc')
plt.scatter(re.linklen.values,re.localcc.values,color='green',s=10)
plt.legend(loc=4)
plt.subplot(324)
plt.axis([0,lenmax,0,40])
#plt.title("Diameter")
#plt.yscale('log')
#plt.xscale('log')
plt.xlabel(r'Linking Length [Degree]')
plt.ylabel(r'Clique Number')
plt.scatter(re.linklen.values,re.numclique.values,s=10)
plt.plot(re.linklen.values,re.numclique.values)
plt.subplot(325)
plt.axis([0,lenmax,0,25000])
#plt.title("Diameter")
#plt.yscale('log')
#plt.xscale('log')
plt.xlabel(r'Linking Length [Degree]')
plt.ylabel(r'Number of Subcomponents')
plt.scatter(re.linklen.values,re.subcomp.values,s=10)
plt.plot(re.linklen.values,re.subcomp.values)
plt.subplot(326)
plt.axis([0,lenmax,0,0.0012])
#plt.title("Diameter")
#plt.yscale('log')
#plt.xscale('log')
plt.xlabel(r'Linking Length [Degree]')
plt.ylabel(r'Degree Centralization')
plt.scatter(re.linklen.values,re.dgcent.values,s=10)
plt.plot(re.linklen.values,re.dgcent.values)
plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0)
fig.savefig("graphSparkSDSS.eps")
plt.show()
