10691143420

Committed 03 Sep 2024 09:48PM UTC coverage: 76.813% (+0.8%) from 76.01%

Build # 10691143420

Build Type

Pull #35

github

Committed by

web-flow

Commit Message

Merge 1ebfe99fb into 81a2299af

Pull Request Pull Request #35: Bump actions/download-artifact from 2 to 4.1.7 in /.github/workflows

Run Details

1928 of 2510 relevant lines covered (76.81%)

10.59 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

85.2

/src/welltestpy/tools/trilib.py

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

import numpy as np

__all__ = ["triangulate", "sym"]


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)[0]

    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[0]]

    for i in range(1, len(res)):
        addable = True
        for j in sol:
            addable &= not _solequal(res[i], j, prec)
        if addable:
            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)[0]

    res = [pntscount * [0]]

    dis = distances[sp1, sp2]

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

    for i in range(pntscount - 2):
        print("add point {}".format(i))
        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


def _addpoints(sol, distances, prec):
    """
    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


def _addpoint(sol, i, distances, prec):
    """
    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)[0]

    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[0], rot_e[1]))
        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[0] - a[0]
    A[1, 1] = b[0] - a[0]
    A[1, 0] = -(b[1] - a[1])
    A[0, 1] = +(b[1] - a[1])
    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))


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

1	"""welltestpy subpackage providing routines for triangulation."""
2	# pylint: disable=C0103
3	from copy import deepcopy as dcopy	14✔
4
5	import numpy as np	14✔
6
7	__all__ = ["triangulate", "sym"]	14✔
8
9
10	def triangulate(distances, prec, all_pos=False):	14✔
11	"""Triangulate points by given distances.
12
13	try to triangulate points by given distances within a symmetric matrix
14	'distances' with ``distances[i,j] = \|pi-pj\|``
15
16	thereby ``p0`` will be set to the origin ``(0,0)`` and ``p1`` to
17	``(\|p0-p1\|,0)``
18
19	Parameters
20	----------
21	distances : :class:`numpy.ndarray`
22	Given distances among the point to be triangulated.
23	It hat to be a symmetric matrix with a vanishing diagonal and
24
25	``distances[i,j] = \|pi-pj\|``
26
27	If a distance is unknown, you can set it to ``-1``.
28	prec : :class:`float`
29	Given Precision to be used within the algorithm. This can be used to
30	smooth away measure errors
31	all_pos : :class:`bool`, optional
32	If `True` all possible constellations will be calculated. Otherwise,
33	the first possibility will be returned.
34	Default: False
35	"""
36	if not _distvalid(distances, prec / 3.0):	14✔
37	raise ValueError("Given distances are not valid")	×
38
39	pntscount = np.shape(distances)[0]	14✔
40
41	res = []	14✔
42
43	# try to triangulate with all posible starting-point constellations
44	for n in range(pntscount - 1):	14✔
45	for m in range(n + 1, pntscount):	14✔
46	print("")	14✔
47	print("Startingconstelation {} {}".format(n, m))	14✔
48	tmpres = _triangulatesgl(distances, n, m, prec)	14✔
49
50	for i in tmpres:	14✔
51	res.append(i)	14✔
52
53	if res and not all_pos:	14✔
54	break	14✔
55	if res and not all_pos:	14✔
56	break	14✔
57
58	if res == []:	14✔
59	print("no possible or unique constellation")	×
60	return []	×
61
62	res = _rotate(res)	14✔
63
64	print("number of overall results: {}".format(len(res)))	14✔
65	sol = [res[0]]	14✔
66
67	for i in range(1, len(res)):	14✔
68	addable = True	14✔
69	for j in sol:	14✔
70	addable &= not _solequal(res[i], j, prec)	14✔
71	if addable:	14✔
72	sol.append(res[i])	14✔
73
74	return sol	14✔
75
76
77	def _triangulatesgl(distances, sp1, sp2, prec):	14✔
78	"""
79	Try to triangulate points.
80
81	With starting points sp1 and sp2 and a given precision.
82	Thereby sp1 will be at the origin (0,0) and sp2 will be at (\|sp2-sp1\|,0).
83	"""
84	res = []	14✔
85
86	if distances[sp1, sp2] < -0.5:	14✔
87	return res	×
88
89	pntscount = np.shape(distances)[0]	14✔
90
91	res = [pntscount * [0]]	14✔
92
93	dis = distances[sp1, sp2]	14✔
94
95	res[0][sp1] = np.array([0.0, 0.0])	14✔
96	res[0][sp2] = np.array([dis, 0.0])	14✔
97
98	for i in range(pntscount - 2):	14✔
99	print("add point {}".format(i))	14✔
100	iterres = []	14✔
101	for sglres in res:	14✔
102	tmpres, state = _addpoints(sglres, distances, prec)	14✔
103
104	if state == 0:	14✔
105	for k in tmpres:	14✔
106	iterres.append(dcopy(k))	14✔
107
108	if iterres == []:	14✔
109	return []	×
110	res = dcopy(iterres)	14✔
111
112	print("number of temporal results: {}".format(len(res)))	14✔
113
114	return res	14✔
115
116
117	def _addpoints(sol, distances, prec):	14✔
118	"""
119	Try for each point to add it to a given solution-approach.
120
121	gives all possibilities and a status about the solution:
122	state = 0: possibilities found
123	state = 1: no possibilities
124	state = 2: solution-approach has a contradiction with a point
125	"""
126	res = []	14✔
127
128	posfound = False	14✔
129
130	# generate all missing points in the solution approach
131	pleft = []	14✔
132	for n, m in enumerate(sol):	14✔
133	if np.ndim(m) == 0:	14✔
134	pleft.append(n)	14✔
135
136	# try each point to add to the given solution-approach
137	for i in pleft:	14✔
138	ires, state = _addpoint(sol, i, distances, prec)	14✔
139
140	# if a point is addable, add new solution-approach to the result
141	if state == 0:	14✔
142	posfound = True	14✔
143	for j in ires:	14✔
144	res.append(dcopy(j))	14✔
145	# if one point gives a contradiction, return empty result and state 2
146	elif state == 2:	×
147	return [], 2	×
148
149	# if no point is addable, return empty result and state 1
150	if posfound:	14✔
151	return res, 0	14✔
152
153	return res, 1	×
154
155
156	def _addpoint(sol, i, distances, prec):	14✔
157	"""
158	Try to add point i to a given solution-approach.
159
160	gives all possibilities and a status about the solution:
161	state = 0: possibilities found
162	state = 1: no possibilities but no contradiction
163	state = 2: solution-approach has a contradiction with point i
164	"""
165	res = []	14✔
166
167	# if i is already part of the solution return it
168	if np.ndim(sol[i]) != 0:	14✔
169	return [sol], 0	×
170
171	# collect the points already present in the solution
172	solpnts = []	14✔
173	for n, m in enumerate(sol):	14✔
174	if np.ndim(m) != 0:	14✔
175	solpnts.append(n)	14✔
176
177	# number of present points
178	pntscount = len(solpnts)	14✔
179
180	# try to add the point in all possible constellations
181	for n in range(pntscount - 1):	14✔
182	for m in range(n + 1, pntscount):	14✔
183	tmppnt, state = _pntcoord(	14✔
184	sol, i, solpnts[n], solpnts[m], distances, prec
185	)
186
187	# if possiblities are found, add them (at most 2! (think about))
188	if state == 0:	14✔
189	for pnt in tmppnt:	14✔
190	res.append(dcopy(sol))	14✔
191	res[-1][i] = pnt	14✔
192
193	# if one possiblity or a contradiction is found, return the result
194	if state != 1:	14✔
195	return res, state	14✔
196
197	# if the state remaind 1, return empty result and no contradiction
198	return res, state	×
199
200
201	def _pntcoord(sol, i, n, m, distances, prec):	14✔
202	"""
203	Generate coordinates for point i in constellation to points m and n.
204
205	Check if these coordinates are valid with all other points in the solution.
206	"""
207	tmppnt = []	14✔
208
209	state = 1	14✔
210
211	pntscount = len(sol)	14✔
212
213	# if no distances known, return empty result and the unknown-state
214	if distances[i, n] < -0.5 or distances[i, m] < -0.5:	14✔
215	return tmppnt, state	×
216
217	# if the Triangle inequality is not fulfilled give a contradiction
218	if distances[i, n] + distances[i, m] < _dist(sol[n], sol[m]):	14✔
219	state = 2	×
220	return tmppnt, state	×
221
222	# generate the affine rotation to bring the points in the right place
223	g = _affinef(_invtranmat(_tranmat(sol[n], sol[m])))	14✔
224
225	# generate the coordinates
226	x = _xvalue(distances[i, n], distances[i, m], _dist(sol[n], sol[m]))	14✔
227	y1, y2 = _yvalue(distances[i, n], distances[i, m], _dist(sol[n], sol[m]))	14✔
228
229	# generate the possible positions
230	pos1 = g(np.array([x, y1]))	14✔
231	pos2 = g(np.array([x, y2]))	14✔
232
233	valid1 = True	14✔
234	valid2 = True	14✔
235
236	# check if the possible positions are valid
237	for k in range(pntscount):	14✔
238	if np.ndim(sol[k]) != 0 and distances[i, k] > -0.5:	14✔
239	valid1 &= abs(_dist(sol[k], pos1) - distances[i, k]) < prec	14✔
240	valid2 &= abs(_dist(sol[k], pos2) - distances[i, k]) < prec	14✔
241
242	# if any position is valid, add it to the result
243	if valid1 or valid2:	14✔
244	state = 0	14✔
245	same = abs(y1 - y2) < prec / 4.0	14✔
246	if valid1:	14✔
247	tmppnt.append(dcopy(pos1))	14✔
248	if valid2 and not same:	14✔
249	tmppnt.append(dcopy(pos2))	14✔
250	# if the positions are not valid, give a contradiction
251	else:
252	state = 2	×
253
254	return tmppnt, state	14✔
255
256
257	def _solequal(sol1, sol2, prec):	14✔
258	"""
259	Compare two different solutions with a given precision.
260
261	Return True if they equal.
262	"""
263	res = True	14✔
264
265	for sol_1, sol_2 in zip(sol1, sol2):	14✔
266	if np.ndim(sol_1) != 0 and np.ndim(sol_2) != 0:	14✔
267	res &= _dist(sol_1, sol_2) < prec	14✔
268	elif np.ndim(sol_1) != 0 and np.ndim(sol_2) == 0:	×
269	return False	×
270	elif np.ndim(sol_1) == 0 and np.ndim(sol_2) != 0:	×
271	return False	×
272
273	return res	14✔
274
275
276	def _distvalid(dis, err=0.0, verbose=True):	14✔
277	"""
278	Check if the given distances between the points are valid.
279
280	I.e. if they fulfill the triangle-equation.
281	"""
282	valid = True	14✔
283	valid &= np.all(dis == dis.T)	14✔
284	valid &= np.all(dis.diagonal() == 0.0)	14✔
285
286	pntscount = np.shape(dis)[0]	14✔
287
288	for i in range(pntscount - 2):	14✔
289	for j in range(i + 1, pntscount - 1):	14✔
290	for k in range(j + 1, pntscount):	14✔
291	if dis[i, j] > -0.5 and dis[i, k] > -0.5 and dis[j, k] > -0.5:	14✔
292	valid &= dis[i, j] + dis[j, k] >= dis[i, k] - abs(err)	14✔
293	valid &= dis[i, k] + dis[k, j] >= dis[i, j] - abs(err)	14✔
294	valid &= dis[j, i] + dis[i, k] >= dis[j, k] - abs(err)	14✔
295
296	if verbose and not dis[i, j] + dis[j, k] >= dis[i, k]:	14✔
297	print("{} {} {} for {}{}".format(i, j, k, i, k))	×
298	if verbose and not dis[i, k] + dis[k, j] >= dis[i, j]:	14✔
299	print("{} {} {} for {}{}".format(i, j, k, i, j))	×
300	if verbose and not dis[j, i] + dis[i, k] >= dis[j, k]:	14✔
301	print("{} {} {} for {}{}".format(i, j, k, j, k))	×
302
303	return valid	14✔
304
305
306	def _xvalue(a, b, c):	14✔
307	"""
308	Get the x-value for the upper point of a triangle.
309
310	where c is the length of the down side starting in the origin and
311	lying on the x-axes, a is the distance of the unknown point to the origen
312	and b is the distance of the unknown point to the righter given point
313	"""
314	return (a2 + c2 - b*2) / (2 c)	14✔
315
316
317	def _yvalue(b, a, c):	14✔
318	"""
319	Get the two possible y-values for the upper point of a triangle.
320
321	where c is the length of the down side starting in the origin and
322	lying on the x-axes, a is the distance of the unknown point to the origen
323	and b is the distance of the unknown point to the righter given point
324	"""
325	# check flatness to eliminate numerical errors when the triangle is flat
326	if a + b <= c or a + c <= b or b + c <= a:	14✔
327	return 0.0, -0.0	×
328
329	res = 2 * ((a * b) ** 2 + (a * c) ** 2 + (b * c) ** 2)	14✔
330	res -= a4 + b4 + c**4	14✔
331	# in case of numerical errors set res to 0 (hope you check validity before)
332	res = max(res, 0.0)	14✔
333	res = np.sqrt(res)	14✔
334	res /= 2 * c	14✔
335	return res, -res	14✔
336
337
338	def _rotate(res):	14✔
339	"""
340	Rotate all solutions in res.
341
342	So that p0 is at the origin and p1 is on the positive x-axes.
343	"""
344	rotres = dcopy(res)	14✔
345
346	for rot_i, rot_e in enumerate(rotres):	14✔
347	g = _affinef(*_tranmat(rot_e[0], rot_e[1]))	14✔
348	for rot_e_j, rot_e_e in enumerate(rot_e):	14✔
349	rotres[rot_i][rot_e_j] = g(rot_e_e)	14✔
350
351	return rotres	14✔
352
353
354	def _tranmat(a, b):	14✔
355	"""
356	Get the coefficients for the affine-linear function f(x)=Ax+s.
357
358	Which fulfills that A is a rotation-matrix,
359	f(a) = [0,0] and f(b) = [\|b-a\|,0].
360	"""
361	A = np.zeros((2, 2))	14✔
362	A[0, 0] = b[0] - a[0]	14✔
363	A[1, 1] = b[0] - a[0]	14✔
364	A[1, 0] = -(b[1] - a[1])	14✔
365	A[0, 1] = +(b[1] - a[1])	14✔
366	A /= _dist(a, b)	14✔
367	s = -np.dot(A, a)	14✔
368	return A, s	14✔
369
370
371	def _invtranmat(A, s):	14✔
372	"""
373	Get the coefficients for the affine-linear function g(x)=Bx+t.
374
375	which is inverse to f(x)=Ax+s
376	"""
377	B = np.linalg.inv(A)	14✔
378	t = -np.dot(B, s)	14✔
379	return B, t	14✔
380
381
382	def _affinef(A, s):	14✔
383	"""Get an affine-linear function f(x) = Ax+s."""
384
385	def func(x):	14✔
386	"""Affine-linear function func(x) = Ax+s."""
387	return np.dot(A, x) + s	14✔
388
389	return func	14✔
390
391
392	def _affinef_pnt(a1, a2, b1, b2, prec=0.01):	14✔
393	"""
394	Get an affine-linear function that maps f(ai) = bi.
395
396	if \|a2-a1\| == \|b2-b1\| with respect to the given precision
397	"""
398	if not abs(_dist(a1, a2) - _dist(b1, b2)) < prec:	×
399	raise ValueError("Points are not in isometric relation")	×
400
401	func_a = _affinef(*_tranmat(a1, a2))	×
402	func_b = _affinef(_invtranmat(_tranmat(b1, b2)))	×
403
404	def func(x):	×
405	"""Affine-linear function func(ai) = bi."""
406	return func_b(func_a(x))	×
407
408	return func	×
409
410
411	def _dist(v, w):	14✔
412	"""Get the distance between two given point vectors v and w."""
413	return np.linalg.norm(np.array(v) - np.array(w))	14✔
414
415
416	def sym(A):	14✔
417	"""Get the symmetrized version of a lower or upper triangle-matrix A."""
418	return A + A.T - np.diag(A.diagonal())	14✔

GeoStat-Framework / welltestpy / 10691143420

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous