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

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)

##                      form of differential evolution)
import pygmo as pg
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

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.