Curve_fit 2D-функция, содержащая пустые массивы -> фигуры (3,3,9) и (3,1) не выровнены: 9 (dim 2)! = 3 (dim 0)

Question

Это мой первый вопрос здесь.Я надеюсь, что предоставлю достаточно информации для вас, ребята.

Я столкнулся с проблемами при использовании curve_fit с функцией, которая использует массивы numpy.Использование функции с фиксированными параметрами работает нормально.

Задача:

Я хочу найти два вектора.Что я знаю, так это как векторы вращаются вокруг друг друга, системы координат и результирующего угла между одним вектором и осью Y системы координат.

Проблема:

my_funcфункция, которую я хочу использовать, но она выдает ошибку только при помещении ее в curve_fit.

my_func2 должен быть разбивкой проблемы, которая работает.Насколько я вижу, функции в принципе одинаковы, за исключением того, что my_func2 имеет все менее и менее сложные формулы.

Это код, который я написал:

import numpy as np
from scipy.optimize import curve_fit

'''The function I want to optimize. 2D-xdata with 6 free parameters.'''
def my_func(X, hx, hy, hz, lx, ly, lz):

    # These are the two independently measured angles (2D-xdata).
    phi, alpha = X

    # y-axis of the coordinate system.
    yaxis = np.array([[0],
                      [1],
                      [0]])

    # First wanted vector h (first three parameters).
    h = np.array([[hx],
                  [hy],
                  [hz]])

    # Second wanted vector l (last three parameters).
    l = np.array([[lx],
                  [ly],
                  [lz]])

    # Projection matrix.
    Pxy = np.array([[1, 0, 0],
                    [0, 1, 0],
                    [0, 0, 0]])

    # Needed to generate the rotation matrix around the unknown vector h.
    h_norm = h / np.linalg.norm(h)
    n1, n2, n3 = h_norm[0][0], h_norm[1][0], h_norm[2][0]

    # Rotation matrix for rotation around the vector h by angle alpha.
    R_h = np.array([[n1 ** 2 * (1 - np.cos(alpha)) + np.cos(alpha), n1 * n2 * (1 - np.cos(alpha)) - n3 * np.sin(alpha), n1 * n3 * (1 - np.cos(alpha)) + n2 * np.sin(alpha)],
                    [n1 * n2 * (1 - np.cos(alpha)) + n3 * np.sin(alpha), n2 ** 2 * (1 - np.cos(alpha)) + np.cos(alpha), n2 * n3 * (1 - np.cos(alpha)) - n1 * np.sin(alpha)],
                    [n1 * n3 * (1 - np.cos(alpha)) - n2 * np.sin(alpha), n2 * n3 * (1 - np.cos(alpha)) + n1 * np.sin(alpha), n3 ** 2 * (1 - np.cos(alpha)) + np.cos(alpha)]])

    # Rotate the vector l around the vector h by angle alpha.
    l_rot = np.dot(R_h, l)

    # Rotation matrix for rotation around x-axis by angle phi.
    R_x = np.array([[1, 0, 0],
                    [0, np.cos(phi), -np.sin(phi)],
                    [0, np.sin(phi), np.cos(phi)]])

    # Rotate the vector l_rot around the x-axis by angle phi.
    l_final = np.dot(R_x, l_rot)

    # Project the vector l_final into the xy-plane.
    l_final_xy = np.dot(Pxy, l_final)

    # Get the angle between the projected vector l_final_xy and the y-axis.
    angle = np.arccos(np.vdot(l_final_xy, yaxis) / (np.linalg.norm(l_final_xy)))

    # Return angle in degree.
    return angle * 180 / np.pi


'''A simplified version of the function above with less complex formulas.'''
def my_func2(X, a1, a2, a3, b1, b2, b3):

    # Represents phi and alpha of my_func.
    x1, x2 = X

    # Represents the first wanted vector of my_func.
    va = np.array([[a1],
                   [a2],
                   [a3]])

    # Represents the second wanted vector of my_func.
    vb = np.array([[b1],
                   [b2],
                   [b3]])

    # Represents the rotation matrix of my_func. It depends on the x-data and the parameters.
    M1 = np.array([[x1 * a1, x2 * b1, 0],
                   [0, x1 * a2, x2 * b2],
                   [x2 * b3, 0, x1 * a3]])

    # Some simplified math with the wanted vectors and the generated matrix.
    v_new = np.vdot(np.dot(M1, va), vb)

    return v_new

Вот некоторые испытания.

# Some x-data: phi and alpha.
xdata = [[0, 0, 0, 30, 30, 30, 60, 60, 60],
         [0, 90, 180, 0, 90, 180, 0, 90, 180]]

# Some y-data.
ydata = [10, 11, 12, 13, 14, 15, 16, 17, 18]

# Test if my_func works as expected.
print(my_func([np.pi / 4, np.pi / 4], 1, 0, 0, 1, 1, 1))

В этой строке выводится 135.0, что правильно.Я также проверил другие значения, и результат всегда выглядит правильно.

print(curve_fit(my_func2, xdata, ydata)[0])

Эта строка печатает [-0.88635298 2.75337506 0.66050304 0.13882423 0.01404608 0.02166652].Итак, подгонка упрощенной задачи работает.

print(curve_fit(my_func, xdata, ydata)[0])

В этой строке выдается следующая ошибка:

l_rot = np.dot(R_h, l) ValueError: shapes (3,3,9) and (3,1) not aligned: 9 (dim 2) != 3 (dim 0)

Итак, последний вопрос: почему я запускаюв размерные проблемы только при использовании curve_fit и как мне обойтись?

Answer 1

Благодаря Mr.T Я заметил, что мое понимание подгонки кривой было неверным.Я думал, что он будет перебирать переданные x- и y-данные и вводить каждый набор значений в функцию.

В действительности, аппроксимация кривой вводит дырку x-data, которая представляет собой список / массив.Функция должна иметь возможность обрабатывать этот список / массив самостоятельно, как объяснил Mr.T во втором комментарии.

Чтобы решить мою проблему, я просто добавил в свою функцию цикл for, который перебирает список / массив x-data.Теперь возвращаемое значение представляет собой список значений вместо одного значения.Одно значение для каждого набора х-данных.

Я не уверен, что это лучшее решение, но рабочая программа ниже.

import numpy as np
from scipy.optimize import curve_fit

'''The function I want to optimize. 2D-xdata with 6 free parameters.'''
def my_func(X, hx, hy, hz, lx, ly, lz):

    angles = []

    for index in range(len(X)):

        # These are the two independently measured angles (2D-xdata).
        phi, alpha = X[index]

        # y-axis of the coordinate system.
        yaxis = np.array([[0],
                          [1],
                          [0]])

        # First wanted vector h (first three parameters).
        h = np.array([[hx],
                      [hy],
                      [hz]])

        # Second wanted vector l (last three parameters).
        l = np.array([[lx],
                      [ly],
                      [lz]])

        # Projection matrix.
        Pxy = np.array([[1, 0, 0],
                        [0, 1, 0],
                        [0, 0, 0]])

        # Needed to generate the rotation matrix around the unknown vector h.
        h_norm = h / np.linalg.norm(h)
        n1, n2, n3 = h_norm[0][0], h_norm[1][0], h_norm[2][0]

        # Rotation matrix for rotation around the vector h by angle alpha.
        R_h = np.array([[n1 ** 2 * (1 - np.cos(alpha)) + np.cos(alpha), n1 * n2 * (1 - np.cos(alpha)) - n3 * np.sin(alpha), n1 * n3 * (1 - np.cos(alpha)) + n2 * np.sin(alpha)],
                        [n1 * n2 * (1 - np.cos(alpha)) + n3 * np.sin(alpha), n2 ** 2 * (1 - np.cos(alpha)) + np.cos(alpha), n2 * n3 * (1 - np.cos(alpha)) - n1 * np.sin(alpha)],
                        [n1 * n3 * (1 - np.cos(alpha)) - n2 * np.sin(alpha), n2 * n3 * (1 - np.cos(alpha)) + n1 * np.sin(alpha), n3 ** 2 * (1 - np.cos(alpha)) + np.cos(alpha)]])

        # Rotate the vector l around the vector h by angle alpha.
        l_rot = np.dot(R_h, l)

        # Rotation matrix for rotation around x-axis by angle phi.
        R_x = np.array([[1, 0, 0],
                        [0, np.cos(phi), -np.sin(phi)],
                        [0, np.sin(phi), np.cos(phi)]])

        # Rotate the vector l_rot around the x-axis by angle phi.
        l_final = np.dot(R_x, l_rot)

        # Project the vector l_final into the xy-plane.
        l_final_xy = np.dot(Pxy, l_final)

        # Get the angle between the projected vector l_final_xy and the y-axis.
        angle = np.arccos(np.vdot(l_final_xy, yaxis) / (np.linalg.norm(l_final_xy)))

        angles.append(angle * 180 / np.pi)

    # Return angle in degree.
    return angles

# Some x-data: phi and alpha.
xdata = [[0, 0],
         [0, 90],
         [0, 180],
         [30, 0],
         [30, 90],
         [30, 180],
         [60, 0],
         [60, 90],
         [60, 180]]

# Some y-data.
ydata = [10, 11, 12, 13, 14, 15, 16, 17, 18]

print(curve_fit(my_func, xdata, ydata)[0])