[SciPy-User] numeric deformation when using minimize() - looking for good constrains, boundaries or just a clever way to express the problem

Wed Oct 7 05:28:42 EDT 2020

import numpy as np
from scipy.optimize import minimize

def find_ellipse_from_chunk(data):
    n = data.shape[0]  # number of points
    assert data.shape == (n, 3)
    guess_center_point = data.mean(axis=0)
    guess_a_axis_vector = np.array([0.0, 0.0, 1.0])
    guess_b_axis_vector = np.array([0.0, 1.0])
    guess_phases = np.array([0.0, 0.0])
    p0 = np.hstack([guess_center_point, guess_a_axis_vector, guess_b_axis_vector, guess_phases])

    def ellipse_func(x, data):
        center_point = x[0:3]
        a_axis_vector = x[3:6]
        b_axis_vector = np.hstack((x[6], x[7], (x[3] * x[6] + x[4] * x[7]) / (-x[5])))
        a_phase, b_phase = x[8:10]
        t = np.linspace(0, (n / 20) * 2 * np.pi, n)[np.newaxis].T  # this must be transposed, so the math adds up
        error = center_point + a_axis_vector * np.cos(t * a_phase) + b_axis_vector * np.sin(t + b_phase) - data
        error_sum = np.sum(error ** 2) + np.sum(error ** 2) * np.dot(a_axis_vector, b_axis_vector)
        return error_sum

    res = minimize(ellipse_func, p0, args=data)
    x = res.x
    result = np.hstack((x[0:8], (x[3] * x[6] + x[4] * x[7]) / (-x[5]), x[8:10]))
    return result

def main():
    data = np.array([[-4.62767933, -4.6275775, -4.62735346, -4.62719652, -4.62711625, -4.62717975,
                      -4.62723845, -4.62722407, -4.62713901, -4.62708749, -4.62703238, -4.62689101,
                      -4.62687185, -4.62694013, -4.62701082, -4.62700483, -4.62697488, -4.62686825,
                      -4.62675683, -4.62675204],
                     [-1.58625998, -1.58625039, -1.58619648, -1.58617611, -1.58620606, -1.5861833,
                      -1.5861821, -1.58619169, -1.58615814, -1.58616893, -1.58613179, -1.58615934,
                      -1.58611262, -1.58610782, -1.58613179, -1.58614017, -1.58613059, -1.58612699,
                      -1.58607428, -1.58610183],
                     [-0.96714786, -0.96713827, -0.96715984, -0.96715145, -0.96716703, -0.96712869,
                      -0.96716104, -0.96713228, -0.96719698, -0.9671838, -0.96717062, -0.96717062,
                      -0.96715744, -0.96707717, -0.96709275, -0.96706519, -0.96715026, -0.96711791,
                      -0.96713588, -0.96714786]])

    x = find_ellipse_from_chunk(data.T)
    center_point = x[0:3]
    a_axis_vector = x[3:6]
    b_axis_vector = x[6:9]
    a_phase, b_phase = x[9:11]
    print("a: " + str(a_axis_vector) + " b: " + str(b_axis_vector))

if __name__ == '__main__':
    main()

This is the algorithm that displays the weirdness.

I am trying to fit my data to ellipses. The weired form of the b_axis_vector is a mathematical attempt to make the a_axis_vector and b_axis_vector, which should represent the major and minor axis of the ellipse, stand orthogonal on each other. But apparently this form of expressing the vector has a bad effect on the optimization algorithm, because the resulting vectors for the major and minor axis are rather unbalanced:

a: [-4.93203143e-07 -2.14349197e-05  5.00000007e-01] b: [-1.77464077e-04 -4.05250071e-05 -1.91235220e-09]

always, the third dimensions of the vectors is much bigger (for vector a) or smaller (for vector b) then the other components and the resulting ellipse ends up as a line.

What are ways that are more accomodating to numerical optimization for orthogonal vectors? How could constrictions look that would help? I really would appriciate helpful pointers, since this is clearly not my core expertise!

This is a version of the code that has a simple 3d plot that is helpful to view the data and the fitted ellipse:

import numpy as np
from scipy.optimize import minimize
from matplotlib.pyplot import cm
import matplotlib.pyplot as plt
import matplotlib as mpl

def find_ellipse_from_chunk(data):
    n = data.shape[0]  # number of points
    assert data.shape == (n, 3)
    guess_center_point = data.mean(axis=0)
    guess_a_axis_vector = np.array([0.0, 0.0, 1.0])
    guess_b_axis_vector = np.array([0.0, 1.0])
    guess_phases = np.array([0.0, 0.0])
    p0 = np.hstack([guess_center_point, guess_a_axis_vector, guess_b_axis_vector, guess_phases])

    def ellipse_func(x, data):
        center_point = x[0:3]
        a_axis_vector = x[3:6]
        b_axis_vector = np.hstack((x[6], x[7], (x[3] * x[6] + x[4] * x[7]) / (-x[5])))
        a_phase, b_phase = x[8:10]
        t = np.linspace(0, (n / 20) * 2 * np.pi, n)[np.newaxis].T  # this must be transposed, so the math adds up
        error = center_point + a_axis_vector * np.cos(t * a_phase) + b_axis_vector * np.sin(t + b_phase) - data
        error_sum = np.sum(error ** 2) + np.sum(error ** 2) * np.dot(a_axis_vector, b_axis_vector)
        return error_sum

    res = minimize(ellipse_func, p0, args=data)
    x = res.x
    result = np.hstack((x[0:8], (x[3] * x[6] + x[4] * x[7]) / (-x[5]), x[8:10]))
    return result

def calculate_ellipses(ellipse, points):
    calculated_ellipse = np.zeros((points, 3))

    center_point = ellipse[0:3]
    a_axis_vector = ellipse[3:6]
    b_axis_vector = ellipse[6:9]
    a_phase, b_phase = ellipse[9:11]

    print("a: " + str(a_axis_vector) + " b: " + str(b_axis_vector))
    t = np.linspace(0, (points / 20) * 2 * np.pi, points)[np.newaxis].T
    calculated_ellipse = center_point + a_axis_vector * np.cos(
        t * a_phase) + b_axis_vector * np.sin(t + b_phase)

    return calculated_ellipse.transpose()

def plot_sensor_in_3d(data, calculated_ellipse):
    fig_3d = plt.figure(num=None, figsize=(10, 10), dpi=80, facecolor='w', edgecolor='k')
    mpl.rcParams['legend.fontsize'] = 10
    ax_3d = fig_3d.gca(projection='3d')
    title = "3D representation of magnetic field vector over 20ms \nfor "  # + self.filestem
    fig_3d.suptitle(title)
    plt.xlabel('x-direction of magnetic flux density [10uT]')
    plt.ylabel('y-direction of magnetic flux density [10uT]')
    color = iter(cm.rainbow(np.linspace(0, 1, 2)))
    c = next(color)
    ax_3d.scatter(data[0, :], data[1, :], data[2, :], s=10, color=c)
    c = next(color)
    ax_3d.plot(calculated_ellipse[0, :], calculated_ellipse[1, :], calculated_ellipse[2, :], color=c)
    plt.show()

def main():
    data = np.array([[-4.62767933, -4.6275775, -4.62735346, -4.62719652, -4.62711625, -4.62717975,
                      -4.62723845, -4.62722407, -4.62713901, -4.62708749, -4.62703238, -4.62689101,
                      -4.62687185, -4.62694013, -4.62701082, -4.62700483, -4.62697488, -4.62686825,
                      -4.62675683, -4.62675204],
                     [-1.58625998, -1.58625039, -1.58619648, -1.58617611, -1.58620606, -1.5861833,
                      -1.5861821, -1.58619169, -1.58615814, -1.58616893, -1.58613179, -1.58615934,
                      -1.58611262, -1.58610782, -1.58613179, -1.58614017, -1.58613059, -1.58612699,
                      -1.58607428, -1.58610183],
                     [-0.96714786, -0.96713827, -0.96715984, -0.96715145, -0.96716703, -0.96712869,
                      -0.96716104, -0.96713228, -0.96719698, -0.9671838, -0.96717062, -0.96717062,
                      -0.96715744, -0.96707717, -0.96709275, -0.96706519, -0.96715026, -0.96711791,
                      -0.96713588, -0.96714786]])

    x = find_ellipse_from_chunk(data.T)
    print(str(x))
    center_point = x[0:3]
    a_axis_vector = x[3:6]
    b_axis_vector = x[6:9]
    a_phase, b_phase = x[9:11]
    print("a: " + str(a_axis_vector) + " b: " + str(b_axis_vector))
    calculated_ellipse = calculate_ellipses(x, data.shape[1])
    plot_sensor_in_3d(data, calculated_ellipse)

if __name__ == '__main__':
    main()

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-user/attachments/20201007/48864c0b/attachment-0001.html>