# Source code for welltestpy.tools.trilib

```"""welltestpy subpackage providing routines for triangulation."""
# pylint: disable=C0103
from copy import deepcopy as dcopy

import numpy as np

__all__ = ["triangulate", "sym"]

[docs]def triangulate(distances, prec, all_pos=False):
"""Triangulate points by given distances.

try to triangulate points by given distances within a symmetric matrix
'distances' with ``distances[i,j] = |pi-pj|``

thereby ``p0`` will be set to the origin ``(0,0)`` and ``p1`` to
``(|p0-p1|,0)``

Parameters
----------
distances : :class:`numpy.ndarray`
Given distances among the point to be triangulated.
It hat to be a symmetric matrix with a vanishing diagonal and

``distances[i,j] = |pi-pj|``

If a distance is unknown, you can set it to ``-1``.
prec : :class:`float`
Given Precision to be used within the algorithm. This can be used to
smooth away measure errors
all_pos : :class:`bool`, optional
If `True` all possible constellations will be calculated. Otherwise,
the first possibility will be returned.
Default: False
"""
if not _distvalid(distances, prec / 3.0):
raise ValueError("Given distances are not valid")

pntscount = np.shape(distances)

res = []

# try to triangulate with all posible starting-point constellations
for n in range(pntscount - 1):
for m in range(n + 1, pntscount):
print("")
print("Startingconstelation {} {}".format(n, m))
tmpres = _triangulatesgl(distances, n, m, prec)

for i in tmpres:
res.append(i)

if res and not all_pos:
break
if res and not all_pos:
break

if res == []:
print("no possible or unique constellation")
return []

res = _rotate(res)

print("number of overall results: {}".format(len(res)))
sol = [res]

for i in range(1, len(res)):
for j in sol:
addable &= not _solequal(res[i], j, prec)
sol.append(res[i])

return sol

def _triangulatesgl(distances, sp1, sp2, prec):
"""
Try to triangulate points.

With starting points sp1 and sp2 and a given precision.
Thereby sp1 will be at the origin (0,0) and sp2 will be at (|sp2-sp1|,0).
"""
res = []

if distances[sp1, sp2] < -0.5:
return res

pntscount = np.shape(distances)

res = [pntscount * ]

dis = distances[sp1, sp2]

res[sp1] = np.array([0.0, 0.0])
res[sp2] = np.array([dis, 0.0])

for i in range(pntscount - 2):
iterres = []
for sglres in res:
tmpres, state = _addpoints(sglres, distances, prec)

if state == 0:
for k in tmpres:
iterres.append(dcopy(k))

if iterres == []:
return []
res = dcopy(iterres)

print("number of temporal results: {}".format(len(res)))

return res

"""
Try for each point to add it to a given solution-approach.

gives all possibilities and a status about the solution:
state = 0: possibilities found
state = 1: no possibilities
state = 2: solution-approach has a contradiction with a point
"""
res = []

posfound = False

# generate all missing points in the solution approach
pleft = []
for n, m in enumerate(sol):
if np.ndim(m) == 0:
pleft.append(n)

# try each point to add to the given solution-approach
for i in pleft:
ires, state = _addpoint(sol, i, distances, prec)

# if a point is addable, add new solution-approach to the result
if state == 0:
posfound = True
for j in ires:
res.append(dcopy(j))
# if one point gives a contradiction, return empty result and state 2
elif state == 2:
return [], 2

# if no point is addable, return empty result and state 1
if posfound:
return res, 0

return res, 1

"""
Try to add point i to a given solution-approach.

gives all possibilities and a status about the solution:
state = 0: possibilities found
state = 1: no possibilities but no contradiction
state = 2: solution-approach has a contradiction with point i
"""
res = []

# if i is already part of the solution return it
if np.ndim(sol[i]) != 0:
return [sol], 0

# collect the points already present in the solution
solpnts = []
for n, m in enumerate(sol):
if np.ndim(m) != 0:
solpnts.append(n)

# number of present points
pntscount = len(solpnts)

# try to add the point in all possible constellations
for n in range(pntscount - 1):
for m in range(n + 1, pntscount):
tmppnt, state = _pntcoord(
sol, i, solpnts[n], solpnts[m], distances, prec
)

# if possiblities are found, add them (at most 2! (think about))
if state == 0:
for pnt in tmppnt:
res.append(dcopy(sol))
res[-1][i] = pnt

# if one possiblity or a contradiction is found, return the result
if state != 1:
return res, state

# if the state remaind 1, return empty result and no contradiction
return res, state

def _pntcoord(sol, i, n, m, distances, prec):
"""
Generate coordinates for point i in constellation to points m and n.

Check if these coordinates are valid with all other points in the solution.
"""
tmppnt = []

state = 1

pntscount = len(sol)

# if no distances known, return empty result and the unknown-state
if distances[i, n] < -0.5 or distances[i, m] < -0.5:
return tmppnt, state

# if the Triangle inequality is not fulfilled give a contradiction
if distances[i, n] + distances[i, m] < _dist(sol[n], sol[m]):
state = 2
return tmppnt, state

# generate the affine rotation to bring the points in the right place
g = _affinef(*_invtranmat(*_tranmat(sol[n], sol[m])))

# generate the coordinates
x = _xvalue(distances[i, n], distances[i, m], _dist(sol[n], sol[m]))
y1, y2 = _yvalue(distances[i, n], distances[i, m], _dist(sol[n], sol[m]))

# generate the possible positions
pos1 = g(np.array([x, y1]))
pos2 = g(np.array([x, y2]))

valid1 = True
valid2 = True

# check if the possible positions are valid
for k in range(pntscount):
if np.ndim(sol[k]) != 0 and distances[i, k] > -0.5:
valid1 &= abs(_dist(sol[k], pos1) - distances[i, k]) < prec
valid2 &= abs(_dist(sol[k], pos2) - distances[i, k]) < prec

# if any position is valid, add it to the result
if valid1 or valid2:
state = 0
same = abs(y1 - y2) < prec / 4.0
if valid1:
tmppnt.append(dcopy(pos1))
if valid2 and not same:
tmppnt.append(dcopy(pos2))
# if the positions are not valid, give a contradiction
else:
state = 2

return tmppnt, state

def _solequal(sol1, sol2, prec):
"""
Compare two different solutions with a given precision.

Return True if they equal.
"""
res = True

for sol_1, sol_2 in zip(sol1, sol2):
if np.ndim(sol_1) != 0 and np.ndim(sol_2) != 0:
res &= _dist(sol_1, sol_2) < prec
elif np.ndim(sol_1) != 0 and np.ndim(sol_2) == 0:
return False
elif np.ndim(sol_1) == 0 and np.ndim(sol_2) != 0:
return False

return res

def _distvalid(dis, err=0.0, verbose=True):
"""
Check if the given distances between the points are valid.

I.e. if they fulfill the triangle-equation.
"""
valid = True
valid &= np.all(dis == dis.T)
valid &= np.all(dis.diagonal() == 0.0)

pntscount = np.shape(dis)

for i in range(pntscount - 2):
for j in range(i + 1, pntscount - 1):
for k in range(j + 1, pntscount):
if dis[i, j] > -0.5 and dis[i, k] > -0.5 and dis[j, k] > -0.5:
valid &= dis[i, j] + dis[j, k] >= dis[i, k] - abs(err)
valid &= dis[i, k] + dis[k, j] >= dis[i, j] - abs(err)
valid &= dis[j, i] + dis[i, k] >= dis[j, k] - abs(err)

if verbose and not dis[i, j] + dis[j, k] >= dis[i, k]:
print("{} {} {} for {}{}".format(i, j, k, i, k))
if verbose and not dis[i, k] + dis[k, j] >= dis[i, j]:
print("{} {} {} for {}{}".format(i, j, k, i, j))
if verbose and not dis[j, i] + dis[i, k] >= dis[j, k]:
print("{} {} {} for {}{}".format(i, j, k, j, k))

return valid

def _xvalue(a, b, c):
"""
Get the x-value for the upper point of a triangle.

where c is the length of the down side starting in the origin and
lying on the x-axes, a is the distance of the unknown point to the origen
and b is the distance of the unknown point to the righter given point
"""
return (a**2 + c**2 - b**2) / (2 * c)

def _yvalue(b, a, c):
"""
Get the two possible y-values for the upper point of a triangle.

where c is the length of the down side starting in the origin and
lying on the x-axes, a is the distance of the unknown point to the origen
and b is the distance of the unknown point to the righter given point
"""
# check flatness to eliminate numerical errors when the triangle is flat
if a + b <= c or a + c <= b or b + c <= a:
return 0.0, -0.0

res = 2 * ((a * b) ** 2 + (a * c) ** 2 + (b * c) ** 2)
res -= a**4 + b**4 + c**4
# in case of numerical errors set res to 0 (hope you check validity before)
res = max(res, 0.0)
res = np.sqrt(res)
res /= 2 * c
return res, -res

def _rotate(res):
"""
Rotate all solutions in res.

So that p0 is at the origin and p1 is on the positive x-axes.
"""
rotres = dcopy(res)

for rot_i, rot_e in enumerate(rotres):
g = _affinef(*_tranmat(rot_e, rot_e))
for rot_e_j, rot_e_e in enumerate(rot_e):
rotres[rot_i][rot_e_j] = g(rot_e_e)

return rotres

def _tranmat(a, b):
"""
Get the coefficients for the affine-linear function f(x)=Ax+s.

Which fulfills that A is a rotation-matrix,
f(a) = [0,0] and f(b) = [|b-a|,0].
"""
A = np.zeros((2, 2))
A[0, 0] = b - a
A[1, 1] = b - a
A[1, 0] = -(b - a)
A[0, 1] = +(b - a)
A /= _dist(a, b)
s = -np.dot(A, a)
return A, s

def _invtranmat(A, s):
"""
Get the coefficients for the affine-linear function g(x)=Bx+t.

which is inverse to f(x)=Ax+s
"""
B = np.linalg.inv(A)
t = -np.dot(B, s)
return B, t

def _affinef(A, s):
"""Get an affine-linear function f(x) = Ax+s."""

def func(x):
"""Affine-linear function func(x) = Ax+s."""
return np.dot(A, x) + s

return func

def _affinef_pnt(a1, a2, b1, b2, prec=0.01):
"""
Get an affine-linear function that maps f(ai) = bi.

if |a2-a1| == |b2-b1| with respect to the given precision
"""
if not abs(_dist(a1, a2) - _dist(b1, b2)) < prec:
raise ValueError("Points are not in isometric relation")

func_a = _affinef(*_tranmat(a1, a2))
func_b = _affinef(*_invtranmat(*_tranmat(b1, b2)))

def func(x):
"""Affine-linear function func(ai) = bi."""
return func_b(func_a(x))

return func

def _dist(v, w):
"""Get the distance between two given point vectors v and w."""
return np.linalg.norm(np.array(v) - np.array(w))

[docs]def sym(A):
"""Get the symmetrized version of a lower or upper triangle-matrix A."""
return A + A.T - np.diag(A.diagonal())
```