Skip to content

Commit

Permalink
Several updates, functionality improvements, and fixed a few issues
Browse files Browse the repository at this point in the history
  • Loading branch information
@tomasstolker
tomasstolker committed Sep 18, 2024
1 parent ec8f9c4 commit c37667c
Show file tree
Hide file tree
Showing 16 changed files with 636 additions and 429 deletions.
10 changes: 5 additions & 5 deletions docs/tutorials/read_isochrone.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@
"id": "7f6f8c20",
"metadata": {},
"source": [
"In this tutorial, we will work with data from an evolutionary model. We will extract an isochrone and cooling curve by interpolating the data at a fixed age and mass, respectively. We will also compute synthetic photometry for a given age and mass by using the associated grid with model spectra."
"In this tutorial, we will work with data from an evolutionary model. We will extract an isochrone and cooling track by interpolating the data at a fixed age and mass, respectively. We will also compute synthetic photometry for a given age and mass by using the associated grid with model spectra."
]
},
{
Expand Down Expand Up @@ -493,15 +493,15 @@
"id": "65c4ab2e",
"metadata": {},
"source": [
"## Extracting a cooling curve"
"## Extracting a cooling track"
]
},
{
"cell_type": "markdown",
"id": "7c6d9faa",
"metadata": {},
"source": [
"Similarly to extracting an isochrone, we can extract a cooling curve with the [get_cooling_curve](https://species.readthedocs.io/en/latest/species.read.html#species.read.read_isochrone.ReadIsochrone.get_cooling_curve) method, which will interpolate the evolutionary data at a fixed mass and a range of ages. Instead of providing an `numpy` array with ages, we can also set the argument of `ages` to `None`. In that case it use the ages that are available in the original data."
"Similarly to extracting an isochrone, we can extract a cooling track with the [get_cooling_track](https://species.readthedocs.io/en/latest/species.read.html#species.read.read_isochrone.ReadIsochrone.get_cooling_track) method, which will interpolate the evolutionary data at a fixed mass and a range of ages. Instead of providing an `numpy` array with ages, we can also set the argument of `ages` to `None`. In that case it use the ages that are available in the original data."
]
},
{
Expand All @@ -511,7 +511,7 @@
"metadata": {},
"outputs": [],
"source": [
"cooling_box = read_iso.get_cooling_curve(mass=10.,\n",
"cooling_box = read_iso.get_cooling_track(mass=10.,\n",
" ages=None,\n",
" filters_color=None,\n",
" filter_mag=None)"
Expand All @@ -522,7 +522,7 @@
"id": "a08d4f32",
"metadata": {},
"source": [
"The cooling curve that is returned by [get_cooling_curve](https://species.readthedocs.io/en/latest/species.read.html#species.read.read_isochrone.ReadIsochrone.get_cooling_curve) is stored in a [CoolingBox](https://species.readthedocs.io/en/latest/species.core.html#species.core.box.CoolingBox). Let's have a look at the content by again using the [open_box](https://species.readthedocs.io/en/latest/species.core.html#species.core.box.Box.open_box) method."
"The cooling track that is returned by [get_cooling_track](https://species.readthedocs.io/en/latest/species.read.html#species.read.read_isochrone.ReadIsochrone.get_cooling_track) is stored in a [CoolingBox](https://species.readthedocs.io/en/latest/species.core.html#species.core.box.CoolingBox). Let's have a look at the content by again using the [open_box](https://species.readthedocs.io/en/latest/species.core.html#species.core.box.Box.open_box) method."
]
},
{
Expand Down
10 changes: 8 additions & 2 deletions species/core/box.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -546,13 +546,19 @@ def create_box(boxtype, **kwargs):
box = ResidualsBox()
box.name = kwargs["name"]
box.photometry = kwargs["photometry"]
box.spectrum = kwargs["spectrum"]
if "spectrum" in kwargs:
box.spectrum = kwargs["spectrum"]
if "model_name" in kwargs:
box.model_name = kwargs["model_name"]
if "chi2_red" in kwargs:
box.chi2_red = kwargs["chi2_red"]

elif boxtype == "samples":
box = SamplesBox()
box.spectrum = kwargs["spectrum"]
if "spectrum" in kwargs:
box.spectrum = kwargs["spectrum"]
if "model_name" in kwargs:
box.model_name = kwargs["model_name"]
box.parameters = kwargs["parameters"]
box.samples = kwargs["samples"]
box.ln_prob = kwargs["ln_prob"]
Expand Down
17 changes: 9 additions & 8 deletions species/data/companion_data/companion_data.json
View file
Original file line number Diff line number Diff line change
Expand Up @@ -112,7 +112,7 @@
"Currie et al. 2013",
"Gaia Early Data Release 3",
"Males et al. 2014",
"Mirek Brandt et al. 2021",
"Brandt et al. 2021",
"Miret-Roig et al. 2020",
"Nowak et al. 2020",
"Stolker et al. 2019",
Expand Down Expand Up @@ -152,7 +152,7 @@
],
"references": [
"Gaia Early Data Release 3",
"Mirek Brandt et al. 2021",
"Brandt et al. 2021",
"Miret-Roig et al. 2020",
"Nowak et al. 2020"
]
Expand Down Expand Up @@ -369,7 +369,7 @@
"Gaia Early Data Release 3",
"Galicher et al. 2011",
"Marois et al. 2010",
"Mirek Brandt et al. 2021",
"Brandt et al. 2021",
"Wang et al. 2018",
"Zurlo et al. 2016"
]
Expand Down Expand Up @@ -448,7 +448,7 @@
"Gaia Early Data Release 3",
"Galicher et al. 2011",
"Marois et al. 2010",
"Mirek Brandt et al. 2021",
"Brandt et al. 2021",
"Wang et al. 2018",
"Zurlo et al. 2016"
]
Expand Down Expand Up @@ -527,7 +527,7 @@
"Gaia Early Data Release 3",
"Galicher et al. 2011",
"Marois et al. 2010",
"Mirek Brandt et al. 2021",
"Brandt et al. 2021",
"Wang et al. 2018",
"Zurlo et al. 2016"
]
Expand Down Expand Up @@ -596,7 +596,7 @@
"Currie et al. 2014",
"Gaia Early Data Release 3",
"Marois et al. 2010",
"Mirek Brandt et al. 2021",
"Brandt et al. 2021",
"Wang et al. 2018",
"Zurlo et al. 2016"
]
Expand Down Expand Up @@ -947,8 +947,9 @@
0.03
],
"mass_companion": [
20.0,
20.0
12.7,
-1.0,
1.2
],
"accretion": false,
"age": [
Expand Down
100 changes: 60 additions & 40 deletions species/data/database.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2214,11 +2214,15 @@ def add_samples(
group_path = f"results/fit/{tag}/fixed_param/{key}"
hdf5_file.create_dataset(group_path, data=value)

if "spec_type" in attr_dict:
dset.attrs["type"] = attr_dict["spec_type"]
if "model_type" in attr_dict:
dset.attrs["model_type"] = attr_dict["model_type"]
elif "spec_type" in attr_dict:
dset.attrs["model_type"] = attr_dict["spec_type"]

if "spec_name" in attr_dict:
dset.attrs["spectrum"] = attr_dict["spec_name"]
if "model_name" in attr_dict:
dset.attrs["model_name"] = attr_dict["model_name"]
elif "spec_name" in attr_dict:
dset.attrs["model_name"] = attr_dict["spec_name"]

dset.attrs["n_param"] = int(len(modelpar))
dset.attrs["sampler"] = str(sampler)
Expand Down Expand Up @@ -2590,8 +2594,15 @@ def get_mcmc_spectra(
hdf5_file = h5py.File(self.database, "r")
dset = hdf5_file[f"results/fit/{tag}/samples"]

spectrum_type = dset.attrs["type"]
spectrum_name = dset.attrs["spectrum"]
if "model_type" in dset.attrs:
model_type = dset.attrs["model_type"]
else:
model_type = dset.attrs["type"]

if "model_name" in dset.attrs:
model_name = dset.attrs["model_name"]
else:
model_name = dset.attrs["spectrum"]

if "n_param" in dset.attrs:
n_param = dset.attrs["n_param"]
Expand Down Expand Up @@ -2635,7 +2646,7 @@ def get_mcmc_spectra(
elif dset.attrs[f"parameter{i}"][:9] == "corr_amp_":
ignore_param.append(dset.attrs[f"parameter{i}"])

if spec_res is not None and spectrum_type == "calibration":
if spec_res is not None and model_type == "calibration":
warnings.warn(
"Smoothing of the spectral resolution is not "
"implemented for calibration spectra."
Expand Down Expand Up @@ -2678,24 +2689,24 @@ def get_mcmc_spectra(

hdf5_file.close()

if spectrum_type == "model":
if spectrum_name == "planck":
if model_type in ["model", "atmosphere"]:
if model_name == "planck":
from species.read.read_planck import ReadPlanck

readmodel = ReadPlanck(wavel_range)

elif spectrum_name == "powerlaw":
elif model_name == "powerlaw":
pass

else:
from species.read.read_model import ReadModel

readmodel = ReadModel(spectrum_name, wavel_range=wavel_range)
readmodel = ReadModel(model_name, wavel_range=wavel_range)

elif spectrum_type == "calibration":
elif model_type == "calibration":
from species.read.read_calibration import ReadCalibration

readcalib = ReadCalibration(spectrum_name, filter_name=None)
readcalib = ReadCalibration(model_name, filter_name=None)

boxes = []

Expand Down Expand Up @@ -2724,15 +2735,15 @@ def get_mcmc_spectra(
elif "distance" not in model_param and distance is not None:
model_param["distance"] = distance

if spectrum_type == "model":
if spectrum_name == "planck":
if model_type in ["model", "atmosphere"]:
if model_name == "planck":
specbox = readmodel.get_spectrum(
model_param,
spec_res,
wavel_resample=wavel_resample,
)

elif spectrum_name == "powerlaw":
elif model_name == "powerlaw":
if wavel_resample is not None:
warnings.warn(
"The 'wavel_resample' parameter is not support by the "
Expand Down Expand Up @@ -2779,7 +2790,7 @@ def get_mcmc_spectra(

specbox = create_box(
boxtype="model",
model=spectrum_name,
model=model_name,
wavelength=specbox_0.wavelength,
flux=flux_comb,
parameters=model_param,
Expand All @@ -2794,7 +2805,7 @@ def get_mcmc_spectra(
ext_filter=ext_filter,
)

elif spectrum_type == "calibration":
elif model_type == "calibration":
specbox = readcalib.get_spectrum(model_param)

boxes.append(specbox)
Expand Down Expand Up @@ -2864,8 +2875,15 @@ def get_mcmc_photometry(
elif "nparam" in dset.attrs:
n_param = dset.attrs["nparam"]

spectrum_type = dset.attrs["type"]
spectrum_name = dset.attrs["spectrum"]
if "model_type" in dset.attrs:
model_type = dset.attrs["model_type"]
else:
model_type = dset.attrs["type"]

if "model_name" in dset.attrs:
model_name = dset.attrs["model_name"]
else:
model_name = dset.attrs["spectrum"]

if "binary" in dset.attrs:
binary = dset.attrs["binary"]
Expand Down Expand Up @@ -2901,8 +2919,8 @@ def get_mcmc_photometry(
for i in range(n_param):
param.append(dset.attrs[f"parameter{i}"])

if spectrum_type == "model":
if spectrum_name == "powerlaw":
if model_type in ["model", "atmosphere"]:
if model_name == "powerlaw":
from species.phot.syn_phot import SyntheticPhotometry

synphot = SyntheticPhotometry(filter_name)
Expand All @@ -2911,12 +2929,12 @@ def get_mcmc_photometry(
else:
from species.read.read_model import ReadModel

readmodel = ReadModel(spectrum_name, filter_name=filter_name)
readmodel = ReadModel(model_name, filter_name=filter_name)

elif spectrum_type == "calibration":
elif model_type == "calibration":
from species.read.read_calibration import ReadCalibration

readcalib = ReadCalibration(spectrum_name, filter_name=filter_name)
readcalib = ReadCalibration(model_name, filter_name=filter_name)

mcmc_phot = np.zeros((samples.shape[0]))

Expand All @@ -2936,8 +2954,8 @@ def get_mcmc_photometry(
elif "distance" not in model_param and distance is not None:
model_param["distance"] = distance

if spectrum_type == "model":
if spectrum_name == "powerlaw":
if model_type in ["model", "atmosphere"]:
if model_name == "powerlaw":
from species.util.model_util import powerlaw_spectrum

pl_box = powerlaw_spectrum(synphot.wavel_range, model_param)
Expand Down Expand Up @@ -3004,7 +3022,7 @@ def get_mcmc_photometry(
else:
mcmc_phot[i], _ = readmodel.get_flux(model_param)

elif spectrum_type == "calibration":
elif model_type == "calibration":
if phot_type == "magnitude":
app_mag, _ = readcalib.get_magnitude(model_param=model_param)
mcmc_phot[i] = app_mag[0]
Expand Down Expand Up @@ -3229,7 +3247,10 @@ def get_samples(
for item in dset.attrs:
attributes[item] = dset.attrs[item]

spectrum = dset.attrs["spectrum"]
if "model_name" in dset.attrs:
model_name = dset.attrs["model_name"]
else:
model_name = dset.attrs["spectrum"]

if "n_param" in dset.attrs:
n_param = dset.attrs["n_param"]
Expand All @@ -3239,7 +3260,6 @@ def get_samples(
if "ln_evidence" in dset.attrs:
ln_evidence = dset.attrs["ln_evidence"]
else:
# For backward compatibility
ln_evidence = None

param = []
Expand All @@ -3249,8 +3269,6 @@ def get_samples(
print(f" - {param[-1]}")

# Printing uniform and normal priors
# Check if attributes are present for
# backward compatibility

uniform_priors = {}
normal_priors = {}
Expand Down Expand Up @@ -3337,7 +3355,7 @@ def get_samples(

return create_box(
"samples",
spectrum=spectrum,
model_name=model_name,
parameters=param,
samples=samples,
ln_prob=ln_prob,
Expand Down Expand Up @@ -3377,7 +3395,6 @@ def get_evidence(self, tag: str) -> Tuple[float, float]:
if "ln_evidence" in dset.attrs:
ln_evidence = dset.attrs["ln_evidence"]
else:
# For backward compatibility
ln_evidence = (None, None)

return ln_evidence[0], ln_evidence[1]
Expand Down Expand Up @@ -3427,12 +3444,16 @@ def get_pt_profiles(
hdf5_file = h5py.File(self.database, "r")
dset = hdf5_file[f"results/fit/{tag}/samples"]

spectrum = dset.attrs["spectrum"]
if "model_name" in dset.attrs:
model_name = dset.attrs["model_name"]
else:
model_name = dset.attrs["spectrum"]

pt_profile = dset.attrs["pt_profile"]

if spectrum != "petitradtrans":
if model_name != "petitradtrans":
raise ValueError(
f"The model spectrum of the posterior samples is '{spectrum}' "
f"The model spectrum of the posterior samples is '{model_name}' "
f"instead of 'petitradtrans'. Extracting P-T profiles is "
f"therefore not possible."
)
Expand Down Expand Up @@ -3866,8 +3887,8 @@ def add_retrieval(

dset = hdf5_file.create_dataset(f"results/fit/{tag}/samples", data=samples)

dset.attrs["type"] = "model"
dset.attrs["spectrum"] = "petitradtrans"
dset.attrs["model_type"] = "retrieval"
dset.attrs["model_name"] = "petitradtrans"
dset.attrs["n_param"] = len(parameters)

if "parallax" in radtrans:
Expand Down Expand Up @@ -4356,7 +4377,6 @@ def get_retrieval_spectra(
temp_nodes = dset.attrs["temp_nodes"]

else:
# For backward compatibility
temp_nodes = None

# Get distance
Expand Down
Loading

1 comment on commit c37667c

@gabrielastro
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That looks nice! Thank you. Small ones:

  • "has to to be build manually" → "has to to be built manually" in fit_evolution.py
  • "and radius parameters. Unless" → "and radius parameters, unless" in read_model.py

Please sign in to comment.