""" .. module:: linalg :platform: Unix, Windows :synopsis: Provides linear algebra utility functions .. moduleauthor:: Onur Rauf Bingol <orbingol@gmail.com> """ import os import math from copy import deepcopy from functools import reduce from .exceptions import GeomdlException from . import _linalg try: from functools import lru_cache except ImportError: from .functools_lru_cache import lru_cache def vector_cross(vector1, vector2): """ Computes the cross-product of the input vectors. :param vector1: input vector 1 :type vector1: list, tuple :param vector2: input vector 2 :type vector2: list, tuple :return: result of the cross product :rtype: tuple """ try: if vector1 is None or len(vector1) == 0 or vector2 is None or len(vector2) == 0: raise ValueError("Input vectors cannot be empty") except TypeError as e: print("An error occurred: {}".format(e.args[-1])) raise TypeError("Input must be a list or tuple") except Exception: raise if not 1 < len(vector1) <= 3 or not 1 < len(vector2) <= 3: raise ValueError("The input vectors should contain 2 or 3 elements") # Convert 2-D to 3-D, if necessary if len(vector1) == 2: v1 = [float(v) for v in vector1] + [0.0] else: v1 = vector1 if len(vector2) == 2: v2 = [float(v) for v in vector2] + [0.0] else: v2 = vector2 # Compute cross product vector_out = [(v1[1] * v2[2]) - (v1[2] * v2[1]), (v1[2] * v2[0]) - (v1[0] * v2[2]), (v1[0] * v2[1]) - (v1[1] * v2[0])] # Return the cross product of the input vectors return vector_out def vector_dot(vector1, vector2): """ Computes the dot-product of the input vectors. :param vector1: input vector 1 :type vector1: list, tuple :param vector2: input vector 2 :type vector2: list, tuple :return: result of the dot product :rtype: float """ try: if vector1 is None or len(vector1) == 0 or vector2 is None or len(vector2) == 0: raise ValueError("Input vectors cannot be empty") except TypeError as e: print("An error occurred: {}".format(e.args[-1])) raise TypeError("Input must be a list or tuple") except Exception: raise # Compute dot product prod = 0.0 for v1, v2 in zip(vector1, vector2): prod += v1 * v2 # Return the dot product of the input vectors return prod def vector_multiply(vector_in, scalar): """ Multiplies the vector with a scalar value. This operation is also called *vector scaling*. :param vector_in: vector :type vector_in: list, tuple :param scalar: scalar value :type scalar: int, float :return: updated vector :rtype: tuple """ scaled_vector = [v * scalar for v in vector_in] return scaled_vector def vector_sum(vector1, vector2, coeff=1.0): """ Sums the vectors. This function computes the result of the vector operation :math:`\\overline{v}_{1} + c * \\overline{v}_{2}`, where :math:`\\overline{v}_{1}` is ``vector1``, :math:`\\overline{v}_{2}` is ``vector2`` and :math:`c` is ``coeff``. :param vector1: vector 1 :type vector1: list, tuple :param vector2: vector 2 :type vector2: list, tuple :param coeff: multiplier for vector 2 :type coeff: float :return: updated vector :rtype: list """ summed_vector = [v1 + (coeff * v2) for v1, v2 in zip(vector1, vector2)] return summed_vector def vector_normalize(vector_in, decimals=18): """ Generates a unit vector from the input. :param vector_in: vector to be normalized :type vector_in: list, tuple :param decimals: number of significands :type decimals: int :return: the normalized vector (i.e. the unit vector) :rtype: list """ try: if vector_in is None or len(vector_in) == 0: raise ValueError("Input vector cannot be empty") except TypeError as e: print("An error occurred: {}".format(e.args[-1])) raise TypeError("Input must be a list or tuple") except Exception: raise # Calculate magnitude of the vector magnitude = vector_magnitude(vector_in) # Normalize the vector if magnitude > 0: vector_out = [] for vin in vector_in: vector_out.append(vin / magnitude) # Return the normalized vector and consider the number of significands return [float(("{:." + str(decimals) + "f}").format(vout)) for vout in vector_out] else: raise ValueError("The magnitude of the vector is zero") def vector_generate(start_pt, end_pt, normalize=False): """ Generates a vector from 2 input points. :param start_pt: start point of the vector :type start_pt: list, tuple :param end_pt: end point of the vector :type end_pt: list, tuple :param normalize: if True, the generated vector is normalized :type normalize: bool :return: a vector from start_pt to end_pt :rtype: list """ try: if start_pt is None or len(start_pt) == 0 or end_pt is None or len(end_pt) == 0: raise ValueError("Input points cannot be empty") except TypeError as e: print("An error occurred: {}".format(e.args[-1])) raise TypeError("Input must be a list or tuple") except Exception: raise ret_vec = [] for sp, ep in zip(start_pt, end_pt): ret_vec.append(ep - sp) if normalize: ret_vec = vector_normalize(ret_vec) return ret_vec def vector_mean(*args): """ Computes the mean (average) of a list of vectors. The function computes the arithmetic mean of a list of vectors, which are also organized as a list of integers or floating point numbers. .. code-block:: python :linenos: # Import geomdl.utilities module from geomdl import utilities # Create a list of vectors as an example vector_list = [[1, 2, 3], [4, 5, 6], [7, 8, 9]] # Compute mean vector mean_vector = utilities.vector_mean(*vector_list) # Alternative usage example (same as above): mean_vector = utilities.vector_mean([1, 2, 3], [4, 5, 6], [7, 8, 9]) :param args: list of vectors :type args: list, tuple :return: mean vector :rtype: list """ sz = len(args) mean_vector = [0.0 for _ in range(len(args[0]))] for input_vector in args: mean_vector = [a+b for a, b in zip(mean_vector, input_vector)] mean_vector = [a / sz for a in mean_vector] return mean_vector def vector_magnitude(vector_in): """ Computes the magnitude of the input vector. :param vector_in: input vector :type vector_in: list, tuple :return: magnitude of the vector :rtype: float """ sq_sum = 0.0 for vin in vector_in: sq_sum += vin**2 return math.sqrt(sq_sum) def vector_angle_between(vector1, vector2, **kwargs): """ Computes the angle between the two input vectors. If the keyword argument ``degrees`` is set to *True*, then the angle will be in degrees. Otherwise, it will be in radians. By default, ``degrees`` is set to *True*. :param vector1: vector :type vector1: list, tuple :param vector2: vector :type vector2: list, tuple :return: angle between the vectors :rtype: float """ degrees = kwargs.get('degrees', True) magn1 = vector_magnitude(vector1) magn2 = vector_magnitude(vector2) acos_val = vector_dot(vector1, vector2) / (magn1 * magn2) angle_radians = math.acos(acos_val) if degrees: return math.degrees(angle_radians) else: return angle_radians def vector_is_zero(vector_in, tol=10e-8): """ Checks if the input vector is a zero vector. :param vector_in: input vector :type vector_in: list, tuple :param tol: tolerance value :type tol: float :return: True if the input vector is zero, False otherwise :rtype: bool """ if not isinstance(vector_in, (list, tuple)): raise TypeError("Input vector must be a list or a tuple") res = [False for _ in range(len(vector_in))] for idx in range(len(vector_in)): if abs(vector_in[idx]) < tol: res[idx] = True return all(res) def point_translate(point_in, vector_in): """ Translates the input points using the input vector. :param point_in: input point :type point_in: list, tuple :param vector_in: input vector :type vector_in: list, tuple :return: translated point :rtype: list """ try: if point_in is None or len(point_in) == 0 or vector_in is None or len(vector_in) == 0: raise ValueError("Input arguments cannot be empty") except TypeError as e: print("An error occurred: {}".format(e.args[-1])) raise TypeError("Input must be a list or tuple") except Exception: raise # Translate the point using the input vector point_out = [coord + comp for coord, comp in zip(point_in, vector_in)] return point_out def point_distance(pt1, pt2): """ Computes distance between two points. :param pt1: point 1 :type pt1: list, tuple :param pt2: point 2 :type pt2: list, tuple :return: distance between input points :rtype: float """ if len(pt1) != len(pt2): raise ValueError("The input points should have the same dimension") dist_vector = vector_generate(pt1, pt2, normalize=False) distance = vector_magnitude(dist_vector) return distance def point_mid(pt1, pt2): """ Computes the midpoint of the input points. :param pt1: point 1 :type pt1: list, tuple :param pt2: point 2 :type pt2: list, tuple :return: midpoint :rtype: list """ if len(pt1) != len(pt2): raise ValueError("The input points should have the same dimension") dist_vector = vector_generate(pt1, pt2, normalize=False) half_dist_vector = vector_multiply(dist_vector, 0.5) return point_translate(pt1, half_dist_vector) @lru_cache(maxsize=os.environ['GEOMDL_CACHE_SIZE'] if "GEOMDL_CACHE_SIZE" in os.environ else 16) def matrix_identity(n): """ Generates a :math:`N \\times N` identity matrix. :param n: size of the matrix :type n: int :return: identity matrix :rtype: list """ imat = [[1.0 if i == j else 0.0 for i in range(n)] for j in range(n)] return imat def matrix_pivot(m, sign=False): """ Computes the pivot matrix for M, a square matrix. This function computes * the permutation matrix, :math:`P` * the product of M and P, :math:`M \\times P` * determinant of P, :math:`det(P)` if ``sign = True`` :param m: input matrix :type m: list, tuple :param sign: flag to return the determinant of the permutation matrix, P :type sign: bool :return: a tuple containing the matrix product of M x P, P and det(P) :rtype: tuple """ mp = deepcopy(m) n = len(mp) p = matrix_identity(n) # permutation matrix num_rowswap = 0 for j in range(0, n): row = j a_max = 0.0 for i in range(j, n): a_abs = abs(mp[i][j]) if a_abs > a_max: a_max = a_abs row = i if j != row: num_rowswap += 1 for q in range(0, n): # Swap rows p[j][q], p[row][q] = p[row][q], p[j][q] mp[j][q], mp[row][q] = mp[row][q], mp[j][q] if sign: return mp, p, math.pow(-1, num_rowswap) return mp, p def matrix_inverse(m): """ Computes the inverse of the matrix via LUP decomposition. :param m: input matrix :type m: list, tuple :return: inverse of the matrix :rtype: list """ mp, p = matrix_pivot(m) m_inv = lu_solve(mp, p) return m_inv def matrix_determinant(m): """ Computes the determinant of the square matrix :math:`M` via LUP decomposition. :param m: input matrix :type m: list, tuple :return: determinant of the matrix :rtype: float """ mp, p, sign = matrix_pivot(m, sign=True) m_l, m_u = lu_decomposition(mp) det = 1.0 for i in range(len(m)): det *= m_l[i][i] * m_u[i][i] det *= sign return det def matrix_transpose(m): """ Transposes the input matrix. The input matrix :math:`m` is a 2-dimensional array. :param m: input matrix with dimensions :math:`(n \\times m)` :type m: list, tuple :return: transpose matrix with dimensions :math:`(m \\times n)` :rtype: list """ num_cols = len(m) num_rows = len(m[0]) m_t = [] for i in range(num_rows): temp = [] for j in range(num_cols): temp.append(m[j][i]) m_t.append(temp) return m_t def matrix_multiply(mat1, mat2): """ Matrix multiplication (iterative algorithm). The running time of the iterative matrix multiplication algorithm is :math:`O(n^{3})`. :param mat1: 1st matrix with dimensions :math:`(n \\times p)` :type mat1: list, tuple :param mat2: 2nd matrix with dimensions :math:`(p \\times m)` :type mat2: list, tuple :return: resultant matrix with dimensions :math:`(n \\times m)` :rtype: list """ n = len(mat1) p1 = len(mat1[0]) p2 = len(mat2) if p1 != p2: raise GeomdlException("Column - row size mismatch") try: # Matrix - matrix multiplication m = len(mat2[0]) mat3 = [[0.0 for _ in range(m)] for _ in range(n)] for i in range(n): for j in range(m): for k in range(p2): mat3[i][j] += float(mat1[i][k] * mat2[k][j]) except TypeError: # Matrix - vector multiplication mat3 = [0.0 for _ in range(n)] for i in range(n): for k in range(p2): mat3[i] += float(mat1[i][k] * mat2[k]) return mat3 def matrix_scalar(m, sc): """ Matrix multiplication by a scalar value (iterative algorithm). The running time of the iterative matrix multiplication algorithm is :math:`O(n^{2})`. :param m: input matrix :type m: list, tuple :param sc: scalar value :type sc: int, float :return: resultant matrix :rtype: list """ mm = [[0.0 for _ in range(len(m[0]))] for _ in range(len(m))] for i in range(len(m)): for j in range(len(m[0])): mm[i][j] = float(m[i][j] * sc) return mm def triangle_normal(tri): """ Computes the (approximate) normal vector of the input triangle. :param tri: triangle object :type tri: elements.Triangle :return: normal vector of the triangle :rtype: tuple """ vec1 = vector_generate(tri.vertices[0].data, tri.vertices[1].data) vec2 = vector_generate(tri.vertices[1].data, tri.vertices[2].data) return vector_cross(vec1, vec2) def triangle_center(tri, uv=False): """ Computes the center of mass of the input triangle. :param tri: triangle object :type tri: elements.Triangle :param uv: if True, then finds parametric position of the center of mass :type uv: bool :return: center of mass of the triangle :rtype: tuple """ if uv: data = [t.uv for t in tri] mid = [0.0, 0.0] else: data = tri.vertices mid = [0.0, 0.0, 0.0] for vert in data: mid = [m + v for m, v in zip(mid, vert)] mid = [float(m) / 3.0 for m in mid] return tuple(mid) @lru_cache(maxsize=os.environ['GEOMDL_CACHE_SIZE'] if "GEOMDL_CACHE_SIZE" in os.environ else 128) def binomial_coefficient(k, i): """ Computes the binomial coefficient (denoted by *k choose i*). Please see the following website for details: http://mathworld.wolfram.com/BinomialCoefficient.html :param k: size of the set of distinct elements :type k: int :param i: size of the subsets :type i: int :return: combination of *k* and *i* :rtype: float """ # Special case if i > k: return float(0) # Compute binomial coefficient k_fact = math.factorial(k) i_fact = math.factorial(i) k_i_fact = math.factorial(k - i) return float(k_fact / (k_i_fact * i_fact)) def lu_decomposition(matrix_a): """ LU-Factorization method using Doolittle's Method for solution of linear systems. Decomposes the matrix :math:`A` such that :math:`A = LU`. The input matrix is represented by a list or a tuple. The input matrix is **2-dimensional**, i.e. list of lists of integers and/or floats. :param matrix_a: Input matrix (must be a square matrix) :type matrix_a: list, tuple :return: a tuple containing matrices L and U :rtype: tuple """ # Check if the 2-dimensional input matrix is a square matrix q = len(matrix_a) for idx, m_a in enumerate(matrix_a): if len(m_a) != q: raise ValueError("The input must be a square matrix. " + "Row " + str(idx + 1) + " has a size of " + str(len(m_a)) + ".") # Return L and U matrices return _linalg.doolittle(matrix_a) def forward_substitution(matrix_l, matrix_b): """ Forward substitution method for the solution of linear systems. Solves the equation :math:`Ly = b` using forward substitution method where :math:`L` is a lower triangular matrix and :math:`b` is a column matrix. :param matrix_l: L, lower triangular matrix :type matrix_l: list, tuple :param matrix_b: b, column matrix :type matrix_b: list, tuple :return: y, column matrix :rtype: list """ q = len(matrix_b) matrix_y = [0.0 for _ in range(q)] matrix_y[0] = float(matrix_b[0]) / float(matrix_l[0][0]) for i in range(1, q): matrix_y[i] = float(matrix_b[i]) - sum([matrix_l[i][j] * matrix_y[j] for j in range(0, i)]) matrix_y[i] /= float(matrix_l[i][i]) return matrix_y def backward_substitution(matrix_u, matrix_y): """ Backward substitution method for the solution of linear systems. Solves the equation :math:`Ux = y` using backward substitution method where :math:`U` is a upper triangular matrix and :math:`y` is a column matrix. :param matrix_u: U, upper triangular matrix :type matrix_u: list, tuple :param matrix_y: y, column matrix :type matrix_y: list, tuple :return: x, column matrix :rtype: list """ q = len(matrix_y) matrix_x = [0.0 for _ in range(q)] matrix_x[q - 1] = float(matrix_y[q - 1]) / float(matrix_u[q - 1][q - 1]) for i in range(q - 2, -1, -1): matrix_x[i] = float(matrix_y[i]) - sum([matrix_u[i][j] * matrix_x[j] for j in range(i, q)]) matrix_x[i] /= float(matrix_u[i][i]) return matrix_x def lu_solve(matrix_a, b): """ Computes the solution to a system of linear equations. This function solves :math:`Ax = b` using LU decomposition. :math:`A` is a :math:`N \\times N` matrix, :math:`b` is :math:`N \\times M` matrix of :math:`M` column vectors. Each column of :math:`x` is a solution for corresponding column of :math:`b`. :param matrix_a: matrix A :type m_l: list :param b: matrix of M column vectors :type b: list :return: x, the solution matrix :rtype: list """ # Variable initialization dim = len(b[0]) num_x = len(b) x = [[0.0 for _ in range(dim)] for _ in range(num_x)] # LU decomposition m_l, m_u = lu_decomposition(matrix_a) # Solve the system of linear equations for i in range(dim): bt = [b1[i] for b1 in b] y = forward_substitution(m_l, bt) xt = backward_substitution(m_u, y) for j in range(num_x): x[j][i] = xt[j] # Return the solution return x def lu_factor(matrix_a, b): """ Computes the solution to a system of linear equations with partial pivoting. This function solves :math:`Ax = b` using LUP decomposition. :math:`A` is a :math:`N \\times N` matrix, :math:`b` is :math:`N \\times M` matrix of :math:`M` column vectors. Each column of :math:`x` is a solution for corresponding column of :math:`b`. :param matrix_a: matrix A :type m_l: list :param b: matrix of M column vectors :type b: list :return: x, the solution matrix :rtype: list """ # Variable initialization dim = len(b[0]) num_x = len(b) x = [[0.0 for _ in range(dim)] for _ in range(num_x)] # LUP decomposition mp, p = matrix_pivot(matrix_a) m_l, m_u = lu_decomposition(mp) # Solve the system of linear equations for i in range(dim): bt = [b1[i] for b1 in b] y = forward_substitution(m_l, bt) xt = backward_substitution(m_u, y) for j in range(num_x): x[j][i] = xt[j] # Return the solution return x def linspace(start, stop, num, decimals=18): """ Returns a list of evenly spaced numbers over a specified interval. Inspired from Numpy's linspace function: https://github.com/numpy/numpy/blob/master/numpy/core/function_base.py :param start: starting value :type start: float :param stop: end value :type stop: float :param num: number of samples to generate :type num: int :param decimals: number of significands :type decimals: int :return: a list of equally spaced numbers :rtype: list """ start = float(start) stop = float(stop) if abs(start - stop) <= 10e-8: return [start] num = int(num) if num > 1: div = num - 1 delta = stop - start return [float(("{:." + str(decimals) + "f}").format((start + (float(x) * float(delta) / float(div))))) for x in range(num)] return [float(("{:." + str(decimals) + "f}").format(start))] def frange(start, stop, step=1.0): """ Implementation of Python's ``range()`` function which works with floats. Reference to this implementation: https://stackoverflow.com/a/36091634 :param start: start value :type start: float :param stop: end value :type stop: float :param step: increment :type step: float :return: float :rtype: generator """ i = 0.0 x = float(start) # Prevent yielding integers. x0 = x epsilon = step / 2.0 yield x # always yield first value while x + epsilon < stop: i += 1.0 x = x0 + i * step yield x if stop > x: yield stop # for yielding last value of the knot vector if the step is a large value, like 0.1 def convex_hull(points): """ Returns points on convex hull in counterclockwise order according to Graham's scan algorithm. Reference: https://gist.github.com/arthur-e/5cf52962341310f438e96c1f3c3398b8 .. note:: This implementation only works in 2-dimensional space. :param points: list of 2-dimensional points :type points: list, tuple :return: convex hull of the input points :rtype: list """ turn_left, turn_right, turn_none = (1, -1, 0) def cmp(a, b): return (a > b) - (a < b) def turn(p, q, r): return cmp((q[0] - p[0])*(r[1] - p[1]) - (r[0] - p[0])*(q[1] - p[1]), 0) def keep_left(hull, r): while len(hull) > 1 and turn(hull[-2], hull[-1], r) != turn_left: hull.pop() if not len(hull) or hull[-1] != r: hull.append(r) return hull points = sorted(points) l = reduce(keep_left, points, []) u = reduce(keep_left, reversed(points), []) return l.extend(u[i] for i in range(1, len(u) - 1)) or l def is_left(point0, point1, point2): """ Tests if a point is Left|On|Right of an infinite line. Ported from the C++ version: on http://geomalgorithms.com/a03-_inclusion.html .. note:: This implementation only works in 2-dimensional space. :param point0: Point P0 :param point1: Point P1 :param point2: Point P2 :return: >0 for P2 left of the line through P0 and P1 =0 for P2 on the line <0 for P2 right of the line """ return ((point1[0] - point0[0]) * (point2[1] - point0[1])) - ((point2[0] - point0[0]) * (point1[1] - point0[1])) def wn_poly(point, vertices): """ Winding number test for a point in a polygon. Ported from the C++ version: http://geomalgorithms.com/a03-_inclusion.html .. note:: This implementation only works in 2-dimensional space. :param point: point to be tested :type point: list, tuple :param vertices: vertex points of a polygon vertices[n+1] with vertices[n] = vertices[0] :type vertices: list, tuple :return: True if the point is inside the input polygon, False otherwise :rtype: bool """ wn = 0 # the winding number counter v_size = len(vertices) - 1 # loop through all edges of the polygon for i in range(v_size): # edge from V[i] to V[i+1] if vertices[i][1] <= point[1]: # start y <= P.y if vertices[i + 1][1] > point[1]: # an upward crossing if is_left(vertices[i], vertices[i + 1], point) > 0: # P left of edge wn += 1 # have a valid up intersect else: # start y > P.y (no test needed) if vertices[i + 1][1] <= point[1]: # a downward crossing if is_left(vertices[i], vertices[i + 1], point) < 0: # P right of edge wn -= 1 # have a valid down intersect # return wn return bool(wn)