API tools faq
paste
Advertisement
dailerfm

gsTools

dailerfm
Feb 10th, 2020
139
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 11.44 KB | None | 0 0
raw download clone embed print report
  1. import numpy as np
  2. import matplotlib.pyplot as plt
  3. import os
  4. import sys
  5. import corner
  6.  
  7. paramname = {r'rho0':r'$\log_{10} \rho_0$','gamma0':r'$\gamma_0$','gamma1':r'$\gamma_1$','gamma2':r'$\gamma_2$','gamma3':r'$\gamma_3$','gamma4':r'$\gamma_4$',r'rhos':r'$\log_{10} \rho_s$', r'rs':r'$\log_{10}r_s$', r'alpha':r'$\alpha$', r'beta':r'$\beta$', r'gamma':r'$\gamma$', r'beta0':r'$\beta_0$', r'betainf':r'$\beta_{\infty}$', r'ra':r'$\log_{10}r_a$', r'eta':r'$\eta$', r'm1':r'$\log_{10}m_1$', r'a1':r'$a_1$', r'm2':r'$\log_{10}m_2$',r'a2':r'$a_2$', r'm3':r'$\log_{10}m_3$', r'a3':r'$a_3$', r'mstar':r'$\log_{10}M_*$' }
  8.  
  9.  
  10. def check_galdata(all_gals, workdir):
  11.  
  12. for gal in all_gals:
  13. all_good = True
  14.  
  15. if os.path.exists(workdir + '/GalaxyData/' + gal + '_VSPs.txt') == False:
  16. print("No VSP file in GalaxyData for ", gal)
  17. all_good = False
  18. if os.path.exists(workdir + '/GalaxyData/' + gal + '_Mstar.txt') == False:
  19. print("No Mstar file in GalaxyData for ",gal)
  20. all_good = False
  21. if os.path.exists(workdir + '/GalaxyData/' + gal + '_PlumParam.txt') == False:
  22. print("No PlumParam file in GalaxyData for ",gal)
  23. all_good = False
  24. if os.path.exists(workdir + '/GalaxyData/' + gal + '_Rhalf.txt') == False:
  25. print("No Rhalf file in GalaxyData for ",gal)
  26. all_good = False
  27. if os.path.exists(workdir + '/GalaxyData/' + gal + '_SurfDen.txt') == False:
  28. print("No SurfDen file in GalaxyData for ",gal)
  29. all_good = False
  30. if os.path.exists(workdir + '/GalaxyData/' + gal + '_KinDat.txt') == False:
  31. print("No KinDat file in GalaxyData for ",gal)
  32. all_good = False
  33. if all_good == False:
  34. print("Please pre-process the data for galaxy ", gal)
  35. print("Goodbye!")
  36. sys.exit()
  37. return 0
  38.  
  39.  
  40. def check_float(text):
  41. isfloat = False
  42.  
  43. while isfloat == False:
  44. try:
  45. var = float(input(text))
  46. isfloat = True
  47. except ValueError:
  48. print("Must be a float/integer")
  49. #isfloat = False
  50. return float(var)
  51.  
  52. def checkdirs(workdir, codedir):
  53. kinphot = os.path.isdir(workdir + "/KinPhotDat")
  54. if kinphot == False:
  55. os.system("mkdir " + workdir + "/KinPhotDat")
  56. galdat = os.path.isdir(workdir + "/GalaxyData")
  57. if galdat == False:
  58. os.system("mkdir " + workdir + "/GalaxyData")
  59. gal_list = os.path.exists(workdir + "/galaxy_list.txt")
  60. if gal_list == False:
  61. os.system("touch " + workdir + "/galaxy_list.txt")
  62. projects = os.path.exists(workdir + "/projects.txt")
  63. if projects == False:
  64. os.system("touch " + workdir + "/projects.txt")
  65.  
  66. code = os.path.exists(workdir + "/pygravsphere.py")
  67. if code == False:
  68. os.system("cp " + codedir + "/pygravsphere.py " + workdir + "/pygravsphere.py")
  69. return 0
  70.  
  71. def get_cmap(n, name='hsv'):
  72.  
  73. return plt.cm.get_cmap(name, n)
  74.  
  75. def plot_chains(samples, workdir, project_name, galaxy, walkers_plot):
  76. foptions = open(workdir + project_name + '/options.txt', 'r')
  77. input_opt = (foptions.readline()).split()
  78.  
  79. steps = int(input_opt[3])
  80. burn_in = int(input_opt[4])
  81. nwalkers = int(input_opt[5])
  82.  
  83. ll = len(samples)/(100*nwalkers)
  84. sample_split = np.array_split(samples, ll)
  85. sample_new = [[] for i in range(nwalkers)]
  86.  
  87.  
  88. params = np.loadtxt(workdir + project_name + '/Submissions/priors.txt', dtype = 'str')
  89. valid, = np.where(params[:,6] == 'False')
  90. param_names = params[valid, 0]
  91.  
  92. for s in range(0, ll):
  93. second_split = np.array_split(sample_split[s], nwalkers)
  94. for s2 in range(nwalkers):
  95. sample_new[s2].append(second_split[s2])
  96.  
  97. sample_new = np.concatenate(np.concatenate(sample_new))
  98.  
  99. mod = int(float(nwalkers)/float(walkers_plot))
  100.  
  101. for p in range(0, len(param_names)):
  102.  
  103. fig, ax1 = plt.subplots(1,1, figsize = (16,4))
  104.  
  105. cmap = get_cmap(nwalkers)
  106.  
  107. plot_sample = np.array_split(sample_new, nwalkers)
  108.  
  109. for i in range(nwalkers):
  110. if i%mod == 0:
  111. ax1.plot(plot_sample[i][:,p], color =cmap(i), lw = 1.0, alpha = 0.5)
  112.  
  113.  
  114. ax1.set_ylabel(paramname.get(param_names[p]), fontsize = 12)
  115. ax1.set_xlabel('Step', fontsize = 12)
  116.  
  117.  
  118. fig.savefig(workdir + project_name + '/%s/' % galaxy + 'Param_%d.png' % p, bbox_inches = "tight")
  119. plt.close()
  120. return 0
  121.  
  122.  
  123. def plot_triangle(samples, workdir, project_name, galaxy):
  124.  
  125. params = np.loadtxt(workdir + project_name + '/Submissions/priors.txt', dtype = 'str')
  126. valid, = np.where(params[:,6] == 'False')
  127. param_names = params[valid, 0]
  128.  
  129. #print np.vectorize(paramname.get)(param_names)
  130.  
  131. fig = corner.corner(samples[:,:-1], labels=np.vectorize(paramname.get)(param_names), quantiles=[0.16, 0.5, 0.84],title_kwargs={"fontsize": 14})
  132. fig.savefig(workdir + project_name + '/%s/' % galaxy + 'Triangle.png')
  133. plt.close()
  134. return 0
  135.  
  136.  
  137.  
  138.  
  139. def get_lims(data,r):
  140.  
  141.  
  142.  
  143.  
  144. out = np.zeros((len(r), 6))
  145.  
  146.  
  147. out[:,0] = r
  148. out[:,1] = np.median(data, axis = 0)
  149. out[:,2] = np.percentile(data, 16, axis = 0)
  150. out[:,3]= np.percentile(data, 84, axis = 0)
  151. out[:,4]= np.percentile(data, 2.5, axis = 0)
  152. out[:,5]= np.percentile(data, 97.5, axis = 0)
  153.  
  154. return out
  155.  
  156.  
  157. def get_lims_all(data,r):
  158.  
  159.  
  160.  
  161.  
  162. out = np.zeros((len(r), 102))
  163. out[:,0] = r
  164. for i in range(0, 101):
  165.  
  166. out[:,i+1] = np.percentile(data, i, axis = 0)
  167.  
  168.  
  169.  
  170.  
  171. return out
  172.  
  173.  
  174. def get_lims_loglog(data, tot_bins):
  175.  
  176.  
  177.  
  178. out = np.zeros((len(tot_bins), 6))
  179.  
  180. out[:,0] = tot_bins
  181. out[:,1] = np.median(data, axis = 0)
  182. out[:,2] = np.percentile(data, 16, axis = 0)
  183. out[:,3]= np.percentile(data, 84, axis = 0)
  184. out[:,4]= np.percentile(data, 2.5, axis = 0)
  185. out[:,5]= np.percentile(data, 97.5, axis = 0)
  186.  
  187. return out
  188.  
  189.  
  190. def get_lims_sig(data, workdir, galaxy):
  191. pos = np.loadtxt(workdir + '/GalaxyData/%s_KinDat'% (galaxy) + '.txt')
  192. pos = pos[:,0]
  193.  
  194. out = np.zeros((len(pos), 6))
  195.  
  196. out[:,0] = pos
  197. out[:,1] = np.median(data, axis = 0)
  198. out[:,2] = np.percentile(data, 16, axis = 0)
  199. out[:,3]= np.percentile(data, 84, axis = 0)
  200. out[:,4]= np.percentile(data, 2.5, axis = 0)
  201. out[:,5]= np.percentile(data, 97.5, axis = 0)
  202.  
  203. return out
  204.  
  205. def create_sub(project_name, num_cores, timevar, workdir,codedir, anis, darkmatter, vsps, plummer, num_walkers, burn_in, steps, int_points, mpi_opt):
  206. if mpi_opt == '3':
  207. core = int(num_cores)
  208. galaxies = np.loadtxt(workdir + '/galaxy_list.txt', ndmin = 1,dtype = 'str')
  209. prestr1 = workdir + project_name + '/Submissions/'
  210. prestr2 = workdir + project_name + '/OutErr/'
  211.  
  212. for galaxy in galaxies:
  213. with open(codedir + "/sub_script.txt") as forg:
  214. newText=forg.read()
  215.  
  216. newText = newText.replace('CORENUM', "%d" %core)
  217. newText = newText.replace('-J GALID', '-J %s_' % galaxy + '%s' % project_name)
  218. newText = newText.replace('GALID.err', prestr2 + "%s_" %galaxy + "%s.err" % project_name)
  219. newText = newText.replace('GALID.out', prestr2 + "%s_" %galaxy + "%s.out" % project_name )
  220. newText = newText.replace('TIME', "%s" %timevar )
  221.  
  222. forg.close()
  223.  
  224. with open(prestr1 + '%s.sh' %galaxy, "w") as f:
  225. f.write(newText + '\n')
  226. f.write("python " + codedir + '/write_script_mpi.py' + " %s"%(workdir) + " %s"%(codedir) + " %s" %(project_name) + " %s"%(galaxy) + " %s"%(core) + " %s"%(num_walkers) + " %s"%(burn_in) + " %s"%(steps) + " %s"%(int_points) + " %s"%(darkmatter) + " %s"%(anis) + " %s"%(vsps) + " %s"%(plummer) + " standard" )
  227. f.close()
  228. os.system("chmod u+x " + prestr1 + '%s.sh' % galaxy)
  229.  
  230. elif mpi_opt == '2':
  231. core = int(num_cores)
  232. galaxies = np.loadtxt(workdir + '/galaxy_list.txt', ndmin = 1,dtype = 'str')
  233. prestr1 = workdir + project_name + '/Submissions/'
  234. prestr2 = workdir + project_name + '/OutErr/'
  235. for galaxy in galaxies:
  236. with open(prestr1 + '%s.sh' %galaxy, "w") as f:
  237. f.write("python " + codedir + '/write_script_multi.py' + " %s"%(workdir) + " %s"%(codedir) + " %s" %(project_name) + " %s"%(galaxy) + " %s"%(core) + " %s"%(num_walkers) + " %s"%(burn_in) + " %s"%(steps) + " %s"%(int_points) + " %s"%(darkmatter) + " %s"%(anis) + " %s"%(vsps) + " %s"%(plummer) + " standard" )
  238. f.close()
  239. os.system("chmod u+x " + prestr1 + '%s.sh' % galaxy)
  240.  
  241.  
  242.  
  243. else:
  244. galaxies = np.loadtxt(workdir + '/galaxy_list.txt', ndmin = 1,dtype = 'str')
  245. prestr1 = workdir + project_name + '/Submissions/'
  246. prestr2 = workdir + project_name + '/OutErr/'
  247. for galaxy in galaxies:
  248. with open(prestr1 + '%s.sh' %galaxy, "w") as f:
  249. f.write("#!/bin/sh" + '\n')
  250. f.write("python " + codedir + '/write_script.py' + " %s"%(workdir) + " %s"%(codedir) + " %s" %(project_name) + " %s"%(galaxy) + " %s"%(num_walkers) + " %s"%(burn_in) + " %s"%(steps) + " %s"%(int_points) + " %s"%(darkmatter) + " %s"%(anis) + " %s"%(vsps) + " %s"%(plummer) + " standard" )
  251. f.close()
  252. os.system("chmod u+x " + prestr1 + '%s.sh' % galaxy)
  253.  
  254. return 0
  255.  
  256. def create_ana_sub(project_name, workdir, codedir, dm_option, beta_option, plummer_option, samples, mpi_opt, cut_off, chi_cut, min_rad, max_rad, points,timevar):
  257.  
  258.  
  259. core = 1
  260. galaxies = np.loadtxt(workdir + '/galaxy_list.txt', ndmin = 1,dtype = 'str')
  261.  
  262.  
  263.  
  264. prestr1 = workdir + project_name + '/Analysis/Submissions/'
  265. prestr2 = workdir + project_name + '/Analysis/OutErr/'
  266.  
  267. if mpi_opt == 'y':
  268. for galaxy in galaxies:
  269. with open(codedir + "/sub_script.txt") as forg:
  270. newText=forg.read()
  271. newText = newText.replace('CORENUM', "%d" %core)
  272. newText = newText.replace('-J GALID', '-J %s_' % galaxy + '%s' % project_name)
  273. newText = newText.replace('GALID.err', prestr2 + "%s_" %galaxy + "%s.err" % project_name)
  274. newText = newText.replace('GALID.out', prestr2 + "%s_" %galaxy + "%s.out" % project_name )
  275. newText = newText.replace('TIME', "%s" %timevar )
  276.  
  277. forg.close()
  278. with open(prestr1 + '%s.sh' %galaxy, "w") as f:
  279. f.write(newText + '\n')
  280.  
  281. f.write("python " + codedir + '/AnalysisCode/' + 'gs_analysis.py ' + workdir +' ' +codedir + ' '+ project_name + ' ' + str(galaxy) + ' ' + str(int(samples)) + ' ' + str(int(cut_off)) + ' ' + str(float(chi_cut)) + ' ' + str(float(min_rad)) + ' ' + str(float(max_rad)) + ' ' + str(int(points)) + '\n')
  282.  
  283. f.close()
  284.  
  285. else:
  286. for galaxy in galaxies:
  287. with open(prestr1 + '%s.sh' %galaxy, "w") as f:
  288. f.write("python " + codedir + '/AnalysisCode/' + 'gs_analysis.py ' + workdir +' ' +codedir + ' '+ project_name + ' ' + str(galaxy) + ' ' + str(int(samples)) + ' ' + str(int(cut_off)) + ' ' + str(float(chi_cut)) + ' ' + str(float(min_rad)) + ' ' + str(float(max_rad)) + ' ' + str(int(points)) + '\n')
  289. f.close()
  290. os.system("chmod u+x " + prestr1 + '%s.sh' %galaxy)
  291.  
  292. return 0
  293.  
  294. def submit_jobs(workdir, project_name, all_gals, sub_com):
  295. sub_com = sub_com.strip()
  296. for galaxy in all_gals:
  297. if not sub_com:
  298. sub_script = open(workdir + '/' + project_name + '/Submissions/' + '%s.sh' %galaxy, 'r')
  299. first = sub_script.readline()
  300. first = first.split()
  301. if first[0] == 'python':
  302. os.system('. ' + workdir + '/' + project_name + '/Submissions/' + '%s.sh' % galaxy)
  303. else:
  304. sub_script = open(workdir + '/' + project_name + '/Submissions/' + '%s.sh' % galaxy, 'r')
  305. first = sub_script.readline()
  306. first = first.split()
  307. if first[0] == 'python':
  308. os.system('. ' + workdir + '/' + project_name + '/Submissions/' + '%s.sh' % galaxy)
  309. else:
  310. os.system(sub_com.strip() + ' ' + workdir + '/' + project_name + '/Submissions/' + '%s.sh' % galaxy)
  311.  
  312. return 0
  313.  
  314. def preprocess(workdir, codedir, all_gals):
  315. sys.path.append(codedir)
  316. from GSpro import gal_input
  317. for galaxy in all_gals:
  318. gal_input.galaxy_data(galaxy, workdir + '/GalaxyData/', workdir + '/KinPhotDat/')
  319.  
  320. return 0
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!