import escape as esc
import numpy as np

esc.require("0.9.8")

Loading material database from C:\dev\escape-core\python\src\escape\scattering\..\data\mdb\materials.db

Models and Optimizers

In the previous notebooks we introduced parameters and functors. Here we connect those pieces to the fitting workflow: a model compares a theoretical functor against experimental data and adjusts parameters to reduce the mismatch.

Functors and data

Functors are used directly: a model takes a functor and experimental data. The functor is evaluated over the data coordinates (vectorized, with optional multithreading; see the functors notebook for num_threads and array evaluation). In practice, this is the step where a mathematical description of the sample or instrument is turned into numbers that can be compared against measured data. Below we create two functors with shared parameters and use them to generate "experimental" data that we will fit later.

# Functors and coordinates: create two functors with shared parameters, then evaluate over x to get "experimental" curves
X = esc.var("X")
a = esc.par("a", 100, userlim=[50, 150])
b = esc.par("b", 4, userlim=[0, 10])
c = esc.par("c", 0.5, userlim=[0.1, 2])
d = esc.par("d", 0.5, userlim=[0, 2])

f1 = a * esc.exp(-((X - b) ** 2) / c) + d
f2 = a * esc.exp(-((X - b) ** 2) / (2 * c)) + d

x = np.linspace(0, 10, 200)
exp1 = f1(x)
exp2 = f2(x)

C:\Users\User\AppData\Local\Temp\ipykernel_6028\754473634.py:12: UserWarning: Input array does not own its data (e.g. it is a view or slice); data will be copied
  exp1 = f1(x)
C:\Users\User\AppData\Local\Temp\ipykernel_6028\754473634.py:13: UserWarning: Input array does not own its data (e.g. it is a view or slice); data will be copied
  exp2 = f2(x)

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(x, exp1, copy=True)
dt2 = esc.data(x, exp2, copy=True)

w1 = dt1.show().config(ylog=True)
w2 = dt2.show().config(ylog=True)

esc.show(w1, w2)

Model

Model is a container for the functor 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(f1, dt1, weight_type="data", residuals_scale="none", name="Model1")
mdl2 = esc.model(f2, dt2, weight_type="data", residuals_scale="none", name="Model2")

w1 = mdl1.show().config(ylog=True)
w2 = mdl2.show().config(ylog=True)

esc.show(w1, w2)

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([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()
opt.show().config_model(ylog=True)

opt

Name: Levenberg-Marquardt 	Parameters number:          4
Parameter           	Value          	+-	Error     	Units     	Fixed
a                   	         100.43	+-	1.6326    	          	    0
b                   	         4.0034	+-	0.0077272 	          	    0
c                   	        0.45586	+-	0.0092425 	          	    0
d                   	        0.38051	+-	0.064225  	          	    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(
    [mdl1, mdl2],
    maxiter=100,
    polish_final_maxiter=50,
    polish_candidate_maxiter=0,
    ftol=0.5,
    nupdate=1,
)

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