Positionssytem
==============
Zum Bestimmen der Position mittles Schall wird eine Lösung einer Gleichung 4. Grades verlangt.
Deswegen gibt es im nächsten Abschnitt einen Solver für Gleichungen 4. Grades.
Gleichungen 5 und höheren Grades werden mit Newtons Tangentenverfahren gelöst.
Das ist nicht immer zuverlässig.
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from math import *
class R:
def __contains__(self, other):
return True
def __eq__(self, other):
return True
R = R()
def solve(coefficients):
assert coefficients, 'need something to solve'
order = len(coefficients) - 1
name = "solve_{}".format(order)
if not name in globals():
return solve_n(coefficients)
function = globals()[name]
return function(*coefficients)
def get_coefficients(*zeros, multiplier = 1, add = 0):
l = [1]
for zero in zeros:
l = [1] + [l[i] * -zero + l[i + 1] for i in range(len(l) - 1)] + [l[-1] * -zero]
l = [a * multiplier for a in l]
l[-1] += add
return l
assert get_coefficients(5, -2) == [1, -3, -10], "1x² - 3x - 10"
def call(coefficients, x):
return sum(x ** (len(coefficients) - i - 1) * c for i, c in enumerate(coefficients))
assert call([1, -3, -10], -2) == 0
assert call([1, -3, -10], 5) == 0
assert call([1, -3, -10], 0) != 0
def approximates(a, b, epsilon = 1e-10):
return abs(a - b) < epsilon
def test_solver(*zeros, epsilon = 1e-10, debug = False):
if debug: print("------ Test Solver ------")
coefficients = get_coefficients(*zeros)
for zero in zeros:
assert approximates(call(coefficients, zero), 0, epsilon = epsilon), \
"call({}, {}) == {} == {}".format(coefficients, zero, call(coefficients, zero), 0)
zeros = list(set(zeros))
zeros.sort()
solutions = list(solve(coefficients))
solutions.sort()
differences_are_very_small = [approximates(a, b, epsilon = epsilon) for a, b in zip(zeros, solutions)]
assert len(zeros) == len(solutions), "{} == {}".format(zeros, solutions)
assert all(differences_are_very_small), "{} == {}".format(zeros, solutions)
def test_zeros(*zeros, epsilon = 1e-10, expected_number_of_zeros = 1, add = 0, debug = False):
if debug: print("------ Test Zeros ------")
c = get_coefficients(*zeros, add = add)
if debug:
print("coefficients: {}".format(c))
zeros = solve(c)
if debug:
print("zeros: {}".format(zeros))
assert len(zeros) == expected_number_of_zeros
for zero in zeros:
assert approximates(call(c, zero), 0), "f({}) == 0 <=> {} == 0".format(zero, call(c, zero), 0)
def test_polynom(*coefficients, epsilon = 1e-10):
zeros = list(solve(coefficients))
zeros.sort()
for zero in zeros:
assert approximates(call(coefficients, zero), 0, epsilon = epsilon), \
"call({}, {}) == {} == {}".format(coefficients, zero, call(coefficients, zero), 0)
def coefficients_to_string(coefficients):
if not coefficients:
return ''
if len(coefficients) == 1:
return "{}".format(coefficients[0])
s = "{} * x + {}".format(coefficients[-2], coefficients[-1])
for exponent, coefficient in enumerate(reversed(coefficients[:-2]), 2):
s = s = "{} * x**{} + ".format(coefficient, exponent) + s
return s
assert coefficients_to_string([3, 4, 5, 6]) == "3 * x**3 + 4 * x**2 + 5 * x + 6", coefficients_to_string([3, 4, 5, 6])
def plot_function(coefficients, from_x = -2, to_x = 4, samples = 2000):
xs = np.linspace(from_x, to_x, samples)
plt.plot(xs, [call(coefficients, x )for x in (s / samples * (to_x - from_x) + from_x for s in range(samples))])
return plt.title("y = f(x) = " + coefficients_to_string(coefficients))
plot_function([1, -12, 3, -8, 5])
<matplotlib.text.Text at 0x589a5f0>
def solve_0(a):
if a == 0:
return [0]
return []
test_solver()
def solve_1(a, b):
"x | ax + b = 0"
if a == 0:
return solve_0(b)
return [- b / a]
test_solver(3)
test_solver(-2)
def solve_2(a, b, c):
if a == 0:
return solve_1(b, c)
p = b / a
q = c / a
D = p * p / 4 - q
if D < 0:
return []
if D == 0:
return [-p / 2]
return [- p / 2 + D ** 0.5, - p / 2 - D ** 0.5]
test_solver(3, 2.3)
test_solver(-2, 100)
def differenciate(coefficients):
return [c * (len(coefficients) - i) for i, c in enumerate(coefficients[:-1])]
def newton(coefficients, epsilon = 1e-15, max_iterations = 1000000, start = 1):
differenciation = differenciate(coefficients)
differenciation2 = differenciate(differenciation)
zero = start
last_zero = None
for iteration in range(1, max_iterations + 1):
y = call(coefficients, zero)
if abs(y) < epsilon:
break
last_zero = zero
slope = call(differenciation, zero)
#print("slope:", slope)
# optimized newton
# http://www.tavrodir.lima-city.de/Arbeiten/Newton-%20und%20Sekantenverfahren%20Beleg.pdf
assert slope
slope2 = call(differenciation2, zero)
#print("slope2:", slope2)
#print("x:", zero)
zero -= y / slope / (1 - y * slope2 / slope / slope)
if last_zero == zero:
break
print(iteration)
if iteration == max_iterations:
raise ValueError("Could not complete newton({}, epsilon = {}, max_iterations = {}) at x = {} f(x) = {}"\
.format(coefficients, epsilon, max_iterations, zero, call(coefficients, zero)))
return zero
def test_newton(*zeros, epsilon = 1e-15):
zero = newton(get_coefficients(*zeros), epsilon = epsilon)
assert any(approximates(zero, x) for x in zeros), "{} in {}".format(zero, zeros)
#test_newton(2)
#test_newton(-5, 2)
#test_newton(10.3, -0.0001, 2)
#test_newton(-33, 100, 2, -3, epsilon = 0.000001)
test_newton(21, 22, 23)
6651
c = get_coefficients(88, 8, 888)
newton(c, start = 887.9993)
5
888.0000000000001
def solve_3_newton(p, q, r, s):
if p == 0:
return solve_2(q, r, s)
zero = newton([p, q, r, s])
a = p
b = q + p * zero
c = - s / zero
assert approximates(c - b * zero, r), "r = {} = c - bz = {} - {} = {}".format(r, c, b * zero, c - b * zero)
return [zero] + solve([a, b, c])
def curt(x):
"cubic root"
if x == 0:
return 0
sign = abs(x) / x
x = abs(x)
return (x ** (1./3.)) * sign
def solve_3(A, B, C, D):
if A == 0:
return solve_2(B, C, D)
# x³ + ax² + bx + c = 0
a = B / A
b = C / A
c = D / A
# x = z - a / 3
# z³ + pz + q = 0
p = b - a*a / 3
q = 2 * a*a*a / 27 - a * b / 3 + c
D = q*q / 4 + p*p*p / 27
if D > 0:
#print("{3} -> D > 0")
#print("D: {}".format(D))
#print("sqrt_of_D: {}".format(sqrt(D)))
#print("minus_q_half: {}".format(- q / 2))
u = curt(- q / 2 + sqrt(D))
#print("u: {}".format(u))
v = curt(- q / 2 - sqrt(D))
#print("v: {}".format(v))
z = u + v
#print("z: {}".format(z))
return [z - a/3]
elif p == 0 and q == 0:
# print("{3} -> p == 0 and q == 0")
# z = 0
return [- a/3]
# elif p == 0: # and q != 0
# # impossible case since it implies D < 0
# print("-> p == 0")
# z = (-q) ** (1.0 / 3.0)
# return [z - a/3]
elif D == 0: # and (p != 0 or q != 0)
# print("{3} -> D == 0 and (p != 0 or q != 0)")
u = curt(- q / 2)
z1 = 2 * u
z23 = - u
return [z1 - a/3, z23 - a/3]
else: # D < 0
# print("{3} -> D < 0")
#print("a: {}".format(a))
#print("b: {}".format(b))
#print("c: {}".format(c))
#print("p: {}".format(p))
#print("q: {}".format(q))
#print("sq3p: {}".format(sqrt(-3 / p)))
m = sqrt(-4/3 * p)
#print("m: {}".format(m))
n = 1 / 3 * acos(-q/2 * sqrt(-27 / p/p/p))
#print("n: {}".format(n))
z2 = -m * cos(n + pi / 3)
z1 = m * cos(n)
z3 = -m * cos(n - pi / 3)
#print("z2: {}".format(z2))
#print("z1: {}".format(z1))
#print("z3: {}".format(z3))
assert z1.imag == 0
assert z2.imag == 0
assert z3.imag == 0
return [z1 - a/3, z2 - a/3, z3 - a/3]
if 1:
print(" D > 0 expected")
test_zeros(1, 2, 3, add = 100)
print(" p == 0 and q == 0 expected")
test_solver(0,0,0)
print(" D == 0 and (p != 0 or q != 0) expected")
test_zeros(0,0,3,expected_number_of_zeros = 2)
print(" D < 0 expected")
test_solver(1,2,3)
test_solver(3, 5, -2222)
test_solver(0.0001,-10.12333,0)
test_solver(-10,-10.12333,0)
test_solver(-10,0.0001,0)
test_solver(-10,0.0001,-10.12333)
solve_3(1, 0, 1, 30)
D > 0 expected p == 0 and q == 0 expected D == 0 and (p != 0 or q != 0) expected D < 0 expected
[-2.9999999999999742]
# 1.00*x^3 + -380.00*x^2 + 22647.00*x + -369468.00
#solve_3(1, -380, 22647, -369468)
32.999996185302734 * 0.0000001 -33
-32.99999670000038
def sqrt_2(a):
if a == 0:
return 0
b = 1
last_b = 0
while last_b != b:
last_b = b
b = (b + a / b) / 2
return b
assert approximates(sqrt_2(2), sqrt(2)), "{} == {} <=> sqrt_2(2) == sqrt(2)".format(sqrt_2(2), sqrt(2))
def curt_2(a):
if a == 0:
return 0
b = a
last_b_1 = 0;
last_b_2 = 0;
while (last_b_1 != b and last_b_2 != b):
last_b_1 = b;
b = (b + a / b / b) / 2;
last_b_2 = b;
b = (b + a / b / b) / 2;
return b
assert approximates(curt_2(2), curt(2)), "{} == {} <=> sqrt_2(2) == sqrt(2)".format(curt_2(2), curt(2))
def sqrt_2(a):
if a == 0:
return 0
b = 1
for i in range(10):
b = (b + a / b) / 2
return b
sqrt_2(225.03703703703704), sqrt(225.03703703703704)
(15.001234517100151, 15.001234517100153)
def curt_3(a):
if a == 0:
return 0
b = 1
for i in range(50):
b = (b + a / b / b) / 2;
return b
curt_3(15.001234517100153), 15.001234517100153** (1/3)
(2.466279729829537, 2.466279729829537)
def solve_4(e, f, g, h, i):
if e == 0:
return solve_3(f, g, h, i)
#print("Solving: {} = 0".format(coefficients_to_string([e, f, g, h, i])))
a = f / e
b = g / e
c = h / e
d = i / e
#print("Solving: {} = 0".format(coefficients_to_string([1, a, b, c, d])))
# x = z - a/4
p = b - 3 / 8 * a*a
q = a*a*a / 8 - a * b / 2 + c
r = -3/256 * a*a*a*a + a*a * b / 16 - a * c / 4 + d
#print(-3/256 * a*a*a*a)
#print(-3/256 * a*a*a*a + a*a * b / 16 )
#print(-3/256 * a*a*a*a + a*a * b / 16 - a * c / 4)
#print(-3/256 * a*a*a*a + a*a * b / 16 - a * c / 4 + d)
#print("p: {} q: {} r: {}".format(p, q, r))
if approximates(p, 0) and approximates(q, 0):
# print("{4} p == 0 and q == 0")
if approximates(r, 0):
return [0]
elif r > 0:
return []
return (-r) ** 0.5
elif approximates(q, 0):
# print("{4} p != 0 and q == 0")
# v = z²
vs = solve([1, p, r])
xs = []
#print('vs: {}'.format(vs))
for v in vs:
if approximates(v, 0):
v = 0
elif v < 0:
continue
z1 = v ** 0.5
z2 = -z1
xs.append(z1 - a / 4)
xs.append(z2 - a / 4)
return xs
#print("{{4}} solve({})".format([8, -4*p, -8*r, 4 * p * r - q*q]))
Ps = solve([8, -4*p, -8*r, 4 * p * r - q*q])
#print("Ps: {}".format(Ps))
has_solution = False
for P in Ps:
Q_squared = 2 * P - p
R_squared = P * P - r
#print("{{4}} P: {} Q_squared: {} R_squared: {}".format(P, Q_squared, R_squared))
if approximates(Q_squared, 0):
Q_squared = 0
elif Q_squared < 0:
continue
if approximates(R_squared, 0):
R_squared = 0
elif R_squared < 0:
continue
Q1 = sqrt(Q_squared)
Q2 = -Q1
R1 = sqrt(R_squared)
R2 = -R1
for R, Q in ((R1, Q1), (R1, Q2), (R2, Q1), (R2, Q2)):
# print("{{4}} P: {} approximately 0: {} | {}".format(P, Q * R + -q / 2, Q * R - -q / 2))
if approximates(Q * R, -q / 2, epsilon = 1e-9):
# print("approximation of 0: {}".format(Q * R +q / 2))
has_solution = True
break
break
if not has_solution:
#print('{4} no solution found')
return []
zs = solve([1, Q, P + R]) + solve([1, -Q, P - R])
# x = z - a/4
return [z - a / 4 for z in zs]
test_solver(40,42,43,44)
test_solver(1,2,3,4)
test_solver(0.018466806307818034, -0.02575054181030545, -18.59653702522963, -19.600040761851144)
test_solver(0.002078698483378716, -0.013961191268707651, -18.91277956509588, -19.279199464702053)
test_solver(--10.12333, --10, -0, -0.01, epsilon = 1e-6)
test_solver(-10.12333, -10, 0, 0.0001, epsilon = 1e-6)
def solve_n(coefficients):
"does not work if newton failes"
if coefficients[0] == 0:
return solve(coefficients[1:])
zero = newton(coefficients)
if zero == 0:
assert coefficients[-1] == 0
return [0] + solve(coefficients[:-1])
less_coefficients = [coefficients[0]]
for coefficient in coefficients[1:-1]:
less_coefficients.append(coefficient + zero * less_coefficients[-1])
assert approximates(less_coefficients[-1], -coefficients[-1] / zero), \
"{} = {} = {} / z | z = {}".format(less_coefficients[-1], -coefficients[-1] / zero, -coefficients[-1], zero)
return [zero] + solve(less_coefficients)
test_solver(1, 2, 3, 4, 5)
test_solver(0, 1, 2, 3, 4, 5)
#test_solver(70, 5, 30, 1, -1, -8) # should choose another epsilon
1 1 1000000
--------------------------------------------------------------------------- ValueError Traceback (most recent call last) <ipython-input-22-a36831337ae9> in <module>() 14 15 test_solver(1, 2, 3, 4, 5) ---> 16 test_solver(0, 1, 2, 3, 4, 5) 17 #test_solver(70, 5, 30, 1, -1, -8) # should choose another epsilon <ipython-input-8-0af21d08908c> in test_solver(epsilon, debug, *zeros) 26 zeros = list(set(zeros)) 27 zeros.sort() ---> 28 solutions = list(solve(coefficients)) 29 solutions.sort() 30 differences_are_very_small = [approximates(a, b, epsilon = epsilon) for a, b in zip(zeros, solutions)] <ipython-input-7-15a417d63dc0> in solve(coefficients) 4 name = "solve_{}".format(order) 5 if not name in globals(): ----> 6 return solve_n(coefficients) 7 function = globals()[name] 8 return function(*coefficients) <ipython-input-22-a36831337ae9> in solve_n(coefficients) 11 less_coefficients.append(coefficient + zero * less_coefficients[-1]) 12 assert approximates(less_coefficients[-1], -coefficients[-1] / zero), "{} = {} = {} / z | z = {}".format(less_coefficients[-1], -coefficients[-1] / zero, -coefficients[-1], zero) ---> 13 return [zero] + solve(less_coefficients) 14 15 test_solver(1, 2, 3, 4, 5) <ipython-input-7-15a417d63dc0> in solve(coefficients) 4 name = "solve_{}".format(order) 5 if not name in globals(): ----> 6 return solve_n(coefficients) 7 function = globals()[name] 8 return function(*coefficients) <ipython-input-22-a36831337ae9> in solve_n(coefficients) 3 if coefficients[0] == 0: 4 return solve(coefficients[1:]) ----> 5 zero = newton(coefficients) 6 if zero == 0: 7 assert coefficients[-1] == 0 <ipython-input-13-359be5fcc843> in newton(coefficients, epsilon, max_iterations, start) 29 print(iteration) 30 if iteration == max_iterations: ---> 31 raise ValueError("Could not complete newton({}, epsilon = {}, max_iterations = {}) at x = {} f(x) = {}" .format(coefficients, epsilon, max_iterations, zero, call(coefficients, zero))) 32 return zero 33 ValueError: Could not complete newton([1, -14, 71, -154, 120, 0], epsilon = 1e-15, max_iterations = 1000000) at x = 2.0000000000000044 f(x) = -8.526512829121202e-14
Die Grundlage des Positionsystems ist ein Rechtwinkliges Dreieck ⊿ABC, in dem die Lautsprecher angeordnet sind.
a ⊥ b <=> ᵧ = 90°
Durch die Schalldifferenz der Lautsprecher, entstehen 3 Hyperbeln, die die Position P(x, y) genau bestimmen.
Strecke CA ∥ CX = h2 ⊥ h1 = Strecke XP ∥ BC.
def get_dt(a, b, x, y):
u = sqrt(x*x + y*y)
w = sqrt(y*y + (b - x)*(b - x))
v = sqrt((a + y)*(a + y) + x*x)
#print("{{p}} u, v, w = {}, {}, {}".format(u, v, w))
g1 = v - u
g2 = v - w
g3 = u - w
return g1, g2, g3
def check(dt1, dt2, a, b, x, y):
"=> the error"
dt3 = dt2 - dt1
g1, g2, g3 = get_dt(a, b, x, y)
# if g1 == dt1 and g2 == dt2: print('{p} g1 == dt1 and g2 == dt2')
error = (g1 - dt1)**2 + (g2 - dt2)**2 + (g3 - dt3)**2
return error
def position(dt1, dt2, a, b):
"""
dt1 = PB - CP
dt2 = BP - AP
"""
g1 = dt1
print("g1: {}".format(g1))
g2 = dt2
print("g2: {}".format(g2))
g3 = g2 - g1 # CP - AP
print("g3: {}".format(g3))
n = 4 * g1*g1
o = 4 * a*a
p = 4 * a*a*a - 4 * a * g1*g1
q = a**4 - 2 * a*a * g1*g1 + g1**4
print("n: {}".format(n))
print("p: {}".format(p))
print("q: {}".format(q))
r = 4 * g3*g3
s = 4 * b*b
t = -4 * b**3 + 4 * b * g3*g3
v = b**4 - 2 * b*b * g3*g3 + g3**4
print("r: {}".format(r))
print("s: {}".format(s))
print("t: {}".format(t))
print("v: {}".format(v))
d = o - n
e = s - r
print("d: {}".format(d))
print("e: {}".format(e))
Ps = []
if r == 0:
# print("{p} r == 0")
xs = solve([e, t, v])
for x in xs:
ys = solve([d, p, q - n * x*x])
Ps.extend([(x, y) for y in ys])
elif n == 0:
print("{p} n == 0")
print("solving {} {} {}".format(d, p, q))
ys = solve([d, p, q])
for y in ys:
print("solving {} {} {}".format(e, t, v - r * y*y))
xs = solve([e, t, v - r * y*y])
Ps.extend([(x, y) for x in xs])
else:
F = -2 * e * r
G = 2 * e * v
H = e*e / n/n
K = F / n
L = G / n
M = t*t / n
A = H * d*d + K * d + r*r
B = 2 * H * d * p + K * p
C = 2 * H * d * q + H * p*p + K * q + L * d - 2 * r * v - M * d
D = 2 * H * p * q + L * p - M * p
E = H * q*q + L * q + v*v - M * q
print("A:", A);
print("B:", B);
print("C:", C);
print("D:", D);
print("E:", E);
ys = solve([A, B, C, D, E])
# print("{{p}} ys: {} == solve({})".format(ys, [A, B, C, D, E]))
for y in ys:
xs = solve([e, t, v - r * y*y])
# print("{{p}} y: {} xs: {}".format(y, xs))
Ps.extend([(x, y) for x in xs])
# check the results
choose = lambda xy: check(dt1, dt2, a, b, xy[0], xy[1])
print("{{p}} PS: {}".format("\n{p} ".join(map(lambda xy: "{}\t{}\t{}".format(xy[0], xy[1], choose(xy)), Ps))))
x, y = min(Ps, key = choose)
return x, y
def assert_approximates(a, b, epsilon = 1e-15):
assert approximates(a, b), "{} == {} but difference is {}".format(a, b, a-b)
def test_position(x0, y0, a = 1, b = 2, debug = False):
if debug: print("------ Test x: {} y: {} ------".format(x0, y0))
dt1, dt2, dt3 = get_dt(a, b, x0, y0)
if debug: print("{{test}} dt1, dt2, dt3 = {}, {}, {}".format(dt1, dt2, dt3))
x, y = position(dt1, dt2, a, b)
if debug: print("{{test}} x0: {} y0: {} \n{{test}} x: {} y: {}".format(x0, y0, x, y))
assert approximates(x, x0, epsilon = 1e-8)
assert approximates(y, y0, epsilon = 1e-8)
if 0:
# r != 0 and n != 0
test_position(2, -1)
test_position(100, -12)
test_position(123.12312321321, -100)
test_position(-333.213, -12)
# r == 0
test_position(1, 1)
test_position(1, 0)
test_position(1, -1)
test_position(0, 0)
test_position(1, -0.5)
# n == 0
test_position(12, -0.5)
test_position(-333.213, -12, 20, 10)
g1: -0.11998715618670985 g2: -10.113697356467924 g3: -9.993710200281214 n: 0.05758767059909562 p: 31998.84824658802 q: 159988.48267315145 r: 399.49697426881914 s: 400 t: -5.030257311808782 v: 0.015814680389667046 d: 1599.942412329401 e: 0.5030257311808555 A: 184305836.10572746 B: 7589175139.439026 C: 116068814847.07921 D: 781207904655.6042 E: 1952914341350.1167 {p} PS: 245.3490205162642 -8.52681486022444 808.7020588120417 {p} -235.34902051626375 -8.52681486022444 0.12011212203363518 {p} 246.79390565496567 -8.57809686477929 808.5034277289541 {p} -236.79390565496522 -8.57809686477929 0.11517534120497079 {p} 343.2129999962256 -11.999999999915298 799.133599738705 {p} -333.21299999622516 -11.999999999915298 2.858296757282845e-23 {p} 345.24642238055645 -12.072162749735583 798.9939485374928 {p} -335.246422380556 -12.072162749735583 2.557274145726705e-05
l = [0.0332622528, -0.0665245056, -3.991614818]
solve(l)
[12.0001974401655, -10.0001974401655]
call(differenciate(l), 1)
-0.0332622528
Hier wird das Minimum der Abstandswerte bestimmt.
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from math import *
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d, Axes3D #<-- Note the capitalization!
fig = plt.figure()
ax = Axes3D(fig) #<-- Note the difference from your original code...
X = np.linspace(-5, 5, 200)
Y = np.linspace(-5, 5, 300)
X, Y = np.meshgrid(X, Y)
Z = X**2 + Y** 2
cset = ax.plot_surface(X, Y, Z,cmap=plt.cm.coolwarm, linewidth=0, antialiased=False)
ax.clabel(cset, fontsize=9, inline=1)
plt.show()
fig, ax = plt.subplots()
p = ax.pcolor(X/(2*pi), Y/(2*pi), Z, cmap=plt.cm.RdBu, vmin=abs(Z).min(), vmax=abs(Z).max())
cb = fig.colorbar(p, ax=ax)
def get_dt(a, b, x, y):
u = sqrt(x*x + y*y)
w = sqrt(y*y + (b - x)*(b - x))
v = sqrt((a + y)*(a + y) + x*x)
#print("{{p}} u, v, w = {}, {}, {}".format(u, v, w))
g1 = v - u
g2 = v - w
g3 = u - w
return g1, g2, g3
def check(dt1, dt2, a, b, x, y):
"=> the error"
dt3 = dt2 - dt1
g1, g2, g3 = get_dt(a, b, x, y)
# if g1 == dt1 and g2 == dt2: print('{p} g1 == dt1 and g2 == dt2')
error = abs(g1 - dt1) #+ (g2 - dt2)**2 + (g3 - dt3)**2
return error
a = 2
b = 1
x0 = -3
y0 = -3
dt1, dt2, dt3 = get_dt(a, b, x0, y0)
X = np.linspace(-30, 30, 200)
Y = np.linspace(-30, 30, 200)
Z = [[check(dt1, dt2, a, b, x, y) for x in X] for y in Y]
fig = plt.figure()
ax = Axes3D(fig) #<-- Note the difference from your original code...
X, Y = np.meshgrid(X, Y)
cset = ax.plot_surface(X, Y, Z,cmap=plt.cm.coolwarm, linewidth=0, antialiased=False)
plt.show()
The figure above shows the hyperbol form of the heard position
def get_dt(a, b, x, y):
u = sqrt(x*x + y*y)
w = sqrt(y*y + (b - x)**2)
#v = sqrt(y*y + (a + x)**2) # in straight line
v = sqrt((a + y)*(a + y) + x*x) # rectangular
#print("{{p}} u, v, w = {}, {}, {}".format(u, v, w))
g1 = v - u
g2 = v - w
g3 = u - w
return g1, g2, g3
def check(dt1, dt2, a, b, x, y):
"=> the error"
dt3 = dt2 - dt1
g1, g2, g3 = get_dt(a, b, x, y)
# if g1 == dt1 and g2 == dt2: print('{p} g1 == dt1 and g2 == dt2')
error = abs(g1 - dt1)**2 + abs(g2 - dt2)**2# + abs(dt3 - g3)
return (error)
def plot_dts(x0, y0, a, b, dim = 30):
dt1, dt2, dt3 = get_dt(a, b, x0, y0)
X = np.linspace(-dim, dim, 500)
Y = np.linspace(dim, -dim, 500)
Z = [[check(dt1, dt2, a, b, x, y) for x in X] for y in Y]
fig = plt.figure()
ax = Axes3D(fig) #<-- Note the difference from your original code...
X, Y = np.meshgrid(X, Y)
cset = ax.plot_surface(X, Y, Z,cmap=plt.cm.coolwarm, linewidth=0, antialiased=False)
ax.view_init(50, 200)
plt.show()
plot_dts(-3, -3, 2, 1, 10)
plot_dts(5, 5, 4, 4)
from math import sqrt
def get_dt(x, y, p):
return sqrt((x - p[0])**2 + (y - p[1]) ** 2)
def get_dts(x, y, Ps):
A, B, C, D = Ps
AP = get_dt(x, y, A)
BP = get_dt(x, y, B)
CP = get_dt(x, y, C)
DP = get_dt(x, y, D)
return (BP - CP), (BP - AP), (BP - DP)
def error_f(x):
return x * x
error_f = abs
use_4_points = True
def get_error(x, y, dts, Ps):
dts2 = get_dts(x, y, Ps)
return error_f(dts2[0] - dts[0]) + error_f(dts2[1] - dts[1]) + error_f((dts2[2] - dts[2])) * use_4_points
def plot_dts(x0, y0, Ps, dim = 30, middle = None):
dts = get_dts(x0, y0, Ps)
if middle is None:
middle = x0, y0
X = np.linspace(-dim + middle[0], dim + middle[0], 100)
Y = np.linspace(middle[1] - dim, middle[1] + dim, 100)
Z = [[get_error(x, y, dts, Ps) for x in X] for y in Y]
X, Y = np.meshgrid(X, Y)
Z = np.array(Z)
fig, ax = plt.subplots()
p = ax.pcolor(X, Y, Z, cmap=plt.cm.RdBu, vmin=min(map(min, Z)), vmax=max(map(max, Z)))
cb = fig.colorbar(p, ax=ax)
#plot_dts(0, 0.75, [(0,0), (0,1), (1,0)], 3)
use_4_points = True
plot_dts(12, -.5, [(1,0), (0,2),(0,0), (1, 2) ], 30)
The piture above shows the error function for 4 points and 3 dts. The valley on the left does decrease with distance to middle.
use_4_points = False
plot_dts(2, -1, [(1,0), (0,2),(0,0), (1, 2) ], 2)
plot_dts(0.5, 0.5, [(1,0), (0,2),(0,0), (1, 2) ],2)
This hole above seems like the only minimum. We can use a simple minimization function to get there. The valley on the left creates problems. It decreses with distance but never reaches an optimum. The error will fall until float precision does not capture it any more. we need good startingpoints to deal with this.
The following grafics suggests that for points not in the rectangle we get out of it easily.
plot_dts(12, -.5, [(1,0), (0,2),(0,0), (1, 2) ], 20)
def reflect(x, x0, m):
return [x0[i] + m * (x0[i] - x[i]) for i in range(len(x))]
def minimize(function, points, alpha = 1, gamma = 2, peta = -0.5, delta = 0.5, error = 1e-15):
# Nelder Mead in Python
# https://en.wikipedia.org/wiki/Nelder%E2%80%93Mead_method
dim = len(points[0])
x0 = None
ii = 0
for i in range(200):
ii += 1
# 1
values = [(function(point), point) for point in points]
values.sort()
values, points = zip(*values)
points = list(points)
print(points)
print('error:', values)
# 2
x0 = [sum([point[i] for point in points[:-1]]) / (len(points) - 1) for i in range(dim)]
if values[0] < error:
break
print(" ----- x0:", x0, "-----")
# 3 reflection
xr = reflect(points[-1], x0, alpha)
vr = function(xr)
print("xr:", xr, vr)
if vr < values[-2] and vr >= values[0]:
points[-1] = xr
print('#reflection')
continue
# 4 expansion
if values[0]> vr:
xe = reflect(points[-1], x0, gamma)
ve = function(xe)
print("xe:", xe, ve)
print('#expansion')
if ve < vr:
points[-1] = xe
continue
else:
points[-1] = xr
continue
# 5 contraction
assert vr >= values[-2]
xc = reflect(points[-1], x0, peta)
vc = function(xc)
print("xc:", xc, vc)
if vc < vr:
print("#contraction")
points[-1] = xc
continue
# 6 reduction
points[1:] = [reflect(points[i], points[0], delta) for i in range(1, len(points))]
print("#reduction")
continue
print('iterations:', ii)
return x0
def minimization_function(dts, Ps):
return lambda p: get_error(p[0], p[1], dts, Ps)
def minimization_function_form_xy(x0, y0, Ps):
dts = get_dts(x0, y0, Ps)
print("---- minimize ----")
return minimization_function(dts, Ps)
middle_influence_on_starting_points = 0.99
def mix(v1, v2, v3):
middle = (v1 + v2 + v3) / 3
return v1 * (1 - middle_influence_on_starting_points) + middle * middle_influence_on_starting_points
def mixPoint(P1, P2, P3):
return mix(P1[0], P2[0], P3[0]), mix(P1[1], P2[1], P3[1]),
def test_minimize(x0, y0, Ps):
starting_points = [mixPoint(Ps[0], Ps[1], Ps[2]), mixPoint(Ps[1], Ps[0], Ps[2]), mixPoint(Ps[2], Ps[1], Ps[0])]
return minimize(minimization_function_form_xy(x0, y0, Ps), starting_points)
_error_f = error_f
error_f = lambda x: x*x
#test_minimize(12, -.5, [(1,0), (0,2),(0,0), (1,2)])
test_minimize(2, -1, [(1,0), (0,2),(0,0), (1,2) ])
#test_minimize(1, 0.0000001, [(1,0), (0,2),(0,0), (1,2) ])
error_f = _error_f
---- minimize ---- [(0.33999999999999997, 0.6599999999999999), (0.32999999999999996, 0.6599999999999999), (0.32999999999999996, 0.6799999999999999)] error: (3.567557327949578, 3.597770753945049, 3.7721843050231763) ----- x0: [0.33499999999999996, 0.6599999999999999] ----- xr: [0.33999999999999997, 0.6399999999999999] 3.3994245774059744 xe: [0.345, 0.6199999999999999] 3.222502484797208 #expansion [[0.345, 0.6199999999999999], (0.33999999999999997, 0.6599999999999999), (0.32999999999999996, 0.6599999999999999)] error: (3.222502484797208, 3.567557327949578, 3.597770753945049) ----- x0: [0.34249999999999997, 0.6399999999999999] ----- xr: [0.355, 0.6199999999999999] 3.193435137478141 xe: [0.3675, 0.5999999999999999] 3.0019912306588785 #expansion [[0.3675, 0.5999999999999999], [0.345, 0.6199999999999999], (0.33999999999999997, 0.6599999999999999)] error: (3.0019912306588785, 3.222502484797208, 3.567557327949578) ----- x0: [0.35624999999999996, 0.6099999999999999] ----- xr: [0.37249999999999994, 0.5599999999999998] 2.694619416189266 xe: [0.38874999999999993, 0.5099999999999998] 2.317060859564891 #expansion [[0.38874999999999993, 0.5099999999999998], [0.3675, 0.5999999999999999], [0.345, 0.6199999999999999]] error: (2.317060859564891, 3.0019912306588785, 3.222502484797208) ----- x0: [0.37812499999999993, 0.5549999999999998] ----- xr: [0.4112499999999999, 0.48999999999999977] 2.1358550395770415 xe: [0.44437499999999985, 0.4249999999999997] 1.6944974509897435 #expansion [[0.44437499999999985, 0.4249999999999997], [0.38874999999999993, 0.5099999999999998], [0.3675, 0.5999999999999999]] error: (1.6944974509897435, 2.317060859564891, 3.0019912306588785) ----- x0: [0.4165624999999999, 0.46749999999999975] ----- xr: [0.4656249999999998, 0.33499999999999963] 1.2408766244833394 xe: [0.5146874999999997, 0.2024999999999995] 0.7209940214393626 #expansion [[0.5146874999999997, 0.2024999999999995], [0.44437499999999985, 0.4249999999999997], [0.38874999999999993, 0.5099999999999998]] error: (0.7209940214393626, 1.6944974509897435, 2.317060859564891) ----- x0: [0.4795312499999998, 0.31374999999999964] ----- xr: [0.5703124999999997, 0.1174999999999995] 0.4488740800735879 xe: [0.6610937499999996, -0.07875000000000065] 0.14945787723676793 #expansion [[0.6610937499999996, -0.07875000000000065], [0.5146874999999997, 0.2024999999999995], [0.44437499999999985, 0.4249999999999997]] error: (0.14945787723676793, 0.7209940214393626, 1.6944974509897435) ----- x0: [0.5878906249999997, 0.06187499999999943] ----- xr: [0.7314062499999995, -0.3012500000000008] 0.09710192885603265 xe: [0.8749218749999993, -0.664375000000001] 0.11706103121048016 #expansion [[0.7314062499999995, -0.3012500000000008], [0.6610937499999996, -0.07875000000000065], [0.5146874999999997, 0.2024999999999995]] error: (0.09710192885603265, 0.14945787723676793, 0.7209940214393626) ----- x0: [0.6962499999999996, -0.19000000000000072] ----- xr: [0.8778124999999994, -0.5825000000000009] 0.09626800658359042 xe: [1.0593749999999993, -0.9750000000000012] 0.12172936726622967 #expansion [[0.8778124999999994, -0.5825000000000009], [0.7314062499999995, -0.3012500000000008], [0.6610937499999996, -0.07875000000000065]] error: (0.09626800658359042, 0.09710192885603265, 0.14945787723676793) ----- x0: [0.8046093749999994, -0.44187500000000085] ----- xr: [0.9481249999999992, -0.805000000000001] 0.12223847654324427 xc: [0.7328515624999995, -0.2603125000000007] 0.0878363120685515 #contraction [[0.7328515624999995, -0.2603125000000007], [0.8778124999999994, -0.5825000000000009], [0.7314062499999995, -0.3012500000000008]] error: (0.0878363120685515, 0.09626800658359042, 0.09710192885603265) ----- x0: [0.8053320312499994, -0.4214062500000008] ----- xr: [0.8792578124999993, -0.5415625000000008] 0.08572256624024802 xe: [0.9531835937499993, -0.6617187500000008] 0.08920249877862851 #expansion [[0.8792578124999993, -0.5415625000000008], [0.7328515624999995, -0.2603125000000007], [0.8778124999999994, -0.5825000000000009]] error: (0.08572256624024802, 0.0878363120685515, 0.09626800658359042) ----- x0: [0.8060546874999994, -0.4009375000000008] ----- xr: [0.7342968749999994, -0.21937500000000065] 0.08002700501870935 xe: [0.6625390624999994, -0.037812500000000526] 0.16265776007358052 #expansion [[0.7342968749999994, -0.21937500000000065], [0.8792578124999993, -0.5415625000000008], [0.7328515624999995, -0.2603125000000007]] error: (0.08002700501870935, 0.08572256624024802, 0.0878363120685515) ----- x0: [0.8067773437499994, -0.38046875000000074] ----- xr: [0.8807031249999994, -0.5006250000000008] 0.07517661773640709 xe: [0.9546289062499993, -0.6207812500000007] 0.07957846110793954 #expansion [[0.8807031249999994, -0.5006250000000008], [0.7342968749999994, -0.21937500000000065], [0.8792578124999993, -0.5415625000000008]] error: (0.07517661773640709, 0.08002700501870935, 0.08572256624024802) ----- x0: [0.8074999999999994, -0.3600000000000007] ----- xr: [0.7357421874999995, -0.17843750000000058] 0.07431711799297423 xe: [0.6639843749999996, 0.003124999999999545] 0.18395130038836038 #expansion [[0.7357421874999995, -0.17843750000000058], [0.8807031249999994, -0.5006250000000008], [0.7342968749999994, -0.21937500000000065]] error: (0.07431711799297423, 0.07517661773640709, 0.08002700501870935) ----- x0: [0.8082226562499994, -0.3395312500000007] ----- xr: [0.8821484374999995, -0.4596875000000007] 0.06472565295904391 xe: [0.9560742187499995, -0.5798437500000007] 0.06997719802727304 #expansion [[0.8821484374999995, -0.4596875000000007], [0.7357421874999995, -0.17843750000000058], [0.8807031249999994, -0.5006250000000008]] error: (0.06472565295904391, 0.07431711799297423, 0.07517661773640709) ----- x0: [0.8089453124999995, -0.31906250000000064] ----- xr: [0.7371874999999996, -0.1375000000000005] 0.07155090672410658 #reflection [[0.8821484374999995, -0.4596875000000007], [0.7371874999999996, -0.1375000000000005], [0.7357421874999995, -0.17843750000000058]] error: (0.06472565295904391, 0.07155090672410658, 0.07431711799297423) ----- x0: [0.8096679687499995, -0.2985937500000006] ----- xr: [0.8835937499999994, -0.4187500000000006] 0.05448513775830273 xe: [0.9575195312499993, -0.5389062500000006] 0.06047601897708571 #expansion [[0.8835937499999994, -0.4187500000000006], [0.8821484374999995, -0.4596875000000007], [0.7371874999999996, -0.1375000000000005]] error: (0.05448513775830273, 0.06472565295904391, 0.07155090672410658) ----- x0: [0.8828710937499995, -0.43921875000000066] ----- xr: [1.0285546874999993, -0.7409375000000008] 0.08438616605064611 xc: [0.8100292968749996, -0.2883593750000006] 0.051192852392358386 #contraction [[0.8100292968749996, -0.2883593750000006], [0.8835937499999994, -0.4187500000000006], [0.8821484374999995, -0.4596875000000007]] error: (0.051192852392358386, 0.05448513775830273, 0.06472565295904391) ----- x0: [0.8468115234374995, -0.35355468750000063] ----- xr: [0.8114746093749996, -0.24742187500000057] 0.04220582101191696 xe: [0.7761376953124997, -0.1412890625000005] 0.04606037575818465 #expansion [[0.8114746093749996, -0.24742187500000057], [0.8100292968749996, -0.2883593750000006], [0.8835937499999994, -0.4187500000000006]] error: (0.04220582101191696, 0.051192852392358386, 0.05448513775830273) ----- x0: [0.8107519531249996, -0.2678906250000006] ----- xr: [0.7379101562499998, -0.11703125000000059] 0.07160319543248124 xc: [0.8471728515624994, -0.3433203125000006] 0.04902020635992824 #contraction [[0.8114746093749996, -0.24742187500000057], [0.8471728515624994, -0.3433203125000006], [0.8100292968749996, -0.2883593750000006]] error: (0.04220582101191696, 0.04902020635992824, 0.051192852392358386) ----- x0: [0.8293237304687495, -0.2953710937500006] ----- xr: [0.8486181640624995, -0.3023828125000006] 0.0393777385679667 xe: [0.8679125976562494, -0.3093945312500006] 0.03472651880783037 #expansion [[0.8679125976562494, -0.3093945312500006], [0.8114746093749996, -0.24742187500000057], [0.8471728515624994, -0.3433203125000006]] error: (0.03472651880783037, 0.04220582101191696, 0.04902020635992824) ----- x0: [0.8396936035156245, -0.2784082031250006] ----- xr: [0.8322143554687496, -0.21349609375000056] 0.028154388589202668 xe: [0.8247351074218746, -0.14858398437500053] 0.023146505943798155 #expansion [[0.8247351074218746, -0.14858398437500053], [0.8679125976562494, -0.3093945312500006], [0.8114746093749996, -0.24742187500000057]] error: (0.023146505943798155, 0.03472651880783037, 0.04220582101191696) ----- x0: [0.846323852539062, -0.22898925781250057] ----- xr: [0.8811730957031243, -0.21055664062500057] 0.013719063340603 xe: [0.9160223388671866, -0.19212402343750057] 0.005588226500145451 #expansion [[0.9160223388671866, -0.19212402343750057], [0.8247351074218746, -0.14858398437500053], [0.8679125976562494, -0.3093945312500006]] error: (0.005588226500145451, 0.023146505943798155, 0.03472651880783037) ----- x0: [0.8703787231445306, -0.17035400390625055] ----- xr: [0.8728448486328118, -0.0313134765625005] 0.013430573199333183 #reflection [[0.9160223388671866, -0.19212402343750057], [0.8728448486328118, -0.0313134765625005], [0.8247351074218746, -0.14858398437500053]] error: (0.005588226500145451, 0.013430573199333183, 0.023146505943798155) ----- x0: [0.8944335937499992, -0.11171875000000053] ----- xr: [0.9641320800781238, -0.07485351562500053] 0.0025461878657244814 xe: [1.0338305664062484, -0.03798828125000053] 0.015962560482479092 #expansion [[0.9641320800781238, -0.07485351562500053], [0.9160223388671866, -0.19212402343750057], [0.8728448486328118, -0.0313134765625005]] error: (0.0025461878657244814, 0.005588226500145451, 0.013430573199333183) ----- x0: [0.9400772094726553, -0.13348876953125055] ----- xr: [1.0073095703124988, -0.2356640625000006] 0.0029406831103553077 #reflection [[0.9641320800781238, -0.07485351562500053], [1.0073095703124988, -0.2356640625000006], [0.9160223388671866, -0.19212402343750057]] error: (0.0025461878657244814, 0.0029406831103553077, 0.005588226500145451) ----- x0: [0.9857208251953113, -0.15525878906250057] ----- xr: [1.055419311523436, -0.11839355468750057] 0.006192685178671064 xc: [0.950871582031249, -0.17369140625000057] 0.0013314754964706024 #contraction [[0.950871582031249, -0.17369140625000057], [0.9641320800781238, -0.07485351562500053], [1.0073095703124988, -0.2356640625000006]] error: (0.0013314754964706024, 0.0025461878657244814, 0.0029406831103553077) ----- x0: [0.9575018310546863, -0.12427246093750055] ----- xr: [0.9076940917968739, -0.012880859375000497] 0.01063419364684844 xc: [0.9824057006835926, -0.17996826171875058] 0.0008816854447736241 #contraction [[0.9824057006835926, -0.17996826171875058], [0.950871582031249, -0.17369140625000057], [0.9641320800781238, -0.07485351562500053]] error: (0.0008816854447736241, 0.0013314754964706024, 0.0025461878657244814) ----- x0: [0.9666386413574208, -0.17682983398437557] ----- xr: [0.9691452026367179, -0.2788061523437506] 0.00970120704095365 xc: [0.9653853607177723, -0.12584167480468805] 0.00023521243282375378 #contraction [[0.9653853607177723, -0.12584167480468805], [0.9824057006835926, -0.17996826171875058], [0.950871582031249, -0.17369140625000057]] error: (0.00023521243282375378, 0.0008816854447736241, 0.0013314754964706024) ----- x0: [0.9738955307006825, -0.15290496826171932] ----- xr: [0.9969194793701159, -0.13211853027343806] 0.001370831850993627 xc: [0.9623835563659657, -0.16329818725585993] 0.0005651156058341884 #contraction [[0.9653853607177723, -0.12584167480468805], [0.9623835563659657, -0.16329818725585993], [0.9824057006835926, -0.17996826171875058]] error: (0.00023521243282375378, 0.0005651156058341884, 0.0008816854447736241) ----- x0: [0.963884458541869, -0.144569931030274] ----- xr: [0.9453632164001454, -0.1091716003417974] 0.00010002573316609938 xe: [0.9268419742584219, -0.07377326965332082] 0.001336960914102282 #expansion [[0.9453632164001454, -0.1091716003417974], [0.9653853607177723, -0.12584167480468805], [0.9623835563659657, -0.16329818725585993]] error: (0.00010002573316609938, 0.00023521243282375378, 0.0005651156058341884) ----- x0: [0.9553742885589589, -0.11750663757324273] ----- xr: [0.9483650207519521, -0.07171508789062553] 0.0017732550506173023 xc: [0.9588789224624623, -0.14040241241455131] 0.0001058621761042048 #contraction [[0.9453632164001454, -0.1091716003417974], [0.9588789224624623, -0.14040241241455131], [0.9653853607177723, -0.12584167480468805]] error: (0.00010002573316609938, 0.0001058621761042048, 0.00023521243282375378) ----- x0: [0.9521210694313038, -0.12478700637817436] ----- xr: [0.9388567781448354, -0.12373233795166066] 8.476917952664252e-05 xe: [0.9255924868583669, -0.12267766952514697] 0.0005224089432561776 #expansion [[0.9388567781448354, -0.12373233795166066], [0.9453632164001454, -0.1091716003417974], [0.9588789224624623, -0.14040241241455131]] error: (8.476917952664252e-05, 0.00010002573316609938, 0.0001058621761042048) ----- x0: [0.9421099972724905, -0.11645196914672903] ----- xr: [0.9253410720825186, -0.09250152587890675] 0.0006258923776300859 xc: [0.9504944598674764, -0.1284271907806402] 9.938912407551278e-06 #contraction [[0.9504944598674764, -0.1284271907806402], [0.9388567781448354, -0.12373233795166066], [0.9453632164001454, -0.1091716003417974]] error: (9.938912407551278e-06, 8.476917952664252e-05, 0.00010002573316609938) ----- x0: [0.9446756190061558, -0.12607976436615043] ----- xr: [0.9439880216121662, -0.14298792839050345] 0.000314855778943972 xc: [0.9450194177031506, -0.11762568235397391] 1.3344122577957499e-05 #contraction [[0.9504944598674764, -0.1284271907806402], [0.9450194177031506, -0.11762568235397391], [0.9388567781448354, -0.12373233795166066]] error: (9.938912407551278e-06, 1.3344122577957499e-05, 8.476917952664252e-05) ----- x0: [0.9477569387853135, -0.12302643656730705] ----- xr: [0.9566570994257917, -0.12232053518295344] 8.256021822722524e-05 xc: [0.9433068584650744, -0.12337938725948386] 2.029051885732362e-05 #contraction [[0.9504944598674764, -0.1284271907806402], [0.9450194177031506, -0.11762568235397391], [0.9433068584650744, -0.12337938725948386]] error: (9.938912407551278e-06, 1.3344122577957499e-05, 2.029051885732362e-05) ----- x0: [0.9477569387853135, -0.12302643656730705] ----- xr: [0.9522070191055526, -0.12267348587513024] 2.163740575220506e-05 xc: [0.9455318986251939, -0.12320291191339545] 4.789101510540189e-06 #contraction [[0.9455318986251939, -0.12320291191339545], [0.9504944598674764, -0.1284271907806402], [0.9450194177031506, -0.11762568235397391]] error: (4.789101510540189e-06, 9.938912407551278e-06, 1.3344122577957499e-05) ----- x0: [0.9480131792463351, -0.12581505134701781] ----- xr: [0.9510069407895196, -0.13400442034006171] 4.6703536355394394e-05 xc: [0.9465162984747428, -0.12172036685049586] 1.3107298585212362e-06 #contraction [[0.9465162984747428, -0.12172036685049586], [0.9455318986251939, -0.12320291191339545], [0.9504944598674764, -0.1284271907806402]] error: (1.3107298585212362e-06, 4.789101510540189e-06, 9.938912407551278e-06) ----- x0: [0.9460240985499684, -0.12246163938194565] ----- xr: [0.9415537372324604, -0.11649608798325112] 2.99849042204015e-05 xc: [0.9482592792087223, -0.12544441508129292] 2.021998866328808e-06 #contraction [[0.9465162984747428, -0.12172036685049586], [0.9482592792087223, -0.12544441508129292], [0.9455318986251939, -0.12320291191339545]] error: (1.3107298585212362e-06, 2.021998866328808e-06, 4.789101510540189e-06) ----- x0: [0.9473877888417326, -0.12358239096589439] ----- xr: [0.9492436790582712, -0.12396187001839333] 1.4312675353914908e-06 #reflection [[0.9465162984747428, -0.12172036685049586], [0.9492436790582712, -0.12396187001839333], [0.9482592792087223, -0.12544441508129292]] error: (1.3107298585212362e-06, 1.4312675353914908e-06, 2.021998866328808e-06) ----- x0: [0.947879988766507, -0.12284111843444459] ----- xr: [0.9475006983242917, -0.12023782178759626] 5.171088804960252e-06 xc: [0.9480696339876147, -0.12414276675786876] 2.663552747567006e-07 #contraction [[0.9480696339876147, -0.12414276675786876], [0.9465162984747428, -0.12172036685049586], [0.9492436790582712, -0.12396187001839333]] error: (2.663552747567006e-07, 1.3107298585212362e-06, 1.4312675353914908e-06) ----- x0: [0.9472929662311788, -0.12293156680418231] ----- xr: [0.9453422534040863, -0.12190126358997129] 3.986996212918106e-06 xc: [0.948268322644725, -0.12344671841128782] 1.781976795813939e-07 #contraction [[0.948268322644725, -0.12344671841128782], [0.9480696339876147, -0.12414276675786876], [0.9465162984747428, -0.12172036685049586]] error: (1.781976795813939e-07, 2.663552747567006e-07, 1.3107298585212362e-06) ----- x0: [0.9481689783161699, -0.12379474258457829] ----- xr: [0.9498216581575969, -0.1258691183186607] 3.094025973576465e-06 xc: [0.9473426383954564, -0.12275755471753708] 1.6398049776816287e-07 #contraction [[0.9473426383954564, -0.12275755471753708], [0.948268322644725, -0.12344671841128782], [0.9480696339876147, -0.12414276675786876]] error: (1.6398049776816287e-07, 1.781976795813939e-07, 2.663552747567006e-07) ----- x0: [0.9478054805200907, -0.12310213656441245] ----- xr: [0.9475413270525668, -0.12206150637095614] 7.658755199645084e-07 xc: [0.9479375572538526, -0.12362245166114061] 3.2085520650762965e-08 #contraction [[0.9479375572538526, -0.12362245166114061], [0.9473426383954564, -0.12275755471753708], [0.948268322644725, -0.12344671841128782]] error: (3.2085520650762965e-08, 1.6398049776816287e-07, 1.781976795813939e-07) ----- x0: [0.9476400978246545, -0.12319000318933884] ----- xr: [0.947011873004584, -0.12293328796738986] 4.083719235201221e-07 xc: [0.9479542102346898, -0.12331836080031333] 2.8031353127652525e-08 #contraction [[0.9479542102346898, -0.12331836080031333], [0.9479375572538526, -0.12362245166114061], [0.9473426383954564, -0.12275755471753708]] error: (2.8031353127652525e-08, 3.2085520650762965e-08, 1.6398049776816287e-07) ----- x0: [0.9479458837442712, -0.12347040623072697] ----- xr: [0.9485491290930861, -0.12418325774391686] 3.9933880012379706e-07 xc: [0.9476442610698639, -0.12311398047413202] 2.158317282627181e-08 #contraction [[0.9476442610698639, -0.12311398047413202], [0.9479542102346898, -0.12331836080031333], [0.9479375572538526, -0.12362245166114061]] error: (2.158317282627181e-08, 2.8031353127652525e-08, 3.2085520650762965e-08) ----- x0: [0.9477992356522769, -0.12321617063722268] ----- xr: [0.9476609140507011, -0.12280988961330475] 1.232902674476558e-07 xc: [0.9478683964530648, -0.12341931114918164] 3.866170423707493e-09 #contraction [[0.9478683964530648, -0.12341931114918164], [0.9476442610698639, -0.12311398047413202], [0.9479542102346898, -0.12331836080031333]] error: (3.866170423707493e-09, 2.158317282627181e-08, 2.8031353127652525e-08) ----- x0: [0.9477563287614643, -0.12326664581165683] ----- xr: [0.9475584472882388, -0.12321493082300033] 3.681197069200203e-08 xc: [0.9478552694980771, -0.12329250330598508] 7.4161378276023025e-09 #contraction [[0.9478683964530648, -0.12341931114918164], [0.9478552694980771, -0.12329250330598508], [0.9476442610698639, -0.12311398047413202]] error: (3.866170423707493e-09, 7.4161378276023025e-09, 2.158317282627181e-08) ----- x0: [0.9478618329755709, -0.12335590722758336] ----- xr: [0.948079404881278, -0.12359783398103469] 5.2025732712747136e-08 xc: [0.9477530470227173, -0.12323494385085769] 4.183530411143511e-09 #contraction [[0.9478683964530648, -0.12341931114918164], [0.9477530470227173, -0.12323494385085769], [0.9478552694980771, -0.12329250330598508]] error: (3.866170423707493e-09, 4.183530411143511e-09, 7.4161378276023025e-09) ----- x0: [0.947810721737891, -0.12332712750001967] ----- xr: [0.947766173977705, -0.12336175169405426] 1.7714252354728633e-09 xe: [0.947721626217519, -0.12339637588808886] 1.127234122457603e-08 #expansion [[0.947766173977705, -0.12336175169405426], [0.9478683964530648, -0.12341931114918164], [0.9477530470227173, -0.12323494385085769]] error: (1.7714252354728633e-09, 3.866170423707493e-09, 4.183530411143511e-09) ----- x0: [0.9478172852153849, -0.12339053142161796] ----- xr: [0.9478815234080524, -0.12354611899237823] 1.732865052902221e-08 xc: [0.9477851661190511, -0.12331273763623782] 2.4355355149778405e-10 #contraction [[0.9477851661190511, -0.12331273763623782], [0.947766173977705, -0.12336175169405426], [0.9478683964530648, -0.12341931114918164]] error: (2.4355355149778405e-10, 1.7714252354728633e-09, 3.866170423707493e-09) ----- x0: [0.947775670048378, -0.12333724466514603] ----- xr: [0.9476829436436913, -0.12325517818111043] 7.764191586911505e-09 xc: [0.9478220332507215, -0.12337827790716384] 7.326586182320066e-10 #contraction [[0.9477851661190511, -0.12331273763623782], [0.9478220332507215, -0.12337827790716384], [0.947766173977705, -0.12336175169405426]] error: (2.4355355149778405e-10, 7.326586182320066e-10, 1.7714252354728633e-09) ----- x0: [0.9478035996848863, -0.12334550777170084] ----- xr: [0.9478410253920675, -0.1233292638493474] 2.5314471955483703e-09 xc: [0.9477848868312957, -0.12335362973287756] 3.9029951800590476e-10 #contraction [[0.9477851661190511, -0.12331273763623782], [0.9477848868312957, -0.12335362973287756], [0.9478220332507215, -0.12337827790716384]] error: (2.4355355149778405e-10, 3.9029951800590476e-10, 7.326586182320066e-10) ----- x0: [0.9477850264751734, -0.12333318368455769] ----- xr: [0.9477480196996253, -0.12328808946195154] 1.460997075996243e-09 xc: [0.9478035298629475, -0.12335573079586076] 1.3524777282385514e-10 #contraction [[0.9478035298629475, -0.12335573079586076], [0.9477851661190511, -0.12331273763623782], [0.9477848868312957, -0.12335362973287756]] error: (1.3524777282385514e-10, 2.4355355149778405e-10, 3.9029951800590476e-10) ----- x0: [0.9477943479909993, -0.12333423421604929] ----- xr: [0.9478038091507028, -0.12331483869922102] 6.001208059786977e-10 xc: [0.9477896174111475, -0.12334393197446342] 7.592841497902612e-11 #contraction [[0.9477896174111475, -0.12334393197446342], [0.9478035298629475, -0.12335573079586076], [0.9477851661190511, -0.12331273763623782]] error: (7.592841497902612e-11, 1.3524777282385514e-10, 2.4355355149778405e-10) ----- x0: [0.9477965736370475, -0.12334983138516209] ----- xr: [0.9478079811550439, -0.12338692513408636] 1.0671689820533232e-09 xc: [0.9477908698780493, -0.12333128451069995] 1.4170588853635575e-11 #contraction [[0.9477908698780493, -0.12333128451069995], [0.9477896174111475, -0.12334393197446342], [0.9478035298629475, -0.12335573079586076]] error: (1.4170588853635575e-11, 7.592841497902612e-11, 1.3524777282385514e-10) ----- x0: [0.9477902436445984, -0.12333760824258169] ----- xr: [0.9477769574262493, -0.12331948568930262] 1.9938817426710398e-10 xc: [0.9477968867537729, -0.12334666951922123] 3.810668682245979e-11 #contraction [[0.9477908698780493, -0.12333128451069995], [0.9477968867537729, -0.12334666951922123], [0.9477896174111475, -0.12334393197446342]] error: (1.4170588853635575e-11, 3.810668682245979e-11, 7.592841497902612e-11) ----- x0: [0.9477938783159111, -0.12333897701496059] ----- xr: [0.9477981392206747, -0.12333402205545776] 3.371842908474786e-11 #reflection [[0.9477908698780493, -0.12333128451069995], [0.9477981392206747, -0.12333402205545776], [0.9477968867537729, -0.12334666951922123]] error: (1.4170588853635575e-11, 3.371842908474786e-11, 3.810668682245979e-11) ----- x0: [0.947794504549362, -0.12333265328307885] ----- xr: [0.947792122344951, -0.12331863704693648] 1.8160497282002648e-10 xc: [0.9477956956515674, -0.12333966140115005] 2.4684073843588673e-12 #contraction [[0.9477956956515674, -0.12333966140115005], [0.9477908698780493, -0.12333128451069995], [0.9477981392206747, -0.12333402205545776]] error: (2.4684073843588673e-12, 1.4170588853635575e-11, 3.371842908474786e-11) ----- x0: [0.9477932827648083, -0.12333547295592501] ----- xr: [0.947788426308942, -0.12333692385639225] 3.0486084803828616e-11 xc: [0.9477957109927415, -0.12333474750569139] 9.752816051518152e-12 #contraction [[0.9477956956515674, -0.12333966140115005], [0.9477957109927415, -0.12333474750569139], [0.9477908698780493, -0.12333128451069995]] error: (2.4684073843588673e-12, 9.752816051518152e-12, 1.4170588853635575e-11) ----- x0: [0.9477957033221545, -0.12333720445342072] ----- xr: [0.9478005367662596, -0.12334312439614149] 2.6214895141661657e-11 xc: [0.9477932866001019, -0.12333424448206035] 3.822494180290415e-12 #contraction [[0.9477956956515674, -0.12333966140115005], [0.9477932866001019, -0.12333424448206035], [0.9477957109927415, -0.12333474750569139]] error: (2.4684073843588673e-12, 3.822494180290415e-12, 9.752816051518152e-12) ----- x0: [0.9477944911258347, -0.1233369529416052] ----- xr: [0.9477932712589279, -0.123339158377519] 4.62221298468443e-12 xc: [0.9477951010592881, -0.12333585022364829] 3.2637415544279835e-12 #contraction [[0.9477956956515674, -0.12333966140115005], [0.9477951010592881, -0.12333585022364829], [0.9477932866001019, -0.12333424448206035]] error: (2.4684073843588673e-12, 3.2637415544279835e-12, 3.822494180290415e-12) ----- x0: [0.9477953983554277, -0.12333775581239917] ----- xr: [0.9477975101107535, -0.12334126714273799] 8.432669494668112e-12 xc: [0.9477943424777648, -0.12333600014722976] 1.183865709625325e-12 #contraction [[0.9477943424777648, -0.12333600014722976], [0.9477956956515674, -0.12333966140115005], [0.9477951010592881, -0.12333585022364829]] error: (1.183865709625325e-12, 2.4684073843588673e-12, 3.2637415544279835e-12) ----- x0: [0.9477950190646661, -0.1233378307741899] ----- xr: [0.9477949370700441, -0.12333981132473151] 2.8108506408508435e-12 xc: [0.9477950600619771, -0.1233368404989191] 1.2235527852869303e-12 #contraction [[0.9477943424777648, -0.12333600014722976], [0.9477950600619771, -0.1233368404989191], [0.9477956956515674, -0.12333966140115005]] error: (1.183865709625325e-12, 1.2235527852869303e-12, 2.4684073843588673e-12) ----- x0: [0.947794701269871, -0.12333642032307443] ----- xr: [0.9477937068881745, -0.12333317924499881] 8.617469748954597e-12 xc: [0.9477951984607191, -0.12333804086211224] 6.804812089619629e-13 #contraction [[0.9477951984607191, -0.12333804086211224], [0.9477943424777648, -0.12333600014722976], [0.9477950600619771, -0.1233368404989191]] error: (6.804812089619629e-13, 1.183865709625325e-12, 1.2235527852869303e-12) ----- x0: [0.947794770469242, -0.12333702050467099] ----- xr: [0.9477944808765069, -0.12333720051042288] 1.1374314805336478e-13 xe: [0.9477941912837717, -0.12333738051617477] 9.564553125849579e-15 #expansion [[0.9477941912837717, -0.12333738051617477], [0.9477951984607191, -0.12333804086211224], [0.9477943424777648, -0.12333600014722976]] error: (9.564553125849579e-15, 6.804812089619629e-13, 1.183865709625325e-12) ----- x0: [0.9477946948722454, -0.12333771068914351] ----- xr: [0.947795047266726, -0.12333942123105726] 1.9071653455676448e-12 xc: [0.9477945186750052, -0.12333685541818663] 3.439870877160463e-13 #contraction [[0.9477941912837717, -0.12333738051617477], [0.9477945186750052, -0.12333685541818663], [0.9477951984607191, -0.12333804086211224]] error: (9.564553125849579e-15, 3.439870877160463e-13, 6.804812089619629e-13) ----- x0: [0.9477943549793885, -0.12333711796718069] ----- xr: [0.9477935114980578, -0.12333619507224915] 4.4854386841214583e-13 xc: [0.9477947767200539, -0.12333757941464646] 2.5035267748156887e-13 #contraction [[0.9477941912837717, -0.12333738051617477], [0.9477947767200539, -0.12333757941464646], [0.9477945186750052, -0.12333685541818663]] error: (9.564553125849579e-15, 2.5035267748156887e-13, 3.439870877160463e-13) ----- x0: [0.9477944840019128, -0.12333747996541061] ----- xr: [0.9477944493288204, -0.12333810451263459] 3.1263800239102516e-13 xc: [0.947794501338459, -0.12333716769179862] 1.4030867200641002e-13 #contraction [[0.9477941912837717, -0.12333738051617477], [0.947794501338459, -0.12333716769179862], [0.9477947767200539, -0.12333757941464646]] error: (9.564553125849579e-15, 1.4030867200641002e-13, 2.5035267748156887e-13) ----- x0: [0.9477943463111154, -0.12333727410398669] ----- xr: [0.9477939159021769, -0.12333696879332692] 4.064608890586861e-14 #reflection [[0.9477941912837717, -0.12333738051617477], [0.9477939159021769, -0.12333696879332692], [0.947794501338459, -0.12333716769179862]] error: (9.564553125849579e-15, 4.064608890586861e-14, 1.4030867200641002e-13) ----- x0: [0.9477940535929743, -0.12333717465475084] ----- xr: [0.9477936058474896, -0.12333718161770307] 2.7092035817505186e-13 xc: [0.9477942774657167, -0.12333717117327472] 2.3878099250732096e-14 #contraction [[0.9477941912837717, -0.12333738051617477], [0.9477942774657167, -0.12333717117327472], [0.9477939159021769, -0.12333696879332692]] error: (9.564553125849579e-15, 2.3878099250732096e-14, 4.064608890586861e-14) ----- x0: [0.9477942343747442, -0.12333727584472474] ----- xr: [0.9477945528473115, -0.12333758289612257] 9.989754828103379e-14 xc: [0.9477940751384606, -0.12333712231902583] 5.742723938490118e-15 #contraction [[0.9477940751384606, -0.12333712231902583], [0.9477941912837717, -0.12333738051617477], [0.9477942774657167, -0.12333717117327472]] error: (5.742723938490118e-15, 9.564553125849579e-15, 2.3878099250732096e-14) ----- x0: [0.9477941332111162, -0.1233372514176003] ----- xr: [0.9477939889565157, -0.12333733166192587] 4.833000190070062e-14 xc: [0.9477942053384164, -0.12333721129543751] 3.715522082636485e-15 #contraction [[0.9477942053384164, -0.12333721129543751], [0.9477940751384606, -0.12333712231902583], [0.9477941912837717, -0.12333738051617477]] error: (3.715522082636485e-15, 5.742723938490118e-15, 9.564553125849579e-15) ----- x0: [0.9477941402384384, -0.12333716680723167] ----- xr: [0.9477940891931051, -0.12333695309828857] 3.556123521039264e-14 xc: [0.9477941657611051, -0.12333727366170322] 6.866967183731575e-16 #contraction [[0.9477941657611051, -0.12333727366170322], [0.9477942053384164, -0.12333721129543751], [0.9477940751384606, -0.12333712231902583]] error: (6.866967183731575e-16, 3.715522082636485e-15, 5.742723938490118e-15) iterations: 74