Source code for pyrigi.misc

"""
This module provides various following miscellaneous functions.
"""

from math import isclose, log10

import numpy as np
import sympy as sp
from sympy import Matrix, MatrixBase

from pyrigi.data_type import Sequence, Number, InfFlex, Vertex, Point


try:
    from IPython.core.magic import register_cell_magic

    @register_cell_magic
    def skip_execution(line, cell):
        print(
            "This cell was marked to be skipped (probably due to long execution time)."
        )
        print("Remove the cell magic `%%skip_execution` to run it.")
        return

except NameError:
    pass


def _doc_category(category):
    """
    Decorator for doc categories.
    """

    def decorator_doc_category(func):
        setattr(func, "_doc_category", category)
        return func

    return decorator_doc_category


def _generate_category_tables(
    cls, tabs, cat_order=None, include_all=False, add_attributes=True
) -> str:
    """
    Generate a formatted string that categorizes methods of a given class.

    Parameters
    ----------
    cls:
        A class.
    tabs:
        The number of indentation levels that are applied to the output.
    cat_order:
        Optional list specifying the order in which categories appear
        in the output.
    include_all:
        Optional boolean determining whether methods without a specific category
        should be included.
    add_attributes:
        Optional boolean determining whether the public attributes should
        be listed among attribute getters.
    """
    if cat_order is None:
        cat_order = []
    categories = {}
    for func in dir(cls):
        if callable(getattr(cls, func)) and func[:2] != "__":
            f = getattr(cls, func)
            if hasattr(f, "_doc_category"):
                category = f._doc_category
                if category in categories:
                    categories[category].append(func)
                else:
                    categories[category] = [func]
            elif include_all:
                if "Not categorized" in categories:
                    categories["Not categorized"].append(func)
                else:
                    categories["Not categorized"] = [func]
        elif isinstance(getattr(cls, func), property) and add_attributes:
            category = "Attribute getters"
            if category in categories:
                categories[category].append(func)
            else:
                categories[category] = [func]

    for category in categories:
        if category not in cat_order:
            cat_order.append(category)

    res = "Methods\n-------\n"
    for category, functions in sorted(
        categories.items(), key=lambda t: cat_order.index(t[0])
    ):
        res += f"**{category}**"
        res += "\n\n.. autosummary::\n\n    "
        res += "\n    ".join(functions)
        res += "\n\n"
    indent = "".join(["    " for _ in range(tabs)])
    return ("\n" + indent).join(res.splitlines())


def _generate_two_orthonormal_vectors(dim: int, random_seed: int = None) -> Matrix:
    """
    Generate two random numeric orthonormal vectors in the given dimension.

    The vectors are in the columns of the returned matrix.

    Parameters
    ----------
    dim:
        The dimension in which the vectors are generated.
    random_seed:
        Seed for generating random vectors.
        When the same value is provided, the same vectors are generated.
    """

    if random_seed is not None:
        np.random.seed(random_seed)

    matrix = np.random.randn(dim, 2)

    # for numerical stability regenerate some elements
    tmp = np.random.randint(0, dim - 1)
    while abs(matrix[tmp, 1]) < 1e-6:
        matrix[tmp, 1] = np.random.randn(1, 1)

    while abs(matrix[-1, 0]) < 1e-6:
        matrix[-1, 0] = np.random.randn(1, 1)

    tmp = np.dot(matrix[:-1, 0], matrix[:-1, 1]) * -1
    matrix[-1, 1] = tmp / matrix[-1, 0]

    # normalize
    matrix[:, 0] = matrix[:, 0] / np.linalg.norm(matrix[:, 0])
    matrix[:, 1] = matrix[:, 1] / np.linalg.norm(matrix[:, 1])
    return matrix


def _generate_three_orthonormal_vectors(dim: int, random_seed: int = None) -> Matrix:
    """
    Generate three random numeric orthonormal vectors in the given dimension.

    Parameters
    ----------
    dim:
        The dimension in which the vectors are generated.
    random_seed:
        Seed for generating random vectors.
        When the same value is provided, the same vectors are generated.

    Notes
    -----
    The vectors are in the columns of the returned matrix. To ensure that the
    vectors are uniformly distributed over the Stiefel manifold, we need to
    ensure that the triangular matrix `R` has positive diagonal elements.
    """

    if random_seed is not None:
        np.random.seed(random_seed)

    matrix = np.random.randn(dim, 3)
    Q, R = np.linalg.qr(matrix)
    return Q @ np.diag(np.sign(np.diag(R)))


[docs] def is_zero(expr: Number, numerical: bool = False, tolerance: float = 1e-9) -> bool: """ Return if the given expression is zero. Parameters ---------- expr: Expression that is checked. numerical: If ``True``, then the check is done only numerically with the given tolerance. If ``False`` (default), the check is done symbolically, ``sympy`` method ``equals`` is used. tolerance: The tolerance that is used in the numerical check coordinate-wise. """ if not numerical: zero_bool = sp.cancel(sp.sympify(expr)).equals(0) if zero_bool is None: raise RuntimeError( "It could not be determined by sympy " + "whether the given sympy expression is zero." + "Please report this as an issue on Github " + "(https://github.com/PyRigi/PyRigi/issues)." ) return zero_bool else: return isclose( sympy_expr_to_float(expr, tolerance=tolerance), 0, abs_tol=tolerance, )
[docs] def is_zero_vector( vector: Sequence[Number], numerical: bool = False, tolerance: float = 1e-9 ) -> bool: """ Return if the given vector is zero. Parameters ---------- vector: Vector that is checked. numerical: If ``True``, then the check is done only numerically with the given tolerance. If ``False`` (default), the check is done symbolically, ``sympy`` attribute ``is_zero`` is used. tolerance: The tolerance that is used in the numerical check coordinate-wise. """ if not isinstance(vector, Matrix): vector = point_to_vector(vector) return all( [is_zero(coord, numerical=numerical, tolerance=tolerance) for coord in vector] )
[docs] def sympy_expr_to_float( expression: Sequence[Number] | Matrix | Number, tolerance: float = 1e-9 ) -> list[float] | float: """ Convert a sympy expression to (numerical) floats. If the given ``expression`` is a ``Sequence`` of ``Numbers`` or a ``Matrix``, then each individual element is evaluated and a list of ``float`` is returned. If the input is just a single sympy expression, it is evaluated and returned as a ``float``. Parameters ---------- expression: The sympy expression. tolerance: Intended level of numerical accuracy. Notes ----- The method :func:`.data_type.point_to_vector` is used to ensure that the input is consistent with the sympy format. """ try: if isinstance(expression, list | tuple | Matrix): return [ float( sp.sympify(coord).evalf( int(round(2.5 * log10(tolerance ** (-1) + 1))) ) ) for coord in point_to_vector(expression) ] return float( sp.sympify(expression).evalf(int(round(2.5 * log10(tolerance ** (-1) + 1)))) ) except sp.SympifyError: raise ValueError(f"The expression `{expression}` could not be parsed by sympy.")
def _normalize_flex( inf_flex: InfFlex, numerical: bool = False, tolerance: float = 1e-9 ) -> InfFlex: """ Divide a vector by its Euclidean norm. Parameters ---------- inf_flex: The infinitesimal flex that is supposed to be normalized. numerical: Boolean determining whether a numerical or symbolic normalization is performed. tolerance: Intended level of numerical accuracy. """ if isinstance(inf_flex, dict): if numerical: _inf_flex = { v: sympy_expr_to_float(flex, tolerance=tolerance) for v, flex in inf_flex.items() } flex_norm = np.linalg.norm(sum(_inf_flex.values(), [])) if isclose(flex_norm, 0, abs_tol=tolerance): raise ValueError("The norm of this flex is almost zero.") return { v: tuple([q / flex_norm for q in flex]) for v, flex in _inf_flex.items() } flex_norm = sp.sqrt(sum([q**2 for flex in inf_flex.values() for q in flex])) if is_zero(flex_norm, numerical=numerical, tolerance=tolerance): raise ValueError("The norm of this flex is zero.") return {v: tuple([q / flex_norm for q in flex]) for v, flex in inf_flex.items()} elif isinstance(inf_flex, Sequence): if numerical: _inf_flex = [ sympy_expr_to_float(flex, tolerance=tolerance) for flex in inf_flex ] flex_norm = np.linalg.norm(_inf_flex) if isclose(flex_norm, 0, abs_tol=tolerance): raise ValueError("The norm of this flex is almost zero.") return [flex / flex_norm for flex in _inf_flex] flex_norm = sp.sqrt(sum([flex**2 for flex in inf_flex])) if is_zero(flex_norm, numerical=numerical, tolerance=tolerance): raise ValueError("The norm of this flex is zero.") return [flex / flex_norm for flex in inf_flex] else: raise TypeError("`inf_flex` does not have the correct type.") def _vector_distance_pointwise( dict1: dict[Vertex, Sequence[Number]], dict2: dict[Vertex, Sequence[Number]], numerical: bool = False, ) -> float: """ Compute the Euclidean distance between two realizations or pointwise vectors. This method computes the Euclidean distance from the realization ``dict_1`` to ``dict2`` considering them as vectors. The keys of ``dict1`` and ``dict2`` must be the same. Parameters ---------- dict1, dict2: The dictionaries that are used for the distance computation. numerical: Boolean determining whether a numerical or symbolic normalization is performed. """ if not set(dict1.keys()) == set(dict2.keys()) or not len(dict1) == len(dict2): raise ValueError("`dict1` and `dict2` are not based on the same vertex set.") if numerical: return float( np.linalg.norm( [ x - y for v in dict1.keys() for x, y in zip( dict1[v], dict2[v], ) ] ) ) return sp.sqrt( sum( [ (x - y) ** 2 for v in dict1.keys() for x, y in zip( dict1[v], dict2[v], ) ] ) )
[docs] def is_isomorphic_graph_list(list1, list2) -> bool: """ Return whether two lists of graphs are the same up to graph isomorphism. """ if len(list1) != len(list2): return False for graph1 in list1: count_copies = 0 for grapht in list1: if graph1.is_isomorphic(grapht): count_copies += 1 count_found = 0 for graph2 in list2: if graph1.is_isomorphic(graph2): count_found += 1 if count_found == count_copies: break else: return False return True
[docs] def point_to_vector(point: Point) -> Matrix: """ Return point as single column sympy Matrix. """ if isinstance(point, MatrixBase) or isinstance(point, np.ndarray): if ( len(point.shape) > 1 and point.shape[0] != 1 and point.shape[1] != 1 ) or len(point.shape) > 2: raise ValueError("Point could not be interpreted as column vector.") if isinstance(point, np.ndarray): point = np.array([point]) if len(point.shape) == 1 else point point = Matrix( [ [float(point[i, j]) for i in range(point.shape[0])] for j in range(point.shape[1]) ] ) return point if (point.shape[1] == 1) else point.transpose() if not isinstance(point, Sequence) or isinstance(point, str): raise TypeError("The point must be a Sequence of Numbers.") try: res = Matrix(point) except Exception as e: raise ValueError("A coordinate could not be interpreted by sympify:\n" + str(e)) if res.shape[0] != 1 and res.shape[1] != 1: raise ValueError("Point could not be interpreted as column vector.") return res if (res.shape[1] == 1) else res.transpose()
def _null_space(A: np.array, tolerance: float = 1e-8) -> np.array: """ Compute the kernel of a numpy matrix. Parameters ---------- tolerance: Used tolerance for the selection of the vectors in the kernel of the numerical matrix. """ _, s, vh = np.linalg.svd(A, full_matrices=True) tol = np.amax(s) * tolerance num = np.sum(s > tol, dtype=int) Q = vh[num:, :].T.conj() return Q