CH82 – optimization¶
Input: Defines the model and constraints
from HJCFIT import read_idealized_bursts
from HJCFIT.likelihood import QMatrix
name = "CH82.scn"
tau = 1e-4
tcrit = 4e-3
graph = [["V", "V", "V", 0, 0],
["V", "V", 0, "V", 0],
["V", 0, "V", "V", "V"],
[ 0, "V", "V", "V", 0],
[ 0, 0, "V", 0, "V"]]
nopen = 2
qmatrix = QMatrix([[ -3050, 50, 3000, 0, 0 ],
[ 2./3., -1502./3., 0, 500, 0 ],
[ 15, 0, -2065, 50, 2000 ],
[ 0, 15000, 4000, -19000, 0 ],
[ 0, 0, 10, 0, -10 ] ], 2)
bursts = read_idealized_bursts(name, tau=tau, tcrit=tcrit)
Creates the constraints, the likelihood function, as well as a function to create random Q-matrix.
from scipy.optimize import minimize
from numpy import nan, zeros, arange
import numpy as np
from HJCFIT.likelihood.random import qmatrix as random_qmatrix
from HJCFIT.likelihood import QMatrix, Log10Likelihood
from HJCFIT.likelihood.optimization import reduce_likelihood
likelihood = Log10Likelihood(bursts, nopen, tau, tcrit)
reduced = reduce_likelihood(likelihood, graph)
x = reduced.to_reduced_coords( random_qmatrix(5).matrix )
constraints = []
def create_inequality_constraints(i, value=0e0, sign=1e0):
f = lambda x: sign * (x[i] - value)
def df(x):
a = zeros(x.shape)
a[i] = sign
return a
return f, df
for i in range(len(x)):
f, df = create_inequality_constraints(i)
constraints.append({'type': 'ineq', 'fun': f, 'jac': df})
f, df = create_inequality_constraints(i, 1e4, -1)
constraints.append({'type': 'ineq', 'fun': f, 'jac': df})
def random_starting_point():
from numpy import inf, nan
from HJCFIT.likelihood.random import rate_matrix as random_rate_matrix
for i in range(100):
matrix = random_rate_matrix(N=len(qmatrix.matrix), zeroprob=0)
x = reduced.to_reduced_coords( matrix )
try:
result = reduced(x)
print(result, reduced.to_full_coords(x))
except:
pass
else:
# `result != nan` was always True -- NaN compares unequal to
# everything, itself included -- so a NaN likelihood used to be
# accepted as a valid starting point.
if np.isfinite(result): break
else: raise RuntimeError("Could not create random matrix")
return x
def does_not_throw(x):
try: return -reduced(x)
except: return nan
Performs the minimization
import math
methods = ['COBYLA', 'SLSQP']
x = random_starting_point()
print ('x=', x)
# `does_not_throw` returns -log10L, so smaller is better. The running best has
# to be held in that same convention: it used to be seeded with reduced(x),
# a positive log10 likelihood, and then compared against result.fun, its
# negation, so the comparison was between incompatible quantities and was
# very nearly always true.
best = (x.copy(), does_not_throw(x))
for i in range(len(methods)):
result = minimize(does_not_throw,
x,
method=methods[i],
constraints=constraints,
options={'maxiter': 1000, 'disp':True})
print(result)
if not math.isnan(result.fun):
# result.x, not x: x is where this attempt started, not where it
# finished. The best parameters reported used to be a starting point.
if result.fun < best[1]: best = (result.x.copy(), result.fun)
# `and i > 4` could never be true with two methods in the list.
if result.success: break
x = x + random_starting_point() * 1e-2
if np.all(np.isnan(x)): x = random_starting_point()
print(best[0])
print('log10 likelihood =', -best[1])
-640.830069172272 [[ -6.21619693e-01 3.80010323e-01 2.41609369e-01 0.00000000e+00
0.00000000e+00]
[ 3.10903942e+03 -3.10919359e+03 0.00000000e+00 1.54171557e-01
0.00000000e+00]
[ 1.59277846e-01 0.00000000e+00 -9.08913654e+03 9.08867581e+03
3.01453712e-01]
[ 0.00000000e+00 4.39644066e-01 2.78736992e-01 -7.18381058e-01
0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 1.54808746e-01 0.00000000e+00
-1.54808746e-01]]
x= [ 3.80010323e-01 2.41609369e-01 3.10903942e+03 1.54171557e-01
1.59277846e-01 9.08867581e+03 3.01453712e-01 4.39644066e-01
2.78736992e-01 1.54808746e-01]
fun: -2062.8258070089187
maxcv: 8.7670065147940707e-16
message: 'Maximum number of function evaluations has been exceeded.'
nfev: 1000
status: 2
success: False
x: array([ -8.76700651e-16, 1.75097884e+02, 3.11563723e+03,
2.68478978e+02, 6.06596893e+02, 9.06857871e+03,
1.73472348e-18, 1.77190965e+01, -3.03804457e-16,
7.92641819e-01])
-697.0699667052597 [[ -2.18554574e-01 8.96065006e-02 1.28948074e-01 0.00000000e+00
0.00000000e+00]
[ 6.92249289e+02 -6.92781143e+02 0.00000000e+00 5.31853771e-01
0.00000000e+00]
[ 8.75734586e-01 0.00000000e+00 -1.91477903e+00 6.89175466e-01
3.49868979e-01]
[ 0.00000000e+00 5.54235637e-01 4.39234514e-02 -5.98159089e-01
0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 9.79442657e+02 0.00000000e+00
-9.79442657e+02]]
Inequality constraints incompatible (Exit mode 4)
Current function value: -2284.629372492554
Iterations: 177
Function evaluations: 2189
Gradient evaluations: 177
fun: -2284.629372492554
jac: array([ -2.08709717e-01, -5.42224910e+08, -2.52990723e-02,
2.75032878e+05, 1.75594303e+09, -1.57243136e+09,
-5.42756597e+08, -2.75028975e+05, 5.62684071e+08,
-6.88560304e+07, 0.00000000e+00])
message: 'Inequality constraints incompatible'
nfev: 2189
nit: 177
njev: 177
status: 4
success: False
x: array([ 4.64593912e-07, 3.41051917e+02, 2.65544934e+03,
1.38876283e+03, 9.99999918e+03, 4.32621889e+03,
2.19612726e+02, 2.99939102e+00, 2.05470191e+00,
-1.68337691e-14])
-447.89702673666727 [[ -2.96248722e+03 2.73428480e-01 2.96221379e+03 0.00000000e+00
0.00000000e+00]
[ 9.99336618e-01 -1.00764635e+00 0.00000000e+00 8.30973004e-03
0.00000000e+00]
[ 5.19874297e-03 0.00000000e+00 -5.01074858e+03 3.97100105e-01
5.01034628e+03]
[ 0.00000000e+00 8.74501090e-01 2.32778358e+03 -2.32865808e+03
0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 6.91322340e-03 0.00000000e+00
-6.91322340e-03]]
[ 3.80906388e-01 2.42898850e-01 3.11596192e+03 1.59490095e-01
1.68035192e-01 9.08868270e+03 3.04952401e-01 4.45186422e-01
2.79176226e-01 9.94923531e+00]
-2284.629372492554