-
Notifications
You must be signed in to change notification settings - Fork 0
/
minimize.py
82 lines (72 loc) · 2.5 KB
/
minimize.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
import numpy as np
from scipy.optimize import minimize
from plot import plot
# Array of (x,y) points
xdata = []
with open("fit01.dat") as f:
for line in f:
line = line.split()
xdata.append((float(line[0]), float(line[1])))
# Array of values at (x,y)
ydata = [0.0 for x in range(len(xdata))]
# Array of experimental error at (x,y)
sigma = [1.0 for x in range(len(xdata))]
# Initial guess of fit parameters
x0 = [2.29663137951, 0.00702885711719, 2.09847731881e-07, -8.84759981594e-08, -2.99660935836e-07, -0.000417687931814]
# Define the objective function to be minimised
# params ... array holding the values of the fit parameters
# X ... array holding inputs of observed data
# Y ... array holding values of observed data
# Err ... array holding errors of observed data
def func(params, X, Y, Err):
# Extract current values of fit parameters
t1 = 1
t2, mu, x1, x2, epsf, V = params
# Compute sum of chi-squared values
chi2a = 0.0
chi2b = 0.0
for n in range(len(X)):
k = X[n]
epsk = -2*t1*(np.cos(k[0])+np.cos(k[1]))-4*t2*np.cos(k[0])*np.cos(k[1])-mu
x0 = -2*x1*(np.cos(k[0])+np.cos(k[1]))-4*x2*np.cos(k[0])*np.cos(k[1])-epsf
det = np.sqrt(0.25*(epsk-x0)*(epsk-x0)+V*V)
Ea = 0.5*(epsk+x0)+det
Eb = 0.5*(epsk+x0)-det
chi2a = chi2a + (Y[n] - Ea)*(Y[n] - Ea)#/(Err[n]*Err[n]) # Uncomment to include error
chi2b = chi2b + (Y[n] - Eb)*(Y[n] - Eb)#/(Err[n]*Err[n]) # Uncomment to include error
return min(chi2a,chi2b)
# Optimize using one of the methods listed below
methods = [
# 'Newton-CG',
# 'dogleg',
# 'trust-ncg',
'Nelder-Mead',
'Powell',
'CG',
'BFGS',
'L-BFGS-B',
'TNC',
'COBYLA',
'SLSQP'
]
result = minimize(func, x0, args=(xdata, ydata, sigma), method=methods[0])
params = result.x
output = "[" + str(params[0]) + ", " \
+ str(params[1]) + ", " \
+ str(params[2]) + ", " \
+ str(params[3]) + ", " \
+ str(params[4]) + ", " \
+ str(params[5]) + "]\n"
with open("data.txt", "a") as f:
f.write(output)
output = "fitt1=" + str(1) + ";\n" + \
"fitt2=" + str(params[0]) + ";\n" + \
"fitmu=" + str(params[1]) + ";\n" + \
"fitx1=" + str(params[2]) + ";\n" + \
"fitx2=" + str(params[3]) + ";\n" + \
"fitf=" + str(params[4]) + ";\n" + \
"fitV=" + str(params[5]) + ";"
output = output.replace("e","*10^")
print output
# Plot energy dispersion
plot(params)