import escape as esc
import numpy as np
esc.require("0.9.7")
from escape.utils.widgets import show

Loading material database from /home/dkor/Data/Development/workspace_escape/escape/python/src/escape/scattering/../data/mdb/materials.db

Models and Optimizers

In the previous notebooks we have demonstrated two major building blocks of ESCAPE package: parameters and functors. In this notebook we discuss objects which are relevant to optimization.

Kernel

Kernel does very little but very important job, it vectorizes a given functor, i.e. it evaluates functor over successive elements of the input array. And it supports multithreading, making all CPU cores busy.

Below we create two functors with shared parameters and corresponding kernels. This kernels we use to create "experimental data" which we will fit later.

#identity factor
X=esc.var("X")

c=esc.par("c", 0.5)
b=esc.par("b", 5)
d=esc.par("d", 0.5 )
a=esc.par("a", 100 )

#a and b are shared parameters
expr1=a*esc.pow(esc.cos(X), 2.0)*esc.exp(-d*X)+b
expr2=a*esc.pow(esc.sin(X), 2.0)*esc.exp(-c*X)+b

k1=esc.kernel("Expr1 kernel", expr1, multithreaded=True)
k2=esc.kernel("Expr2 kernel", expr2, multithreaded=True)

x=np.arange(-10, 10, 0.01, dtype=float)

exp1=np.zeros(x.shape, dtype=float)
exp2=np.zeros(x.shape, dtype=float)

k1(x, exp1)
k2(x, exp2)


show(k1, coordinates=x, ylog=True)
show(k2, coordinates=x, ylog=True)

Next we create data objects using the generated arrays with poisson noise.

#Data

#adding some noise to data
exp1 = np.random.poisson(exp1)
exp2 = np.random.poisson(exp2)

dt1=esc.data("Data1", x, exp1, copy=True)
dt2=esc.data("Data2", x, exp2, copy=True)

show(dt1, ylog=True)
show(dt2, ylog=True)

Model

Model is a container for the kernel and the experimental data. It is responsible for calculating residuals and cost, which are necessary to optimize all parameters provided by the user. The model has two settings, weight_type=("none", "data") and residuals_scale=("none", "log"). The latter is used to calculate residuals and the former indicates which weights will be used to calculate cost.

mdl1=esc.model("Model1", k1, dt1, weight_type="none", residuals_scale="none")
mdl2=esc.model("Model2", k2, dt2, weight_type="none", residuals_scale="none")
show(mdl1, ylog=True)
show(mdl2, ylog=True)

Optimizer

Currently there are two optimization algorithms supported by ESCAPE.

The Levenberg-Marquardt algorithm

It has been quite a long period of time the standard method for non-linear data fitting and still is very usefull for many applications. It starts at the initial value of the cost function and, as a gradient descent trust region method, steps in the direction of the derivative until it reaches the local minimum. The cost function is the sum of square differences between theory and model. This algorithm uses a numerical approximation of the Jacobian matrix to set the step direction and an adaptive algorithm for the size of the trust region.

One should use this method when there is a reasonable fit exist near the minimum and one would like to get the best value for the optimized parameters. Compared to stochastic methods, LM converges much faster.

opt=esc.levmar("LM optimizer", [mdl1, mdl2], xtol = 1e-10, nupdate=1)

Before optimization we will "shake" values of our parameters. That means that the starting values of the parameters will be chosen randomly but within their limits. They have to be explicitly selected by the user as shown below.

opt.shake()
opt()
show(opt, ylog=True)

opt

Name: LM optimizer        	Parameters number:          4
Parameter           	Value          	+-	Error     	Units     	Fixed
a                   	         100.21	+-	0.0090089 	          	    0
d                   	        0.49972	+-	9.5397e-06	          	    0
b                   	         5.1368	+-	0.018523  	          	    0
c                   	         0.4996	+-	1.0996e-05	          	    0

Differential evolution

Differential evolution is a stochastic population based algorithm which uses differences between solutions points as a guide to selecting new points for new solution. A difference vector is computed for a pair of randomly chosen points for each member of the population. Based on crossover setting a random subset of vector's components are added to the current point. If the cost of this new point is smaller than the cost of the current point, it replaces the current point. Convergence with differential evolution will be slower, but more robust.

ESCAPE user can boost the convergence speed of DE algorithm using the following polishing settings:

- polish_final_maxiter - maximum number of iterations for LM algorithm after DE
- polish_candidate_maxiter - maximum number of iterations for LM algorithm for each candidate

Each time the candidate is chosen by the DE algorithm, LM method is started to slightly improve its value in the direction of the derivative. We recommend to use small number of iterations for the candidate polishing.

opt=esc.diffevol("DE optimizer", [mdl1, mdl2], maxiter=100,
                  polish_final_maxiter=50, 
                  polish_candidate_maxiter=2, minconv=1e-3, nupdate=1)

opt.shake()
opt()
show(opt, ylog=True)

opt

Name: DE optimizer        	Parameters number:          4
Parameter           	Value          	+-	Error     	Units     	Fixed
a                   	         100.21	+-	0.0090089 	          	    0
d                   	        0.49972	+-	9.5397e-06	          	    0
b                   	         5.1368	+-	0.018523  	          	    0
c                   	         0.4996	+-	1.0996e-05	          	    0