Maximum Likelihood Fit ( i.e., a \(\chi^2\)-fit)
The first example will be a simple \(\chi^2\)-fit, where the optical model parameters will be adjusted to reproduce observed elastic scattering data. In this case I am going to look at the \(^{48}\) Ca \((p, p)\) at \(25\) MeV. This data comes from Ref. [Mc1986]. Just selecting the … optical potential, we first have to create a FRESCO file that is compatible with pfunk. In particular, pfunk swaps variable values based on line numbers. While this can seem cumbersome, the upshot is that swapping based on line number eliminates any computational overhead with doing find and replace. Below the relevant lines are highlighted.
0 48Ca(p,p) 25.0 MeV
1 NAMELIST
2 &fresco
3 hcm = 0.1
4 rmatch = 40.0
5 jtmax = 20
6 absend = 0
7 jump(1:6) = 0, 0, 0, 0, 0, 0
8 jbord(1:6) = 0, 0, 0, 0, 0.0, 0.0
9 kqmax = 1
10 thmin = 0.0
11 thmax = 180.0
12 thinc = 1.0
13 it0 = 1
14 smats = 2
15 xstabl = 1
16 elab(1:4) = 25.0, 0, 0, 0
17 nlab(1:3) = 0, 0, 0
18 /
19
20 &partition
21 namep = 'p'
22 massp = 1
23 zp = 1
24 namet = 'Ca'
25 masst = 48
26 zt = 20
27 qval = 0.0
28 pwf = .false.
29 nex = 1
30 /
31
32
33 &states
34 jp = 0.5
35 bandp = 1
36 cpot = 1
37 jt = 0.0
38 bandt = 1
39 /
40
41
42 &partition
43 /
44
45 &pot
46 kp = 1
47 at = 48
48 rc = 1.25
49 /
50
51 &pot
52 kp = 1
53 type = 1
54 p1 = 52.2
55 p2 = 1.17
56 p3 = 0.75
57 p4 = 2.8
58 p5 = 1.32
59 p6 = 0.63
60 /
61
62 &pot
63 kp = 1
64 type = 2
65 p4 = 7.6
66 p5 = 1.32
67 p6 = 0.63
68 /
69
70 &pot
71 kp = 1
72 type = 3
73 shape = 0
74 p1 = 6.2
75 p2 = 1.01
76 p3 = 0.75
77 /
78
79 &pot
80 /
81
82 &overlap
83 /
84
85 &coupling
86 /
In order for pfunk to swap these values we need to specify the input file path, line number, and variable name. The most convenient way to do this for the line number and names is using lists.
fresco_path = './48Ca_elastic_new.in'
fresco_names = ['p1', 'p2', 'p3', 'p4', ('p5', 'p5'), ('p6', 'p6')]
fresco_positions = [54, 55, 56, 57, (58, 66), (59, 67)]
Note that the line numbers above start at zero. The names and numbers in parentheses tell pfunk to swap one number in at two positions in the input file. This is done because often global optical potentials use the same radius and diffuseness parameters for the imaginary surface and volume potentials:
The parenthesis enforce this condition by swapping the enclosed parameters simultaneously.
Now it is time to initialize our model. The first step is to create an instance
of a pfunk.model.Model object:
model = pfunk.model.Model(fresco_path, fresco_names, fresco_positions)
Since this is a \(\chi^2\)-fit the only remaining task is to assign a likelihood function. The trick is that a \(\chi^2\)-fit is simply the log-likelihood of a normal distribution. There are several methods in the model object to assist in the creation of the necessary probability functions. For a likelihood function all that needs to be specified is the path of the fresco output, the path to the data that is to be fit, and any additional model parameters (which I will come back to later). The data we want to fit should be in a text file with labeled columns of angle, cross section, and uncertainty. As an example here are the first four lines of the data file.
theta sigma erry
10.19108 0.80674 0.080674
12.95117 0.89819 0.089819
15.28662 1.012 0.1012
Next, by default, the FRESCO elastic scattering cross section is written to fort.201.
To add a likelihood function based on elastic scattering call the method
pfunk.model.Model.create_elastic_likelihood():
model.create_elastic_likelihood('fort.201', elastic_data_path)
Finally, pfunk needs to know that we are done adding likelihood functions so that it can create the total likelihood function.
model.create_likelihood()
We now what to fit the optical model parameters to the data. Minimization is
provided by the module pfunk.model_fit , and specifically a model can be
minimized using the class pfunk.model_fit.MAPFit. First the model needs
to have starting parameters. The easiest way to do this if we just have potential
parameters is to pull them from the pfunk.fresco_classes.NamelistInput object which is generated automatically by the model. This object stores the initial parameters, and this list can be copied to the model objects initial parameters with:
model.x0 = model.fresco.x0[:]
We now initialize the fit object:
fit = model.funk.model_fit.MAPFit(model, lnlike=True, percent_range=5.0)
where we have told MAPFit object to minimize just the likelihood via lnlike=True and to
consider a global minimization in a range of 5.0*model.x0 up and down with percent_range=5.0.
Because model.x0 controls the search range of the global minimization, different starting values
can be chosen by changing the array fit.x0_custom. With this we simply need to run the minimization, which
can be done using basin hopping, pfunk.model_fit.MAPFit.run(), differential evolution,
pfunk.model_fit.MAPFit.run_differential(), or dual annealing, pfunk.model_fit.MAPFit.run_anneal().
Annealing usually provides the best compromise between speed and search range:
fit.run_anneal(max_iter=1000)
This will print out a message after each iteration, but the message will not make too much sense. The first value
is the current value of likelihood function, which you should check to make sure everything is still going smoothly.
After this finishes the minimized parameter values can be found in the array fit.results.x. An example of this
calculation can be found in this example.