diff --git a/examples/plot.py b/examples/plot.py index c59569cf99aec2fb4d2521ed4406d35baeb779f2..941a33ff03841be073f9f6f800f15faf2786bb98 100644 --- a/examples/plot.py +++ b/examples/plot.py @@ -66,19 +66,19 @@ accx_new = accx_new[rank] accy_new = accy_new[rank] accz_new = accz_new[rank] -# Read Gadget accelerations -data=loadtxt("gadget_dump.dat") -id = data[:,0] -accx_g=data[:,5] -accy_g=data[:,6] -accz_g=data[:,7] +# # Read Gadget accelerations +# data=loadtxt("gadget_dump.dat") +# id = data[:,0] +# accx_g=data[:,5] +# accy_g=data[:,6] +# accz_g=data[:,7] -# Sort accelerations -rank = argsort(id) -id = id[rank] -accx_g = accx_g[rank] -accy_g = accy_g[rank] -accz_g = accz_g[rank] +# # Sort accelerations +# rank = argsort(id) +# id = id[rank] +# accx_g = accx_g[rank] +# accy_g = accy_g[rank] +# accz_g = accz_g[rank] # Build error ------------------------------------------------ @@ -91,9 +91,9 @@ errx_new = (accx_new - accx_e )/abs(accx_e) erry_new = (accy_new - accy_e )/abs(accy_e) errz_new = (accz_new - accz_e )/abs(accz_e) -errx_g = (accx_g - accx_e )/abs(accx_e) -erry_g = (accy_g - accy_e )/abs(accy_e) -errz_g = (accz_g - accz_e )/abs(accz_e) +# errx_g = (accx_g - accx_e )/abs(accx_e) +# erry_g = (accy_g - accy_e )/abs(accy_e) +# errz_g = (accz_g - accz_e )/abs(accz_e) # Statistics meanx_bh = mean(errx_bh[abs(errx_bh) < 0.1]) @@ -121,7 +121,7 @@ stdz_new = std(errz_new[abs(errz_new) < 0.1]) figure(frameon=True) subplot(311, title="Acceleration along X") -plot(id, errx_g , 'gs') +#plot(id, errx_g , 'gs') plot(id, errx_bh , 'rx') plot(id, errx_new , 'b.') text(id[-1], 0.18, "B-H: $%5.3f\\pm%5.3f$\n QuickShed: $%5.3f\\pm%5.3f$"%(meanx_bh, stdx_bh, meanx_new, stdx_new), backgroundcolor="w", va="top", ha="right" ) @@ -132,7 +132,7 @@ xlim(0,id[-1]) grid() subplot(312, title="Acceleration along Y") -plot(id, erry_g , 'gs') +#plot(id, erry_g , 'gs') plot(id, erry_bh , 'rx') plot(id, erry_new , 'b.') text(id[-1], 0.18, "B-H: $%5.3f\\pm%5.3f$\n QuickShed: $%5.3f\\pm%5.3f$"%(meany_bh, stdy_bh, meany_new, stdy_new), backgroundcolor="w", va="top", ha="right" ) @@ -144,7 +144,7 @@ xlim(0,id[-1]) grid() subplot(313, title="Acceleration along Z") -plot(id, errz_g , 'gs', label="Gadget") +#plot(id, errz_g , 'gs', label="Gadget") plot(id, errz_new , 'b.', label="QuickShed") plot(id, errz_bh , 'rx', label="Legacy") #text(id[-1], 0.18, "B-H: $%5.3f\\pm%5.3f$\n QuickShed: $%5.3f\\pm%5.3f$"%(meanz_bh, stdz_bh, meanz_new, stdz_new), backgroundcolor="w", va="top", ha="right" ) @@ -165,20 +165,20 @@ bins = linspace(-3, 3, 10000) figure(frameon=True) subplot(311, title="Acceleration along X")#, yscale='log') -hist(errx_g, bins=bins, normed=1, histtype='step', rwidth=0.01, color='g', label="Gadget") +#hist(errx_g, bins=bins, normed=1, histtype='step', rwidth=0.01, color='g', label="Gadget") hist(errx_bh, bins=bins, normed=1, histtype='step', rwidth=0.01, color='r', label="Legacy") hist(errx_new, bins=bins, normed=1, histtype='step', rwidth=0.01, color='b', label="QuickShed") legend(loc="upper right") xlim(-0.03, 0.03) subplot(312, title="Acceleration along Y") -hist(erry_g, bins=bins, normed=1, histtype='step', rwidth=0.01, color='g') +#hist(erry_g, bins=bins, normed=1, histtype='step', rwidth=0.01, color='g') hist(erry_bh, bins=bins, normed=1, histtype='step', rwidth=0.01, color='r') hist(erry_new, bins=bins, normed=1, histtype='step', rwidth=0.01, color='b') xlim(-0.03, 0.03) subplot(313, title="Acceleration along Z") -hist(errz_g, bins=bins, normed=1, histtype='step', rwidth=0.01, color='g') +#hist(errz_g, bins=bins, normed=1, histtype='step', rwidth=0.01, color='g') hist(errz_bh, bins=bins, normed=1, histtype='step', rwidth=0.01, color='r') hist(errz_new, bins=bins, normed=1, histtype='step', rwidth=0.01, color='b') xlim(-0.03, 0.03) @@ -187,9 +187,9 @@ savefig("histogram.png") # Print error -errx_g = errx_g[ (errx_g > -0.1) & (errx_g < 0.1) ] -erry_g = erry_g[ (erry_g > -0.1) & (erry_g < 0.1) ] -errz_g = errz_g[ (errz_g > -0.1) & (errz_g < 0.1) ] +# errx_g = errx_g[ (errx_g > -0.1) & (errx_g < 0.1) ] +# erry_g = erry_g[ (erry_g > -0.1) & (erry_g < 0.1) ] +# errz_g = errz_g[ (errz_g > -0.1) & (errz_g < 0.1) ] errx_bh = errx_bh[ (errx_bh > -0.1) & (errx_bh < 0.1) ] erry_bh = erry_bh[ (erry_bh > -0.1) & (erry_bh < 0.1) ] @@ -200,6 +200,6 @@ erry_new = erry_new[ (erry_new > -0.1) & (erry_new < 0.1) ] errz_new = errz_new[ (errz_new > -0.1) & (errz_new < 0.1) ] -print "Error for Gadget: ( x= %4.3e"%std(errx_g), "y= %4.3e"%std(erry_g), "z= %4.3e"%std(errz_g), ") Combined= %4.3e"%(( std(errx_g) + std(erry_g) + std(errz_g) )/3.) +#print "Error for Gadget: ( x= %4.3e"%std(errx_g), "y= %4.3e"%std(erry_g), "z= %4.3e"%std(errz_g), ") Combined= %4.3e"%(( std(errx_g) + std(erry_g) + std(errz_g) )/3.) print "Error for Legacy: ( x= %4.3e"%std(errx_bh), "y= %4.3e"%std(erry_bh), "z= %4.3e"%std(errz_bh), ") Combined= %4.3e"%(( std(errx_bh) + std(erry_bh) + std(errz_bh) )/3.) print "Error for QuShed: ( x= %4.3e"%std(errx_new), "y= %4.3e"%std(erry_new), "z= %4.3e"%std(errz_new), ") Combined= %4.3e"%(( std(errx_new) + std(erry_new) + std(errz_new) )/3.) diff --git a/examples/plot_gadget.py b/examples/plot_gadget.py new file mode 100644 index 0000000000000000000000000000000000000000..c59569cf99aec2fb4d2521ed4406d35baeb779f2 --- /dev/null +++ b/examples/plot_gadget.py @@ -0,0 +1,205 @@ +#!/usr/bin/env python2 +# -*- coding: utf-8 -*- +import matplotlib +matplotlib.use("Agg") +from pylab import * +import numpy +# tex stuff +#rc('font',**{'family':'serif','serif':['Palatino']}) +params = {'axes.labelsize': 16, + 'axes.titlesize': 20, + 'text.fontsize': 14, + 'legend.fontsize': 14, + 'xtick.labelsize': 14, + 'ytick.labelsize': 14, + 'text.usetex': True, + + 'figure.figsize' : (12,9), + 'figure.subplot.left' : 0.07, # the left side of the subplots of the figure + 'figure.subplot.right' : 0.98 , # the right side of the subplots of the figure + 'figure.subplot.bottom' : 0.09 , # the bottom of the subplots of the figure + 'figure.subplot.top' : 0.9 , # the top of the subplots of the figure + 'figure.subplot.wspace' : 0.26 , # the amount of width reserved for blank space between subplots + 'figure.subplot.hspace' : 0.22 , # the amount of height reserved for white space between subplots + + 'lines.markersize' : 6, + 'lines.linewidth' : 2, + +# 'axes.formatter.limits' : (-2, 1), + + 'text.latex.unicode': True + } +rcParams.update(params) +rc('font', family='serif') +import sys +import os + +print "Plotting..." + +# Read Quickshed accelerations +data=loadtxt("particle_dump.dat") +id = data[:,0] +accx_e=data[:,5] +accy_e=data[:,6] +accz_e=data[:,7] + +accx_bh=data[:,8] +accy_bh=data[:,9] +accz_bh=data[:,10] + +accx_new=data[:,11] +accy_new=data[:,12] +accz_new=data[:,13] + +# Sort accelerations +rank = argsort(id) +id = id[rank] +accx_e = accx_e[rank] +accy_e = accy_e[rank] +accz_e = accz_e[rank] + +accx_bh = accx_bh[rank] +accy_bh = accy_bh[rank] +accz_bh = accz_bh[rank] + +accx_new = accx_new[rank] +accy_new = accy_new[rank] +accz_new = accz_new[rank] + +# Read Gadget accelerations +data=loadtxt("gadget_dump.dat") +id = data[:,0] +accx_g=data[:,5] +accy_g=data[:,6] +accz_g=data[:,7] + +# Sort accelerations +rank = argsort(id) +id = id[rank] +accx_g = accx_g[rank] +accy_g = accy_g[rank] +accz_g = accz_g[rank] + + +# Build error ------------------------------------------------ + +errx_bh = (accx_bh - accx_e )/abs(accx_e) +erry_bh = (accy_bh - accy_e )/abs(accy_e) +errz_bh = (accz_bh - accz_e )/abs(accz_e) + +errx_new = (accx_new - accx_e )/abs(accx_e) +erry_new = (accy_new - accy_e )/abs(accy_e) +errz_new = (accz_new - accz_e )/abs(accz_e) + +errx_g = (accx_g - accx_e )/abs(accx_e) +erry_g = (accy_g - accy_e )/abs(accy_e) +errz_g = (accz_g - accz_e )/abs(accz_e) + +# Statistics +meanx_bh = mean(errx_bh[abs(errx_bh) < 0.1]) +stdx_bh = std(errx_bh[abs(errx_bh) < 0.1]) +meany_bh = mean(erry_bh[abs(erry_bh) < 0.1]) +stdy_bh = std(erry_bh[abs(erry_bh) < 0.1]) +meanz_bh = mean(errz_bh[abs(errz_bh) < 0.1]) +stdz_bh = std(errz_bh[abs(errz_bh) < 0.1]) + +meanx_new = mean(errx_new[abs(errx_new) < 0.1]) +stdx_new = std(errx_new[abs(errx_new) < 0.1]) +meany_new = mean(erry_new[abs(erry_new) < 0.1]) +stdy_new = std(erry_new[abs(erry_new) < 0.1]) +meanz_new = mean(errz_new[abs(errz_new) < 0.1]) +stdz_new = std(errz_new[abs(errz_new) < 0.1]) + +# meanx_g = mean(errx_g) +# stdx_g = std(errx_g) +# meany_g = mean(erry_g) +# stdy_g = std(erry_g) +# meanz_g = mean(errz_g) +# stdz_g = std(errz_g) + +# Plot ------------------------------------------------------- +figure(frameon=True) + +subplot(311, title="Acceleration along X") +plot(id, errx_g , 'gs') +plot(id, errx_bh , 'rx') +plot(id, errx_new , 'b.') +text(id[-1], 0.18, "B-H: $%5.3f\\pm%5.3f$\n QuickShed: $%5.3f\\pm%5.3f$"%(meanx_bh, stdx_bh, meanx_new, stdx_new), backgroundcolor="w", va="top", ha="right" ) + + +ylim(-0.2, 0.2) +xlim(0,id[-1]) +grid() + +subplot(312, title="Acceleration along Y") +plot(id, erry_g , 'gs') +plot(id, erry_bh , 'rx') +plot(id, erry_new , 'b.') +text(id[-1], 0.18, "B-H: $%5.3f\\pm%5.3f$\n QuickShed: $%5.3f\\pm%5.3f$"%(meany_bh, stdy_bh, meany_new, stdy_new), backgroundcolor="w", va="top", ha="right" ) + + +ylim(-0.2, 0.2) +xlim(0,id[-1]) + +grid() + +subplot(313, title="Acceleration along Z") +plot(id, errz_g , 'gs', label="Gadget") +plot(id, errz_new , 'b.', label="QuickShed") +plot(id, errz_bh , 'rx', label="Legacy") +#text(id[-1], 0.18, "B-H: $%5.3f\\pm%5.3f$\n QuickShed: $%5.3f\\pm%5.3f$"%(meanz_bh, stdz_bh, meanz_new, stdz_new), backgroundcolor="w", va="top", ha="right" ) +legend(loc="upper right") + +ylim(-0.2, 0.2) +xlim(0,id[-1]) +grid() + +savefig("accelerations.png") + + + + +# Plot ------------------------------------------------------- +bins = linspace(-3, 3, 10000) + + +figure(frameon=True) +subplot(311, title="Acceleration along X")#, yscale='log') +hist(errx_g, bins=bins, normed=1, histtype='step', rwidth=0.01, color='g', label="Gadget") +hist(errx_bh, bins=bins, normed=1, histtype='step', rwidth=0.01, color='r', label="Legacy") +hist(errx_new, bins=bins, normed=1, histtype='step', rwidth=0.01, color='b', label="QuickShed") +legend(loc="upper right") +xlim(-0.03, 0.03) + +subplot(312, title="Acceleration along Y") +hist(erry_g, bins=bins, normed=1, histtype='step', rwidth=0.01, color='g') +hist(erry_bh, bins=bins, normed=1, histtype='step', rwidth=0.01, color='r') +hist(erry_new, bins=bins, normed=1, histtype='step', rwidth=0.01, color='b') +xlim(-0.03, 0.03) + +subplot(313, title="Acceleration along Z") +hist(errz_g, bins=bins, normed=1, histtype='step', rwidth=0.01, color='g') +hist(errz_bh, bins=bins, normed=1, histtype='step', rwidth=0.01, color='r') +hist(errz_new, bins=bins, normed=1, histtype='step', rwidth=0.01, color='b') +xlim(-0.03, 0.03) + +savefig("histogram.png") + + +# Print error +errx_g = errx_g[ (errx_g > -0.1) & (errx_g < 0.1) ] +erry_g = erry_g[ (erry_g > -0.1) & (erry_g < 0.1) ] +errz_g = errz_g[ (errz_g > -0.1) & (errz_g < 0.1) ] + +errx_bh = errx_bh[ (errx_bh > -0.1) & (errx_bh < 0.1) ] +erry_bh = erry_bh[ (erry_bh > -0.1) & (erry_bh < 0.1) ] +errz_bh = errz_bh[ (errz_bh > -0.1) & (errz_bh < 0.1) ] + +errx_new = errx_new[ (errx_new > -0.1) & (errx_new < 0.1) ] +erry_new = erry_new[ (erry_new > -0.1) & (erry_new < 0.1) ] +errz_new = errz_new[ (errz_new > -0.1) & (errz_new < 0.1) ] + + +print "Error for Gadget: ( x= %4.3e"%std(errx_g), "y= %4.3e"%std(erry_g), "z= %4.3e"%std(errz_g), ") Combined= %4.3e"%(( std(errx_g) + std(erry_g) + std(errz_g) )/3.) +print "Error for Legacy: ( x= %4.3e"%std(errx_bh), "y= %4.3e"%std(erry_bh), "z= %4.3e"%std(errz_bh), ") Combined= %4.3e"%(( std(errx_bh) + std(erry_bh) + std(errz_bh) )/3.) +print "Error for QuShed: ( x= %4.3e"%std(errx_new), "y= %4.3e"%std(erry_new), "z= %4.3e"%std(errz_new), ") Combined= %4.3e"%(( std(errx_new) + std(erry_new) + std(errz_new) )/3.)