In [1]:
%matplotlib inline
from galaxy_analysis.plot.plot_styles import *
import matplotlib.pyplot as plt
import yt
import numpy as np
import deepdish as dd
from scipy.optimize import brentq
from galaxy_analysis.utilities import utilities
from galaxy_analysis.analysis import Galaxy
from galaxy_analysis.utilities import functions
from galaxy_analysis.particle_analysis import particle_types as ptypes
import time
WDIR = '/home/aemerick/work/enzo_runs/sn_H2atten_H2sh/'
In [2]:
reload(functions)
Out[2]:
In [2]:
gal = Galaxy('DD0400', wdir = WDIR)
In [3]:
# stars_lifetime = gal.df[('io','particle_model_lifetime')]
In [4]:
data = dd.io.load(WDIR + 'gas_abundances.h5')
In [30]:
nrow = 2
ncol = 3
fig, all_ax = plt.subplots(nrow,ncol)
fig.set_size_inches(8*ncol,8*nrow)
fbins = data['DD0200']['Disk']['mass_fraction']['bins']
binsize = (fbins[1:]) - (fbins[:-1])
centers = 0.5 * (fbins[1:] + fbins[:-1])
phase = 'Molecular'
field = 'O_Fraction'
ylim = [0.0,6]
ds_list = np.arange(120.0, 401, 10)
min_median = 1.0E80
max_median = -1.0E80
axi = 0
axj = 0
phases = ['CNM','WNM','WIM','HIM','Disk']
phase_labels = ['Cold Neutral','Warm Neutral','Warm Ionized','Hot Ionized','All ISM']
mean_phase = {}
std_phase = {}
actual_mean_phase = {}
actual_median_phase = {}
actual_std_phase = {}
time = np.zeros(np.size(ds_list))
for k in phases:
mean_phase[k] = np.zeros(np.size(ds_list))
std_phase[k] = np.zeros(np.size(ds_list))
actual_mean_phase[k] = np.zeros(np.size(ds_list))
actual_median_phase[k] = np.zeros(np.size(ds_list))
actual_std_phase[k] = np.zeros(np.size(ds_list))
xmin = 10000
xmax = -1000000
for i,j in enumerate(ds_list):
dsname = 'DD%0004i'%(j)
axi = 0
axj = 0
for ip, phase in enumerate(phases):
# print axi, axj, phase
ax = all_ax[(axi, axj)]
y = data[dsname][phase]['mass_fraction'][field]['hist']
mean = data[dsname][phase]['mass_fraction'][field]['mean']
std = data[dsname][phase]['mass_fraction'][field]['std']
median = np.log10(np.interp(0.5, np.cumsum(y)/(1.0*np.sum(y)), centers))
label = data[dsname]['general']['Time'] - 46.0
min_median = np.min([median, min_median])
max_median = np.max([median, max_median])
frac = 0.8
alpha = ((i / (np.size(ds_list)*1.0))) * frac
halpha = ((i / (np.size(ds_list)*1.0))) * frac * 0.5
color = magma(alpha/frac) #'black' # viridis(alpha)
lw = 4.0
hlw = 4.0
if j == np.max(ds_list):
color = 'black'
alpha = 1.0
lw = 5.0
hlw = 5.0
norm_y = y / binsize
plot_bins = np.log10(fbins)
plot_histogram(ax, plot_bins, norm_y/np.max(norm_y), lw = hlw, color = color, alpha = halpha)
selection = (y > 0) * (centers > 10**(median - 4)) * (centers < 10**(median+4))
fit_x = centers[selection]
y_to_fit = norm_y[selection]
# take initial guess parameters from distribution values - compute logged dist values
try:
u_guess = np.log( mean / (np.sqrt(1.0 + std*std/(mean*mean))))
std_guess = np.sqrt(np.log(1.0 + std*std/(mean*mean)))
xplot = np.logspace(np.log10(np.min(fit_x)), np.log10(np.max(fit_x)),4000)
lognorm = functions.lognormal()
popt, pcov = lognorm.fit_function(fit_x, y_to_fit, p0 = [u_guess, std_guess])
xplot = np.logspace(np.log10(np.min(fit_x)), np.log10(np.max(fit_x)),4000)
yplot = lognorm._f(xplot, *popt)
ax.plot(np.log10(xplot), yplot/np.max(yplot), lw = lw, color = color , ls = '-', alpha = alpha) #, label = 'Fit')
except:
print phase, dsname, 'failure'
popt = [None, None]
mean_phase[phase][i] = popt[0]
actual_mean_phase[phase][i] = mean
actual_median_phase[phase][i] = median
std_phase[phase][i] = popt[1]
actual_std_phase[phase][i] = std
time[i] = label
axj = axj + 1
if axj >= ncol:
axi = axi + 1
axj = 0
xmin = np.min([np.log10(xplot[0]),xmin])
xmax = np.max([np.log10(xplot[-1]),xmax])
for i in np.arange(2):
for j in np.arange(3):
all_ax[(i,j)].set_xlim(xmax - 5, xmax)
all_ax[(i,j)].semilogy()
all_ax[(i,j)].set_ylim(1.0E-2,1.0)
ba = 8
ba = 0
all_ax[(i,j)].set_xlim(-7 - ba, -2.5 - ba)
xy = (np.min(all_ax[(i,j)].get_xlim())+ 0.2,np.max(all_ax[(i,j)].get_ylim())*0.95)
all_ax[(i,j)].annotate(phase_labels[(3)*(i) + j], xy = xy, xytext=xy)
all_ax[(i,j)].set_xlabel(r'log(Oxygen Mass Fraction)')
all_ax[(i,j)].set_ylabel(r'Peak Normalized PDF')
#ax.set_xlim( np.log10(min_median) - 2, np.log10(max_median) + 2)
#fig.savefig(field + '_phase_evolution.png')
fig.savefig('log_panel_phase_evolution.png')
In [ ]:
In [7]:
print data['DD0500']['CNM']['mass_fraction']['O_Fraction'].keys()
print np.sort(data.keys())
In [54]:
data[dsname][phase]['mass_fraction']['O_Fraction'].keys()
Out[54]:
In [5]:
from scipy import integrate
from scipy.stats import distributions
from scipy import stats
def random_sample_from_cdf(npoints, xval, CDF):
u = np.random.rand(npoints)
rsample = np.interp(u, CDF, xval)
return rsample
def _load_data(ldata, dsname, phase, field, centers = None, IQR = False, return_dict = False):
y = ldata[dsname][phase]['mass_fraction'][field]['hist']
mean = ldata[dsname][phase]['mass_fraction'][field]['mean']
std = ldata[dsname][phase]['mass_fraction'][field]['std']
q1 = ldata[dsname][phase]['mass_fraction'][field]['Q1']
q3 = ldata[dsname][phase]['mass_fraction'][field]['Q3']
q10 = ldata[dsname][phase]['mass_fraction'][field]['decile_1']
q90 = ldata[dsname][phase]['mass_fraction'][field]['decile_9']
if q3 is None or q1 is None:
iqr = None
q90_q10_range = None
else:
iqr = np.log10(q3) - np.log10(q1) # ldata[dsname][phase]['mass_fraction'][field]['inner_quartile_range']
q90_q10_range = np.log10(q90) - np.log10(q10)
label = ldata[dsname]['general']['Time'] - 46.0
#if (not ('median' in ldata[dsname][phase]['mass_fraction'][field].keys())):
median = np.log10(np.interp(0.5, np.cumsum(y)/(1.0*np.sum(y)), centers))
#else:
# print 'loading median from file'
# median = ldata[dsname][phase]['mass_fraction'][field]['median']
if return_dict:
rdict = {'y' : y, 'mean' : mean, 'median' : median, 'std' : std, 'label' : label, 'q1' : q1, 'q3' : q3,
'iqr' : iqr, 'q10' : q10, 'q90' : q90, 'q90_q10_range' : q90_q10_range}
return rdict
else:
if IQR:
return y, mean, median, std, label, iqr
else:
return y, mean, median, std, label
def _poisson_error(y):
N = np.sum(y) * 1.0
error = np.sqrt(y*y*(N+y)/(N**3))
return error
def _compute_fit_and_plot(ax, bins, y, mean=None, median=None, std=None, show_hist = True, plot_PDF = True,
plot_CDF = False, alpha = 0.5, halpha = 0.25, color = 'black', ls = '-', Mtot = 1.0E4,
fit_only = False, method = 'curve_fit', x_shift = 0.0, peak_normalize = True):
D_KS = None
D_KS_alpha = None
CDF_diff = None
KS_alpha = None
KS_result = None
KS_pval = None
if plot_CDF:
plot_PDF = False
binsize = (bins[1:]) - (bins[:-1])
centers = 0.5 * (bins[1:] + bins[:-1])
frac = 0.8
lw = 4.0
hlw = 4.0
#if j == np.max(ds_list):
# color = 'black'
# alpha = 1.0
# lw = 5.0
# hlw = 5.0
norm_y = y / binsize
plot_bins = np.log10(bins)
if show_hist and (not fit_only):
if plot_PDF:
if peak_normalize:
_y = np.max(norm_y)
else:
_y = 1.0
plot_histogram(ax, plot_bins + x_shift, norm_y / _y, lw = hlw, color = color, alpha = halpha, ls = ls)
elif plot_CDF:
plot_histogram(ax, plot_bins + x_shift, np.cumsum(y), lw = hlw, color = color, alpha = halpha, ls = ls)
else:
plot_histogram(ax, plot_bins + x_shift, y, lw = hlw, color = color, alpha = halpha, ls = ls)
if ((mean is not None) and (median is not None) and (std is not None)):
selection = (y > 0) * (centers > 10**(median - 4)) * (centers < 10**(median+4))
selection = (y > np.max(y)*1.0E-3) # * (centers > 10**(median - 4)) * (centers < 10**(median + 4))
if np.size(centers[selection]) < 3:
selection = (y > 0) * (centers > 10**(median-4)) * (centers < 10**(median+4))
if np.size(centers[selection]) < 3:
print y
fit_x = centers[selection]
y_to_fit = norm_y[selection]
# print fit_x
# take initial guess parameters from distribution values - compute logged dist values
if True:
#error = _poisson_error(y[selection]*Mtot) / Mtot
#norm_error = error / binsize[selection]
norm_error = np.sqrt(y[selection]*Mtot) / binsize[selection] / Mtot
#norm_error = 0.01 * norm_y[selection]
#if np.size(norm_error[ norm_error > y_to_fit]) > 0:
# norm_error[norm_error > y_to_fit] = 0.95*y_to_fit[norm_error>y_to_fit]
if np.size(norm_error[norm_error == 0.0]) > 0:
print 'getting zero error'
# norm_error[norm_error == 0.0] = np.median(norm_error[norm_error>0])
u_guess = np.log( mean / (np.sqrt(1.0 + std*std/(mean*mean))))
std_guess = np.sqrt(np.log(1.0 + std*std/(mean*mean)))
xplot = np.logspace(np.log10(np.min(fit_x)), np.log10(np.max(fit_x)),4000)
lognorm = functions.lognormal()
# print mean, std, u_guess, std_guess
#popt, pcov = lognorm.fit_function(fit_x, y_to_fit, p0 = [u_guess, std_guess],
# # sigma = norm_error, absolute_sigma = False,
# bounds = ([-60,0],[0,10]))
popt, pcov = lognorm.fit_function(fit_x, y_to_fit, p0 = [u_guess, std_guess,],
# sigma = norm_error, absolute_sigma = False,
bounds = ([u_guess - 10, u_guess + 10], [-np.inf, np.inf]),
method = method, data_cdf = np.cumsum(y[selection]))
xplot = np.logspace(np.log10(np.min(fit_x))-3, np.log10(np.max(fit_x))+3,4000)
yplot = lognorm._f(xplot, *popt)
#
# Testing soem things here for goodness of fit estimation
#
# 1) chi squared - errors are hard
chisqr = np.sum( (lognorm._f(fit_x, *popt) - y_to_fit)**2 / (norm_error)**2)
reduced_chisqr = chisqr / (np.size(y_to_fit) - 2.0)
# 2) Cumulative distribution and computing D between the two distributions
mf = np.zeros(np.size(bins)-1)
func = lambda x : lognorm._f(x, *popt)
for i in np.arange(0, np.size(bins)-1):
mf[i] = integrate.quad(func, bins[i], bins[i+1])[0]
mf = np.cumsum(mf)
CDF_diff = np.cumsum(y) - mf
D_KS = np.max(np.abs(CDF_diff))
D_KS_alpha = np.sqrt(2.0/ (1.0*np.size( y[y>0]) )*(np.size(y[y>0])))
KS_alpha = 2.0 * np.log(- 2.0 *(D_KS / D_KS_alpha)**2)
# try this way
#fit_sample = random_sample_from_CDF(10000, x[y>0], lognorm._CDF(x[y>0], *popt))
#nsample = 20
#dat_sample = random_sample_from_cdf(nsample, centers[y>0], np.cumsum(y[y>0]))
#func = lambda x : lognorm._CDF(x, *popt)
#KS_result = stats.kstest(dat_sample, func)
# compute p-value using scipy
N = np.size(y[y>0])
pval_two = distributions.kstwobign.sf(D_KS*np.sqrt(N))
if N > 2666 or pval_two > 0.80 - N*0.3/1000.0:
KS_pval = distributions.kstwobign.sf(D_KS*np.sqrt(N))
else:
KS_pval = distributions.ksone.sf(D_KS,N)*2
if not fit_only:
if plot_PDF:
if peak_normalize:
_y = np.max(yplot)
else:
_y = 1.0
ax.plot(np.log10(xplot) + x_shift, yplot/_y, lw = lw, color = color , ls = ls, alpha = alpha) #, label = 'Fit')
elif plot_CDF:
cplot = (bins[1:] + bins[:-1])*0.5
ax.plot(np.log10(cplot) + x_shift, mf, lw = lw, color = color, ls = ls, alpha = alpha)
ax.plot([np.log10(cplot)[-1] + x_shift, 0 + x_shift], [mf[-1]]*2, lw = lw, color = color, ls = ls, alpha = alpha)
else:
mf = np.zeros(np.size(bins)-1)
func = lambda x : lognorm._f(x, *popt)
for i in np.arange(0, np.size(bins)-1):
mf[i] = integrate.quad(func, bins[i], bins[i+1])[0]
if not fit_only:
cplot = (bins[1:] + bins[:-1])*0.5
ax.plot(np.log10(cplot) + x_shift, mf, lw = lw, color = color, ls = ls, alpha = alpha)
else:
print phase, dsname, 'failure'
popt = [None, None]
# print np.trapz(yplot,xplot)
pdf_mean = np.trapz(centers * norm_y, x = centers)
pdf_mf = np.zeros(np.size(bins)-2)
func = lambda x : np.interp(x, centers, norm_y)
for i in np.arange(0, np.size(bins)-2):
pdf_mf[i] = integrate.quad(func, bins[i+1], bins[i+2])[0]
pdf_mf = np.cumsum(pdf_mf)
pdf_median = np.interp(0.5, pdf_mf, centers[1:])
pdf_mode = np.argmax( norm_y )
pdf_std = np.sqrt( np.trapz(centers*centers*norm_y, x = centers) - pdf_mean**2)
# return x values of fit, mean of fit, and std of fit
mean_fit = popt[0]
std_fit = popt[1]
output_dict = {'mean' : mean_fit, 'std' : std_fit, 'chisqr' : chisqr, 'norm_error' : norm_error,
'rchisqr' : reduced_chisqr, 'npoints' : np.size(y_to_fit), 'N' : np.size(y_to_fit),
'fit_x' : fit_x, 'rel_error' : np.abs(norm_error/norm_y[selection]),
'CDF_diff' : CDF_diff, 'D_KS' : D_KS, 'D_KS_alpha' : D_KS_alpha, 'KS_alpha' : KS_alpha,
'num_pdf_mean' : pdf_mean, 'num_pdf_median' : pdf_median, 'num_pdf_mode' : pdf_mode,
'num_pdf_std' : pdf_std, 'KS_pval' : KS_pval, 'KS_result' : KS_result, 'og_mean' : mean}
return xplot, mean_fit, std_fit, output_dict
else:
if ((mean is not None) or (median is not None) or (std is not None)):
print "Must provide mean, median, and std to enable a lognormal fit"
raise ValueError
# return nothing if just plotting the line and doing no fitting
return
# mean_phase[phase][i] = popt[0]
# actual_mean_phase[phase][i] = mean
# std_phase[phase][i] = popt[1]
# actual_std_phase[phase][i] = std
# time[i] = label
In [11]:
fbins = data['DD0200']['Disk']['mass_fraction']['bins']
binsize = (fbins[1:]) - (fbins[:-1])
centers = 0.5 * (fbins[1:] + fbins[:-1])
phases = ['CNM','WNM','WIM','HIM','Disk']
ds_list = np.arange(120, 401, 10) # [ 100, 150, 200, 250, 300, 350]
ele = ['O','Fe','N','Ba']
count = 0
output_stats = {}
for dsi in ds_list:
dsname = 'DD%0004i'%(dsi)
output_stats[dsi] = {}
for phase in phases:
output_stats[dsi][phase] = {}
for e in ele:
field = e + '_Fraction'
rdict = _load_data(data, dsname, phase, field, centers = centers, return_dict = True)
y = rdict['y']
mean = rdict['mean']
median = rdict['median']
std = rdict['std']
iqr = rdict['iqr']
Mtot = data[dsname][phase]['general']['total_mass']
if (mean is None) or (np.size(y[y > 0]) == 0): # phase not present
for k in output:
output[k] = None
output_stats[dsi][phase][e] = output
continue
xplot, mean_fit, std_fit, output = _compute_fit_and_plot(None, fbins, y, mean = mean, median = median, std = std,
Mtot=Mtot, method = 'KS',
fit_only = True)
output['actual_mean'] = mean
output['actual_std'] = std
output['actual_median'] = median
output['actual_IQR'] = iqr
output['actual_q90q10'] = rdict['q90_q10_range']
output_stats[dsi][phase][e] = output
count += 1
if count >= 9:
print dsi
count = 0
# convert to an array we can plot as a function of time
output_stats_data = {}
for phase in phases:
output_stats_data[phase] = {}
for e in ele:
output_stats_data[phase][e] = {}
for k in output_stats[ds_list[0]][phase][e].keys():
output_stats_data[phase][e][k] = np.array( [output_stats[x][phase][e][k] for x in ds_list])
In [10]:
fig, ax = plt.subplots(2,2)
fig.set_size_inches(16,16)
for i, phase in enumerate(['Molecular','CNM','WNM','WIM']): #,'HIM']):
if phase == 'Disk':
continue
mean = output_stats_data[phase]['O']['actual_median']
std = output_stats_data[phase]['O']['actual_IQR']
# mean = np.exp(M + 0.5*S**2)
# std = np.sqrt(np.exp(S**2 + 2.0 * M)*(np.exp(S*S) - 1.0))
#select =
ax[(0,0)].plot(ds_list, mean, lw = 3, ls = '-', label = phase)
ax[(0,1)].plot(ds_list, std , lw = 3, ls = '-', label = phase)
ax[(0,0)].set_ylabel(r'Oxygen - Median')
ax[(0,1)].set_ylabel(r'Oxygen - IQR')
mean = output_stats_data[phase]['N']['actual_median']
std = output_stats_data[phase]['N']['actual_IQR']
# mean = np.exp(M + 0.5*S**2)
# std = np.sqrt(np.exp(S**2 + 2.0 * M)*(np.exp(S*S) - 1.0))
#select =
ax[(1,0)].plot(ds_list, mean, lw = 3, ls = '-', label = phase)
ax[(1,1)].plot(ds_list, std , lw = 3, ls = '-', label = phase)
ax[(1,0)].set_ylabel(r'Nitrogen - Median')
ax[(1,1)].set_ylabel(r'Nitrogen - IQR')
#ax[(0,1)].set_ylim(0,2)
#ax[(1,1)].set_ylim(0,2)
for a in [(0,0),(0,1),(1,0),(1,1)]:
# ax[a].set_xlim(0.0, np.max(time))
ax[a].set_xlabel('Time (Myr)')
#ax[1].semilogy()
#ax[0].set_ylim(-8, -2)
#ax[1].set_ylim(1.0E-1, 100.0)
ax[(0,0)].set_ylim(-7, -2)
ax[(1,0)].set_ylim(-7, -2)
ax[(0,1)].set_ylim(0.0, 1.0)
ax[(1,1)].set_ylim(0.0, 1.0)
ax[(0,0)].legend(loc='best')
print output_stats_data['CNM']['O']['actual_IQR']
In [16]:
fig, ax = plt.subplots(2,2)
fig.set_size_inches(16,16)
for i, phase in enumerate(['Molecular','CNM','WNM','WIM','HIM']):
if phase == 'Disk':
continue
mean = output_stats_data[phase]['O']['mean']
std = output_stats_data[phase]['O']['std']
# mean = np.exp(M + 0.5*S**2)
# std = np.sqrt(np.exp(S**2 + 2.0 * M)*(np.exp(S*S) - 1.0))
#select =
ax[(0,0)].plot(ds_list, mean, lw = 5, ls = '-', label = phase)
ax[(0,1)].plot(ds_list, std , lw = 5, ls = '-', label = phase)
ax[(0,0)].set_ylabel(r'O $\mu$')
ax[(0,1)].set_ylabel(r'O $\sigma$')
mean = output_stats_data[phase]['N']['mean']
std = output_stats_data[phase]['N']['std']
# mean = np.exp(M + 0.5*S**2)
# std = np.sqrt(np.exp(S**2 + 2.0 * M)*(np.exp(S*S) - 1.0))
#select =
ax[(1,0)].plot(ds_list, mean, lw = 5, ls = '-', label = phase)
ax[(1,1)].plot(ds_list, std , lw = 5, ls = '-', label = phase)
ax[(1,0)].set_ylabel(r'N $\mu$')
ax[(1,1)].set_ylabel(r'N $\sigma$')
ax[(0,1)].set_ylim(0,2)
ax[(1,1)].set_ylim(0,2)
for a in [(0,0),(0,1),(1,0),(1,1)]:
# ax[a].set_xlim(0.0, np.max(time))
ax[a].set_xlabel('Time (Myr)')
#ax[1].semilogy()
#ax[0].set_ylim(-8, -2)
#ax[1].set_ylim(1.0E-1, 100.0)
ax[(0,0)].set_ylim(-12,-5)
ax[(1,0)].set_ylim(-17,-10)
ax[(0,0)].legend(loc='best')
Out[16]:
In [18]:
fig, ax = plt.subplots(2,2)
fig.set_size_inches(16,16)
for i, phase in enumerate(['Molecular','CNM','WNM','WIM','HIM']):
if phase == 'Disk':
continue
mean = output_stats_data[phase]['O']['num_pdf_mean']
std = output_stats_data[phase]['O']['num_pdf_std']
# mean = np.exp(M + 0.5*S**2)
# std = np.sqrt(np.exp(S**2 + 2.0 * M)*(np.exp(S*S) - 1.0))
#select =
ln_mean = np.log( mean / np.sqrt(1.0 + std**2/mean**2))
ln_std = np.sqrt( np.log(1.0 + std**2/mean**2))
ax[(0,0)].plot(ds_list, ln_mean, lw = 5, ls = '-', label = phase)
ax[(0,1)].plot(ds_list, ln_std , lw = 5, ls = '-', label = phase)
ax[(0,0)].set_ylabel(r'O $\mu$')
ax[(0,1)].set_ylabel(r'O $\sigma$')
mean = output_stats_data[phase]['N']['num_pdf_mean']
std = output_stats_data[phase]['N']['num_pdf_std']
# mean = np.exp(M + 0.5*S**2)
# std = np.sqrt(np.exp(S**2 + 2.0 * M)*(np.exp(S*S) - 1.0))
#select =
ln_mean = np.log( mean / np.sqrt(1.0 + std**2/mean**2))
ln_std = np.sqrt( np.log(1.0 + std**2/mean**2))
ax[(1,0)].plot(ds_list, ln_mean, lw = 5, ls = '-', label = phase)
ax[(1,1)].plot(ds_list, ln_std , lw = 5, ls = '-', label = phase)
ax[(1,0)].set_ylabel(r'N $\mu$')
ax[(1,1)].set_ylabel(r'N $\sigma$')
ax[(0,1)].set_ylim(0,2)
ax[(1,1)].set_ylim(0,2)
for a in [(0,0),(0,1),(1,0),(1,1)]:
# ax[a].set_xlim(0.0, np.max(time))
ax[a].set_xlabel('Time (Myr)')
#ax[1].semilogy()
#ax[0].set_ylim(-8, -2)
#ax[1].set_ylim(1.0E-1, 100.0)
ax[(0,0)].set_ylim(-12,-5)
ax[(1,0)].set_ylim(-17,-10)
ax[(0,0)].legend(loc='best')
In [23]:
fig, ax = plt.subplots(2,2)
fig.set_size_inches(16,16)
for i, phase in enumerate(['Molecular','CNM','WNM','WIM','HIM']):
if phase == 'Disk':
continue
pdf_mean = output_stats_data[phase]['O']['num_pdf_mean']
mean = output_stats_data[phase]['O']['mean']
pdf_std = output_stats_data[phase]['O']['num_pdf_std']
std = output_stats_data[phase]['O']['std']
ln_pdf_mean = np.log( pdf_mean / np.sqrt(1.0 + pdf_std**2/pdf_mean**2))
ln_pdf_std = np.sqrt( np.log(1.0 + pdf_std**2/pdf_mean**2))
mean_err = np.abs((mean - ln_pdf_mean)/ln_pdf_mean)
std_err = np.abs((std - ln_pdf_std)/ln_pdf_std)
# mean = np.exp(M + 0.5*S**2)
# std = np.sqrt(np.exp(S**2 + 2.0 * M)*(np.exp(S*S) - 1.0))
#select =
ax[(0,0)].plot(ds_list, mean_err, lw = 5, ls = '-', label = phase)
ax[(0,1)].plot(ds_list, std_err , lw = 5, ls = '-', label = phase)
ax[(0,0)].set_ylabel(r'O $\mu$')
ax[(0,1)].set_ylabel(r'O $\sigma$')
pdf_mean = output_stats_data[phase]['N']['num_pdf_mean']
mean = output_stats_data[phase]['N']['mean']
pdf_std = output_stats_data[phase]['N']['num_pdf_std']
std = output_stats_data[phase]['N']['std']
ln_pdf_mean = np.log( pdf_mean / np.sqrt(1.0 + pdf_std**2/pdf_mean**2))
ln_pdf_std = np.sqrt( np.log(1.0 + pdf_std**2/pdf_mean**2))
mean_err = np.abs((mean - ln_pdf_mean)/ln_pdf_mean)
std_err = np.abs((std - ln_pdf_std)/ln_pdf_std)
# mean = np.exp(M + 0.5*S**2)
# std = np.sqrt(np.exp(S**2 + 2.0 * M)*(np.exp(S*S) - 1.0))
#select =
ax[(1,0)].plot(ds_list, mean_err, lw = 5, ls = '-', label = phase)
ax[(1,1)].plot(ds_list, std_err , lw = 5, ls = '-', label = phase)
ax[(1,0)].set_ylabel(r'N $\mu$')
ax[(1,1)].set_ylabel(r'N $\sigma$')
ax[(0,0)].semilogy()
ax[(0,1)].semilogy()
ax[(1,0)].semilogy()
ax[(1,1)].semilogy()
#ax[(0,1)].set_ylim(0,2)
#ax[(1,1)].set_ylim(0,2)
for a in [(0,0),(0,1),(1,0),(1,1)]:
# ax[a].set_xlim(0.0, np.max(time))
ax[a].set_xlabel('Time (Myr)')
#ax[1].semilogy()
#ax[0].set_ylim(-8, -2)
#ax[1].set_ylim(1.0E-1, 100.0)
#ax[(0,0)].set_ylim(-12,-5)
#ax[(1,0)].set_ylim(-17,-10)
ax[(0,0)].legend(loc='best')
In [24]:
fig, ax = plt.subplots(2,2)
fig.set_size_inches(16,16)
for i, phase in enumerate(['Disk']):# enumerate(['Molecular','CNM','WNM','WIM']):#'WNM','WIM','HIM']):
#if phase == 'Disk':
continue
# mean = np.exp(M + 0.5*S**2)
# std = np.sqrt(np.exp(S**2 + 2.0 * M)*(np.exp(S*S) - 1.0))
y = output_stats_data[phase]['O']['KS_pval']
ax[(0,0)].plot(ds_list, y, lw = 5, ls = '-', label = phase)
#print phase, 'O', output_stats_data[phase]['O']['KS_alpha']
ax[(1,0)].plot(ds_list, y , lw = 5, ls = '-', label = phase)
ax[(0,0)].set_ylabel(r'O D KS')
ax[(1,0)].set_ylabel(r'O D KS')
#mean = output_stats_data[phase]['N']['mean']
#std = output_stats_data[phase]['N']['std']
# mean = np.exp(M + 0.5*S**2)
# std = np.sqrt(np.exp(S**2 + 2.0 * M)*(np.exp(S*S) - 1.0))
#select =
y = output_stats_data[phase]['N']['KS_pval']
ax[(0,1)].plot(ds_list, y, lw = 5, ls = '-', label = phase)
#print phase, 'O', output_stats_data[phase]['N']['KS_alpha']
ax[(1,1)].plot(ds_list, y , lw = 5, ls = '-', label = phase)
ax[(0,1)].set_ylabel(r'N D KS')
ax[(1,1)].set_ylabel(r'N D KS')
#ax[(1,0)].set_ylim(0,2)
#ax[(1,1)].set_ylim(0,2)
for a in [(0,0),(0,1),(1,0),(1,1)]:
ax[a].set_ylim(0.0,1.0)
ax[a].set_xlabel('Time (Myr)')
plt.minorticks_on()
#ax[1].semilogy()
#ax[0].set_ylim(-8, -2)
#ax[1].set_ylim(1.0E-1, 100.0)
ax[(0,0)].legend(loc='best')
Out[24]:
In [20]:
np.log10(10.0)
Out[20]:
In [28]:
fig, ax = plt.subplots(2,3)
fig.set_size_inches(15,10)
pc = {'Molecular' : 'C1', 'CNM' : 'C0', 'WNM' : 'C2', 'WIM':'C4','HIM' : 'C3'}
x = ds_list - ds_list[0]
for i, phase in enumerate(['HIM','WIM','WNM','CNM']):
if phase == 'Disk':
continue
mean = output_stats_data[phase]['O']['actual_mean']
iqr = output_stats_data[phase]['O']['actual_IQR']
median = output_stats_data[phase]['O']['actual_median']
std = output_stats_data[phase]['O']['actual_std']
q90q10 = output_stats_data[phase]['O']['actual_q90q10']
select = std != np.array(None)
ds_list = np.array(ds_list)
# mean = np.exp(M + 0.5*S**2)
# std = np.sqrt(np.exp(S**2 + 2.0 * M)*(np.exp(S*S) - 1.0))
ls = '-'
if phase == 'HIM':
ls = ':'
#select =
#mean = np.log10(mean)
#median = np.log10(median)
#std = np.log10(std)
log_mean = np.log10((mean[select]).astype(np.float64))
ax[(0,0)].plot(x, median, lw = 3, ls = ls, label = phase, color = pc[phase])
ax[(0,1)].plot(x[select], log_mean - median[select], lw = 3, ls = ls, label = phase, color = pc[phase])
ax[(0,2)].plot(x[select], q90q10[select], lw = 3, ls = ls, label = phase, color = pc[phase])
ax[(0,0)].set_ylabel(r'log(Median)')
ax[(0,1)].set_ylabel(r'log(Mean) - log(Median) (dex)')
ax[(0,2)].set_ylabel(r'log(90%) - log(10%) (dex)')
ele2 = 'Ba'
mean = output_stats_data[phase][ele2]['actual_mean']
iqr = output_stats_data[phase][ele2]['actual_IQR']
median = output_stats_data[phase][ele2]['actual_median']
std = output_stats_data[phase][ele2]['actual_std']
q90q10 = output_stats_data[phase][ele2]['actual_q90q10']
select = std != np.array(None)
# mean = np.exp(M + 0.5*S**2)
# std = np.sqrt(np.exp(S**2 + 2.0 * M)*(np.exp(S*S) - 1.0))
#select =
#mean = np.log10(mean)
#median = np.log10(median)
#std = np.log10(std)
log_mean = np.log10((mean[select]).astype(np.float64))
ax[(1,0)].plot(x, median, lw = 3, ls = ls, label = phase,color = pc[phase])
ax[(1,1)].plot(x[select], log_mean - median[select] , lw = 3, ls = ls, label = phase,color = pc[phase])
ax[(1,2)].plot(x[select], q90q10[select], lw = 3, ls = ls, label = phase,color = pc[phase])
ax[(1,0)].set_ylabel(r'log(Median)')
ax[(1,1)].set_ylabel(r'log(Mean) - log(Median) (dex)')
ax[(1,2)].set_ylabel(r'log(90%) - log(10%) (dex)')
plt.minorticks_on()
#ax[(0,2)].set_ylim(0,2)
#ax[(1,2)].set_ylim(0,2)
for a in [(0,0),(0,1),(1,0),(1,1), (0,2),(1,2)]:
ax[a].set_xlim(0.0, 500.0)
ax[a].set_xlabel(r'Time (Myr)')
#for a in [(0,0),(1,0)]:
# ax[a].semilogy()
ax[(0,0)].set_ylim(-8, -2)
ax[(0,1)].set_ylim(0, 1.5)
ax[(0,2)].set_ylim(0.25, 3.0)
xy = (10, -2.5)
ax[(0,0)].annotate('Oxygen',xy,xy)
xy = (10, -10.5)
ax[(1,0)].annotate('Barium',xy,xy)
#ax[(0,2)].set_ylim(1.0E-7, 10**(-1.5))
#ax[1].semilogy()
#ax[0].set_ylim(-8, -2)
#ax[1].set_ylim(1.0E-1, 100.0)
ax[(1,0)].set_ylim( ax[(0,0)].get_ylim())
ax[(1,0)].set_ylim( -16, -10)
ax[(1,1)].set_ylim( ax[(0,1)].get_ylim())
#ax[(1,2)].set_ylim(1.0E-8, 10**(-2.5))
#ax[(0,0)].set_ylim(-12,-5)
ax[(1,2)].set_ylim( ax[(0,2)].get_ylim())
ax[(0,0)].legend(loc='best')
plt.tight_layout()
times = 1.0 * x
fig.savefig('oxygen_barium_phases_CDF_evolution.png')
In [24]:
# convert to an array we can plot as a function of time
output_stats_data = {}
for phase in phases:
output_stats_data[phase] = {}
for e in ele:
output_stats_data[phase][e] = {}
for k in output_stats[ds_list[0]][phase][e].keys():
output_stats_data[phase][e][k] = np.array( [output_stats[x][phase][e][k] for x in ds_list])
In [19]:
output_stats_data['CNM']['O']['mean']
Out[19]:
In [113]:
nrow = 2
ncol = 2
ninch = 8
fig, all_ax = plt.subplots(nrow,ncol)
fig.set_size_inches(ninch*ncol,ninch*nrow)
fbins = data['DD0100']['Disk']['mass_fraction']['bins']
binsize = (fbins[1:]) - (fbins[:-1])
centers = 0.5 * (fbins[1:] + fbins[:-1])
# phase = 'Molecular'
# field = 'O_Fraction'
ylim = [0.0,6]
ds_list = np.arange(80, 551, 10) # [ 100, 150, 200, 250, 300, 350]
ele = ['O','Fe','N','Ba']
phase = 'HIM'
min_median = 1.0E80
max_median = -1.0E80
axi = 0
axj = 0
phases = ['Molecular','CNM','WNM','WIM','HIM','Disk']
phase_labels = ['Molecular','Cold Neutral','Warm Neutral','Warm Ionized','Hot Ionized','All ISM']
mean_phase = {}
std_phase = {}
actual_mean_phase = {}
actual_std_phase = {}
time = np.zeros(np.size(ds_list))
#for k in phases:
# mean_phase[k] = np.zeros(np.size(ds_list))
# std_phase[k] = np.zeros(np.size(ds_list))
# actual_mean_phase[k] = np.zeros(np.size(ds_list))
# actual_std_phase[k] = np.zeros(np.size(ds_list))
def _load_data(ldata, dsname, phase, field, centers = None):
y = ldata[dsname][phase]['mass_fraction'][field]['hist']
mean = ldata[dsname][phase]['mass_fraction'][field]['mean']
std = ldata[dsname][phase]['mass_fraction'][field]['std']
label = ldata[dsname]['general']['Time'] - 46.0
#if (not ('median' in ldata[dsname][phase]['mass_fraction'][field].keys())):
median = np.log10(np.interp(0.5, np.cumsum(y)/(1.0*np.sum(y)), centers))
#else:
# print 'loading median from file'
# median = ldata[dsname][phase]['mass_fraction'][field]['median']
return y, mean, median, std, label
xmin = 10000
xmax = -1000000
plot_median = {}
D_KS_dict = {}
for e in ele:
D_KS_dict[e] = np.zeros(np.size(ds_list))
for i,j in enumerate(ds_list):
dsname = 'DD%0004i'%(j)
axi = 0
axj = 0
for ip, element in enumerate(ele):
# print axi, axj, phase
if nrow == 1:
axindex = axj
else:
axindex = (axi,axj)
ax = all_ax[axindex]
field = element + '_Fraction'
centers = 0.5 * (fbins[1:] + fbins[:-1])
y, mean, median, std, label = _load_data(data, dsname, phase, field, centers = centers)
Mtot = data[dsname][phase]['general']['total_mass']
frac = 0.8
alpha = ((i / (np.size(ds_list)*1.0))) * frac
halpha = ((i / (np.size(ds_list)*1.0))) * frac * 0.5
color = magma(alpha/frac) #'black' # viridis(alpha)
if j == ds_list[-1]:
color = 'black'
if mean is None:
continue
xplot, mean_fit, std_fit, output = _compute_fit_and_plot(ax, fbins, y, mean = mean, median = median, std = std,
alpha = alpha, halpha = halpha, color = color, Mtot=Mtot)
print "%i %s %f %f %f"%(j, element, output['mean'], output['std'], output['D_KS'])
D_KS_dict[element][i] = output['D_KS']
axj = axj + 1
if axj >= ncol:
axi = axi + 1
axj = 0
min_median = np.min([median, min_median])
max_median = np.max([median, max_median])
xmin = np.min([np.log10(xplot[0]),xmin])
xmax = np.max([np.log10(xplot[-1]),xmax])
plot_median[element] = median
for i in np.arange(nrow):
for j in np.arange(ncol):
if nrow == 1:
axindex = j
else:
axindex = (i,j)
all_ax[axindex].set_xlim(xmax - 5, xmax)
all_ax[axindex].semilogy()
all_ax[axindex].set_ylim(1.0E-2,1.0)
ba = 8
ba = 0
# all_ax[axindex].set_xlim(-7 - ba, -2.5 - ba)
# xy = (np.min(all_ax[axindex].get_xlim())+ 0.2,np.max(all_ax[axindex].get_ylim())*0.95)
# all_ax[axindex].annotate(phase_labels[(3)*(i) + j], xy = xy, xytext=xy)
all_ax[axindex].set_ylabel(r'Peak Normalized PDF')
#ax.set_xlim( np.log10(min_median) - 2, np.log10(max_median) + 2)
count = 0
for i in np.arange(nrow):
for j in np.arange(ncol):
if nrow == 1:
axindex = j
else:
axindex = (i,j)
e = ele[count]
all_ax[axindex].set_xlim( plot_median[e] - 3, plot_median[e] + 3)
all_ax[axindex].set_xlabel(r'log(' + e +' Mass Fraction)')
count = count + 1
#ax[(0,0)].set_xlim(
#ax[(0,0)].
#fig.savefig(field + '_phase_evolution.png')
plt.tight_layout()
fig.savefig('lognorm_O_Fe_N_Ba.png')
In [114]:
fig, ax = plt.subplots()
fig.set_size_inches(8,8)
for e in ele:
ax.plot(ds_list, D_KS_dict[e], lw = 3, label = e)
ax.set_xlabel('Time (Myr)')
ax.set_ylabel(r'D$_{\rm CDF}$')
plt.tight_layout()
ax.legend(loc = 'best')
Out[114]:
In [42]:
fig, ax = plt.subplots(2)
fig.set_size_inches(24,16)
ds_num = 400
dsname = 'DD%0004i'%(ds_num)
fbins = data[dsname]['Disk']['mass_fraction']['bins']
binsize = (fbins[1:]) - (fbins[:-1])
centers = 0.5 * (fbins[1:] + fbins[:-1])
# phase = 'Molecular'
# field = 'O_Fraction'
keys = data[dsname]['Disk']['mass_fraction'].keys()
keys = [k.split('_F')[0] for k in keys if '_Fraction' in k]
keys = [k for k in keys if len(k) <= 2]
elements = utilities.sort_by_anum([k for k in keys if not k in ['H','He','H1','H2']])
print elements
elements = ['Ba','Y','As','Sr','Mn','Na','Ca','N','Ni','Mg','S','Si','Fe','C','O']
print elements
def _load_data(ldata, dsname, ephase, field, centers = None):
y = ldata[dsname][phase]['mass_fraction'][field]['hist']
mean = ldata[dsname][phase]['mass_fraction'][field]['mean']
std = ldata[dsname][phase]['mass_fraction'][field]['std']
label = ldata[dsname]['general']['Time'] - 46.0
#if (not ('median' in ldata[dsname][phase]['mass_fraction'][field].keys())):
median = np.log10(np.interp(0.5, np.cumsum(y)/(1.0*np.sum(y)), centers))
#else:
# print 'loading median from file'
# median = ldata[dsname][phase]['mass_fraction'][field]['median']
return y, mean, median, std, label
plot_log = False # True
xmin = 10000
xmax = -1000000
colors = ['C' + str(i) for i in np.arange(9)]
lss = ['-','--']
ci = li = 0
xprev = -1000
plot_median = {}
ymax_order = np.zeros(np.size(elements))
phase = 'Disk'
for i, e in enumerate(elements):
field = e + '_Fraction'
y, mean, median, std, label = _load_data(data, dsname, phase, field, centers = centers)
ymax_order[i] = np.argmax(y/binsize)
sort = np.argsort(ymax_order)
sort_e = np.array(elements)[sort]
label_pos = np.empty((np.size(elements),))
label_pos[::2] = 1
label_pos[1::2] = -1
label_pos = label_pos[sort]
#print sort_e
xychi_y = 1.1
for i, e in enumerate(elements):
dsname = 'DD%0004i'%(ds_num)
field = e + '_Fraction'
centers = 0.5 * (fbins[1:] + fbins[:-1])
y, mean, median, std, label = _load_data(data, dsname, phase, field, centers = centers)
frac = 0.8
alpha = halpha = 1.0
# alpha = ((i / (np.size(elements)*1.0))) * frac
# halpha = ((i / (np.size(elements)*1.0))) * frac * 0.5
# color = magma(i/(1.0*np.size(elements))) #'black' # viridis(alpha)
color = colors[ci]
ls = lss[li]
Mtot = data[dsname][phase]['general']['total_mass']
_compute_fit_and_plot(ax[0], fbins, y,
alpha = alpha, halpha = halpha, color = color, Mtot=Mtot)
xplot, mean_fit, std_fit, output = _compute_fit_and_plot(ax[1], fbins, y, mean = mean, median = median, std = std,
alpha = alpha, halpha = halpha, color = color, show_hist = False, ls = ls, Mtot = Mtot, method = "KS")
#print "%s %i %f %f %f %f"%(e, output['N'], output['mean'], output['std'], output['chisqr'], output['rchisqr'])
print e, output['KS_result'][0], output['KS_result'][1], output['D_KS'], output['KS_pval']
#ax[2].plot(np.log10(output['fit_x']), output['rel_error'], lw = 3, color = color)
#if e == 'Mn':
# print y
if plot_log:
ytext = 4.0
else:
ytext = 1.05
xtext = np.log10(centers[np.argmax(y/binsize)]) - 0.1 - 0.05
xa = np.log10(centers[np.argmax(y/binsize)])
ya = 1.0
pos = label_pos[i]
#if np.abs(xtext -xprev) < 0.2:
if pos > 0:
xtext = xtext + 0.05
if plot_log:
ytext = 2.0
else:
ytext = 1.125
# pos = pos * -1
xy = (xtext, ytext)
xya = (xa,ya)
ax[0].annotate(e, xy = xya, xytext=xy, color = color,
arrowprops=dict(arrowstyle="-", connectionstyle="arc3"))
ax[1].annotate(e, xy = xya, xytext=xy, color = color,
arrowprops=dict(arrowstyle="-",connectionstyle="arc3"))
xychi = (-4, xychi_y)
ax[1].annotate(r'D$_{\rm CDF}$ ' + e + ' = %0.3f'%(output['D_KS']),
xy = xychi, xytext=xychi, color = color)
xychi_y = xychi_y - 0.065 # *10.0**(-0.34)
xprev = xy[0]
ci = ci + 1
if ci >= np.size(colors):
ci = 0
li = li + 1
for i in [0,1]:
ax[i].set_xlim(xmax - 5, xmax)
if plot_log:
ax[i].semilogy()
ax[i].set_ylim(1.0E-5,8.0)
ax[i].set_xlim(-14,-1.5)
ax[i].set_xticks(np.arange(-14,-1))
else:
ax[i].set_ylim(0.0,1.2)
ax[i].set_xlim(-14,-2)
ax[i].set_xticks(np.arange(-14,-1))
ax[i].set_yticks(np.arange(0,11,2)/10.0)
ax[i].set_xlabel(r'log(Mass Fraction)')
#ax[2].set_ylabel('Normalized Error')
#ax[2].set_xlim(ax[1].get_xlim())
#ax[2].semilogy()
ax[0].set_ylabel('Peak Normalized PDF from Data')
ax[1].set_ylabel('Peak Normalized PDF from Fit')
fig.subplots_adjust(hspace=0)
plt.setp([a.get_xticklabels() for a in fig.axes[:-1]], visible=False)
plt.minorticks_on()
#ax[(0,0)].set_xlim(
#ax[(0,0)].
#fig.savefig(field + '_phase_evolution.png')
#plt.tight_layout()
outname = 'lognorm_all_elements_' + phase + '_400.png'
if plot_log:
outname = 'log_' + outname
fig.savefig(outname)
In [14]:
fig, ax = plt.subplots(2)
fig.set_size_inches(24,16)
ds_num = 400
dsname = 'DD%0004i'%(ds_num)
fbins = data[dsname]['Disk']['mass_fraction']['bins']
binsize = (fbins[1:]) - (fbins[:-1])
centers = 0.5 * (fbins[1:] + fbins[:-1])
# phase = 'Molecular'
# field = 'O_Fraction'
keys = data[dsname]['Disk']['mass_fraction'].keys()
keys = [k.split('_F')[0] for k in keys if '_Fraction' in k]
keys = [k for k in keys if len(k) <= 2]
elements = utilities.sort_by_anum([k for k in keys if not k in ['H','He','H1','H2']])
print elements
elements = ['Ba','Y','As','Sr','Mn','Na','Ca','N','Ni','Mg','S','Si','Fe','C','O']
print elements
def _load_data(ldata, dsname, phase, field, centers = None):
y = ldata[dsname][phase]['mass_fraction'][field]['hist']
mean = ldata[dsname][phase]['mass_fraction'][field]['mean']
std = ldata[dsname][phase]['mass_fraction'][field]['std']
label = ldata[dsname]['general']['Time'] - 46.0
#if (not ('median' in ldata[dsname][phase]['mass_fraction'][field].keys())):
median = np.log10(np.interp(0.5, np.cumsum(y)/(1.0*np.sum(y)), centers))
#else:
# print 'loading median from file'
# median = ldata[dsname][phase]['mass_fraction'][field]['median']
return y, mean, median, std, label
plot_log = True
xmin = 10000
xmax = -1000000
colors = ['C' + str(i) for i in np.arange(9)]
lss = ['-','--']
ci = li = 0
xprev = -1000
plot_median = {}
ymax_order = np.zeros(np.size(elements))
phase = 'Molecular'
for i, e in enumerate(elements):
field = e + '_Fraction'
y, mean, median, std, label = _load_data(data, dsname, phase, field, centers = centers)
ymax_order[i] = np.argmax(y/binsize)
sort = np.argsort(ymax_order)
sort_e = np.array(elements)[sort]
label_pos = np.empty((np.size(elements),))
label_pos[::2] = 1
label_pos[1::2] = -1
label_pos = label_pos[sort]
#print sort_e
xychi_y = 1.1
for i, e in enumerate(elements):
dsname = 'DD%0004i'%(ds_num)
field = e + '_Fraction'
centers = 0.5 * (fbins[1:] + fbins[:-1])
y, mean, median, std, label = _load_data(data, dsname, phase, field, centers = centers)
frac = 0.8
alpha = halpha = 1.0
# alpha = ((i / (np.size(elements)*1.0))) * frac
# halpha = ((i / (np.size(elements)*1.0))) * frac * 0.5
# color = magma(i/(1.0*np.size(elements))) #'black' # viridis(alpha)
color = colors[ci]
ls = lss[li]
Mtot = data[dsname][phase]['general']['total_mass']
_compute_fit_and_plot(ax[0], fbins, y,
alpha = alpha, halpha = halpha, color = color, Mtot=Mtot)
xplot, mean_fit, std_fit, output = _compute_fit_and_plot(ax[1], fbins, y, mean = mean, median = median, std = std,
alpha = alpha, halpha = halpha, color = color, show_hist = False, ls = ls, Mtot = Mtot, method = "KS")
#print "%s %i %f %f %f %f"%(e, output['N'], output['mean'], output['std'], output['chisqr'], output['rchisqr'])
#print e, output['KS_result'][0], output['KS_result'][1], output['D_KS'], output['KS_pval']
#ax[2].plot(np.log10(output['fit_x']), output['rel_error'], lw = 3, color = color)
print e, np.log10(output['num_pdf_mean']), np.log10(output['og_mean']), np.log10(output['num_pdf_median'])
#if e == 'Mn':
# print y
if plot_log:
ytext = 4.0
else:
ytext = 1.05
xtext = np.log10(centers[np.argmax(y/binsize)]) - 0.1 - 0.05
xa = np.log10(centers[np.argmax(y/binsize)])
ya = 1.0
pos = label_pos[i]
#if np.abs(xtext -xprev) < 0.2:
if pos > 0:
xtext = xtext + 0.05
if plot_log:
ytext = 2.0
else:
ytext = 1.125
# pos = pos * -1
xy = (xtext, ytext)
xya = (xa,ya)
ax[0].annotate(e, xy = xya, xytext=xy, color = color,
arrowprops=dict(arrowstyle="-", connectionstyle="arc3"))
ax[1].annotate(e, xy = xya, xytext=xy, color = color,
arrowprops=dict(arrowstyle="-",connectionstyle="arc3"))
xychi = (-4, xychi_y)
ax[1].annotate(r'D$_{\rm CDF}$ ' + e + ' = %0.3f'%(output['D_KS']),
xy = xychi, xytext=xychi, color = color)
xychi_y = xychi_y - 0.065 # *10.0**(-0.34)
xprev = xy[0]
ci = ci + 1
if ci >= np.size(colors):
ci = 0
li = li + 1
for i in [0,1]:
ax[i].set_xlim(xmax - 5, xmax)
if plot_log:
ax[i].semilogy()
ax[i].set_ylim(1.0E-5,8.0)
ax[i].set_xlim(-14,-1.5)
ax[i].set_xticks(np.arange(-14,-1))
else:
ax[i].set_ylim(0.0,1.2)
ax[i].set_xlim(-14,-2)
ax[i].set_xticks(np.arange(-14,-1))
ax[i].set_yticks(np.arange(0,11,2)/10.0)
ax[i].set_xlabel(r'log(Mass Fraction)')
#ax[2].set_ylabel('Normalized Error')
#ax[2].set_xlim(ax[1].get_xlim())
#ax[2].semilogy()
ax[0].set_ylabel('Peak Normalized PDF from Data')
ax[1].set_ylabel('Peak Normalized PDF from Fit')
fig.subplots_adjust(hspace=0)
plt.setp([a.get_xticklabels() for a in fig.axes[:-1]], visible=False)
plt.minorticks_on()
#ax[(0,0)].set_xlim(
#ax[(0,0)].
#fig.savefig(field + '_phase_evolution.png')
#plt.tight_layout()
outname = 'loglog_lognorm_all_elements_' + phase + '_400.png'
if plot_log:
outname = 'log_' + outname
fig.savefig(outname)
In [16]:
fig, ax = plt.subplots(2,2, sharey=True) #, sharex=True)
fig.set_size_inches(16,16)
ds_num = 400
dsname = 'DD%0004i'%(ds_num)
fbins = data[dsname]['Disk']['mass_fraction']['bins']
binsize = (fbins[1:]) - (fbins[:-1])
centers = 0.5 * (fbins[1:] + fbins[:-1])
# phase = 'Molecular'
# field = 'O_Fraction'
keys = data[dsname]['Disk']['mass_fraction'].keys()
keys = [k.split('_F')[0] for k in keys if '_Fraction' in k]
keys = [k for k in keys if len(k) <= 2]
elements = utilities.sort_by_anum([k for k in keys if not k in ['H','He','H1','H2']])
print elements
elements = ['Ba','Y','As','Sr','Mn','Na','Ca','N','Ni','Mg','S','Si','Fe','C','O']
print elements
def _load_data(ldata, dsname, phase, field, centers = None):
y = ldata[dsname][phase]['mass_fraction'][field]['hist']
mean = ldata[dsname][phase]['mass_fraction'][field]['mean']
std = ldata[dsname][phase]['mass_fraction'][field]['std']
label = ldata[dsname]['general']['Time'] - 46.0
#if (not ('median' in ldata[dsname][phase]['mass_fraction'][field].keys())):
median = np.log10(np.interp(0.5, np.cumsum(y)/(1.0*np.sum(y)), centers))
#else:
# print 'loading median from file'
# median = ldata[dsname][phase]['mass_fraction'][field]['median']
return y, mean, median, std, label
plot_log = True
xmin = 10000
xmax = -1000000
colors = ['C' + str(i) for i in np.arange(9)]
lss = ['-','--']
ci = li = 0
xprev = -1000
plot_median = {}
ymax_order = np.zeros(np.size(elements))
phase = 'Disk'
for i, e in enumerate(elements):
field = e + '_Fraction'
y, mean, median, std, label = _load_data(data, dsname, phase, field, centers = centers)
ymax_order[i] = np.argmax(y/binsize)
sort = np.argsort(ymax_order)
sort_e = np.array(elements)[sort]
label_pos = np.empty((np.size(elements),))
label_pos[::2] = 1
label_pos[1::2] = -1
label_pos = label_pos[sort]
#print sort_e
xychi_y = 1.1
colors = {'Molecular':'C0','CNM':'C1','WNM':'C2','WIM':'C3','HIM':'C4'}
lss = {'O' : '-', 'N' : '-'}
xs = {'O': 0.0, 'N' : 0.0}
for i, e in enumerate(['N','O']):
x_shift = xs[e]
for phase in ['Molecular','CNM','WNM','WIM','HIM']:
dsname = 'DD%0004i'%(ds_num)
field = e + '_Fraction'
centers = 0.5 * (fbins[1:] + fbins[:-1])
y, mean, median, std, label = _load_data(data, dsname, phase, field, centers = centers)
frac = 0.8
alpha = halpha = 1.0
# alpha = ((i / (np.size(elements)*1.0))) * frac
# halpha = ((i / (np.size(elements)*1.0))) * frac * 0.5
# color = magma(i/(1.0*np.size(elements))) #'black' # viridis(alpha)
color = colors[phase]
ls = lss[e]
Mtot = data[dsname][phase]['general']['total_mass']
_compute_fit_and_plot(ax[(0,i)], fbins, y, x_shift = x_shift,
alpha = alpha, halpha = halpha, color = color, Mtot=Mtot, peak_normalize = False)
xplot, mean_fit, std_fit, output = _compute_fit_and_plot(ax[(1,i)], fbins, y, mean = mean, median = median, std = std,
alpha = alpha, halpha = halpha, color = color, show_hist = False, ls = ls,
Mtot = Mtot, method = "KS", x_shift = x_shift, peak_normalize = False)
#print "%s %i %f %f %f %f"%(e, output['N'], output['mean'], output['std'], output['chisqr'], output['rchisqr'])
print e, output['D_KS'], output['KS_pval'], np.log10(output['num_pdf_mean'])
#ax[2].plot(np.log10(output['fit_x']), output['rel_error'], lw = 3, color = color)
#if e == 'Mn':
# print y
if plot_log:
ytext = 4.0
else:
ytext = 1.05
xtext = np.log10(centers[np.argmax(y/binsize)]) - 0.1 - 0.05
xa = np.log10(centers[np.argmax(y/binsize)])
ya = 1.0
pos = label_pos[i]
#if np.abs(xtext -xprev) < 0.2:
if pos > 0:
xtext = xtext + 0.05
if plot_log:
ytext = 2.0
else:
ytext = 1.125
# pos = pos * -1
xy = (xtext, ytext)
xya = (xa,ya)
# ax[0].annotate(e, xy = xya, xytext=xy, color = color,
# arrowprops=dict(arrowstyle="-", connectionstyle="arc3"))
# ax[1].annotate(e, xy = xya, xytext=xy, color = color,
# arrowprops=dict(arrowstyle="-",connectionstyle="arc3"))
xychi = (-4, xychi_y)
# ax[1].annotate(r'D$_{\rm CDF}$ ' + e + ' = %0.3f'%(output['D_KS']),
# xy = xychi, xytext=xychi, color = color)
xychi_y = xychi_y - 0.065 # *10.0**(-0.34)
xprev = xy[0]
ci = ci + 1
if ci >= np.size(colors):
ci = 0
li = li + 1
for i in [(0,0),(0,1),(1,0),(1,1)]:
ax[i].set_xlim(xmax - 5, xmax)
if plot_log:
ax[i].semilogy()
ax[i].set_ylim(5.0E-2, 2.0E6)
xmin = -9.5
xmax = -0.9
ax[i].set_xlim(xmin,xmax)
ax[i].set_xticks(np.arange(np.floor(xmin),np.ceil(xmax)+0.5))
else:
ax[i].set_ylim(0.0,1.2)
ax[i].set_xlim(-14.1,-1.9)
ax[i].set_xticks(np.arange(-14,-1))
ax[i].set_yticks(np.arange(0,11,2)/10.0)
ax[i].set_xlabel(r'log(Mass Fraction)')
labels = {'Molecular' : 'Molecular', 'CNM' : "Cold Neutral", 'WNM' :"Warm Neutral", "WIM" : "Warm Ionized", "HIM" : 'Hot Ionized'}
for phase in phases:
if phase == 'Disk':
continue
ax[(0,0)].plot([-1000,-999],[1,2], lw = 3, color = colors[phase], ls = '-', label = labels[phase])
#ax[2].set_ylabel('Normalized Error')
#ax[2].set_xlim(ax[1].get_xlim())
#ax[2].semilogy()
#ax[0].set_ylabel('Peak Normalized PDF from Data')
#ax[1].set_ylabel('Peak Normalized PDF from Fit')
#plt.subplots_adjust(wspace=None, hspace=None)
fig.subplots_adjust(wspace = 0, hspace=0)
#plt.setp([a.get_xticklabels() for a in fig.axes[:-1]], visible=False)
plt.minorticks_on()
#ax[(0,0)].set_xlim(
#ax[(0,0)].
ax[(0,0)].legend(loc='best')
#fig.savefig(field + '_phase_evolution.png')
#plt.tight_layout()
outname = 'lognorm_select_elements_all_phases_400.png'
if plot_log:
outname = 'log_' + outname
fig.savefig(outname)
In [49]:
fig, ax = plt.subplots(2,2) #, sharey=True) #, sharex=True)
fig.set_size_inches(12,12)
ds_num = 400
dsname = 'DD%0004i'%(ds_num)
fbins = data[dsname]['Disk']['mass_fraction']['bins']
binsize = (fbins[1:]) - (fbins[:-1])
centers = 0.5 * (fbins[1:] + fbins[:-1])
# phase = 'Molecular'
# field = 'O_Fraction'
keys = data[dsname]['Disk']['mass_fraction'].keys()
keys = [k.split('_F')[0] for k in keys if '_Fraction' in k]
keys = [k for k in keys if len(k) <= 2]
elements = utilities.sort_by_anum([k for k in keys if not k in ['H','He','H1','H2']])
print elements
elements = ['Ba','Y','As','Sr','Mn','Na','Ca','N','Ni','Mg','S','Si','Fe','C','O']
print elements
def _load_data(ldata, dsname, phase, field, centers = None):
y = ldata[dsname][phase]['mass_fraction'][field]['hist']
mean = ldata[dsname][phase]['mass_fraction'][field]['mean']
std = ldata[dsname][phase]['mass_fraction'][field]['std']
label = ldata[dsname]['general']['Time'] - 46.0
#if (not ('median' in ldata[dsname][phase]['mass_fraction'][field].keys())):
median = np.log10(np.interp(0.5, np.cumsum(y)/(1.0*np.sum(y)), centers))
#else:
# print 'loading median from file'
# median = ldata[dsname][phase]['mass_fraction'][field]['median']
return y, mean, median, std, label
plot_log = True
xmin = 10000
xmax = -1000000
colors = ['C' + str(i) for i in np.arange(9)]
lss = ['-','--']
ci = li = 0
xprev = -1000
plot_median = {}
ymax_order = np.zeros(np.size(elements))
phase = 'Disk'
for i, e in enumerate(elements):
field = e + '_Fraction'
y, mean, median, std, label = _load_data(data, dsname, phase, field, centers = centers)
ymax_order[i] = np.argmax(y/binsize)
sort = np.argsort(ymax_order)
sort_e = np.array(elements)[sort]
label_pos = np.empty((np.size(elements),))
label_pos[::2] = 1
label_pos[1::2] = -1
label_pos = label_pos[sort]
#print sort_e
xychi_y = 1.1
colors = {'Molecular':'C0','CNM':'C1','WNM':'C2','WIM':'C4','HIM':'C3'}
lss = {'O' : '-', 'Ba' : '-'}
xs = {'O': 0.0, 'Ba' : 0.0}
for i, e in enumerate(['O','Ba']):
x_shift = xs[e]
for phase in ['Molecular','CNM','WNM','WIM','HIM']:
dsname = 'DD%0004i'%(ds_num)
field = e + '_Fraction'
centers = 0.5 * (fbins[1:] + fbins[:-1])
y, mean, median, std, label = _load_data(data, dsname, phase, field, centers = centers)
frac = 0.8
alpha = halpha = 1.0
# alpha = ((i / (np.size(elements)*1.0))) * frac
# halpha = ((i / (np.size(elements)*1.0))) * frac * 0.5
# color = magma(i/(1.0*np.size(elements))) #'black' # viridis(alpha)
color = colors[phase]
ls = lss[e]
Mtot = data[dsname][phase]['general']['total_mass']
_compute_fit_and_plot(ax[(0,i)], fbins, y, x_shift = x_shift,
alpha = alpha, halpha = halpha, color = color, Mtot=Mtot, peak_normalize = False, plot_CDF=False)
xplot, mean_fit, std_fit, output = _compute_fit_and_plot(ax[(1,i)], fbins, y, mean = mean, median = median, std = std,
alpha = alpha, halpha = halpha, color = color, show_hist = False, ls = ls,
Mtot = Mtot, method = "KS", x_shift = x_shift, peak_normalize = False, plot_CDF = False)
#print "%s %i %f %f %f %f"%(e, output['N'], output['mean'], output['std'], output['chisqr'], output['rchisqr'])
print e, output['D_KS'], output['KS_pval'], np.log10(output['num_pdf_mean'])
#ax[2].plot(np.log10(output['fit_x']), output['rel_error'], lw = 3, color = color)
#if e == 'Mn':
# print y
if plot_log:
ytext = 4.0
else:
ytext = 1.05
xtext = np.log10(centers[np.argmax(y/binsize)]) - 0.1 - 0.05
xa = np.log10(centers[np.argmax(y/binsize)])
ya = 1.0
pos = label_pos[i]
#if np.abs(xtext -xprev) < 0.2:
if pos > 0:
xtext = xtext + 0.05
if plot_log:
ytext = 2.0
else:
ytext = 1.125
# pos = pos * -1
xy = (xtext, ytext)
xya = (xa,ya)
# ax[0].annotate(e, xy = xya, xytext=xy, color = color,
# arrowprops=dict(arrowstyle="-", connectionstyle="arc3"))
# ax[1].annotate(e, xy = xya, xytext=xy, color = color,
# arrowprops=dict(arrowstyle="-",connectionstyle="arc3"))
xychi = (-4, xychi_y)
# ax[1].annotate(r'D$_{\rm CDF}$ ' + e + ' = %0.3f'%(output['D_KS']),
# xy = xychi, xytext=xychi, color = color)
xychi_y = xychi_y - 0.065 # *10.0**(-0.34)
xprev = xy[0]
ci = ci + 1
if ci >= np.size(colors):
ci = 0
li = li + 1
for i in [(0,0),(0,1),(1,0),(1,1)]:
ax[i].set_xlim(xmax - 5, xmax)
# ax[i].set_ylim(, 1.0)
#xmin = -7.5
#xmax = -1.8
#ax[i].set_xlim(xmin,xmax)
##ax[i].set_xticks(np.arange(np.floor(xmin),np.ceil(xmax)+0.5))
ax[i].set_xlabel(r'log(Mass Fraction)')
if plot_log:
ax[i].semilogy()
for i in [0,1]:
ax[(i,0)].set_xlim(-6.5 , -0.8)
ax[(i,1)].set_xlim(-14.5, -8.8)
ax[(i,0)].set_ylim(1.0E-5, 1.0E5)
ax[(i,1)].set_ylim(1.0E3, 1.0E13)
labels = {'Molecular' : 'Molecular', 'CNM' : "Cold Neutral", 'WNM' :"Warm Neutral", "WIM" : "Warm Ionized", "HIM" : 'Hot Ionized'}
for phase in phases:
if phase == 'Disk':
continue
ax[(0,0)].plot([-1000,-999],[1,2], lw = 3, color = colors[phase], ls = '-', label = labels[phase])
ax[(0,1)].plot([-1000,-999],[1,2], lw = 3, color = colors[phase], ls = '-', label = labels[phase])
#ax[2].set_ylabel('Normalized Error')
#ax[2].set_xlim(ax[1].get_xlim())
ax[(0,1)].legend(loc='lower left')
xfrac = 0.8
yfrac = 0.95
for i in [0,1]:
xlim = ax[(i,0)].get_xlim()
ylim = ax[(i,0)].get_ylim()
xy = ( xfrac*(xlim[1]-xlim[0])+xlim[0], 10.0**(yfrac*(np.log10(ylim[1])-np.log10(ylim[0])) + np.log10(ylim[0])) )
ax[(i,0)].annotate('Oxygen', xy, xy)
xlim = ax[(i,1)].get_xlim()
ylim = ax[(i,1)].get_ylim()
xy = ( xfrac*(xlim[1]-xlim[0])+xlim[0], 10.0**(yfrac*(np.log10(ylim[1])-np.log10(ylim[0])) + np.log10(ylim[0])) )
ax[(i,1)].annotate('Barium', xy, xy)
ax[(0,0)].set_ylabel('PDF from Data')
ax[(1,0)].set_ylabel('Log-normal PDF Fit')
ax[(0,1)].set_ylabel('PDF from Data')
ax[(1,1)].set_ylabel('Log-normal PDF Fit')
# plt.subplots_adjust(wspace=None, hspace=None)
# fig.subplots_adjust(wspace = 0, hspace=0)
#plt.setp([a.get_xticklabels() for a in fig.axes[:-1]], visible=False)
plt.minorticks_on()
#ax[(0,0)].set_xlim(
#ax[(0,0)].
plt.tight_layout()
#fig.savefig(field + '_phase_evolution.png')
#plt.tight_layout()
outname = 'lognorm_select_elements_all_phases_400.png'
if plot_log:
outname = 'log_' + outname
fig.savefig(outname)
In [14]:
fig, ax = plt.subplots(2)
fig.set_size_inches(24,16)
fbins = data[dsname]['Disk']['mass_fraction']['bins']
binsize = (fbins[1:]) - (fbins[:-1])
centers = 0.5 * (fbins[1:] + fbins[:-1])
# phase = 'Molecular'
# field = 'O_Fraction'
keys = data[dsname]['Disk']['mass_fraction'].keys()
keys = [k.split('_F')[0] for k in keys if '_Fraction' in k]
keys = [k for k in keys if len(k) <= 2]
elements = utilities.sort_by_anum([k for k in keys if not k in ['H','He','H1','H2']])
print elements
elements = ['Ba','Y','As','Sr','Mn','Na','Ca','N','Ni','Mg','S','Si','Fe','C','O']
print elements
def _load_data(ldata, dsname, phase, field, centers = None):
y = ldata[dsname][phase]['mass_fraction'][field]['hist']
mean = ldata[dsname][phase]['mass_fraction'][field]['mean']
std = ldata[dsname][phase]['mass_fraction'][field]['std']
label = ldata[dsname]['general']['Time'] - 46.0
#if (not ('median' in ldata[dsname][phase]['mass_fraction'][field].keys())):
median = np.log10(np.interp(0.5, np.cumsum(y)/(1.0*np.sum(y)), centers))
#else:
# print 'loading median from file'
# median = ldata[dsname][phase]['mass_fraction'][field]['median']
return y, mean, median, std, label
plot_log = True
xmin = 10000
xmax = -1000000
colors = ['C' + str(i) for i in np.arange(9)]
lss = ['-','--']
ci = li = 0
xprev = -1000
plot_median = {}
ymax_order = np.zeros(np.size(elements))
phase = 'CNM'
for i, e in enumerate(elements):
field = e + '_Fraction'
y, mean, median, std, label = _load_data(data, dsname, phase, field, centers = centers)
ymax_order[i] = np.argmax(y/binsize)
sort = np.argsort(ymax_order)
sort_e = np.array(elements)[sort]
label_pos = np.empty((np.size(elements),))
label_pos[::2] = 1
label_pos[1::2] = -1
label_pos = label_pos[sort]
#print sort_e
for i, e in enumerate(elements):
dsname = 'DD%0004i'%(ds_num)
field = e + '_Fraction'
centers = 0.5 * (fbins[1:] + fbins[:-1])
y, mean, median, std, label = _load_data(data, dsname, phase, field, centers = centers)
frac = 0.8
alpha = halpha = 1.0
# alpha = ((i / (np.size(elements)*1.0))) * frac
# halpha = ((i / (np.size(elements)*1.0))) * frac * 0.5
# color = magma(i/(1.0*np.size(elements))) #'black' # viridis(alpha)
color = colors[ci]
ls = lss[li]
_compute_fit_and_plot(ax[0], fbins, y, plot_PDF = False,
alpha = alpha, halpha = halpha, color = color)
xplot, mean_fit, std_fit, output = _compute_fit_and_plot(ax[1], fbins, y, mean = mean, median = median, std = std, plot_PDF = False,
alpha = alpha, halpha = halpha, color = color, show_hist = False, ls = ls)
if plot_log:
ytext = 4.0
else:
ytext = 1.05
xtext = np.log10(centers[np.argmax(y)]) - 0.1 - 0.05
xa = np.log10(centers[np.argmax(y)])
ya = 1.0
pos = label_pos[i]
#if np.abs(xtext -xprev) < 0.2:
if pos > 0:
xtext = xtext + 0.05
if plot_log:
ytext = 2.0
else:
ytext = 1.125
# pos = pos * -1
xy = (xtext, ytext)
xya = (xa,ya)
ax[0].annotate(e, xy = xya, xytext=xy, color = color,
arrowprops=dict(arrowstyle="-", connectionstyle="arc3"))
ax[1].annotate(e, xy = xya, xytext=xy, color = color,
arrowprops=dict(arrowstyle="-",connectionstyle="arc3"))
xprev = xy[0]
ci = ci + 1
if ci >= np.size(colors):
ci = 0
li = li + 1
for i in [0,1]:
ax[i].set_xlim(xmax - 5, xmax)
if plot_log:
ax[i].semilogy()
ax[i].set_ylim(1.0E-5,8.0)
ax[i].set_xlim(-14,-1.5)
ax[i].set_xticks(np.arange(-14,-1))
else:
ax[i].set_ylim(0.0,0.15)
ax[i].set_xlim(-14,-2)
ax[i].set_xticks(np.arange(-14,-1))
ax[i].set_yticks(np.arange(0,1.5,0.2)/10.0)
ax[i].set_xlabel(r'log(Mass Fraction)')
ax[0].set_ylabel('Fraction of Mass from Data')
ax[1].set_ylabel('Fraction of Mass from Fit')
fig.subplots_adjust(hspace=0)
plt.setp([a.get_xticklabels() for a in fig.axes[:-1]], visible=False)
plt.minorticks_on()
#ax[(0,0)].set_xlim(
#ax[(0,0)].
#fig.savefig(field + '_phase_evolution.png')
#plt.tight_layout()
outname = 'lognorm_int_all_elements_' + phase + '.png'
if plot_log:
outname = 'log_' + outname
fig.savefig(outname)
In [18]:
fig, ax = plt.subplots(2)
fig.set_size_inches(24,16)
ds_num = 300
dsname = 'DD%0004i'%(ds_num)
fbins = data[dsname]['Disk']['mass_fraction']['bins']
binsize = (fbins[1:]) - (fbins[:-1])
centers = 0.5 * (fbins[1:] + fbins[:-1])
# phase = 'Molecular'
# field = 'O_Fraction'
keys = data[dsname]['Disk']['mass_fraction'].keys()
keys = [k.split('_F')[0] for k in keys if '_Fraction' in k]
keys = [k for k in keys if len(k) <= 2]
elements = utilities.sort_by_anum([k for k in keys if not k in ['H','He','H1','H2']])
print elements
elements = ['Ba','Y','As','Sr','Mn','Na','Ca','N','Ni','Mg','S','Si','Fe','C','O']
print elements
def _load_data(ldata, dsname, phase, field, centers = None):
y = ldata[dsname][phase]['mass_fraction'][field]['hist']
mean = ldata[dsname][phase]['mass_fraction'][field]['mean']
std = ldata[dsname][phase]['mass_fraction'][field]['std']
label = ldata[dsname]['general']['Time'] - 46.0
#if (not ('median' in ldata[dsname][phase]['mass_fraction'][field].keys())):
median = np.log10(np.interp(0.5, np.cumsum(y)/(1.0*np.sum(y)), centers))
#else:
# print 'loading median from file'
# median = ldata[dsname][phase]['mass_fraction'][field]['median']
return y, mean, median, std, label
plot_log = False
xmin = 10000
xmax = -1000000
colors = ['C' + str(i) for i in np.arange(9)]
lss = ['-','--']
ci = li = 0
xprev = -1000
plot_median = {}
ymax_order = np.zeros(np.size(elements))
phase = 'HIM'
for i, e in enumerate(elements):
field = e + '_Fraction'
y, mean, median, std, label = _load_data(data, dsname, phase, field, centers = centers)
ymax_order[i] = np.argmax(y/binsize)
sort = np.argsort(ymax_order)
sort_e = np.array(elements)[sort]
label_pos = np.empty((np.size(elements),))
label_pos[::2] = 1
label_pos[1::2] = -1
label_pos = label_pos[sort]
#print sort_e
xychi_y = 0.9
for i, e in enumerate(elements):
dsname = 'DD%0004i'%(ds_num)
field = e + '_Fraction'
centers = 0.5 * (fbins[1:] + fbins[:-1])
y, mean, median, std, label = _load_data(data, dsname, phase, field, centers = centers)
frac = 0.8
alpha = halpha = 1.0
# alpha = ((i / (np.size(elements)*1.0))) * frac
# halpha = ((i / (np.size(elements)*1.0))) * frac * 0.5
# color = magma(i/(1.0*np.size(elements))) #'black' # viridis(alpha)
color = colors[ci]
ls = lss[li]
_compute_fit_and_plot(ax[0], fbins, y, plot_CDF = True,
alpha = alpha, halpha = halpha, color = color)
xplot, mean_fit, std_fit, output = _compute_fit_and_plot(ax[1], fbins, y, mean = mean, median = median, std = std, plot_CDF = True,
alpha = alpha, halpha = halpha, color = color, show_hist = False, ls = ls, method = 'KS')
p = np.exp( -7.01256*output['D_KS']**2 * (output['npoints'] + 2.78019) + 2.99587*output['D_KS']*(output['npoints']+2.78019)**(0.5) - 0.122119 + 0.97498/output['npoints']**(0.5) + 1.67997/output['npoints'] )
print e, output['D_KS'], output['mean'], output['std'], output['KS_pval']
if plot_log:
ytext = 4.0
else:
ytext = 1.05
xtext = np.log10(centers[np.argmax(y)]) - 0.1 - 0.05
xa = np.log10(centers[np.argmax(y)])
ya = 1.0
pos = label_pos[i]
#if np.abs(xtext -xprev) < 0.2:
if pos > 0:
xtext = xtext + 0.05
if plot_log:
ytext = 2.0
else:
ytext = 1.125
# pos = pos * -1
xy = (xtext, ytext)
xya = (xa,ya)
ax[0].annotate(e, xy = xya, xytext=xy, color = color,
arrowprops=dict(arrowstyle="-", connectionstyle="arc3"))
ax[1].annotate(e, xy = xya, xytext=xy, color = color,
arrowprops=dict(arrowstyle="-",connectionstyle="arc3"))
xprev = xy[0]
xychi = (-4, xychi_y)
xychi_y = xychi_y - 0.05
ax[1].annotate(r'p$_{\rm KS}$ %2s = %0.3f'%(e,output['KS_pval']),
xy = xychi, xytext=xychi, color = color)
ci = ci + 1
if ci >= np.size(colors):
ci = 0
li = li + 1
for i in [0,1]:
ax[i].set_xlim(xmax - 5, xmax)
if plot_log:
ax[i].semilogy()
ax[i].set_ylim(1.0E-5,8.0)
ax[i].set_xlim(-14,-1.5)
ax[i].set_xticks(np.arange(-14,-1))
else:
ax[i].set_ylim(0.0,1.2)
ax[i].set_xlim(-14,-2)
ax[i].set_xticks(np.arange(-14,-1))
ax[i].set_yticks(np.arange(0,11,2)/10.0)
ax[i].set_xlabel(r'log(Mass Fraction)')
ax[0].set_ylabel('Fraction of Mass from Data')
ax[1].set_ylabel('Fraction of Mass from Fit')
fig.subplots_adjust(hspace=0)
plt.setp([a.get_xticklabels() for a in fig.axes[:-1]], visible=False)
plt.minorticks_on()
#ax[(0,0)].set_xlim(
#ax[(0,0)].
#fig.savefig(field + '_phase_evolution.png')
#plt.tight_layout()
outname = 'lognorm_CDF_all_elements_' + phase + '.png'
if plot_log:
outname = 'log_' + outname
fig.savefig(outname)
In [16]:
nrow = 2
ncol = 3
fig, all_ax = plt.subplots(nrow,ncol)
fig.set_size_inches(8*ncol,8*nrow)
fbins = data['DD0100']['Disk']['mass_fraction']['bins']
binsize = (fbins[1:]) - (fbins[:-1])
centers = 0.5 * (fbins[1:] + fbins[:-1])
phase = 'Molecular'
field = 'O_Fraction'
ylim = [0.0,6]
ds_list = np.arange(60, 120, 5)
min_median = 1.0E80
max_median = -1.0E80
axi = 0
axj = 0
phases = ['Molecular','CNM','WNM','WIM','HIM','Disk']
mean_phase = {}
std_phase = {}
actual_mean_phase = {}
actual_std_phase = {}
time = np.zeros(np.size(ds_list))
for k in phases:
mean_phase[k] = np.zeros(np.size(ds_list))
std_phase[k] = np.zeros(np.size(ds_list))
actual_mean_phase[k] = np.zeros(np.size(ds_list))
actual_std_phase[k] = np.zeros(np.size(ds_list))
xmin = 10000
xmax = -1000000
for i,j in enumerate(ds_list):
dsname = 'DD%0004i'%(j)
axi = 0
axj = 0
for ip, phase in enumerate(phases):
# print axi, axj, phase
ax = all_ax[(axi, axj)]
y = data[dsname][phase]['mass_fraction'][field]['hist']
mean = data[dsname][phase]['mass_fraction'][field]['mean']
std = data[dsname][phase]['mass_fraction'][field]['std']
median = np.log10(np.interp(0.5, np.cumsum(y)/(1.0*np.sum(y)), centers))
label = data[dsname]['general']['Time'] - 46.0
# centers = 0.5 * (fbins[1:] + fbins[:-1])
# y, mean, median, std, label = _load_data(data, dsname, phase, field, centers = centers)
min_median = np.min([median, min_median])
max_median = np.max([median, max_median])
frac = 0.8
alpha = ((i / (np.size(ds_list)*1.0))) * frac
halpha = ((i / (np.size(ds_list)*1.0))) * frac * 0.5
color = magma(alpha/frac) #'black' # viridis(alpha)
lw = 4.0
hlw = 4.0
if j == np.max(ds_list):
color = 'black'
alpha = 1.0
lw = 5.0
hlw = 5.0
norm_y = y / binsize
plot_bins = np.log10(fbins)
norm = np.max(norm_y)
norm = 1.0
plot_histogram(ax, plot_bins, norm_y/norm, lw = hlw, color = color, alpha = halpha)
selection = (y > 0) * (centers > 10**(median - 4)) * (centers < 10**(median+4))
fit_x = centers[selection]
y_to_fit = norm_y[selection]
# take initial guess parameters from distribution values - compute logged dist values
u_guess = np.log( mean / (np.sqrt(1.0 + std*std/(mean*mean))))
std_guess = np.sqrt(np.log(1.0 + std*std/(mean*mean)))
xplot = np.logspace(np.log10(np.min(fit_x)), np.log10(np.max(fit_x)),4000)
try:
lognorm = functions.lognormal()
popt, pcov = lognorm.fit_function(fit_x, y_to_fit, p0 = [u_guess, std_guess])
xplot = np.logspace(np.log10(np.min(fit_x)), np.log10(np.max(fit_x)),4000)
yplot = lognorm._f(xplot, *popt)
norm = np.max(yplot)
norm = 1.0
ax.plot(np.log10(xplot), yplot/norm, lw = lw, color = color , ls = '-', alpha = alpha) #, label = 'Fit')
except:
print phase, dsname, 'failure'
popt = [None, None]
mean_phase[phase][i] = popt[0]
actual_mean_phase[phase][i] = mean
std_phase[phase][i] = popt[1]
actual_std_phase[phase][i] = std
time[i] = label
axj = axj + 1
if axj >= ncol:
axi = axi + 1
axj = 0
xmin = np.min([np.log10(xplot[0]),xmin])
xmax = np.max([np.log10(xplot[-1]),xmax])
for i in np.arange(2):
for j in np.arange(3):
all_ax[(i,j)].set_xlim(xmax - 5, xmax)
#all_ax[(i,j)].set_ylim(1.0E-4,1.0)
all_ax[(i,j)].semilogy()
ba = 8
ba = 0
all_ax[(i,j)].set_xlim(-10 - ba, -2.5 - ba)
xy = (np.min(all_ax[(i,j)].get_xlim())+ 0.2,np.max(all_ax[(i,j)].get_ylim())*0.95)
all_ax[(i,j)].annotate(phases[(3)*(i) + j], xy = xy, xytext=xy)
all_ax[(i,j)].set_xlabel(r'log(Mass Fraction)')
all_ax[(i,j)].set_ylabel(r'Peak Normalized PDF')
#ax.set_xlim( np.log10(min_median) - 2, np.log10(max_median) + 2)
fig.savefig(field + '_phase_evolution.png')
In [ ]:
In [17]:
time
Out[17]:
In [80]:
fig, ax = plt.subplots(1,2)
fig.set_size_inches(16,8)
for i, phase in enumerate(phases):
if phase == 'Disk':
continue
mean = actual_mean_phase[phase]
std = actual_std_phase[phase]
# mean = np.exp(M + 0.5*S**2)
# std = np.sqrt(np.exp(S**2 + 2.0 * M)*(np.exp(S*S) - 1.0))
#select =
ax[0].plot(time, np.log10(mean), lw = 5, ls = '-')
ax[1].plot(time, std / mean, lw = 5, ls = '-', label = phase)
ax[0].set_ylabel(r'Mean')
ax[1].set_ylabel(r'Coefficient of Variation')
for a in [ax[0],ax[1]]:
a.set_xlim(0.0, np.max(time))
a.set_xlabel('Time (Myr)')
ax[1].semilogy()
ax[0].set_ylim(-8, -2)
ax[1].set_ylim(1.0E-1, 100.0)
ax[1].legend(loc='best')
fig.savefig('phase_mean_var_from_distribution.png')
In [39]:
fig, ax = plt.subplots(1,2)
fig.set_size_inches(16,8)
for i, phase in enumerate(phases):
if phase == 'Disk':
continue
M = mean_phase[phase]
S = std_phase[phase]
mean = np.exp(M + 0.5*S**2)
std = np.sqrt(np.exp(S**2 + 2.0 * M)*(np.exp(S*S) - 1.0))
#select =
ax[0].plot(time, mean, lw = 5, ls = '-')
ax[1].plot(time, std / mean, lw = 5, ls = '-', label = phase)
ax[0].set_ylabel(r'Mean (from fit)')
ax[1].set_ylabel(r'Coefficient of Variation (from fit)')
for a in [ax[0],ax[1]]:
a.set_xlim(0.0, np.max(time))
a.set_xlabel('Time (Myr)')
a.semilogy()
ax[0].set_ylim(1.0E-6, 1.0E-2)
ax[1].set_ylim(0.1, 100.0)
ax[1].legend(loc='best')
fig.savefig('phase_mean_var_from_fit.png')
In [ ]:
In [17]:
ptimes = ds_list # in Myr
dt = 100.0
bins = np.arange(-15,-2, 0.2)
all_hist = {}
for i, t in enumerate(ptimes):
selection = ptypes.select_formed_stars(gal.ds, gal.df, np.max([t - dt, 0.0]),
np.min([t + dt, np.max(ptimes)]))
ydata = np.log10(gal.df['particle_O_fraction'][selection])
hist, bins = np.histogram(ydata, bins = bins)
hist = hist / (10**(bins[1:]) - 10**(bins[:-1])) / (1.0 * np.size(ydata)) # normalize to PDF
all_hist[t] = hist
In [26]:
Out[26]:
In [39]:
t_o = gal.df['creation_time'].convert_to_units('Myr')
xdata = times*1.0
fig, all_axes = plt.subplots(2,3)
fig.set_size_inches(24,16)
axi = axj = 0
for i,phase in enumerate(['CNM','WNM','WIM','HIM','Disk']):
ax = all_axes[(axi,axj)]
xdata = times + 119.0
mean = actual_mean_phase[phase]
std = actual_std_phase[phase]
median = actual_median_phase[phase]
#mean = np.exp(mol_mean + 0.5*mol_std**2)
#std = np.sqrt(np.exp(mol_std**2 + 2.0 * mol_mean)*(np.exp(mol_std*mol_std) - 1.0))
xdata = xdata[std>0]
mean = mean[std>0]
median = median[std>0]
std = std[std>0]
# number of standard deviations away from mean for each star
sigma_stars = (np.log10(gal.df['particle_O_fraction']) - np.interp(t_o, xdata, median)) #/ np.interp(t_o, xdata, std)
x = gal.df[('io','particle_Fe_over_H')]
ax.scatter(t_o - np.min(t_o), sigma_stars, alpha = 0.5, color = 'black')
ax.set_xlim(np.min(t_o-np.min(t_o)), np.max(t_o-np.min(t_o)))
#ax.set_xlim(-8,0)
ax.set_ylim(-2,2)
ax.plot(ax.get_xlim(), [0.0,0.0], color = 'black', lw = 4, ls = '--')
#diff = median - np.log10(mean)
#print diff
#ax.plot(xdata-np.min(t_o.value), diff, color = 'red', lw = 3, ls = '--')
ax.set_ylabel('Distance From Median O Fraction (dex)')
ax.set_xlabel('Time (Myr)')
xy = (10.0, ax.get_ylim()[1] - 0.2)
ax.annotate(phase_labels[i], xy=xy,xytext=xy)
axj = axj + 1
if axj >= 3:
axj = 0
axi = axi + 1
plt.tight_layout()
fig.savefig('stellar_distance_to_median.png')
In [40]:
nrow = 2
ncol = 3
fig, all_ax = plt.subplots(nrow,ncol)
fig.set_size_inches(8*ncol,8*nrow)
fbins = data['DD0200']['Disk']['mass_fraction']['bins']
binsize = (fbins[1:]) - (fbins[:-1])
centers = 0.5 * (fbins[1:] + fbins[:-1])
phase = 'Molecular'
field = 'Ba_Fraction'
ylim = [0.0,6]
ds_list = np.arange(120.0, 401, 10)
min_median = 1.0E80
max_median = -1.0E80
axi = 0
axj = 0
phases = ['CNM','WNM','WIM','HIM','Disk']
phase_labels = ['Cold Neutral','Warm Neutral','Warm Ionized','Hot Ionized','All ISM']
mean_phase = {}
std_phase = {}
actual_mean_phase = {}
actual_median_phase = {}
actual_std_phase = {}
time = np.zeros(np.size(ds_list))
for k in phases:
mean_phase[k] = np.zeros(np.size(ds_list))
std_phase[k] = np.zeros(np.size(ds_list))
actual_mean_phase[k] = np.zeros(np.size(ds_list))
actual_median_phase[k] = np.zeros(np.size(ds_list))
actual_std_phase[k] = np.zeros(np.size(ds_list))
xmin = 10000
xmax = -1000000
for i,j in enumerate(ds_list):
dsname = 'DD%0004i'%(j)
axi = 0
axj = 0
for ip, phase in enumerate(phases):
# print axi, axj, phase
ax = all_ax[(axi, axj)]
y = data[dsname][phase]['mass_fraction'][field]['hist']
mean = data[dsname][phase]['mass_fraction'][field]['mean']
std = data[dsname][phase]['mass_fraction'][field]['std']
median = np.log10(np.interp(0.5, np.cumsum(y)/(1.0*np.sum(y)), centers))
label = data[dsname]['general']['Time'] - 46.0
min_median = np.min([median, min_median])
max_median = np.max([median, max_median])
frac = 0.8
alpha = ((i / (np.size(ds_list)*1.0))) * frac
halpha = ((i / (np.size(ds_list)*1.0))) * frac * 0.5
color = magma(alpha/frac) #'black' # viridis(alpha)
lw = 4.0
hlw = 4.0
if j == np.max(ds_list):
color = 'black'
alpha = 1.0
lw = 5.0
hlw = 5.0
norm_y = y / binsize
plot_bins = np.log10(fbins)
plot_histogram(ax, plot_bins, norm_y/np.max(norm_y), lw = hlw, color = color, alpha = halpha)
selection = (y > 0) * (centers > 10**(median - 4)) * (centers < 10**(median+4))
fit_x = centers[selection]
y_to_fit = norm_y[selection]
# take initial guess parameters from distribution values - compute logged dist values
try:
u_guess = np.log( mean / (np.sqrt(1.0 + std*std/(mean*mean))))
std_guess = np.sqrt(np.log(1.0 + std*std/(mean*mean)))
xplot = np.logspace(np.log10(np.min(fit_x)), np.log10(np.max(fit_x)),4000)
lognorm = functions.lognormal()
popt, pcov = lognorm.fit_function(fit_x, y_to_fit, p0 = [u_guess, std_guess])
xplot = np.logspace(np.log10(np.min(fit_x)), np.log10(np.max(fit_x)),4000)
yplot = lognorm._f(xplot, *popt)
ax.plot(np.log10(xplot), yplot/np.max(yplot), lw = lw, color = color , ls = '-', alpha = alpha) #, label = 'Fit')
except:
print phase, dsname, 'failure'
popt = [None, None]
mean_phase[phase][i] = popt[0]
actual_mean_phase[phase][i] = mean
actual_median_phase[phase][i] = median
std_phase[phase][i] = popt[1]
actual_std_phase[phase][i] = std
time[i] = label
axj = axj + 1
if axj >= ncol:
axi = axi + 1
axj = 0
xmin = np.min([np.log10(xplot[0]),xmin])
xmax = np.max([np.log10(xplot[-1]),xmax])
for i in np.arange(2):
for j in np.arange(3):
all_ax[(i,j)].set_xlim(xmax - 5, xmax)
all_ax[(i,j)].semilogy()
all_ax[(i,j)].set_ylim(1.0E-2,1.0)
ba = 8
ba = 0
all_ax[(i,j)].set_xlim(-7 - ba, -2.5 - ba)
xy = (np.min(all_ax[(i,j)].get_xlim())+ 0.2,np.max(all_ax[(i,j)].get_ylim())*0.95)
all_ax[(i,j)].annotate(phase_labels[(3)*(i) + j], xy = xy, xytext=xy)
all_ax[(i,j)].set_xlabel(r'log(Oxygen Mass Fraction)')
all_ax[(i,j)].set_ylabel(r'Peak Normalized PDF')
#ax.set_xlim( np.log10(min_median) - 2, np.log10(max_median) + 2)
#fig.savefig(field + '_phase_evolution.png')
fig.savefig('log_panel_phase_evolution.png')
In [41]:
t_o = gal.df['creation_time'].convert_to_units('Myr')
xdata = times*1.0
fig, all_axes = plt.subplots(2,3)
fig.set_size_inches(24,16)
axi = axj = 0
for i,phase in enumerate(['CNM','WNM','WIM','HIM','Disk']):
ax = all_axes[(axi,axj)]
xdata = times + 119.0
mean = actual_mean_phase[phase]
std = actual_std_phase[phase]
median = actual_median_phase[phase]
#mean = np.exp(mol_mean + 0.5*mol_std**2)
#std = np.sqrt(np.exp(mol_std**2 + 2.0 * mol_mean)*(np.exp(mol_std*mol_std) - 1.0))
xdata = xdata[std>0]
mean = mean[std>0]
median = median[std>0]
std = std[std>0]
# number of standard deviations away from mean for each star
sigma_stars = (np.log10(gal.df['particle_Ba_fraction']) - np.interp(t_o, xdata, median)) #/ np.interp(t_o, xdata, std)
x = gal.df[('io','particle_Fe_over_H')]
ax.scatter(t_o - np.min(t_o), sigma_stars, alpha = 0.5, color = 'black')
ax.set_xlim(np.min(t_o-np.min(t_o)), np.max(t_o-np.min(t_o)))
#ax.set_xlim(-8,0)
ax.set_ylim(-2,2)
ax.plot(ax.get_xlim(), [0.0,0.0], color = 'black', lw = 4, ls = '--')
#diff = median - np.log10(mean)
#print diff
#ax.plot(xdata-np.min(t_o.value), diff, color = 'red', lw = 3, ls = '--')
ax.set_ylabel('Distance From Median Ba Fraction (dex)')
ax.set_xlabel('Time (Myr)')
xy = (10.0, ax.get_ylim()[1] - 0.2)
ax.annotate(phase_labels[i], xy=xy,xytext=xy)
axj = axj + 1
if axj >= 3:
axj = 0
axi = axi + 1
plt.tight_layout()
fig.savefig('stellar_distance_to_median_Ba.png')
In [ ]:
fig, ax = plt.subplots()
fig.set_size_inches(8,8)
for i in [150, 200, 250, 300]:
plot_histogram(ax, bins, all_hist[i] , lw = 3)
# print i, all_hist[i]
ax.set_xlabel('log(O Fraction)')
ax.set_ylabel('PDF')
ax.set_xlim(-8,-3)
ax.semilogy()
In [ ]:
lognorm?
In [10]:
In [ ]: