Skip to content

Commit

Permalink
feature: batch processing of lists of objects is now working.
Browse files Browse the repository at this point in the history 
Some further refactoring of the speedyfit.py main module.
  • Loading branch information
@vosjo
vosjo committed Dec 2, 2020
1 parent 11cae62 commit c508710
Show file tree
Hide file tree
Showing 3 changed files with 105 additions and 68 deletions.
2 changes: 1 addition & 1 deletion setup.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@
test_suite='pytest.collector',
tests_require=['pytest'],
entry_points={
'console_scripts': ['speedyfit=speedyfit.speedyfit:main'],
'console_scripts': ['speedyfit=speedyfit.speedyfit:main', 'speedyfit-batch=speedyfit.speedyfit_batch:main'],
},
classifiers=[
"Development Status :: 3 - Alpha",
Expand Down
127 changes: 69 additions & 58 deletions speedyfit/speedyfit.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -203,70 +203,17 @@ def fit_sed(setup, photbands, obs, obs_err):
if 'g2' in names: names.remove('g2')
samples = samples[names]

return results, samples, constraints, gridnames

def main():
parser = argparse.ArgumentParser()
parser.add_argument('filename', action="store", type=str, help='use setup given in this file')
parser.add_argument("-empty", dest='empty', type=str, default=None,
help="Create empty setup file ('single' or 'double')")
parser.add_argument('--phot', dest='photometry', action='store_true',
help='When creating a new setupfile, use this option to also download photometry from Vizier and Tap archives.')
parser.add_argument('--noplot', dest='noplot', action='store_true',
help="Don't show any plots, only store to disk.")
args, variables = parser.parse_known_args()

if args.empty is not None:

objectname = args.filename
filename = objectname + '_single.yaml' if args.empty == 'single' else objectname + '_binary.yaml'

plx, e_plx = photometry_query.get_parallax(objectname)

out = default_single if args.empty == 'single' else default_double
out = out.replace('<photfilename>', objectname + '.phot')
out = out.replace('<objectname>', objectname)
out = out.replace('<plx>', str(plx))
out = out.replace('<e_plx>', str(e_plx))

ofile = open(filename, 'w')
ofile.write(out)
ofile.close()

if args.photometry:
photometry = photometry_query.get_photometry(objectname, filename=objectname + '.phot')

sys.exit()

if args.filename is None:
print("Nothing to do")
sys.exit()

# -- load the setup file
setupfile = open(args.filename)
setup = yaml.safe_load(setupfile)
setupfile.close()

# -- obtain the observations
photbands, obs, obs_err = get_observations(setup)

# -- perform the SED fit
results, samples, constraints, gridnames = fit_sed(setup, photbands, obs, obs_err)

print("================================================================================")
print("")
print("Resulting parameter values and errors:")

percentiles = setup.get('percentiles', [16, 50, 84])
pc = np.percentile(samples.view(np.float64).reshape(samples.shape + (-1,)), percentiles, axis=0)
pars = {}
for p, v, e1, e2 in zip(samples.dtype.names, pc[1], pc[1] - pc[0], pc[2] - pc[1]):
results[p] = [results[p], v, e1, e2]
pars[p] = v

print(" Par Best Pc emin emax")
for p in samples.dtype.names:
print(" {:10s} = {} {} -{} +{}".format(p, *plotting.format_parameter(p, results[p])))
return results, samples, constraints, gridnames


def write_results(setup, results, samples, obs, obs_err, photbands):

outpars, outvals = [], []
for par in ['teff', 'logg', 'L', 'rad']:
Expand Down Expand Up @@ -302,8 +249,10 @@ def main():
fileio.write_summary2hdf5(setup['objectname'], samples, obs, obs_err, photbands, pars=results,
grids=setup['grids'], filename=h5file)

# -- Plotting

def plot_results(setup, results, samples, constraints, gridnames, obs, obs_err, photbands):

# check for 10 possible plots. Should be enough for now.
for i in range(10):

pindex = 'plot' + str(i)
Expand Down Expand Up @@ -354,6 +303,68 @@ def main():
if not setup[pindex].get('path', None) is None:
pl.savefig(setup[pindex].get('path', 'distribution.png'))


def main():
parser = argparse.ArgumentParser()
parser.add_argument('filename', action="store", type=str, help='use setup given in this file')
parser.add_argument("-empty", dest='empty', type=str, default=None,
help="Create empty setup file ('single' or 'double')")
parser.add_argument('--phot', dest='photometry', action='store_true',
help='When creating a new setupfile, use this option to also download photometry from Vizier and Tap archives.')
parser.add_argument('--noplot', dest='noplot', action='store_true',
help="Don't show any plots, only store to disk.")
args, variables = parser.parse_known_args()

if args.empty is not None:

objectname = args.filename
filename = objectname + '_single.yaml' if args.empty == 'single' else objectname + '_binary.yaml'

plx, e_plx = photometry_query.get_parallax(objectname)

out = default_single if args.empty == 'single' else default_double
out = out.replace('<photfilename>', objectname + '.phot')
out = out.replace('<objectname>', objectname)
out = out.replace('<plx>', str(plx))
out = out.replace('<e_plx>', str(e_plx))

ofile = open(filename, 'w')
ofile.write(out)
ofile.close()

if args.photometry:
photometry = photometry_query.get_photometry(objectname, filename=objectname + '.phot')

sys.exit()

if args.filename is None:
print("Nothing to do")
sys.exit()

# -- load the setup file
setupfile = open(args.filename)
setup = yaml.safe_load(setupfile)
setupfile.close()

# -- obtain the observations
photbands, obs, obs_err = get_observations(setup)

# -- perform the SED fit
results, samples, constraints, gridnames = fit_sed(setup, photbands, obs, obs_err)

# -- write the results
write_results(setup, results, samples, obs, obs_err, photbands)

# -- create plots
plot_results(setup, results, samples, constraints, gridnames, obs, obs_err, photbands)

print("================================================================================")
print("")
print("Resulting parameter values and errors:")
print(" Par Best Pc emin emax")
for p in samples.dtype.names:
print(" {:10s} = {} {} -{} +{}".format(p, *plotting.format_parameter(p, results[p])))

if not args.noplot:
pl.show()

Expand Down
44 changes: 35 additions & 9 deletions speedyfit/speedyfit_batch.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@
from astropy.coordinates import Angle

from speedyfit import photometry_query
from speedyfit.speedyfit import fit_sed, get_observations
from speedyfit.speedyfit import fit_sed, get_observations, write_results, plot_results
from speedyfit.default_setup import default_double, default_single


Expand Down Expand Up @@ -52,10 +52,10 @@ def convert_dec(dec):
def read_setup(setup):

if os.path.isfile(setup):
infile = os.open(setup, 'r')
setup = infile.read_lines()
infile = open(setup, 'r')
setup = infile.readlines()
infile.close()
setup = '/n'.join(setup)
setup = ''.join(setup)

elif setup == 'single':
setup = default_single
Expand All @@ -79,9 +79,11 @@ def prepare_setup(object_list, basedir, default_setup):

out = copy.copy(default_setup)
if '<photfilename>' in out:
out = out.replace('<photfilename>', objectname + '.phot')
photfilename_ = basedir + '/' + objectname + '/' + objectname + '.phot'
out = out.replace('<photfilename>', photfilename_)
if '<objectname>' in out:
out = out.replace('<objectname>', objectname)
objectname_ = basedir + '/' + objectname + '/' + objectname
out = out.replace('<objectname>', objectname_)
if '<plx>' in out:
out = out.replace('<plx>', str(plx))
out = out.replace('<e_plx>', str(e_plx))
Expand All @@ -104,13 +106,17 @@ def prepare_photometry(object_list, basedir, skip_existing=True):
os.mkdir(basedir+'/'+objectname)

try:
photometry = photometry_query.get_photometry(objectname, filename=objectname + '.phot')
filename = basedir + '/' + objectname + '/' + objectname + '.phot'
photometry = photometry_query.get_photometry(objectname, filename=filename)
except Exception as e:
print("Failed to obtain photometry for {}".format(objectname))
print(e)


def fit_seds(object_list, basedir):

all_results = {}

for objectname in object_list:

# read the setup
Expand All @@ -125,12 +131,32 @@ def fit_seds(object_list, basedir):
continue

try:
# get the photometry
photbands, obs, obs_err = get_observations(setup)
fit_sed(setup, photbands, obs, obs_err)

# do the SED fit
results, samples, constraints, gridnames = fit_sed(setup, photbands, obs, obs_err)

# add results to results dic
all_results.setdefault('objectname', []).append(objectname)
for k, v in results.items():
all_results.setdefault(k, []).append(v[0])
all_results.setdefault(k+'_err', []).append(( v[2] + v[3] ) / 2.0)

# write the results
write_results(setup, results, samples, obs, obs_err, photbands)

# make any requested plots
plot_results(setup, results, samples, constraints, gridnames, obs, obs_err, photbands)

except Exception as e:
print("Could not fit {}".format(objectname))
print(e)

all_results = pd.DataFrame(data=all_results)
all_results.to_csv('results.csv')
print(all_results)



def main():
Expand Down Expand Up @@ -168,4 +194,4 @@ def main():

# fit all systems if requested
if args.fit:
pass
fit_seds(object_list, basedir)

0 comments on commit c508710

Please sign in to comment.