I'm trying to make a code to solve a overdetermined system of non-linear equations.
It is to be used to locate positions of instruments on sea bottom: for each instrument there are several surface sources with known coordinates (Xs, Ys, Zs) and the time the signal propagated to the instrument Ts, so I want to find location of the instrument (Xi, Yi, Zi) and water velocity (Vw).
I have a large number of points for which I need to run it, but I'm stuck at modifying the code from single to multiple points. Basically I want to be able to set equation coefficients from the array and then loop solution over it.
Any advice will be appreciated.
here is my code for single instrument and 6 hard-coded sources:
import numpy as np
import scipy.optimize
# observed values
Ts = np.array([0.0025, .002, .0039, .0041, .0036, .0015])
def f6(p): # 6 equations to solve with hard coded sources coordinates (Xs, Ys, Zs=0) ; p[0:2] - Xi, Yi, Zi, p[3] - Vw
return np.array([
((p[0] - 0) ** 2 + (p[1] - 0) ** 2 + (p[2] - 0) ** 2) / p[3]**2,
((p[0] - 0) ** 2 + (p[1] - 5) ** 2 + (p[2] - 0) ** 2) / p[3]**2,
((p[0] - 6) ** 2 + (p[1] - 5) ** 2 + (p[2] - 0) ** 2) / p[3]**2,
((p[0] - 6) ** 2 + (p[1] - 0) ** 2 + (p[2] - 0) ** 2) / p[3]**2,
((p[0] - 6) ** 2 + (p[1] - 3) ** 2 + (p[2] - 0) ** 2) / p[3]**2,
((p[0] - 0) ** 2 + (p[1] - 3) ** 2 + (p[2] - 0) ** 2) / p[3]**2])
def system(p): # Returns the residuals
return f6(p) - Ts**2
res = scipy.optimize.least_squares(system, [1, 1, 1, 1450],
jac='2-point', loss='soft_l1', gtol=5e-16, verbose=2)
This code works as expected. But when I try to move to parameter definition of system of equations is fails.
import numpy as np
import scipy.optimize
# observed values
Ts = np.array([0.0025, .002, .0039, .0041, .0036, .0015])
# array with coefficients for equations scoor[0] - Xs, scoor[1] - Ys, Zs=0
scoor = np.array([[0, 0, 6, 6, 6, 0],
[0, 5, 5, 0, 3, 3]])
def f6s(p): # making 6 equations to solve, using scoor coefficients; p[0:2] - Xi, Yi, Zi, p[3] - Vw
return np.asarray([((p[0] - scoor[0]) ** 2 + (p[1] - scoor[1]) ** 2 + (p[2] - 0) ** 2) / p[3]**2])
def system1(p): Returns the residuals
return f6s(p) - Ts**2
res2 = scipy.optimize.least_squares(system1, [1, 1, 1, 1450],
jac='2-point', loss='soft_l1', gtol=5e-16, verbose=2)
This does not work (( and I can't figure out what is wrong ... here is error message:
Traceback (most recent call last):
File "C:\ProgramData\Anaconda3\envs\tomo\lib\site-packages\IPython\core\interactiveshell.py", line 3444, in run_code
exec(code_obj, self.user_global_ns, self.user_ns)
File "<ipython-input-19-979c61a9e4d5>", line 29, in <module>
res2 = scipy.optimize.least_squares(system1, [1, 1, 1, 1450],
File "C:\ProgramData\Anaconda3\envs\tomo\lib\site-packages\scipy\optimize\_lsq\least_squares.py", line 823, in least_squares
raise ValueError("`fun` must return at most 1-d array_like. "
ValueError: `fun` must return at most 1-d array_like. f0.shape: (1, 6)
What do I have wrong?
SOLVED: definition of equations in def f6s() should be flatted.
def f6s(p): # making 6 equations to solve, using scoor coefficients; p[0:2] - Xi, Yi, Zi, p[3] - Vw
return np.asarray([((p[0] - scoor[0]) ** 2 + (p[1] - scoor[1]) ** 2 + (p[2] - 0) ** 2) / p[3]**2]).flatten()