<jupyter_start><jupyter_text>Tutorial Outline **Tutorial contents**:- Adding new yield tables- Choosing element set- Training of a neural network- Running MCMC analysis- Computing Bayes/LOO-CV scores- Plotting LOO-CV element predictionsThe above are based on the Philcox & Rybizki (2017, submitted) paper which should be cited when using this code. This is based on the $\mathit{Chempy}$ software, described in Rybizki et al. (2017, [arXiv:1702.08729](https://arxiv.org/abs/1702.08729)) and full tutorials for this can be found [here](https://github.com/jan-rybizki/Chempy/tree/master/tutorials).** Requirements**:Before running this tutorial, the $\mathit{ChempyScoring}$ code and its dependencies must be installed as described [here](https://github.com/oliverphilcox/ChempyScoring/blob/master/requirements.txt)The authors [<NAME>](mailto:<EMAIL>) and [<NAME>](<EMAIL>) are happy to assist with any problems which may arise  Step 1: Load Yield TablesFirst we must load in the Nucleosynthetic yield table to be tested. Here we will test the SN2 net yields of Portinari et al. (1998,  [arXiv:astro-ph/9711337](https://arxiv.org/abs/astro-ph/9711337)). These include elements for stars of mass 6-120Msun, as used in the Illustris simulation. These are constructed from the text files in the `Chempy/input/yields/Portinari_1998/` directory.The yield tables provide data for 11 masses and  four metalicities. The following are important considerations:- Yield tables must be **net** yields such that they sum to 0. These are obtained either directly from the input tables (as below) or from gross yields by subtracting the initial abundances (e.g. `net_yield` = [`gross_yield` - (`initial_mass_fraction`)$*$(`ejected_mass`)]/'`initial_mass`)- Output yields must be in fractions of the initial stellar mass (not the ejected mass)- `mass_in_remnants` field is the fraction of the initial stellar mass left as a remnant- `unprocessed_mass_in_winds` field is $1$-`mass_in_remnants-`$\Sigma$`(elemental yields)`. The last factor is to account for the yields not quite summing to zero due to missing elements etc.- The output table must be in the form shown below for compatibility with $\mathit{Chempy}$.To add this into *Chempy* we add a `Portinari_net()` function to the `SN2_feedback()` class in `yields.py` as shown below:<jupyter_code>## NB: The Portinari_net definition should be inserted into the yields.py file

from Chempy import localpath # For file locations
import numpy as np

class SN2_feedback(object):
    def __init__(self):   
        """
        This is the object that holds the feedback table for SN2 stars.
        The different methods load different tables from the literature. 
        They are in the input/yields/ folder.
        """

    def Portinari_net(self):
        '''
        Loading the yield table from Portinari1998. 
        These are presented as net yields in fractions of initial stellar mass.
        '''
        
        # Define metallicities in table
        self.metallicities = [0.0004,0.004,0.008,0.02,0.05] 
        
        # Load one table
        x = np.genfromtxt(localpath + 'input/yields/Portinari_1998/0.02.txt',names=True)
        
        # Define masses and elements in yield tables
        self.masses = list(x['Mass']) # In solar masses
        self.elements = list(x.dtype.names[3:])

        self.table = {} # Output dictionary for yield tables
        
        for metallicity in self.metallicities:
            additional_keys = ['Mass', 'mass_in_remnants','unprocessed_mass_in_winds'] 
            names = additional_keys + self.elements # These are fields in dictionary
            
            # Create empty record array of correct size
            base = np.zeros(len(self.masses))
            list_of_arrays = []
            for i in range(len(names)):
                list_of_arrays.append(base)
            yield_subtable = np.core.records.fromarrays(list_of_arrays,names=names)
            
            # Add mass field to subtable (in solar masses)
            yield_subtable['Mass'] = np.array(self.masses)
            
            # Read in yield tbale
            x = np.genfromtxt(localpath + 'input/yields/Portinari_1998/%s.txt' %(metallicity),names=True)
            
            # Read in element yields
            for item in self.elements:
                yield_subtable[item] = np.divide(x[item],x['Mass']) # Yields must be in mass fraction
            
            # Add fractional mass in remnants
            yield_subtable['mass_in_remnants'] = np.divide(x['Mass'] - x['ejected_mass'], x['Mass'])

            # Add unprocessed mass as 1-remnants (with correction if summed net yields are not exactly zero)
            for i,item in enumerate(self.masses):
                yield_subtable['unprocessed_mass_in_winds'][i] = 1. - (yield_subtable['mass_in_remnants'][i] + sum(list(yield_subtable[self.elements][i])))

            # Add subtable to output table
            self.table[metallicity] = yield_subtable<jupyter_output><empty_output><jupyter_text>We can now test the new yield table (using C as an example element).<jupyter_code># Define correct yield table
from Chempy.wrapper import SN2_feedback
basic_sn2 = SN2_feedback()
getattr(basic_sn2, 'Portinari_net')()

print("C Yields")
for metallicity in basic_sn2.metallicities:
    print("\n Metallicity = %.2e" %(metallicity))
    for i in range(len(basic_sn2.masses)):
        print("Mass = %d, Yield  = %.6e" %(basic_sn2.masses[i],basic_sn2.table[metallicity]['C'][i]))<jupyter_output>C Yields

 Metallicity = 4.00e-04
Mass = 6, Yield  = -3.710133e-05
Mass = 7, Yield  = -4.233929e-05
Mass = 9, Yield  = 9.881267e-03
Mass = 12, Yield  = 1.408225e-02
Mass = 15, Yield  = 1.415953e-02
Mass = 20, Yield  = 1.410735e-02
Mass = 30, Yield  = 1.160323e-02
Mass = 40, Yield  = 8.404800e-03
Mass = 60, Yield  = 8.438067e-03
Mass = 100, Yield  = 7.125290e-03
Mass = 120, Yield  = 1.630867e-02

 Metallicity = 4.00e-03
Mass = 6, Yield  = -3.069767e-04
Mass = 7, Yield  = -3.448071e-04
Mass = 9, Yield  = 9.192600e-03
Mass = 12, Yield  = 1.338642e-02
Mass = 15, Yield  = 1.512660e-02
Mass = 20, Yield  = 1.265245e-02
Mass = 30, Yield  = 9.859900e-03
Mass = 40, Yield  = 8.351425e-03
Mass = 60, Yield  = 6.609650e-03
Mass = 100, Yield  = 3.126120e-03
Mass = 120, Yield  = 3.591392e-04

 Metallicity = 8.00e-03
Mass = 6, Yield  = -5.534767e-04
Mass = 7, Yield  = -6.486429e-04
Mass = 9, Yield  = 2.709333e-02
Mass = 12, Yield  = 1.224050e-02
Mass = 15, Yield  = 1.476993e-02
Mass = 20, Yield  = 1.17[...]<jupyter_text>Pre-implemented Yield SetsIn the code, we have implemented the (net) tables from the following papers:CC-SN: - Chieffi & Limongi 2004 (`chieffi04_net`)- Ritter et al. 2017 (in prep.) (`NuGrid_net`)- Nomoto et al. 2013 (`Nomoto2013_net`)- Portinari 1998 (`Portinari_net`)- West & Heger 2017 (in prep.) (`West17_net`)- TNG tables (Pillepich et al. 2017) (`TNG_net`)SN Ia (Gross yields):- Iwamoto et al. 1999 (`Iwamoto`)- Theilemann et al. 2003 (`Thielemann`)- Seitenzahl et al. 2013 (`Seitenzahl`)- TNG tables (Pillepich et al. 2017) (`TNG`)AGB:- Karakas 2010 (`Karakas_net_yield`)- Karakas 2016 (`Karakas16_net`)- Ventura 2013 (`Ventura_net`)- TNG tables (Pillepich et al. 2017) (`TNG_net`)  Step 2: Choice of ElementsWe are free to use any set of chemical elements in this analysis. The set should be chosen to match the simulation. 28 elements up to Ge (excluding Li, B, Be) and the IllustrisTNG elements (Pillepich et al. 2017) were used in the Philcox & Rybizki (2017).To select the required elements we modify the `Chempy/parameter.py` file.Both `elements_to_trace` **and** `initial_neural_names` must be changed here.- `elements_to_trace` contains elements in the proto-solar data-file, including B, Be, Li and H which are not predicted directly by the neural network. - `initial_neural_names` is the elements predicted by the network only (as [X/Fe] or [Fe/H] abundances).In addition, if extra elements are added, it should be checked that they are predicted by the yield tables and feature in the `Chempy/input/stars/proto_sun_all.npy` observational data-set.<jupyter_code>## Modify these lines in Chempy/parameter.py

# This field should contain all required elements  (and B,Be,Li,H) in alphabetical order
elements_to_trace = ['Al', 'Ar', 'B', 'Be', 'C', 'Ca', 'Cl', 'Co', 'Cr', 'Cu', 'F', 'Fe', 'Ga', 'Ge', 'H', 'He', 'K', 'Li', 'Mg', 'Mn', 'N', 'Na', 'Ne', 'Ni', 'O', 'P', 'S', 'Sc', 'Si', 'Ti', 'V', 'Zn']

# This field contains names of elements predicted by neural network
initial_neural_names = ['Al', 'Ar', 'C', 'Ca', 'Cl', 'Co', 'Cr', 'Cu', 'F', 'Fe', 'Ga', 'Ge', 'He', 'K', 'Mg', 'Mn', 'N', 'Na', 'Ne', 'Ni', 'O', 'P', 'S', 'Sc', 'Si', 'Ti', 'V', 'Zn']<jupyter_output><empty_output><jupyter_text>Step 3: Create Neural Network DatasetNow that the yield set has been implemented, we must create a training data-set for the neural network.Firstly it is important to change the `parameter.py` file such that *Chempy* uses the correct yields. Here we add the new SN2 yield table name and set *Chempy* to use it by default.<jupyter_code>## Modify these lines in Chempy/parameter.py to add new yield set
yield_table_name_sn2_list = ['chieffi04','Nugrid','Nomoto2013','Portinari', 'chieffi04_net', 'Nomoto2013_net','NuGrid_net','West17_net','TNG_net']#'Frischknecht16_net'
yield_table_name_sn2_index = 3
yield_table_name_sn2 = yield_table_name_sn2_list[yield_table_name_sn2_index]<jupyter_output><empty_output><jupyter_text>We can test this as follows:<jupyter_code># Load parameter file
from Chempy.parameter import ModelParameters
a = ModelParameters()

# Print new SN2 yield table name
print(a.yield_table_name_sn2)<jupyter_output>Portinari_net<jupyter_text>We must also set the list of parameters to optimize over to include only the 5 free *Chempy* parameters (i.e. not $\beta$). This will be changed later in the analysis, but must be done at this point, else the `training_data()` routine will fail. **Here we can also change the free parameters (e.g. to fix $\alpha_\mathrm{IMF}$) **For compatibility reasons $\beta$ (when later added) is included in the SSP_parameters definition.<jupyter_code># Modify these lines in Chempy/parameter.py file

SSP_parameters =  [-2.29,-2.75] # Prior values
SSP_parameters_to_optimize = ['high_mass_slope','log10_N_0'] # Parameter names
assert len(SSP_parameters) == len(SSP_parameters_to_optimize)

ISM_parameters =  [-0.3,0.55,0.5]
ISM_parameters_to_optimize = ['log10_starformation_efficiency', 'log10_sfr_scale', 'outflow_feedback_fraction']
assert len(ISM_parameters) == len(ISM_parameters_to_optimize)<jupyter_output><empty_output><jupyter_text>We must also turn OFF the neural network predictions:<jupyter_code># In Chempy/parameter.py
UseNeural=False

# To test
a.UseNeural<jupyter_output><empty_output><jupyter_text>Data-sets can be created using the `Chempy.neural` module as follows and are saved in the `Neural/` directory. This uses multiprocessing to create a training data-set using 10 values of each of the 5 free *Chempy* parameters. (This value can be changed using `training_size` in `parameter.py`). These are written to file as `Neural/training_abundances.npy` (abundance output) and `Neural/training_norm_grid.npy` (normalised input)<jupyter_code>from Chempy.neural import training_data
#training_data()<jupyter_output><empty_output><jupyter_text>The above was run on a 64-core machine, taking 45 minutes. This is the most computationally expensive section, as $\mathit{Chempy}$ must be run $10^5$ times so multiprocessing is used here with as many cores as possible to speed this up.   Step 4: Train Neural Network Next we must train the neural network using the previously constructed data-sets. This can be simply done with the `Chempy/neural.py` `create_network()` function. This creates and trains a 30-neuron network over 1000 training epochs, using a learning rate of 0.007 by default (optimised via a validation data-set).<jupyter_code>from Chempy.neural import create_network
#create_network(Plot=True)<jupyter_output><empty_output><jupyter_text>The above was run on an 8-core machine taking 20 minutes. The trained network is saved as `Neural/neural_model.npz`. A loss plot is also produced, showing the network loss function against training epoch (if `Plot=True`).Using `Chempy/neural.py`'s `neural_output()` function we may simulate the output of *Chempy* for any set of input parameters;<jupyter_code>from Chempy.parameter import ModelParameters
a = ModelParameters()

from Chempy.neural import neural_output

# Here we compute the neural network predictions for the prior values of the free parameters (a.p0) 
# Here the output is for the Chieffi & Limongi yield tables and all 28 elements
output = neural_output(a.p0)

print("Element \t Predicted abundance")
print("------------------------------------")
for i in range(len(output)):
    print(a.initial_neural_names[i],"\t\t",output[i])<jupyter_output>Element 	 Predicted abundance
------------------------------------
Al 		 0.0393322117908
Ar 		 -0.0297753069018
C 		 -0.13662989818
Ca 		 -0.0828822019434
Cl 		 -0.557419770709
Co 		 -0.132588300145
Cr 		 0.0641349769778
Cu 		 -0.144458035335
F 		 -0.569666503541
Fe 		 0.256995791139
Ga 		 -0.0674121953433
Ge 		 -0.346140289736
He 		 -0.21917714508
K 		 -0.695086484994
Mg 		 -0.116784140666
Mn 		 0.200105075396
N 		 -0.0497552835193
Na 		 0.183123889506
Ne 		 0.169685670872
Ni 		 0.233162051334
O 		 0.00395271527296
P 		 -0.0362743018159
S 		 0.119659840693
Sc 		 -0.384143635707
Si 		 0.1511760692
Ti 		 -0.271989694516
V 		 -0.164913581524
Zn 		 -0.328304738689<jupyter_text>Step 5: Run MCMC analysisWe can now run MCMC on the yield set and produce optimal parameters and a corner plot. This is done as shown below, using a trained neural network.Whichever neural network is in the `Neural/` directory will be used.Before running this analysis we must add in the beta parameter into the `parameter.py` script and turn **on**  the neural network. Here we optimize for $\log_{10}\beta$ using a prior of $1.0\pm0.5$, but this can be changed into linear space using the `beta_param` identifier instead<jupyter_code># Modify these lines in Chempy/parameter.py file

SSP_parameters =  [1.0, -2.29,-2.75] # Prior values
SSP_parameters_to_optimize = ['log10_beta','high_mass_slope','log10_N_0'] # Parameter names (or use beta_param)
assert len(SSP_parameters) == len(SSP_parameters_to_optimize)

# These lines set the priors:
priors = {
    'beta_param' : (5.0,2.5,0),
    'log10_beta' : (1.0,0.5,0)} # (continued for all parameters)

UseNeural = True # Use the neural network
# This is a wrapper for MCMC (with initial global optimisation to speed up convergence)
from Chempy.wrapper import single_star_optimization
single_star_optimization() # NB: if not using a neural network, multi_star_optimization() should be used<jupyter_output>['Proto-sun_all']
first minimization for each star separately took:  0 seconds
step  1 of  1000
1.73569619634 1.7353424002
calculation so far took 1.6636462211608887  seconds
step  2 of  1000
1.73569619634 1.73446346283
calculation so far took 1.9948980808258057  seconds
step  3 of  1000
1.73569619634 1.73111670318
calculation so far took 2.3175220489501953  seconds
step  4 of  1000
1.73569619634 1.72968438599
calculation so far took 2.7558858394622803  seconds
step  5 of  1000
1.73569619634 1.72378296961
calculation so far took 3.1382975578308105  seconds
step  6 of  1000
1.73569619634 1.70758419918
calculation so far took 3.51373291015625  seconds
step  7 of  1000
1.73569619634 1.69546069415
calculation so far took 3.8425889015197754  seconds
step  8 of  1000
1.73569619634 1.65553515173
calculation so far took 4.18914532661438  seconds
step  9 of  1000
1.73569619634 1.6015698814
calculation so far took 4.543027877807617  seconds
step  10 of  1000
1.73569619634 1.46826507008
calculati[...]<jupyter_text>Now we restructure the MCMC chain found in the `mcmc/` directory and print the posterior predictions. The parameter names can be taken from the `mcmc/parameter_names.npy` file.<jupyter_code>from Chempy.plot_mcmc import restructure_chain
restructure_chain('mcmc/') # Restructure the MCMC chain

# Print median and 1-sigma values of the posterior PDF
positions = np.load('mcmc/posteriorPDF.npy')
med = np.median(positions,axis=0) # Median
up = np.percentile(positions,100-15.865,axis=0) # Upper 1 sigma bound
low = np.percentile(positions,15.865,axis=0) # Lower 1 sigma bound

print('\nMCMC Parameters \n-----------------')
for i in range(6):
    print('%.2f +%.2f -%.2f' %(med[i],up[i]-med[i],med[i]-low[i]))
# Now create the corner plot (saved as mcmc/parameter_space_sorted.png)
%pylab inline
from Chempy.plot_mcmc import plot_mcmc_chain_with_prior
plot_mcmc_chain_with_prior('mcmc/',use_prior=True,only_first_star=False,plot_true_parameters=False,plot_only_SSP_parameter=False)<jupyter_output>Populating the interactive namespace from numpy and matplotlib<jupyter_text>Step 6: Compute Overall Scores We can now compute both the Bayes and LOO-CV scores for the yield table. Before computation, we must have the `UseNeural` switch set to **on** and have $\beta$ listed in the parameters to optimize as in Step 4.The code for scoring is in the `overall_scores.py` file.Firstly, the Bayes score. This outputs the Bayesian score (in linear space) along with the estimate of its error deriving from the Monte Carlo integration only. Scores are saved in the `OverallScores/` directory as `Bayes_score_[YIELD SET].npz`. Here we use the Chieffi & Limongi 2004 yields scores for illustration.<jupyter_code>from Chempy.overall_scores import overall_Bayes
#overall_Bayes() # This takes around 10 minutes to run on an 8-core machine

# To see the output:
dat=np.load('OverallScores/Bayes_score - chieffi04_net, Karakas_net_yield, Seitenzahl.npz')
print('Log10 Bayes Score: %.3f +/- %.3f' %(np.log10(dat['score']),dat['score_err']/(dat['score']*np.log(10))))<jupyter_output>Log10 Bayes Score: 1.611 +/- 0.002<jupyter_text>We now calculate the CV score. This is computed ten times by default to account for the aforementioned scatter. This takes around 2 hours on an 8-core machine. Scores are saved as `OverallScores/CV_normalised_element_predictions_[YIELD TABLE].npz`<jupyter_code>from Chempy.overall_scores import CV_element_predictions
#CV_element_predictions() # Run on a faster machine

# To see the output:
dat = np.load('OverallScores/CV_normalised_element_predictions_West17_net.npz')

# Load normalised scores (already in log-space)
scores = dat['normalised_scores']

# Print score and error
print('Normalised West & Heger (2017, in prep.) LOO-CV score = %.2f + %.2f - %.2f' %(np.median(scores),np.percentile(scores,100-15.865)-np.median(scores),np.median(scores)-np.percentile(scores,15.865)))<jupyter_output>Normalised West & Heger (2017, in prep.) LOO-CV score = -1.01 + 0.02 - 0.03<jupyter_text>Step 7: Visualising Best Element PredictionsWe can now plot the elemental predictions from the LOO-CV scores (using the element mean and width found from MCMC with that element excluded - see paper for further details).<jupyter_code># Load mean + sigma + element names
pred_mean = np.mean(dat['mean'],axis=0)
pred_sigma = np.mean(dat['sigma'],axis=0) # Taking mean over all 10 iterations
elements = dat['elements']

# Load in proto-solar data in same element order
ps_dat = np.load('Chempy/input/stars/Proto-sun_all.npy')
ps_nam = ps_dat.dtype.names
ps_abun = []; ps_err = []
for el in elements:
    for i in range(len(ps_nam)):
        if ps_nam[i] == el:
            ps_abun.append(ps_dat[0][i])
            ps_err.append(ps_dat[1][i])

# Configure plot
plt.figure(figsize=(20,10))
large_text = 28
text_size = 22
marker_size= 10
small_text = 16
plt.rc('font', family='serif',size = large_text)
plt.rc('xtick', labelsize=small_text)
plt.rc('ytick', labelsize=text_size)
plt.rc('axes', labelsize=text_size, lw=1.0)
plt.rc('lines', linewidth = 2)
plt.rcParams['ytick.major.pad']='8'
plt.rcParams['text.latex.preamble']=[r"\usepackage{libertine}"]
params = {'text.usetex' : True,
    'font.family' : 'libertine',
    'text.latex.unicode': True,
    }
plt.rcParams.update(params)
            
            
# Plot both sets of data
plt.errorbar(np.arange(len(elements)),pred_mean,yerr=pred_sigma,fmt='x',c='k',alpha=0.5,label='LOO-CV Predictions')
plt.errorbar(np.arange(len(elements)),ps_abun,yerr=ps_err,fmt='o',c='b',label='Observations')
plt.legend()
plt.xlabel('Element')
plt.ylabel('[X/Fe] abundance')

# Change x-axis to element names
ax = plt.gca()
la=plt.setp(ax,xticks=np.arange(len(elements)), xticklabels=elements)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11122953202209088, 'fraction_numerical': 0.07387850014533476, 'mean_word_length': 4.2930769230769235, 'pattern_counts': {'":': 0, '<': 58, '<?xml version=': 0, '>': 58, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '42164', 'license': 'MIT', 'n_tokens_mistral': 8069, 'n_tokens_neox': 6867, 'n_words': 2191, 'path': 'Scoring & Parameter Choice Tutorial.ipynb', 'repo_name': 'oliverphilcox/ChempyMulti', 'size': '20671'}