Permalink
Browse files

Added a plotting script. Modded some of the shell script cheat sheets.

  • Loading branch information...
1 parent 54f064c commit c8af0a5c132c148a91f7a62d5e93c2c70015b4b7 @matthewbellis matthewbellis committed Apr 6, 2012
@@ -1,3 +1,7 @@
time ../../CUDA_code/Calculate_arc_length_two_datasets 2MASS_40k_RA_Dec.dat 2MASS_40k_RA_Dec.dat
time ../../CUDA_code/Calculate_arc_length_two_datasets 2MASS_40k_RA_Dec.dat 2MASS_40k_RA_Dec.dat -o logbinning_test.dat -w 0.2 -L 0.001 -N15 -l 2
time ../../CUDA_code/Calculate_arc_length_two_datasets 2MASS_40k_RA_Dec.dat 2MASS_40k_RA_Dec.dat -o logbinning_test.dat -w 0.5 -L 0.001 -N15 -l 1
+
+time ../../CUDA_code/Calculate_arc_length_two_datasets ~/Work/Astronomy/catalogs/Wechsler/wg100k_ra_dec.dat ~/Work/Astronomy/catalogs/Wechsler/wg100k_ra_dec.dat -o logbinning_test.dat -w 0.5 -L 0.001 -N25 -l 1
+time ../../CUDA_code/Calculate_arc_length_two_datasets ../2MASS_study/mc100k_ra_dec.dat ../2MASS_study/mc100k_ra_dec.dat -o logbinning_test_mcmc.dat -w 0.5 -L 0.001 -N25 -l 1
+time ../../CUDA_code/Calculate_arc_length_two_datasets ~/Work/Astronomy/catalogs/Wechsler/wg100k_ra_dec.dat ../2MASS_study/mc100k_ra_dec.dat -o logbinning_test_mc_data.dat -w 0.5 -L 0.001 -N25 -l 1
@@ -1 +1,2 @@
python generate_random_RA_Decl_in_radians.py --mass-range 6.28 --pos-range 1 -n 40000 > MonteCarlo_40k_RA_Dec.csv
+python generate_random_RA_Decl_in_radians.py --mass-range 6.28 --pos-range 1 -n 200000 > mc200k_ra_dec.dat
@@ -93,7 +93,7 @@ def main():
nparticles = int(options.nparticles)
############################################################################
############################################################################
- print "Right_Ascension , Declination"
+ #print "Right_Ascension , Declination"
# Generate the config file
print nparticles
@@ -104,34 +104,32 @@ def main():
# Generate the mass either from a flat or Gaussian distribution.
mass = [0.0]
if mass_flat:
- #value = 2.0*mass_range*np.random.random_sample() - mass_range
mass[0] = mass_range*np.random.random_sample()
- #mass[0] = 2.0*acos(value)
- #print value
- output += "%0.4f , " % (mass[0])
+ output += "%-7.4f " % (mass[0])
else:
# Make sure value is positive
while mass[0]<=0.0:
mass = np.random.normal(mass_mean, mass_width,1)
- output += "%0.4f , " % (mass[0])
+ output += "%-7.4f " % (mass[0])
# Generate the initial position either from a flat or Gaussian distribution.
pos = [0.0]
if pos_flat:
- #pos[0] = 2.0*pos_range*np.random.random_sample() - pos_range
pos[0] = acos(2.0*pos_range*np.random.random_sample() - pos_range) - 1.57
- output += "%0.4f , " % (pos[0])
+ output += "%-7.4f " % (pos[0])
else:
# Make sure value is positive
while pos[0]<=0.0:
pos = np.random.normal(pos_mean, pos_width,1)
- output += "%0.4f , " % (pos[0])
+ output += "%-7.4f " % (pos[0])
############### Do some acceptance correction by cutting out some stars
############### in a swath.
value = True
+ '''
if abs(pos[0] - cos(mass[0])) < 0.15:
value = False
+ '''
if value:
print output
@@ -161,7 +161,7 @@ def main():
subplots[0].set_xlabel(r"$\theta$ (degrees)", fontsize=24, weight='bold')
subplots[0].set_ylabel(r"w($\theta$)", fontsize=24, weight='bold')
- #subplots[0].set_xscale('log')
+ subplots[0].set_xscale('log')
#subplots[0].set_yscale('log')
#subplots[0].set_xlim(0.01,100)
@@ -0,0 +1,128 @@
+#!/usr/bin/env python
+"""
+Make a plot of some numbers read in from a .csv file with a header.
+"""
+
+import sys
+from math import *
+import numpy as np
+import csv
+import matplotlib.pyplot as plt
+import matplotlib.mlab as mlab
+import matplotlib.font_manager as fm
+
+from optparse import OptionParser
+from pylab import *
+
+
+################################################################################
+# main
+################################################################################
+def main():
+
+ # Parse the command line options
+ myusage = "\nusage: %prog [options] <file1.csv>"
+
+ parser = OptionParser(usage = myusage)
+ parser.add_option("-x", "--x-axis", dest="x_index", default=0,
+ help="Column in file to use as x-axis values.")
+ parser.add_option("-y", "--y-axis", dest="y_index", default=1,
+ help="Column in file to use as y-axis values.")
+ parser.add_option("-n", "--plot_every_nth_point", dest="plot_every_nth_point",
+ default=1,
+ help="If you're plotting a really big file, then only plot every nth point. [default=1]")
+
+ # Parse the options
+ (options, args) = parser.parse_args()
+
+ plot_every_nth_point = int(options.plot_every_nth_point)
+
+ ################################################################################
+ ################################################################################
+ # Make a figure on which to plot stuff.
+
+ figs = []
+ for i in xrange(2):
+ figs.append(plt.figure(figsize=(12, 8), dpi=90, facecolor='w', edgecolor='k'))
+# fig2 = plt.figure(figsize=(12, 8), dpi=90, facecolor='w', edgecolor='k')
+ #
+ # Usage is XYZ: X=how many rows to divide.
+ # Y=how many columns to divide.
+ # Z=which plot to plot based on the first being '1'.
+ # So '111' is just one plot on the main figure.
+ ################################################################################
+ subplots = []
+ for i in range(0,2):
+ subplots.append(figs[i].add_subplot(1,1,1))
+
+ ################################################################################
+ ################################################################################
+
+ ############################################################################
+ # Open the file, assuming that it is .csv format and checking to see
+ # that it exists.
+ ############################################################################
+ filename = args[0]
+
+ infile = open(filename,'r')
+
+ ############################################################################
+
+# formatter = ScalarFormatter()
+# formatter.set_scientific(True)
+# formatter.set_powerlimits((-4,4))
+
+ plt.xticks(fontsize=14, weight='bold')
+ plt.yticks(fontsize=14, weight='bold')
+
+ x_index = int(options.x_index)
+ y_index = int(options.y_index)
+
+
+ xpts = []
+ ypts = []
+ zpts = []
+ xaxis_title = ""
+ yaxis_title = ""
+
+
+ i=0
+ for row in infile:
+ #print row
+ if i%10000==0:
+ print i
+
+ if i>=20000000:
+ break
+
+ vals = row.split()
+ if i>0:
+ if len(vals)>=3:
+ if i%plot_every_nth_point==0:
+ xpts.append(float(vals[0]))
+ ypts.append(float(vals[1]))
+ zpts.append(float(vals[2]))
+
+ i += 1
+
+ #print xpts[i]
+ #print ypts[i]
+ #myplots = subplots[0].plot(xpts, ypts,'o',markersize=1)
+ h = subplots[1].hist(zpts,bins=500)
+
+ subplots[0].set_xlabel(xaxis_title, fontsize=14, weight='bold')
+ subplots[0].set_ylabel(yaxis_title, fontsize=14, weight='bold')
+
+ subplots[1].set_xlabel("z", fontsize=24, weight='bold')
+ subplots[1].set_ylabel("# galaxies", fontsize=24, weight='bold')
+
+ #infile_basename = filename[0].split('/')[-1].split('.')[0]
+ #output_file_name = "plot_together_%s_x%d_y%d.png" % (infile_basename,x_index,y_index)
+ #plt.savefig(output_file_name)
+ plt.show()
+
+################################################################################
+# Top-level script evironment
+################################################################################
+if __name__ == "__main__":
+ main()

0 comments on commit c8af0a5

Please sign in to comment.