I am trying to adapt this code to use in a program:
#! 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)
I have the problem in the function def Df_2b(c)
in df2b_dc = empty((len(c), x.size))
. I've been looking at the definition of empty
and I don't see two things make sense:
- The double parenthesis that there is to pass the arguments.
- The second argument is an integer that gives the number of elements in x (I don't see the point).
The error it gives me is:
d_diffs_R_c = np.empty(leng_c, size_x) TypeError: data type not understood
My version of the code is this:
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
numpy.empty
receives three parameters:shape
: specifies the "shape" of the array, its dimensions. It is a required parameter and can be an integer or an iterable containing integers and allowing indexing (tuple, list, another NumPy array, etc) if the array has more than one dimension. Namely:numpy.empty(3)
creates a one-dimensional array of 3 elements:numpy.empty((2, 3))
creates a matrix with two rows and three columns:numpy.empty((2, 3), dtype='i')
creates an array of two rows and three columns but with data of type 32-bit signed integer:So we can continue, for example
numpy.empty((10, 10, 4))
create a three-dimensional array, typical for storing RGBA images.dtype
: Specifies the type of data to be stored by the array. For more information about supported types and how to specify them, see the documentation: Data type objects .order
: can have two values‘C’
, to store the arrays in the C style ( row-major order ) or‘F’
to use the Fortran style ( column-major order ). It is also optional, by default the C style is followed. An array is characterized by having all its elements in contiguous positions in memory, so the main difference lies in the order in which the elements of the array are placed when it is not one-dimensional:This has some implications regarding the efficiency or convenience of one over the other in very specific cases, but this is another topic outside the scope of this question.
The output is a NumPy array with the specified dimensions and type but with the elements uninitialized (they will contain arbitrary data, garbage resulting from what other programs or the OS previously wrote to those memory addresses) unlike what it does
numpy.zeros
(initializes all elements to zero),numpy.ones
(initializes all elements to one) ornumpy.full
(initializes all elements to a value passed as an argument. That is, it limits itself to reserving enough memory to contain all the elements of the array, depending on the data type .Only if the type is an object (for example Python objects like
str
,list
, etc) are the a elements initializedNone
. It is important to take this into account because before operating with the values you have to make sure that you have initialized them to a valid value , otherwise the results will be indeterminate.Your error is because you are passing the dimensions of the array as separate arguments and not as the first parameter, so you are passing the number of columns in the array as the second parameter, which is
dtype
when you should actually do:where
(leng_c, size_x)
is simply a tuple, equivalent to doing:The above creates an array of
leng_c
xsize_x
withdtype
default ( float64 ) and row-major order with the elements uninitialized. I suppose you are trying to fit the circle by least squares, so you need a matrix with two rows (one for the coordinatesx
and one for the coordinatesy
) and as many columns as there are points (size_x
).