Estoy intentando adaptar este código para usarlo en un programa:
#! python
# == METHOD 2b ==
method_2b = "leastsq with jacobian"
def calc_R(xc, yc):
""" calculate the distance of each data points from the center (xc, yc) """
return sqrt((x-xc)**2 + (y-yc)**2)
def f_2b(c):
""" calculate the algebraic distance between the 2D points and the mean circle centered at c=(xc, yc) """
Ri = calc_R(*c)
return Ri - Ri.mean()
def Df_2b(c):
""" Jacobian of f_2b
The axis corresponding to derivatives must be coherent with the col_deriv option of leastsq"""
xc, yc = c
df2b_dc = empty((len(c), x.size))
Ri = calc_R(xc, yc)
df2b_dc[0] = (xc - x)/Ri # dR/dxc
df2b_dc[1] = (yc - y)/Ri # dR/dyc
df2b_dc = df2b_dc - df2b_dc.mean(axis=1)[:, newaxis]
return df2b_dc
center_estimate = x_m, y_m
center_2b, ier = optimize.leastsq(f_2b, center_estimate, Dfun=Df_2b, col_deriv=True)
xc_2b, yc_2b = center_2b
Ri_2b = calc_R(*center_2b)
R_2b = Ri_2b.mean()
residu_2b = sum((Ri_2b - R_2b)**2)
El problema lo tengo en la función def Df_2b(c)
en df2b_dc = empty((len(c), x.size))
. He estado mirando la definición de empty
y no le encuentro sentido a dos cosas:
- El doble paréntesis que hay para pasarle los argumentos.
- El segundo argumento es un entero que da el numero de elementos que hay en x (no le veo el sentido).
El error que me da es:
d_diffs_R_c = np.empty(leng_c, size_x) TypeError: data type not understood
Mi versión del código es esta:
import numpy as np
def calc_R(x, y, xc, yc):
radii = np.sqrt((x - xc)**2 + (y - yc)**2)
return radii
def diffs_R(c, x, y):
Ri = calc_R(x, y, *c)
diffs = Ri - Ri.mean()
return diffs
def jacobian(c, x, y):
""" Jacobian of diffs_R
The axis corresponding to derivatives must be coherent with the
col_deriv option of leastsq"""
xc, yc = c
size_x = np.size(x)
leng_c = len(c)
print(size_x)
print(leng_c)
d_diffs_R_c = np.empty(leng_c, size_x)
Ri = calc_R(xc, yc)
d_diffs_R_c[0] = (xc - x) / Ri # dR/dxc
d_diffs_R_c[1] = (yc - y) / Ri # dR/dyc
d_diffs_R_c = d_diffs_R_c - d_diffs_R_c.mean(axis=1)[:, np.newaxis]
return d_diffs_R_c