Bayesian Optimization#

# Execute once during the start
!git clone https://github.com/cfteach/modules.git
!pip install scikit-optimize 

%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

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 = 0.5 #.5  
pitch = 4.0 #10.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)


detector2.geometry_display(Z, Y, R, y_min=y_min, y_max=y_max,block=False,pause=5) #5

print("# of wires: ", num_wires)

#------------- TRACK GENERATION -----------#
print(":::: TRACK GENERATION ::::")
t = detector2.Tracks(b_min= y_min, b_max=y_max, alpha_mean=0, alpha_std=0.2) #b_min=-100, b_max=100
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]
print("fraction of tracks detected: ",frac_detected)

::::: 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
:::: TRACK GENERATION ::::

fraction of tracks detected:  0.264

Define Objectives#

Defines the objective function of the problem that can be used in the Bayesian Optimization.

#------------- OBJECTIVE FUNCTION ---------------#

class Obj:

  def __init__(self,R=0.5,pitch=4.0):
    self.R = R 
    self.pitch = pitch  

  def update_pars(self,r,p):
    self.R = r 
    self.pitch = p  
    

  def objective(self,x):

    #you have to upack R, pitch if included in the BO
    y1, y2, y3, z1, z2, z3 = x 
    Z, Y = detector2.Tracker(self.R, self.pitch, y1, y2, y3, z1, z2, z3).create_geometry()
    #a track is detected if at least two wires have been hit
    val = detector2.get_score(Z, Y, tracks, self.R)[0] 
    #print(self.R, self.pitch)

    return 1. - val  # the smaller, the better.

Bayesian Optimization of the Detector Geometry#

You can find different options for the optimization.

We refer to sklearn optimization package: https://scikit-optimize.github.io/stable/

Different approaches for the regression part.

from skopt import gp_minimize, gbrt_minimize, dummy_minimize

# You can change values within pre-defined ranges 
# Below we are considering only the shifts along y and z 
# (of each plane with respect to the previous plane)

y1_min, y1_max = [0., 4.]
y2_min, y2_max = [0., 4.]
y3_min, y3_max = [0., 4.]
z1_min, z1_max = [2., 10.]
z2_min, z2_max = [2., 10.]
z3_min, z3_max = [2., 10.]

dims = [(y1_min, y1_max),(y2_min, y2_max), (y3_min, y3_max), 
        (z1_min, z1_max), (z2_min, z2_max),(z3_min, z3_max)]

ncalls = 100
rand_st = 1717 #for reproducibility