MOBO#
# run only once during the notebook execution
!git clone https://github.com/cfteach/modules.git &> /dev/null
!pip install ax-platform &> /dev/null
!pip install ipyvolume &> /dev/null
!pip install plotly
Requirement already satisfied: plotly in /home/ksuresh/miniconda3/envs/env_AID2E/lib/python3.12/site-packages (5.22.0)
Requirement already satisfied: tenacity>=6.2.0 in /home/ksuresh/miniconda3/envs/env_AID2E/lib/python3.12/site-packages (from plotly) (8.3.0)
Requirement already satisfied: packaging in /home/ksuresh/miniconda3/envs/env_AID2E/lib/python3.12/site-packages (from plotly) (24.0)
import torch
torch.cuda.is_available()
True
%load_ext autoreload
%autoreload 2
from IPython.display import display, Math, Latex
import os
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
#import AI4NP_detector_opt.sol2.detector2 as detector2
import modules.detector2 as detector2
import re
import pickle
#import dill
import torch
from ax.metrics.noisy_function import GenericNoisyFunctionMetric
from ax.service.utils.report_utils import exp_to_df #https://ax.dev/api/service.html#ax.service.utils.report_utils.exp_to_df
from ax.runners.synthetic import SyntheticRunner
# Plotting imports and initialization
from ax.utils.notebook.plotting import render, init_notebook_plotting
from ax.plot.contour import plot_contour
from ax.plot.pareto_utils import compute_posterior_pareto_frontier
from ax.plot.pareto_frontier import plot_pareto_frontier
#init_notebook_plotting()
# Model registry for creating multi-objective optimization models.
from ax.modelbridge.registry import Models
# Analysis utilities, including a method to evaluate hypervolumes
from ax.modelbridge.modelbridge_utils import observed_hypervolume
from ax import SumConstraint
from ax import OrderConstraint
from ax import ParameterConstraint
from ax.core.search_space import SearchSpace
from ax.core.parameter import RangeParameter,ParameterType
from ax.core.objective import MultiObjective, Objective, ScalarizedObjective
from ax.core.optimization_config import ObjectiveThreshold, MultiObjectiveOptimizationConfig
from ax.core.experiment import Experiment
from botorch.utils.multi_objective.box_decompositions.dominated import DominatedPartitioning
from ax.core.data import Data
from ax.core.types import ComparisonOp
from sklearn.utils import shuffle
from functools import wraps
Create detector geometry and simulate tracks#
The module detector creates a simple 2D geometry of a wire based tracker made by 4 planes.
The adjustable parameters are the radius of each wire, the pitch (along the y axis), and the shift along y and z of a plane with respect to the previous one.
A total of 8 parameters can be tuned.
The goal of this toy model, is to tune the detector design so to optimize the efficiency (fraction of tracks which are detected) as well as the cost for its realization. As a proxy for the cost, we use the material/volume (the surface in 2D) of the detector. For a track to be detetected, in the efficiency definition we require at least two wires hit by the track.
So we want to maximize the efficiency (defined in detector.py) and minimize the cost.
LIST OF PARAMETERS#
(baseline values)
R = .5 [cm]
pitch = 4.0 [cm]
y1 = 0.0, y2 = 0.0, y3 = 0.0, z1 = 2.0, z2 = 4.0, z3 = 6.0 [cm]
# CONSTANT PARAMETERS
#------ define mother region ------#
y_min=-10.1
y_max=10.1
N_tracks = 1000
print("::::: BASELINE PARAMETERS :::::")
R = .5
pitch = 4.0
y1 = 0.0
y2 = 0.0
y3 = 0.0
z1 = 2.0
z2 = 4.0
z3 = 6.0
print("R, pitch, y1, y2, y3, z1, z2, z3: ", R, pitch, y1, y2, y3, z1, z2, z3,"\n")
#------------- GEOMETRY ---------------#
print(":::: INITIAL GEOMETRY ::::")
tr = detector2.Tracker(R, pitch, y1, y2, y3, z1, z2, z3)
Z, Y = tr.create_geometry()
num_wires = detector2.calculate_wires(Y, y_min, y_max)
volume = detector2.wires_volume(Y, y_min, y_max,R)
detector2.geometry_display(Z, Y, R, y_min=y_min, y_max=y_max,block=False,pause=5) #5
print("# of wires: ", num_wires, ", volume: ", volume)
#------------- TRACK GENERATION -----------#
print(":::: TRACK GENERATION ::::")
t = detector2.Tracks(b_min=y_min, b_max=y_max, alpha_mean=0, alpha_std=0.2)
tracks = t.generate(N_tracks)
detector2.geometry_display(Z, Y, R, y_min=y_min, y_max=y_max,block=False, pause=-1)
detector2.tracks_display(tracks, Z,block=False,pause=-1)
#a track is detected if at least two wires have been hit
score = detector2.get_score(Z, Y, tracks, R)
frac_detected = score[0]
resolution = score[1]
print("fraction of tracks detected: ",frac_detected)
print("resolution: ",resolution)
::::: BASELINE PARAMETERS :::::
R, pitch, y1, y2, y3, z1, z2, z3: 0.5 4.0 0.0 0.0 0.0 2.0 4.0 6.0
:::: INITIAL GEOMETRY ::::
# of wires: 20 , volume: 62.800000000000004
:::: TRACK GENERATION ::::
fraction of tracks detected: 0.264
resolution: 0.24613882479204957
Define Objectives#
Defines a class for the objectives of the problem that can be used in the MOO.
class objectives():
def __init__(self,tracks,y_min,y_max):
self.tracks = tracks
self.y_min = y_min
self.y_max = y_max
def wrapper_geometry(fun):
def inner(self):
R, pitch, y1, y2, y3, z1, z2, z3 = self.X
self.geometry(R, pitch, y1, y2, y3, z1, z2, z3)
return fun(self)
return inner
def update_tracks(self, new_tracks):
self.tracks = new_tracks
def update_design_point(self,X):
self.X = X
def geometry(self,R, pitch, y1, y2, y3, z1, z2, z3):
tr = detector2.Tracker(R, pitch, y1, y2, y3, z1, z2, z3)
self.R = R
self.Z, self.Y = tr.create_geometry()
@wrapper_geometry
def calc_score(self):
res = detector2.get_score(self.Z, self.Y, self.tracks, self.R)
assert res[0] >= 0 and res[1] >= 0,"Fraction or Resolution negative."
return res
def get_score(self,X):
R, pitch, y1, y2, y3, z1, z2, z3 = X
self.geometry(R, pitch, y1, y2, y3, z1, z2, z3)
res = detector2.get_score(self.Z, self.Y, self.tracks, self.R)
return res
def get_volume(self):
volume = detector2.wires_volume(self.Y, self.y_min, self.y_max,self.R)
return volume
res = objectives(tracks,y_min,y_max)
#res.geometry(R, pitch, y1, y2, y3, z1, z2, z3)
X = R, pitch, y1, y2, y3, z1, z2, z3
#fscore = res.get_score(X)
res.update_design_point(X)
fscore = res.calc_score()[0]
fvolume = res.get_volume()
print("...check: ", fvolume, fscore)
...check: 62.800000000000004 0.264
Multi-Objective Optimization#
We will be using ax-platform (https://ax.dev).
In this example we will be using Multi-Objective Bayesian Optimization (MOBO) using qNEHVI + SAASBO
Notice that every function is minimized. Our efficiency is defined as an tracking inefficiency = 1 - efficiency
We add the resolution as a third objective. The average residual of the track hit from the wire centre is used as a proxy for the resolution for this toy-model
#---------------------- BOTORCH FUNCTIONS ------------------------#
def build_experiment(search_space,optimization_config):
experiment = Experiment(
name="pareto_experiment",
search_space=search_space,
optimization_config=optimization_config,
runner=SyntheticRunner(),
)
return experiment
def glob_fun(loc_fun):
@wraps(loc_fun)
def inner(xdic):
x_sorted = [xdic[p_name] for p_name in xdic.keys()] #it assumes x will be given as, e.g., dictionary
res = list(loc_fun(x_sorted))
return res
return inner
def initialize_experiment(experiment,N_INIT):
sobol = Models.SOBOL(search_space=experiment.search_space)
experiment.new_batch_trial(sobol.gen(N_INIT)).run()
return experiment.fetch_data()
@glob_fun
def ftot(xdic):
return (1- res.get_score(xdic)[0], res.get_volume(), res.get_score(xdic)[1])
def f1(xdic):
return ftot(xdic)[0] #obj1
def f2(xdic):
return ftot(xdic)[1] #obj2
def f3(xdic):
return ftot(xdic)[2] #obj3
tkwargs = {
"dtype": torch.double,
"device": torch.device("cuda" if torch.cuda.is_available() else "cpu"),
}
# Define Hyper-parameters for the optimization
N_BATCH = 5
BATCH_SIZE = 5
dim_space = 8 # len(X)
N_INIT = 2 * (dim_space + 1) #
lowerv = np.array([0.5,2.5,0.,0.,0.,2.,2.,2.])
upperv = np.array([1.0,5.0,4.,4.,4.,10.,10.,10.])
from ax import (
Data,
Experiment,
Metric,
Objective,
OptimizationConfig,
ParameterType,
RangeParameter,
Runner,
SearchSpace,
)
# defining the search space one can also include constraints in this function
search_space = SearchSpace(
parameters=
[RangeParameter(name=f"x{i}", lower=lowerv[i], upper=upperv[i],
parameter_type=ParameterType.FLOAT) for i in range(dim_space)]
)
print (search_space)
# define the metrics for optimization
metric_a = GenericNoisyFunctionMetric("a", f=f1, noise_sd=0.0, lower_is_better=True)
metric_b = GenericNoisyFunctionMetric("b", f=f2, noise_sd=0.0, lower_is_better=True)
metric_c = GenericNoisyFunctionMetric("c", f=f3, noise_sd=0.0, lower_is_better=True)
mo = MultiObjective(objectives=[Objective(metric=metric_a),
Objective(metric=metric_b),
Objective(metric=metric_c)
]
)
ref_point = [-1.1]*len(mo.metrics)
refpoints = torch.Tensor(ref_point).to(**tkwargs) # [1.1, 1.1, 1.1] for 3 objs
objective_thresholds = [ObjectiveThreshold(metric=metric, bound=val, relative=False, op=ComparisonOp.LEQ)
for metric, val in zip(mo.metrics, refpoints) #---> this requires defining a torch.float64 object --- by default is (-)1.1 for DTLZ
]
optimization_config = MultiObjectiveOptimizationConfig(
objective=mo,
objective_thresholds=objective_thresholds
)
SearchSpace(parameters=[RangeParameter(name='x0', parameter_type=FLOAT, range=[0.5, 1.0]), RangeParameter(name='x1', parameter_type=FLOAT, range=[2.5, 5.0]), RangeParameter(name='x2', parameter_type=FLOAT, range=[0.0, 4.0]), RangeParameter(name='x3', parameter_type=FLOAT, range=[0.0, 4.0]), RangeParameter(name='x4', parameter_type=FLOAT, range=[0.0, 4.0]), RangeParameter(name='x5', parameter_type=FLOAT, range=[2.0, 10.0]), RangeParameter(name='x6', parameter_type=FLOAT, range=[2.0, 10.0]), RangeParameter(name='x7', parameter_type=FLOAT, range=[2.0, 10.0])], parameter_constraints=[])
SearchSpace(parameters=[RangeParameter(name='x0', parameter_type=ParameterType.FLOAT, lower=0.5, upper= 1.0), RangeParameter(name='x1', parameter_type=ParameterType.FLOAT, lower=2.5, upper=5.0), RangeParameter(name='x2', parameter_type=ParameterType.FLOAT, lower=0.0, upper=4.0), RangeParameter(name='x3', parameter_type=ParameterType.FLOAT, lower=0.0, upper=4.0), RangeParameter(name='x4', parameter_type=ParameterType.FLOAT, lower=0.0, upper= 4.0), RangeParameter(name='x5', parameter_type=ParameterType.FLOAT, lower = 2.0, upper = 10.0), RangeParameter(name='x6', parameter_type=ParameterType.FLOAT, lower =2.0, upper = 10.0), RangeParameter(name='x7', parameter_type=ParameterType.FLOAT, lower = 2.0, upper = 10.0)], parameter_constraints=[])
SearchSpace(parameters=[RangeParameter(name='x0', parameter_type=FLOAT, range=[0.5, 1.0]), RangeParameter(name='x1', parameter_type=FLOAT, range=[2.5, 5.0]), RangeParameter(name='x2', parameter_type=FLOAT, range=[0.0, 4.0]), RangeParameter(name='x3', parameter_type=FLOAT, range=[0.0, 4.0]), RangeParameter(name='x4', parameter_type=FLOAT, range=[0.0, 4.0]), RangeParameter(name='x5', parameter_type=FLOAT, range=[2.0, 10.0]), RangeParameter(name='x6', parameter_type=FLOAT, range=[2.0, 10.0]), RangeParameter(name='x7', parameter_type=FLOAT, range=[2.0, 10.0])], parameter_constraints=[])
experiment = build_experiment(search_space,optimization_config)
data = initialize_experiment(experiment,N_INIT)
hv_list = []
for i in range(N_BATCH):
print("\n\n...PROCESSING BATCH n.: {}\n\n".format(i+1))
model = Models.FULLYBAYESIANMOO(
experiment=experiment,
data=data, # tell the data
# use fewer num_samples and warmup_steps to speed up this tutorial
num_samples=32,#256
warmup_steps=64,#512
torch_device=tkwargs["device"],
verbose=False, # Set to True to print stats from MCMC
disable_progbar=False, # Set to False to print a progress bar from MCMC
)
generator_run = model.gen(BATCH_SIZE) #ask BATCH_SIZE points
trial = experiment.new_batch_trial(generator_run=generator_run)
trial.run()
data = Data.from_multiple_data([data, trial.fetch_data()]) #https://ax.dev/api/core.html#ax.Data.from_multiple_data
print("\n\n\n...calculate df via exp_to_df (i.e., global dataframe so far):\n\n")
metric_names = {index: i for index, i in enumerate(mo.metric_names)}
N_METRICS = len(metric_names)
df = exp_to_df(experiment).sort_values(by=["trial_index"])
outcomes = torch.tensor(df[mo.metric_names].values)
#outcomes, _ = data_to_outcomes(data, N_INIT, i+1, BATCH_SIZE, N_METRICS, metric_names)
partitioning = DominatedPartitioning(ref_point=refpoints, Y=outcomes.to(**tkwargs))
try:
hv = partitioning.compute_hypervolume().item()
except:
hv = 0
print("Failed to compute hv")
hv_list.append(hv)
print(f"Iteration: {i+1}, HV: {hv}")
...PROCESSING BATCH n.: 1
Sample: 100%|██████████| 96/96 [00:22, 4.24it/s, step size=6.02e-01, acc. prob=0.813]
Sample: 100%|██████████| 96/96 [00:27, 3.52it/s, step size=3.24e-01, acc. prob=0.932]
Sample: 100%|██████████| 96/96 [00:24, 3.87it/s, step size=4.99e-01, acc. prob=0.887]
[WARNING 05-31 11:53:19] ax.service.utils.report_utils: Column reason missing for all trials. Not appending column.
...calculate df via exp_to_df (i.e., global dataframe so far):
Iteration: 1, HV: 804.0200360632218
...PROCESSING BATCH n.: 2
Sample: 100%|██████████| 96/96 [00:20, 4.76it/s, step size=7.79e-01, acc. prob=0.756]
Sample: 100%|██████████| 96/96 [00:26, 3.63it/s, step size=3.59e-01, acc. prob=0.900]
Sample: 100%|██████████| 96/96 [00:20, 4.74it/s, step size=6.42e-01, acc. prob=0.740]
/home/ksuresh/miniconda3/envs/env_AID2E/lib/python3.12/site-packages/linear_operator/utils/cholesky.py:40: NumericalWarning:
A not p.d., added jitter of 1.0e-08 to the diagonal
/home/ksuresh/miniconda3/envs/env_AID2E/lib/python3.12/site-packages/linear_operator/utils/cholesky.py:40: NumericalWarning:
A not p.d., added jitter of 1.0e-08 to the diagonal
[WARNING 05-31 11:54:45] ax.service.utils.report_utils: Column reason missing for all trials. Not appending column.
...calculate df via exp_to_df (i.e., global dataframe so far):
Iteration: 2, HV: 804.023028474837
...PROCESSING BATCH n.: 3
Sample: 100%|██████████| 96/96 [00:24, 3.98it/s, step size=6.14e-01, acc. prob=0.848]
Sample: 100%|██████████| 96/96 [00:30, 3.12it/s, step size=4.38e-01, acc. prob=0.935]
Sample: 100%|██████████| 96/96 [00:25, 3.78it/s, step size=6.12e-01, acc. prob=0.652]
/home/ksuresh/miniconda3/envs/env_AID2E/lib/python3.12/site-packages/linear_operator/utils/cholesky.py:40: NumericalWarning:
A not p.d., added jitter of 1.0e-08 to the diagonal
/home/ksuresh/miniconda3/envs/env_AID2E/lib/python3.12/site-packages/botorch/optim/optimize.py:251: RuntimeWarning:
Optimization failed in `gen_candidates_scipy` with the following warning(s):
[NumericalWarning('A not p.d., added jitter of 1.0e-08 to the diagonal'), NumericalWarning('A not p.d., added jitter of 1.0e-07 to the diagonal'), NumericalWarning('A not p.d., added jitter of 1.0e-06 to the diagonal'), NumericalWarning('A not p.d., added jitter of 1.0e-05 to the diagonal'), NumericalWarning('A not p.d., added jitter of 1.0e-04 to the diagonal'), NumericalWarning('A not p.d., added jitter of 1.0e-03 to the diagonal'), NumericalWarning('A not p.d., added jitter of 1.0e-08 to the diagonal'), NumericalWarning('A not p.d., added jitter of 1.0e-07 to the diagonal'), NumericalWarning('A not p.d., added jitter of 1.0e-06 to the diagonal'), NumericalWarning('A not p.d., added jitter of 1.0e-05 to the diagonal'), NumericalWarning('A not p.d., added jitter of 1.0e-04 to the diagonal'), NumericalWarning('A not p.d., added jitter of 1.0e-03 to the diagonal'), BotorchWarning('Low-rank cholesky updates failed due NaNs or due to an ill-conditioned covariance matrix. Falling back to standard sampling.'), OptimizationWarning('Optimization failed within `scipy.optimize.minimize` with status 2 and message ABNORMAL_TERMINATION_IN_LNSRCH.'), NumericalWarning('A not p.d., added jitter of 1.0e-08 to the diagonal'), NumericalWarning('A not p.d., added jitter of 1.0e-07 to the diagonal'), NumericalWarning('A not p.d., added jitter of 1.0e-06 to the diagonal'), NumericalWarning('A not p.d., added jitter of 1.0e-05 to the diagonal'), NumericalWarning('A not p.d., added jitter of 1.0e-04 to the diagonal'), NumericalWarning('A not p.d., added jitter of 1.0e-03 to the diagonal'), BotorchWarning('Low-rank cholesky updates failed due NaNs or due to an ill-conditioned covariance matrix. Falling back to standard sampling.'), NumericalWarning('A not p.d., added jitter of 1.0e-08 to the diagonal'), NumericalWarning('A not p.d., added jitter of 1.0e-07 to the diagonal'), NumericalWarning('A not p.d., added jitter of 1.0e-06 to the diagonal'), NumericalWarning('A not p.d., added jitter of 1.0e-05 to the diagonal'), NumericalWarning('A not p.d., added jitter of 1.0e-04 to the diagonal'), NumericalWarning('A not p.d., added jitter of 1.0e-03 to the diagonal'), BotorchWarning('Low-rank cholesky updates failed due NaNs or due to an ill-conditioned covariance matrix. Falling back to standard sampling.'), OptimizationWarning('Optimization failed within `scipy.optimize.minimize` with status 2 and message ABNORMAL_TERMINATION_IN_LNSRCH.'), NumericalWarning('A not p.d., added jitter of 1.0e-08 to the diagonal'), NumericalWarning('A not p.d., added jitter of 1.0e-07 to the diagonal'), NumericalWarning('A not p.d., added jitter of 1.0e-06 to the diagonal'), NumericalWarning('A not p.d., added jitter of 1.0e-05 to the diagonal'), NumericalWarning('A not p.d., added jitter of 1.0e-04 to the diagonal'), NumericalWarning('A not p.d., added jitter of 1.0e-03 to the diagonal'), BotorchWarning('Low-rank cholesky updates failed due NaNs or due to an ill-conditioned covariance matrix. Falling back to standard sampling.'), NumericalWarning('A not p.d., added jitter of 1.0e-08 to the diagonal'), NumericalWarning('A not p.d., added jitter of 1.0e-07 to the diagonal'), NumericalWarning('A not p.d., added jitter of 1.0e-06 to the diagonal'), NumericalWarning('A not p.d., added jitter of 1.0e-05 to the diagonal'), NumericalWarning('A not p.d., added jitter of 1.0e-04 to the diagonal'), NumericalWarning('A not p.d., added jitter of 1.0e-03 to the diagonal')]
Trying again with a new set of initial conditions.
/home/ksuresh/miniconda3/envs/env_AID2E/lib/python3.12/site-packages/botorch/optim/optimize.py:251: RuntimeWarning:
Optimization failed on the second try, after generating a new set of initial conditions.
/home/ksuresh/miniconda3/envs/env_AID2E/lib/python3.12/site-packages/linear_operator/utils/cholesky.py:40: NumericalWarning:
A not p.d., added jitter of 1.0e-08 to the diagonal
/home/ksuresh/miniconda3/envs/env_AID2E/lib/python3.12/site-packages/linear_operator/utils/cholesky.py:40: NumericalWarning:
A not p.d., added jitter of 1.0e-07 to the diagonal
/home/ksuresh/miniconda3/envs/env_AID2E/lib/python3.12/site-packages/linear_operator/utils/cholesky.py:40: NumericalWarning:
A not p.d., added jitter of 1.0e-06 to the diagonal
/home/ksuresh/miniconda3/envs/env_AID2E/lib/python3.12/site-packages/linear_operator/utils/cholesky.py:40: NumericalWarning:
A not p.d., added jitter of 1.0e-05 to the diagonal
/home/ksuresh/miniconda3/envs/env_AID2E/lib/python3.12/site-packages/linear_operator/utils/cholesky.py:40: NumericalWarning:
A not p.d., added jitter of 1.0e-04 to the diagonal
/home/ksuresh/miniconda3/envs/env_AID2E/lib/python3.12/site-packages/linear_operator/utils/cholesky.py:40: NumericalWarning:
A not p.d., added jitter of 1.0e-03 to the diagonal
/home/ksuresh/miniconda3/envs/env_AID2E/lib/python3.12/site-packages/botorch/acquisition/cached_cholesky.py:150: BotorchWarning:
Low-rank cholesky updates failed due NaNs or due to an ill-conditioned covariance matrix. Falling back to standard sampling.
[WARNING 05-31 11:57:15] ax.service.utils.report_utils: Column reason missing for all trials. Not appending column.
...calculate df via exp_to_df (i.e., global dataframe so far):
Iteration: 3, HV: 804.023028474837
...PROCESSING BATCH n.: 4
Sample: 100%|██████████| 96/96 [00:26, 3.58it/s, step size=3.34e-01, acc. prob=0.922]
Sample: 100%|██████████| 96/96 [00:31, 3.01it/s, step size=7.06e-01, acc. prob=0.808]
Sample: 100%|██████████| 96/96 [00:28, 3.36it/s, step size=7.25e-01, acc. prob=0.784]
/home/ksuresh/miniconda3/envs/env_AID2E/lib/python3.12/site-packages/linear_operator/utils/cholesky.py:40: NumericalWarning:
A not p.d., added jitter of 1.0e-08 to the diagonal
/home/ksuresh/miniconda3/envs/env_AID2E/lib/python3.12/site-packages/linear_operator/utils/cholesky.py:40: NumericalWarning:
A not p.d., added jitter of 1.0e-08 to the diagonal
/home/ksuresh/miniconda3/envs/env_AID2E/lib/python3.12/site-packages/linear_operator/utils/cholesky.py:40: NumericalWarning:
A not p.d., added jitter of 1.0e-08 to the diagonal
/home/ksuresh/miniconda3/envs/env_AID2E/lib/python3.12/site-packages/linear_operator/utils/cholesky.py:40: NumericalWarning:
A not p.d., added jitter of 1.0e-08 to the diagonal
/home/ksuresh/miniconda3/envs/env_AID2E/lib/python3.12/site-packages/linear_operator/utils/cholesky.py:40: NumericalWarning:
A not p.d., added jitter of 1.0e-08 to the diagonal
[WARNING 05-31 11:59:26] ax.modelbridge.base: TorchModelBridge(model=FullyBayesianMOOBotorchModel) was not able to generate 5 unique candidates. Generated arms have the following weights, as there are repeats:
[3.0, 1.0, 1.0]
[WARNING 05-31 11:59:26] ax.service.utils.report_utils: Column reason missing for all trials. Not appending column.
...calculate df via exp_to_df (i.e., global dataframe so far):
Iteration: 4, HV: 804.5932903331675
...PROCESSING BATCH n.: 5
Sample: 100%|██████████| 96/96 [00:28, 3.40it/s, step size=5.58e-01, acc. prob=0.864]
Sample: 100%|██████████| 96/96 [00:30, 3.11it/s, step size=5.71e-01, acc. prob=0.892]
Sample: 100%|██████████| 96/96 [00:49, 1.94it/s, step size=1.15e-01, acc. prob=0.970]
/home/ksuresh/miniconda3/envs/env_AID2E/lib/python3.12/site-packages/linear_operator/utils/cholesky.py:40: NumericalWarning:
A not p.d., added jitter of 1.0e-08 to the diagonal
[WARNING 05-31 12:01:59] ax.service.utils.report_utils: Column reason missing for all trials. Not appending column.
...calculate df via exp_to_df (i.e., global dataframe so far):
Iteration: 5, HV: 805.555527745325
Analysis of Results#
Inspecting the Hyper volume statistics#
import plotly.express as px
fig = px.scatter(x = np.arange(N_BATCH) + 1, y = hv_list,
labels={"x": "N_BATCHES",
"y": "Hyper Volume"},
width = 800, height = 800,
title = "HyperVolume Improvement", )
fig.update_traces(marker=dict(size=8,
line=dict(width=2,
color='DarkSlateGrey')),
selector=dict(mode='marker+line'))
fig.data[0].update(mode = "markers+lines")
fig.show()
import matplotlib.pyplot as plt
plt.figure(figsize = (10, 7))
plt.plot(np.arange(N_BATCH) + 1 , hv_list, "ro-")
plt.xlabel("N_BATCHES", fontsize = 12)
plt.ylabel("Hyper Volume", fontsize = 12)
plt.show()
Overall Performance in the Objective space.#
fig1 = px.scatter_3d(df, x="a", y="b", z = "c", color = "trial_index",
labels = { "a": "InEfficiency",
"b": "Volume",
"c": "Resolution"
}, hover_data = df.columns,
height = 800, width = 800)
fig1.show()
Exploration as a function of Iteration number#
obj_fig = px.scatter_3d(df, x="a", y="b", z = "c", animation_frame="trial_index", color="trial_index",
range_x=[0., 0.6], range_y=[0. , 400.], range_z=[0., 0.6],
labels = { "a": "InEfficiency",
"b": "Volume",
"c": "Resolution"}, hover_data = df.columns,
width = 800, height = 800)
obj_fig.update(layout_coloraxis_showscale=False)
obj_fig.show()
Computing posterior pareto frontiers.#
Once can sample expected approximate pareto front solution from the built surrogate model.
from ax.core import metric
# https://ax.dev/api/plot.html#ax.plot.pareto_utils.compute_posterior_pareto_frontier
# absolute_metrics – List of outcome metrics that should NOT be relativized w.r.t. the status quo
# (all other outcomes will be in % relative to status_quo).
# Note that approximated pareto frontier is can be visualized only against 2 objectives.
# So one can try to make mixed plots, to see the ``
n_points_surrogate = 25
frontier = [] #(a,b), (a,c), (b,c)
metric_combos = [(metric_a, metric_b), (metric_a, metric_c), (metric_b, metric_c)]
for combo in metric_combos:
print ("computing pareto frontier : ", combo)
frontier.append(compute_posterior_pareto_frontier(
experiment=experiment,
data=experiment.fetch_data(),
primary_objective=combo[0], #_b
secondary_objective=combo[1], #_a
absolute_metrics=["a", "b", "c"],
num_points=n_points_surrogate,
))
#render(plot_pareto_frontier(frontier, CI_level=0.9))
#res_front = plot_pareto_frontier(frontier, CI_level=0.8)
computing pareto frontier : (GenericNoisyFunctionMetric('a'), GenericNoisyFunctionMetric('b'))
/home/ksuresh/miniconda3/envs/env_AID2E/lib/python3.12/site-packages/botorch/optim/fit.py:102: OptimizationWarning:
`scipy_minimize` terminated with status OptimizationStatus.FAILURE, displaying original message from `scipy.optimize.minimize`: ABNORMAL_TERMINATION_IN_LNSRCH
computing pareto frontier : (GenericNoisyFunctionMetric('a'), GenericNoisyFunctionMetric('c'))
/home/ksuresh/miniconda3/envs/env_AID2E/lib/python3.12/site-packages/botorch/optim/fit.py:102: OptimizationWarning:
`scipy_minimize` terminated with status OptimizationStatus.FAILURE, displaying original message from `scipy.optimize.minimize`: ABNORMAL_TERMINATION_IN_LNSRCH
computing pareto frontier : (GenericNoisyFunctionMetric('b'), GenericNoisyFunctionMetric('c'))
/home/ksuresh/miniconda3/envs/env_AID2E/lib/python3.12/site-packages/botorch/optim/fit.py:102: OptimizationWarning:
`scipy_minimize` terminated with status OptimizationStatus.FAILURE, displaying original message from `scipy.optimize.minimize`: ABNORMAL_TERMINATION_IN_LNSRCH
print ("Metric_a, Metric_b")
render(plot_pareto_frontier(frontier[0], CI_level=0.8))
Metric_a, Metric_b
Validating the computed pareto front performance#
Since the model is trained on objectives, One can perform k-fold validation to see the performance of the surrgoate model’s prediction
from ax.modelbridge.cross_validation import cross_validate
from ax.plot.diagnostic import tile_cross_validation
#https://ax.dev/api/_modules/ax/modelbridge/cross_validation.html
cv = cross_validate(model, folds = 5)
render(tile_cross_validation(cv))
Exercise 3#
Determine the Pareto set from the 3D front and choose an optimal point
Plot the optimal configuration of the tracker corresponding to that point
Do analysis of convergence
Visualize the point with a radar or petal diagram, following https://pymoo.org/visualization/index.html