# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---
# # Tutorial 3: PRISM class properties
# In this tutorial, we use the knowledge gained in the previous two tutorials to explore how the instance/class properties of the `Pipeline`, `Emulator` and `ModelLink` classes can be used for various tasks.
# This includes inspecting the state of the emulator; modifying the pipeline operations and more.
# It is assumed here that the reader has successfully completed the first two tutorials ([Basic usage](1_basic_usage.ipynb) and [ModelLink subclasses](2_modellink_subclasses.ipynb)) and has a basic understanding of Python instance/class properties.
# For this tutorial, we will use the `LineLink` class definition that was introduced in the previous tutorial.
# Therefore, before we can begin, we have to define this class again, initialize the required classes and construct the first iteration:
# +
# Imports
import numpy as np
from prism import Pipeline
from prism.emulator import Emulator
from prism.modellink import ModelLink
# LineLink class definition
class LineLink(ModelLink):
    # Define default model parameters (optional)
    def get_default_model_parameters(self):
        par_dict = {
            'A': [-10, 10, 3],  # Intercept in [-10, 10] with estimate of 3
            'B': [0, 5, 1.5]}   # Slope in [0, 5] with estimate of 1.5
        return(par_dict)
    # Override call_model abstract method
    def call_model(self, emul_i, par_set, data_idx):
        # Calculate the value on a straight line for requested data_idx
        vals = par_set['A']+np.array(data_idx)*par_set['B']
        return(vals)
    # Override get_md_var abstract method
    def get_md_var(self, emul_i, par_set, data_idx):
        # Calculate the model discrepancy variance
        # For a straight line, this value can be set to a constant
        return(1e-4*np.ones_like(data_idx))
# LineLink initialization
model_data = {1: [4.5, 0.1],    # f(1) = 4.5 +- 0.1
              2.5: [6.8, 0.1],  # f(2.5) = 6.8 +- 0.1
              -2: [0, 0.1]}     # f(-2) = 0 +- 0.1
modellink_obj = LineLink(model_data=model_data)
# Pipeline initialization
pipe = Pipeline(modellink_obj, working_dir='prism_line')
# Construction of first iteration
pipe.construct(1)
# -
# ## Pipeline properties
# The `Pipeline` class holds all information that is required to perform all operations in *PRISM* that do not directly modify the emulator (which is contained in the `Emulator` class).
# This includes making projection figures; analyzing an emulator iteration (which does not modify the emulator itself); telling the `Emulator` how to construct an iteration; handling system paths; and many other things.
# As such, we are able to change most of the underlying `Pipeline` properties at any given moment.
# This allows us to modify the operations in the `Pipeline` without making direct changes to the `Emulator` (which could render specific results invalid).
# Note that many parameters of the `Pipeline` can be set during initialization by using the *prism_par* argument (see [PRISM parameters](https://prism-tool.readthedocs.io/en/latest/user/descriptions.html#prism-parameters) for their descriptions).
# When a new emulator is being constructed, the `Pipeline` object will generate a Latin-Hypercube design of `n_sam_init` samples, which by default is set to:
pipe.n_sam_init
# These samples are then evaluated in the model and used to construct an emulator.
# However, for a model like the one in our `LineLink`, this is probably a bit overkill.
# So, we could set it to, let's say, $100$ and reconstruct the first iteration of the emulator:
pipe.n_sam_init = 100
pipe.construct(1, force=True)
# Here, we reconstructed the first iteration of the emulator by first setting `n_sam_init` to $100$ and then calling the `construct()`-method with *force=True* (which overrides the default behavior of skipping construction if already finished).
# Although the part of parameter space that is still remaining is fairly similar as before, we can see that the evaluation rate of the emulator has significantly increased.
# This could potentially be very beneficial to us if we were to evaluate the emulator many times in later iterations.
#
# However, we can also see that the `Pipeline` is warning us that there are probably not enough plausible samples to construct a more accurate emulator iteration.
# The reason for this is because there is such a small part of parameter space still remaining as plausible, that we have to evaluate much more samples to obtain a decent number of plausible samples.
# The number of samples that are evaluated in the emulator during an analysis is influenced by the iteration number `emul_i`, the number of model parameters `n_par` and the base evaluation number `base_eval_sam`.
# While we obviously cannot change the iteration number or the number of model parameters, we can change the base evaluation number.
# At this moment, the base evaluation number is set to:
pipe.base_eval_sam
# As we used $100$ samples for the first iteration and we want to have at least this many plausible samples for the next iteration, let's set `base_eval_sam` to $16,000$ and reanalyze the iteration:
pipe.base_eval_sam = 16000
pipe.analyze()
# And now we have enough plausible samples for the next iteration.
# We can also see that due to evaluating much more samples, that the percent of parameter space remaining is more accurate.
# If we wanted to, we could check how many samples have been evaluated, how many are plausible and what they are:
print("Number of evaluated samples in iteration 1: %i" % (pipe.n_eval_sam[1]))
print("Number of plausible samples in iteration 1: %i" % (pipe.n_impl_sam[1]))
print("Plausible samples: %s" % (pipe.impl_sam))
# The two numbers in the output show us the same thing as the `details()` overview, but we can also see that all the plausible samples are quite clustered.
# This would mean that the next emulator iteration is going to be much more accurate than the current one.
# Note that another way to increase the number of plausible samples in an iteration, is by changing the implausibility parameters, as shown in the first tutorial.
# A few other useful properties are:
print("Path to emulator HDF5-file: %r" % (pipe.hdf5_file))
print("Is file-logging enabled? %s" % (pipe.do_logging))
print("Are parameters split into active and passive parameters? %s" % (pipe.do_active_anal))
print("Bound ModelLink object: %s" % (pipe.modellink))
print("Current implausibility parameters: %s" % (pipe.impl_cut[1]))
print("Current # of implausibility wildcards: %i" % (pipe.cut_idx[1]))
# As a general rule when using *PRISM*: All class attributes that hold data specific to an emulator iteration, start with index 1 instead of index 0, as we have shown above for some class properties.
# This was done to make it more intuitive to retrieve iteration specific data, as you can simply use the iteration number as the index.
# ## Emulator properties
# Whereas it can be useful to look at the `Pipeline` properties, it does not tell us anything about the state of the emulator itself.
# Although they cannot be modified, looking at the `Emulator` properties can give us a lot of information.
# As the underlying algorithms of *PRISM* involve a lot of math, they can however be a bit harder to understand.
# Keep in mind that only some `Emulator` properties are converted to a (more) human-readable representation from their underlying counterparts.
# All `Emulator` properties can be accessed with `pipe.emulator.xxx`.
# Probably the most interesting and useful property to look at, is the `poly_terms` property (which combines several other properties together and solely exists for the user):
pipe.emulator.poly_terms[1]
# This dictionary shows us what the polynomial terms are for calculating every data point in the first emulator iteration.
# So, if we remember that our model is given by the function $f(x) = A+Bx$, then we can immediately see that the polynomial terms are very close to what the real terms should look like, where $A$ is multiplied by unity and $B$ is multiplied by whatever the value of $x$ is for that data point (with $x$ being the key in the outer dictionary).
# We can also check the performance of the emulator iteration by looking at the amount of variance in the model outputs that the emulator could not explain with the polynomial functions.
# This variance is called the *residual variance*, and a general rule is that the lower it is, the more accurate is the emulator:
pipe.emulator.rsdl_var[1]
# These values tell us basically the same thing as the `poly_terms` property, just in a different way: the emulator is capturing almost all model behavior.
# The difference between these two properties however, is that we usually do not know what the polynomial terms should look like, while we do know that a low residual variance is a good thing.
# Therefore, it is usually best to view both at the same time, where the `rsdl_var` tells us how much we should believe that the `poly_terms` are accurate.
# There are also a few other `Emulator` properties we could look at.
# For example, we might be interested in checking which parameters are considered 'active'.
# Although the `details()` overview tells us which parameters are considered active at least once over all data points combined (this is indicated by the asterisk in front of the parameter name), we would like to know the active parameters for each data point individually.
# This can be done in the following way:
pipe.emulator.active_par_data[1]
# Of course, we already knew that both parameters are active for all data points, as otherwise their regression functions would not contain any polynomial terms that contains either parameter.
# However, when emulating models that have a large number of model parameters, checking the `active_par_data` property for this information is much easier than the `poly_terms` property.
# (In case you would like to see this for yourself, go back to the emulator made in the first tutorial with the `GaussianLink` class and take a look at its `poly_terms`, `rsdl_var` and `active_par_data` properties.)