I want to minimise the objective function: objective_function_2() subject to the inequality constraints 69 < T(x) < 71. The objective of the optimisation is to find the set of Fourier series parameters [c0, c1, c2, wavelength] that minimises the objective function (curve fitting problem). fx() is the function that calculates the residuals between the data and the estimated values of the Fourier series. T(x) is the function calculate_tightness() that calculates the range of the fourier series estimated at each iteration.
The problem is that I don't understand why I get this error message: ValueError: all the input arrays must have same number of dimensions, but the array at index 0 has 2 dimension(s) and the array at index 1 has 1 dimension(s)
the objective_function_2() produce an n-vector of residuals and calculate_tightness() one scalar value.
Does anyone have any idea why this is happening?
import numpy as np
from scipy.optimize import LinearConstraint, NonlinearConstraint, BFGS, minimize
ref_fld = np.array([-479.9024323 , -469.80142114, -459.70040999, -449.59939883,
-439.49838768, -429.39737652, -419.29636536, -409.19535421,
-399.09434305, -388.9933319 , -378.89232074, -368.79130958,
-358.69029843, -348.58928727, -338.48827611, -328.38726496,
-318.2862538 , -308.18524265, -298.08423149, -287.98322033,
-277.88220918, -267.78119802, -257.68018686, -247.57917571,
-237.47816455, -227.3771534 , -217.27614224, -207.17513108,
-197.07411993, -186.97310877, -176.87209762, -166.77108646,
-156.6700753 , -146.56906415, -136.46805299, -126.36704183,
-116.26603068, -106.16501952, -96.06400837, -85.96299721,
-75.86198605, -65.7609749 , -55.65996374, -45.55895258,
-35.45794143, -25.35693027, -15.25591912, -5.15490796,
4.9461032 , 15.04711435, 25.14812551, 35.24913667,
45.35014782, 55.45115898, 65.55217013, 75.65318129,
85.75419245, 95.8552036 , 105.95621476, 116.05722591,
126.15823707, 136.25924823, 146.36025938, 156.46127054,
166.5622817 , 176.66329285, 186.76430401, 196.86531516,
206.96632632, 217.06733748, 227.16834863, 237.26935979,
247.37037095, 257.4713821 , 267.57239326, 277.67340441,
287.77441557, 297.87542673, 307.97643788, 318.07744904,
328.1784602 , 338.27947135, 348.38048251, 358.48149366,
368.58250482, 378.68351598, 388.78452713, 398.88553829,
408.98654944, 419.0875606 , 429.18857176, 439.28958291,
449.39059407, 459.49160523, 469.59261638, 479.69362754,
489.79463869, 499.89564985, 509.99666101, 520.09767216])
fld = np.array([-300.41522506, -120.9280477 , -274.77413647, 494.45656622,
-44.00495929, -479.90233432, 58.55913797, -326.056248 ,
84.20018256, 443.17449743])
flr = np.array([-13.20752855, 38.56985419, 44.28484794, -51.64708478,
-10.50558888, -49.95878419, -53.88137785, -12.73304144,
-54.2792669 , -7.59544309])
def fourier_series(x, c0, c1, c2, w):
"""
Parameters
----------
x
c0
c1
c2
w
Returns
-------
"""
v = np.array(x.astype(float))
# v.fill(c0)
v = c0 + c1 * np.cos(2 * np.pi / w * x) + c2 * np.sin(2 * np.pi / w * x)
return np.rad2deg(np.arctan(v))
def calculate_splot(ref_fold_frame, popt):
return np.rad2deg(np.arctan(fourier_series(ref_fold_frame, *popt)))
def calculate_tightness(theta):
curve = calculate_splot(ref_fld, theta)
amax = np.arctan(np.deg2rad(curve.max()))
amin = np.arctan(np.deg2rad(curve.min()))
return 180 - np.rad2deg(2*np.tan((amax - amin) / 2))
def fx(theta, x, y):
# function to calculate residuals for optimisation (least squares)
return np.tan(np.deg2rad(y)) - fourier_series(x, *theta)
def objective_function_2(theta):
x = fld
y = flr
return fx(theta, x, y) + calculate_tightness(theta)
Cx = NonlinearConstraint(calculate_tightness, 69, 71, jac='2-point', hess=BFGS())
x0 = [0, 1, 1, 500]
res = minimize(objective_function_2, x0, constraints=[Cx], method='trust-constr',
options={'verbose': 1}, )
Error message
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
/tmp/ipykernel_91/1210121354.py in <module>
1 Cx = NonlinearConstraint(calculate_tightness, 69, 71, jac='2-point', hess=BFGS()) #setup_optimisation_constraints(constraints_matrix, [60, -1e-2], [71, 1e-4], linear=False)
2 x0 = theta[recovery<1][0]
----> 3 res = minimize(objective_function_2, x0, constraints=[Cx], method='trust-constr',
4
5 options={'verbose': 1}, )
/mnt/c/Users/rcha0044/LoopPhD/LOOP_ENV/env_18/lib/python3.8/site-packages/scipy/optimize/_minimize.py in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
632 constraints, callback=callback, **options)
633 elif meth == 'trust-constr':
--> 634 return _minimize_trustregion_constr(fun, x0, args, jac, hess, hessp,
635 bounds, constraints,
636 callback=callback, **options)
/mnt/c/Users/rcha0044/LoopPhD/LOOP_ENV/env_18/lib/python3.8/site-packages/scipy/optimize/_trustregion_constr/minimize_trustregion_constr.py in _minimize_trustregion_constr(fun, x0, args, grad, hess, hessp, bounds, constraints, xtol, gtol, barrier_tol, sparse_jacobian, callback, maxiter, verbose, finite_diff_rel_step, initial_constr_penalty, initial_tr_radius, initial_barrier_parameter, initial_barrier_tolerance, factorization_method, disp)
507
508 elif method == 'tr_interior_point':
--> 509 _, result = tr_interior_point(
510 objective.fun, objective.grad, lagrangian_hess,
511 n_vars, canonical.n_ineq, canonical.n_eq,
/mnt/c/Users/rcha0044/LoopPhD/LOOP_ENV/env_18/lib/python3.8/site-packages/scipy/optimize/_trustregion_constr/tr_interior_point.py in tr_interior_point(fun, grad, lagr_hess, n_vars, n_ineq, n_eq, constr, jac, x0, fun0, grad0, constr_ineq0, jac_ineq0, constr_eq0, jac_eq0, stop_criteria, enforce_feasibility, xtol, state, initial_barrier_parameter, initial_tolerance, initial_penalty, initial_trust_radius, factorization_method)
302 s0 = np.maximum(-1.5*constr_ineq0, np.ones(n_ineq))
303 # Define barrier subproblem
--> 304 subprob = BarrierSubproblem(
305 x0, s0, fun, grad, lagr_hess, n_vars, n_ineq, n_eq, constr, jac,
306 barrier_parameter, tolerance, enforce_feasibility,
/mnt/c/Users/rcha0044/LoopPhD/LOOP_ENV/env_18/lib/python3.8/site-packages/scipy/optimize/_trustregion_constr/tr_interior_point.py in __init__(self, x0, s0, fun, grad, lagr_hess, n_vars, n_ineq, n_eq, constr, jac, barrier_parameter, tolerance, enforce_feasibility, global_stop_criteria, xtol, fun0, grad0, constr_ineq0, jac_ineq0, constr_eq0, jac_eq0)
51 self.xtol = xtol
52 self.fun0 = self._compute_function(fun0, constr_ineq0, s0)
---> 53 self.grad0 = self._compute_gradient(grad0)
54 self.constr0 = self._compute_constr(constr_ineq0, constr_eq0, s0)
55 self.jac0 = self._compute_jacobian(jac_eq0, jac_ineq0, s0)
/mnt/c/Users/rcha0044/LoopPhD/LOOP_ENV/env_18/lib/python3.8/site-packages/scipy/optimize/_trustregion_constr/tr_interior_point.py in _compute_gradient(self, g)
137
138 def _compute_gradient(self, g):
--> 139 return np.hstack((g, -self.barrier_parameter*np.ones(self.n_ineq)))
140
141 def _compute_jacobian(self, J_eq, J_ineq, s):
<__array_function__ internals> in hstack(*args, **kwargs)
/mnt/c/Users/rcha0044/LoopPhD/LOOP_ENV/env_18/lib/python3.8/site-packages/numpy/core/shape_base.py in hstack(tup)
343 return _nx.concatenate(arrs, 0)
344 else:
--> 345 return _nx.concatenate(arrs, 1)
346
347
<__array_function__ internals> in concatenate(*args, **kwargs)
ValueError: all the input arrays must have same number of dimensions, but the array at index 0 has 2 dimension(s) and the array at index 1 has 1 dimension(s)
SciPy/NumPy/Python version information
1.7.1 1.21.3 sys.version_info(major=3, minor=8, micro=10, releaselevel='final', serial=0)
