In [11]:
import numpy as np
dat = np.load('./data/HeCS_pared_small.npy')
# for n in range(len(dat)):
# print('')
# Ngal = dat['Ngal'][n]
# print('cluster index:', n)
# print('velocity standard deviation:', np.std(dat['vlos'][n][0:Ngal]))
# print ('cluster M500 from SZ:', dat['MSZ'][n])
# print ('cluster M200c:', dat['Mtot'][n])
# print (' ')
In [2]:
import pandas as pd
datasum = pd.DataFrame(columns = ['velsigma','M500_SZ','M200c'], index = range(len(dat)))
for n in range(len(dat)):
Ngal = dat['Ngal'][n]
velsigma = np.std(dat['vlos'][n][0:Ngal])
m500 = dat['MSZ'][n]
m200c = dat['Mtot'][n]
datasum.loc[n] = pd.Series({'velsigma':velsigma ,
'M500_SZ' : m500 ,
'M200c' : m200c })
In [3]:
## PURE POWER LAW
alpha = 0.382
sig15 = 1244
## CONTAMINATED POWER LAW
# alpha = 0.359
# sig15 = 753
M_pow = (datasum['velsigma']/sig15)**(1/alpha) * 1.e15
M_pow.name = 'M_powlaw'
ep_m200c = (M_pow - datasum['M200c'])/datasum['M200c']
ep_m200c.name = 'ep_M200c'
datacalc = pd.concat([M_pow,ep_m200c], axis=1)
In [4]:
datares = datasum.join(datacalc)
datares.head(15)
Out[4]:
In [5]:
import matplotlib as mpl
% matplotlib inline
mpl.rcParams['figure.figsize'] = [10.,10.]
plt = mpl.pyplot
In [6]:
def plothist(X,Y, n=10):
increment = (X.max() - X.min())/n
i = X.min()
y = []
err = []
while i < X.max():
y.append(Y[(X>=i) & (X<(i+increment))].mean())
err.append(Y[(X>=i) & (X<(i+increment))].std())
i += increment
plt.errorbar(np.array(range(n))*increment + X.min() + increment/2.,
y,
err,
)
In [13]:
def binnedplot(X,Y, n=10, percentiles = [15,50,85], mean=True, label='', ax = None):
if n=='auto':
n = len(X)/1000
increment = (X.max() - X.min())/n
i = X.min()
x = []
y = []
y_percent = np.ndarray((0,len(percentiles)))
err = []
while i < X.max():
bindata = Y[(X>=i) & (X<(i+increment))]
i+=increment
if len(bindata) == 0: continue
x.append(i-increment/2)
y.append(bindata.mean())
y_p = np.percentile(bindata,percentiles)
y_percent = np.append(y_percent, [y_p],axis=0)
err.append(bindata.std())
if mean:
if ax is None: plt.plot(x,y, label=label+'mean')
else: ax.plot(x,y, label=label+'mean')
y_percent = np.swapaxes(y_percent,0,1)
for i in range(len(percentiles)):
if ax is None:
plt.plot(x, y_percent[i],label=label+str(percentiles[i]))
else:
ax.plot(x, y_percent[i],label=label+str(percentiles[i]))
return x,y,err
In [14]:
X = np.log10(datares['M200c'].astype('float64'))
Y = np.log10(datares['M_powlaw'].astype('float64'))
plt.plot(X,Y,'.')
plt.xlabel('log[ M200c ]', fontsize=20)
plt.ylabel('log[ M_powlaw ]', fontsize=20)
binnedplot(X,Y,10)
Out[14]:
In [8]:
X = np.log10(datares['M200c'].astype('float64'))
Y = (datares['ep_M200c'].astype('float64'))
plt.plot(X,Y,'.')
plt.xlabel('log[ M200c ]', fontsize=20)
plt.ylabel(r'$\epsilon_{M200c}$', fontsize=20)
plothist(X,Y,10)
In [9]:
n = 20
Ngal = dat['Ngal'][n]
vlos = dat['vlos'][n][0:Ngal]
h = plt.hist(np.absolute(vlos), 12)
plt.xlabel(r'$v_{los}$', fontsize=20)
plt.ylabel(r'$dN/dv$', fontsize=20)
print('cluster index:', n)
print('velocity standard deviation:', np.std(dat['vlos'][n][0:Ngal]))
print ('cluster M500 from SZ:', dat['MSZ'][n])
print ('cluster M200c:', dat['Mtot'][n])
In [10]:
type(dat)
for n in range(len(dat)):
Ngal = dat['Ngal'][n]
velsigma = np.std(dat['vlos'][n][0:Ngal])
m500 = dat['MSZ'][n]
m200c = dat['Mtot'][n]
datasum.loc[n] = pd.Series({'velsigma':velsigma ,
'M500_SZ' : m500 ,
'M200c' : m200c })
In [11]:
dat = np.load('./data/HeCS_pared_small.npy')
featuresList = []
massList = []
numHalos = len(dat)
#loop through all of the halos
for h in range(numHalos):
numSubs = dat['Ngal'][h]
#make a numpy array that will hold, in this case, 3 features for each of the numSubs subhalos. Be aware that SDM requires at least 2 dimensions, so if you want to use line of sight velocity, you'll have to fill in a dummy value for the 2nd dimension. I used an array of 1.0's
subfeat = np.zeros((numSubs, 2), 'f')
#loop through the subhalos/galaxies in this halo:
for sh in range (0, numSubs):
#and fill in each of the features
subfeat[sh][0] = np.abs(dat['vlos'][h][sh])
subfeat[sh][1] = 1.0
#append the numpy array of features to the featuresList
featuresList.append(subfeat)
#and the mass to the massList
massList.append(dat['Mtot'])
#turn this into a features object:
# features = sdm.Features(featuresList, mass=massList)
#model.crossvalidate is the workhorse. It will do the 10-fold crossvalidation to predict answers for all of your data. It outputs root mean squared error (rmse) and mass predictions:
# rmse, predictions = model.crossvalidate(trainData, trainData.mass)
In [19]:
tempres = np.load('./tempresults.npy', encoding='latin1').item()
tempres.keys()
Out[19]:
In [20]:
X = np.log10(tempres['mass']).astype('float64')
Y = np.log10(tempres['predictions']).astype('float64')
plt.plot(X,Y,'.')
plt.xlabel('log[ M200c ]', fontsize=20)
plt.ylabel('log[ M_SDM ]', fontsize=20)
plothist(X,Y,10)
In [14]:
X.astype
Out[14]:
In [4]:
x = np.array([4,5,6,67])
x = np.append(x,124)
x[x<15]
Out[4]:
In [ ]:
In [10]:
comadat = np.load('./data/coma.p')[0]
In [ ]:
vlos = []
for i in comadat:
vlos.append(i[0])
h = plt.hist(vlos, 15)
plt.xlabel(r'$v_{los}$', fontsize=20)
plt.ylabel(r'$dN/dv$', fontsize=20)
velsigma_coma = np.std(vlos)
M_pow_coma = (velsigma_coma/sig15)**(1/alpha) * 1.e15
print ('cluster: COMA')
print ('velocity standard deviation:', np.std(vlos))
print ('mass from power law:', "%.4g" % M_pow_coma)
In [ ]: