In [284]:
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
from importlib import reload
import mathutils
import fitsio
import setcover
import archespec
import scipy.optimize as op
from scipy.ndimage import gaussian_filter1d
import cookb_signalsmooth
from matplotlib.colors import LogNorm
In [322]:
plt.rcParams["figure.figsize"]=(10,8)
plt.rcParams['axes.linewidth'] = 2
plt.rcParams['xtick.major.size'] = 15
plt.rcParams['xtick.major.width'] = 2
plt.rcParams['xtick.minor.size'] = 10
plt.rcParams['xtick.minor.width'] = 2
plt.rcParams['xtick.labelsize'] = 25
plt.rcParams['ytick.major.size'] = 15
plt.rcParams['ytick.major.width'] = 2
plt.rcParams['ytick.minor.size'] = 10
plt.rcParams['ytick.minor.width'] = 2
plt.rcParams['ytick.labelsize'] = 25
In [330]:
masterwave, allflux, allivar = archespec.rest_allspec_readin()
objs_ori = archespec.arche_readin()
nobj = objs_ori.size
z = objs_ori['Z']
mass = objs_ori['MASS']
sfr = objs_ori['SFR']
In [331]:
index_wave_all = np.searchsorted(masterwave, [3605., 6997.])
tmpflux = allflux[index_wave_all[0]:index_wave_all[1],:]
tmpivar = allivar[index_wave_all[0]:index_wave_all[1],:]
tmpwave = masterwave[index_wave_all[0]:index_wave_all[1]]
tmploglam = np.log10(tmpwave)
print(tmploglam.size, tmploglam.size//15*15)
In [332]:
itmp = (np.where(np.logical_and(np.logical_and(\
np.logical_and(np.logical_and(mass>8, mass<12), sfr>1e-4),\
z<0.3), z>0.05)))[0]
iuse = np.random.choice(itmp, 1000)
print(iuse.shape)
chi2 = np.zeros((iuse.size, iuse.size))
A = np.zeros((iuse.size, iuse.size))
In [333]:
newwave = np.median(tmpwave.reshape(tmploglam.size//15, 15), axis=1)
newflux = np.sum(tmpflux.reshape(tmploglam.size//15, 15, tmpflux.shape[1]), axis=1)
newivar = 1./np.sum(1./tmpivar.reshape(tmploglam.size//15, 15, tmpflux.shape[1]), axis=1)
In [334]:
tmp_yerr = 1./np.sqrt(newivar[:, iuse].T.reshape(iuse.size, newwave.size))
tmp_y = newflux[:,iuse].T
for i in np.arange(iuse.size):
tmp_x = newflux[:, iuse[i]].T.reshape(1,newwave.size)
tmp_xerr = 1./np.sqrt(newivar[:, iuse[i]].T.reshape(1,newwave.size))
A_tmp, chi2_tmp = mathutils.quick_amplitude(tmp_x, tmp_y, tmp_xerr, tmp_yerr)
A[i,:] = A_tmp
chi2[i,:] = chi2_tmp
In [399]:
a_matrix = chi2/newwave.size<3.5
cost = np.ones(iuse.size)
In [374]:
plt.hist(np.ravel(chi2/newwave.size), 40, range=[0,100])
Out[374]:
In [400]:
g = setcover.SetCover(a_matrix, cost)
In [401]:
g.CFT()
In [402]:
iarchetype = np.nonzero(g.s)[0]
print(iarchetype.shape)
tmpmedian = tmpflux[:,iuse[iarchetype]]
In [403]:
# Calculate [O II] EW
wave_oii = 3728.
index_oii = np.searchsorted(tmpwave, wave_oii)
flux_oii = np.sum(tmpmedian[index_oii-5:index_oii+5,:], axis=0)
flux_left = np.sum(tmpmedian[index_oii-15:index_oii-5,:], axis=0)
flux_right = np.sum(tmpmedian[index_oii+5:index_oii+15,:], axis=0)
ew_oii = flux_oii/(flux_left+flux_right)*2.
isort_oii = np.argsort(ew_oii)[::-1]
wave_ha = 6563.
index_ha = np.searchsorted(tmpwave, wave_ha)
flux_ha = np.sum(tmpmedian[index_ha-5:index_ha+5,:], axis=0)
flux_haleft = np.sum(tmpmedian[index_ha-115:index_ha-105,:], axis=0)
flux_haright = np.sum(tmpmedian[index_ha+105:index_ha+115,:], axis=0)
ew_ha = flux_ha/(flux_haleft+flux_haright)*2.
isort_ha = np.argsort(ew_ha)[::-1]
wave_norm = 4600.
index_norm = np.searchsorted(tmpwave, wave_norm)
flux_norm = np.median(tmpmedian[index_norm-5:index_norm+5,:], axis=0)
wave_left = 4200.
index_left = np.searchsorted(tmpwave, wave_left)
flux_left = np.sum(tmpmedian[index_left-15:index_left+15,:], axis=0)
wave_right = 5500.
index_right = np.searchsorted(tmpwave, wave_right)
flux_right = np.sum(tmpmedian[index_right-15:index_right+15,:], axis=0)
color_cont = flux_right/flux_left
isort_cont = np.argsort(color_cont)
In [404]:
#i = -1
i += 1
fig1 = plt.figure(figsize=(20,10))
ax1 = fig1.add_subplot(111)
fig2 = plt.figure(figsize=(20,10))
ax2 = fig2.add_subplot(111)
#x = gaussian_filter1d(tmpwave, 1.5)
x = tmpwave
for i in np.arange(iarchetype.size):
#y = gaussian_filter1d(tmpmedian[:,isort_oii[i]]/flux_norm[isort_oii[i]], 1.5)
#y = cookb_signalsmooth.smooth(tmpmedian[:,isort_oii[i]]/flux_norm[isort_oii[i]], 20)
#iwave_use = np.where(tmpmedian[:,isort_cont[i]]>1.E-5)[0]
iwave_use = np.where(tmpmedian[:,isort_ha[i]]>1.E-5)[0]
y = cookb_signalsmooth.smooth(tmpmedian[iwave_use,isort_ha[i]]/flux_norm[isort_ha[i]], 10)
#ax.plot(x, y)
ax1.plot(x[iwave_use], y)
ax2.loglog(x[iwave_use], y)
ax1.set_xlim(3600, 7200)
ax2.set_xlim(3600, 7200)
ax1.set_xticks((np.arange(4)+1.)*1000.+3000.)
ax2.set_xticks((np.arange(4)+1.)*1000.+3000.)
ax1.set_ylim(1E-1, 5E0)
ax2.set_ylim(1E-1, 2E1)
filefile1 = "/Users/Benjamin/Dropbox/Zhu_Projects/Archetype/Movie/Full50_"+'%03d'%i+".jpg"
filefile2 = "/Users/Benjamin/Dropbox/Zhu_Projects/Archetype/Movie/Loglog_Full50_"+'%03d'%i+".jpg"
fig1.savefig(filefile1)
fig2.savefig(filefile2)
ax1.clear()
ax2.clear()
#print(objs_ori['RA'][iuse[iarchetype[isort_oii[i]]]], objs_ori['DEC'][iuse[iarchetype[isort_oii[i]]]])
#print(objs_ori['RA'][iuse[iarchetype]])
for i in (np.arange(iarchetype.size))[::-1]:
print(objs_ori['RA'][iuse[iarchetype[isort_ha[i]]]], objs_ori['DEC'][iuse[iarchetype[isort_ha[i]]]])
In [395]:
for i in (np.arange(iarchetype.size))[::-1]:
print(objs_ori['RA'][iuse[iarchetype[isort_oii[i]]]], objs_ori['DEC'][iuse[iarchetype[isort_oii[i]]]])
In [106]:
150./360.*24
Out[106]:
In [107]:
print(isort_oii)
In [221]:
a_school = np.zeros((11,11), dtype=bool)
cost_school = np.ones(11)
In [222]:
# 1
a_school[0:4, 0] = True
# 2
a_school[0:3, 1] = True
a_school[4, 1] = True
# 3
a_school[0:6, 2] = True
# 4
a_school[0, 3] = True
a_school[2:4, 3] = True
a_school[5:7, 3] = True
# 5
a_school[1:3, 4] = True
a_school[4:6, 4] = True
a_school[7:9, 4] = True
# 6
a_school[2:8, 5] = True
# 7
a_school[3, 6] = True
a_school[5:8, 6] = True
# 8
a_school[4:10, 7] = True
# 9
a_school[4, 8] = True
a_school[7:10, 8] = True
# 10
a_school[7:11, 9] = True
# 11
a_school[9:11, 10] = True
In [223]:
print(a_school*0.1)
In [224]:
fig = plt.figure(figsize=(11,11))
ax = fig.add_subplot(111)
mat = ax.imshow(a_school, cmap='hot_r', interpolation='nearest', vmax=3., vmin=0)
#plt.yticks(np.arange(a_school.shape[0]-1)+0.5, np.arange(a_school.shape[0]-1)+1)#+1)
#plt.xticks(np.arange(a_school.shape[1]-1)+0.5, np.arange(a_school.shape[0]-1)+1)#+1)
plt.yticks(np.arange(a_school.shape[0]-1)+0.5, [])#+1)
plt.xticks(np.arange(a_school.shape[1]-1)+0.5, [])#+1)
#plt.xticks(rotation=30)
#plt.xlabel('Concert Dates')
# this places 0 or 1 centered in the individual squares
for x in np.arange(a_school.shape[0]):
for y in np.arange(a_school.shape[1]):
ax.annotate(str(a_school[x, y].astype(int))[0], xy=(y, x),
horizontalalignment='center', verticalalignment='center', fontsize=20)
ax.grid(True)
#plt.show()
In [256]:
fig = plt.figure(figsize=(11,11))
ax = fig.add_subplot(111)
#ax.hist2d(dr7_mass['MEDIAN'][iarche], np.log10(dr7_sfr['MEDIAN'][iarche]), bins=160)
#ax.plot(dr7_mass['MEDIAN'][iarche], np.log10(dr7_sfr['MEDIAN'][iarche]), '.')
ax.plot(dr7_mass['MEDIAN'][iarche], dr7_sfr['MEDIAN'][iarche], '.')
ax.set_ylim(2, -2)
Out[256]:
In [278]:
dr7_info = fitsio.read('/Users/Benjamin/AstroData/Garching/gal_info_dr7_v5_2.fit.gz', ext=1)
dr7_ssfr = fitsio.read('/Users/Benjamin/AstroData/Garching/gal_totspecsfr_dr7_v5_2.fits.gz', ext=1)
dr7_mass = fitsio.read('/Users/Benjamin/AstroData/Garching/totlgm_dr7_v5_2.fit.gz', ext=1)
In [369]:
idr7 = np.nonzero(np.logical_and(np.logical_and(np.logical_and(np.logical_and(np.logical_and(np.logical_and(
dr7_info['Z']>0.02, dr7_info['Z']<0.30),
dr7_mass['MEDIAN']>7.), dr7_mass['MEDIAN']<13.),
dr7_ssfr['MEDIAN']<-1), dr7_ssfr['MEDIAN']>-20),
dr7_info['SN_MEDIAN']>3.))[0]
iarche_dr7 = np.random.choice(idr7, 10000)
print(iarche_dr7.shape)
In [296]:
print(dr7_info.size, dr7_mass.size, dr7_sfr.size)
In [248]:
dr7_sfr.dtype
Out[248]:
In [405]:
fig = plt.figure(figsize=(11,11))
ax = fig.add_subplot(111)
ax.hist2d(dr7_mass['MEDIAN'][idr7], dr7_ssfr['MEDIAN'][idr7], bins=80, cmin=3, cmap='Greys', norm=LogNorm())
#ax.plot(dr7_mass['MEDIAN'][iarche], np.log10(dr7_sfr['MEDIAN'][iarche]), '.')
#ax.plot(dr7_mass[1000:11000]['MEDIAN'], dr7_sfr[1000:11000]['MEDIAN'], '.')
ax.plot(mass[iarchetype], sfr[iarchetype]-mass[iarchetype], '*', ms=15)
ax.set_ylim(-7, -14)
ax.set_xlim(7, 13.)
ax.set_xlabel(r'$\log_{10}$ $M_*$ [M$_\odot$]', fontsize=22)
ax.set_ylabel(r'sSFR [yr$^{-1}$]', fontsize=22)
print(iarchetype.size)
In [410]:
fig = plt.figure(figsize=(11,11))
ax = fig.add_subplot(111)
mat = ax.imshow(a_matrix[0:100, 0:100], cmap='hot_r', interpolation='nearest', vmax=3., vmin=0)
#plt.yticks(np.arange(a_school.shape[0]-1)+0.5, np.arange(a_school.shape[0]-1)+1)#+1)
#plt.xticks(np.arange(a_school.shape[1]-1)+0.5, np.arange(a_school.shape[0]-1)+1)#+1)
#plt.yticks(np.arange(a_matrix.shape[0]-1)+0.5, [])#+1)
#plt.xticks(np.arange(a_matrix.shape[1]-1)+0.5, [])#+1)
ax.set_xticklabels([])
ax.set_yticklabels([])
ax.set_xticks([])
ax.set_yticks([])
#plt.xticks(rotation=30)
#plt.xlabel('Concert Dates')
# this places 0 or 1 centered in the individual squares
#for x in np.arange(a_school.shape[0]):
# for y in np.arange(a_school.shape[1]):
# ax.annotate(str(a_school[x, y].astype(int))[0], xy=(y, x),
# horizontalalignment='center', verticalalignment='center', fontsize=20)
#ax.grid(True)
#plt.show()
Out[410]:
In [445]:
fig = plt.figure(figsize=(11,11))
ax = fig.add_subplot(111)
mat = ax.imshow(a_matrix, cmap='hot_r', interpolation='nearest', vmax=2.6, vmin=0)
#plt.yticks(np.arange(a_school.shape[0]-1)+0.5, np.arange(a_school.shape[0]-1)+1)#+1)
#plt.xticks(np.arange(a_school.shape[1]-1)+0.5, np.arange(a_school.shape[0]-1)+1)#+1)
#plt.yticks(np.arange(a_matrix.shape[0]-1)+0.5, [])#+1)
#plt.xticks(np.arange(a_matrix.shape[1]-1)+0.5, [])#+1)
ax.set_xticklabels([])
ax.set_yticklabels([])
ax.set_xticks([])
ax.set_yticks([])
#plt.xticks(rotation=30)
#plt.xlabel('Concert Dates')
# this places 0 or 1 centered in the individual squares
#for x in np.arange(a_school.shape[0]):
# for y in np.arange(a_school.shape[1]):
# ax.annotate(str(a_school[x, y].astype(int))[0], xy=(y, x),
# horizontalalignment='center', verticalalignment='center', fontsize=20)
#ax.grid(True)
#plt.show()
Out[445]:
In [461]:
fig = plt.figure(figsize=(11,11))
ax = fig.add_subplot(111)
mat = ax.imshow(chi2[0:11, 0:11]/newwave.size, cmap='hot', interpolation='nearest', vmax=8., vmin=-2)
#plt.yticks(np.arange(a_school.shape[0]-1)+0.5, np.arange(a_school.shape[0]-1)+1)#+1)
#plt.xticks(np.arange(a_school.shape[1]-1)+0.5, np.arange(a_school.shape[0]-1)+1)#+1)
plt.yticks(np.arange(11-1)+0.5, [])#+1)
plt.xticks(np.arange(11-1)+0.5, [])#+1)
#ax.set_xticklabels([])
#ax.set_yticklabels([])
#ax.set_xticks([])
#ax.set_yticks([])
#plt.xticks(rotation=30)
#plt.xlabel('Concert Dates')
# this places 0 or 1 centered in the individual squares
#for x in np.arange(a_school.shape[0]):
# for y in np.arange(a_school.shape[1]):
# ax.annotate(str(a_school[x, y].astype(int))[0], xy=(y, x),
# horizontalalignment='center', verticalalignment='center', fontsize=20)
ax.grid(True)
#plt.show()
In [438]:
fig = plt.figure(figsize=(11,11))
ax = fig.add_subplot(111)
mat = ax.imshow(a_matrix[0:11, 0:11], cmap='hot_r', interpolation='nearest', vmax=3., vmin=0)
#plt.yticks(np.arange(a_school.shape[0]-1)+0.5, np.arange(a_school.shape[0]-1)+1)#+1)
#plt.xticks(np.arange(a_school.shape[1]-1)+0.5, np.arange(a_school.shape[0]-1)+1)#+1)
plt.yticks(np.arange(11-1)+0.5, [])#+1)
plt.xticks(np.arange(11-1)+0.5, [])#+1)
ax.set_xticklabels([])
ax.set_yticklabels([])
#ax.set_xticks([])
#ax.set_yticks([])
#plt.xticks(rotation=30)
#plt.xlabel('Concert Dates')
# this places 0 or 1 centered in the individual squares
#for x in np.arange(a_school.shape[0]):
# for y in np.arange(a_school.shape[1]):
# ax.annotate(str(a_school[x, y].astype(int))[0], xy=(y, x),
# horizontalalignment='center', verticalalignment='center', fontsize=20)
ax.grid(True)
#plt.show()
In [458]:
fig = plt.figure(figsize=(12,6))
ax = fig.add_subplot(111)
x = tmpwave
iwave_use = np.where(tmpmedian[:,4]>1.E-5)[0]
y = cookb_signalsmooth.smooth(tmpmedian[iwave_use,4], 10)
ax.plot(x[iwave_use], y)
ax.set_xticks([])
ax.set_yticks([])
Out[458]:
In [455]:
fig = plt.figure(figsize=(12,6))
ax = fig.add_subplot(111)
x = tmpwave
iwave_use = np.where(tmpmedian[:,8]>1.E-5)[0]
y = cookb_signalsmooth.smooth(tmpmedian[iwave_use,8], 10)
ax.plot(x[iwave_use], y)
ax.set_xticks([])
ax.set_yticks([])
Out[455]:
In [462]:
g_school = setcover.SetCover(a_school, cost_school)
In [467]:
g_school.CFT()
In [468]:
np.nonzero(g_school.s)
Out[468]:
In [470]:
print(a_school)
In [ ]: