2

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

1 respuesta 1

3

numpy.empty recibe tres parámetros:

  • shape: especifica la "forma" de la matriz, sus dimensiones. Es un parámetro obligatorio y puede ser un entero o un iterable que contenga enteros y que permita indexado (tupla, lista, otro array de NumPy, etc) si la matriz tiene más de una dimensión. Es decir:

    • numpy.empty(3) crea un array unidimensional de 3 elementos:

      >>> numpy.empty()
      array([  1.33734888e+165,   9.04783241e-043,   7.11699256e-091])
      
    • numpy.empty((2, 3)) crea una matriz de dos filas y tres columnas:

      >>> numpy.empty((2, 3))
      array([[  1.22383476e-307,   8.01063937e-307,   8.90071135e-308],
             [  8.06612616e-308,   1.16820876e-307,   7.56590980e-308]])
      
    • numpy.empty((2, 3), dtype='i') crea una matriz de dos filas y tres columnas pero con datos de tipo entero con signo de 32 bits:

      >>> numpy.empty((2, 3), dtype='i')
      array([[ 909403705, 1647786293,  808531757],
             [ 758592052,  892679218, 1633771875]], dtype=int32)
      

    Así podemos seguir, por ejemplo numpy.empty((10, 10, 4)) crea una array tridimencional, típico para almacenar imágenes RGBA.

  • dtype : especifica el tipo de los datos que va a almacenar el array. Para más información sobre los tipos soportados y como especificarlos podemos ver la documentación: Data type objects.

  • order: puede tener dos valores ‘C’, para almacenar los arrays en el estilo de C (row-major order) o ‘F’ para usar el estilo de Fortran (column-major order). Es también opcional, por defecto se sigue el estilo de C. Un array se caracteriza por tener todos sus elementos en posiciones contiguas en memoria por lo que la diferencia principal estriba en el orden en el que se colocan los elementos del array cuando este no es unidimensional:

    introducir la descripción de la imagen aquí

    Esto tiene algunas implicaciones en cuanto a eficiencia o conveniencia de uno sobre otro en casos muy concretos, pero este es otro tema fuera del alcance de esta pregunta.

La salida es un array NumPy con las dimensiones y tipo especificadas pero con los elementos sin inicializar (contendrán datos arbitrarios, basura resultado de lo que otros programas o el SO escribieran anteriormente en esas direcciones de memoria) a diferencia de lo que hace numpy.zeros (inicializa todos los elementos a cero), numpy.ones (inicializa todos los elementos a uno) o numpy.full(inicializa todos los elementos a un valor pasado como argumento. Es decir, se limita a reservar memoria suficiente para contener todos los elementos del array, en función del tipo de dato.

Solo si el tipo es un objeto (por ejemplo objetos Python como str, list, etc) se inicializan los elementos a None. Es importante tenerlo en cuenta ya que antes de operar con los valores hay que asegurarse de haberlos inicializado a un valor válido, de lo contrario los resultados serán indeterminados.

Tu error se debe a que estas pasando las dimensiones de la matriz como argumentos independientes y no como primer parámetro, por lo que pasas el número de columnas de la matriz como segundo parámetro, que es dtype, cuando efectivamente debes hacer:

d_diffs_R_c = np.empty((leng_c, size_x))

dónde (leng_c, size_x) es una tupla simplemente, equivalente a hacer:

s = (leng_c, size_x)
d_diffs_R_c = np.empty(s)

Lo anterior crea una matriz de leng_c x size_x con dtype por defecto (float64) y con row-major order con los elementos sin inicializar. Supongo que tratas de hacer un ajuste del círculo por mínimos cuadrados, por tanto necesitas una matriz de dos filas (una para las coordenadas x y otra para las coordenada y) y tantas columnas como puntos tienes (size_x).

1
  • Efectivamente, lo voy a usar para ajustar una nube de datos a una circunferencia.
    – Zhisi
    Commented el 23 may. 2018 a las 12:17

Tu Respuesta

By clicking “Publica tu respuesta”, you agree to our terms of service and acknowledge you have read our privacy policy.

¿No es la respuesta que buscas? Examina otras preguntas con la etiqueta o formula tu propia pregunta.