This page will guide you on the different options the NEST package provides for dating stars.

You can download these docs as a jupyter notebook here tutorial.ipynb

Start by importing the package:

import NEST

List of available models

Each model is a neural network trained on a different evolutionary stellar grid.

NEST.available_models()

BaSTIModel (http://basti-iac.oa-abruzzo.inaf.it/)
PARSECModel (https://stev.oapd.inaf.it/PARSEC/)
MISTModel (https://waps.cfa.harvard.edu/MIST/)
GenevaModel (https://www.unige.ch/sciences/astro/evolution/en/database/syclist/)
DartmouthModel (https://rcweb.dartmouth.edu/stellar/)
YaPSIModel (http://www.astro.yale.edu/yapsi/)
BaSTI_HSTModel (http://basti-iac.oa-abruzzo.inaf.it/)

Loading a model

By default, the age models use sklearn to run the neural networks. If you do not have it installed, or do not want to use it, you can pass the argument use_sklearn=False to switch to a numpy version.

You can also pass a verbose=False argument if you want the model to run quietly.

Note

Be aware though, in this mode it becomes significantly slower once you run the models on more than ~100 stars

age_model = NEST.BaSTIModel()

Tip

The color and magnitudes are, for most models, in Gaia passbands. Magnitude is \(M_G\) and color is \((G_{BP}-G_{RP})\). The only exception is the BaSTI_HSTModel(), which uses \(F814W\) as magnitude and \((F606W - F814W)\) as color.

If in doubt, print age_model to get a description of it.

Show code cell content

Hide code cell content

age_model

BaSTI Age Model, Gaia photometry: mag=MG, col=(GBP-GRP).

Estimating the age of a star

To estimate the age of a star, use the ages_prediction() method, and pass it a metallicity met, magnitude mag and color col. Note that, with only these values provided, the return value is a np.array of shape (1,1). The reason will be clear if you check later uses.

age_model.ages_prediction(met=0.0,mag=2.0,col=1.0)

array([[1.81505156]])

You can use the check_domain() method to verify that the star is within the bounds of the isochrones points used for training the neural networks. The function only checks for this in the HR Diagram, so you pass it a metallicity met, magnitude mag and color col.

age_model.check_domain(met=0.0,mag=2.0,col=1.0)

array([ True])

Estimating the age of a star with GBP & GRP

If you are using Gaia magnitudes and colors, you can also pass the red GRP and blue GBP magnitudes separately as additional inputs. This generally results in better estimations.

age_model.ages_prediction(met=0.0,mag=2.0,col=1.0,GBP=2.0,GRP=1.0)

Show code cell output

Hide code cell output

array([[1.34888727]])

Estimating the age of multiple stars

If you pass lists of values of size N as inputs instead of single values, the return value of the ages_prediction() method will be a np.array of shape (N,1). Again, the reason for the 2D shape will be clear by checking later uses.

Tip

If you have tqdm installed, you can display the progress of the age predictions. By default, this display is on, but you can turn it off when you initiate the age_model by passing the argument use_tqdm=False.

age_model.ages_prediction(
    met=[0.0,-1.0],
    mag=[2.0,3.0],
    col=[1.0,0.5]
)

Show code cell output

Hide code cell output

array([[1.81505156],
       [5.83193499]])

age_model.ages_prediction(
    met=[0.0,-1.0],
    mag=[2.0,3.0],
    col=[1.0,0.5],
    GBP=[2.0,1.0],
    GRP=[1.0,0.5]
)

Show code cell output

Hide code cell output

array([[1.34888727],
       [2.20902813]])

You can also check if the stars are within the training domain all at once by passing lists of metallicity met, magnitude mag & color values col.

age_model.check_domain(
    met=[0.0,-1.0],
    mag=[2.0,3.0],
    col=[1.0,0.5]
)

Show code cell output

Hide code cell output

array([ True,  True])

That way, you can crossmatch both outputs to only keep stars whose age estimates we can reliably trust:

met = [0.0,-1.0,-2.0,-3.0]
mag = [2.0,3.0,3.0,5.0]
col = [1.0,0.5,0.5,0.5]
GBP = [2.0,1.0,1.0,1.0]
GRP = [1.0,0.5,0.5,0.5]
ages = age_model.ages_prediction(met=met,mag=mag,col=col,GBP=GBP,GRP=GRP)
print(ages)
inside_domain = age_model.check_domain(met=met,mag=mag,col=col)
print(inside_domain)
ages = ages[inside_domain]#Only keep ages of stars within the training domain of the NN
print(ages)

[[1.34888727]
 [2.20902813]
 [3.38402904]
 [6.08718654]]
[ True  True  True False]
[[1.34888727]
 [2.20902813]
 [3.38402904]]

Estimating the age of a star with uncertainties

If you have uncertainty values for a star, you can pass them to the ages_prediction() method by adding values for emet, emag & ecol. The idea is that multiple age predictions will be made by offsetting the star’s metallicity met, magnitude mag & color col by adding a random gaussian offset of width (emet,emag,ecol). The number of realisations made is controlled by the parameter n. The output of the method will then be a np.array of shape (1,n).

age_model.ages_prediction(met=0.0,mag=2.0,col=1.0,emet=0.1,emag=0.05,ecol=0.05,n=10)

Show code cell output

Hide code cell output

array([[6.60308496, 1.72800752, 1.66448607, 1.74689925, 2.08590546,
        1.68951072, 1.87893172, 1.79672179, 1.75266411, 1.80040514]])

Once the ages have been predicted for each star n times, you can call statistical methods model.mean_ages(), model.median_ages(), model.mode_ages() and model.std_ages()

print(age_model.mean_ages())
print(age_model.median_ages())
print(age_model.mode_ages())
print(age_model.std_ages())

[2.27466167]
[1.77469295]
[1.68651377]
[1.4472521]

Estimating the age of multiple stars with uncertainties

Finally, you can combine multiple (N) stars age estimations and added uncertainties, with n Monte Carlo realisations. The output of the ages_prediction() method will then be a np.array of shape (N,n). A value of n=100 seems to be a good middle ground between speed and accurate predictions. If the number of stars is low, n>=1000 gives enough realisations to draw a PDF.

age_model.ages_prediction(
    met=[0.0,-1.0],
    mag=[2.0,3.0],
    col=[1.0,0.5],
    emet=[0.1,0.1],
    emag=[0.05,0.05],
    ecol=[0.05,0.05],
    n=10
)

Show code cell output

Hide code cell output

array([[3.80336953, 1.69983455, 4.90082245, 1.69778142, 1.59160446,
        1.71340115, 2.05606171, 1.7911002 , 1.78933451, 2.02847367],
       [4.86888566, 6.07496224, 5.08471375, 6.947272  , 5.2717213 ,
        5.54440037, 5.73574561, 5.86507676, 5.72814486, 5.99924428]])

print(age_model.mean_ages())
print(age_model.median_ages())
print(age_model.mode_ages())
print(age_model.std_ages())

[2.30717837 5.71201668]
[1.79021736 5.73194523]
[1.70239371 5.73512226]
[1.06045561 0.55593791]

Working with big datasets

In the event you are working with big datasets, you might not want to store every sampled age for every star, as you might run out of memory to store all this information. In this case, you can pass the argument store_samples=False to the age_prediction() method. Instead of returning a (N,n) np.array, the method will return a dict with keys mean, median, mode and std, which each contain a (N) np.array with the given statistical property.

Note

Keep in mind this will cause a significant slowdown, as the statistical properties have to be computed for each star individually before discarding all its samples.

age_model.ages_prediction(
    met=[0.0,-1.0],
    mag=[2.0,3.0],
    col=[1.0,0.5],
    emet=[0.1,0.1],
    emag=[0.05,0.05],
    ecol=[0.05,0.05],
    n=100,
    store_samples=False
)

Show code cell output

Hide code cell output

{'mean': array([2.068395  , 5.92732228]),
 'median': array([1.89050387, 5.87078113]),
 'mode': array([1.775, 5.825]),
 'std': array([0.56649702, 0.79865838])}

Coeval (cluster) population age

If you are dating a stellar population that is coeval, you might want to estimate the overall age of your population. To this end, after running the age prediction (with store_samples=True), you can call the population_age(). This works by multiplying individual age probability distribution functions (pdfs) together and returning the maximum value of this global pdf. By default, it uses 280 bins between 0 and 14 Gyr, but you can use the nbins, min_age and max_age arguments to change this.

age_model.ages_prediction(
    met=[0.0,0.0],
    mag=[2.0,3.0],
    col=[1.0,0.5],
    emet=[0.05,0.05],
    emag=[0.05,0.05],
    ecol=[0.05,0.05],
    n=10_000
)

age_model.population_age()

\[\tau=1.73_{-0.29}^{+-0.06}\,\mathrm{Gyr}\]

The population_age() method accepts the following arguments :

• nbins=280 Number of bins used to compute the population age.

• min_age=0 and max_age=14 Limits of the bins used to compute the population age.

• n_mc=100 Number of Monte Carlo samples for computing the population age.

• check_domain=True Whether to remove stars outside of the training domain.

• epsilon=None Hardcodes the ε value used in the computation of the population age.

Age model properties

• age_model.samples returns an (N,n,m) np.array of the latest samples used for the age predictions, if they were stored.

• age_model.ages returns an (N,n) np.array of the latest age predictions ran. N : number of stars, n : number of samples, m : inputs dimension (3 if using metallicity,magnitude & color, 5 if using GBP & GRP in addition).# Age model properties

• age_model.samples returns an (N,n,m) np.array of the latest samples used for the age predictions, if they were stored.

• age_model.ages returns an (N,n) np.array of the latest age predictions ran. N : number of stars, n : number of samples, m : inputs dimension (3 if using metallicity,magnitude & color, 5 if using GBP & GRP in addition).

#Here is an example on how you can retrieve the samples
#of the age predictions to investigate individual predictions
age_model.ages_prediction(
    met=[0.0,-1.0],
    mag=[2.0,3.0],
    col=[1.0,0.5],
    emet=[0.1,0.1],
    emag=[0.1,0.1],
    ecol=[0.1,0.1],
    n=100
)

star_i = 1#Select the second star
magnitudes = age_model.samples[star_i,:,1]
colors = age_model.samples[star_i,:,2]
ages = age_model.ages[star_i]

import matplotlib.pyplot as plt
plt.scatter(colors,magnitudes,c=ages)
plt.colorbar(label='Age [Gyr]')
plt.xlabel(r'$G_{BP}$ - $G_{RP}$ [mag]')
plt.ylabel(r'$M_G$ [mag]')
plt.ylim(plt.ylim()[::-1]);

Plotting an HR Diagram

Warning

You need the matplotlib package for this function.

To visualize your star sample in an HR Diagram, you can use the HR_diagram() method. It plots your stars, colored by their estimated age. It also plots isochrones from the evolutionary grid the age model Neural Network was trained on.

Note

By default, the curve data is not included in the NEST package, so the first time, you will be prompted to download them (27.9Mb).

The isochrone_met option lets you choose which set of isochrones will be plotted. You can choose between [M/H] = {-2,1,0}. This is a restricted range, so the displayed isochrones might not line up perfectly with your data points. Of course, the derived age was computed using the whole set of isochrones.

If you already ran ages_prediction(), the method will by default plot the stars colored by age. If you only want to plot isochrones, pass the plot_stars=False argument. The age_type option lets you choose the type of ages used to color the stars, between 'median', 'mean' and 'mode'. The check_domain option restricts the plotted stars to stars being in the domain of the training grid (true by default).

Tip

In an interactive matplotlib backend, you can hover over the isochrones to get their age.

import numpy as np
import matplotlib.pyplot as plt

met = np.random.uniform(-2, 0, size=100)
mag = np.random.uniform(-2, 4, size=100)
col = np.random.uniform( 0, 2, size=100)
age_model.ages_prediction(met=met,mag=mag,col=col)
age_model.HR_diagram();

Finally, if you are plotting a star population, the method will by default plot the isochrone with the same age as the overall age of the population (computed using age_model.population_age() if you stored samples, otherwise using the median of the type of individual ages you chose). You can disable this by setting plot_isochrone to False. If you want to plot a specific isochrone(s), you can set isochrone_ages to an age value in Gyr or a list of age values, as well as isochrone_ages_std for their associated uncertainties.

Warning

Currently, plotting exact individual isochrones is restricted to the BaSTI model, as it involves interpolation between a precomputed isochrone set, which is as of yet not fully implemented for other grids. For other models, the method will plot the closest isochrone in the available set of isochrones.

age_model = NEST.BaSTIModel()
fig,ax = age_model.HR_diagram(isochrone_ages=[0.5,1,3,5,10,12])

If you want to plot isochrones on one of your existing figures, you can pass the method your fig and ax objects, and nothing else. Otherwise, the HR_diagram() method returns a (fig,ax) tuple, so you can plot more objects on top by using the outputed figure and axis.

from matplotlib.colors import LogNorm

fig,ax = plt.subplots(figsize=(5,5))

age_model = NEST.BaSTIModel()
age_model.HR_diagram(fig=fig,ax=ax)

mag = np.random.normal(1  , .5  , size=10_000)
col = np.random.normal(.75, .125, size=10_000)

ax.hist2d(col,mag,bins=200,range=((-1,3),(-5,10)),norm=LogNorm())

plt.xlim(-1,3)
plt.ylim(10,-5)

plt.xlabel(r'$(G_{BP} - G_{RP})_0$ [mag]')
plt.ylabel(r'$M_G$ [mag]');

There are a lot of parameters to customize the look of the plotted HR Diagram. Here is a list :

• isochrone_met=0. The metallicity (between -2,1 and 0) of the plotted set of isochrones.

• plot_isochrone=True Whether to plot the fitted isochrone of the stellar population.

• plot_stars=True Whether to plot individual stars.

• n_mc=100 Number of Monte Carlo samples for computing the population age.

• epsilon=None Hardcodes the ε value used in the computation of the population age.

• plot_isochrone_uncertainty=True Whether to plot uncertainties in the isochrones as shaded areas.

• isochrone_ages=None, An age value or list of age values of isochrones to be displayed.

• isochrone_ages_std=None The associated uncertainties.

• age_type='median' Type of age to be used for coloring stars.

• check_domain=True Whether to remove stars outside of training domain.

• fig=None Existing figure to be used for plotting.

• ax=None Existing axis to be used for plotting.

• star_cmap=cmaps['tempo_R'] Colormap to be used for stars.

• isochrone_cmap=cmaps['acton'] Colormap to be used for isochrones.

• loc_legend='best' Position of the legend.

• axis_fontsize=12 Font size of the axes.

• legend_fontsize=10 Font size of the legend.

• colorbar_fontsize=10 Font size of the colorbar.

• selected_isochrone_linewidth=2 Line width of the colored isochrones.

• uncertainty_alpha=0.25 Opacity of the uncertainty shaded areas.

• colorbar=True Enables or disables the colorbar.

• colorbar_lims=None Limits of the colorbar for the stars scatter.

• colorbar_loc=right Position of the colorbar.

• colorbar_pad=0.02 Padding of the colorbar.

Any other keyword argument will be passed to the function plotting the set of isochrones, so for example, if you only want the best fitted isochrone but no other, you can pass the argument alpha=0.