https://github.com/miguelzuma/hi_class_public
Tip revision: 16ae0f6ccfcee513146ec36b690678f34fb687f4 authored by Miguel Zumalacarregui on 23 May 2020, 00:31:57 UTC
Merge pull request #8 from emiliobellini/hi_class
Merge pull request #8 from emiliobellini/hi_class
Tip revision: 16ae0f6
CPU
###################################################################
#
# CPU, a CLASS Plotting Utility
# v1.3
# Benjamin Audren, 07.11.2011
#
# This is a small python program aimed to gain time when comparing two
# spectra, i.e. from CAMB and CLASS, or a non-linear spectrum to a
# linear one. It is designed to be used in a command line fashion, not
# being restricted to your CLASS directory, though it recognized mainly
# CLASS output format. Far from perfect, or complete, it could use any
# suggestion for enhancing it, just to avoid losing time on useless
# matters for others. Be warned that, when comparing with other
# format, the following is assumed: there are no empty line (especially
# at the end of file). Gnuplot comment lines (starting with a # are
# allowed). This issue will cause a non-very descriptive error in CPU,
# any suggestion for testing it is welcome. Example of use: To
# superimpose two different spectra and see their global shape : python
# CPU output/lcdm_z2_pk.dat output/lncdm_z2_pk.dat To precisely compare
# the non-linear contribution of one-loop w.r.t. to TRG method: python
# CPU -i output/lcdm_1l_z1_pk.dat output/lcdm_1l_z1_pk_nl_density.dat
# output/lcdm_trg_pk_nl_density.dat The documentation available with
# --help should cover any question.
###################################################################
from scipy import *
from scipy import interpolate
import numpy as np
import os,sys,string,io,subprocess
import argparse
# parser for input arguments
parser = argparse.ArgumentParser(description='CPU, a CLASS Plotting Utility, specify wether you want to superimpose, divide or interpolate different files, and behold!',
epilog='''A standard usage would be, for instance :
python CPU output/test_1l_pk.dat output/test_1l_pk_nl_density.dat -i 0
python CPU output/wmap_cl.dat output/planck_cl.dat''',
formatter_class=argparse.RawDescriptionHelpFormatter)
parser.add_argument('files', metavar='Files',
type=str, nargs='*',
help='the relative path of the desired file to plot\nfrom the root directory of CLASS')
parser.add_argument('-d, --divide',
dest='merging',action='store_const',const='blend_against',default='blend_together',
help='divide the spectra from different files. The k-values must be strictly identical (default: plot every graph on the same plot)')
parser.add_argument('-i, --interpolate',metavar='index',
nargs='?',dest='interp',const=True,default=False,
help='interpolate the spectra on each set of k-values if different. Use to compare two (or more) graphs with different k-values. If nothing more is specified, it will plot both index in gnuplot. Otherwise, just add -i 0 to plot only the first one.')
parser.add_argument('-t', metavar='pk, cl_lin,etc...',
help='specify the file type (whether pk or cl (if not specified, cl==cl_lin, other choices are cl_log and cl_ll for log linear), only needed if your file name does not already contain one of these keywords...')
parser.add_argument('-colnum',metavar='2, 3, ...',
help='specify the column you want to plot. By default, column is 2: you are plotting TT spectrum')
parser.add_argument('-x, --x11',metavar='gnuplot terminal',
dest='term',action='store_const',const='x11',default='aqua',
help='if your version of gnuplot does not support the aqua terminal, it will pick instead the x11 one, for crappier displays')
parser.add_argument('-p, --print',
dest='printfile',action='store_true',default=False,
help='print the graph directly in a .png file')
parser.add_argument('-r, --repeat',
dest='repeat',action='store_true',default=False,
help='repeat the step for all redshifts with same base name, whatever the desired action')
parser.add_argument('-c, --cleaning',metavar='path',
nargs='?',dest='cleaning',const=True,default=False,
help='remove all .plt files in the current directory (if none specified) : keep your folders clean !')
# Remove all .plt, interp and divided files generated in output/
# By default, will clean current directory.
def clean(path):
pattern1,pattern2,pattern3='.plt','_interp.dat','_divided.dat'
path+="/"
print 'Cleaning {0} directory from .plt files ...'.format(path)
i=0
for each in os.listdir(path):
if ( (each.rfind(pattern1)!=-1) or (each.rfind(pattern2)!=-1) or (each.rfind(pattern3)!=-1) ):
name= "{0}{1}".format(path,each)
print ' deleting {0}'.format(name)
if os.remove(name)!= None:
break
i+=1
if i==0:
print ' ==> Already as clean as possible'
else:
print ' ==> Done'
# Errors ###########################
def error_format():
print " Hum... have you really provided a .dat file ?"
print " We apologise for the inconvenience, but CPU is exiting now"
exit()
def error_type():
print " Spectrum type unrecognized or unsupported yet, sorry!"
print " Please try the -t command in addition to your previous call"
print " We apologise for the inconvenience"
exit()
def error_number_of_files():
print " You specified a wrong number of files for the operation you demanded, maybe you do not want to divide only one file, for instance ?"
print " We apologise for the inconvenience"
exit()
####################################
# Subfunction called to write on the gnuplot script file the proper file names.
def printstring(names):
print_file=names[0].rstrip(".dat")
print_file+=".ps"
print_string="set terminal po enhanced color\nset output '{0}'\nreplot".format(print_file)
print '{0} has been generated'.format(print_file)
return print_string
# Subfunction called to write the most of the gnuplot script file according to options.
def headers_plot_file(spectrum_type,names,plot_line,term):
tmp=open(names[0],"r")
# printing option
if ((args.printfile is True) and (plot_line is not False)):
print_string=printstring(names)
else:
print_string=""
# for interpolation: index choice option
if plot_line is True:
if args.interp is not True:
if args.interp is not False:
index_string=" index {0} ".format(args.interp)
else:
index_string=" "
else:
index_string=" "
# if args colnum not specified, put it to 2,
if args.colnum is None:
args.colnum=2
if spectrum_type=='cl_lin':
if plot_line is True:
plot_string="plot "
for name in names:
plot_string+="'{0}'".format(name)+"{0}u 1:{1} w l,".format(index_string,args.colnum)
plot_string=plot_string.rstrip(",")
plot_string+="\n"
else:
plot_string=''
return "set terminal {0} enhanced\n".format(term)+"set xlabel 'l'\nset ylabel 'l(l+1)C_l / 2{/Symbol p}'\nset xr [2:]\nset format y '%.0t*10^{%T}'\nset key right\nset title 'CLASS output'\n"+plot_string+print_string
elif spectrum_type=='cl_log':
for line in tmp:
if line.rfind('multipoles')!=-1:
lmax= line.split(None)[-1]
break
if plot_line is True:
plot_string="plot "
for name in names:
plot_string+="'{0}'".format(name)+"{0}u 1:{1} w l,".format(index_string,args.colnum)
plot_string=plot_string.rstrip(",")
plot_string+="\n"
else:
plot_string=''
return "set terminal {0} enhanced\n".format(term)+"set logscale x\nset xr [2:{0}]\n".format(lmax)+"set xlabel 'l'\nset format y '%.0t*10^{%T}'\n"+"set ylabel 'l(l+1)C_l / 2{/Symbol p}'\nset key left\nset title 'CLASS output'\n"+plot_string+print_string
elif spectrum_type=='cl_ll':
for line in tmp:
if line.rfind('multipoles')!=-1:
lmax= line.split()[-1]
break
if plot_line is True:
plot_string="plot "
for name in names:
plot_string+="'{0}'".format(name)+"{0}u (sqrt($1)):{1} w l,".format(index_string,args.colnum)
plot_string=plot_string.rstrip(",")
plot_string=plot_string+"\n"
else:
plot_string=''
return "set terminal {0} enhanced\n".format(term)+"set xlabel 'l'\nset ylabel 'l(l+1)C_l / 2{/Symbol p}'\nset key right\nset title 'CLASS output'\nset xtics ('2' (2),'100' (100), '500' (500), '1000' (1000), '1500' (1500), '2000' (2000),'2500' (2500))"+"\nset format y '%.0t*10^{%T}'\n"+"set xr [2:{0}]\n".format(lmax)+plot_string+print_string
elif spectrum_type=='pk':
if plot_line is True:
plot_string="plot "
for name in names:
plot_string+="'{0}'".format(name)+"{0}w l,".format(index_string)
plot_string=plot_string.rstrip(",")
plot_string+="\n"
else:
plot_string=""
for line in tmp:
if line.rfind('redshift')!=-1:
z= line.split("=")[1]
z= z.rstrip("\n")
break
else:
z= 'I have no idea'
return "set terminal {0} enhanced\n".format(term)+"set logscale\nset xlabel 'k (h/Mpc)'\nset ylabel 'P_k (Mpc/h)^3'\nset key right\nset title 'Power spectrum at z={0}'\n".format(z)+plot_string+print_string
else:
return None
tmp.close()
# Print all asked files one next to the other (default operation if files with different k values)
def blend_together(names,spectrum_type,term):
gnuplot_file=names[0].replace(".dat",".plt")
if gnuplot_file==names[0]:
error_format()
print ' creating {0}'.format(gnuplot_file)
plotfile = open(gnuplot_file, "w")
plotfile.write(headers_plot_file(spectrum_type, names,True,term))
plotfile.close()
_plot(gnuplot_file)
# Print all asked files with respect to the same k (or anything) values, all y-data being divided by the y data of the first file.
# Only valid if the k values are exactly the same (values and number of values).
def blend_against(names,spectrum_type,term):
gnuplot_file=names[0].replace(".dat",".plt")
data_file=names[0].replace(".dat","_divided.dat")
if gnuplot_file==names[0]:
error_format()
print ' creating {0} and {1}'.format(gnuplot_file,data_file)
string=['' for rows in range(10000)]
imax=0
datafile=open(data_file,"w")
for name in names:
i=0
tmp=open(name,"r")
for line in tmp:
if ((line.find('#')==-1) and (line.rfind('000000')==-1)):
string[i]=string[i]+line.rstrip("\n")+"\t"
if i>imax:
imax=i
i+=1
tmp.close()
for i in range(imax):
datafile.write(string[i]+"\n")
plotfile = open(gnuplot_file,"w")
plotfile.write(headers_plot_file(spectrum_type,names,False,term))
plot_string='unset logscale\nplot '
x_base=1
field_base=2
field=2
size=len(names)
for i in range(size):
field+=2
plot_string+="'{0}' u {1}:(${2}/${3}) w l,".format(data_file,x_base,field,field_base)
plot_string=plot_string.rstrip(',')
plot_string+="\n"
plotfile.write(plot_string)
plotfile.close()
datafile.close()
_plot(gnuplot_file)
# If k values are different, an interpolation is done and outputs a data file. For a two files case:
# File 1 contains ( k1 | P1(k1) ), File 2 ( k2 | P2(k2) ). The data file created will contain:
# k1 | P1(k1) | P2(k1)_interp
#(blank space for gnuplot using index 0, etc)
# k2 | P1(k2)_interp | P2(k2)
def blend_against_interp(names,spectrum_type,term):
gnuplot_file=names[0].replace(".dat",".plt")
data_file=names[0].replace(".dat","_interp.dat")
if gnuplot_file==names[0]:
error_format()
print ' creating {0} and {1}'.format(gnuplot_file,data_file)
plotfile = open(gnuplot_file,"w")
l=0
jmax=[0 for col in range(10)]
lmax=0
kmax=10000000
kmin=0
spam = np.array([[[0 for col in range(10)] for row in range(10000)] for depth in range(2)],dtype=float)
for name in names:
currentfile = open(name,"r")
print ' reading {0}..'.format(name)
j=0
for line in currentfile:
if (line.find('#')==-1):
line=line.split()
spam[0,j,l]=float(line[0])
spam[1,j,l]=float(line[1])
j+=1
jmax[l]=j
if spam[0,jmax[l]-1,l]<=kmax:
kmax=spam[0,jmax[l]-1,l]
if spam[0,0,l]>=kmin:
kmin=spam[0,0,l]
l+=1
currentfile.close()
lmax=l
lower_bound=[0 for col in range(10)]
upper_bound=[0 for col in range(10)]
#determining the upper and lower bound for each file, to only do interpolation
for l in range (lmax):
for j in range (jmax[l]):
if ((spam[0,j,l]<kmin) and (spam[0,j+1,l]>=kmin)):
lower_bound[l]=j+1
if (spam[0,j,l]<=kmax):
upper_bound[l]=j
print ' -> done'
#creating the new data file
curves=[np.array for row in range(lmax)]
interpolated=[np.array for row in range(lmax)]
for l in range (lmax):
x=spam[0,lower_bound[l]:upper_bound[l],l]
y=spam[1,lower_bound[l]:upper_bound[l],l]
curves[l]=interpolate.splrep(x,y)
datafile = open(data_file,"w")
for l in range (lmax):
for ll in range(lmax):
if ll==l:
interpolated[ll]=spam[1,lower_bound[l]:upper_bound[l],l]
else:
x2=spam[0,lower_bound[l]:upper_bound[l],l]
interpolated[ll]=interpolate.splev(x2,curves[ll],der=0)
for j in range(upper_bound[l]-lower_bound[l]-2):
data_string=''
for ll in range(lmax):
data_string=data_string+str((interpolated[ll])[j+1])+"\t"
datafile.write(str(spam[0,j+lower_bound[l]+1,l])+"\t"+data_string+"\n")
datafile.write("\n\n")
datafile.close()
# if index chosen
if args.interp is not True:
index_string_interp=" index {0}".format(args.interp)
else:
index_string_interp=""
plotfile = open(gnuplot_file,"w")
plotfile.write(headers_plot_file(spectrum_type,names,False,term))
plotfile.write("unset logscale \n")
if args.t in ['cl_log', 'pk']:
plotfile.write("set logscale x\n")
plotfile.write("set xr [{0}:{1}]\n".format(kmin,kmax))
#plotfile.write("set yr [-0.001:0.005]\n")
plot_string="plot "
for l in range(lmax-1):
plot_string+="'{0}' {1} u 1:(${2}/$2-1) w l,".format(data_file,index_string_interp,l+3)
plot_string=plot_string.rstrip(",")
plot_string=plot_string+"\n"
plotfile.write(plot_string)
if (args.printfile is True):
plotfile.write(printstring(names))
plotfile.close()
_plot(gnuplot_file)
# Launch session of gnuplot with generated gnuplot script file
def _plot(gnuplot_file):
os.system("gnuplot -persist '{0}'".format(gnuplot_file))
#######################################################
################## MAIN PART ##########################
#######################################################
print '~~~ CPU, a CLASS Plotting Utility ~~~'
args = parser.parse_args()
# check if the user want to clean its directory first
if args.cleaning is not False:
if args.cleaning is not True:
clean(args.cleaning)
else:
clean(os.getcwd())
exit()
# if there are no argument in the input, print usage
if len(args.files)==0:
parser.print_usage()
exit()
# if the first file name contains cl or pk, infer the type of desired spectrum
if ((args.files[0].rfind('cl')!=-1) and (args.t is None)):
spectrum_type='cl_lin'
args.t = 'cl_lin'
elif args.files[0].rfind('pk')!=-1:
spectrum_type='pk'
args.t = 'pk'
elif args.t is not None:
spectrum_type=args.t
else:
error_type()
# repeater
temp_path=args.files[0].split("/")
if len(temp_path)> 1:
path=temp_path[-2]
else:
path=os.getcwd()
list_file=['' for i in range(len(args.files)*10)]
if (temp_path[-1].rfind("z")!=-1 and args.repeat is True):
root_name=temp_path[-1].split("z")[-2]
for any in range (0,len(args.files)):
extension_name=args.files[any].split("/")[-1].split("_")[-1]
i=0
for each in os.listdir(path):
if (("z" in each) and (each.split("z")[-2]==root_name) and (each.split("_")[-1]==extension_name)):
if len(temp_path)>1:
list_file[any+len(args.files)*i]=temp_path[-2]+"/"+each
else:
list_file[any+len(args.files)*i]=each
i+=1
if args.repeat is True:
repeat_len=i-1
else:
repeat_len=0
else:
for any in range(0,len(args.files)):
list_file[any]=args.files[any]
repeat_len=0
# actual computation
for i in range(0,repeat_len+1):
local_files=['' for k in range(len(args.files))]
for j in range(0,len(args.files)):
local_files[j]=list_file[j+len(args.files)*i]
if args.interp is False:
if args.merging=='blend_together':
blend_together(local_files,spectrum_type,args.term)
else:
if len(local_files)<2:
error_number_of_files()
else:
blend_against(local_files,spectrum_type,args.term)
else:
print '**interpolating (please wait)'
blend_against_interp(local_files,spectrum_type,args.term)
exit()