In [1]:
%matplotlib inline
import astropy.stats
import matplotlib.pyplot as plt
from matplotlib import ticker
from matplotlib.ticker import FormatStrFormatter
import numpy as np
import os
import pandas as pd
from scipy import spatial
In [2]:
#! /usr/bin/env python
# Victor Calderon
# June 29, 2016
# Vanderbilt University
from __future__ import print_function, division, absolute_import
__author__ =['Victor Calderon']
__copyright__ =["Copyright 2016 Victor Calderon, Index function"]
__email__ =['victor.calderon@vanderbilt.edu']
__maintainer__ =['Victor Calderon']
import glob
"""
"""
def Index(directory, datatype):
"""
Indexes the files in a directory `directory' with a
specific data type.
Parameters
----------
directory: str
Absolute path to the folder that is indexed.
datatype: str
Data type of the files to be indexed in the folder.
Returns
-------
file_array: array_like
num.array of indexed files in the folder 'directory'
with specific datatype.
Examples
--------
>>> Index('~/data', '.txt')
>>> array(['A.txt', 'Z'.txt', ...])
"""
assert(os.path.exists(directory))
files = np.array(glob.glob('{0}/*{1}'.format(directory, datatype)))
return files
In [3]:
def sph_to_cart(ra,dec,cz):
"""
Converts spherical coordinates to Cartesian coordinates.
Parameters
----------
ra: array-like
right-ascension of galaxies in degrees
dec: array-like
declination of galaxies in degrees
cz: array-like
velocity of galaxies in km/s
Returns
-------
coords: array-like, shape = N by 3
x, y, and z coordinates
"""
cz_dist = cz/70. #converts velocity into distance
x_arr = cz_dist*np.cos(np.radians(ra))*np.cos(np.radians(dec))
y_arr = cz_dist*np.sin(np.radians(ra))*np.cos(np.radians(dec))
z_arr = cz_dist*np.sin(np.radians(dec))
coords = np.column_stack((x_arr,y_arr,z_arr))
return coords
In [4]:
def calc_dens(n_val,r_val):
"""
Returns densities of spheres with radius being the distance to the
nth nearest neighbor.
Parameters
----------
n_val = integer
The 'N' from Nth nearest neighbor
r_val = array-like
An array with the distances to the Nth nearest neighbor for
each galaxy
Returns
-------
dens: array-like
An array with the densities of the spheres created with radii
to the Nth nearest neighbor.
"""
dens = np.array([(3.*(n_val+1)/(4.*np.pi*r_val[hh]**3)) \
for hh in range(len(r_val))])
return dens
In [5]:
def plot_calcs(mass,mass_err=False,ratio_err=False):
"""
Returns values for plotting the stellar mass function and
mass ratios
Parameters
----------
mass: array-like
A 1D array with mass values
Returns
-------
bin_centers: array-like
An array with the medians of the mass bins
"""
bins = np.linspace(9.2,11.8,14)
dlogM = 0.2
mass_counts, edges = np.histogram(mass,bins)
bin_centers = 0.5*(edges[:-1]+edges[1:])
mass_freq = mass_counts/float(len(mass))/dlogM
ratio_dict = {}
frac_val = [2,4,10]
for ii in frac_val:
ratio_dict[ii] = {}
# Calculations for the lower density cut
frac_data = int(len(mass)/ii)
frac_mass = mass[0:frac_data]
counts, edges = np.histogram(frac_mass,bins)
# Calculations for the higher density cut
frac_mass_2 = mass[-frac_data:]
counts_2, edges_2 = np.histogram(frac_mass_2,bins)
if ratio_err == True:
yerr = []
yerr.append(((counts_2*1.)/(counts*1.)*np.sqrt(1./counts + 1./counts_2)))
# Ratio determination
ratio_counts = (1.*counts_2)/(1.*counts)
ratio_dict[ii] = ratio_counts
if ratio_err == True:
ratio_dict_list = [[] for xx in xrange(2)]
ratio_dict_list[0] = ratio_dict
ratio_dict_list[1] = yerr
ratio_dict = ratio_dict_list
if mass_err == True:
mass_freq_list = [[] for xx in xrange(2)]
mass_freq_list[0] = mass_freq
mass_freq_list[1] = np.sqrt(mass_counts)/float(len(mass))/dlogM
mass_freq = np.array(mass_freq_list)
return bin_centers, mass_freq, ratio_dict
In [6]:
def plot_ratios(data,nn_val,ax,ii,zz):
if zz==16 :
ax.set_xlabel('$\log\ M_{*}$',fontsize=18)
if zz%3 != 0:
# ax.get_yaxis().set_visible(False)
ax.yaxis.set_ticklabels([])
# ax.get_yaxis().set_ticks([])
bins = np.linspace(9.2,11.8,14)
if ii==0:
title_label = 'Mass Ratio 50/50, {0} NN'.format(nn_val)
frac_val = 2
ax.text(0.05, 0.95, title_label,horizontalalignment='left',verticalalignment='top',transform=ax.transAxes,fontsize=12)
elif ii==1:
title_label = 'Mass Ratio 25/75, {0} NN'.format(nn_val)
frac_val = 4
ax.text(0.05, 0.95, title_label,horizontalalignment='left',verticalalignment='top',transform=ax.transAxes,fontsize=12)
elif ii==2:
title_label = 'Mass Ratio 10/90, {0} NN'.format(nn_val)
frac_val = 10
ax.text(0.05, 0.95, title_label,horizontalalignment='left',verticalalignment='top',transform=ax.transAxes,fontsize=12)
# Calculations
frac_data = (len(data)/frac_val)
frac_mass = data[0:frac_data]
counts, edges = np.histogram(frac_mass,bins)
bins_cens = 0.5*(edges[:-1]+edges[1:])
#
frac_mass_2 = data[-frac_data:]
counts_2, edges_2 = np.histogram(frac_mass_2,bins)
ratio_counts = (1.*counts_2)/(1.*counts)
x_plots = bins_cens
yerr = ((counts_2*1.)/(counts*1.)*np.sqrt(1./counts + 1./counts_2))
# return yerr
ax.axhline(y=1,c="blue",linewidth=0.5,zorder=0)
ax.set_ylim([0,5])
ax.set_xlim([9.2,11.8])
ax.tick_params(axis='both', labelsize=12)
ax.set_xticks(np.arange(9., 12., 0.5))
ax.set_yticks([1.,3.])
ax.plot(x_plots,ratio_counts)
ax.errorbar(x_plots,ratio_counts,yerr = yerr)
In [8]:
neigh_vals = [1,2,3,5,10,20]
eco_path = r"C:\Users\Hannah\Desktop\Vanderbilt_REU\Stellar_mass_env_density\Catalogs\ECO_true"
ECO_true = (Index(eco_path,'.txt'))
names = ['ra','dec','cz','logMstar']
PD_eco = pd.read_csv(ECO_true[0],sep="\s+", usecols=(0,1,2,4),header=None,\
skiprows=1,names=names)
eco_comp = PD_eco[PD_eco.logMstar >= 9.3]
ra_eco = (np.array(eco_comp)).T[0]
dec_eco = (np.array(eco_comp)).T[1]
cz_eco = (np.array(eco_comp)).T[2]
mass_eco = (np.array(eco_comp)).T[3]
coords_eco = sph_to_cart(ra_eco,dec_eco,cz_eco)
eco_neighbor_tree = spatial.cKDTree(coords_eco)
eco_tree_dist = np.array(eco_neighbor_tree.query(coords_eco,21)[0])
eco_mass_dist = np.column_stack((mass_eco,eco_tree_dist.T[neigh_vals].T))
eco_dens = ([calc_dens(neigh_vals[jj],\
(eco_mass_dist.T[range(1,7)[jj]])) for jj in range(len(neigh_vals))])
eco_mass_dens = [(np.column_stack((mass_eco,eco_dens[ii]))) for ii in range(len(neigh_vals))]
eco_idx = [(eco_mass_dens[jj].T[1].argsort()) for jj in range(len(neigh_vals))]
eco_mass_dat = [(eco_mass_dens[jj][eco_idx[jj]].T[0]) for jj in range(len(neigh_vals))]
# eco_ratio_info = [[] for xx in xrange(len(eco_mass_dat))]
# for ii in range(len(eco_mass_dat)):
# bin_centers, eco_freq, eco_ratio_info[ii] = plot_calcs(eco_mass_dat[ii],mass_err=True,ratio_err=True)
# for jj in (range(len(eco_mass_dat))):
# eco_medians = np.array(bin_func(eco_mass_dist,(jj+1),bootstrap=True))
In [9]:
np.seterr(divide='ignore',invalid='ignore')
nn_val_arr = [1,2,3,5,10,20]
nrow_num = int(6)
ncol_num = int(3)
%matplotlib qt
fig, axes = plt.subplots(nrows=nrow_num, ncols=ncol_num, figsize=(100,200), sharex= True)
plt.subplots_adjust(left=0.02, bottom=0.09, right=1.00, top=1.00,hspace=0,wspace=0)
# figtitle = fig.suptitle("Mass Distributions by NN Density", fontsize=18)
# figtitle.set_y(1.04)
axes_flat = axes.flatten()
zz = int(0)
yerr = [[[] for ii in range(0,3)] for yy in range(len(eco_mass_dat))]
for kk in range(len(eco_mass_dat)):
for ii in range(0,3):
yerr[kk][ii] = plot_ratios(eco_mass_dat[kk], nn_val_arr[kk], axes_flat[zz],ii,zz)
zz += 1
In [10]:
# print (yerr)
plt.show()