<jupyter_start><jupyter_text>Task 14 - Parameter Study OptimisationThe previous task sampled from the available parameters and aimed to cover the parameter space efficiently. This task uses Scikit-opt Gaussian Processing to home-in on the optimal solution. Optimisation is even more efficient as it removes the need to sample the entire parameter space and, instead, focuses on the area of interest.In optimisation algorithms it is common to see a combination of exploration and exploitation to find the optimal value.<jupyter_code># Install dependencies

import json
import numpy as np
import pandas as pd
import adaptive
import holoviews
import ipywidgets
import nest_asyncio
import plotly.graph_objects as go

from tqdm import tqdm
from pathlib import Path
from skopt import gp_minimize
from skopt.utils import dump, load
from scipy.interpolate import griddata

from openmc_model import objective

adaptive.notebook_extension()
nest_asyncio.apply()

# method for saving results in json file
def output_result(filepath, result):
    filename = filepath
    Path(filename).parent.mkdir(parents=True, exist_ok=True)
    with open(filename, mode="w", encoding="utf-8") as f:
        json.dump(result, f, indent=4)<jupyter_output><empty_output><jupyter_text>The following codes run neutronics simulations using a simple pre-defined model. We will compare optimised simulations with simple parameter sweeps.  1D OptimisationThe code below runs a simple parameter sweep to obtain TBR as a function of breeder to multiplier ratio in a 1D parameter sweep. These results are the 'true' TBR values across the parameter space that we will compare our optimised results with.<jupyter_code># get_true_values_1D

tbr_values = []
for breeder_percent_in_breeder_plus_multiplier in tqdm(np.linspace(0, 100, 101)):
    tbr_values.append({'breeder_percent_in_breeder_plus_multiplier':breeder_percent_in_breeder_plus_multiplier,
                       'tbr':-objective([breeder_percent_in_breeder_plus_multiplier])})

# results saved in json file
output_result("outputs/1d_tbr_values.json", tbr_values)<jupyter_output><empty_output><jupyter_text>The next code block runs an optimised simulation using the same model, but will search the parameter space for where TBR is maximum. It does this by sampling the parameter space, fitting the results using Gaussian Processing and running a new simulation at the point where TBR is maximum according to the fitted data. If this process is iterated sufficiently, the simulations performed get closer and closer to the point across the parameter space where TBR is maximum.<jupyter_code># get_optimised_values_1d

learner = adaptive.SKOptLearner(objective,
                                dimensions=[(0., 100.)],
                                base_estimator="GP",
                                acq_func="gp_hedge",
                                acq_optimizer="lbfgs",
                               )
runner = adaptive.Runner(learner, ntasks=1, goal=lambda l: l.npoints > 30)

runner.live_info()

runner.ioloop.run_until_complete(runner.task)

# results saved in json file
output_result("outputs/1d_optimised_values.json", dict(learner.data))<jupyter_output><empty_output><jupyter_text>The next code block plots the 'true' simulation data and optimisation data on the same graph. This allows us to see how close the optimisation got to the true maximum TBR across the parameter space.<jupyter_code># 1_plot_1d_optimisation
# Note, optimisation functions tend to minimise the value therefore there are a few negative signs in these scripts

# Loads true data for comparison
data = pd.read_json('outputs/1d_tbr_values.json')

x_data=data['breeder_percent_in_breeder_plus_multiplier']
fx=-data['tbr']


# Load optimisation data

with open('outputs/1d_optimised_values.json', 'r') as f: data = json.load(f).items()

x_vals = [i[0] for i in data]
tbr_vals = [-i[1] for i in data]

# Print max TBR from optimisation data
print('Maximum TBR of ', tbr_vals[-1], 'found with a breeder percent in breeder plus multiplier of ', x_vals[-1])


fig = go.Figure()

# Plot samples from optimsation points
fig.add_trace(go.Scatter(x = x_vals,
                         y = tbr_vals,
                         name="Samples from optimisation",
                         mode='markers',
                         marker=dict(color='red', size=10)
                        )
                )

# Plot true function.
fig.add_trace(go.Scatter(name="True value (unknown)",
                          x = x_data,
                          y = [-i for i in fx],
                          mode='lines',
                          line = {'shape': 'spline'},
                          marker=dict(color='green')
                         )
              )

fig.update_layout(title='Optimal breeder percent in breeder plus multiplier',
                  xaxis={'title': 'breeder percent in breeder plus multiplier', 'range': [0, 100]},
                  yaxis={'title': 'TBR', 'range': [0.1, 2]}
                 )

fig.show()<jupyter_output><empty_output><jupyter_text>As shown, the optimisation samples are spread across the parameter space but are more dense towards the true TBR maximum . This shows how the optimisation homes-in on this point by repeatedly simulating and fitting data.To reach the true maximum TBR value, sufficient simulations must be performed so that the data trend across the parameter space evaluated to a sufficient accuracy. However, optimisation achieved this using fewer samples than the sweep of the entire parameter space as it focused on sampling the important areas of the space (101 sweep samples vs 30 optimised samples).This was a 1D problem, however, the same techniques can be applied to N-dimension problems but the number of simulations required increases. The next example is a 2D dimensional problem where the optimal breeder to multiplier ratio and enrichment are being found.  2D OptimisationThe code below runs a simple parameter sweep to obtain TBR as a function of breeder to multiplier ratio and enrichment in a 2D parameter sweep. These results are the 'true' TBR values across the parameter space that we will compare our optimised results with.<jupyter_code># get_true_values_2D

tbr_values = []
for breeder_percent_in_breeder_plus_multiplier in tqdm(np.linspace(0, 100, 20)):
    for blanket_breeder_li6_enrichment in np.linspace(0, 100, 20):
        tbr_values.append({'breeder_percent_in_breeder_plus_multiplier': breeder_percent_in_breeder_plus_multiplier,
                           'blanket_breeder_li6_enrichment': blanket_breeder_li6_enrichment,
                           'tbr': -objective([breeder_percent_in_breeder_plus_multiplier,
                                              blanket_breeder_li6_enrichment])
                          })

# results saved in json file
output_result("outputs/2d_tbr_values.json", tbr_values)<jupyter_output><empty_output><jupyter_text>The next code block runs an optimised simulation but searches the 2D parameter space for where TBR is maximum.<jupyter_code># get_optimised_values_2d

# Uses adaptive sampling methods from task X to obtain starting points for the optimiser
learner = adaptive.Learner2D(objective, bounds=[(0, 100), (0, 100)])
runner = adaptive.Runner(learner, ntasks=1, goal=lambda l: l.npoints > 40)
runner.live_info()
runner.ioloop.run_until_complete(runner.task)


# Gaussian Processes based optimisation that returns an SciPy optimisation object
res = gp_minimize(objective,          # the function to minimize
                  dimensions=[(0., 100.), (0., 100.)],       # the bounds on each dimension of x
                  n_calls=40,         # the number of evaluations of f
                  n_random_starts=0,  # the number of random initialization points
                  verbose=True,
                  x0=[i for i in list(learner.data.keys())], # initial data from the adaptive sampling method
                  y0=list(learner.data.values()) # initial data from the adaptive sampling method
                  )

# saves 2d optimisation results in .dat file
dump(res, "outputs/2d_optimised_values.dat")<jupyter_output><empty_output><jupyter_text>The next code block plots the true results and optimised results on the same 2D scatter graph.<jupyter_code># 2d_plot_2d_optimisation_scatter.py

# load true data for comparison
data = pd.read_json('outputs/2d_tbr_values.json')
x=data['breeder_percent_in_breeder_plus_multiplier']
y=data['blanket_breeder_li6_enrichment']
z=data['tbr']

# Print max TBR from optimisation data
print('Optimal breeder_percent_in_breeder_plus_multiplier_ratio = ', res.x[0])
print('Optimal Li6 enrichment = ', res.x[1])
print('Maximum TBR = ', -res.fun)

fig = go.Figure()

fig.add_trace(go.Scatter3d(name='TBR values found during optimisation',
                           x=[x[0] for x in res.x_iters],
                           y=[x[1] for x in res.x_iters],
                           z=-res.func_vals,
                           mode='markers',
                           marker=dict(size=7)
                          )
             )

fig.add_trace(go.Scatter3d(name='True values',
                           x=x,
                           y=y,
                           z=z,
                           mode='markers',
                           marker=dict(size=7)
                          )
             )


fig.add_trace(go.Scatter3d(name='Maximum TBR value found',
                           x=[res.x[0]],
                           y=[res.x[1]],
                           z=[-res.fun],
                           mode='markers',
                           marker=dict(size=7)
                          )
             )

fig.update_layout(title='Optimal Li6 enrichment and breeder percent in breeder plus multiplier',
                  scene={'yaxis': {'title': 'Li6 enrichment percent'},
                         'zaxis': {'title': 'breeder percent in breeder plus multiplier'},
                         'zaxis': {'title': 'TBR'}
                        }
                 )

fig.show()<jupyter_output><empty_output><jupyter_text>As shown, the optimisation samples are spread across the parameter space but are more dense towards the true TBR maximum . This shows how the optimisation homes-in on this point by repeatedly simulating and fitting data. In this case, this is a 2D fitting.We can also produce a contour graph to show similar results.<jupyter_code># 2_plot_2d_optimisation_contour

# Print max TBR from optimisation data
print('Optimal Li6 enrichment = ', res.x[0])
print('Optimal breeder percent in breeder plus multiplier = ', res.x[1])
print('Maximum TBR = ', -res.fun)


# creates a grid and interploates values on it
xi = np.linspace(0, 100, 100)
yi = np.linspace(0, 100, 100)
zi = griddata((x, y), z, (xi[None,:], yi[:,None]), method='linear')


fig = go.Figure()

# plots interpolated values as colour map plot
fig.add_trace(trace = go.Contour(
                z=zi,
                x=yi,
                y=xi,
        colorscale="Viridis",
        opacity=0.9,
        line=dict(width=0, smoothing=0.85),
        contours=dict(
            showlines=False,
            showlabels=False,
            size=0,
            labelfont=dict(size=15,),
        ),
    ))

fig.add_trace(go.Scatter(name='TBR values found during optimisation',
                         x=[x[0] for x in res.x_iters],
                         y=[x[1] for x in res.x_iters],
                         hovertext=-res.func_vals,
                         hoverinfo="text",
                         marker={"size": 8},
                         mode='markers'
                        )
             )


# This add the final optimal value found during the optimisation as a seperate scatter point on the graph
fig.add_trace(go.Scatter(name='Maximum TBR value found',
                         x=[res.x[0]],
                         y=[res.x[1]],
                         hovertext=[-res.fun],
                         hoverinfo="text",
                         marker={"size": 8},
                         mode='markers'
                        )
             )

fig.update_layout(title='',
                  xaxis={'title': 'breeder percent in breeder plus multiplier', 'range':(-1, 101)},
                  yaxis={'title': 'blanket breeder li6 enrichment', 'range':(-1, 101)},
                  legend_orientation="h"
                 )

fig.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08270313757039421, 'fraction_numerical': 0.012952534191472244, 'mean_word_length': 2.563933486238532, 'pattern_counts': {'":': 2, '<': 33, '<?xml version=': 0, '>': 35, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, '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': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '491384', 'license': 'MIT', 'n_tokens_mistral': 3580, 'n_tokens_neox': 3395, 'n_words': 1113, 'path': 'tasks/task_14_parameter_study_optimisation/parameter_study_optimisation.ipynb', 'repo_name': 'pshriwise/openmc_workshop', 'size': '12430'}