<jupyter_start><jupyter_text>Optimising Model Parameters Using the Evolutionary Optimizers<jupyter_code>from pprint import pprint
import matplotlib.pyplot as plt
import nglview as nv

%matplotlib inline
import isambard
import isambard.optimisation.evo_optimizers as ev_opts
from isambard.optimisation.evo_optimizers import Parameter
def show_ball_and_stick(ampal):
    view = nv.show_text(ampal.pdb)
    view.add_ball_and_stick()
    view.remove_cartoon()
    return view<jupyter_output><empty_output><jupyter_text>Getting Ready Three things are required to run an optimisation:1. A specification1. A list of amino-acid sequences to be packed onto the model1. A list of parametersThe first two are pretty straight forward, chances are you know the specification you want to use and the sequence you want to optimize. In this example I'll use the `CoiledCoil` specification, using the `from_parameters` class method and I'll use the sequence for the [basis-set dimer](http://www.rcsb.org/pdb/explore/explore.do?structureId=4DZM).<jupyter_code>specification = isambard.specifications.CoiledCoil.from_parameters
sequences = ['EIAALKQEIAALKKENAALKWEIAALKQ', 'EIAALKQEIAALKKENAALKWEIAALKQ']<jupyter_output><empty_output><jupyter_text>So far, so simple, but parameters require a bit more effort to set up.  Parameters If we look at the signature of `CoiledCoil.from_parameters` it looks like this:```from_parameters(    n,    aa=28,    major_radius=None,    major_pitch=None,    phi_c_alpha=26.42,    minor_helix_type='alpha',    auto_build=True    )``` It is expecting arguments in this order, so we need to supply parameters to the optimizer that match this format. We can make a list of parameters using the `Parameter` class:<jupyter_code>parameters = [
    Parameter.static('Oligomeric State', 2),
    Parameter.static('Helix Length', 28),
    Parameter.dynamic('Radius', 5.0, 1.0),
    Parameter.dynamic('Pitch', 200, 60),
    Parameter.dynamic('PhiCA', 283, 27),  # 283 is equivalent to a g position
]
parameters<jupyter_output><empty_output><jupyter_text>>  Note> "Parameters" in the `from_parameters` class method name is referring to the geometric parameters that describe a coiled-coil, not the `Parameter` class in the `evo_optimizer` module. `Parameter`s can be created using the `static` or `dynamic` class method. `Parameter`s created with `Parameter.static` will have fixed values during optimisation, whereas those created with `Parameter.dynamic` will be modified by the optimizer.The parameters are created using the following arguments:```pythonParameter.static('Oligomeric State', 2)                ^                   ^___ The static value                |___ a human readable labelParameter.dynamic('Radius', 5.0, 1.0)                 ^         ^    ^___ Value range                 |         |___ Mean value                 |___ a human readable label```Both parameter types have a human-readable label. This is for your own reference really, you can use anything you like, as long as you know what parameter it corresponds to. Static parameters have a single value supplied, while the dynamic value has a mean value and a range. For example, in the radius parameter, we've given a mean value of `5.0`, with a range of `1.0`, this means that the optimizer can create values for radius of between `4.0` and `6.0`.Once you've made your parameter list you can test it to see if it can be used to generate a model correctly:<jupyter_code>default_values = [x.default_value for x in parameters]
print(default_values)<jupyter_output><empty_output><jupyter_text>The `default_value` property returns either the static value or the mean of a set of dynamic values. We can unpack this list in order to supply arguments to the specification.<jupyter_code>test_model = specification(*default_values)
print(test_model)
print(test_model.major_radii)
print(test_model.major_pitches)
print(test_model.phi_c_alphas)<jupyter_output><empty_output><jupyter_text>As you can see, the model was created as expected using the values we supplied for the parameters. Now we're ready to run an optimisation.  Running an Optimisation Once you have all these things in place, you're ready to initialise the optimizer. There are 4 different evolutionary optimizers in the `evo_optimizers` module:1. `GA` (Genetic Algorithm) - Good at eliminating unfavourable regions of parameter space.1. `DE` (Differential Evolution) - A simple and versatile algorithm that quickly finds the minima in smooth energy landscapes.1. `PSO` - (Particle Swarm Optimisation) - A noisy algorithm that explores a lot of the parameter space.1. `CMAES` - (Covariance Matrix Adaptive Evolutionary Strategy) - Thoroughly explores regions of parameter space that are likely to contain the minima. >  Note> It's difficult to know which algorithm is going to be best for your problem, so it's a good idea to try out all of them on a small scale. If you're really unsure about which one to use, the `GA` provides a good balance of parameter-space coverage and thoroughness. Let's take a look at the signature for the `GA`:```python init signatureGA(specification, sequences, parameters, build_fn, eval_fn,   cxpb=0.5, mutpb=0.2, **kwargs)``` Here you can see that the `GA`, and in fact all of the evolutionary optimizers, have 5 required arguments to initialize them. We've already sorted out the first three (specification, sequence and parameters), what about `build_fn` and `eval_fn`? These arguments are functions that the optimizer will use to build a model from the parameters it generates and then to evaluate those models. This means that the optimizers are highly customisable (ask someone if you'd like to know more about this), but to get started, you might want to use the optimizer with predefined build and evaluation functions.<jupyter_code>opt_ga = ev_opts.GA.buff_internal_eval(specification, sequences, parameters)<jupyter_output><empty_output><jupyter_text>The `buff_internal_eval` class method uses `ev_opts.default_build` as the build method and uses the BUFF internal energy as the evaluation function. **This is a good place to start, but this might not be appropriate for all problems.** BUFF is a very soft force field, which is good for parametric modelling as it tolerates small clashes that can occur as a result of using a geometric description with restricted degrees of freedom. For your problem you might need a harder force field, or even some completely different way of evaluating models, in which case you can use another evaluation function or make your own.Now we can start the optimizer using the `run_opt` method. `run_opt` has two required arguments: `pop_size` and `generations`. `pop_size` defines the number of children (i.e. sets of parameters) to be evaluated per generation, while `generations` defines the total number of generations you want to run. We can also supply an additional `cores` argument, which defines the number of CPU cores you wish to use for the optimisation. **Running the next command might take a while! ~2 mins.**<jupyter_code>opt_ga.run_opt(50, 5)<jupyter_output><empty_output><jupyter_text>The optimizer prints various statistics about the optimisation both as it's running and after it finishes. After each generation is complete, a range of information about that generation is printed:* **gen** - the generation number* **evals** - the number of models built and evaluated in that generation. Even though a specific `pop_size` is specified, not all models are built, as models with the same or very similar parameters are not recreated. This is most obvious in the `GA` where parameters with a high "fitness" are retained in the subsequent generation.* **avg** - the average value of the evaluation metric for the generation.* **std** - the standard deviation of the evaluation metric for the generation.* **min** - the lowest (best) score in the generation.* **max** - the highest (worst) score in the generation. >  Note> This logging information can be written to disk using the `log`, `log_path` and `run_id` keyword arguments on `run_opt`. Once the optimisation has finished, the parameters and score for the best model are printed. You can get the AMPAL object for the best model by using the `best_model` property:<jupyter_code>optimized_model = opt_ga.best_model
print(optimized_model)
print(optimized_model.major_radii)
print(optimized_model.major_pitches)
print(optimized_model.phi_c_alphas)<jupyter_output><empty_output><jupyter_text>You can then write the PDB to a file and take a look at your optimized model:<jupyter_code>show_ball_and_stick(optimized_model)<jupyter_output><empty_output><jupyter_text>Of course, this was a very quick optimisation, and to find something closer to the global minima you may need larger values for `pop_size` and `generations`. The std column of output printed by the optimizer gives a good indication of the degree of convergence on a single solution during optimisation, but you should be critical of the optimisation, running it multiple times and evaluating the models produced thoroughly.  Evaluating an Optimisation There are some tools in the optimizer that can give information about the optimisation run that we have performed. The most useful of these is the `make_energy_funnel_data` function, which compares all of the models created during the optimisation to the best model produced.<jupyter_code>funnel_data = opt_ga.make_energy_funnel_data()
pprint(funnel_data[:5])<jupyter_output><empty_output><jupyter_text>`make_energy_funnel_data` returns a list of tuples containing the RMSD between a model and the best model from the optimisation, the BUFF internal energy for that model and the generation that the model came from. We can use these data to make a plot showing gross properties of the optimisation:<jupyter_code>xs = [x[0] for x in funnel_data]
ys = [x[1] for x in funnel_data]
zs = [x[2] for x in funnel_data]
fig, ax = plt.subplots(1, 1)
fig.set_size_inches(4, 4)
ax.grid()
ax.scatter(xs, ys, c=zs)
ax.set_xlabel('RMSD (Å)')
ax.set_ylabel('BUFF Internal Energy')<jupyter_output><empty_output><jupyter_text>In this plot, every model produced by the optimizer is shown as a point, with its BUFF internal energy as the Y axis and the RMSD from the best scoring model as the X axis. The colour represents the generation, purple is the first generation and yellow is the last generation.You can see that there is some sort of structure to the data points, which makes sense; if models have similar RMSDs, they should probably have a similar BUFF internal energy.Let's take a closer look at the best scoring models:<jupyter_code>fig, ax = plt.subplots(1, 1)
fig.set_size_inches(4, 4)
ax.grid()
ax.scatter(xs, ys, c=zs)
ax.set_ylim(top=-700)
ax.set_xlabel('RMSD (Å)')
ax.set_ylabel('BUFF Internal Energy')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.06704220149598301, 'fraction_numerical': 0.006648813371502447, 'mean_word_length': 5.185037121644775, 'pattern_counts': {'":': 0, '<': 49, '<?xml version=': 0, '>': 55, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '157637', 'license': 'MIT', 'n_tokens_mistral': 3036, 'n_tokens_neox': 2895, 'n_words': 1379, 'path': 'tutorial3/003Tutorial.ipynb', 'repo_name': 'CCPBioSim/isambard_workshop', 'size': '10829'}