Tutorial: Own Problems and Models for Evolutionary Optimization¶
This tutorial demonstrates how to define own geometry models and optimization problems using the EO module of pyGDM. The evolutionary optimization module interfaces pagmo/pygmo.
- We will create a very simple model with only two free parameters: A cuboidal antenna of fixed height. The free parameters shall be the width and length of the rectangular cross-section.
- The problem we’ll implement below will be the maximization of the scattering cross-section.
We start by loading the modules we’ll use
import numpy as np from pyGDM2 import structures from pyGDM2 import materials from pyGDM2 import fields from pyGDM2 import core from pyGDM2 import linear from pyGDM2 import tools from pyGDM2 import visu ## --- the main EO routine from pyGDM2.EO.core import run_eo
Implementing a Geometry Model¶
We will implement a very simple model based on the
geometry in the
structures module. The free parameters will be
length and width of the rod, the height shall be fixed.
from pyGDM2.EO.models import BaseModel class ownModel(BaseModel): """optimization-model for a rectangular antenna Notes ----- The purpose of this model is only for demonstration. A problem with only two free parameters can more easily be solved by systematic analysis of all possible permutations (brute-force) A user-defined model must at least implement the following three methods: - `__init__(self, sim, **kwargs)` (constructor) - `get_bounds(self)` (optimization parameter box-boundaries) - `generate_structure(self, params)` (geometry generator) """ def __init__(self, sim, limits_W, limits_L, height): """ Mandatory actions in the constructor are: - the constructor of the parent class "BaseModel" must be called with the `simulation` object defining the pyGDM simulation part - call `self.random_structure()`, which is defined in "BaseModel" """ ## --- call BaseModel constructor, pass at least the "simulation" instance super(self.__class__, self).__init__(sim) ## --- width and length limits for rectrangular antenna self.limits_W = limits_W self.limits_L = limits_L ## --- fixed parameters self.height = height # height of rectangular antenna ## --- init with random values. `random_structure` is defined in `BaseModel` self.random_structure() def get_bounds(self): """Return lower and upper boundaries for parameters Both of the lower/upper boundary lists must contain a value for each free parameter. In our case (2 free parameters), thus the two lists contain 2 elements each: The box-boundaries for the width and length of the rectangle. """ self.lbnds = [self.limits_W, self.limits_L] self.ubnds = [self.limits_W, self.limits_L] return self.lbnds, self.ubnds def generate_structure(self, params): """generate the structure This method implements the actual structure generation as function of the free parameters, passed as numpy array `params`. The geometry is generated as a list of mesh-point positions (list/array of 3D-coordinate tuples (x,y,z)). The structure must be internally set using `self.set_structure()`, the latter being defined in `BaseModel`. """ ## --- the order of `params` must correspond to ## the boundary order, used in `get_bounds` W, L = params, params ## --- generate the structure geometry: rectangular nano-rod from pyGDM2 import structures geometry = structures.rect_wire(self.sim.struct.step, L, self.height, W) ## --- mandatory call: set new structure-geometry ## (`set_structure` is a `BaseModel` method) self.set_structure(geometry)
Implementing an Optimization Problem¶
We will also implement a simple problem for optimization. We want to find a structure geometry (in this demonstration: The optimum dimensions of the rectangular rod), maximizing the scattering cross-section under plane wave illumination for a specific wavelength and polarization. You will see below, that this is very simple. The problem of maximizing the scattering cross-section is defined by only 8 lines of actual code.
*Note:* pyGDM assumes that the objective function is to be *maximized*. The pygmo package assumes minimization, thus internally the objective function will be multiplied by -1 to effectively render the maximization to a minimization task. If you want to implement a minimization, you might apply the same trick: use the negative of the target value(s) as objective function.
from pyGDM2.EO.problems import BaseProblem class ownProblem(BaseProblem): """Problem to maximize scattering cross-section (SCS) of nano-structure Notes ----- A user-defined problem must at least implement the following methods: - `__init__(self, model, **kwargs)` (constructor) - `objective_function(self, params)` (optimization target) """ def __init__(self, model): """ In the Problem constructor, the constructor of the parent class "BaseProblem" must be called, passing the `model` object. """ super(self.__class__, self).__init__(model) def objective_function(self, params): """calculate the scattering cross-section This method implements the actual optimization problem by calculating the objective function, also often called "fitness". The objective function returns the optimization target as function of the parameter-set `params`. The fitness will usually be calculated as follows: - (1) the structure generator is called with the parameter- vector `params` (which is provided by the optimization algorithm) - (2) a pyGDM simulation is performed using the model-geometry defined by the parameters `params` - (3) the pyGDM results are evaluated and the result returned. Below, we calculate the scattering cross-section of the considered particle. """ ## --- generate structure definded by `params` self.model.generate_structure(params) ## --- GDM sim. / cross-section calculation (field_index=0) core.scatter(self.model.sim, verbose=False) ext_cs, sca_cs, abs_cs = linear.extinct(self.model.sim, 0) return float(sca_cs)
Running the Optimization¶
Now, let’s test our user-defined model and problem. Let’s search for the rectangular antenna geometry that optimally scatters light of \(1\,\mu\)m wavelength.
Setup the pyGDM-Simulation¶
We setup a pyGDM simulation just as usual, with the only difference, that the geometry is left blank, since this will be varied and determined during the optimization by the model-class, we defined above.
## --- Setup structure mesh = 'cube' step = 20 material = materials.gold() # particle material n1, n2 = 1.0, 1.0 # in vacuum ## Empty dummy-geometry, will be replaced on run-time by EO-model geometry =  struct = structures.struct(step, geometry, material, n1,n2, structures.get_normalization(mesh)) ## --- Setup incident field field_generator = fields.planewave # planwave excitation kwargs = dict(theta = [0.0]) # single (linear) polarization wavelengths =  # single wavelength efield = fields.efield(field_generator, wavelengths=wavelengths, kwargs=kwargs) ## --- Simulation initialization sim = core.simulation(struct, efield)
/home/wiecha/.local/lib/python2.7/site-packages/pyGDM2/structures.py:119: UserWarning: Emtpy structure geometry. warnings.warn("Emtpy structure geometry.")
Setup the Evolutionary Optimization¶
After having defined the simulation conditions, we need to setup the particle geometry, the target problem and evolutionary optimization algorithm.
- The geometry will be provided by our model class
- The problem is defined in our class
- The evolutionary optimization algorithm will be one of the pygmo package
# ============================================================================= # structure model # ============================================================================= model = ownModel(sim, limits_W=[2,20], limits_L=[2,20], height=1) # ============================================================================= # optimization problem # ============================================================================= problem = ownProblem(model) # ============================================================================= # setup optimization algorithm # ============================================================================= ## --- size of population population = 15 # Nr of individuals ## --- stop criteria max_time = 60 # seconds max_iter = 20 # max. iterations max_nonsuccess = 5 # max. consecutive iterations without improvement ## --- other config generations = 1 # Nr of generations to evolve between status reports save_all_generations = False filename = "own_problem.eo" # path to file to save results (optional) ## Use algorithm "sade" (jDE variant, a self-adaptive ## form of differential evolution) import pygmo as pg algorithm = pg.sade algorithm_kwargs = dict() # optional kwargs, to be passed to the algorithm # ============================================================================= # run the optimization # ============================================================================= eo_dict = run_eo(problem, population=population, algorithm=algorithm, generations=generations, max_time=max_time, max_iter=max_iter, max_nonsuccess=max_nonsuccess, filename=filename)
/home/wiecha/.local/lib/python2.7/site-packages/pyGDM2/structures.py:104: UserWarning: Minimum structure Z-value lies below substrate level! Shifting structure bottom to Z=step/2. " Shifting structure bottom to Z=step/2.")
---------------------------------------------- Starting new optimization ---------------------------------------------- iter # 1, time: 2.4s, progress # 0, f_evals: 30 (non-success: 1) - champion fitness: [-289300.0] iter # 2, time: 4.4s, progress # 1, f_evals: 45 - champion fitness: [-355440.0] iter # 3, time: 6.5s, progress # 1, f_evals: 60 (non-success: 1) iter # 4, time: 8.5s, progress # 2, f_evals: 75 - champion fitness: [-437960.0] iter # 5, time: 10.6s, progress # 2, f_evals: 90 (non-success: 1) iter # 6, time: 12.6s, progress # 2, f_evals: 105 (non-success: 2) iter # 7, time: 14.7s, progress # 3, f_evals: 120 - champion fitness: [-481860.0] iter # 8, time: 16.9s, progress # 3, f_evals: 135 (non-success: 1) iter # 9, time: 19.1s, progress # 3, f_evals: 150 (non-success: 2) iter # 10, time: 21.3s, progress # 3, f_evals: 165 (non-success: 3) iter # 11, time: 23.4s, progress # 3, f_evals: 180 (non-success: 4) iter # 12, time: 25.7s, progress # 3, f_evals: 195 (non-success: 5) -------- maximum non-successful iterations reached
Analyze and plot the results¶
pyGDM.EO.tools offers several tools for rapid and simple reloading
of the optimization results. We will use
obtain the final best solution of the optimization. This function
returns the pyGDM
simulation object containing the optimum geometry:
## --- some more imports import copy import matplotlib.pyplot as plt from pyGDM2.EO import tools as EOtools ## --- get the best candidate from the optimization sim = EOtools.get_best_candidate(eo_dict, verbose=True) print sim visu.structure(sim)
Best candidate after 3 iterations with improvement: fitness = ['-4.8186e+05'] Testing: recalculating fitness... Done. Everything OK. =============== GDM Simulation Information =============== precision: <type 'numpy.float32'> / <type 'numpy.complex64'> ------ nano-object ------- Homogeneous object. material: "Gold, Johnson/Christy" mesh type: cubic nominal stepsize: 20nm nr. of meshpoints: 209 ----- incident field ----- field generator: "planewave" 1 wavelengths between 1000.0 and 1000.0nm 1 incident field configurations per wavelength ------ environment ------- n3 = (1+0j) <-- top n2 = (1+0j) <-- structure zone (height "spacing" = 5000.0nm) n1 = (1+0j) <-- substrate ===== *core.scatter* ====== self-consistent fields are available
<matplotlib.collections.PathCollection at 0x7fb83469b4d0>
Calculate a Scattering Spectrum for the Resulting Structure¶
We hope to have obtained a nano-antenna, resonant at the target
wavelength (in our case \(1\mu\)m). To check this, we would like
to compute the scattering spectrum of the structure. Since the
simulation object was defined for only one single
illumination wavelength, we need to create a new pyGDM simulation
#============================================================================== # setup new simulation to calculate spectrum #============================================================================== ## --- structure struct = copy.deepcopy(sim.struct) ## --- incident field field_generator = fields.planewave # planwave excitation wavelengths = np.arange(500, 1500, 20) # spectrum kwargs = dict(theta = [0.0, 90.0]) # 0/90 deg polarizations efield = fields.efield(field_generator, wavelengths=wavelengths, kwargs=kwargs) ## --- spectrally resolved simulation sim_spectrum = core.simulation(struct, efield) #============================================================================== # run simulation for the spectrum #============================================================================== core.scatter(sim_spectrum, verbose=0) ## --- scattering spectrum field_kwargs = tools.get_possible_field_params_spectra(sim_spectrum) wl, spec_ext0 = tools.calculate_spectrum(sim_spectrum, field_kwargs, linear.extinct) asca0 = spec_ext0.T ## --- plot plt.figure(figsize=(10,5)) plt.subplot2grid((1,5), (0,0), colspan=3) plt.title("scattering spectrum") plt.plot(wl, asca0, label="plane wave, pol. || X") plt.legend(loc='best', fontsize=10) plt.xlabel("wavelength (nm)") plt.ylabel("scat.section (nm^2)") plt.subplot2grid((1,5), (0,3), colspan=2, aspect="equal") plt.title('structure geometry') visu.structure(sim_spectrum, show=False) plt.show()
The structure found by evolutionary optimization has indeed a scattering resonance at \(1000\,\)nm wavelength. Furthermore, since we aimed to maximize the cross-section, the \(Y\)-dimension of the rectangle is chosen as large as possible (within the allowed limits), to increase the cross-section. Maximizing the scattering efficiency instead, would lead to a very narrow antenna of same length, in order to reduce the geometric cross-section.
Note: Technical Detail¶
Most EO-tools can be called either with the
eo_dict as obtained by
run_eo or with the filename to the stored optimizaiton results. For
sim = EOtools.get_best_candidate(eo_dict, verbose=True)
Best candidate after 3 iterations with improvement: fitness = ['-4.8186e+05'] Testing: recalculating fitness... Done. Everything OK.
is equivalent to calling
get_best_candidate, giving the path to the
filename = "own_problem.eo" # path to saved EO results sim = EOtools.get_best_candidate(filename, verbose=True)
Best candidate after 3 iterations with improvement: fitness = ['-4.8186e+05'] Testing: recalculating fitness... Done. Everything OK.