From 617b685bfb10eaea7aed80a57e241577eb7ac1c6 Mon Sep 17 00:00:00 2001 From: healther Date: Mon, 8 Apr 2013 14:41:49 +0200 Subject: [PATCH] docstrings added docstrings to analysis, distribution, execute, setup and starformation fixed a bug in analysis (wrong detection area due to magnetudes) fixed a bug in plot due to inf-values if no stars are detected --- analysis.py | 28 ++--- distribution.py | 21 +++- execute.py | 15 ++- plot.py | 37 ++++--- setup.py | 11 +- starformation.py | 43 +++++--- starformationtest.py | 256 +++++++++++++++++++++++++++++++++++++++++++ 7 files changed, 360 insertions(+), 51 deletions(-) create mode 100755 starformationtest.py diff --git a/analysis.py b/analysis.py index 71f2d5c..9977786 100755 --- a/analysis.py +++ b/analysis.py @@ -11,9 +11,19 @@ def main(folder, quiet=0): - '''This script gives the number of "observed" stars from the sampled datafiles in "folder"/ + '''analysis.main(folder, quiet=0) + +This script gives the number of "observed" stars from the sampled datafiles in "folder"/ according to the selection criteria from Yusef-Zadeh et al +input: +folder String Specifiecs the folder where the files are +quiet boolean =1 surpresses all standard output + +output: +returns a file named __expected_number in which contains a numpy array of the simulation parameters +number, A_v, Aperature_size, Age and the expected 'detected' number +The first line of the file is an ','-seperated head of the contained informations ''' if quiet: @@ -36,21 +46,12 @@ def main(folder, quiet=0): #only using files created by the automated simulation if fil.startswith('sim_') and not 'settings' in fil.encode("ascii"): print ("%s/%s" % (folder,fil.encode("ascii")), file=output_stream) - #f = open("%s/%s" % (folder,fil.encode("ascii")), 'r') - #headers = f.readline().strip().split(',') - #data = np.loadtxt(f) - #f.close() # Read in hdulist = fits.open('%s/%s' %(folder,fil)) -# av = hdulist[1].header['age'] data = hdulist[1].data - ##selecting the relevant columns - #c1 = headers.index('corrected_flux %s' % color1) - #c2 = headers.index('corrected_flux %s' % color2) - #calculating magnitudes from fluxes and converting to CMD-data x = -2.5*(np.log10(data['cflux %s' % color1]/64130) - np.log10(data['cflux %s' % color2]/7140)) y = -2.5*(np.log10(data['cflux %s' % color2]/7140)) @@ -59,9 +60,10 @@ def main(folder, quiet=0): # efficiency? accuracy? n=0 #selecting "observed" stars - for i in range(len(x)): - if (y[i] < -10./3. * (x[i]-1.) + 10.) and (max_mag < y[i] < min_mag): - n = n+1 + #for i in range(len(x)): + #if (y[i] > -10./3. * (x[i]-1.) + 10.) and (max_mag < y[i] < min_mag): + #n = n+1 + n = sum(np.logical_and( (y > -10./3. * (x-1.) + 10.), np.logical_and(max_mag < y, y < min_mag))) tmp, av, apera, age = fil.split('_') out.append([Decimal(av), Decimal(apera), Decimal(age), n]) diff --git a/distribution.py b/distribution.py index 1d74855..753bc04 100755 --- a/distribution.py +++ b/distribution.py @@ -9,8 +9,27 @@ class distribution( object ): - '''to be added + '''distribution(func, a, b) +A distribution object represents a probability distribution, the provided function is taken as the +probability density function and an cumulative density function and its inverse a calculated. a and +b are taken as the lower and the upper bound respectively. +func need not to be normalized + +input: +func function specifies the probability density function of the distribution +a float specifies the lower bound of the distribution +b float specifies the upper bound of the distribution + +members: +pdf() returns the probability density function, takes a float as an argument +cdf() returns the cumulative probability density function, takes a float as an argument +ppf() returns the inverse cumulative probability funciton, takes a float as an argument +norm() returns the normalisation factor (integral over pdf) +sample(int n) returns an array of n sampled values +mean() returns the mean value of the distribution +_lowerbound lowerbound of the distribution +_upperbound upperbound of the distribution ''' def __init__(self, func = lambda x: 1, a = 0., b = 1000.): self._pdf = func diff --git a/execute.py b/execute.py index e1320f0..be27d6f 100644 --- a/execute.py +++ b/execute.py @@ -11,7 +11,9 @@ def main(quiet = False): - ''' This script produces a grid of expected numbers of stars according to the selection + '''execute.main(quiet = False) + +This script produces a grid of expected numbers of stars according to the selection criteria of Yusef-Zedah et al. 2009, 702,178-225 The Astrophysical Journal. The grid is in av for visual extinction, apera for aperature size and age for the maxage of the starformation size @@ -25,7 +27,7 @@ def main(quiet = False): else: output_stream = sys.stdout - sfr = .08 + sfr = .01 # star mass function def f(x): # Kouper IMF #http://adsabs.harvard.edu/abs/2001MNRAS.322..231K @@ -61,8 +63,9 @@ def g(x): for i in range(len(avs)): for j in range(len(aperas)): for k in range(len(ages)): - starformation.main(massfunction = mf, starformationhistory = sf[k], A_v = avs[i], sfr = .01, \ - apera = aperas[j], maxage = ages[k], appendix = "%s_%03d_%06d_%07d" % ('sim',avs[i],aperas[j],ages[k]), quiet=True) + starformation.main(massfunction = mf, starformationhistory = sf[k], \ + A_v = avs[i], sfr = sfr, apera = aperas[j], maxage = ages[k], \ + appendix = "%s_%03d_%06d_%07d" % ('sim',avs[i],aperas[j],ages[k]), quiet=True) print(avs[i],aperas[j],ages[k], l/len(avs)/len(aperas)/len(ages), file=output_stream) l = l+1 parameters.append([avs[i],aperas[j],ages[k]]) @@ -70,7 +73,7 @@ def g(x): t2 = time() # end of simulation print(t2, t1, t2-t1) - print ('number of simulations run: %s' %k , file=output_stream) + print ('number of simulations run: %s' %l , file=output_stream) head = ['AV', 'Aperature_size', 'Age'] f = open('out/__head', 'w') f.write( ','.join(head)+'\n' ) @@ -79,7 +82,7 @@ def g(x): t3 = time() # end of saving data - analysis.main('out', quiet=True) + analysis.main('out') print ('analysis complete' , file=output_stream) t4 = time() # end of analysing data diff --git a/plot.py b/plot.py index 5a10ab8..2a27984 100644 --- a/plot.py +++ b/plot.py @@ -1,5 +1,5 @@ # -*- coding: utf-8 -*- -#test commentary +from __future__ import print_function, division from mpl_toolkits.mplot3d import Axes3D from matplotlib import cm from matplotlib.mlab import griddata @@ -11,7 +11,7 @@ def main(): aperature = 1 - visual_ex = 1 + visual_ex = 0 f = open('out/__expected_number', 'r') @@ -23,14 +23,13 @@ def main(): fig = plt.figure() ax = fig.gca(projection='3d') avs = np.linspace(10.0, 50.0, 5) - aperas = np.linspace(10000, 50000, 5) + aperas = np.logspace(2, 5, 4) ages = np.linspace(500000, 2000000, 11) - numbers = data[:,3].reshape(5, 5, 11) - numbers = np.roll(numbers,4,2) + numbers = np.clip(data[:,3].reshape(5, 4, 11),5.,100000.) #prevent inf if no stars are detected X, Y = np.meshgrid(avs, ages) - Z = 559.*.01/np.transpose(numbers[:,aperature,:]) + Z = 559.*.08/np.transpose(numbers[:,aperature,:]) surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm) fig.colorbar(surf) @@ -44,30 +43,34 @@ def main(): ax.set_xlabel('av') ax.set_ylabel('age') ax.set_zlabel('starformation rate in M_sun/year') + #ax.set_zlim(0.,1.) + ax.set_title('Starformation as a function of\n visual extinction Av and age at apperaturesize=1kpc') - plt.savefig('plot/av-age.svg') + #plt.savefig('plot/av-age.svg') fig2 = plt.figure() ax2 = fig2.gca(projection='3d') - X2, Y2 = np.meshgrid(ages, aperas) - Z2 = 559.*.01/numbers[visual_ex,:,:] + X2, Y2 = np.meshgrid(ages, np.log10(aperas)) + Z2 = 559.*.08/numbers[visual_ex,:,:] surf2 = ax2.plot_surface(X2, Y2, Z2, rstride=1, cstride=1, cmap=cm.coolwarm) fig2.colorbar(surf2) - x2 = data[:,2].reshape(5, 5, 11)[visual_ex,:,:].reshape(55) - y2 = data[:,1].reshape(5, 5, 11)[visual_ex,:,:].reshape(55) - z2 = 559.*.01/data[:,3].reshape(5, 5, 11)[visual_ex,:,:].reshape(55) + x2 = data[:,2].reshape(5, 4, 11)[visual_ex,:,:].reshape(44) + y2 = np.log10(data[:,1].reshape(5, 4, 11)[visual_ex,:,:].reshape(44)) + z2 = 559.*.08/data[:,3].reshape(5, 4, 11)[visual_ex,:,:].reshape(44) ax2.scatter(x2,y2,z2) ax2.set_xlabel('age') - ax2.set_ylabel('apera') + ax2.set_ylabel('log(apera)') ax2.set_zlabel('starformation rate in M_sun/year') + # ax2.set_zlim(0.,1.) + ax2.set_title('Starformation as a function of\n visual extinction apperaturesize and age at Av=10') - plt.savefig('plot/age-apera.svg') + # plt.savefig('plot/age-apera.png') @@ -78,7 +81,7 @@ def cmd(folder, av, apera, age): xmin = -1. xmax = 8. ymin = 15. - ymax = 0. + ymax = -1. selectmax = 0. selectmin = 8. @@ -86,8 +89,8 @@ def cmd(folder, av, apera, age): #headers = f.readline().strip().split(',') #data = np.loadtxt(f) #f.close() - - hdulist = fits.open('%s/%s' %(folder,'sim_%s_%s_%s'%(av,apera,age))) + + hdulist = fits.open('%s/%s' %(folder,'sim_%03d_%06d_%07d'%(av,apera,age))) # av = hdulist[1].header['age'] data = hdulist[1].data diff --git a/setup.py b/setup.py index 253c865..bf19c80 100644 --- a/setup.py +++ b/setup.py @@ -9,8 +9,17 @@ def main(quiet=False): - '''Pulling all necessary files for using the starformation script + '''setup.main(quiet=False) +Pulling all necessary files for using the starformation script + +This script pulls the fits-files for the radiation models from the MPIA-folder of Thomas +Robitaille, if the files have been moved feel free to contact robitaille@mpia.de + +It also pulls the extinction law from +http://caravan.astro.wisc.edu/protostars/info.php?topic=sedfitter_results +if you have problems with them please contact +http://caravan.astro.wisc.edu/protostars/contact.php ''' if quiet: output_stream = StringIO() diff --git a/starformation.py b/starformation.py index 81de9c7..efb79dc 100755 --- a/starformation.py +++ b/starformation.py @@ -11,23 +11,40 @@ from astropy.io import fits -def main(massfunction = 0, starformationhistory = 0, A_v = 10.0, sfr = .01, apera = 24000, maxage = 2000000., distance = 8.0, appendix='default', quiet=0, precise=0): - '''Creates a sample of stars +def main(massfunction = 0, starformationhistory = 0, A_v = 10.0, sfr = .01, apera = 24000,\ + maxage = 2000000., distance = 8.0, appendix='default', quiet=0, precise=0): + '''main(massfunction = 0, starformationhistory = 0, A_v = 10.0, sfr = .01, apera = 24000,\ + maxage = 2000000., distance = 8.0, appendix='default', quiet=0, precise=0) + +Creates a sample of stars input: -A_v float value for the visual extinction -sfr float Star formation rate in M_sun/year, is assumed to be constant over maxage -apera float used aperature size for selecting the fluxes of the protostars -maxage float age of the star formation site, sfr is assumed to be constant -appendix String sets the outputfilename, default is the starting time (via time.time()) -quiet boolean if true (=1) surpresses all output -precise boolean if true (=1) sample single star till expected mass reached, else calculate - expected number (viaand sample as an array +massfunction distribution +starformation history distribution +A_v float value for the visual extinction +sfr float Star formation rate in M_sun/year, is assumed to be constant over maxage +apera float used aperature size for selecting the fluxes of the protostars +maxage float age of the star formation site, sfr is assumed to be constant +distance float distance to the simulated starformation site +appendix String sets the outputfilename, default is the starting time (via time.time()) +quiet boolean if true (=1) surpresses all standard output +precise boolean if true (=1) sample single star till expected mass reached, else calculate + expected number and sample as an array + +The distributions must provide an object which has the following members: + float cdf(float x) returns the integrated distribution up to x, is used to calculate + the expected mass + float _upperbound returns the upper limit of the distribution, is used to calculate + the expected mass + float[] sample(int n) returns an array of n floats, sampled from the distribution + float mean() returns the mean value of the distribution + output: -returns two files in the folder 'out/' the _settings file contains the used values of - A_v, sfr, apera, maxage, number of sampled stars, their cumulated mass and - the expected mass +returns a fits file in the out-folder, either using the appendix as filename or the time of the + starting of the script in order to prevent overwriting existing files + In the header of the fits-file are the values: A_v, sfr, apera, maxage and distance recorded + In the data part are the age, mass, modelnumber and the uncorrected and corrected fluxes ''' if quiet: diff --git a/starformationtest.py b/starformationtest.py new file mode 100755 index 0000000..aaf65d1 --- /dev/null +++ b/starformationtest.py @@ -0,0 +1,256 @@ +# -*- coding: utf-8 -*- +from __future__ import print_function +from time import time +from StringIO import StringIO +import sys +import distribution as dist +import numpy as np +import matplotlib.pyplot as plt +import scipy.spatial +from astropy.table import Table, Column +from astropy.io import fits + +@profile +def main(massfunction = 0, starformationhistory = 0, A_v = 10.0, sfr = .01, apera = 24000,\ + maxage = 2000000., distance = 8.0, appendix='default', quiet=0, precise=0): + '''Creates a sample of stars + +input: +massfunction distribution +starformation history distribution +A_v float value for the visual extinction +sfr float Star formation rate in M_sun/year, is assumed to be constant over maxage +apera float used aperature size for selecting the fluxes of the protostars +maxage float age of the star formation site, sfr is assumed to be constant +distance float distance to the simulated starformation site +appendix String sets the outputfilename, default is the starting time (via time.time()) +quiet boolean if true (=1) surpresses all standard output +precise boolean if true (=1) sample single star till expected mass reached, else calculate + expected number and sample as an array + +The distributions must provide an object which has the following members: + float cdf(float x) returns the integrated distribution up to x, is used to calculate + the expected mass + float _upperbound returns the upper limit of the distribution, is used to calculate + the expected mass + float[] sample(int n) returns an array of n floats, sampled from the distribution + float mean() returns the mean value of the distribution + + +output: +returns a fits file in the out-folder, either using the appendix as filename or the time of the + starting of the script in order to prevent overwriting existing files + In the header of the fits-file are the values: A_v, sfr, apera, maxage and distance recorded + In the data part are the age, mass, modelnumber and the uncorrected and corrected fluxes +''' + + if quiet: + output_stream = StringIO() + else: + output_stream = sys.stdout + + t0 = time() + if appendix=='default': # making sure not to overwrite former output + appendix=t0 # by using the starting time as an unique id + #parameter settings + k_v = 211.4 # opacity in v_band in cm^2/g + # wavelength of the corresponding filterband in microns + wavelength = [1.235, 1.662, 2.159, 3.550, 4.493, 5.731, 7.872, 23.68, 71.42, 155.9] + models = ['2H', '2J', '2K', 'I1', 'I2', 'I3', 'I4', 'M1', 'M2', 'M3'] + + + if massfunction == 0 and starformationhistory == 0: + # star mass function + def f(x): # Kouper IMF + #http://adsabs.harvard.edu/abs/2001MNRAS.322..231K + if x<.01: + return 0 + elif x < .08: + return 62.46192*x**-.3 + elif x < .5: + return 1.413352*x**-1.8 + elif x < 1.: + return x**-2.3 # also value 2.7 given eq.(6) or eq (2) + else: + return x**-2.3 + f = np.vectorize(f) + massfunction = dist.distribution(f, .1, 50.) + + # star formation history + def g(x): + return sfr + g = np.vectorize(g) + starformationhistory = dist.distribution(g, 10000., maxage) + + + + cumass = 0. #sampled mass + exmass = starformationhistory.cdf()(starformationhistory._upperbound) #expected mass formed + stars = [] #storing for the sample + n = 0 + + t1 = time() # startup completed + + if precise == True: + while cumass < exmass: + mass, age = massfunction.sample(1)[0], starformationhistory.sample(1)[0] + cumass = cumass + mass + stars.append([n, age, mass]) + if n % 10000 == 0: + print (n, cumass, file=output_stream) #reporting progress + n = n+1 + else: + n = int(exmass/ massfunction.mean()) + mass, age = massfunction.sample(n), starformationhistory.sample(n) + cumass = np.sum(mass) + stars = [[i, age[i], mass[i]] for i in range(n)] + + print ('number of sampled stars: %s' %n , file=output_stream) + print ('mass of sampled stars: %s' % cumass , file=output_stream) + print ('mean mass: %s' % (cumass/n), file=output_stream) + print ('expected mass of stars: %s' % exmass , file=output_stream) + t2 = time() # sampleing completed + + + # python code for model contact + #initial parameters + model = [ fits.open('models/%s.fits' % mod) for mod in models ] # fits-data for the model + param = fits.open('models/parameters.fits.gz') # modelparameter + app_num = [ np.interp(apera, model[i][2].data.field(0), range(model[i][2].data.field(0).size)) for i in range(len(models)) ] + + # to do: + # check for interpolation of aperature size + + # sampling viewing angle + angle = np.random.random_integers(0,9,len(stars)) + #reading model grid + mass = param[1].data['MASSC'][::10] + age = param[1].data['TIME'][::10] + grid = np.vstack([age, mass]).transpose() + + #converting to logspace + stars = np.asarray(stars) + grid = np.log10(grid) + stars[:,1:] = np.log10(stars[:,1:]) + + output = stars.tolist() #creating output + + #normalizing for nearest neighbor search + grid[0,:] = grid[0,:]/(grid[0,:].max() - grid[0,:].min()) + grid[1,:] = grid[1,:]/(grid[1,:].max() - grid[1,:].min()) + stars[1,:] = stars[1,:]/(grid[0,:].max() - grid[0,:].min()) + stars[2,:] = stars[2,:]/(grid[1,:].max() - grid[1,:].min()) + + t3 = time() #model data load complete + + tree = scipy.spatial.cKDTree(grid,leafsize=10) #search tree + matches = [tree.query(star[1:] , k=1)[1] for star in stars] #saves matches with (dist, index) + + t4 = time() #matching sample to data complete + + # extracting fluxes + fluxes = [0 for j in range(len(models)) ] + indices = 10*np.asarray(matches) + angle + for j in range(len(models)): + fluxes[j] = model[j][1].data[indices]['TOTAL_FLUX'][:,app_num[j]] + + + + # applying extinction + extinction = np.loadtxt('models/extinction_law.ascii') + k_lambda = np.interp(wavelength, extinction[:,0], extinction[:,1]) + correctionfactor = 10.**(-.4 * A_v * k_lambda / k_v) + + newfluxes = [0 for j in range(len(models)) ] + for j in range(len(models)): + newfluxes[j] = np.asarray(fluxes[j]) * correctionfactor[j] * (1./distance)**2 + + + t5 = time() #extracting fluxes complete + + # saving data to output: #, log10(age), log10(mass), modelmatch, (flux, flux_error, corrected_flux, corrected_flux_error) for each model + fluxes = np.asarray(fluxes) + newfluxes = np.asarray(newfluxes) + output = np.vstack([np.asarray(output).transpose(), matches, fluxes, newfluxes]).transpose() + + # creating the output file + #head = ['#', 'age', 'mass', 'model'] + #for mod in models: + #head.append('flux %s' % mod) + #for mod in models: + #head.append('corrected_flux %s' % mod) + #f = open('out/%s' % appendix, 'w') + #f.write( ','.join(head)+'\n' ) + ##np.savetxt(f, output) + #np.savetxt(f, output) + #f.close() + + # create table + + t = Table() + t.add_column(Column(name='age', data=output[:,1])) + t.add_column(Column(name='mass', data=output[:,2])) + t.add_column(Column(name='model', data=output[:,3])) + for i in range(len(models)): + t.add_column(Column(name='flux %s' % models[i], data=output[:,4+i])) + for i in range(len(models)): + t.add_column(Column(name='cflux %s' % models[i], data=output[:,4+len(models)+i])) + + header = fits.Header() + header['AV'] = A_v + header['SFR'] = sfr + header['APPERA'] = apera + header['MAXAGE'] = maxage + header['DIST'] = distance + + fits.writeto('out/%s' % appendix, np.array(t), header, clobber=True) + + + ## creating the settings file + #f = open('out/%s_settings' % appendix, 'w') + #settings = '%s #visual extinction A_v \n' % A_v + #settings = settings + '%s #star formation rate sfr \n' % sfr + #settings = settings + '%s #star formation time time \n' % maxage + #settings = settings + '%s #aperature size apera\n' % apera + #settings = settings + '%s #number of sampled stars \n' % len(stars) + #settings = settings + '%s #cumulated sampled mass \n' % cumass + #settings = settings + '%s #expected mass \n' % exmass + #f.write(settings) + #f.close() + + t6 = time() #saving complete + +# timing possibility for optimization efforts + + print( 'starting script at %f' %(t0), file=output_stream) + print( 'initializing %f' %(t1-t0), file=output_stream) + print( "sampleing %f" %(t2-t1), file=output_stream) + print( "model data load %f" %(t3-t2), file=output_stream) + print( "matching model %f" %(t4-t3), file=output_stream) + print( "extracting fluxes %f" %(t5-t4), file=output_stream) + print( "saving %f" %(t6-t5), file=output_stream) + print( "________________________", file=output_stream) + print( "total runtime %f" %(t6-t0), file=output_stream) + print( "finishing script %f" %t6, file=output_stream) + +def f(x): # Kouper IMF + #http://adsabs.harvard.edu/abs/2001MNRAS.322..231K + if x<.01: + return 0 + elif x < .08: + return 62.46192*x**-.3 + elif x < .5: + return 1.413352*x**-1.8 + elif x < 1.: + return x**-2.3 # also value 2.7 given eq.(6) or eq (2) + else: + return x**-2.3 +f = np.vectorize(f) +massfunction = dist.distribution(f, .1, 50.) + + # star formation history +def g(x): + return .08 +g = np.vectorize(g) +starformationhistory = dist.distribution(g, 10000., 2000000.) +main(sfr = .08, massfunction = massfunction, starformationhistory = starformationhistory)