# 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

```
In [1]:
```

```
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 `rect_wire`

geometry in the `structures`

module. The free parameters will be
length and width of the rod, the height shall be fixed.

```
In [2]:
```

```
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[0], self.limits_L[0]]
self.ubnds = [self.limits_W[1], self.limits_L[1]]
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[0], params[1]
## --- 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.

```
In [3]:
```

```
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.

```
In [4]:
```

```
## --- 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 = [1000] # 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
`ownModel`

- The problem is defined in our class
`ownProblem`

- The evolutionary optimization algorithm will be one of the pygmo package

```
In [5]:
```

```
# =============================================================================
# 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 `get_best_candidate`

to
obtain the final best solution of the optimization. This function
returns the pyGDM `simulation`

object containing the optimum geometry:

```
In [6]:
```

```
## --- 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
```

```
Out[6]:
```

```
<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
optimization `simulation`

object was defined for only one single
illumination wavelength, we need to create a new pyGDM simulation
instance:

```
In [7]:
```

```
#==============================================================================
# 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)[0]
wl, spec_ext0 = tools.calculate_spectrum(sim_spectrum, field_kwargs, linear.extinct)
asca0 = spec_ext0.T[1]
## --- 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
instance,

```
In [8]:
```

```
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
stored results:

```
In [9]:
```

```
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.
```