I'm writing code to fit parameters of related pairs of functions from a family using the scipy.optimize.least_squares() function and it seems that the test parameters are not being passed to the function correctly. Here's a simplified illustration of the problem.
import numpy as np
import numpy.linalg as la
import scipy.optimize
import math
Ha = {'H': lambda x, a, b : a + b*x, 'nParams': 2}
Hb = {'H': lambda x, a, b, c : a + b*x + c*x*x, 'nParams': 3}
Hc = {'H': lambda x, a, b, c, d : a*math.cos(b*(x-c)) + d, 'nParams': 4}
points1 = [(1,0), (2,4), (3,5), (7,8)]
points2 = [(1,1), (2,6), (4,1), (6,6)]
def coupled_1 (x, g, *params, **kwargs):
H1, H2 = kwargs.values() # H1 and H2 are dictionaries containing the two curves
# to be fit, 'H', and their respective numbers of
# parameters, 'nParams'
matrix = np.diag([H1['H'](x, *params[:H1['nParams']]), H2['H'](x, *params[H1['nParams']:])]) \
+ np.diag([g],k=1) + np.diag([g],k=-1)
return la.eigh(matrix)['eigenvalues'][0]
def coupled_2 (x, g, *params, **kwargs):
H1, H2 = kwargs.values()
matrix = np.diag([H1['H'](x, *params[:H1['nParams']]), H2['H'](x, *params[H1['nParams']:])]) \
+ np.diag([g],k=1) + np.diag([g],k=-1)
return la.eigh(matrix)['eigenvalues'][1]
def getParameters(H1, H2, pts1, pts2):
nParams = H1['nParams'] + H2['nParams']
# These points sometimes come in pairs and sometimes don't so I've found no better way to zip them
pts = [(p1, p2) for p1 in pts1 for p2 in pts2 if p1[0] == p2[0]]
res = lambda *params, **kwargs: \
np.sqrt(sum( [( p[0][1]-coupled_1(p[0][0], *params, **kwargs) )**2 \
+ (p[1][1]-coupled_2(p[0][0], *params, **kwargs) )**2 for p in pts] ))
result = scipy.optimize.least_squares(res,[1] + [0]*nParams, kwargs={'H1':H1,'H2':H2})
return result['x']
params = getParameters(Ha, Hc, points1, points2)
I'm receiving the following error from Ha()
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Cell In[5], line 37
34 result = scipy.optimize.least_squares(res,[1] + [0]*nParams, kwargs={'H1':H1,'H2':H2})
35 return result['x']
---> 37 params = getParameters(Ha, Hc, points1, points2)
38 print(params)
Cell In[5], line 34, in getParameters(H1, H2, pts1, pts2)
30 pts = [(p1, p2) for p1 in pts1 for p2 in pts2 if p1[0] == p2[0]]
31 res = lambda *params, **kwargs: \
32 np.sqrt(sum( [( p[0][1]-coupled_1(p[0][0], *params, **kwargs) )**2 \
33 + (p[1][1]-coupled_2(p[0][0], *params, **kwargs) )**2 for p in pts] ))
---> 34 result = scipy.optimize.least_squares(res,[1] + [0]*nParams, kwargs={'H1':H1,'H2':H2})
35 return result['x']
File ~\anaconda3\lib\site-packages\scipy\optimize\_lsq\least_squares.py:830, in least_squares(fun, x0, jac, bounds, method, ftol, xtol, gtol, x_scale, loss, f_scale, diff_step, tr_solver, tr_options, jac_sparsity, max_nfev, verbose, args, kwargs)
827 if method == 'trf':
828 x0 = make_strictly_feasible(x0, lb, ub)
--> 830 f0 = fun_wrapped(x0)
832 if f0.ndim != 1:
833 raise ValueError("`fun` must return at most 1-d array_like. "
834 "f0.shape: {0}".format(f0.shape))
File ~\anaconda3\lib\site-packages\scipy\optimize\_lsq\least_squares.py:825, in least_squares.<locals>.fun_wrapped(x)
824 def fun_wrapped(x):
--> 825 return np.atleast_1d(fun(x, *args, **kwargs))
Cell In[5], line 32, in getParameters.<locals>.<lambda>(*params, **kwargs)
29 # These points sometimes come in pairs and sometimes don't so I've found no better way to zip them
30 pts = [(p1, p2) for p1 in pts1 for p2 in pts2 if p1[0] == p2[0]]
31 res = lambda *params, **kwargs: \
---> 32 np.sqrt(sum( [( p[0][1]-coupled_1(p[0][0], *params, **kwargs) )**2 \
33 + (p[1][1]-coupled_2(p[0][0], *params, **kwargs) )**2 for p in pts] ))
34 result = scipy.optimize.least_squares(res,[1] + [0]*nParams, kwargs={'H1':H1,'H2':H2})
35 return result['x']
Cell In[5], line 32, in <listcomp>(.0)
29 # These points sometimes come in pairs and sometimes don't so I've found no better way to zip them
30 pts = [(p1, p2) for p1 in pts1 for p2 in pts2 if p1[0] == p2[0]]
31 res = lambda *params, **kwargs: \
---> 32 np.sqrt(sum( [( p[0][1]-coupled_1(p[0][0], *params, **kwargs) )**2 \
33 + (p[1][1]-coupled_2(p[0][0], *params, **kwargs) )**2 for p in pts] ))
34 result = scipy.optimize.least_squares(res,[1] + [0]*nParams, kwargs={'H1':H1,'H2':H2})
35 return result['x']
Cell In[5], line 17, in coupled_1(x, g, *params, **kwargs)
14 H1, H2 = kwargs.values() # H1 and H2 are dictionaries containing the two curves
15 # to be fit, 'H', and their respective numbers of
16 # parameters, 'nParams'
---> 17 matrix = np.diag([H1['H'](x, *params[:H1['nParams']]), H2['H'](x, *params[H1['nParams']:])]) \
18 + np.diag([g],k=1) + np.diag([g],k=-1)
19 return la.eigh(matrix)['eigenvalues'][0]
TypeError: <lambda>() missing 2 required positional arguments: 'a' and 'b'
It seems that the *params list is not making it to Ha(), but I'm unsure why.
Here is a different approach for the res function, this works for me.
import numpy as np
import numpy.linalg as la
import scipy.optimize
import math
# Define the dictionaries for the functions
Ha = {'H': lambda x, a, b: a + b*x, 'nParams': 2}
Hb = {'H': lambda x, a, b, c: a + b*x + c*x*x, 'nParams': 3}
Hc = {'H': lambda x, a, b, c, d: a*math.cos(b*(x-c)) + d, 'nParams': 4}
points1 = [(1,0), (2,4), (3,5), (7,8)]
points2 = [(1,1), (2,6), (4,1), (6,6)]
def coupled_1(x, g, *params, **kwargs):
H1, H2 = kwargs['H1'], kwargs['H2']
matrix = np.diag([H1['H'](x, *params[:H1['nParams']]), H2['H'](x, *params[H1['nParams']:])]) \
+ np.diag([g], k=1) + np.diag([g], k=-1)
return la.eigh(matrix)[0][0] # Get the smallest eigenvalue
def coupled_2(x, g, *params, **kwargs):
H1, H2 = kwargs['H1'], kwargs['H2']
matrix = np.diag([H1['H'](x, *params[:H1['nParams']]), H2['H'](x, *params[H1['nParams']:])]) \
+ np.diag([g], k=1) + np.diag([g], k=-1)
return la.eigh(matrix)[0][1] # Get the second smallest eigenvalue
def getParameters(H1, H2, pts1, pts2):
nParams = H1['nParams'] + H2['nParams']
# Filter points with matching x-values
pts = [(p1, p2) for p1 in pts1 for p2 in pts2 if p1[0] == p2[0]]
def res(params, **kwargs):
error_sum = 0
for p1, p2 in pts:
x_val = p1[0]
y1 = p1[1]
y2 = p2[1]
error_sum += (y1 - coupled_1(x_val, *params, **kwargs))**2
error_sum += (y2 - coupled_2(x_val, *params, **kwargs))**2
return np.sqrt(error_sum)
# Initial guess for parameters
initial_guess = [1] + [0]*nParams
# Optimize parameters using least squares method
result = scipy.optimize.least_squares(res, initial_guess, kwargs={'H1': H1, 'H2': H2})
return result['x']
# Obtain parameters
params = getParameters(Ha, Hc, points1, points2)
print(params)